1 #include "mrilib.h"
2
3 #ifdef USE_OMP
4 #include <omp.h>
5 #endif
6
7 /*** NOT 7D SAFE ***/
8
9 /****************************************************************************/
10 /** functions to "warp" 3D images -- not very efficient, but quite general **/
11 /****************************************************************************/
12
13 /*--------------------------------------------------------------------------*/
14
15 static int wtype = MRI_LINEAR ;
mri_warp3D_method(int mode)16 void mri_warp3D_method( int mode ){ wtype = mode ; } /** set interpolation **/
17
18 /*--------------------------------------------------------------------------*/
19
20 static byte *womask = NULL ; /* 19 Nov 2004 */
21
mri_warp3D_set_womask(MRI_IMAGE * wim)22 void mri_warp3D_set_womask( MRI_IMAGE *wim )
23 {
24 womask = (wim == NULL || wim->kind != MRI_byte) ? NULL
25 : MRI_BYTE_PTR(wim) ;
26 }
27
28 #undef SKIP
29 #define SKIP(i,j,k) (womask!=NULL && womask[(i)+(j)*nxnew+(k)*nxynew]==0)
30
31 /*--------------------------------------------------------------------------*/
32
33 static int zout = 1 ;
mri_warp3D_zerout(int zzzz)34 void mri_warp3D_zerout( int zzzz ){ zout = zzzz ; } /** outside = 0? **/
35
36 /*--------------------------------------------------------------------------*/
37 /* Macros for interpolation and access to data arrays. */
38
39 #undef FAR
40 #define FAR(i,j,k) far[(i)+(j)*nx+(k)*nxy] /* input point */
41 #undef NAR
42 #define NAR(i,j,k) nar[(i)+(j)*nxnew+(k)*nxynew] /* output point */
43
44 /** clip input mm to range [0..nn] (for indexing) **/
45
46 #undef CLIP
47 #define CLIP(mm,nn) if(mm < 0)mm=0; else if(mm > nn)mm=nn
48
49 /* define Lagrange cubic interpolation polynomials */
50
51 #undef P_M1
52 #undef P_00
53 #undef P_P1
54 #undef P_P2
55 #undef P_FACTOR
56 #define P_M1(x) (-(x)*((x)-1)*((x)-2))
57 #define P_00(x) (3*((x)+1)*((x)-1)*((x)-2))
58 #define P_P1(x) (-3*(x)*((x)+1)*((x)-2))
59 #define P_P2(x) ((x)*((x)+1)*((x)-1))
60 #define P_FACTOR 4.62962963e-3 /* 1/216 = final scaling factor */
61
62 /*---------------------------------------------------------------------------*/
63 /*! Transform a 3D image geometrically, using cubic interpolation.
64 See comments for mri_warp3D() for details about the parameters.
65 -----------------------------------------------------------------------------*/
66
mri_warp3D_cubic(MRI_IMAGE * im,int nxnew,int nynew,int nznew,void wf (float,float,float,float *,float *,float *))67 MRI_IMAGE *mri_warp3D_cubic(
68 MRI_IMAGE *im, int nxnew, int nynew, int nznew,
69 void wf(float,float,float,float *,float *,float *) )
70 {
71 MRI_IMAGE *imfl , *newImg ;
72 float *ffar ;
73 int nx,ny,nz,nxy, nxynew , nxyznew , pp ;
74 float bot,top ;
75
76 ENTRY("mri_warp3D_cubic") ;
77
78 /*--- sanity check inputs ---*/
79
80 if( im == NULL || wf == NULL ) RETURN(NULL) ;
81
82 /*-- dimensional analysis --*/
83
84 nx = im->nx ; /* input image dimensions */
85 ny = im->ny ;
86 nz = im->nz ; nxy = nx*ny ;
87
88 nxnew = (nxnew > 0) ? nxnew : nx ; /* output image dimensions */
89 nynew = (nynew > 0) ? nynew : ny ;
90 nznew = (nznew > 0) ? nznew : nz ;
91 nxynew = nxnew*nynew ; nxyznew = nxynew * nznew ;
92
93 /*----- Allow for different input image types, by breaking them into
94 components, doing each one separately, and then reassembling.
95 This is needed since we only warp float-valued images below. -----*/
96
97 switch( im->kind ){
98
99 case MRI_float: /* use input directly */
100 imfl = im ; break ;
101
102 default:{ /* floatize input */
103 imfl = mri_to_float(im) ;
104 newImg = mri_warp3D_cubic( imfl , nxnew,nynew,nznew , wf ) ;
105 mri_free(imfl) ;
106 imfl = mri_to_mri(im->kind,newImg) ;
107 if( imfl != NULL ){ mri_free(newImg); newImg = imfl; }
108 RETURN(newImg) ;
109 }
110
111 case MRI_rgb:{ /* break into 3 pieces */
112 MRI_IMARR *imar = mri_rgb_to_3float(im) ;
113 MRI_IMAGE *rim,*gim,*bim ;
114 rim = mri_warp3D_cubic( IMARR_SUBIM(imar,0), nxnew,nynew,nznew, wf ) ;
115 mri_free( IMARR_SUBIM(imar,0) ) ;
116 gim = mri_warp3D_cubic( IMARR_SUBIM(imar,1), nxnew,nynew,nznew, wf ) ;
117 mri_free( IMARR_SUBIM(imar,1) ) ;
118 bim = mri_warp3D_cubic( IMARR_SUBIM(imar,2), nxnew,nynew,nznew, wf ) ;
119 mri_free( IMARR_SUBIM(imar,2) ) ;
120 FREE_IMARR(imar) ;
121 newImg = mri_3to_rgb( rim,gim,bim ) ;
122 mri_free(rim); mri_free(gim); mri_free(bim); RETURN(newImg);
123 }
124
125 case MRI_complex:{ /* break into 2 pieces */
126 MRI_IMARR *imar = mri_complex_to_pair(im) ;
127 MRI_IMAGE *rim, *iim ;
128 rim = mri_warp3D_cubic( IMARR_SUBIM(imar,0), nxnew,nynew,nznew, wf ) ;
129 mri_free( IMARR_SUBIM(imar,0) ) ;
130 iim = mri_warp3D_cubic( IMARR_SUBIM(imar,1), nxnew,nynew,nznew, wf ) ;
131 mri_free( IMARR_SUBIM(imar,1) ) ;
132 FREE_IMARR(imar) ;
133 newImg = mri_pair_to_complex( rim , iim ) ;
134 mri_free(rim); mri_free(iim); RETURN(newImg);
135 }
136
137 } /* end of special cases of input datum */
138
139 /*----- at this point, imfl is in float format -----*/
140
141 ffar = MRI_FLOAT_PTR( imfl ) ; /* input image data */
142
143 newImg = mri_new_vol( nxnew,nynew,nznew, MRI_float ) ; /* make output image */
144
145 bot = top = ffar[0] ; /* find input data range */
146 for( pp=1 ; pp < imfl->nvox ; pp++ ){
147 if( ffar[pp] > top ) top = ffar[pp] ;
148 else if( ffar[pp] < bot ) bot = ffar[pp] ;
149 }
150
151 /*** loop over output points and warp to them ***/
152
153 AFNI_OMP_START ;
154 #pragma omp parallel if( nxyznew > 99999 )
155 {
156 float *far , *nar ;
157 float nxh,nyh,nzh , val ;
158 int ix_m1,ix_00,ix_p1,ix_p2 ; /* interpolation indices */
159 int jy_m1,jy_00,jy_p1,jy_p2 ; /* (input image) */
160 int kz_m1,kz_00,kz_p1,kz_p2 ;
161 float wt_m1,wt_00,wt_p1,wt_p2 ; /* interpolation weights */
162 float f_jm1_km1, f_j00_km1, f_jp1_km1, f_jp2_km1, /* interpolants */
163 f_jm1_k00, f_j00_k00, f_jp1_k00, f_jp2_k00,
164 f_jm1_kp1, f_j00_kp1, f_jp1_kp1, f_jp2_kp1,
165 f_jm1_kp2, f_j00_kp2, f_jp1_kp2, f_jp2_kp2,
166 f_km1 , f_k00 , f_kp1 , f_kp2 ;
167 int nx1,ny1,nz1 , ii,jj,kk , ix,jy,kz , qq , do_zout=zout ;
168 float xpr,ypr,zpr, xx,yy,zz, fx,fy,fz ;
169
170 far = ffar ; /* input image data */
171 nar = MRI_FLOAT_PTR( newImg ) ; /* output image data */
172
173 nxh = nx-0.5 ; nyh = ny-0.5 ; nzh = nz-0.5 ;
174 nx1 = nx-1 ; ny1 = ny-1 ; nz1 = nz-1 ;
175
176 #pragma omp for
177 for( qq=0 ; qq < nxyznew ; qq++ ){ /* voxel loop */
178 xpr = ii = qq % nxnew ; /* break 1D index qq into 3D index (ii,jj,kk) */
179 zpr = kk = qq / nxynew ;
180 ypr = jj = (qq-kk*nxynew) / nxnew ;
181 if( SKIP(ii,jj,kk) ) continue ; /* 19 Nov 2004 */
182 wf( xpr,ypr,zpr , &xx,&yy,&zz ) ; /* get xx,yy,zz in original image */
183
184 if( do_zout &&
185 (xx < -0.5 || xx > nxh ||
186 yy < -0.5 || yy > nyh || zz < -0.5 || zz > nzh ) ){ nar[qq] = 0.0f; continue; }
187
188 ix = floor(xx) ; fx = xx - ix ; /* integer and */
189 jy = floor(yy) ; fy = yy - jy ; /* fractional coords */
190 kz = floor(zz) ; fz = zz - kz ;
191
192 ix_m1 = ix-1 ; ix_00 = ix ; ix_p1 = ix+1 ; ix_p2 = ix+2 ;
193 CLIP(ix_m1,nx1) ; CLIP(ix_00,nx1) ; CLIP(ix_p1,nx1) ; CLIP(ix_p2,nx1) ;
194
195 jy_m1 = jy-1 ; jy_00 = jy ; jy_p1 = jy+1 ; jy_p2 = jy+2 ;
196 CLIP(jy_m1,ny1) ; CLIP(jy_00,ny1) ; CLIP(jy_p1,ny1) ; CLIP(jy_p2,ny1) ;
197
198 kz_m1 = kz-1 ; kz_00 = kz ; kz_p1 = kz+1 ; kz_p2 = kz+2 ;
199 CLIP(kz_m1,nz1) ; CLIP(kz_00,nz1) ; CLIP(kz_p1,nz1) ; CLIP(kz_p2,nz1) ;
200
201 wt_m1 = P_M1(fx) ; wt_00 = P_00(fx) ; /* interpolation weights */
202 wt_p1 = P_P1(fx) ; wt_p2 = P_P2(fx) ;
203
204 #undef XINT
205 #define XINT(j,k) wt_m1*FAR(ix_m1,j,k)+wt_00*FAR(ix_00,j,k) \
206 +wt_p1*FAR(ix_p1,j,k)+wt_p2*FAR(ix_p2,j,k)
207
208 /* interpolate to location ix+fx at each jy,kz level */
209
210 f_jm1_km1 = XINT(jy_m1,kz_m1) ; f_j00_km1 = XINT(jy_00,kz_m1) ;
211 f_jp1_km1 = XINT(jy_p1,kz_m1) ; f_jp2_km1 = XINT(jy_p2,kz_m1) ;
212 f_jm1_k00 = XINT(jy_m1,kz_00) ; f_j00_k00 = XINT(jy_00,kz_00) ;
213 f_jp1_k00 = XINT(jy_p1,kz_00) ; f_jp2_k00 = XINT(jy_p2,kz_00) ;
214 f_jm1_kp1 = XINT(jy_m1,kz_p1) ; f_j00_kp1 = XINT(jy_00,kz_p1) ;
215 f_jp1_kp1 = XINT(jy_p1,kz_p1) ; f_jp2_kp1 = XINT(jy_p2,kz_p1) ;
216 f_jm1_kp2 = XINT(jy_m1,kz_p2) ; f_j00_kp2 = XINT(jy_00,kz_p2) ;
217 f_jp1_kp2 = XINT(jy_p1,kz_p2) ; f_jp2_kp2 = XINT(jy_p2,kz_p2) ;
218
219 /* interpolate to jy+fy at each kz level */
220
221 wt_m1 = P_M1(fy) ; wt_00 = P_00(fy) ;
222 wt_p1 = P_P1(fy) ; wt_p2 = P_P2(fy) ;
223
224 f_km1 = wt_m1 * f_jm1_km1 + wt_00 * f_j00_km1
225 + wt_p1 * f_jp1_km1 + wt_p2 * f_jp2_km1 ;
226
227 f_k00 = wt_m1 * f_jm1_k00 + wt_00 * f_j00_k00
228 + wt_p1 * f_jp1_k00 + wt_p2 * f_jp2_k00 ;
229
230 f_kp1 = wt_m1 * f_jm1_kp1 + wt_00 * f_j00_kp1
231 + wt_p1 * f_jp1_kp1 + wt_p2 * f_jp2_kp1 ;
232
233 f_kp2 = wt_m1 * f_jm1_kp2 + wt_00 * f_j00_kp2
234 + wt_p1 * f_jp1_kp2 + wt_p2 * f_jp2_kp2 ;
235
236 /* interpolate to kz+fz to get output */
237
238 wt_m1 = P_M1(fz) ; wt_00 = P_00(fz) ;
239 wt_p1 = P_P1(fz) ; wt_p2 = P_P2(fz) ;
240
241 val = P_FACTOR * ( wt_m1 * f_km1 + wt_00 * f_k00
242 + wt_p1 * f_kp1 + wt_p2 * f_kp2 ) ;
243
244 if( val > top ) val = top ; /* clip to input data range */
245 else if( val < bot ) val = bot ;
246
247 nar[qq] = val ;
248 }
249
250 } /* end OpenMP */
251 AFNI_OMP_END ;
252
253 /*** cleanup and return ***/
254
255 if( im != imfl ) mri_free(imfl) ; /* throw away unneeded workspace */
256 RETURN(newImg);
257 }
258
259 /*---------------------------------------------------------------------------*/
260 /*! Transform a 3D image geometrically, using linear interpolation.
261 See comments for mri_warp3D() for details about the parameters.
262 -----------------------------------------------------------------------------*/
263
mri_warp3D_linear(MRI_IMAGE * im,int nxnew,int nynew,int nznew,void wf (float,float,float,float *,float *,float *))264 MRI_IMAGE *mri_warp3D_linear(
265 MRI_IMAGE *im, int nxnew, int nynew, int nznew,
266 void wf(float,float,float,float *,float *,float *) )
267 {
268 MRI_IMAGE *imfl , *newImg ;
269 int nx,ny,nz,nxy , nxynew,nxyznew ;
270
271 #if 0
272 int nzset ; float zzsum ;
273 #endif
274
275 ENTRY("mri_warp3D_linear") ;
276
277 /*--- sanity check inputs ---*/
278
279 if( im == NULL || wf == NULL ) RETURN(NULL) ;
280
281 /*-- dimensional analysis --*/
282
283 nx = im->nx ; /* input image dimensions */
284 ny = im->ny ;
285 nz = im->nz ; nxy = nx*ny ;
286
287 #if 0
288 fprintf(stderr,"mri_warp3D_linear: nx=%d ny=%d nz=%d\n",nx,ny,nz) ;
289 nzset = 0 ; zzsum = 0.0 ;
290 #endif
291
292 nxnew = (nxnew > 0) ? nxnew : nx ; /* output image dimensions */
293 nynew = (nynew > 0) ? nynew : ny ;
294 nznew = (nznew > 0) ? nznew : nz ;
295 nxynew = nxnew*nynew ; nxyznew = nxynew*nznew ;
296
297 /*----- allow for different input image types, by breaking them into
298 components, doing each one separately, and then reassembling -----*/
299
300 switch( im->kind ){
301
302 case MRI_float: /* use input directly */
303 imfl = im ; break ;
304
305 default:{ /* floatize-input */
306 imfl = mri_to_float(im) ;
307 newImg = mri_warp3D_linear( imfl , nxnew,nynew,nznew , wf ) ;
308 mri_free(imfl) ;
309 imfl = mri_to_mri(im->kind,newImg) ;
310 if( imfl != NULL ){ mri_free(newImg); newImg = imfl; }
311 RETURN(newImg) ;
312 }
313
314 case MRI_rgb:{
315 MRI_IMARR *imar = mri_rgb_to_3float(im) ;
316 MRI_IMAGE *rim,*gim,*bim ;
317 rim = mri_warp3D_linear( IMARR_SUBIM(imar,0), nxnew,nynew,nznew, wf ) ;
318 mri_free( IMARR_SUBIM(imar,0) ) ;
319 gim = mri_warp3D_linear( IMARR_SUBIM(imar,1), nxnew,nynew,nznew, wf ) ;
320 mri_free( IMARR_SUBIM(imar,1) ) ;
321 bim = mri_warp3D_linear( IMARR_SUBIM(imar,2), nxnew,nynew,nznew, wf ) ;
322 mri_free( IMARR_SUBIM(imar,2) ) ;
323 FREE_IMARR(imar) ;
324 newImg = mri_3to_rgb( rim,gim,bim ) ;
325 mri_free(rim); mri_free(gim); mri_free(bim); RETURN(newImg);
326 }
327
328 case MRI_complex:{
329 MRI_IMARR *imar = mri_complex_to_pair(im) ;
330 MRI_IMAGE *rim, *iim ;
331 rim = mri_warp3D_linear( IMARR_SUBIM(imar,0), nxnew,nynew,nznew, wf ) ;
332 mri_free( IMARR_SUBIM(imar,0) ) ;
333 iim = mri_warp3D_linear( IMARR_SUBIM(imar,1), nxnew,nynew,nznew, wf ) ;
334 mri_free( IMARR_SUBIM(imar,1) ) ;
335 FREE_IMARR(imar) ;
336 newImg = mri_pair_to_complex( rim , iim ) ;
337 mri_free(rim); mri_free(iim); RETURN(newImg);
338 }
339
340 } /* end of special cases of input datum */
341
342 /*----- at this point, imfl is in float format -----*/
343
344 newImg = mri_new_vol( nxnew,nynew,nznew, MRI_float ) ; /* make output image */
345
346 AFNI_OMP_START ;
347 #pragma omp parallel if( nxyznew > 99999 )
348 {
349 float *far , *nar ;
350 float nxh,nyh,nzh ;
351 float xpr,ypr,zpr, xx,yy,zz, fx,fy,fz ;
352 int ii,jj,kk , qq , ix,jy,kz , nx1,ny1,nz1 ;
353 int ix_00,ix_p1 ; /* interpolation indices */
354 int jy_00,jy_p1 ; /* (input image) */
355 int kz_00,kz_p1 ;
356 float wt_00,wt_p1 ; /* interpolation weights */
357 float f_j00_k00, f_jp1_k00, f_j00_kp1, f_jp1_kp1, f_k00, f_kp1 ;
358 int do_zout = zout ;
359
360 far = MRI_FLOAT_PTR( imfl ) ; /* input image data */
361 nar = MRI_FLOAT_PTR( newImg ) ; /* output image data */
362
363 /*** loop over output points and warp to them ***/
364
365 nxh = nx-0.5 ; nyh = ny-0.5 ; nzh = nz-0.5 ;
366 nx1 = nx-1 ; ny1 = ny-1 ; nz1 = nz-1 ;
367
368 #pragma omp for
369 for( qq=0 ; qq < nxyznew ; qq++ ){ /* voxel loop */
370 xpr = ii = qq % nxnew ; /* break 1D index qq into 3D index (ii,jj,kk) */
371 zpr = kk = qq / nxynew ;
372 ypr = jj = (qq-kk*nxynew) / nxnew ;
373 if( SKIP(ii,jj,kk) ) continue ; /* 19 Nov 2004 */
374 wf( xpr,ypr,zpr , &xx,&yy,&zz ) ; /* get xx,yy,zz in original image */
375
376 if( do_zout &&
377 (xx < -0.5 || xx > nxh ||
378 yy < -0.5 || yy > nyh || zz < -0.5 || zz > nzh ) ){ nar[qq] = 0.0f ; continue ; }
379
380 ix = floor(xx) ; fx = xx - ix ; /* integer and */
381 jy = floor(yy) ; fy = yy - jy ; /* fractional coords */
382 kz = floor(zz) ; fz = zz - kz ; /* of warped point */
383
384 ix_00 = ix ; ix_p1 = ix+1 ;
385 CLIP(ix_00,nx1) ; CLIP(ix_p1,nx1) ;
386
387 jy_00 = jy ; jy_p1 = jy+1 ;
388 CLIP(jy_00,ny1) ; CLIP(jy_p1,ny1) ;
389
390 kz_00 = kz ; kz_p1 = kz+1 ;
391 CLIP(kz_00,nz1) ; CLIP(kz_p1,nz1) ;
392
393 wt_00 = 1.0-fx ; wt_p1 = fx ;
394
395 #undef XINT
396 #define XINT(j,k) wt_00*FAR(ix_00,j,k)+wt_p1*FAR(ix_p1,j,k)
397
398 /* interpolate to location ix+fx at each jy,kz level */
399
400 f_j00_k00 = XINT(jy_00,kz_00) ; f_jp1_k00 = XINT(jy_p1,kz_00) ;
401 f_j00_kp1 = XINT(jy_00,kz_p1) ; f_jp1_kp1 = XINT(jy_p1,kz_p1) ;
402
403 /* interpolate to jy+fy at each kz level */
404
405 wt_00 = 1.0-fy ; wt_p1 = fy ;
406
407 f_k00 = wt_00 * f_j00_k00 + wt_p1 * f_jp1_k00 ;
408
409 f_kp1 = wt_00 * f_j00_kp1 + wt_p1 * f_jp1_kp1 ;
410
411 /* interpolate to kz+fz to get output */
412
413 nar[qq] = (1.0-fz) * f_k00 + fz * f_kp1 ;
414
415 } /* end of voxel loop */
416
417 } /* end OpenMP */
418 AFNI_OMP_END ;
419
420 /*** cleanup and return ***/
421
422 if( im != imfl ) mri_free(imfl) ; /* throw away unneeded workspace */
423 RETURN(newImg);
424 }
425
426 /*---------------------------------------------------------------------------*/
427 /*! Transform a 3D image geometrically, using NN 'interpolation'.
428 See comments for mri_warp3D() for details about the parameters.
429 -----------------------------------------------------------------------------*/
430
mri_warp3D_NN(MRI_IMAGE * im,int nxnew,int nynew,int nznew,void wf (float,float,float,float *,float *,float *))431 MRI_IMAGE *mri_warp3D_NN(
432 MRI_IMAGE *im, int nxnew, int nynew, int nznew,
433 void wf(float,float,float,float *,float *,float *) )
434 {
435 MRI_IMAGE *imfl , *newImg ;
436 float *far , *nar ;
437 float xpr,ypr,zpr, xx,yy,zz, fx,fy,fz ;
438 int ii,jj,kk, nx,ny,nz,nxy, nx1,ny1,nz1, ix,jy,kz, nxynew ;
439 float nxh,nyh,nzh ;
440
441 ENTRY("mri_warp3D_NN") ;
442
443 /*--- sanity check inputs ---*/
444
445 if( im == NULL || wf == NULL ) RETURN(NULL) ;
446
447 /*-- dimensional analysis --*/
448
449 nx = im->nx ; nx1 = nx-1 ; /* input image dimensions */
450 ny = im->ny ; ny1 = ny-1 ;
451 nz = im->nz ; nz1 = nz-1 ; nxy = nx*ny ;
452
453 nxnew = (nxnew > 0) ? nxnew : nx ; /* output image dimensions */
454 nynew = (nynew > 0) ? nynew : ny ;
455 nznew = (nznew > 0) ? nznew : nz ;
456 nxynew = nxnew*nynew ;
457
458 /*----- allow for different input image types, by breaking them into
459 components, doing each one separately, and then reassembling -----*/
460
461 switch( im->kind ){
462
463 case MRI_float: /* use input directly */
464 imfl = im ; break ;
465
466 default:{ /* floatize-input */
467 imfl = mri_to_float(im) ;
468 newImg = mri_warp3D_NN( imfl , nxnew,nynew,nznew , wf ) ;
469 mri_free(imfl) ;
470 imfl = mri_to_mri(im->kind,newImg) ;
471 if( imfl != NULL ){ mri_free(newImg); newImg = imfl; }
472 RETURN(newImg) ;
473 }
474
475 case MRI_rgb:{
476 MRI_IMARR *imar = mri_rgb_to_3float(im) ;
477 MRI_IMAGE *rim,*gim,*bim ;
478 rim = mri_warp3D_NN( IMARR_SUBIM(imar,0), nxnew,nynew,nznew, wf ) ;
479 mri_free( IMARR_SUBIM(imar,0) ) ;
480 gim = mri_warp3D_NN( IMARR_SUBIM(imar,1), nxnew,nynew,nznew, wf ) ;
481 mri_free( IMARR_SUBIM(imar,1) ) ;
482 bim = mri_warp3D_NN( IMARR_SUBIM(imar,2), nxnew,nynew,nznew, wf ) ;
483 mri_free( IMARR_SUBIM(imar,2) ) ;
484 FREE_IMARR(imar) ;
485 newImg = mri_3to_rgb( rim,gim,bim ) ;
486 mri_free(rim); mri_free(gim); mri_free(bim); RETURN(newImg);
487 }
488
489 case MRI_complex:{
490 MRI_IMARR *imar = mri_complex_to_pair(im) ;
491 MRI_IMAGE *rim, *iim ;
492 rim = mri_warp3D_NN( IMARR_SUBIM(imar,0), nxnew,nynew,nznew, wf ) ;
493 mri_free( IMARR_SUBIM(imar,0) ) ;
494 iim = mri_warp3D_NN( IMARR_SUBIM(imar,1), nxnew,nynew,nznew, wf ) ;
495 mri_free( IMARR_SUBIM(imar,1) ) ;
496 FREE_IMARR(imar) ;
497 newImg = mri_pair_to_complex( rim , iim ) ;
498 mri_free(rim); mri_free(iim); RETURN(newImg);
499 }
500
501 } /* end of special cases of input datum */
502
503 /*----- at this point, imfl is in float format -----*/
504
505 far = MRI_FLOAT_PTR( imfl ) ; /* input image data */
506
507 newImg = mri_new_vol( nxnew,nynew,nznew, MRI_float ) ; /* make output image */
508 nar = MRI_FLOAT_PTR( newImg ) ; /* output image data */
509
510 /*** loop over output points and warp to them ***/
511
512 nxh = nx-0.5 ; nyh = ny-0.5 ; nzh = nz-0.5 ;
513 for( kk=0 ; kk < nznew ; kk++ ){
514 zpr = kk ;
515 for( jj=0 ; jj < nynew ; jj++ ){
516 ypr = jj ;
517 for( ii=0 ; ii < nxnew ; ii++ ){
518 xpr = ii ;
519 if( SKIP(ii,jj,kk) ) continue ; /* 19 Nov 2004 */
520 wf( xpr,ypr,zpr , &xx,&yy,&zz ) ; /* get xx,yy,zz in original image */
521
522 if( zout &&
523 (xx < -0.5 || xx > nxh ||
524 yy < -0.5 || yy > nyh || zz < -0.5 || zz > nzh ) ){
525
526 NAR(ii,jj,kk) = 0.0 ; continue ;
527 }
528
529 ix = (int)rint(xx) ; jy = (int)rint(yy) ; kz = (int)rint(zz) ;
530 CLIP(ix,nx1) ; CLIP(jy,ny1) ; CLIP(kz,nz1) ;
531
532 NAR(ii,jj,kk) = FAR(ix,jy,kz) ;
533 }
534 }
535 }
536
537 /*** cleanup and return ***/
538
539 if( im != imfl ) mri_free(imfl) ; /* throw away unneeded workspace */
540 RETURN(newImg);
541 }
542
543 /* define quintic interpolation polynomials (Lagrange) */
544
545 #undef Q_M2
546 #undef Q_M1
547 #undef Q_00
548 #undef Q_P1
549 #undef Q_P2
550 #undef Q_P3
551 #define Q_M2(x) (x*(x*x-1.0)*(2.0-x)*(x-3.0)*0.008333333)
552 #define Q_M1(x) (x*(x*x-4.0)*(x-1.0)*(x-3.0)*0.041666667)
553 #define Q_00(x) ((x*x-4.0)*(x*x-1.0)*(3.0-x)*0.083333333)
554 #define Q_P1(x) (x*(x*x-4.0)*(x+1.0)*(x-3.0)*0.083333333)
555 #define Q_P2(x) (x*(x*x-1.0)*(x+2.0)*(3.0-x)*0.041666667)
556 #define Q_P3(x) (x*(x*x-1.0)*(x*x-4.0)*0.008333333)
557
558 /*---------------------------------------------------------------------------*/
559 /*! Transform a 3D image geometrically, using quintic interpolation.
560 See comments for mri_warp3D() for details about the parameters.
561 - RWCox - 06 Aug 2003
562 -----------------------------------------------------------------------------*/
563
mri_warp3D_quintic(MRI_IMAGE * im,int nxnew,int nynew,int nznew,void wf (float,float,float,float *,float *,float *))564 MRI_IMAGE *mri_warp3D_quintic(
565 MRI_IMAGE *im, int nxnew, int nynew, int nznew,
566 void wf(float,float,float,float *,float *,float *) )
567 {
568 MRI_IMAGE *imfl , *newImg ;
569 float *far , *nar ;
570 float xpr,ypr,zpr, xx,yy,zz, fx,fy,fz ;
571 int ii,jj,kk, nx,ny,nz,nxy, nx1,ny1,nz1, ix,jy,kz, nxynew ;
572 float bot,top,val ;
573 float nxh,nyh,nzh ;
574
575 int ix_m2,ix_m1,ix_00,ix_p1,ix_p2,ix_p3 ; /* interpolation indices */
576 int jy_m2,jy_m1,jy_00,jy_p1,jy_p2,jy_p3 ; /* (input image) */
577 int kz_m2,kz_m1,kz_00,kz_p1,kz_p2,kz_p3 ;
578
579 float wt_m2,wt_m1,wt_00,wt_p1,wt_p2,wt_p3 ; /* interpolation weights */
580
581 float f_jm2_km2, f_jm1_km2, f_j00_km2, f_jp1_km2, f_jp2_km2, f_jp3_km2,
582 f_jm2_km1, f_jm1_km1, f_j00_km1, f_jp1_km1, f_jp2_km1, f_jp3_km1,
583 f_jm2_k00, f_jm1_k00, f_j00_k00, f_jp1_k00, f_jp2_k00, f_jp3_k00,
584 f_jm2_kp1, f_jm1_kp1, f_j00_kp1, f_jp1_kp1, f_jp2_kp1, f_jp3_kp1,
585 f_jm2_kp2, f_jm1_kp2, f_j00_kp2, f_jp1_kp2, f_jp2_kp2, f_jp3_kp2,
586 f_jm2_kp3, f_jm1_kp3, f_j00_kp3, f_jp1_kp3, f_jp2_kp3, f_jp3_kp3,
587 f_km2 , f_km1 , f_k00 , f_kp1 , f_kp2 , f_kp3 ;
588
589 ENTRY("mri_warp3D_quintic") ;
590
591 /*--- sanity check inputs ---*/
592
593 if( im == NULL || wf == NULL ) RETURN(NULL) ;
594
595 /*-- dimensional analysis --*/
596
597 nx = im->nx ; nx1 = nx-1 ; /* input image dimensions */
598 ny = im->ny ; ny1 = ny-1 ;
599 nz = im->nz ; nz1 = nz-1 ; nxy = nx*ny ;
600
601 nxnew = (nxnew > 0) ? nxnew : nx ; /* output image dimensions */
602 nynew = (nynew > 0) ? nynew : ny ;
603 nznew = (nznew > 0) ? nznew : nz ;
604 nxynew = nxnew*nynew ;
605
606 /*----- Allow for different input image types, by breaking them into
607 components, doing each one separately, and then reassembling.
608 This is needed since we only warp float-valued images below. -----*/
609
610 switch( im->kind ){
611
612 case MRI_float: /* use input directly */
613 imfl = im ; break ;
614
615 default:{ /* floatize input */
616 imfl = mri_to_float(im) ;
617 newImg = mri_warp3D_quintic( imfl , nxnew,nynew,nznew , wf ) ;
618 mri_free(imfl) ;
619 imfl = mri_to_mri(im->kind,newImg) ;
620 if( imfl != NULL ){ mri_free(newImg); newImg = imfl; }
621 RETURN(newImg) ;
622 }
623
624 case MRI_rgb:{ /* break into 3 pieces */
625 MRI_IMARR *imar = mri_rgb_to_3float(im) ;
626 MRI_IMAGE *rim,*gim,*bim ;
627 rim = mri_warp3D_quintic( IMARR_SUBIM(imar,0), nxnew,nynew,nznew, wf ) ;
628 mri_free( IMARR_SUBIM(imar,0) ) ;
629 gim = mri_warp3D_quintic( IMARR_SUBIM(imar,1), nxnew,nynew,nznew, wf ) ;
630 mri_free( IMARR_SUBIM(imar,1) ) ;
631 bim = mri_warp3D_quintic( IMARR_SUBIM(imar,2), nxnew,nynew,nznew, wf ) ;
632 mri_free( IMARR_SUBIM(imar,2) ) ;
633 FREE_IMARR(imar) ;
634 newImg = mri_3to_rgb( rim,gim,bim ) ;
635 mri_free(rim); mri_free(gim); mri_free(bim); RETURN(newImg);
636 }
637
638 case MRI_complex:{ /* break into 2 pieces */
639 MRI_IMARR *imar = mri_complex_to_pair(im) ;
640 MRI_IMAGE *rim, *iim ;
641 rim = mri_warp3D_quintic( IMARR_SUBIM(imar,0), nxnew,nynew,nznew, wf ) ;
642 mri_free( IMARR_SUBIM(imar,0) ) ;
643 iim = mri_warp3D_quintic( IMARR_SUBIM(imar,1), nxnew,nynew,nznew, wf ) ;
644 mri_free( IMARR_SUBIM(imar,1) ) ;
645 FREE_IMARR(imar) ;
646 newImg = mri_pair_to_complex( rim , iim ) ;
647 mri_free(rim); mri_free(iim); RETURN(newImg);
648 }
649
650 } /* end of special cases of input datum */
651
652 /*----- at this point, imfl is in float format -----*/
653
654 far = MRI_FLOAT_PTR( imfl ) ; /* input image data */
655
656 newImg = mri_new_vol( nxnew,nynew,nznew, MRI_float ) ; /* make output image */
657 nar = MRI_FLOAT_PTR( newImg ) ; /* output image data */
658
659 bot = top = far[0] ; /* find input data range */
660 for( ii=1 ; ii < imfl->nvox ; ii++ ){
661 if( far[ii] > top ) top = far[ii] ;
662 else if( far[ii] < bot ) bot = far[ii] ;
663 }
664
665 /*** loop over output points and warp to them ***/
666
667 nxh = nx-0.5 ; nyh = ny-0.5 ; nzh = nz-0.5 ;
668 for( kk=0 ; kk < nznew ; kk++ ){
669 zpr = kk ;
670 for( jj=0 ; jj < nynew ; jj++ ){
671 ypr = jj ;
672 for( ii=0 ; ii < nxnew ; ii++ ){
673 xpr = ii ;
674 if( SKIP(ii,jj,kk) ) continue ; /* 19 Nov 2004 */
675 wf( xpr,ypr,zpr , &xx,&yy,&zz ) ; /* get xx,yy,zz in original image */
676
677 if( zout &&
678 (xx < -0.5 || xx > nxh ||
679 yy < -0.5 || yy > nyh || zz < -0.5 || zz > nzh ) ){
680
681 NAR(ii,jj,kk) = 0.0 ; continue ;
682 }
683
684 ix = floor(xx) ; fx = xx - ix ; /* integer and */
685 jy = floor(yy) ; fy = yy - jy ; /* fractional coords */
686 kz = floor(zz) ; fz = zz - kz ; /* in input image */
687
688 /* compute indexes from which to interpolate (-2,-1,0,+1,+2,+3),
689 but clipped to lie inside input image volume */
690
691 ix_m1 = ix-1 ; ix_00 = ix ; ix_p1 = ix+1 ; ix_p2 = ix+2 ;
692 CLIP(ix_m1,nx1) ; CLIP(ix_00,nx1) ; CLIP(ix_p1,nx1) ; CLIP(ix_p2,nx1) ;
693 ix_m2 = ix-2 ; ix_p3 = ix+3 ;
694 CLIP(ix_m2,nx1) ; CLIP(ix_p3,nx1) ;
695
696 jy_m1 = jy-1 ; jy_00 = jy ; jy_p1 = jy+1 ; jy_p2 = jy+2 ;
697 CLIP(jy_m1,ny1) ; CLIP(jy_00,ny1) ; CLIP(jy_p1,ny1) ; CLIP(jy_p2,ny1) ;
698 jy_m2 = jy-2 ; jy_p3 = jy+3 ;
699 CLIP(jy_m2,ny1) ; CLIP(jy_p3,ny1) ;
700
701 kz_m1 = kz-1 ; kz_00 = kz ; kz_p1 = kz+1 ; kz_p2 = kz+2 ;
702 CLIP(kz_m1,nz1) ; CLIP(kz_00,nz1) ; CLIP(kz_p1,nz1) ; CLIP(kz_p2,nz1) ;
703 kz_m2 = kz-2 ; kz_p3 = kz+3 ;
704 CLIP(kz_m2,nz1) ; CLIP(kz_p3,nz1) ;
705
706 wt_m1 = Q_M1(fx) ; wt_00 = Q_00(fx) ; /* interpolation weights */
707 wt_p1 = Q_P1(fx) ; wt_p2 = Q_P2(fx) ; /* in x-direction */
708 wt_m2 = Q_M2(fx) ; wt_p3 = Q_P3(fx) ;
709
710 #undef XINT
711 #define XINT(j,k) wt_m2*FAR(ix_m2,j,k)+wt_m1*FAR(ix_m1,j,k) \
712 +wt_00*FAR(ix_00,j,k)+wt_p1*FAR(ix_p1,j,k) \
713 +wt_p2*FAR(ix_p2,j,k)+wt_p3*FAR(ix_p3,j,k)
714
715 /* interpolate to location ix+fx at each jy,kz level */
716
717 f_jm2_km2 = XINT(jy_m2,kz_m2) ; f_jm1_km2 = XINT(jy_m1,kz_m2) ;
718 f_j00_km2 = XINT(jy_00,kz_m2) ; f_jp1_km2 = XINT(jy_p1,kz_m2) ;
719 f_jp2_km2 = XINT(jy_p2,kz_m2) ; f_jp3_km2 = XINT(jy_p3,kz_m2) ;
720
721 f_jm2_km1 = XINT(jy_m2,kz_m1) ; f_jm1_km1 = XINT(jy_m1,kz_m1) ;
722 f_j00_km1 = XINT(jy_00,kz_m1) ; f_jp1_km1 = XINT(jy_p1,kz_m1) ;
723 f_jp2_km1 = XINT(jy_p2,kz_m1) ; f_jp3_km1 = XINT(jy_p3,kz_m1) ;
724
725 f_jm2_k00 = XINT(jy_m2,kz_00) ; f_jm1_k00 = XINT(jy_m1,kz_00) ;
726 f_j00_k00 = XINT(jy_00,kz_00) ; f_jp1_k00 = XINT(jy_p1,kz_00) ;
727 f_jp2_k00 = XINT(jy_p2,kz_00) ; f_jp3_k00 = XINT(jy_p3,kz_00) ;
728
729 f_jm2_kp1 = XINT(jy_m2,kz_p1) ; f_jm1_kp1 = XINT(jy_m1,kz_p1) ;
730 f_j00_kp1 = XINT(jy_00,kz_p1) ; f_jp1_kp1 = XINT(jy_p1,kz_p1) ;
731 f_jp2_kp1 = XINT(jy_p2,kz_p1) ; f_jp3_kp1 = XINT(jy_p3,kz_p1) ;
732
733 f_jm2_kp2 = XINT(jy_m2,kz_p2) ; f_jm1_kp2 = XINT(jy_m1,kz_p2) ;
734 f_j00_kp2 = XINT(jy_00,kz_p2) ; f_jp1_kp2 = XINT(jy_p1,kz_p2) ;
735 f_jp2_kp2 = XINT(jy_p2,kz_p2) ; f_jp3_kp2 = XINT(jy_p3,kz_p2) ;
736
737 f_jm2_kp3 = XINT(jy_m2,kz_p3) ; f_jm1_kp3 = XINT(jy_m1,kz_p3) ;
738 f_j00_kp3 = XINT(jy_00,kz_p3) ; f_jp1_kp3 = XINT(jy_p1,kz_p3) ;
739 f_jp2_kp3 = XINT(jy_p2,kz_p3) ; f_jp3_kp3 = XINT(jy_p3,kz_p3) ;
740
741 /* interpolate to jy+fy at each kz level */
742
743 wt_m1 = Q_M1(fy) ; wt_00 = Q_00(fy) ; wt_p1 = Q_P1(fy) ;
744 wt_p2 = Q_P2(fy) ; wt_m2 = Q_M2(fy) ; wt_p3 = Q_P3(fy) ;
745
746 f_km2 = wt_m2 * f_jm2_km2 + wt_m1 * f_jm1_km2 + wt_00 * f_j00_km2
747 + wt_p1 * f_jp1_km2 + wt_p2 * f_jp2_km2 + wt_p3 * f_jp3_km2 ;
748
749 f_km1 = wt_m2 * f_jm2_km1 + wt_m1 * f_jm1_km1 + wt_00 * f_j00_km1
750 + wt_p1 * f_jp1_km1 + wt_p2 * f_jp2_km1 + wt_p3 * f_jp3_km1 ;
751
752 f_k00 = wt_m2 * f_jm2_k00 + wt_m1 * f_jm1_k00 + wt_00 * f_j00_k00
753 + wt_p1 * f_jp1_k00 + wt_p2 * f_jp2_k00 + wt_p3 * f_jp3_k00 ;
754
755 f_kp1 = wt_m2 * f_jm2_kp1 + wt_m1 * f_jm1_kp1 + wt_00 * f_j00_kp1
756 + wt_p1 * f_jp1_kp1 + wt_p2 * f_jp2_kp1 + wt_p3 * f_jp3_kp1 ;
757
758 f_kp2 = wt_m2 * f_jm2_kp2 + wt_m1 * f_jm1_kp2 + wt_00 * f_j00_kp2
759 + wt_p1 * f_jp1_kp2 + wt_p2 * f_jp2_kp2 + wt_p3 * f_jp3_kp2 ;
760
761 f_kp3 = wt_m2 * f_jm2_kp3 + wt_m1 * f_jm1_kp3 + wt_00 * f_j00_kp3
762 + wt_p1 * f_jp1_kp3 + wt_p2 * f_jp2_kp3 + wt_p3 * f_jp3_kp3 ;
763
764 /* interpolate to kz+fz to get output */
765
766 wt_m1 = Q_M1(fz) ; wt_00 = Q_00(fz) ; wt_p1 = Q_P1(fz) ;
767 wt_p2 = Q_P2(fz) ; wt_m2 = Q_M2(fz) ; wt_p3 = Q_P3(fz) ;
768
769 val = wt_m2 * f_km2 + wt_m1 * f_km1 + wt_00 * f_k00
770 + wt_p1 * f_kp1 + wt_p2 * f_kp2 + wt_p3 * f_kp3 ;
771
772 if( val > top ) val = top ; /* clip to input data range */
773 else if( val < bot ) val = bot ;
774
775 NAR(ii,jj,kk) = val ;
776 }
777 }
778 }
779
780 /*** cleanup and return ***/
781
782 if( im != imfl ) mri_free(imfl) ; /* throw away unneeded workspace */
783 RETURN(newImg);
784 }
785
786 /*--------------------------------------------------------------------------*/
787 /*! Generic warping function. Calls the specific one for the desired
788 interpolation method, which was set by mri_warp3D_method().
789 - im = input image (3D).
790 - nxnew,nynew,nznew = size of output image grid.
791 - wf( inew,jnew,knew , iold,jold,kold ) is a function that takes
792 as input 3 indices in [0..nxnew-1,0..nynew-1,0..nznew-1] (the output
793 image index space) and returns indices in the input image space.
794 - indices in and out of wf() are floats.
795 - by default, values outside the original image volume are treated
796 as zero; this mode can be turned off using mri_warp3D_zerout().
797 ----------------------------------------------------------------------------*/
798
mri_warp3D(MRI_IMAGE * im,int nxnew,int nynew,int nznew,void wf (float,float,float,float *,float *,float *))799 MRI_IMAGE *mri_warp3D(
800 MRI_IMAGE *im, int nxnew, int nynew, int nznew,
801 void wf(float,float,float,float *,float *,float *) )
802 {
803 switch( wtype ){
804
805 default:
806 case MRI_CUBIC:
807 return mri_warp3D_cubic ( im , nxnew,nynew,nznew , wf ) ;
808
809 case MRI_LINEAR:
810 return mri_warp3D_linear ( im , nxnew,nynew,nznew , wf ) ;
811
812 case MRI_NN:
813 return mri_warp3D_NN ( im , nxnew,nynew,nznew , wf ) ;
814
815 case MRI_QUINTIC: /* 06 Aug 2003 */
816 return mri_warp3D_quintic( im , nxnew,nynew,nznew , wf ) ;
817 }
818 return NULL ; /* unreachable */
819 }
820
821 /*---------------------------------------------------------------------------*/
822
823 static float sx_scale, sy_scale, sz_scale ; /* global scaler data */
824
825 /*! Internal transform function for resize "warp". */
826
w3dMRI_scaler(float a,float b,float c,float * x,float * y,float * z)827 static INLINE void w3dMRI_scaler( float a,float b,float c,
828 float *x,float *y,float *z )
829 {
830 *x = sx_scale * a; *y = sy_scale * b; *z = sz_scale * c;
831 }
832
833 /*--------------------------------------------------------------------------*/
834 /*! Special image warp just for resizing.
835 ----------------------------------------------------------------------------*/
836
mri_warp3D_resize(MRI_IMAGE * im,int nxnew,int nynew,int nznew)837 MRI_IMAGE *mri_warp3D_resize( MRI_IMAGE *im , int nxnew, int nynew, int nznew )
838 {
839 int nx,ny,nz , nnx,nny,nnz ;
840
841 if( im == NULL ) return NULL ;
842 nx = im->nx ; ny = im->ny ; nz = im->nz ;
843 nnx = nxnew ; nny = nynew ; nnz = nznew ;
844
845 if( nnx <= 0 && nny <= 0 && nnz <= 0 ) return NULL ;
846
847 sx_scale = (nnx > 0) ? ((float)nx)/nnx : 0.0 ; /* global variables */
848 sy_scale = (nny > 0) ? ((float)ny)/nny : 0.0 ;
849 sz_scale = (nnz > 0) ? ((float)nz)/nnz : 0.0 ;
850
851 if( nnx <= 0 ){
852 sx_scale = MAX(sy_scale,sz_scale) ;
853 nnx = (int)rint(sx_scale*nx) ;
854 }
855 if( nny <= 0 ){
856 sy_scale = MAX(sx_scale,sz_scale) ;
857 nny = (int)rint(sy_scale*ny) ;
858 }
859 if( nnz <= 0 ){
860 sz_scale = MAX(sx_scale,sy_scale) ;
861 nnz = (int)rint(sz_scale*nz) ;
862 }
863
864 return mri_warp3D( im , nnx,nny,nnz , w3dMRI_scaler ) ;
865 }
866
867 /*---------------------------------------------------------------------------*/
868
869 static float a11_aff,a12_aff,a13_aff,a14_aff ,
870 a21_aff,a22_aff,a23_aff,a24_aff ,
871 a31_aff,a32_aff,a33_aff,a34_aff ;
872
873 /*! Internal transform function for affine "warp". */
874
w3dMRI_affine(float a,float b,float c,float * x,float * y,float * z)875 static INLINE void w3dMRI_affine( float a,float b,float c,
876 float *x,float *y,float *z )
877 {
878 *x = a11_aff*a + a12_aff*b + a13_aff*c + a14_aff ;
879 *y = a21_aff*a + a22_aff*b + a23_aff*c + a24_aff ;
880 *z = a31_aff*a + a32_aff*b + a33_aff*c + a34_aff ;
881 }
882
883 /*--------------------------------------------------------------------------*/
884 /*! Special image warp for affine transformation of indices.
885 ----------------------------------------------------------------------------*/
886
mri_warp3D_affine(MRI_IMAGE * im,THD_vecmat aff)887 MRI_IMAGE * mri_warp3D_affine( MRI_IMAGE *im , THD_vecmat aff )
888 {
889 if( im == NULL ) return NULL ;
890
891 a11_aff = aff.mm.mat[0][0] ; a12_aff = aff.mm.mat[0][1] ;
892 a13_aff = aff.mm.mat[0][2] ; a14_aff = aff.vv.xyz[0] ;
893
894 a21_aff = aff.mm.mat[1][0] ; a22_aff = aff.mm.mat[1][1] ;
895 a23_aff = aff.mm.mat[1][2] ; a24_aff = aff.vv.xyz[1] ;
896
897 a31_aff = aff.mm.mat[2][0] ; a32_aff = aff.mm.mat[2][1] ;
898 a33_aff = aff.mm.mat[2][2] ; a34_aff = aff.vv.xyz[2] ;
899
900 return mri_warp3D( im , 0,0,0 , w3dMRI_scaler ) ;
901 }
902
903 /*===========================================================================*/
904 /*===========================================================================*/
905 /******** Functions for warping a dataset using coordinate transforms, *******/
906 /******** rather than the index transformation of the functions above. *******/
907 /*===========================================================================*/
908 /*===========================================================================*/
909
910 static THD_vecmat ijk_to_dicom_in, ijk_to_dicom_out, dicom_to_ijk_in ;
911
912 static void (*warp_out_to_in)(float ,float ,float ,
913 float *,float *,float * ) ;
914
915 /*----------------------------------------------------------------------------*/
916 /*! This is the function passed to mri_warp3D() for the ijk_out to ijk_in
917 transformation. It does the ijk_out to dicom_out transform here,
918 then calls the coordinate transform function for dicom_out to dicom_in,
919 then does the dicom_in to ijk_in transform here.
920 ------------------------------------------------------------------------------*/
921
w3d_warp_func(float xout,float yout,float zout,float * xin,float * yin,float * zin)922 static INLINE void w3d_warp_func( float xout, float yout, float zout,
923 float *xin, float *yin, float *zin )
924 {
925 THD_fvec3 xxx,yyy ; float xi,yi,zi ;
926
927 LOAD_FVEC3(xxx,xout,yout,zout) ;
928 yyy = VECMAT_VEC( ijk_to_dicom_out , xxx ) ; /* ijk_out to dicom_out */
929 warp_out_to_in( yyy.xyz[0],yyy.xyz[1],yyy.xyz[2], /* dicom_out to dicom_in */
930 &(xxx.xyz[0]) , &(xxx.xyz[1]), &(xxx.xyz[2]) ) ;
931 yyy = VECMAT_VEC( dicom_to_ijk_in , xxx ) ; /* dicom_in to ijk_in */
932 *xin = yyy.xyz[0] ; *yin = yyy.xyz[1] ; *zin = yyy.xyz[2] ;
933 }
934
935 /*--------------------------------------------------------------------------*/
936 /*! - Find the 8 corners of the input dataset (voxel edges, not centers).
937 - Warp each one using the provided wfunc().
938 - Return the min and max (x,y,z) coordinates of these warped points.
939 ----------------------------------------------------------------------------*/
940
warp_corners(THD_3dim_dataset * inset,void wfunc (float,float,float,float *,float *,float *),float * xb,float * xt,float * yb,float * yt,float * zb,float * zt)941 static void warp_corners( THD_3dim_dataset *inset,
942 void wfunc(float,float,float,float *,float *,float *),
943 float *xb , float *xt ,
944 float *yb , float *yt , float *zb , float *zt )
945 {
946 THD_dataxes *daxes = inset->daxes ;
947 THD_fvec3 corner , wcorn ;
948 float nx0 = -0.5 , ny0 = -0.5 , nz0 = -0.5 ;
949 float nx1 = daxes->nxx-0.5, ny1 = daxes->nyy-0.5, nz1 = daxes->nzz-0.5 ;
950 float xx,yy,zz , xbot,ybot,zbot , xtop,ytop,ztop ;
951
952 LOAD_FVEC3( corner , nx0,ny0,nz0 ) ;
953 wcorn = VECMAT_VEC( ijk_to_dicom_in , corner ) ;
954 wfunc( wcorn.xyz[0],wcorn.xyz[1],wcorn.xyz[2] , &xx,&yy,&zz ) ;
955 xbot = xtop = xx ;
956 ybot = ytop = yy ;
957 zbot = ztop = zz ;
958
959 LOAD_FVEC3( corner , nx1,ny0,nz0 ) ;
960 wcorn = VECMAT_VEC( ijk_to_dicom_in , corner ) ;
961 wfunc( wcorn.xyz[0],wcorn.xyz[1],wcorn.xyz[2] , &xx,&yy,&zz ) ;
962 xbot = MIN(xx,xbot); xtop = MAX(xx,xtop);
963 ybot = MIN(yy,ybot); ytop = MAX(yy,ytop);
964 zbot = MIN(zz,zbot); ztop = MAX(zz,ztop);
965
966 LOAD_FVEC3( corner , nx0,ny1,nz0 ) ;
967 wcorn = VECMAT_VEC( ijk_to_dicom_in , corner ) ;
968 wfunc( wcorn.xyz[0],wcorn.xyz[1],wcorn.xyz[2] , &xx,&yy,&zz ) ;
969 xbot = MIN(xx,xbot); xtop = MAX(xx,xtop);
970 ybot = MIN(yy,ybot); ytop = MAX(yy,ytop);
971 zbot = MIN(zz,zbot); ztop = MAX(zz,ztop);
972
973 LOAD_FVEC3( corner , nx1,ny1,nz0 ) ;
974 wcorn = VECMAT_VEC( ijk_to_dicom_in , corner ) ;
975 wfunc( wcorn.xyz[0],wcorn.xyz[1],wcorn.xyz[2] , &xx,&yy,&zz ) ;
976 xbot = MIN(xx,xbot); xtop = MAX(xx,xtop);
977 ybot = MIN(yy,ybot); ytop = MAX(yy,ytop);
978 zbot = MIN(zz,zbot); ztop = MAX(zz,ztop);
979
980 LOAD_FVEC3( corner , nx0,ny0,nz1 ) ;
981 wcorn = VECMAT_VEC( ijk_to_dicom_in , corner ) ;
982 wfunc( wcorn.xyz[0],wcorn.xyz[1],wcorn.xyz[2] , &xx,&yy,&zz ) ;
983 xbot = MIN(xx,xbot); xtop = MAX(xx,xtop);
984 ybot = MIN(yy,ybot); ytop = MAX(yy,ytop);
985 zbot = MIN(zz,zbot); ztop = MAX(zz,ztop);
986
987 LOAD_FVEC3( corner , nx1,ny0,nz1 ) ;
988 wcorn = VECMAT_VEC( ijk_to_dicom_in , corner ) ;
989 wfunc( wcorn.xyz[0],wcorn.xyz[1],wcorn.xyz[2] , &xx,&yy,&zz ) ;
990 xbot = MIN(xx,xbot); xtop = MAX(xx,xtop);
991 ybot = MIN(yy,ybot); ytop = MAX(yy,ytop);
992 zbot = MIN(zz,zbot); ztop = MAX(zz,ztop);
993
994 LOAD_FVEC3( corner , nx0,ny1,nz1 ) ;
995 wcorn = VECMAT_VEC( ijk_to_dicom_in , corner ) ;
996 wfunc( wcorn.xyz[0],wcorn.xyz[1],wcorn.xyz[2] , &xx,&yy,&zz ) ;
997 xbot = MIN(xx,xbot); xtop = MAX(xx,xtop);
998 ybot = MIN(yy,ybot); ytop = MAX(yy,ytop);
999 zbot = MIN(zz,zbot); ztop = MAX(zz,ztop);
1000
1001 LOAD_FVEC3( corner , nx1,ny1,nz1 ) ;
1002 wcorn = VECMAT_VEC( ijk_to_dicom_in , corner ) ;
1003 wfunc( wcorn.xyz[0],wcorn.xyz[1],wcorn.xyz[2] , &xx,&yy,&zz ) ;
1004 xbot = MIN(xx,xbot); xtop = MAX(xx,xtop);
1005 ybot = MIN(yy,ybot); ytop = MAX(yy,ytop);
1006 zbot = MIN(zz,zbot); ztop = MAX(zz,ztop);
1007
1008 *xb = xbot; *xt = xtop;
1009 *yb = ybot; *yt = ytop; *zb = zbot; *zt = ztop; return;
1010 }
1011
1012 /*--------------------------------------------------------------------------*/
1013 /*! Geometrically transform a 3D dataset, producing a new dataset.
1014 - w_in2out transforms DICOM coords from input grid to output grid
1015 (only needed if newggg != NULL and flag == WARP3D_NEWGRID;
1016 otherwise, w_in2out can be NULL)
1017 - w_out2in is the inverse of w_in2out (cannot be NULL)
1018 - newggg is a pointer to something that determines the new
1019 grid on which the output dataset will be computed
1020 - NULL means compute on the inset grid
1021 - otherwise, flag determines how newggg is interpreted:
1022 - WARP3D_NEWGRID => newggg is a "float *" to a new grid size
1023 - WARP3D_NEWDSET => newggg is a "THD_3dim_dataset *" to a dataset
1024 whose grid will be used for the output grid
1025 - prefix = new dataset prefix (if NULL, then "warped")
1026 - zpad = number of planes to zeropad on each face of inset (>= 0)
1027 - flag = see above
1028 - Interpolation method can be set using mri_warp3D_method().
1029 - At end, input dataset is unloaded from memory.
1030 - If input is 3D+time, the output dataset won't have slice-wise time
1031 offsets, even if the input did.
1032 ----------------------------------------------------------------------------*/
1033
THD_warp3D(THD_3dim_dataset * inset,void w_in2out (float,float,float,float *,float *,float *),void w_out2in (float,float,float,float *,float *,float *),void * newggg,char * prefix,int zpad,int flag)1034 THD_3dim_dataset * THD_warp3D(
1035 THD_3dim_dataset *inset ,
1036 void w_in2out(float,float,float,float *,float *,float *),
1037 void w_out2in(float,float,float,float *,float *,float *),
1038 void *newggg , char *prefix , int zpad , int flag )
1039 {
1040 THD_3dim_dataset *outset , *qset , *newgset=NULL ;
1041 int nxin,nyin,nzin,nvals , ival ;
1042 int nxout,nyout,nzout ;
1043 float xbot,xtop , ybot,ytop , zbot,ztop ;
1044 int use_newgrid , use_oldgrid ;
1045 float ddd_newgrid=0.0 , fac ;
1046 MRI_IMAGE *inim , *outim , *wim ;
1047
1048 ENTRY("THD_warp3D") ;
1049
1050 /*-- decide which grid the output will be computed on --*/
1051
1052 use_newgrid = ((flag & WARP3D_GRIDMASK) == WARP3D_NEWGRID) && (newggg != NULL) ;
1053 if( use_newgrid ){
1054 float *gg = (float *)newggg ;
1055 if( thd_floatscan(1,gg) == 0 && *gg > 0.0 ) ddd_newgrid = *gg ;
1056 else use_newgrid = 0 ;
1057 } else if( newggg != NULL && (flag & WARP3D_GRIDMASK) == WARP3D_NEWDSET ){
1058 newgset = (THD_3dim_dataset *) newggg ;
1059 if( !ISVALID_DSET(newgset) ) newgset = NULL ;
1060 }
1061 use_oldgrid = !use_newgrid && (newgset == NULL) ;
1062
1063 /*-- see if inputs are valid --*/
1064
1065 if( !ISVALID_DSET(inset) ||
1066 w_out2in == NULL ||
1067 (w_in2out == NULL && use_newgrid) ){
1068
1069 fprintf(stderr,"** ERROR: THD_warp3D has bad inputs!\n") ;
1070 RETURN(NULL);
1071 }
1072
1073 /*-- zeropad and replace input dataset, if desired --*/
1074
1075 if( zpad > 0 ){
1076 qset = THD_zeropad( inset , zpad,zpad,zpad,zpad,zpad,zpad ,
1077 "Quetzalcoatl" , ZPAD_PURGE ) ;
1078 DSET_unload(inset) ;
1079 if( qset == NULL ){
1080 fprintf(stderr,"** ERROR: THD_warp3D can't zeropad!\n"); RETURN(NULL);
1081 }
1082 } else {
1083 qset = inset ;
1084 }
1085
1086 /*-- read input data from disk, if not already present --*/
1087
1088 DSET_load(qset) ;
1089 if( !DSET_LOADED(qset) ){
1090 fprintf(stderr,"** ERROR: THD_warp3D can't load input dataset!\n") ;
1091 if( qset != inset ) DSET_delete(qset) ; else DSET_unload(qset) ;
1092 RETURN(NULL);
1093 }
1094
1095 /*-- compute mapping from input dataset (i,j,k) to DICOM coords --*/
1096
1097 { THD_vecmat ijk_to_xyz , xyz_to_dicom ;
1098
1099 LOAD_DIAG_MAT( ijk_to_xyz.mm , qset->daxes->xxdel,
1100 qset->daxes->yydel, qset->daxes->zzdel );
1101 LOAD_FVEC3 ( ijk_to_xyz.vv , qset->daxes->xxorg,
1102 qset->daxes->yyorg, qset->daxes->zzorg );
1103
1104 xyz_to_dicom.mm = qset->daxes->to_dicomm ;
1105 LOAD_FVEC3( xyz_to_dicom.vv , 0.0,0.0,0.0 ) ;
1106
1107 ijk_to_dicom_in = MUL_VECMAT( xyz_to_dicom , ijk_to_xyz ) ;
1108 dicom_to_ijk_in = INV_VECMAT( ijk_to_dicom_in ) ;
1109 }
1110
1111 /*-- make empty output dataset --*/
1112
1113 nxin = DSET_NX(qset) ;
1114 nyin = DSET_NY(qset) ;
1115 nzin = DSET_NZ(qset) ;
1116 nvals = DSET_NVALS(qset) ;
1117
1118 if( nxin*nyin*nzin <= 1 ){
1119 fprintf(stderr,"** ERROR: THD_warp3D has nxin=%d nyin=%d nzin=%d!\n",
1120 nxin,nyin,nzin ) ;
1121 if( qset != inset ) DSET_delete(qset) ; else DSET_unload(qset) ;
1122 RETURN(NULL) ;
1123 }
1124
1125 if( use_oldgrid ){ /*-- output is on same grid as input --*/
1126
1127 outset = EDIT_empty_copy( qset ) ;
1128
1129 } else if( newgset != NULL ){ /*-- output on same grid as newgset --*/
1130
1131 outset = EDIT_empty_copy( newgset ) ;
1132
1133 EDIT_dset_items( outset ,
1134 ADN_nvals , nvals ,
1135 ADN_type , qset->type ,
1136 ADN_view_type , qset->view_type ,
1137 ADN_func_type , qset->func_type ,
1138 ADN_malloc_type , DATABLOCK_MEM_MALLOC ,
1139 ADN_none ) ;
1140
1141 if( DSET_NUM_TIMES(qset) > 1 ) /* and some time structure? */
1142 EDIT_dset_items( outset ,
1143 ADN_ntt , nvals ,
1144 ADN_tunits , DSET_TIMEUNITS(qset) ,
1145 ADN_ttorg , DSET_TIMEORIGIN(qset) ,
1146 ADN_ttdel , DSET_TR(qset) ,
1147 ADN_ttdur , DSET_TIMEDURATION(qset) ,
1148 ADN_none ) ;
1149
1150 THD_copy_datablock_auxdata( qset->dblk , outset->dblk ) ; /* 08 Jun 2004 */
1151 INIT_STAT_AUX( qset , MAX_STAT_AUX , outset->stat_aux ) ;
1152
1153 } else { /*-- output is on new grid --*/
1154
1155 float xmid,ymid,zmid ;
1156 THD_ivec3 nxyz , orixyz ;
1157 THD_fvec3 dxyz , orgxyz ;
1158
1159 /* compute DICOM coordinates of warped corners */
1160
1161 warp_corners( qset, w_in2out, &xbot,&xtop, &ybot,&ytop, &zbot,&ztop ) ;
1162
1163 nxout = (int)( (xtop-xbot)/ddd_newgrid+0.999 ); if( nxout < 1 ) nxout = 1;
1164 nyout = (int)( (ytop-ybot)/ddd_newgrid+0.999 ); if( nyout < 1 ) nyout = 1;
1165 nzout = (int)( (ztop-zbot)/ddd_newgrid+0.999 ); if( nzout < 1 ) nzout = 1;
1166
1167 xmid = 0.5*(xbot+xtop); ymid = 0.5*(ybot+ytop); zmid = 0.5*(zbot+ztop);
1168 xbot = xmid-0.5*(nxout-1)*ddd_newgrid; xtop = xbot+(nxout-1)*ddd_newgrid;
1169 ybot = ymid-0.5*(nyout-1)*ddd_newgrid; ytop = ybot+(nyout-1)*ddd_newgrid;
1170 zbot = zmid-0.5*(nzout-1)*ddd_newgrid; ztop = zbot+(nzout-1)*ddd_newgrid;
1171
1172 #if 0
1173 if( verb )
1174 fprintf(stderr,"++ Transformed grid:\n"
1175 "++ xbot = %10.4g xtop = %10.4g nx = %d\n"
1176 "++ ybot = %10.4g ytop = %10.4g ny = %d\n"
1177 "++ zbot = %10.4g ztop = %10.4g nz = %d\n" ,
1178 xbot,xtop,nxout , ybot,ytop,nyout , zbot,ztop,nzout ) ;
1179 #endif
1180
1181 if( nxout*nyout*nzout <= 1 ){
1182 fprintf(stderr,"** ERROR: THD_warp3D has nxout=%d nyout=%d nzout=%d!\n",
1183 nxout,nyout,nzout ) ;
1184 if( qset != inset ) DSET_delete(qset) ; else DSET_unload(qset) ;
1185 RETURN(NULL) ;
1186 }
1187
1188 nxyz.ijk[0] = nxout ; dxyz.xyz[0] = ddd_newgrid ; /* setup axes */
1189 nxyz.ijk[1] = nyout ; dxyz.xyz[1] = ddd_newgrid ;
1190 nxyz.ijk[2] = nzout ; dxyz.xyz[2] = ddd_newgrid ;
1191
1192 orixyz.ijk[0] = ORI_R2L_TYPE ; orgxyz.xyz[0] = xbot ;
1193 orixyz.ijk[1] = ORI_A2P_TYPE ; orgxyz.xyz[1] = ybot ;
1194 orixyz.ijk[2] = ORI_I2S_TYPE ; orgxyz.xyz[2] = zbot ;
1195
1196 /** create dataset and mangle it into the desired shape **/
1197
1198 outset = EDIT_empty_copy( NULL ) ; /* void and formless dataset */
1199
1200 EDIT_dset_items( outset , /* give it some structure! */
1201 ADN_nxyz , nxyz ,
1202 ADN_xyzdel , dxyz ,
1203 ADN_xyzorg , orgxyz ,
1204 ADN_xyzorient , orixyz ,
1205 ADN_malloc_type , DATABLOCK_MEM_MALLOC ,
1206 ADN_nvals , nvals ,
1207 ADN_type , qset->type ,
1208 ADN_view_type , qset->view_type ,
1209 ADN_func_type , qset->func_type ,
1210 ADN_none ) ;
1211
1212 if( DSET_NUM_TIMES(qset) > 1 ) /* and some time structure? */
1213 EDIT_dset_items( outset ,
1214 ADN_ntt , nvals ,
1215 ADN_tunits , DSET_TIMEUNITS(qset) ,
1216 ADN_ttorg , DSET_TIMEORIGIN(qset) ,
1217 ADN_ttdel , DSET_TR(qset) ,
1218 ADN_ttdur , DSET_TIMEDURATION(qset) ,
1219 ADN_none ) ;
1220
1221 THD_copy_datablock_auxdata( qset->dblk , outset->dblk ) ; /* 08 Jun 2004 */
1222 INIT_STAT_AUX( qset , MAX_STAT_AUX , outset->stat_aux ) ;
1223
1224 } /*-- end of warping to new grid --*/
1225
1226 nxout = DSET_NX(outset) ;
1227 nyout = DSET_NY(outset) ;
1228 nzout = DSET_NZ(outset) ;
1229
1230 /*-- compute mapping from output dataset (i,j,k) to DICOM coords --*/
1231
1232 { THD_vecmat ijk_to_xyz , xyz_to_dicom ;
1233
1234 LOAD_DIAG_MAT( ijk_to_xyz.mm, outset->daxes->xxdel,
1235 outset->daxes->yydel, outset->daxes->zzdel );
1236 LOAD_FVEC3 ( ijk_to_xyz.vv, outset->daxes->xxorg,
1237 outset->daxes->yyorg, outset->daxes->zzorg );
1238
1239 xyz_to_dicom.mm = outset->daxes->to_dicomm ;
1240 LOAD_FVEC3( xyz_to_dicom.vv , 0.0,0.0,0.0 ) ;
1241
1242 ijk_to_dicom_out = MUL_VECMAT( xyz_to_dicom , ijk_to_xyz ) ;
1243 }
1244
1245 /*-- add prefix to new dataset --*/
1246
1247 if( !THD_filename_ok(prefix) ) prefix = "warped" ;
1248 EDIT_dset_items( outset , ADN_prefix,prefix , ADN_none ) ;
1249
1250 /*-- loop over bricks and warp them --*/
1251
1252 warp_out_to_in = w_out2in ; /* for use in w3d_warp_func(), supra */
1253
1254 for( ival=0 ; ival < nvals ; ival++ ){
1255 inim = DSET_BRICK(qset,ival) ;
1256 fac = DSET_BRICK_FACTOR(qset,ival) ;
1257 if( fac > 0.0 && fac != 1.0 ) wim = mri_scale_to_float( fac , inim ) ;
1258 else wim = inim ;
1259 outim = mri_warp3D( wim , nxout,nyout,nzout , w3d_warp_func ) ;
1260 if( outim == NULL ){
1261 fprintf(stderr,"** ERROR: THD_warp3D fails at ival=%d\n",ival);
1262 DSET_delete(outset);
1263 if( qset != inset ) DSET_delete(qset) ; else DSET_unload(qset) ;
1264 RETURN(NULL) ;
1265 }
1266
1267 if( wim != inim ){ /* if we scaled input */
1268 float gtop , fimfac ;
1269 mri_free(wim) ;
1270
1271 /* 20 Oct 2003: rescale if input was an integer type */
1272
1273 if( outim->kind == MRI_float && MRI_IS_INT_TYPE(inim->kind) ){
1274 fimfac = fac ;
1275 wim = mri_scalize( outim , inim->kind , &fimfac ) ;
1276 EDIT_BRICK_FACTOR( outset , ival , fimfac ) ;
1277 mri_free(outim) ; outim = wim ;
1278 } else { /* input not an integer: */
1279 EDIT_BRICK_FACTOR( outset , ival , 0.0 ) ; /* output=unscaled float */
1280 }
1281
1282 }
1283 EDIT_substitute_brick( outset, ival,outim->kind,mri_data_pointer(outim) );
1284 DSET_unload_one( qset , ival ) ;
1285 }
1286
1287 /*-- done!!! --*/
1288
1289 if( qset != inset ) DSET_delete(qset) ; else DSET_unload(qset) ;
1290
1291 THD_load_statistics( outset ) ;
1292 RETURN( outset ) ;
1293 }
1294
1295 /*--------------------------------------------------------------------------*/
1296
1297 static THD_vecmat aff_out2in , aff_in2out ; /* data for afo2i() and afi2o() */
1298
1299 /*! Internal transform for THD_warp3D_affine(). ----------------------------*/
1300
afo2i(float a,float b,float c,float * x,float * y,float * z)1301 static INLINE void afo2i( float a,float b,float c,
1302 float *x,float *y,float *z )
1303 {
1304 THD_fvec3 xxx , yyy ;
1305 LOAD_FVEC3( xxx , a,b,c ) ;
1306 yyy = VECMAT_VEC( aff_out2in , xxx ) ;
1307 UNLOAD_FVEC3( yyy , *x , *y , *z ) ;
1308 }
1309
1310 /*! Internal transform for THD_warp3D_affine(). ----------------------------*/
1311
afi2o(float a,float b,float c,float * x,float * y,float * z)1312 static INLINE void afi2o( float a,float b,float c,
1313 float *x,float *y,float *z )
1314 {
1315 THD_fvec3 xxx , yyy ;
1316 LOAD_FVEC3( xxx , a,b,c ) ;
1317 yyy = VECMAT_VEC( aff_in2out , xxx ) ;
1318 UNLOAD_FVEC3( yyy , *x , *y , *z ) ;
1319 }
1320
1321 /*--------------------------------------------------------------------------*/
1322 /*! A version of THD_warp3D() that uses an affine map, specified in the
1323 out2in parameter. The main advantage over THD_warp3D() is that the
1324 transformation functions internal -- afo2i() and afi2o() -- which
1325 makes it simpler on the caller and also somewhat faster.
1326 ----------------------------------------------------------------------------*/
1327
THD_warp3D_affine(THD_3dim_dataset * inset,THD_vecmat out2in,void * newggg,char * prefix,int zpad,int flag)1328 THD_3dim_dataset * THD_warp3D_affine(
1329 THD_3dim_dataset *inset ,
1330 THD_vecmat out2in ,
1331 void *newggg , char *prefix , int zpad , int flag )
1332 {
1333 aff_out2in = out2in ;
1334 aff_in2out = INV_VECMAT(aff_out2in) ;
1335
1336 return THD_warp3D( inset , afi2o,afo2i , newggg,prefix,zpad,flag ) ;
1337 }
1338
1339 /*--------------------------------------------------------------------------*/
1340 /*! Internal transform functions for TT <-> MNI coords. Both are RAI,
1341 as per AFNI internal logic. 11 Mar 2004 - RW Cox [Jury Duty Day]
1342 cf. https://imaging.mrc-cbu.cam.ac.uk/imaging/MniTalairach
1343 ----------------------------------------------------------------------------*/
1344
w3d_mni2tta(float mx,float my,float mz,float * tx,float * ty,float * tz)1345 static INLINE void w3d_mni2tta( float mx , float my , float mz ,
1346 float *tx, float *ty, float *tz )
1347 {
1348 *tx = 0.99 * mx ;
1349
1350 if( mz > 0.0 ){
1351 *ty = 0.9688 * my + 0.0460 * mz ;
1352 *tz = -0.0485 * my + 0.9189 * mz ;
1353 } else {
1354 *ty = 0.9688 * my + 0.0420 * mz ;
1355 *tz = -0.0485 * my + 0.8390 * mz ;
1356 }
1357 }
1358
1359 /*------- Should be the inverse of the above! -----------------------------*/
1360
w3d_tta2mni(float tx,float ty,float tz,float * mx,float * my,float * mz)1361 static INLINE void w3d_tta2mni( float tx , float ty , float tz ,
1362 float *mx, float *my, float *mz )
1363 {
1364 *mx = 1.01010 * tx ;
1365 *my = 1.02962 * ty - 0.05154 * tz ;
1366 *mz = 0.05434 * ty + 1.08554 * tz ;
1367 if( *mz < 0.0 ) *mz *= 1.09523 ;
1368 }
1369
1370 /*--------------------------------------------------------------------------*/
1371
THD_warp3D_tta2mni(THD_3dim_dataset * inset,void * newggg,char * prefix,int zpad,int flag)1372 THD_3dim_dataset * THD_warp3D_tta2mni( THD_3dim_dataset *inset , void *newggg ,
1373 char *prefix , int zpad , int flag )
1374 {
1375 return THD_warp3D( inset , w3d_tta2mni,w3d_mni2tta , newggg,prefix,zpad,flag ) ;
1376 }
1377
1378 /*------- Should be the inverse of the above! -----------------------------*/
1379
THD_warp3D_mni2tta(THD_3dim_dataset * inset,void * newggg,char * prefix,int zpad,int flag)1380 THD_3dim_dataset * THD_warp3D_mni2tta( THD_3dim_dataset *inset , void *newggg ,
1381 char *prefix , int zpad , int flag )
1382 {
1383 return THD_warp3D( inset , w3d_mni2tta,w3d_tta2mni , newggg,prefix,zpad,flag ) ;
1384 }
1385