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