1 #include "mrilib.h"
2 
3 /*********************************************************************
4   * Wrapper functions to interpolate images in various ways.
5   * Built on top of functions in mri_genalign_util (thus the GA_).
6   * Could be optimized for various sub-cases, but is that
7     really worth the effort? We'll see.
8 *********************************************************************/
9 
10 /*---------------------------------------------------------------------------*/
11 
mri_interp_scalar_to_floats_pointset(MRI_IMAGE * fim,int npp,float * ip,float * jp,float * kp,int code,float * vv)12 void mri_interp_scalar_to_floats_pointset( MRI_IMAGE *fim,
13                                            int npp, float *ip, float *jp, float *kp,
14                                            int code, float *vv )
15 {
16    MRI_IMAGE *inim=fim ;
17 
18 ENTRY("mri_interp_scalar_to_floats_pointset") ;
19 
20    if( inim == NULL || npp <= 0    ||
21        ip   == NULL || jp  == NULL || kp == NULL || vv == NULL ){
22      ERROR_message("NULL inputs to mri_interp_scalar_to_floats_pointset()") ;
23      EXRETURN ;
24    }
25 
26    /* convert to float? */
27 
28    if( inim->kind != MRI_float ) inim = mri_to_float(fim) ;
29 
30    switch( code ){
31 
32      case MRI_NN:
33        GA_interp_NN     ( inim , npp,ip,jp,kp , vv ) ;
34      break ;
35 
36      case MRI_LINEAR:
37        GA_interp_linear ( inim , npp,ip,jp,kp , vv ) ;
38      break ;
39 
40      case MRI_CUBIC:
41        GA_interp_cubic  ( inim , npp,ip,jp,kp , vv ) ;
42      break ;
43 
44      case MRI_QUINTIC:
45        GA_interp_quintic( inim , npp,ip,jp,kp , vv ) ;
46      break ;
47 
48      case MRI_WSINC5:
49        GA_interp_wsinc5 ( inim , npp,ip,jp,kp , vv ) ;
50      break ;
51 
52      default:
53        ERROR_message("unknown interp code = %d in mri_interp_to_floats_pointset()",code) ;
54 
55    }
56 
57    if( inim != fim ) mri_free(inim) ;
58 
59    EXRETURN ;
60 }
61 
62 /*---------------------------------------------------------------------------*/
63 
mri_interp_to_sametyp_pointset(MRI_IMAGE * fim,int npp,float * ip,float * jp,float * kp,int code,void * vv)64 void mri_interp_to_sametyp_pointset( MRI_IMAGE *fim,
65                                      int npp, float *ip, float *jp, float *kp,
66                                      int code, void *vv )
67 {
68    register int ii ;
69 
70 ENTRY("mri_interp_to_sametyp_pointset") ;
71 
72    if( fim == NULL || npp <= 0    ||
73        ip  == NULL || jp  == NULL || kp == NULL || vv == NULL ){
74      ERROR_message("NULL inputs to mri_interp_to_sametyp_pointset()") ;
75      EXRETURN ;
76    }
77 
78    switch( fim->kind ){
79 
80      case MRI_float:
81        mri_interp_scalar_to_floats_pointset( fim,
82                                              npp, ip,jp,kp, code,(float *)vv ) ;
83      break ;
84 
85      case MRI_int:{
86        int *pp = (int *)vv ;
87        float *ff = (float *)malloc(sizeof(float)*npp) ;
88        mri_interp_scalar_to_floats_pointset( fim ,
89                                              npp , ip,jp,kp , code,ff ) ;
90        for( ii=0 ; ii < npp ; ii++ ) pp[ii] = (int)rintf(ff[ii]) ;
91        free(ff) ;
92      }
93      break ;
94 
95      case MRI_short:{
96        short *pp = (short *)vv ;
97        float *ff = (float *)malloc(sizeof(float)*npp) ;
98        mri_interp_scalar_to_floats_pointset( fim ,
99                                              npp , ip,jp,kp , code,ff ) ;
100        for( ii=0 ; ii < npp ; ii++ ) pp[ii] = SHORTIZE(ff[ii]) ;
101        free(ff) ;
102      }
103      break ;
104 
105      case MRI_byte:{
106        byte *pp = (byte *)vv ;
107        float *ff = (float *)malloc(sizeof(float)*npp) ;
108        mri_interp_scalar_to_floats_pointset( fim ,
109                                              npp , ip,jp,kp , code,ff ) ;
110        for( ii=0 ; ii < npp ; ii++ ) pp[ii] = BYTEIZE(ff[ii]) ;
111        free(ff) ;
112      }
113      break ;
114 
115      case MRI_rgb:{
116        byte *pval = (byte *)vv ;
117        MRI_IMARR *ppar ;
118        float *rr = (float *)malloc(sizeof(float)*npp) ;
119        float *gg = (float *)malloc(sizeof(float)*npp) ;
120        float *bb = (float *)malloc(sizeof(float)*npp) ;
121        ppar = mri_rgb_to_3float(fim) ;
122        mri_interp_scalar_to_floats_pointset( IMARR_SUBIM(ppar,0) ,
123                                              npp , ip,jp,kp , code, rr ) ;
124        mri_interp_scalar_to_floats_pointset( IMARR_SUBIM(ppar,1) ,
125                                              npp , ip,jp,kp , code, gg ) ;
126        mri_interp_scalar_to_floats_pointset( IMARR_SUBIM(ppar,2) ,
127                                              npp , ip,jp,kp , code, bb ) ;
128        for( ii=0 ; ii < npp ; ii++ ){
129          pval[3*ii+0] = BYTEIZE(rr[ii]) ;
130          pval[3*ii+1] = BYTEIZE(gg[ii]) ;
131          pval[3*ii+2] = BYTEIZE(bb[ii]) ;
132        }
133        free(bb) ; free(gg) ; free(rr) ;
134      }
135      break ;
136 
137      case MRI_rgba:{
138        rgba *pval = (rgba *)vv ;
139        MRI_IMARR *ppar ;
140        float *rr = (float *)malloc(sizeof(float)*npp) ;
141        float *gg = (float *)malloc(sizeof(float)*npp) ;
142        float *bb = (float *)malloc(sizeof(float)*npp) ;
143        float *aa = (float *)malloc(sizeof(float)*npp) ;
144        ppar = mri_rgba_to_4float(fim) ;
145        mri_interp_scalar_to_floats_pointset( IMARR_SUBIM(ppar,0) ,
146                                              npp , ip,jp,kp , code, rr ) ;
147        mri_interp_scalar_to_floats_pointset( IMARR_SUBIM(ppar,1) ,
148                                              npp , ip,jp,kp , code, gg ) ;
149        mri_interp_scalar_to_floats_pointset( IMARR_SUBIM(ppar,2) ,
150                                              npp , ip,jp,kp , code, bb ) ;
151        mri_interp_scalar_to_floats_pointset( IMARR_SUBIM(ppar,3) ,
152                                              npp , ip,jp,kp , code, aa ) ;
153        for( ii=0 ; ii < npp ; ii++ ){
154          pval[ii].r = BYTEIZE(rr[ii]) ;
155          pval[ii].g = BYTEIZE(gg[ii]) ;
156          pval[ii].b = BYTEIZE(bb[ii]) ;
157          pval[ii].a = BYTEIZE(rr[ii]) ;
158        }
159        free(aa) ; free(bb) ; free(gg) ; free(rr) ;
160      }
161      break ;
162 
163      case MRI_complex:{
164        complex *pval = (complex *)vv ;
165        MRI_IMARR *ppar ;
166        float *xx = (float *)malloc(sizeof(float)*npp) ;
167        float *yy = (float *)malloc(sizeof(float)*npp) ;
168        ppar =  mri_complex_to_pair(fim) ;
169        mri_interp_scalar_to_floats_pointset( IMARR_SUBIM(ppar,0) ,
170                                              npp , ip,jp,kp , code, xx ) ;
171        mri_interp_scalar_to_floats_pointset( IMARR_SUBIM(ppar,1) ,
172                                              npp , ip,jp,kp , code, yy ) ;
173        for( ii=0 ; ii < npp ; ii++ ){
174          pval[ii].r = xx[ii] ;
175          pval[ii].i = yy[ii] ;
176        }
177        free(yy) ; free(xx) ;
178      }
179      break ;
180 
181      case MRI_fvect:{
182        MRI_IMARR *ppar ; int nff, kk,ii,qq ; float **ffar, *pval=(float *)vv ;
183        ppar = mri_fvect_to_imarr(fim) ;
184        nff  = IMARR_COUNT(ppar) ;
185        ffar = (float **)malloc(sizeof(float *)*nff) ;
186        for( kk=0 ; kk < nff ; kk++ ){
187          ffar[kk] = (float *)malloc(sizeof(float)*npp) ;
188          mri_interp_scalar_to_floats_pointset( IMARR_SUBIM(ppar,kk) ,
189                                                npp , ip,jp,kp , code, ffar[kk] ) ;
190        }
191        for( qq=ii=0 ; ii < npp ; ii++ ){
192          for( kk=0 ; kk < nff ; kk++ ) pval[qq++] = ffar[kk][ii] ;
193        }
194        for( kk=0 ; kk < nff ; kk++ ) free(ffar[kk]) ;
195        free(ffar) ;
196      }
197      break ;
198 
199      default:
200        ERROR_message("unknown image kind = %d in mri_interp_to_sametyp_pointset()",fim->kind) ;
201      break ;
202 
203    }
204 
205    EXRETURN ;
206 }
207 
208 /*---------------------------------------------------------------------------*/
209 
mri_interp_to_floats_block(MRI_IMAGE * fim,mat44 ijk_to_ijk,int code,int ibot,int itop,int jbot,int jtop,int kbot,int ktop)210 MRI_IMAGE * mri_interp_to_floats_block( MRI_IMAGE *fim ,
211                                         mat44 ijk_to_ijk , int code ,
212                                         int ibot , int itop ,
213                                         int jbot , int jtop ,
214                                         int kbot , int ktop  )
215 {
216    MRI_IMAGE *outim ; float *outar ;
217    int nxout, nyout, nzout, nxyz , ii,jj,kk,qq ;
218    float *ip , *jp , *kp ;
219 
220 ENTRY("mri_interp_to_floats_block") ;
221 
222    if( fim == NULL || !ISVALID_MAT44(ijk_to_ijk) ||
223        ibot > itop || jbot > jtop || kbot > ktop   ) RETURN(NULL) ;
224 
225    nxout = itop-ibot+1 ;
226    nyout = jtop-jbot+1 ;
227    nzout = ktop-kbot+1 ; nxyz = nxout * nyout * nzout ;
228 
229    outim = mri_new_vol( nxout, nyout, nzout , MRI_float ) ;
230    outar = MRI_FLOAT_PTR(outim) ;
231 
232    ip = (float *)malloc(sizeof(float)*nxyz) ;
233    jp = (float *)malloc(sizeof(float)*nxyz) ;
234    kp = (float *)malloc(sizeof(float)*nxyz) ;
235 
236    for( qq=kk=kbot ; kk <= ktop ; kk++ ){
237     for( jj=jbot ; jj <= jtop ; jj++ ){
238      for( ii=ibot ; ii <= itop ; ii++,qq++ ){
239        MAT44_VEC( ijk_to_ijk , ii,jj,kk , ip[qq],jp[qq],kp[qq] ) ;
240    }}}
241 
242    mri_interp_scalar_to_floats_pointset( fim, nxyz, ip,jp,kp, code, outar ) ;
243 
244    free(kp); free(jp); free(ip) ;
245 
246    RETURN(outim) ;
247 }
248 
249 /*---------------------------------------------------------------------------*/
250 
mri_interp_to_vectyp_block(MRI_IMAGE * fim,mat44 ijk_to_ijk,int code,int ibot,int itop,int jbot,int jtop,int kbot,int ktop)251 MRI_IMAGE * mri_interp_to_vectyp_block( MRI_IMAGE *fim ,
252                                         mat44 ijk_to_ijk , int code ,
253                                         int ibot , int itop ,
254                                         int jbot , int jtop ,
255                                         int kbot , int ktop  )
256 {
257    MRI_IMAGE *outim ; float *outar ;
258    int nxout, nyout, nzout, nxyz , ii,jj,kk,qq ;
259    float *ip , *jp , *kp ;
260 
261 ENTRY("mri_interp_to_vectyp_block") ;
262 
263    if( fim == NULL || !ISVALID_MAT44(ijk_to_ijk) ||
264        ibot > itop || jbot > jtop || kbot > ktop   ) RETURN(NULL) ;
265 
266    if( ! ISVECTIM(fim) ){
267      outim = mri_interp_to_floats_block( fim,ijk_to_ijk,code,ibot,itop,jbot,jtop,kbot,ktop) ;
268      RETURN(outim) ;
269    }
270 
271    nxout = itop-ibot+1 ;
272    nyout = jtop-jbot+1 ;
273    nzout = ktop-kbot+1 ; nxyz = nxout * nyout * nzout ;
274 
275    if( fim->kind == MRI_fvect ){
276      outim = mri_new_fvectim( nxout , nyout , nzout , fim->vdim ) ;
277    } else {
278      outim = mri_new_vol( nxout, nyout, nzout , fim->kind ) ;
279    }
280    outar = mri_data_pointer(outim) ;
281 
282    ip = (float *)malloc(sizeof(float)*nxyz) ;
283    jp = (float *)malloc(sizeof(float)*nxyz) ;
284    kp = (float *)malloc(sizeof(float)*nxyz) ;
285 
286    for( qq=kk=kbot ; kk <= ktop ; kk++ ){
287     for( jj=jbot ; jj <= jtop ; jj++ ){
288      for( ii=ibot ; ii <= itop ; ii++,qq++ ){
289        MAT44_VEC( ijk_to_ijk , ii,jj,kk , ip[qq],jp[qq],kp[qq] ) ;
290    }}}
291 
292    mri_interp_to_sametyp_pointset( fim, nxyz, ip,jp,kp, code, outar ) ;
293 
294    free(kp); free(jp); free(ip) ;
295 
296    RETURN(outim) ;
297 }
298 
299 /*---------------------------------------------------------------------------*/
300 /* Interpolate 1 point from each image, and return a 1D float array.
301    This function is meant to provide a time-series array from an
302    interpolated grid point, for use in the AFNI graph viewer.
303    It is rather brute force, but as computers are so fast these days .... */
304 /*---------------------------------------------------------------------------*/
305 
imarr_interp_scalar_to_floats_onepoint(MRI_IMARR * imar,float iq,float jq,float kq,int code)306 MRI_IMAGE * imarr_interp_scalar_to_floats_onepoint(
307                     MRI_IMARR *imar,
308                     float iq , float jq , float kq , int code )
309 {
310    MRI_IMAGE *inim, *fim, *outim ; int nim, kk ;
311    float *outar , vv[1] , ip,jp,kp ;
312 
313 ENTRY("imarr_interp_scalar_to_floats_onepoint") ;
314 
315    if( imar == NULL || IMARR_COUNT(imar) < 1 ){
316      ERROR_message("NULL inputs to imarr_interp_scalar_to_floats_onepoint()") ;
317      RETURN(NULL) ;
318    }
319 
320    nim = IMARR_COUNT(imar) ;
321    outim = mri_new( nim , 1 , MRI_float ) ;
322    outar = MRI_FLOAT_PTR(outim) ;
323    ip = iq ; jp = jq ; kp = kq ;
324 
325    for( kk=0 ; kk < nim ; kk++ ){
326 
327      inim = IMARR_SUBIM(imar,kk) ;
328      if( inim == NULL ) continue ;
329 
330      /* convert to float? */
331      if( inim->kind != MRI_float ) inim = mri_to_float(fim) ;
332 
333      switch( code ){
334        case MRI_NN:
335          GA_interp_NN     ( inim , 1,&ip,&jp,&kp , vv ) ;
336        break ;
337        case MRI_LINEAR:
338          GA_interp_linear ( inim , 1,&ip,&jp,&kp , vv ) ;
339        break ;
340        case MRI_CUBIC:
341          GA_interp_cubic  ( inim , 1,&ip,&jp,&kp , vv ) ;
342        break ;
343        case MRI_QUINTIC:
344          GA_interp_quintic( inim , 1,&ip,&jp,&kp , vv ) ;
345        break ;
346        case MRI_WSINC5:
347          GA_interp_wsinc5 ( inim , 1,&ip,&jp,&kp , vv ) ;
348        break ;
349      }
350      outar[kk] = vv[0] ;
351      if( inim != fim ) mri_free(inim) ;
352    }
353 
354    RETURN(outim) ;
355 }
356