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 #include "thd_shear3d.h"
8
9 /*-- prototypes for funcs at end of file --*/
10
11 static void rotate_stdin_points( THD_dfvec3, THD_dmat33, int,THD_dfvec3 ) ;
12
13 /*------------------------------------------------------------------------------*/
14
15 #define MATVEC_DICOM 1
16 #define MATVEC_ORDER 2
17
18 static int verb=0 ;
19
20 /*-- 19 Jun 2001 stuff --*/
21
22 #define MODE_DFILE 1
23 #define MODE_1DFILE 2
24 MRI_IMAGE * get_dfile_params( char * fname , int mode ) ;
25
26 /*------------------------------------------------------------------------------*/
27
main(int argc,char * argv[])28 int main( int argc , char * argv[] )
29 {
30 THD_3dim_dataset * dset=NULL ;
31 char * new_prefix = "rota" , * cpt ;
32 float dx=0 , dy=0 , dz=0 ;
33 int ax1=0,ax2=1,ax3=2 , adx,ady,adz ;
34 char cdx='\0',cdy='\0',cdz='\0' ;
35 float th1=0.0,th2=0.0,th3=0.0 ;
36 int iopt , nvox , rotarg=-1 , dcode=-1 , ival,nval , ihand ;
37 float * fvol ;
38 double cputim=0.0 ;
39 int clipit=1 ; /* 11 Apr 2000 and 16 Apr 2002 */
40 float cbot=0.0,ctop=0.0 ;
41
42 int matvec=0 ; /* 19 July 2000 */
43 THD_dmat33 rmat , pp,ppt ;
44 THD_dfvec3 tvec ;
45
46 int dopoints=0 , doorigin=0 ; /* 21 Nov 2000 */
47 double xo=0.0 , yo=0.0 , zo=0.0 ;
48
49 THD_3dim_dataset *rotpar_dset=NULL , *gridpar_dset=NULL ; /* 07 Feb 2001 */
50
51 int zpad=0 ; /* 05 Feb 2001 */
52
53 char *dname=NULL ; /* 19 Jun 2001 */
54 int dmode=0 ;
55 MRI_IMAGE *dim=NULL ;
56 float *dar=NULL ;
57 int ndar=0 ;
58 int skipit=0 ;
59
60 /*-------------------------------*/
61
62 LOAD_DIAG_DMAT(rmat,1.0,1.0,1.0) ;
63 memset(&tvec, 0, sizeof(THD_dfvec3));
64 /*-- read command line arguments --*/
65
66 if( argc < 2 || strncmp(argv[1],"-help",4) == 0 ){
67 printf(
68 "Usage: 3drotate [options] dataset\n"
69 "Rotates and/or translates all bricks from an AFNI dataset.\n"
70 "'dataset' may contain a sub-brick selector list.\n"
71 "\n"
72 "GENERIC OPTIONS:\n"
73 " -prefix fname = Sets the output dataset prefix name to be 'fname'\n"
74 " -verbose = Prints out progress reports (to stderr)\n"
75 "\n"
76 "OPTIONS TO SPECIFY THE ROTATION/TRANSLATION:\n"
77 "-------------------------------------------\n"
78 "*** METHOD 1 = direct specification:\n"
79 "At most one of these shift options can be used:\n"
80 " -ashift dx dy dz = Shifts the dataset 'dx' mm in the x-direction, etc.,\n"
81 " AFTER rotation.\n"
82 " -bshift dx dy dz = Shifts the dataset 'dx' mm in the x-direction, etc.,\n"
83 " BEFORE rotation.\n"
84 " The shift distances by default are along the (x,y,z) axes of the dataset\n"
85 " storage directions (see the output of '3dinfo dataset'). To specify them\n"
86 " anatomically, you can suffix a distance with one of the symbols\n"
87 " 'R', 'L', 'A', 'P', 'I', and 'S', meaning 'Right', 'Left', 'Anterior',\n"
88 " 'Posterior', 'Inferior', and 'Superior', respectively.\n"
89 "\n"
90 " -rotate th1 th2 th3\n"
91 " Specifies the 3D rotation to be composed of 3 planar rotations:\n"
92 " 1) 'th1' degrees about the 1st axis, followed by\n"
93 " 2) 'th2' degrees about the (rotated) 2nd axis, followed by\n"
94 " 3) 'th3' degrees about the (doubly rotated) 3rd axis.\n"
95 " Which axes are used for these rotations is specified by placing\n"
96 " one of the symbols 'R', 'L', 'A', 'P', 'I', and 'S' at the end\n"
97 " of each angle (e.g., '10.7A'). These symbols denote rotation\n"
98 " about the 'Right-to-Left', 'Left-to-Right', 'Anterior-to-Posterior',\n"
99 " 'Posterior-to-Anterior', 'Inferior-to-Superior', and\n"
100 " 'Superior-to-Inferior' axes, respectively. A positive rotation is\n"
101 " defined by the right-hand rule.\n"
102 "\n"
103 "*** METHOD 2 = copy from output of 3dvolreg:\n"
104 " -rotparent rset\n"
105 " Specifies that the rotation and translation should be taken from the\n"
106 " first 3dvolreg transformation found in the header of dataset 'rset'.\n"
107 " -gridparent gset\n"
108 " Specifies that the output dataset of 3drotate should be shifted to\n"
109 " match the grid of dataset 'gset'. Can only be used with -rotparent.\n"
110 " This dataset should be one this is properly aligned with 'rset' when\n"
111 " overlaid in AFNI.\n"
112 " * If -rotparent is used, then don't use -matvec, -rotate, or -[ab]shift.\n"
113 " * If 'gset' has a different number of slices than the input dataset,\n"
114 " then the output dataset will be zero-padded in the slice direction\n"
115 " to match 'gset'.\n"
116 " * These options are intended to be used to align datasets between sessions:\n"
117 " S1 = SPGR from session 1 E1 = EPI from session 1\n"
118 " S2 = SPGR from session 2 E2 = EPI from session 2\n"
119 " 3dvolreg -twopass -twodup -base S1+orig -prefix S2reg S2+orig\n"
120 " 3drotate -rotparent S2reg+orig -gridparent E1+orig -prefix E2reg E2+orig\n"
121 " The result will have E2reg rotated from E2 in the same way that S2reg\n"
122 " was from S2, and also shifted/padded (as needed) to overlap with E1.\n"
123 "\n"
124 "*** METHOD 3 = give the transformation matrix/vector directly:\n"
125 " -matvec_dicom mfile\n"
126 " -matvec_order mfile\n"
127 " Specifies that the rotation and translation should be read from file\n"
128 " 'mfile', which should be in the format\n"
129 " u11 u12 u13 v1\n"
130 " u21 u22 u23 v2\n"
131 " u31 u32 u33 u3\n"
132 " where each 'uij' and 'vi' is a number. The 3x3 matrix [uij] is the\n"
133 " orthogonal matrix of the rotation, and the 3-vector [vi] is the -ashift\n"
134 " vector of the translation.\n"
135 "\n"
136 "*** METHOD 4 = copy the transformation from 3dTagalign:\n"
137 " -matvec_dset mset\n"
138 " Specifies that the rotation and translation should be read from\n"
139 " the .HEAD file of dataset 'mset', which was created by program\n"
140 " 3dTagalign.\n"
141 " * If -matvec_dicom is used, the matrix and vector are given in Dicom\n"
142 " coordinate order (+x=L, +y=P, +z=S). This is the option to use\n"
143 " if mfile is generated using 3dTagalign -matvec mfile.\n"
144 " * If -matvec_order is used, the matrix and vector are given in the\n"
145 " coordinate order of the dataset axes, whatever they may be.\n"
146 " * You can't mix -matvec_* options with -rotate and -*shift.\n"
147 "\n"
148 "*** METHOD 5 = input rotation+shift parameters from an ASCII file:\n"
149 " -dfile dname *OR* -1Dfile dname\n"
150 " With these methods, the movement parameters for each sub-brick\n"
151 " of the input dataset are read from the file 'dname'. This file\n"
152 " should consist of columns of numbers in ASCII format. Six (6)\n"
153 " numbers are read from each line of the input file. If the\n"
154 " '-dfile' option is used, each line of the input should be at\n"
155 " least 7 numbers, and be of the form\n"
156 " ignored roll pitch yaw dS dL dP\n"
157 " If the '-1Dfile' option is used, then each line of the input\n"
158 " should be at least 6 numbers, and be of the form\n"
159 " roll pitch yaw dS dL dP\n"
160 " (These are the forms output by the '-dfile' and\n"
161 " '-1Dfile' options of program 3dvolreg; see that\n"
162 " program's -help output for the hideous details.)\n"
163 " The n-th sub-brick of the input dataset will be transformed\n"
164 " using the parameters from the n-th line of the dname file.\n"
165 " If the dname file doesn't contain as many lines as the\n"
166 " input dataset has sub-bricks, then the last dname line will\n"
167 " be used for all subsequent sub-bricks. Excess columns or\n"
168 " rows will be ignored.\n"
169 " N.B.: Rotation is always about the center of the volume.\n"
170 " If the parameters are derived from a 3dvolreg run\n"
171 " on a dataset with a different center in xyz-space,\n"
172 " the results may not be what you want!\n"
173 " N.B.: You can't use -dfile/-1Dfile with -points (infra).\n"
174 "\n"
175 "POINTS OPTIONS (instead of datasets):\n"
176 "------------------------------------\n"
177 " -points\n"
178 " -origin xo yo zo\n"
179 " These options specify that instead of rotating a dataset, you will\n"
180 " be rotating a set of (x,y,z) points. The points are read from stdin.\n"
181 " * If -origin is given, the point (xo,yo,zo) is used as the center for\n"
182 " the rotation.\n"
183 " * If -origin is NOT given, and a dataset is given at the end of the\n"
184 " command line, then the center of the dataset brick is used as\n"
185 " (xo,yo,zo). The dataset will NOT be rotated if -points is given.\n"
186 " * If -origin is NOT given, and NO dataset is given at the end of the\n"
187 " command line, then xo=yo=zo=0 is assumed. You probably don't\n"
188 " want this.\n"
189 " * (x,y,z) points are read from stdin as 3 ASCII-formatted numbers per\n"
190 " line, as in 3dUndump. Any succeeding numbers on input lines will\n"
191 " be copied to the output, which will be written to stdout.\n"
192 " * The input (x,y,z) coordinates are taken in the same order as the\n"
193 " axes of the input dataset. If there is no input dataset, then\n"
194 " negative x = R positive x = L }\n"
195 " negative y = A positive y = P } e.g., the DICOM order\n"
196 " negative z = I positive z = S }\n"
197 " One way to dump some (x,y,z) coordinates from a dataset is:\n"
198 "\n"
199 " 3dmaskdump -mask something+tlrc -o xyzfilename -noijk\n"
200 " '3dcalc( -a dset+tlrc -expr x -datum float )'\n"
201 " '3dcalc( -a dset+tlrc -expr y -datum float )'\n"
202 " '3dcalc( -a dset+tlrc -expr z -datum float )'\n"
203 "\n"
204 " (All of this should be on one command line.)\n"
205 "============================================================================\n"
206 "\n"
207 "Example: 3drotate -prefix Elvis -bshift 10S 0 0 -rotate 30R 0 0 Sinatra+orig\n"
208 "\n"
209 "This will shift the input 10 mm in the superior direction, followed by a 30\n"
210 "degree rotation about the Right-to-Left axis (i.e., nod the head forward).\n"
211 "\n"
212 "============================================================================\n"
213 "Algorithm: The rotation+shift is decomposed into 4 1D shearing operations\n"
214 " (a 3D generalization of Paeth's algorithm). The interpolation\n"
215 " (i.e., resampling) method used for these shears can be controlled\n"
216 " by the following options:\n"
217 "\n"
218 " -Fourier = Use a Fourier method (the default: most accurate; slowest).\n"
219 " -NN = Use the nearest neighbor method.\n"
220 " -linear = Use linear (1st order polynomial) interpolation (least accurate).\n"
221 " -cubic = Use the cubic (3rd order) Lagrange polynomial method.\n"
222 " -quintic = Use the quintic (5th order) Lagrange polynomial method.\n"
223 " -heptic = Use the heptic (7th order) Lagrange polynomial method.\n"
224 "\n"
225 " -Fourier_nopad = Use the Fourier method WITHOUT padding\n"
226 " * If you don't mind - or even want - the wraparound effect\n"
227 " * Works best if dataset grid size is a power of 2, possibly\n"
228 " times powers of 3 and 5, in all directions being altered.\n"
229 " * The main use would seem to be to un-wraparound poorly\n"
230 " reconstructed images, by using a shift; for example:\n"
231 " 3drotate -ashift 30A 0 0 -Fourier_nopad -prefix Anew A+orig\n"
232 " * This option is also available in the Nudge Dataset plugin.\n"
233 "\n"
234 " -clipit = Clip results to input brick range [now the default].\n"
235 " -noclip = Don't clip results to input brick range.\n"
236 "\n"
237 " -zpad n = Zeropad around the edges by 'n' voxels during rotations\n"
238 " (these edge values will be stripped off in the output)\n"
239 " N.B.: Unlike to3d, in this program '-zpad' adds zeros in\n"
240 " all directions.\n"
241 " N.B.: The environment variable AFNI_ROTA_ZPAD can be used\n"
242 " to set a nonzero default value for this parameter.\n"
243 ) ;
244
245 printf("\n" MASTER_SHORTHELP_STRING ) ;
246 PRINT_COMPILE_DATE ; exit(0) ;
247 }
248
249 mainENTRY("3drotate main"); machdep(); PRINT_VERSION("3drotate"); AUTHOR("RW Cox");
250
251 /*-- 20 Apr 2001: addto the arglist, if user wants to [RWCox] --*/
252
253 { int new_argc ; char ** new_argv ;
254 addto_args( argc , argv , &new_argc , &new_argv ) ;
255 if( new_argv != NULL ){ argc = new_argc ; argv = new_argv ; }
256 }
257
258 AFNI_logger("3drotate",argc,argv) ;
259
260 #define ERREX(str) (fprintf(stderr,"*** %s\n",str),exit(1))
261
262 iopt = 1 ;
263 while( iopt < argc && argv[iopt][0] == '-' ){
264
265 if( strncmp(argv[iopt],"-zpad",5) == 0 ){ /* 05 Feb 2001 */
266 if( zpad > 0 )
267 fprintf(stderr,"+++ WARNING: second -zpad option!\n") ;
268 zpad = (int) strtod( argv[++iopt] , NULL ) ;
269 if( zpad < 0 ){
270 fprintf(stderr,"*** ERROR: Can't use -zpad %d\n",zpad) ;
271 exit(1) ;
272 }
273 THD_rota_setpad(zpad,zpad,zpad) ;
274 iopt++ ; continue ;
275 }
276
277 if( strncmp(argv[iopt],"-points",7) == 0 ){ /* 21 Nov 2000 */
278 dopoints = 1 ;
279 iopt++ ; continue ;
280 }
281
282 if( strncmp(argv[iopt],"-origin",7) == 0 ){ /* 21 Nov 2000 */
283 xo = strtod( argv[++iopt] , NULL ) ;
284 yo = strtod( argv[++iopt] , NULL ) ;
285 zo = strtod( argv[++iopt] , NULL ) ;
286 doorigin = 1 ;
287 iopt++ ; continue ;
288 }
289
290 if( strncmp(argv[iopt],"-rotpar",7) == 0 ){ /* 07 Feb 2001 */
291
292 ATR_float *atr ;
293
294 if( rotpar_dset != NULL )
295 ERREX("*** Can't use 2 -rotparent options!") ;
296 if( matvec )
297 ERREX("*** Can't combine -rotparent with -matvec!") ;
298 if( dcode > 0 || rotarg > 0 )
299 ERREX("*** Can't use -rotparent with -shift or -rotate options!") ;
300 if( dname != NULL )
301 ERREX("*** Can't use -rotparent with -dfile or -1Dfile options!") ;
302
303 rotpar_dset = THD_open_one_dataset( argv[++iopt] ) ;
304 if( rotpar_dset == NULL )
305 ERREX("*** Can't open -rotparent dataset!\n") ;
306 THD_make_cardinal(rotpar_dset); /* deoblique 21 Oct, 2011 [rickr] */
307
308 atr = THD_find_float_atr( rotpar_dset->dblk , "VOLREG_MATVEC_000000" ) ;
309 if( atr == NULL || atr->nfl < 12 )
310 ERREX("*** -rotparent dataset doesn't have VOLREG attributes!?") ;
311
312 iopt++ ; continue ;
313 }
314
315 if( strncmp(argv[iopt],"-gridpar",7) == 0 ){ /* 07 Feb 2001 */
316
317 if( gridpar_dset != NULL )
318 ERREX("*** Can't use -2 -gridparent options!") ;
319
320 gridpar_dset = THD_open_one_dataset( argv[++iopt] ) ;
321 if( gridpar_dset == NULL )
322 ERREX("*** Can't open -gridparent dataset!\n") ;
323 THD_make_cardinal(gridpar_dset); /* deoblique 21 Oct, 2011 [rickr] */
324
325 iopt++ ; continue ;
326 }
327
328 if( strncmp(argv[iopt],"-matvec_",8) == 0 ){ /* 19 Jul 2000 */
329
330 MRI_IMAGE *matim ; float *matar , sum ;
331
332 if( matvec )
333 ERREX("*** Can't use 2 -matvec options!") ;
334 if( dcode > 0 || rotarg > 0 )
335 ERREX("*** Can't use -matvec with -shift or -rotate options!") ;
336 if( rotpar_dset != NULL )
337 ERREX("*** Can't use -matvec with -rotparent option!") ;
338 if( dname != NULL )
339 ERREX("*** Can't use -matvec with -dfile or -1Dfile options!") ;
340
341 if( strcmp(argv[iopt],"-matvec_order") == 0 ) matvec = MATVEC_ORDER ;
342 else matvec = MATVEC_DICOM ;
343
344 if( strcmp(argv[iopt],"-matvec_dset") == 0){ /* 20 July 2000 */
345 THD_3dim_dataset * mvset ; ATR_float * atr ;
346
347 mvset = THD_open_dataset( argv[++iopt] ) ;
348 if( mvset == NULL ) ERREX("*** Can't read -matvec_dset dataset!") ;
349 THD_make_cardinal(mvset); /* deoblique 21 Oct, 2011 [rickr] */
350 atr = THD_find_float_atr( mvset->dblk , "TAGALIGN_MATVEC" ) ;
351 if( atr == NULL || atr->nfl < 12 )
352 ERREX("*** -matvec_dset doesn't have matrix+vector in .HEAD!") ;
353 matar = atr->fl ;
354 LOAD_DMAT(rmat,matar[0],matar[1],matar[2],
355 matar[4],matar[5],matar[6],
356 matar[8],matar[9],matar[10] ) ;
357 LOAD_DFVEC3(tvec,matar[3],matar[7],matar[11]) ;
358 DSET_delete(mvset) ;
359 } else {
360 matim = mri_read_ascii( argv[++iopt] ) ;
361 if( matim == NULL ) ERREX("Can't read -matvec file!") ;
362 if( matim->nx != 4 || matim->ny != 3 ) ERREX("-matvec file not 4x3!") ;
363
364 matar = MRI_FLOAT_PTR(matim) ;
365 LOAD_DMAT(rmat,matar[0],matar[1],matar[2],
366 matar[4],matar[5],matar[6],
367 matar[8],matar[9],matar[10] ) ;
368 LOAD_DFVEC3(tvec,matar[3],matar[7],matar[11]) ;
369
370 mri_free(matim) ;
371 }
372
373 /* check if matrix is approximately orthogonal */
374 /* [will be orthogonalized in rot_to_shear_matvec() in thd_shear3d.c] */
375
376 pp = TRANSPOSE_DMAT(rmat) ; pp = DMAT_MUL(pp,rmat) ;
377 sum = fabs(pp.mat[0][0]-1.0)+fabs(pp.mat[1][0]) +fabs(pp.mat[2][0])
378 +fabs(pp.mat[0][1]) +fabs(pp.mat[1][1]-1.0)+fabs(pp.mat[2][1])
379 +fabs(pp.mat[0][2]) +fabs(pp.mat[1][2]) +fabs(pp.mat[2][2]-1.0);
380 if( sum > 0.01 ) ERREX("-matvec matrix not orthogonal!") ;
381
382 iopt++ ; continue ;
383 }
384
385 if( strncmp(argv[iopt],"-clipit",4) == 0 ){ /* 11 Apr 2000 */
386 fprintf(stderr,"++ Notice: -clipit is now the default\n") ;
387 clipit = 1 ;
388 iopt++ ; continue ;
389 }
390
391 if( strncmp(argv[iopt],"-noclip",4) == 0 ){ /* 16 Apr 2002 */
392 clipit = 0 ;
393 iopt++ ; continue ;
394 }
395
396 if( strncmp(argv[iopt],"-prefix",4) == 0 ){
397 new_prefix = argv[++iopt] ;
398 iopt++ ; continue ;
399 }
400
401 if( strncmp(argv[iopt],"-verbose",5) == 0 ){
402 verb = 1 ;
403 iopt++ ; continue ;
404 }
405
406 if( strcmp(argv[iopt],"-Fourier_nopad") == 0 ){ /* 13 May 2003 */
407 THD_rota_method( MRI_FOURIER_NOPAD ) ;
408 iopt++ ; continue ;
409 }
410 if( strncmp(argv[iopt],"-Fourier",4) == 0 || strncmp(argv[iopt],"-fourier",4) == 0 ){
411 THD_rota_method( MRI_FOURIER ) ;
412 iopt++ ; continue ;
413 }
414
415 if( strncmp(argv[iopt],"-cubic",4) == 0 || strncmp(argv[iopt],"-Cubic",4) == 0 ){
416 THD_rota_method( MRI_CUBIC ) ;
417 iopt++ ; continue ;
418 }
419
420 if( strncmp(argv[iopt],"-quintic",4) == 0 || strncmp(argv[iopt],"-Quintic",4) == 0 ){
421 THD_rota_method( MRI_QUINTIC ) ;
422 iopt++ ; continue ;
423 }
424
425 if( strncmp(argv[iopt],"-heptic",4) == 0 || strncmp(argv[iopt],"-Heptic",4) == 0 ){
426 THD_rota_method( MRI_HEPTIC ) ;
427 iopt++ ; continue ;
428 }
429
430 if( strncmp(argv[iopt],"-linear",4) == 0 || strncmp(argv[iopt],"-Linear",4) == 0 ){
431 THD_rota_method( MRI_LINEAR ) ; clipit = 0 ;
432 iopt++ ; continue ;
433 }
434
435 if( strncmp(argv[iopt],"-nn",3) == 0 || strncmp(argv[iopt],"-NN",4) == 0 ){
436 THD_rota_method( MRI_NN ) ; clipit = 0 ;
437 iopt++ ; continue ;
438 }
439
440 if( strncmp(argv[iopt],"-ashift",4) == 0 ){
441 if( matvec ) ERREX("*** Can't use -ashift with -matvec!") ;
442 if( dcode > 0 ){fprintf(stderr,"*** Can't use 2 shift options!\n");exit(1);}
443 if( rotpar_dset != NULL )
444 ERREX("*** Can't use -ashift with -rotparent!") ;
445 if( dname != NULL )
446 ERREX("*** Can't use -ashift with -dfile or -1Dfile options!") ;
447 dx = strtod( argv[++iopt] , &cpt ) ; cdx = *cpt ;
448 dy = strtod( argv[++iopt] , &cpt ) ; cdy = *cpt ;
449 dz = strtod( argv[++iopt] , &cpt ) ; cdz = *cpt ;
450 dcode = DELTA_AFTER ;
451 iopt++ ; continue ;
452 }
453
454 if( strncmp(argv[iopt],"-bshift",4) == 0 ){
455 if( matvec ) ERREX("*** Can't use -bshift with -matvec!") ;
456 if( dcode > 0 ){fprintf(stderr,"*** Can't use 2 shift options!\n");exit(1);}
457 if( rotpar_dset != NULL )
458 ERREX("*** Can't use -bshift with -rotparent!") ;
459 if( dname != NULL )
460 ERREX("*** Can't use -bshift with -dfile or -1Dfile options!") ;
461 dx = strtod( argv[++iopt] , &cpt ) ; cdx = *cpt ;
462 dy = strtod( argv[++iopt] , &cpt ) ; cdy = *cpt ;
463 dz = strtod( argv[++iopt] , &cpt ) ; cdz = *cpt ;
464 dcode = DELTA_BEFORE ;
465 iopt++ ; continue ;
466 }
467
468 #if 0
469 if( strncmp(argv[iopt],"-cshift",4) == 0 ){
470 if( dcode > 0 ){fprintf(stderr,"*** Can't use 2 shift options!\n");exit(1);}
471 dx = strtod( argv[++iopt] , &cpt ) ; cdx = *cpt ;
472 dy = strtod( argv[++iopt] , &cpt ) ; cdy = *cpt ;
473 dz = strtod( argv[++iopt] , &cpt ) ; cdz = *cpt ;
474 dcode = DELTA_FIXED ;
475 iopt++ ; continue ;
476 }
477 #endif
478
479 if( strncmp(argv[iopt],"-rotate",4) == 0 ){
480 if( matvec ) ERREX("*** Can't use -rotate with -matvec!") ;
481 if( rotarg > 0 ) ERREX("*** Can't have 2 -rotate options!\n") ;
482 if( rotpar_dset != NULL )
483 ERREX("*** Can't use -rotate with -rotparent!") ;
484 if( dname != NULL )
485 ERREX("*** Can't use -rotate with -dfile or -1Dfile options!") ;
486 rotarg = iopt ; /* save and process later */
487 iopt += 4 ; continue ;
488 }
489
490 /*-- 19 Jun 2001 --*/
491
492 if( strcmp(argv[iopt],"-dfile")==0 || strcmp(argv[iopt],"-1Dfile")==0 ){
493
494 if( dname != NULL )
495 ERREX("*** Can't use 2 -dfile/-1Dfile options!") ;
496 if( matvec )
497 ERREX("*** Can't combine -dfile/-1Dfile with -matvec!") ;
498 if( dcode > 0 || rotarg > 0 )
499 ERREX("*** Can't use -dfile/-1Dfile with -shift or -rotate options!") ;
500 if( rotpar_dset != NULL )
501 ERREX("*** Can't use -dfile/-1Dfile with -rotparent!") ;
502
503 if( strcmp(argv[iopt],"-dfile") ==0 ) dmode = MODE_DFILE ;
504 else if( strcmp(argv[iopt],"-1Dfile")==0 ) dmode = MODE_1DFILE ;
505
506 dname = argv[++iopt] ;
507
508 dim = get_dfile_params( dname , dmode ) ;
509 if( dim == NULL )
510 ERREX("*** Error with -dfile or -1Dfile!") ;
511
512 dar = MRI_FLOAT_PTR(dim) ; /* 6 x ndar array */
513 ndar = dim->ny ;
514 iopt++ ; continue ;
515 }
516
517 fprintf(stderr,"*** Unknown option: %s\n",argv[iopt]) ; exit(1) ;
518 }
519
520 /*-- check for legal combinations --*/
521
522 if( gridpar_dset != NULL && rotpar_dset == NULL ){
523 fprintf(stderr,"+++ WARNING: -gridparent means nothing without -rotparent!\n");
524 DSET_delete( gridpar_dset ) ;
525 gridpar_dset = NULL ;
526 }
527
528 if( gridpar_dset != NULL && dopoints ){
529 fprintf(stderr,"+++ WARNING: -gridparent means nothing with -points!\n") ;
530 DSET_delete( gridpar_dset ) ;
531 gridpar_dset = NULL ;
532 }
533
534 if( doorigin && !dopoints ){
535 fprintf(stderr,"+++ WARNING: -origin means nothing without -points!\n") ;
536 }
537
538 if( dar != NULL && dopoints )
539 ERREX("*** Can't combine -dfile/-1Dfile with -points!") ;
540
541 if( matvec==0 && dcode<0 && rotarg<0 && rotpar_dset==NULL && dar==NULL )
542 ERREX("Don't you want to do anything [no -rotate,-shift,-matvec,-rotparent,-dfile]?");
543
544 /** read input dataset */
545
546 if( iopt >= argc && !dopoints ) ERREX("*** No input dataset?") ;
547
548 if( iopt < argc ){
549 dset = THD_open_dataset( argv[iopt] ) ;
550 if( dset == NULL ){
551 fprintf(stderr,"*** Cannot open dataset %s!\n",argv[iopt]); exit(1);
552 }
553 THD_make_cardinal(dset); /* deoblique 21 Oct, 2011 [rickr] */
554 } else {
555 dset = EDIT_empty_copy(NULL) ; /* 21 Nov 2000: need a fake dataset */
556 if( !doorigin )
557 fprintf(stderr,"+++ WARNING: no -origin and no input dataset on command line!\n") ;
558 }
559
560 if( dopoints && !doorigin ){
561 xo = dset->daxes->xxorg + 0.5*(dset->daxes->nxx - 1)*dset->daxes->xxdel ;
562 yo = dset->daxes->yyorg + 0.5*(dset->daxes->nyy - 1)*dset->daxes->yydel ;
563 zo = dset->daxes->zzorg + 0.5*(dset->daxes->nzz - 1)*dset->daxes->zzdel ;
564 fprintf(stderr,"+++ Using -origin %g %g %g\n",xo,yo,zo) ;
565 }
566
567 /* now can process rotation arguments */
568
569 ihand = THD_handedness(dset) ;
570
571 if( rotarg > 0 ){
572 int neg ;
573 iopt = rotarg ;
574
575 th1 = (PI/180.0) * strtod( argv[++iopt] , &cpt ) ;
576 switch( *cpt ){
577 default: ERROR_exit("Illegal code after th1 in -rotate\n");
578
579 case '\0': case 'x': case 'X': ax1 = 0 ; neg = 0 ; break ;
580 case 'y': case 'Y': ax1 = 1 ; neg = 0 ; break ;
581 case 'z': case 'Z': ax1 = 2 ; neg = 0 ; break ;
582
583 case 'A': case 'P':
584 case 'R': case 'L':
585 case 'I': case 'S': ax1 = THD_axcode(dset,*cpt) ;
586 neg = (ax1 < 0) ;
587 ax1 = abs(ax1) - 1 ; break ;
588 }
589 if( neg ) th1 = -th1 ;
590
591 th2 = (PI/180.0) * strtod( argv[++iopt] , &cpt ) ;
592 switch( *cpt ){
593 default: ERROR_exit("Illegal code after th2 in -rotate\n");
594
595 case 'x': case 'X': ax2 = 0 ; break ;
596 case '\0': case 'y': case 'Y': ax2 = 1 ; break ;
597 case 'z': case 'Z': ax2 = 2 ; break ;
598
599 case 'A': case 'P':
600 case 'R': case 'L':
601 case 'I': case 'S': ax2 = THD_axcode(dset,*cpt) ;
602 neg = (ax2 < 0) ;
603 ax2 = abs(ax2) - 1 ; break ;
604 }
605 if( neg ) th2 = -th2 ;
606
607 th3 = (PI/180.0) * strtod( argv[++iopt] , &cpt ) ;
608 switch( *cpt ){
609 default: ERROR_exit("Illegal code after th3 in -rotate\n");
610
611 case 'x': case 'X': ax3 = 0 ; break ;
612 case 'y': case 'Y': ax3 = 1 ; break ;
613 case '\0': case 'z': case 'Z': ax3 = 2 ; break ;
614
615 case 'A': case 'P':
616 case 'R': case 'L':
617 case 'I': case 'S': ax3 = THD_axcode(dset,*cpt) ;
618 neg = (ax3 < 0) ;
619 ax3 = abs(ax3) - 1 ; break ;
620 }
621 if( neg ) th3 = -th3 ;
622
623 if( th1 == 0.0 && th2 == 0.0 && th3 == 0.0 )
624 INFO_message("All angles after -rotate are 0!") ;
625
626 if( ax1 < 0 || ax1 > 2 || ax2 < 0 || ax2 > 2 || ax3 < 0 || ax3 > 2 )
627 ERROR_exit("Can't understand axes codes in -rotate!") ;
628
629 if( ihand < 0 ){ th1 = -th1 ; th2 = -th2 ; th3 = -th3 ; }
630
631 #if 0
632 fprintf(stderr,"ihand=%d th1=%g th2=%g th3=%g\n",ihand,th1,th2,th3);
633 fprintf(stderr,"ax1=%d ax2=%d ax3=%d\n",ax1,ax2,ax3) ;
634 #endif
635
636 /* notice the minus signs on the angles ------+ */
637 /* | | | */
638 if( dopoints ) /* V V V */
639 rmat = rot_to_matrix( ax1,-th1,ax2,-th2,ax3,-th3 ) ; /* 21 Nov 2000 */
640 }
641
642 /* may need to process shift arguments as well */
643
644 if( dcode > 0 && (cdx != '\0' || cdy != '\0' || cdz != '\0') ){
645 float qdx=0 , qdy=0 , qdz=0 ;
646
647 adx = THD_axcode(dset,cdx) ;
648 if( adx > -99 || dx != 0.0 ){ /* 29 Jan 2001: skip if purely 0 */
649 switch( adx ){
650 case 1: qdx = -dx ; break ;
651 default:
652 case -1: qdx = dx ; break ;
653 case 2: qdy = -dx ; break ;
654 case -2: qdy = dx ; break ;
655 case 3: qdz = -dx ; break ;
656 case -3: qdz = dx ; break ;
657 }
658 }
659
660 ady = THD_axcode(dset,cdy) ;
661 if( ady > -99 || dy != 0.0 ){ /* 29 Jan 2001 */
662 switch( ady ){
663 case 1: qdx = -dy ; break ;
664 case -1: qdx = dy ; break ;
665 case 2: qdy = -dy ; break ;
666 default:
667 case -2: qdy = dy ; break ;
668 case 3: qdz = -dy ; break ;
669 case -3: qdz = dy ; break ;
670 }
671 }
672
673 adz = THD_axcode(dset,cdz) ;
674 if( adz > -99 || dz != 0.0 ){ /* 29 Jan 2001 */
675 switch( adz ){
676 case 1: qdx = -dz ; break ;
677 case -1: qdx = dz ; break ;
678 case 2: qdy = -dz ; break ;
679 case -2: qdy = dz ; break ;
680 case 3: qdz = -dz ; break ;
681 default:
682 case -3: qdz = dz ; break ;
683 }
684 }
685
686 if( verb )
687 fprintf(stderr,"++ Shifting parameters:\n"
688 "++ direction codes: adx=%d ady=%d adz=%d\n"
689 "++ input values: dx=%g dy=%g dz=%g\n"
690 "++ output values: qdx=%g qdy=%g qdz=%g\n" ,
691 adx,ady,adz , dx,dy,dz , qdx,qdy,qdz ) ;
692
693 dx = qdx ; dy = qdy ; dz = qdz ;
694 }
695
696 /*- 19 July 2000: now can deal with -matvec_dicom case -*/
697
698 if( matvec == MATVEC_DICOM ){ /* convert matrix/vector */
699 pp = DBLE_mat_to_dicomm( dset ) ; /* to dataset coord order */
700 ppt = TRANSPOSE_DMAT(pp); /* from the DICOM order! */
701 rmat = DMAT_MUL(ppt,rmat); rmat = DMAT_MUL(rmat,pp);
702 tvec = DMATVEC(ppt,tvec);
703 }
704
705 /*-- 07 Feb 2001: deal with -rotparent and -gridparent case
706 [very similar to -matvec_dicom, actually] --*/
707
708 if( gridpar_dset != NULL ){
709 int mm , nz_gp , nz_ds ;
710
711 /* check for compatibility! */
712
713 mm = THD_dataset_mismatch( gridpar_dset , dset ) ;
714 if( mm & (MISMATCH_DELTA | MISMATCH_ORIENT) ){
715 fprintf(stderr,"*** Fatal Error:\n"
716 "*** -gridparent dataset and input dataset don't\n"
717 "*** match in grid spacing and/or orientation!\n" ) ;
718 exit(1) ;
719 }
720
721 if( DSET_NX(gridpar_dset) != DSET_NX(dset) ||
722 DSET_NY(gridpar_dset) != DSET_NY(dset) ){
723
724 fprintf(stderr,"*** Fatal Error:\n"
725 "*** -gridparent and input datasets\n"
726 "*** don't match in x,y dimensions!\n" ) ;
727 exit(1) ;
728 }
729
730 /* check for zero padding */
731
732 nz_gp = DSET_NZ(gridpar_dset) ; nz_ds = DSET_NZ(dset) ;
733
734 if( nz_gp < nz_ds ){
735 fprintf(stderr,"*** Fatal Error:\n"
736 "*** -gridparent has fewer slices than input dataset!\n") ;
737 exit(1) ;
738 }
739
740 if( nz_gp > nz_ds ){
741 int npad1 = (nz_gp - nz_ds) / 2 ;
742 int npad2 = (nz_gp - nz_ds) - npad1 ;
743 int add_I=0, add_S=0, add_A=0, add_P=0, add_L=0, add_R=0 ;
744 THD_3dim_dataset * pset ;
745 char *sp1=NULL,*sp2=NULL ;
746
747 /* where to add slices? and how many? */
748
749 switch( dset->daxes->zzorient ){
750 case ORI_R2L_TYPE:
751 case ORI_L2R_TYPE: add_R=npad1; add_L=npad2; sp1="R"; sp2="L"; break;
752
753 case ORI_P2A_TYPE:
754 case ORI_A2P_TYPE: add_A=npad1; add_P=npad2; sp1="A"; sp2="P"; break;
755
756 case ORI_I2S_TYPE:
757 case ORI_S2I_TYPE: add_I=npad1; add_S=npad2; sp1="I"; sp2="S"; break;
758 }
759
760 /* add them on */
761
762 if( verb )
763 fprintf(stderr,"+++ Zero padding to match -gridparent: -%s %d -%s %d\n",
764 sp1,npad1,sp2,npad2 ) ;
765
766 pset = THD_zeropad( dset,
767 add_I,add_S,add_A,add_P,add_L,add_R,
768 "Elvis" , ZPAD_PURGE ) ;
769
770 if( pset == NULL ){
771 fprintf(stderr,"*** Fatal Error:\n"
772 "*** Can't properly zeropad input dataset!\n" ) ;
773 exit(1) ;
774 }
775
776 /* toss input datset, replace with padded one */
777
778 DSET_delete(dset) ; dset = pset ;
779 }
780 }
781
782 if( rotpar_dset != NULL ){ /* compute transformation from -rotparent */
783 ATR_float * atr ;
784 float * matar , sum ;
785 THD_fvec3 fv ;
786 THD_dfvec3 dv,ev,qv , cv_e2, cv_e1, cv_s1, cv_s2 ;
787
788 /* load transformation from rotparent */
789
790 atr = THD_find_float_atr( rotpar_dset->dblk , "VOLREG_MATVEC_000000" ) ;
791 matar = atr->fl ;
792 LOAD_DMAT(rmat,matar[0],matar[1],matar[2],
793 matar[4],matar[5],matar[6],
794 matar[8],matar[9],matar[10] ) ;
795 LOAD_DFVEC3(tvec,matar[3],matar[7],matar[11]) ;
796
797 /* check if matrix is orthogonal */
798
799 pp = TRANSPOSE_DMAT(rmat) ; pp = DMAT_MUL(pp,rmat) ;
800 sum = fabs(pp.mat[0][0]-1.0)+fabs(pp.mat[1][0]) +fabs(pp.mat[2][0])
801 +fabs(pp.mat[0][1]) +fabs(pp.mat[1][1]-1.0)+fabs(pp.mat[2][1])
802 +fabs(pp.mat[0][2]) +fabs(pp.mat[1][2]) +fabs(pp.mat[2][2]-1.0);
803 if( sum > 0.01 ) ERREX("-rotparent matrix not orthogonal!") ;
804
805 /* must alter shift [tvec] to allow for differing
806 coordinates in the rotparent, gridparent, and input datasets */
807
808 /* cv_e2 = center of input dataset [Dicom coordinates] */
809
810 fv = THD_dataset_center( dset ) ; /* dataset coords */
811 FVEC3_TO_DFVEC3( fv , cv_e2 ) ; /* convert to double */
812
813 /* cv_e1 = center of gridparent */
814
815 if( gridpar_dset != NULL ){
816 fv = THD_dataset_center( gridpar_dset ) ;
817 FVEC3_TO_DFVEC3( fv , cv_e1 ) ;
818 } else {
819 cv_e1 = cv_e2 ; /* what else to do? */
820 }
821
822 /* cv_s2 = center of rotation in rotparent */
823
824 atr = THD_find_float_atr( rotpar_dset->dblk , "VOLREG_CENTER_OLD" ) ;
825 LOAD_DFVEC3( cv_s2 , atr->fl[0] , atr->fl[1] , atr->fl[2] ) ;
826
827 /* cv_s1 = center of base dataset for rotparent */
828
829 atr = THD_find_float_atr( rotpar_dset->dblk , "VOLREG_CENTER_BASE" ) ;
830 LOAD_DFVEC3( cv_s1 , atr->fl[0] , atr->fl[1] , atr->fl[2] ) ;
831
832 /* compute extra shift due to difference in
833 center of rotation between rotparent and input dataset,
834 then add in shifts caused by -twodup for rotparent and input */
835
836 dv = SUB_DFVEC3( cv_e2 , cv_s2 ) ;
837 ev = DMATVEC( rmat , dv ) ; /* R[E2-S2] */
838
839 dv = ev ; /* vestige of a stupid bug, since fixed */
840
841 ev = SUB_DFVEC3( cv_e1 , cv_s1 ) ; /* E1-S1 */
842
843 qv = SUB_DFVEC3( dv , ev ) ; /* R[E2-S2] + S1-E1 */
844
845 tvec = ADD_DFVEC3( tvec , qv ) ; /* shifted translation vector */
846
847 /* convert transformation from Dicom to dataset coords */
848
849 pp = DBLE_mat_to_dicomm( dset ) ;
850 ppt = TRANSPOSE_DMAT(pp);
851 rmat = DMAT_MUL(ppt,rmat); rmat = DMAT_MUL(rmat,pp); tvec = DMATVEC(ppt,tvec);
852
853 /* modify origin of output dataset to match -gridparent */
854
855 if( gridpar_dset != NULL ){
856 dset->daxes->xxorg = gridpar_dset->daxes->xxorg ;
857 dset->daxes->yyorg = gridpar_dset->daxes->yyorg ;
858 dset->daxes->zzorg = gridpar_dset->daxes->zzorg ;
859
860 /* 12 Feb 2001: adjust origin of time-offsets as well */
861
862 if( dset->taxis != NULL && dset->taxis->nsl > 0 ){
863 dset->taxis->zorg_sl = dset->daxes->zzorg ;
864 }
865 }
866
867 matvec = MATVEC_ORDER ; /* flag that transform comes from rmat/tvec */
868 }
869
870 /*-- 21 Nov 2000: read (x,y,z) points from stdin, process them, quit --*/
871
872 if( dopoints ){
873 THD_dfvec3 xyzorg ;
874 if( !matvec ) LOAD_DFVEC3( tvec , dx,dy,dz ) ;
875 if( dcode < 0 ) dcode = DELTA_AFTER ;
876 LOAD_DFVEC3(xyzorg,xo,yo,zo) ;
877 rotate_stdin_points( xyzorg , rmat , dcode,tvec ) ; /* at end of file */
878 exit(0) ;
879 }
880
881 /*-- 12 Feb 2001: adjust time-offsets for slice direction shifts --*/
882
883 if( dset->taxis != NULL && dset->taxis->nsl > 0 ){
884 int ndz ;
885 int kk,jj , nsl = dset->taxis->nsl ;
886
887 if( matvec )
888 ndz = (int) rint( tvec.xyz[2] / fabs(dset->daxes->zzdel) ) ; /* shift */
889 else
890 ndz = (int) rint( dz / fabs(dset->daxes->xxdel) ) ;
891
892 if( ndz != 0 ){
893 float * tsl = (float *)malloc(sizeof(float)*nsl) ;
894 for( kk=0 ; kk < nsl ; kk ++ ){
895 jj = kk - ndz ;
896 if( jj < 0 || jj >= nsl ) tsl[kk] = 0.0 ;
897 else tsl[kk] = dset->taxis->toff_sl[jj] ;
898 }
899 EDIT_dset_items( dset , ADN_toff_sl , tsl , ADN_none ) ;
900 free(tsl) ;
901 if( verb )
902 fprintf(stderr,"+++ adjusting time-offsets by %d slices\n",ndz) ;
903 }
904 }
905
906 /*- read dataset, prepare to process it, write back out (with new name) -*/
907
908 DSET_mallocize(dset) ;
909 DSET_load(dset) ;
910 dset->dblk->diskptr->storage_mode = STORAGE_BY_BRICK ; /* 14 Jan 2004 */
911
912 dset->idcode = MCW_new_idcode() ; /* 08 Jun 1999 - is a new dataset */
913 EDIT_dset_items( dset ,
914 ADN_prefix , new_prefix ,
915 ADN_label1 , new_prefix ,
916 ADN_none ) ;
917 if( THD_deathcon() && THD_is_file(dset->dblk->diskptr->header_name) ){
918 fprintf(stderr,
919 "** ERROR: Output file %s already exists -- cannot continue!\n",
920 dset->dblk->diskptr->header_name ) ;
921 exit(1) ;
922 }
923
924 /* old history is already in the dataset */
925
926 tross_Make_History( "3drotate" , argc,argv , dset ) ;
927
928 nvox = DSET_NVOX(dset) ;
929 fvol = (float *) malloc( sizeof(float) * nvox ) ;
930
931 nval = DSET_NVALS(dset) ;
932 if( verb ){
933 fprintf(stderr,"+++ %d sub-bricks: ",nval) ;
934 cputim = COX_cpu_time() ;
935 }
936
937 /* 03 May 2005: save rotation center */
938
939 { THD_fvec3 cv ;
940 cv = THD_dataset_center( dset ) ;
941 THD_set_float_atr( dset->dblk , "ROTATE_CENTER_OLD" , 3 , cv.xyz ) ;
942 THD_set_float_atr( dset->dblk , "ROTATE_CENTER_BASE" , 3 , cv.xyz ) ;
943 }
944
945 /*-- loop over all sub-bricks: copy into fvol, rotate fvol, copy back --*/
946
947 for( ival=0 ; ival < nval ; ival++ ){
948
949 if( verb ) fprintf(stderr,"%d",ival) ;
950
951 EDIT_coerce_type( nvox ,
952 DSET_BRICK_TYPE(dset,ival),DSET_ARRAY(dset,ival) ,
953 MRI_float,fvol ) ;
954
955 if( verb ) fprintf(stderr,".") ;
956
957 if( clipit ){ /* 11 Apr 2000 */
958 register int ii ; register float bb,tt ;
959 bb = tt = fvol[0] ;
960 for( ii=1 ; ii < nvox ; ii++ ){
961 if( fvol[ii] < bb ) bb = fvol[ii] ;
962 else if( fvol[ii] > tt ) tt = fvol[ii] ;
963 }
964 cbot = bb ; ctop = tt ;
965 }
966
967 /* 19 Jun 2001: get matrix/vector from rot params in file */
968
969 skipit = 0 ;
970 if( dar != NULL ){
971 THD_dvecmat dvm ; char rotcom[256] ; int jj=ival ;
972 static int ndar_over=0 ;
973
974 if( jj >= ndar ){
975 jj = ndar-1 ;
976 if( ndar_over == 0 )
977 WARNING_message("from brick %d on, using last line (%d) from %s\n",
978 jj , ndar-1 , dname ) ;
979 ndar_over++ ;
980 }
981
982 sprintf(rotcom,"-rotate %.4fI %.4fR %.4fA -ashift %.4fS %.4fL %.4fP",
983 dar[0+6*jj] , dar[1+6*jj] , dar[2+6*jj] ,
984 dar[3+6*jj] , dar[4+6*jj] , dar[5+6*jj] ) ;
985
986 dvm = THD_rotcom_to_matvec( dset , rotcom ) ; /* thd_rotangles.c */
987
988 rmat = dvm.mm ; tvec = dvm.vv ; matvec = MATVEC_ORDER ;
989
990 skipit = (dar[0+6*jj]==0.0 && dar[1+6*jj]==0.0 && dar[2+6*jj]==0.0 &&
991 dar[3+6*jj]==0.0 && dar[4+6*jj]==0.0 && dar[5+6*jj]==0.0 );
992
993 if( !skipit ){
994 skipit = ( fabs(rmat.mat[0][0]-1.0) < 0.00001 ) &&
995 ( fabs(rmat.mat[1][1]-1.0) < 0.00001 ) &&
996 ( fabs(rmat.mat[2][2]-1.0) < 0.00001 ) &&
997 ( fabs(tvec.xyz[0]) < 0.001 ) &&
998 ( fabs(tvec.xyz[1]) < 0.001 ) &&
999 ( fabs(tvec.xyz[2]) < 0.001 ) ;
1000 }
1001
1002 if( skipit ) fprintf(stderr,"[Matrix near identity, skipping]");
1003 }
1004
1005 /** carry out the rotation **/
1006
1007 if( !skipit ){
1008 if( matvec ){
1009 THD_rota_vol_matvec( DSET_NX(dset) , DSET_NY(dset) , DSET_NZ(dset) ,
1010 fabs(DSET_DX(dset)) ,
1011 fabs(DSET_DY(dset)) ,
1012 fabs(DSET_DZ(dset)) ,
1013 fvol , rmat , tvec ) ;
1014 } else {
1015 THD_rota_vol( DSET_NX(dset) , DSET_NY(dset) , DSET_NZ(dset) ,
1016 fabs(DSET_DX(dset)),
1017 fabs(DSET_DY(dset)),
1018 fabs(DSET_DZ(dset)), fvol ,
1019 ax1,th1, ax2,th2, ax3,th3, dcode,dx,dy,dz ) ;
1020
1021 /* 02 May 2005: save matrix+vector for recording below */
1022
1023 rmat = rot_to_matrix( ax1,-th1,ax2,-th2,ax3,-th3 ) ;
1024 LOAD_DFVEC3( tvec , dx,dy,dz ) ;
1025 if( dcode == DELTA_BEFORE ) tvec = DMATVEC(rmat,tvec) ;
1026 }
1027 }
1028
1029 /** 02 May 2005: record operation in header **/
1030
1031 { float matar[12] ;
1032 THD_dmat33 drmat ;
1033 THD_dfvec3 dvec ;
1034 char anam[64] ;
1035
1036 /* must transform [rmat,tvec] from dataset to DICOM coords */
1037 /* (the inverse of the DICOM to dataset transforms earlier) */
1038
1039 pp = DBLE_mat_to_dicomm(dset) ; ppt = TRANSPOSE_DMAT(pp) ;
1040 drmat = DMAT_MUL(pp,rmat) ; drmat = DMAT_MUL(drmat,ppt);
1041 dvec = DMATVEC (pp,tvec) ;
1042
1043 UNLOAD_DMAT(drmat,matar[0],matar[1],matar[2],
1044 matar[4],matar[5],matar[6],
1045 matar[8],matar[9],matar[10] ) ;
1046 UNLOAD_DFVEC3(dvec,matar[3],matar[7],matar[11]) ;
1047 sprintf(anam,"ROTATE_MATVEC_%06d",ival) ;
1048 THD_set_float_atr( dset->dblk , anam , 12 , matar ) ;
1049 }
1050
1051 if( verb ) fprintf(stderr,".") ;
1052
1053 if( clipit ){ /* 11 Apr 2000 */
1054 register int ii ; register float bb,tt ;
1055 bb = cbot ; tt = ctop ;
1056 for( ii=0 ; ii < nvox ; ii++ ){
1057 if( fvol[ii] < bb ) fvol[ii] = bb ;
1058 else if( fvol[ii] > tt ) fvol[ii] = tt ;
1059 }
1060 }
1061
1062 EDIT_coerce_type( nvox , MRI_float,fvol ,
1063 DSET_BRICK_TYPE(dset,ival),DSET_ARRAY(dset,ival) ) ;
1064
1065 } /* end of loop over sub-bricks */
1066
1067 if( verb ){
1068 cputim = COX_cpu_time() - cputim ;
1069 fprintf(stderr,"\n+++ CPU time=%10.3g s" , cputim) ;
1070 if( nval > 1 ) fprintf(stderr," [= %10.3g s/sub-brick]" , cputim/nval) ;
1071 fprintf(stderr,"\n+++ Writing dataset to disk in %s", DSET_BRIKNAME(dset) ) ;
1072 }
1073
1074 dset->dblk->master_nvals = 0 ; /* 11 Apr 2000 hack */
1075 DSET_write(dset) ;
1076 if( verb ) fprintf(stderr,"\n") ;
1077 exit(0) ;
1078 }
1079
1080 /*=================================================================================*/
1081
1082 #include <ctype.h>
1083
1084 #define NBUF 1024 /* line buffer size */
1085
1086 #define DUPOUT(n) fprintf(fpout,"%s",linbuf+n)
1087
rotate_stdin_points(THD_dfvec3 xyzorg,THD_dmat33 rmat,int dcode,THD_dfvec3 tvec)1088 static void rotate_stdin_points( THD_dfvec3 xyzorg, THD_dmat33 rmat,
1089 int dcode, THD_dfvec3 tvec )
1090 {
1091 char linbuf[NBUF] , *cp ;
1092 FILE *fpin=stdin , *fpout=stdout ;
1093 int ii , kk , nbuf , ll , nn , ld ;
1094 double xx,yy,zz ;
1095 THD_dfvec3 xyz ;
1096
1097 if( verb ){
1098 DUMP_DMAT33("Rotation",rmat) ;
1099 }
1100
1101 /** rmat = DMAT_INV(rmat) ; **/
1102
1103 /*-- loop over input lines --*/
1104
1105 ll = ld = 0 ;
1106 while(1){
1107 ll++ ; /* line count */
1108 cp = afni_fgets( linbuf , NBUF , fpin ) ; /* read the line */
1109 if( cp == NULL ) break ; /* end of file => end of loop */
1110 kk = strlen(linbuf) ;
1111 if( kk == 0 ) continue ; /* empty line => get next line */
1112
1113 /* find 1st nonblank */
1114
1115 for( ii=0 ; ii < kk && isspace(linbuf[ii]) ; ii++ ) ; /* nada */
1116 if( ii == kk || /* all blanks */
1117 (linbuf[ii] == '/' && linbuf[ii+1] == '/') ){ /* or comment */
1118
1119 DUPOUT(0) ; continue ;
1120 }
1121
1122 /* scan line for data */
1123
1124 nn = sscanf(linbuf+ii , "%lf%lf%lf%n" , &xx,&yy,&zz,&nbuf ) ;
1125 if( nn < 3 ){
1126 fprintf(stderr,"+++ WARNING: input line %d was incomplete\n",ll) ;
1127 continue ;
1128 }
1129 nbuf += ii ; /* position of next character after zz */
1130
1131 /* process vector */
1132
1133 LOAD_DFVEC3(xyz , xx,yy,zz) ;
1134 xyz = SUB_DFVEC3(xyz,xyzorg) ;
1135 if( dcode == DELTA_BEFORE ) xyz = ADD_DFVEC3(xyz,tvec) ;
1136 xyz = DMATVEC(rmat,xyz) ;
1137 if( dcode == DELTA_AFTER ) xyz = ADD_DFVEC3(xyz,tvec) ;
1138 xyz = ADD_DFVEC3(xyz,xyzorg) ;
1139
1140 fprintf(fpout,"%g %g %g%s",xyz.xyz[0],xyz.xyz[1],xyz.xyz[2],linbuf+nbuf) ;
1141 ld++ ;
1142
1143 } /* end of loop over input lines */
1144
1145 if( verb )
1146 fprintf(stderr,"-points: read %d lines, wrote %d lines\n",ll-1,ld) ;
1147 }
1148
1149 /*-----------------------------------------------------------------------
1150 19 Jun 2001:
1151 Get a 6 x N image of the rotation parameters
1152 -------------------------------------------------------------------------*/
1153
get_dfile_params(char * fname,int mode)1154 MRI_IMAGE * get_dfile_params( char * fname , int mode )
1155 {
1156 MRI_IMAGE *outim , *flim ;
1157 float *oar , *far ;
1158 int nx,ny , ii,jj ;
1159
1160 if( fname == NULL || fname[0] == '\0' ) return NULL ;
1161
1162 flim = mri_read_ascii( fname ) ;
1163 if( flim == NULL ) return NULL ;
1164
1165 nx = flim->nx ; ny = flim->ny ; far = MRI_FLOAT_PTR(flim) ;
1166
1167 if( mode == MODE_DFILE ){ /* skip 1st element of each row */
1168 if( nx < 7 ){
1169 fprintf(stderr,"** -dfile %s has too few columns!\n",fname) ;
1170 mri_free(flim) ; return NULL ;
1171 }
1172 outim = mri_new( 6 , ny , MRI_float ) ;
1173 oar = MRI_FLOAT_PTR(outim) ;
1174
1175 for( jj=0 ; jj < ny ; jj++ ){
1176 oar[0+6*jj] = far[1+nx*jj] ;
1177 oar[1+6*jj] = far[2+nx*jj] ;
1178 oar[2+6*jj] = far[3+nx*jj] ;
1179 oar[3+6*jj] = far[4+nx*jj] ;
1180 oar[4+6*jj] = far[5+nx*jj] ;
1181 oar[5+6*jj] = far[6+nx*jj] ;
1182 }
1183
1184 mri_free(flim) ;
1185
1186 } else if( mode == MODE_1DFILE ){ /* first 6 elements of each row */
1187
1188 if( nx < 6 ){
1189 fprintf(stderr,"** -1Dfile %s has too few columns!\n",fname) ;
1190 mri_free(flim) ; return NULL ;
1191 } else if( nx == 6 ){
1192 outim = flim ;
1193 } else {
1194 outim = mri_new( 6 , ny , MRI_float ) ;
1195 oar = MRI_FLOAT_PTR(outim) ;
1196 for( jj=0 ; jj < ny ; jj++ ){
1197 oar[0+6*jj] = far[0+nx*jj] ;
1198 oar[1+6*jj] = far[1+nx*jj] ;
1199 oar[2+6*jj] = far[2+nx*jj] ;
1200 oar[3+6*jj] = far[3+nx*jj] ;
1201 oar[4+6*jj] = far[4+nx*jj] ;
1202 oar[5+6*jj] = far[5+nx*jj] ;
1203 }
1204 mri_free(flim) ;
1205 }
1206
1207 } else {
1208 fprintf(stderr,"** get_dfile_params: illegal mode=%d\n",mode) ;
1209 mri_free(flim) ; return NULL ;
1210 }
1211
1212 return outim ;
1213 }
1214