1 /*****************************************************************************
2 Major portions of this software are copyrighted by the Medical College
3 of Wisconsin, 1994-2000, and are released under the Gnu General Public
4 License, Version 2. See the file README.Copyright for details.
5 ******************************************************************************/
6
7 /*===========================================================================
8 Routines to rotate and shift a 3D volume using a 4 way shear decomposition,
9 coupled with an FFT-based shifting of the rows -- RWCox [October 1998].
10 =============================================================================*/
11
12 #include "thd_shear3d.h" /* 23 Oct 2000: moved shear funcs to thd_shear3d.c */
13
14 /********** 28 Oct 1999: the shifting routines that were here **********
15 ********** have been removed to file thd_shift2.c **********/
16
17 /*---------------------------------------------------------------------------
18 Set the interpolation method for shifting:
19 input is one of MRI_NN, MRI_LINEAR, MRI_CUBIC, or MRI_FOURIER.
20 -----------------------------------------------------------------------------*/
21
22 typedef void (*shift_func)(int,int,float,float *,float,float *) ;
23 static shift_func shifter = fft_shift2 ;
24 static int shift_method = MRI_FOURIER ;
25
THD_rota_method(int mode)26 void THD_rota_method( int mode )
27 {
28 shift_method = mode ;
29 switch( mode ){
30 case MRI_NN: shifter = nn_shift2 ; break ;
31 case MRI_TSSHIFT: shifter = ts_shift2 ; break ; /* Dec 1999 */
32
33 case MRI_LINEAR: shifter = lin_shift2 ; break ;
34 case MRI_FOURIER: shifter = fft_shift2 ; break ;
35 default:
36 case MRI_CUBIC: shifter = cub_shift2 ; break ;
37 case MRI_QUINTIC: shifter = quint_shift2 ; break ; /* Nov 1998 */
38 case MRI_HEPTIC: shifter = hept_shift2 ; break ; /* Nov 1998 */
39
40 case MRI_FOURIER_NOPAD: shifter = fft_shift2 ; break ; /* 13 May 2003 */
41 }
42 return ;
43 }
44
45 /*---------------------------------------------------------------------------
46 Flip a 3D array about the (x,y) axes:
47 i <--> nx-1-i j <--> ny-1-j
48 -----------------------------------------------------------------------------*/
49
50 #define VV(i,j,k) v[(i)+(j)*nx+(k)*nxy]
51 #define SX(i) (nx1-(i))
52 #define SY(j) (ny1-(j))
53 #define SZ(k) (nz1-(k))
54
flip_xy(int nx,int ny,int nz,float * v)55 static void flip_xy( int nx , int ny , int nz , float * v )
56 {
57 int ii,jj,kk ;
58 int nx1=nx-1,nx2=nx/2, ny1=ny-1,ny2=ny/2, nz1=nz-1,nz2=nz/2, nxy=nx*ny ;
59 float * r1 ;
60
61 r1 = (float *) malloc(sizeof(float)*nx) ; /* save 1 row */
62
63 for( kk=0 ; kk < nz ; kk++ ){ /* for each slice */
64 for( jj=0 ; jj < ny2 ; jj++ ){ /* first 1/2 of rows */
65
66 /* swap rows jj and ny1-jj, flipping them in ii as well */
67
68 for( ii=0 ; ii < nx ; ii++ ) r1[ii] = VV(SX(ii),SY(jj),kk) ;
69 for( ii=0 ; ii < nx ; ii++ ) VV(ii,SY(jj),kk) = VV(SX(ii),jj ,kk) ;
70 for( ii=0 ; ii < nx ; ii++ ) VV(ii,jj ,kk) = r1[ii] ;
71 }
72 if( ny%2 == 1 ){ /* central row? */
73 for( ii=0 ; ii < nx ; ii++ ) r1[ii] = VV(SX(ii),jj,kk) ; /* flip it */
74 for( ii=0 ; ii < nx ; ii++ ) VV(ii,jj,kk) = r1[ii] ; /* restore */
75 }
76 }
77
78 free(r1) ; return ;
79 }
80
81 /*---------------------------------------------------------------------------
82 Flip a 3D array about the (y,z) axes:
83 j <--> ny-1-j k <--> nz-1-k
84 -----------------------------------------------------------------------------*/
85
flip_yz(int nx,int ny,int nz,float * v)86 static void flip_yz( int nx , int ny , int nz , float * v )
87 {
88 int ii,jj,kk ;
89 int nx1=nx-1,nx2=nx/2, ny1=ny-1,ny2=ny/2, nz1=nz-1,nz2=nz/2, nxy=nx*ny ;
90 float * r1 ;
91
92 r1 = (float *) malloc(sizeof(float)*ny) ;
93
94 for( ii=0 ; ii < nx ; ii++ ){
95 for( kk=0 ; kk < nz2 ; kk++ ){
96 for( jj=0 ; jj < ny ; jj++ ) r1[jj] = VV(ii,SY(jj),SZ(kk)) ;
97 for( jj=0 ; jj < ny ; jj++ ) VV(ii,jj,SZ(kk)) = VV(ii,SY(jj),kk ) ;
98 for( jj=0 ; jj < ny ; jj++ ) VV(ii,jj,kk ) = r1[jj] ;
99 }
100 if( nz%2 == 1 ){
101 for( jj=0 ; jj < ny ; jj++ ) r1[jj] = VV(ii,SY(jj),kk) ;
102 for( jj=0 ; jj < ny ; jj++ ) VV(ii,jj,kk) = r1[jj] ;
103 }
104 }
105
106 free(r1) ; return ;
107 }
108
109 /*---------------------------------------------------------------------------
110 Flip a 3D array about the (x,z) axes:
111 i <--> nx-1-i k <--> nz-1-k
112 -----------------------------------------------------------------------------*/
113
flip_xz(int nx,int ny,int nz,float * v)114 static void flip_xz( int nx , int ny , int nz , float * v )
115 {
116 int ii,jj,kk ;
117 int nx1=nx-1,nx2=nx/2, ny1=ny-1,ny2=ny/2, nz1=nz-1,nz2=nz/2, nxy=nx*ny ;
118 float * r1 ;
119
120 r1 = (float *) malloc(sizeof(float)*nx) ;
121
122 for( jj=0 ; jj < ny ; jj++ ){
123 for( kk=0 ; kk < nz2 ; kk++ ){
124 for( ii=0 ; ii < nx ; ii++ ) r1[ii] = VV(SX(ii),jj,SZ(kk)) ;
125 for( ii=0 ; ii < nx ; ii++ ) VV(ii,jj,SZ(kk)) = VV(SX(ii),jj,kk ) ;
126 for( ii=0 ; ii < nx ; ii++ ) VV(ii,jj,kk ) = r1[ii] ;
127 }
128 if( nz%2 == 1 ){
129 for( ii=0 ; ii < nx ; ii++ ) r1[ii] = VV(SX(ii),jj,kk) ;
130 for( ii=0 ; ii < nx ; ii++ ) VV(ii,jj,kk) = r1[ii] ;
131 }
132 }
133
134 free(r1) ; return ;
135 }
136
137 /*---------------------------------------------------------------------------
138 Apply an x-axis shear to a 3D array: x -> x + a*y + b*z + s
139 (dilation factor "f" assumed to be 1.0)
140 -----------------------------------------------------------------------------*/
141
apply_xshear(float a,float b,float s,int nx,int ny,int nz,float * v)142 static void apply_xshear( float a , float b , float s ,
143 int nx , int ny , int nz , float * v )
144 {
145 float * fj0 , * fj1 ;
146 int nx1=nx-1 , ny1=ny-1 , nz1=nz-1 , nxy=nx*ny ;
147 float nx2=0.5*nx1 , ny2=0.5*ny1 , nz2=0.5*nz1 ;
148 int ii,jj,kk , nup=0,nst ;
149 float a0 , a1=0 , st ;
150
151 ENTRY("apply_xshear") ;
152
153 /* don't do anything if shift is less than 0.001 pixel */
154
155 st = fabs(a)*ny2 + fabs(b)*nz2 + fabs(s) ; if( st < 1.e-3 ) EXRETURN ;
156
157 if( shift_method == MRI_FOURIER ){
158 nst = nx + 0.5*st ;
159 nup = csfft_nextup_one35(nst) ;
160 } else if( shift_method == MRI_FOURIER_NOPAD ){
161 nup = csfft_nextup_even(nx) ;
162 }
163
164 for( kk=0 ; kk < nz ; kk++ ){
165 for( jj=0 ; jj < ny ; jj+=2 ){
166 fj0 = v + (jj*nx + kk*nxy) ;
167 fj1 = (jj < ny1) ? (fj0 + nx) : NULL ; /* allow for odd ny */
168 a0 = a*(jj-ny2) + b*(kk-nz2) + s ;
169 a1 = a0 + a ;
170 shifter( nx , nup , a0 , fj0 , a1 , fj1 ) ;
171 }
172 }
173
174 EXRETURN ;
175 }
176
177 /*---------------------------------------------------------------------------
178 Apply a y-axis shear to a 3D array: y -> y + a*x + b*z + s
179 -----------------------------------------------------------------------------*/
180
apply_yshear(float a,float b,float s,int nx,int ny,int nz,float * v)181 static void apply_yshear( float a , float b , float s ,
182 int nx , int ny , int nz , float * v )
183 {
184 float * fj0 , * fj1 ;
185 int nx1=nx-1 , ny1=ny-1 , nz1=nz-1 , nxy=nx*ny ;
186 float nx2=0.5*nx1 , ny2=0.5*ny1 , nz2=0.5*nz1 ;
187 int ii,jj,kk , nup=0,nst ;
188 float a0 , a1=0 , st ;
189
190 ENTRY("apply_yshear") ;
191
192 /* don't do anything if shift is less than 0.001 pixel */
193
194 st = fabs(a)*nx2 + fabs(b)*nz2 + fabs(s) ; if( st < 1.e-3 ) EXRETURN ;
195
196 if( shift_method == MRI_FOURIER ){
197 nst = ny + 0.5*st ;
198 nup = csfft_nextup_one35(nst) ;
199 } else if( shift_method == MRI_FOURIER_NOPAD ){
200 nup = csfft_nextup_even(ny) ;
201 }
202
203 fj0 = (float *) malloc( sizeof(float) * 2*ny ) ; fj1 = fj0 + ny ;
204
205 for( kk=0 ; kk < nz ; kk++ ){
206 for( ii=0 ; ii < nx1 ; ii+=2 ){
207 for( jj=0; jj < ny; jj++ ){ fj0[jj] = VV(ii,jj,kk) ; fj1[jj] = VV(ii+1,jj,kk) ; }
208 a0 = a*(ii-nx2) + b*(kk-nz2) + s ;
209 a1 = a0 + a ;
210 shifter( ny , nup , a0 , fj0 , a1 , fj1 ) ;
211 for( jj=0; jj < ny; jj++ ){ VV(ii,jj,kk) = fj0[jj] ; VV(ii+1,jj,kk) = fj1[jj] ; }
212 }
213
214 if( ii == nx1 ){ /* allow for odd nx */
215 for( jj=0; jj < ny; jj++ ) fj0[jj] = VV(ii,jj,kk) ;
216 a0 = a*(ii-nx2) + b*(kk-nz2) + s ;
217 shifter( ny , nup , a0 , fj0 , a1 , NULL ) ;
218 for( jj=0; jj < ny; jj++ ) VV(ii,jj,kk) = fj0[jj] ;
219 }
220 }
221
222 free(fj0) ; EXRETURN ;
223 }
224
225 /*---------------------------------------------------------------------------
226 Apply a z-axis shear to a 3D array: z -> z + a*x + b*y + s
227 -----------------------------------------------------------------------------*/
228
apply_zshear(float a,float b,float s,int nx,int ny,int nz,float * v)229 static void apply_zshear( float a , float b , float s ,
230 int nx , int ny , int nz , float * v )
231 {
232 float * fj0 , * fj1 ;
233 int nx1=nx-1 , ny1=ny-1 , nz1=nz-1 , nxy=nx*ny ;
234 float nx2=0.5*nx1 , ny2=0.5*ny1 , nz2=0.5*nz1 ;
235 int ii,jj,kk , nup=0,nst ;
236 float a0 , a1=0 , st ;
237
238 ENTRY("apply_zshear") ;
239
240 /* don't do anything if shift is less than 0.001 pixel */
241
242 st = fabs(a)*nx2 + fabs(b)*ny2 + fabs(s) ; if( st < 1.e-3 ) EXRETURN ;
243
244 if( shift_method == MRI_FOURIER ){
245 nst = nz + 0.5*st ;
246 nup = csfft_nextup_one35(nst) ;
247 } else if( shift_method == MRI_FOURIER_NOPAD ){
248 nup = csfft_nextup_even(nz) ;
249 }
250
251 fj0 = (float *) malloc( sizeof(float) * 2*nz ) ; fj1 = fj0 + nz ;
252
253 for( jj=0 ; jj < ny ; jj++ ){
254 for( ii=0 ; ii < nx1 ; ii+=2 ){
255 for( kk=0; kk < nz; kk++ ){ fj0[kk] = VV(ii,jj,kk) ; fj1[kk] = VV(ii+1,jj,kk) ; }
256 a0 = a*(ii-nx2) + b*(jj-ny2) + s ;
257 a1 = a0 + a ;
258 shifter( nz , nup , a0 , fj0 , a1 , fj1 ) ;
259 for( kk=0; kk < nz; kk++ ){ VV(ii,jj,kk) = fj0[kk] ; VV(ii+1,jj,kk) = fj1[kk] ; }
260 }
261
262 if( ii == nx1 ){ /* allow for odd nx */
263 for( kk=0; kk < nz; kk++ ) fj0[kk] = VV(ii,jj,kk) ;
264 a0 = a*(ii-nx2) + b*(jj-ny2) + s ;
265 shifter( nz , nup , a0 , fj0 , a1 , NULL ) ;
266 for( kk=0; kk < nz; kk++ ) VV(ii,jj,kk) = fj0[kk] ;
267 }
268 }
269
270 free(fj0) ; EXRETURN ;
271 }
272
273 /*---------------------------------------------------------------------------
274 Apply a set of shears to a 3D array of floats.
275 Note that we assume that the dilation factors ("f") are all 1.
276 -----------------------------------------------------------------------------*/
277
apply_3shear(MCW_3shear shr,int nx,int ny,int nz,float * vol)278 static void apply_3shear( MCW_3shear shr ,
279 int nx , int ny , int nz , float * vol )
280 {
281 int qq ;
282 float a , b , s ;
283
284 ENTRY("apply_3shear") ;
285
286 if( ! ISVALID_3SHEAR(shr) ) EXRETURN ;
287
288 /* carry out a preliminary 180 flippo ? */
289
290 if( shr.flip0 >= 0 ){
291 switch( shr.flip0 + shr.flip1 ){
292 case 1: flip_xy( nx,ny,nz,vol ) ; break ;
293 case 2: flip_xz( nx,ny,nz,vol ) ; break ;
294 case 3: flip_yz( nx,ny,nz,vol ) ; break ;
295 default: EXRETURN ; /* should not occur */
296 }
297 }
298
299 /* apply each shear */
300
301 for( qq=0 ; qq < 4 ; qq++ ){
302 switch( shr.ax[qq] ){
303 case 0:
304 a = shr.scl[qq][1] ;
305 b = shr.scl[qq][2] ;
306 s = shr.sft[qq] ;
307 apply_xshear( a,b,s , nx,ny,nz , vol ) ;
308 break ;
309
310 case 1:
311 a = shr.scl[qq][0] ;
312 b = shr.scl[qq][2] ;
313 s = shr.sft[qq] ;
314 apply_yshear( a,b,s , nx,ny,nz , vol ) ;
315 break ;
316
317 case 2:
318 a = shr.scl[qq][0] ;
319 b = shr.scl[qq][1] ;
320 s = shr.sft[qq] ;
321 apply_zshear( a,b,s , nx,ny,nz , vol ) ;
322 break ;
323 }
324 }
325
326 EXRETURN ;
327 }
328
329 /*---------------------------------------------------------------------------
330 Set zero padding size for rotations:
331 padding is done before rotate, then stripped off afterwards.
332 02 Feb 2001 -- RWCox
333 -----------------------------------------------------------------------------*/
334
335 static int rotpx=0 , rotpy=0 , rotpz = 0 ;
336 static int rotpset=0 ;
337
THD_rota_setpad(int px,int py,int pz)338 void THD_rota_setpad( int px , int py , int pz )
339 {
340 rotpx = (px > 0) ? px : 0 ;
341 rotpy = (py > 0) ? py : 0 ;
342 rotpz = (pz > 0) ? pz : 0 ;
343 rotpset = 1 ; /* 05 Feb 2001 */
344 return ;
345 }
346
347 /*---------------------------------------------------------------------------*/
348
THD_rota_clearpad(void)349 void THD_rota_clearpad(void) /* 05 Feb 2001 */
350 {
351 rotpx=rotpy=rotpz=0; rotpset=1; return;
352 }
353
THD_rota_envpad(void)354 static void THD_rota_envpad(void)
355 {
356 char * eee = my_getenv("AFNI_ROTA_ZPAD") ;
357 int ppp ;
358
359 if( rotpset ) return ;
360 eee = my_getenv("AFNI_ROTA_ZPAD") ;
361 if( eee != NULL ){
362 ppp = strtol( eee , NULL , 10 ) ;
363 THD_rota_setpad(ppp,ppp,ppp) ;
364 }
365 rotpset = 1 ; return ;
366 }
367
368 /*---------------------------------------------------------------------------
369 Rotate and translate a 3D volume.
370 -----------------------------------------------------------------------------*/
371
372 #undef CLIPIT
373
THD_rota_vol(int nx,int ny,int nz,float xdel,float ydel,float zdel,float * vol,int ax1,float th1,int ax2,float th2,int ax3,float th3,int dcode,float dx,float dy,float dz)374 void THD_rota_vol( int nx , int ny , int nz ,
375 float xdel , float ydel , float zdel , float * vol ,
376 int ax1,float th1, int ax2,float th2, int ax3,float th3,
377 int dcode , float dx , float dy , float dz )
378 {
379 MCW_3shear shr ;
380
381 #ifdef CLIPIT
382 register float bot , top ;
383 register int nxyz=nx*ny*nz , ii ;
384 #endif
385
386 ENTRY("THD_rota_vol") ;
387
388 if( nx < 2 || ny < 2 || nz < 2 || vol == NULL ) EXRETURN ;
389
390 if( xdel == 0.0 ) xdel = 1.0 ;
391 if( ydel == 0.0 ) ydel = 1.0 ;
392 if( zdel == 0.0 ) zdel = 1.0 ;
393
394 if( th1 == 0.0 && th2 == 0.0 && th3 == 0.0 ){ /* nudge rotation */
395 th1 = 1.e-6 ; th2 = 1.1e-6 ; th3 = 0.9e-6 ;
396 }
397
398 #if 0
399 fprintf(stderr,"THD_rota_vol:\n") ;
400 fprintf(stderr," th1=%g th2=%g th3=%g\n",th1,th2,th3) ;
401 fprintf(stderr," dx=%g dy=%g dz=%g\n",dx,dy,dz) ;
402 fprintf(stderr," xdel=%g ydel=%g zdel=%g\n",xdel,ydel,zdel) ;
403 #endif
404
405 shr = rot_to_shear( ax1,-th1 , ax2,-th2 , ax3,-th3 ,
406 dcode,dx,dy,dz , xdel,ydel,zdel ) ;
407
408 if( ! ISVALID_3SHEAR(shr) ){
409 fprintf(stderr,"*** THD_rota_vol: can't compute shear transformation!\n") ;
410 EXRETURN ;
411 }
412
413 #if 0
414 if( MRILIB_verbose )
415 DUMP_3SHEAR("Computed shear",shr) ;
416 #endif
417
418 #ifdef CLIPIT
419 bot = top = vol[0] ;
420 for( ii=1 ; ii < nxyz ; ii++ ){
421 if( vol[ii] < bot ) bot = vol[ii] ;
422 else if( vol[ii] > top ) top = vol[ii] ;
423 }
424 if( bot >= top ) EXRETURN ;
425 #endif
426
427 /********************************/
428 /* 02 Feb 2001: include padding */
429
430 { float * vvv , *www ;
431 int nxp , nyp , nzp ;
432
433 THD_rota_envpad() ; /* 05 Feb 2001 */
434
435 nxp=nx+2*rotpx ; nyp=ny+2*rotpy ; nzp=nz+2*rotpz ;
436
437 if( rotpx > 0 && rotpy > 0 && rotpz > 0 )
438 vvv = EDIT_volpad_even( rotpx,rotpy,rotpz , nx,ny,nz , MRI_float,vol ) ;
439 else
440 vvv = vol ;
441
442 apply_3shear( shr , nxp,nyp,nzp , vvv ) ; /*-- do the actual rotation! --*/
443
444 if( vvv != vol ){
445 www = EDIT_volpad_even( -rotpx,-rotpy,-rotpz , nxp,nyp,nzp , MRI_float,vvv ) ;
446 free(vvv) ;
447 memcpy( vol , www , sizeof(float)*nx*ny*nz ) ; free(www) ;
448 }
449 }
450
451 /********************************/
452
453 #ifdef CLIPIT
454 for( ii=0 ; ii < nxyz ; ii++ ){
455 if( vol[ii] < bot ) vol[ii] = bot ;
456 else if( vol[ii] > top ) vol[ii] = top ;
457 }
458 #endif
459
460 EXRETURN ;
461 }
462
463 /*---------------------------------------------------------------------------
464 Like the above, but with geometrical information about the volume
465 given from the image header
466 -----------------------------------------------------------------------------*/
467
THD_rota3D(MRI_IMAGE * im,int ax1,float th1,int ax2,float th2,int ax3,float th3,int dcode,float dx,float dy,float dz)468 MRI_IMAGE * THD_rota3D( MRI_IMAGE * im ,
469 int ax1,float th1, int ax2,float th2, int ax3,float th3,
470 int dcode , float dx , float dy , float dz )
471 {
472 MRI_IMAGE * jm ;
473 float * jvol ;
474
475 if( ! MRI_IS_3D(im) ){
476 fprintf(stderr,"\n*** THD_rota3D: non-3D image input!\n") ;
477 return NULL ;
478 }
479
480 jm = mri_new_vol( im->nx , im->ny , im->nz , MRI_float ) ;
481 MRI_COPY_AUX(jm,im) ;
482 jvol = MRI_FLOAT_PTR(jm) ;
483
484 EDIT_coerce_type( im->nvox ,
485 im->kind , mri_data_pointer(im) , MRI_float , jvol ) ;
486
487 THD_rota_vol( im->nx , im->ny , im->nz ,
488 fabs(im->dx) , fabs(im->dy) , fabs(im->dz) , jvol ,
489 ax1,th1 , ax2,th2 , ax3,th3 , dcode , dx,dy,dz ) ;
490
491 return jm ;
492 }
493
494 /****************************************************************************
495 Alternative entries, with rotation specified via a 3x3 matrix
496 and shift as a 3-vector -- RWCox - 16 July 2000
497 *****************************************************************************/
498
499 /*---------------------------------------------------------------------------
500 Rotate and translate a 3D volume
501 -----------------------------------------------------------------------------*/
502
503 #undef CLIPIT
504
THD_rota_vol_matvec(int nx,int ny,int nz,float xdel,float ydel,float zdel,float * vol,THD_dmat33 rmat,THD_dfvec3 tvec)505 void THD_rota_vol_matvec( int nx , int ny , int nz ,
506 float xdel , float ydel , float zdel , float * vol ,
507 THD_dmat33 rmat , THD_dfvec3 tvec )
508 {
509 MCW_3shear shr ;
510 int dcode ;
511
512 #ifdef CLIPIT
513 register float bot , top ;
514 register int nxyz=nx*ny*nz , ii ;
515 #endif
516
517 if( nx < 2 || ny < 2 || nz < 2 || vol == NULL ) return ;
518
519 if( xdel == 0.0 ) xdel = 1.0 ;
520 if( ydel == 0.0 ) ydel = 1.0 ;
521 if( zdel == 0.0 ) zdel = 1.0 ;
522
523 shr = rot_to_shear_matvec( rmat , tvec , xdel,ydel,zdel ) ;
524
525 if( ! ISVALID_3SHEAR(shr) ){
526 fprintf(stderr,"*** THD_rota_vol: can't compute shear transformation!\n") ;
527 return ;
528 }
529
530 #if 0
531 DUMP_3SHEAR("Computed shear",shr) ;
532 DUMP_MAT33("rmat",rmat) ;
533 DUMP_FVEC3("tvec",tvec) ;
534 fprintf(stderr,"---- xdel = %f, ydel = %f, zdel = %f\n",
535 xdel, ydel, zdel);
536 #endif
537
538 #ifdef CLIPIT
539 bot = top = vol[0] ;
540 for( ii=1 ; ii < nxyz ; ii++ ){
541 if( vol[ii] < bot ) bot = vol[ii] ;
542 else if( vol[ii] > top ) top = vol[ii] ;
543 }
544 if( bot >= top ) return ;
545 #endif
546
547 /********************************/
548 /* 02 Feb 2001: include padding */
549
550 { float * vvv , *www ;
551 int nxp , nyp , nzp ;
552
553 THD_rota_envpad() ; /* 05 Feb 2001 */
554
555 nxp=nx+2*rotpx ; nyp=ny+2*rotpy ; nzp=nz+2*rotpz ;
556
557 if( rotpx > 0 && rotpy > 0 && rotpz > 0 )
558 vvv = EDIT_volpad_even( rotpx,rotpy,rotpz , nx,ny,nz , MRI_float,vol ) ;
559 else
560 vvv = vol ;
561
562 apply_3shear( shr , nxp,nyp,nzp , vvv ) ; /*-- do the actual rotation! --*/
563
564 if( vvv != vol ){
565 www = EDIT_volpad_even( -rotpx,-rotpy,-rotpz , nxp,nyp,nzp , MRI_float,vvv ) ;
566 free(vvv) ;
567 memcpy( vol , www , sizeof(float)*nx*ny*nz ) ; free(www) ;
568 }
569 }
570
571 /********************************/
572
573 #ifdef CLIPIT
574 for( ii=0 ; ii < nxyz ; ii++ ){
575 if( vol[ii] < bot ) vol[ii] = bot ;
576 else if( vol[ii] > top ) vol[ii] = top ;
577 }
578 #endif
579
580 return ;
581 }
582
583 /*------------------------------------------------------------------------------
584 14 Feb 2001:
585 Like the above, but with geometrical information about the volume
586 given from the image header
587 --------------------------------------------------------------------------------*/
588
THD_rota3D_matvec(MRI_IMAGE * im,THD_dmat33 rmat,THD_dfvec3 tvec)589 MRI_IMAGE * THD_rota3D_matvec( MRI_IMAGE * im, THD_dmat33 rmat,THD_dfvec3 tvec )
590 {
591 MRI_IMAGE * jm ;
592 float * jvol ;
593
594 if( ! MRI_IS_3D(im) ){
595 fprintf(stderr,"\n*** THD_rota3D_matvec: non-3D image input!\n") ;
596 return NULL ;
597 }
598
599 jm = mri_new_vol( im->nx , im->ny , im->nz , MRI_float ) ;
600 MRI_COPY_AUX(jm,im) ;
601 jvol = MRI_FLOAT_PTR(jm) ;
602
603 EDIT_coerce_type( im->nvox ,
604 im->kind , mri_data_pointer(im) , MRI_float , jvol ) ;
605
606 THD_rota_vol_matvec( im->nx , im->ny , im->nz ,
607 fabs(im->dx) , fabs(im->dy) , fabs(im->dz) , jvol ,
608 rmat , tvec ) ;
609 return jm ;
610 }
611