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 "vecmat.h"
8 #include "mrilib.h"
9
10 /***************************************************************************/
11
12 #define ERREX(str) (fprintf(stderr,"** %s\n",str),exit(1))
13
14 THD_dmat33 DBLE_mat_to_dicomm ( THD_3dim_dataset * ) ; /* at end of file */
15 static void TAG_read( char * tfile ) ;
16
17 #define ROTATION 1 /* matrix_type options */
18 #define AFFINE 2
19 #define ROTSCL 3
20
21 static THD_usertaglist * mytagset = NULL ; /* to hold a tagset from -tagset */
22
usage_3dTagalign(int detail)23 void usage_3dTagalign(int detail) {
24 printf("Usage: 3dTagalign [options] dset\n"
25 "Rotates/translates dataset 'dset' to be aligned with the master,\n"
26 "using the tagsets embedded in their .HEAD files.\n"
27 "\n"
28 "Options:\n"
29 " -master mset = Use dataset 'mset' as the master dataset\n"
30 " [this is a nonoptional option]\n"
31 #if 0
32 "\n"
33 " -wtval = Use the numerical value attached to each tag\n"
34 " in the master as weighting factor for that tag\n"
35 " -wt1D ff = Read the *.1D file to use as a vector of weights\n"
36 " for each tag\n"
37 " N.B.: The default is to weight each tag the same.\n"
38 " If you use weights, a larger weight means to\n"
39 " count aligning that tag as more important.\n"
40 #endif
41 "\n"
42 " -tagset tfile = Use the tagset in the .tag file instead of dset.\n"
43 "\n"
44 " -nokeeptags = Don't put transformed locations of dset's tags\n"
45 " into the output dataset [default = keep tags]\n"
46 "\n"
47 " -matvec mfile = Write the matrix+vector of the transformation to\n"
48 " file 'mfile'. This can be used as input to the\n"
49 " '-matvec_in2out' option of 3dWarp, if you want\n"
50 " to align other datasets in the same way (e.g.,\n"
51 " functional datasets).\n"
52 #if 0
53 " N.B.: The matrix+vector of the transformation is also\n"
54 " saved in the .HEAD file of the output dataset\n"
55 " (unless -dummy is used). You can use the\n"
56 " output dataset from 3dTagalign as the dataset\n"
57 " for 3drotates's '-matvec_dset' option; this\n"
58 " lets you avoid using an intermediate mfile.\n"
59 "\n"
60 " -3dWarp = Use the '3dWarp' function to transform the dataset,\n"
61 " rather than the '3drotate' function. The output\n"
62 " dataset will be computed on the same grid as the\n"
63 " master dataset.\n"
64 " N.B.: This option is implied by '-affine' and '-rotscl'.\n"
65 #endif
66 "\n"
67 " -rotate = Compute the best transformation as a rotation + shift.\n"
68 " This is the default.\n"
69 "\n"
70 " -affine = Compute the best transformation as a general affine\n"
71 " map rather than just a rotation + shift. In all\n"
72 " cases, the transformation from input to output\n"
73 " coordinates is of the form\n"
74 " [out] = [R] [in] + [V]\n"
75 " where [R] is a 3x3 matrix and [V] is a 3-vector.\n"
76 " By default, [R] is computed as a proper (det=1)\n"
77 " rotation matrix (3 parameters). The '-affine'\n"
78 " option says to fit [R] as a general matrix\n"
79 " (9 parameters).\n"
80 " N.B.: An affine transformation can rotate, rescale, and\n"
81 " shear the volume. Be sure to look at the dataset\n"
82 " before and after to make sure things are OK.\n"
83 "\n"
84 " -rotscl = Compute transformation as a rotation times an isotropic\n"
85 " scaling; that is, [R] is an orthogonal matrix times\n"
86 " a scalar.\n"
87 " N.B.: '-affine' and '-rotscl' do unweighted least squares.\n"
88 "\n"
89 " -prefix pp = Use 'pp' as the prefix for the output dataset.\n"
90 " [default = 'tagalign']\n"
91 " -verb = Print progress reports\n"
92 " -dummy = Don't actually rotate the dataset, just compute\n"
93 " the transformation matrix and vector. If\n"
94 " '-matvec' is used, the mfile will be written.\n"
95 " -linear }\n"
96 " -cubic } = Chooses spatial interpolation method.\n"
97 " -NN } = [default = cubic]\n"
98 " -quintic }\n"
99 "\n"
100 "Nota Bene:\n"
101 "* The transformation is carried out\n"
102 " using the same methods as program 3dWarp.\n"
103 "\n"
104 "Author: RWCox - 16 Jul 2000, etc.\n"
105 ) ;
106 PRINT_COMPILE_DATE ;
107 return;
108 }
109 /*-------------------------------------------------------------------------*/
110
main(int argc,char * argv[])111 int main( int argc , char * argv[] )
112 {
113 THD_dfvec3 *xx , *yy , dv ;
114 int nvec=0 , ii,jj, iarg ;
115 THD_dvecmat rt , rtinv ;
116 THD_dmat33 pp,ppt , rr ;
117 THD_dfvec3 tt ;
118
119 THD_3dim_dataset *mset=NULL , *dset=NULL ;
120 double *ww=NULL ;
121 int nww=0 ;
122 int keeptags=1 , wtval=0 , verb=0 , dummy=0 ;
123 char * prefix = "tagalign" , *mfile=NULL , *tfile=NULL;
124
125 float *fvol , cbot,ctop , dsum ;
126 int nval , nvox , clipit , ival, RMETH=MRI_CUBIC;
127
128 float matar[12] ;
129
130 int use_3dWarp=1 , matrix_type=ROTATION ;
131
132 mainENTRY("3dTagalign main");
133
134 /*--- help? ---*/
135
136 /*- scan args -*/
137 iarg = 1 ; RMETH=MRI_CUBIC;
138 while( iarg < argc && argv[iarg][0] == '-' ){
139
140 /*-----*/
141
142 if( strcmp(argv[iarg],"-h") == 0 ||
143 strcmp(argv[iarg],"-help") == 0){ /* 22 Apr 2003 */
144 usage_3dTagalign(strlen(argv[iarg]) > 3 ? 2:1);
145 exit(0);
146 }
147
148 /*-----*/
149
150 if( strcmp(argv[iarg],"-NN") == 0 ){
151 RMETH = MRI_NN ; iarg++ ; continue ;
152 }
153 if( strcmp(argv[iarg],"-linear") == 0 ){
154 RMETH = MRI_LINEAR ; iarg++ ; continue ;
155 }
156 if( strcmp(argv[iarg],"-cubic") == 0 ){
157 RMETH = MRI_CUBIC ; iarg++ ; continue ;
158 }
159 if( strcmp(argv[iarg],"-quintic") == 0 ){
160 RMETH = MRI_QUINTIC ; iarg++ ; continue ;
161 }
162
163 /*-----*/
164
165 if( strcmp(argv[iarg],"-rotate") == 0 ){ /* 22 Apr 2003 */
166 matrix_type = ROTATION ; use_3dWarp = 1 ;
167 iarg++ ; continue ;
168 }
169
170 /*-----*/
171
172 if( strcmp(argv[iarg],"-affine") == 0 ){ /* 21 Apr 2003 */
173 matrix_type = AFFINE ; use_3dWarp = 1 ;
174 iarg++ ; continue ;
175 }
176
177 /*-----*/
178
179 if( strcmp(argv[iarg],"-rotscl") == 0 ){ /* 22 Apr 2003 */
180 matrix_type = ROTSCL ; use_3dWarp = 1 ;
181 iarg++ ; continue ;
182 }
183
184 #if 0
185 /*-----*/
186
187 if( strcmp(argv[iarg],"-3dWarp") == 0 ){ /* 21 Apr 2003 */
188 use_3dWarp = 1 ;
189 iarg++ ; continue ;
190 }
191 #endif
192
193 /*-----*/
194
195 if( strcmp(argv[iarg],"-master") == 0 ){
196 if( mset != NULL ) ERREX("Can only have one -master option") ;
197 if( ++iarg >= argc ) ERREX("Need an argument after -master") ;
198
199 mset = THD_open_dataset( argv[iarg] ) ;
200
201 if( mset == NULL ) ERREX("Can't open -master dataset") ;
202 if( mset->tagset == NULL ) ERREX("No tags in -master dataset") ;
203 if( TAGLIST_COUNT(mset->tagset) < 3 ) ERREX("Not enough tags in -master dataset") ;
204
205 for( nvec=ii=0 ; ii < TAGLIST_COUNT(mset->tagset) ; ii++ )
206 if( TAG_SET(TAGLIST_SUBTAG(mset->tagset,ii)) ) nvec++ ;
207
208 if( nvec < 3 ) ERREX("Not enough tags set in -master dataset") ;
209
210 if( nvec < TAGLIST_COUNT(mset->tagset) )
211 fprintf(stderr,"++ WARNING: not all tags are set in -master dataset\n") ;
212
213 if( verb ) fprintf(stderr,"++ Found %d tags in -master dataset\n",nvec) ;
214
215 iarg++ ; continue ;
216 }
217
218 #if 0
219 /*-----*/
220
221 if( strcmp(argv[iarg],"-wtval") == 0 ){
222 if( ww != NULL ) ERREX("Can't have -wtval after -wt1D") ;
223 wtval++ ;
224 iarg++ ; continue ;
225 }
226
227 /*-----*/
228
229 if( strcmp(argv[iarg],"-wt1D") == 0 ){
230 MRI_IMAGE * wtim ; float * wtar ;
231
232 if( wtval ) ERREX("Can't have -wt1D after -wtval") ;
233 if( ww != NULL ) ERREX("Can't have two -wt1D options!") ;
234 if( ++iarg >= argc ) ERREX("Need an argument after -wt1D") ;
235
236 wtim = mri_read_1D( argv[iarg] ) ;
237
238 if( wtim == NULL ) ERREX("Can't read -wtim file") ;
239 if( wtim->ny > 1 ) ERREX("-wtim file has more than one columm") ;
240
241 wtar = MRI_FLOAT_PTR(wtim) ;
242 ww = (double *) malloc(sizeof(double)*wtim->nx) ; nww = wtim->nx ;
243 for( ii=0 ; ii < nww ; ii++ ){
244 ww[ii] = (double) wtar[ii] ;
245 if( ww[ii] < 0.0 ) ERREX("Negative value found in -wt1D file") ;
246 }
247
248 mri_free(wtim) ;
249 iarg++ ; continue ;
250 }
251 #endif
252
253 /*-----*/
254
255 if( strcmp(argv[iarg],"-nokeeptags") == 0 ){
256 keeptags = 0 ;
257 iarg++ ; continue ;
258 }
259
260 /*-----*/
261
262 if( strncmp(argv[iarg],"-verb",5) == 0 ){
263 verb++ ;
264 iarg++ ; continue ;
265 }
266
267 /*-----*/
268
269 if( strcmp(argv[iarg],"-dummy") == 0 ){
270 dummy++ ;
271 iarg++ ; continue ;
272 }
273
274 /*-----*/
275
276 if( strcmp(argv[iarg],"-prefix") == 0 ){
277 if( ++iarg >= argc ) ERREX("Need an argument after -prefix") ;
278 prefix = argv[iarg] ;
279 if( !THD_filename_ok(prefix) ) ERREX("-prefix string is illegal") ;
280 iarg++ ; continue ;
281 }
282
283 /*-----*/
284
285 if( strcmp(argv[iarg],"-matvec") == 0 ){
286 if( ++iarg >= argc ) ERREX("Need an argument after -matvec") ;
287 mfile = argv[iarg] ;
288 if( !THD_filename_ok(mfile) ) ERREX("-matvec string is illegal") ;
289 iarg++ ; continue ;
290 }
291
292 /*-----*/
293
294 if( strcmp(argv[iarg],"-tagset") == 0 ){
295 if( ++iarg >= argc ) ERREX("Need an argument after -tagset") ;
296 tfile = argv[iarg] ;
297 if( !THD_filename_ok(tfile) ) ERREX("-tagset string is illegal") ;
298
299 if( mytagset != NULL ) ERREX("Can only have one -tagset option") ;
300 mytagset = (THD_usertaglist *)calloc(1, sizeof(THD_usertaglist)) ;
301 if( !mytagset ) ERREX("insufficient RAM for taglist") ;
302 TAG_read(tfile) ;
303
304 iarg++ ; continue ;
305 }
306
307
308 /*-----*/
309
310 fprintf(stderr,"** Unknown option: %s\n",argv[iarg]) ;
311 suggest_best_prog_option(argv[0], argv[iarg]);
312 exit(1) ;
313
314 } /* end of scanning command line for options */
315
316 if( argc < 2 ){
317 ERROR_message("Too few options");
318 usage_3dTagalign(0);
319 exit(1) ;
320 }
321
322
323 if( mset == NULL ) ERREX("No -master option found on command line") ;
324
325 #if 0
326 if( ww != NULL && nww < nvec ) ERREX("Not enough weights found in -wt1D file") ;
327
328 /*-- if -wtval, setup weights from master tag values --*/
329
330 if( wtval ){
331 ww = (double *) malloc(sizeof(double)*nvec) ; nww = nvec ;
332 for( ii=jj=0 ; ii < TAGLIST_COUNT(mset->tagset) ; ii++ ){
333 if( TAG_SET(TAGLIST_SUBTAG(mset->tagset,ii)) ){
334 ww[jj] = (double) TAG_VAL(TAGLIST_SUBTAG(mset->tagset,ii)) ;
335
336 if( ww[jj] < 0.0 ) ERREX("Negative value found in -master tag values") ;
337 jj++ ;
338 }
339 }
340 }
341 #endif
342
343 /*-- read input dataset (to match to master dataset) --*/
344
345 if( iarg >= argc ) ERREX("No input dataset?") ;
346
347 dset = THD_open_dataset( argv[iarg] ) ;
348
349 if( dset == NULL ) ERREX("Can't open input dataset") ;
350
351 /* if user specified one, unceremoniously overwrite the dset's tagset */
352 if( mytagset != NULL )
353 dset->tagset = mytagset ;
354
355 if( dset->tagset == NULL ) ERREX("No tags in input dataset") ;
356
357 if( TAGLIST_COUNT(dset->tagset) != TAGLIST_COUNT(mset->tagset) )
358 ERREX("Tag counts don't match in -master and input") ;
359
360 /* check if set tags match exactly */
361
362 for( ii=0 ; ii < TAGLIST_COUNT(mset->tagset) ; ii++ ){
363 if( TAG_SET(TAGLIST_SUBTAG(mset->tagset,ii)) !=
364 TAG_SET(TAGLIST_SUBTAG(dset->tagset,ii)) )
365 ERREX("Set tags don't match in -master and input") ;
366 }
367
368 /*-- load vector lists: xx=master, yy=input --*/
369
370 xx = (THD_dfvec3 *) malloc( sizeof(THD_dfvec3) * nvec ) ;
371 yy = (THD_dfvec3 *) malloc( sizeof(THD_dfvec3) * nvec ) ;
372 dsum = 0.0 ;
373 for( ii=jj=0 ; ii < nvec ; ii++ ){
374 if( TAG_SET(TAGLIST_SUBTAG(mset->tagset,ii)) ){
375
376 /* N.B.: these are in Dicom order now */
377 LOAD_DFVEC3( xx[jj] ,
378 TAG_X( TAGLIST_SUBTAG(mset->tagset,ii) ) ,
379 TAG_Y( TAGLIST_SUBTAG(mset->tagset,ii) ) ,
380 TAG_Z( TAGLIST_SUBTAG(mset->tagset,ii) ) ) ;
381
382 LOAD_DFVEC3( yy[jj] ,
383 TAG_X( TAGLIST_SUBTAG(dset->tagset,ii) ) ,
384 TAG_Y( TAGLIST_SUBTAG(dset->tagset,ii) ) ,
385 TAG_Z( TAGLIST_SUBTAG(dset->tagset,ii) ) ) ;
386
387 dv = SUB_DFVEC3( xx[jj] , yy[jj] ) ;
388 dsum += dv.xyz[0]*dv.xyz[0] + dv.xyz[1]*dv.xyz[1] + dv.xyz[2]*dv.xyz[2] ;
389
390 jj++ ;
391 }
392 }
393
394 dsum = sqrt(dsum/nvec) ;
395 fprintf(stderr,"++ RMS distance between tags before = %.2f mm\n" , dsum ) ;
396
397 /*-- compute best transformation from mset to dset coords --*/
398
399 switch( matrix_type ){
400 default:
401 case ROTATION:
402 rt = DLSQ_rot_trans( nvec , yy , xx , ww ) ; /* in thd_rot3d.c */
403 break ;
404
405 case AFFINE:
406 rt = DLSQ_affine ( nvec , yy , xx ) ; /* 21 Apr 2003 */
407 break ;
408
409 case ROTSCL:
410 rt = DLSQ_rotscl ( nvec , yy , xx , (DSET_NZ(dset)==1) ? 2 : 3 ) ;
411 break ;
412 }
413 rtinv = INV_DVECMAT(rt) ;
414
415 /*-- check for floating point legality --*/
416
417 nval = 0 ;
418 for( ii=0 ; ii < 3 ; ii++ ){
419 dsum = rt.vv.xyz[ii] ; nval += thd_floatscan(1,&dsum) ;
420 for( jj=0 ; jj < 3 ; jj++ ){
421 dsum = rt.mm.mat[ii][jj] ; nval += thd_floatscan(1,&dsum) ;
422 }
423 }
424 if( nval > 0 ){
425 fprintf(stderr,"** Floating point errors during calculation\n"
426 "** of transform matrix and translation vector\n" ) ;
427 exit(1) ;
428 }
429
430 /*-- check for rotation matrix legality --*/
431
432 dsum = DMAT_DET(rt.mm) ;
433
434 if( dsum == 0.0 || (matrix_type == ROTATION && fabs(dsum-1.0) > 0.01) ){
435 fprintf(stderr,"** Invalid transform matrix computed: tags dependent?\n"
436 "** computed [matrix] and [vector] follow:\n" ) ;
437
438 for( ii=0 ; ii < 3 ; ii++ )
439 fprintf(stderr," [ %10.5f %10.5f %10.5f ] [ %10.5f ] \n",
440 rt.mm.mat[ii][0],rt.mm.mat[ii][1],rt.mm.mat[ii][2],rt.vv.xyz[ii] );
441
442 exit(1) ;
443 }
444
445 /*-- print summary --*/
446
447 if( verb ){
448 fprintf(stderr,"++ Matrix & Vector [Dicom: x=R-L; y=A-P; z=I-S]\n") ;
449 for( ii=0 ; ii < 3 ; ii++ )
450 fprintf(stderr," %10.5f %10.5f %10.5f %10.5f\n",
451 rt.mm.mat[ii][0],rt.mm.mat[ii][1],rt.mm.mat[ii][2],rt.vv.xyz[ii] );
452 }
453
454 if( matrix_type == ROTATION || matrix_type == ROTSCL ){
455 double theta, costheta , dist , fac=1.0 ;
456
457 if( matrix_type == ROTSCL ){
458 fac = DMAT_DET(rt.mm); fac = fabs(fac);
459 if( DSET_NZ(dset) == 1 ) fac = sqrt(fac) ;
460 else fac = cbrt(fac) ;
461 }
462
463 costheta = 0.5 * sqrt(1.0 + DMAT_TRACE(rt.mm)/fac ) ;
464 theta = 2.0 * acos(costheta) * 180/3.14159265 ;
465 dist = SIZE_DFVEC3(rt.vv) ;
466
467 fprintf(stderr,"++ Total rotation=%.2f degrees; translation=%.2f mm; scaling=%.2f\n",
468 theta,dist,fac) ;
469 }
470
471 if( mfile ){
472 FILE * mp ;
473
474 if( THD_is_file(mfile) )
475 fprintf(stderr,"++ Warning: -matvec will overwrite file %s\n",mfile) ;
476
477 mp = fopen(mfile,"w") ;
478 if( mp == NULL ){
479 fprintf(stderr,"** Can't write to -matvec %s\n",mfile) ;
480 } else {
481 for( ii=0 ; ii < 3 ; ii++ )
482 fprintf(mp," %10.5f %10.5f %10.5f %10.5f\n",
483 rt.mm.mat[ii][0],rt.mm.mat[ii][1],rt.mm.mat[ii][2],rt.vv.xyz[ii] );
484 fclose(mp) ;
485 if( verb ) fprintf(stderr,"++ Wrote matrix+vector to %s\n",mfile) ;
486 }
487 }
488
489 if( dummy ){
490 fprintf(stderr,"++ This was a -dummy run: no output dataset\n") ; exit(0) ;
491 }
492
493 /*-- 21 Apr 2003: transformation can be done the old way (a la 3drotate),
494 or the new way (a la 3dWarp). --*/
495
496 #if 0
497 if( !use_3dWarp ){ /**** the old way ****/
498
499 /*-- now must scramble the rotation matrix and translation
500 vector from Dicom coordinate order to dataset brick order --*/
501
502 pp = DBLE_mat_to_dicomm( dset ) ;
503 ppt = TRANSPOSE_DMAT(pp) ;
504 rr = DMAT_MUL(ppt,rt.mm) ; rr = DMAT_MUL(rr,pp) ; tt = DMATVEC(ppt,rt.vv) ;
505
506 /*-- now create the output dataset by screwing with the input dataset
507 (this code is adapted from 3drotate.c) --*/
508
509 DSET_mallocize(dset) ;
510 DSET_load( dset ) ; CHECK_LOAD_ERROR(dset) ;
511 dset->idcode = MCW_new_idcode() ;
512 dset->dblk->diskptr->storage_mode = STORAGE_BY_BRICK ; /* 14 Jan 2004 */
513 EDIT_dset_items( dset ,
514 ADN_prefix , prefix ,
515 ADN_label1 , prefix ,
516 ADN_none ) ;
517
518 if( !THD_ok_overwrite() &&
519 (THD_deathcon() && THD_is_file(dset->dblk->diskptr->header_name) )){
520 fprintf(stderr,
521 "** Output file %s already exists -- cannot continue!\n",
522 dset->dblk->diskptr->header_name ) ;
523 exit(1) ;
524 }
525
526 tross_Make_History( "3dTagalign" , argc,argv , dset ) ;
527
528 /*-- if desired, keep old tagset --*/
529
530 if( keeptags ){
531 THD_dfvec3 rv ;
532
533 dsum = 0.0 ;
534 for( jj=ii=0 ; ii < TAGLIST_COUNT(dset->tagset) ; ii++ ){
535 if( TAG_SET(TAGLIST_SUBTAG(dset->tagset,ii)) ){
536 rv = DMATVEC( rt.mm , yy[jj] ) ; /* operating on */
537 rv = ADD_DFVEC3( rt.vv , rv ) ; /* Dicom order */
538
539 dv = SUB_DFVEC3( xx[jj] , rv ) ;
540 dsum += dv.xyz[0]*dv.xyz[0] + dv.xyz[1]*dv.xyz[1]
541 + dv.xyz[2]*dv.xyz[2] ;
542
543 UNLOAD_DFVEC3( rv , TAG_X( TAGLIST_SUBTAG(dset->tagset,ii) ) ,
544 TAG_Y( TAGLIST_SUBTAG(dset->tagset,ii) ) ,
545 TAG_Z( TAGLIST_SUBTAG(dset->tagset,ii) ) ) ;
546
547 jj++ ;
548 }
549 }
550 dsum = sqrt(dsum/nvec) ;
551 fprintf(stderr,"++ RMS distance between tags after = %.2f mm\n" , dsum ) ;
552
553 } else {
554 myRwcFree(dset->tagset) ; /* send it to the dustbin */
555 }
556
557 /*-- rotate sub-bricks --*/
558
559 if( verb ) fprintf(stderr,"++ computing output BRIK") ;
560
561 nvox = DSET_NVOX(dset) ;
562 nval = DSET_NVALS(dset) ;
563 fvol = (float *) malloc( sizeof(float) * nvox ) ;
564
565 THD_rota_method( MRI_HEPTIC ) ;
566 clipit = 1 ;
567
568 for( ival=0 ; ival < nval ; ival++ ){
569
570 /*- get sub-brick out of dataset -*/
571
572 EDIT_coerce_type( nvox ,
573 DSET_BRICK_TYPE(dset,ival),DSET_ARRAY(dset,ival) ,
574 MRI_float,fvol ) ;
575
576 if( clipit ){
577 register int ii ; register float bb,tt ;
578 bb = tt = fvol[0] ;
579 for( ii=1 ; ii < nvox ; ii++ ){
580 if( fvol[ii] < bb ) bb = fvol[ii] ;
581 else if( fvol[ii] > tt ) tt = fvol[ii] ;
582 }
583 cbot = bb ; ctop = tt ;
584 }
585
586 if( verb && nval < 5 ) fprintf(stderr,".") ;
587
588 /*- rotate it -*/
589
590 THD_rota_vol_matvec( DSET_NX(dset) , DSET_NY(dset) , DSET_NZ(dset) ,
591 fabs(DSET_DX(dset)) , fabs(DSET_DY(dset)) ,
592 fabs(DSET_DZ(dset)) ,
593 fvol , rr , tt ) ;
594
595 if( verb ) fprintf(stderr,".") ;
596
597 if( clipit ){
598 register int ii ; register float bb,tt ;
599 bb = cbot ; tt = ctop ;
600 for( ii=0 ; ii < nvox ; ii++ ){
601 if( fvol[ii] < bb ) fvol[ii] = bb ;
602 else if( fvol[ii] > tt ) fvol[ii] = tt ;
603 }
604 }
605
606 if( verb && nval < 5 ) fprintf(stderr,".") ;
607
608 /*- put it back into dataset -*/
609
610 EDIT_coerce_type( nvox, MRI_float,fvol ,
611 DSET_BRICK_TYPE(dset,ival),DSET_ARRAY(dset,ival) );
612
613 } /* end of loop over sub-brick index */
614
615 if( verb ) fprintf(stderr,":") ;
616
617 /* save matrix+vector into dataset, too */
618
619 UNLOAD_DMAT(rt.mm,matar[0],matar[1],matar[2],
620 matar[4],matar[5],matar[6],
621 matar[8],matar[9],matar[10] ) ;
622 UNLOAD_DFVEC3(rt.vv,matar[3],matar[7],matar[11]) ;
623 THD_set_atr( dset->dblk, "TAGALIGN_MATVEC", ATR_FLOAT_TYPE, 12, matar ) ;
624
625 /* write dataset to disk */
626
627 dset->dblk->master_nvals = 0 ; /* in case this was a mastered dataset */
628 DSET_write(dset) ;
629
630 if( verb ) fprintf(stderr,"\n") ;
631
632 } else
633 #endif
634 { /**** the new way: use 3dWarp type transformation ****/
635
636 THD_3dim_dataset *oset ;
637 THD_vecmat tran ;
638
639 #if 0
640 DFVEC3_TO_FVEC3( rt.vv , tran.vv ) ;
641 DMAT_TO_MAT ( rt.mm , tran.mm ) ;
642 #else
643 DFVEC3_TO_FVEC3( rtinv.vv , tran.vv ) ;
644 DMAT_TO_MAT ( rtinv.mm , tran.mm ) ;
645 #endif
646
647 mri_warp3D_method( RMETH ) ;
648 oset = THD_warp3D_affine( dset, tran, mset, prefix, 0, WARP3D_NEWDSET ) ;
649 if( oset == NULL ){
650 fprintf(stderr,"** ERROR: THD_warp3D() fails!\n"); exit(1);
651 }
652
653 tross_Copy_History( dset , oset ) ;
654 tross_Make_History( "3dTagalign" , argc,argv , oset ) ;
655
656 UNLOAD_DMAT(rt.mm,matar[0],matar[1],matar[2],
657 matar[4],matar[5],matar[6],
658 matar[8],matar[9],matar[10] ) ;
659 UNLOAD_DFVEC3(rt.vv,matar[3],matar[7],matar[11]) ;
660 THD_set_atr( oset->dblk, "TAGALIGN_MATVEC", ATR_FLOAT_TYPE, 12, matar ) ;
661
662 /*-- if desired, keep old tagset --*/
663
664 if( keeptags ){
665 THD_dfvec3 rv ;
666
667 oset->tagset = (THD_usertaglist *)calloc(1, sizeof(THD_usertaglist)) ;
668 if( !oset->tagset ) ERREX("insufficient RAM to copy taglist") ;
669 *(oset->tagset) = *(dset->tagset) ; /* swell foop */
670
671 dsum = 0.0 ;
672 for( jj=ii=0 ; ii < TAGLIST_COUNT(oset->tagset) ; ii++ ){
673 if( TAG_SET(TAGLIST_SUBTAG(oset->tagset,ii)) ){
674 rv = DMATVEC( rt.mm , yy[jj] ) ;
675 rv = ADD_DFVEC3( rt.vv , rv ) ;
676
677 dv = SUB_DFVEC3( xx[jj] , rv ) ;
678 dsum += dv.xyz[0]*dv.xyz[0] + dv.xyz[1]*dv.xyz[1]
679 + dv.xyz[2]*dv.xyz[2] ;
680
681 UNLOAD_DFVEC3( rv , TAG_X( TAGLIST_SUBTAG(oset->tagset,ii) ) ,
682 TAG_Y( TAGLIST_SUBTAG(oset->tagset,ii) ) ,
683 TAG_Z( TAGLIST_SUBTAG(oset->tagset,ii) ) ) ;
684
685 jj++ ;
686 }
687 }
688 dsum = sqrt(dsum/nvec) ;
689 fprintf(stderr,"++ RMS distance between tags after = %.2f mm\n" , dsum ) ;
690 }
691
692 DSET_write(oset) ;
693
694 } /* end of 3dWarp-like work */
695
696 exit(0) ;
697 }
698
699 /***********************************************************************/
700
701 #include "thd.h"
702
703 /*---------------------------------------------------------------------
704 This produces a permutation-like matrix that transforms from
705 brick axis coordinates to Dicom order coordinates.
706 -----------------------------------------------------------------------*/
707
DBLE_mat_to_dicomm(THD_3dim_dataset * dset)708 THD_dmat33 DBLE_mat_to_dicomm( THD_3dim_dataset * dset )
709 {
710 THD_dmat33 tod ;
711
712 LOAD_ZERO_DMAT(tod) ;
713
714 switch( dset->daxes->xxorient ){
715 case ORI_R2L_TYPE: tod.mat[0][0] = 1.0 ; break ;
716 case ORI_L2R_TYPE: tod.mat[0][0] = -1.0 ; break ;
717 case ORI_P2A_TYPE: tod.mat[1][0] = -1.0 ; break ;
718 case ORI_A2P_TYPE: tod.mat[1][0] = 1.0 ; break ;
719 case ORI_I2S_TYPE: tod.mat[2][0] = 1.0 ; break ;
720 case ORI_S2I_TYPE: tod.mat[2][0] = -1.0 ; break ;
721
722 default: THD_FATAL_ERROR("illegal xxorient code") ;
723 }
724
725 switch( dset->daxes->yyorient ){
726 case ORI_R2L_TYPE: tod.mat[0][1] = 1.0 ; break ;
727 case ORI_L2R_TYPE: tod.mat[0][1] = -1.0 ; break ;
728 case ORI_P2A_TYPE: tod.mat[1][1] = -1.0 ; break ;
729 case ORI_A2P_TYPE: tod.mat[1][1] = 1.0 ; break ;
730 case ORI_I2S_TYPE: tod.mat[2][1] = 1.0 ; break ;
731 case ORI_S2I_TYPE: tod.mat[2][1] = -1.0 ; break ;
732
733 default: THD_FATAL_ERROR("illegal yyorient code") ;
734 }
735
736 switch( dset->daxes->zzorient ){
737 case ORI_R2L_TYPE: tod.mat[0][2] = 1.0 ; break ;
738 case ORI_L2R_TYPE: tod.mat[0][2] = -1.0 ; break ;
739 case ORI_P2A_TYPE: tod.mat[1][2] = -1.0 ; break ;
740 case ORI_A2P_TYPE: tod.mat[1][2] = 1.0 ; break ;
741 case ORI_I2S_TYPE: tod.mat[2][2] = 1.0 ; break ;
742 case ORI_S2I_TYPE: tod.mat[2][2] = -1.0 ; break ;
743
744 default: THD_FATAL_ERROR("illegal zxorient code") ;
745 }
746
747 return tod ;
748 }
749
750 /*--------------------------------------------------------------------------*/
751
TAG_read(char * str)752 static void TAG_read( char * str ) /* copied from plug_tag.c */
753 {
754 char * cpt , buf[256] , quote ;
755 int ii , jj , ntog , lbuf , kk ;
756 FILE * fp ;
757 char new_label[MAX_TAG_LABEL] ;
758 float new_x , new_y , new_z , new_val ;
759 int new_ti , new_set ;
760
761 /*-- check for suffix on filename --*/
762
763 cpt = strstr( str , ".tag" ) ;
764 if( cpt == NULL ){
765 ii = strlen(str) ;
766 cpt = (char *)malloc(ii+8) ;
767 if( !cpt ) ERREX("insufficient RAM to append .tag to tag file name") ;
768 strcpy(cpt,str) ;
769 if( cpt[ii-1] != '.' ) strcat( cpt , "." ) ;
770 strcat( cpt , "tag" ) ;
771 str = cpt ;
772 }
773
774 /*-- open file --*/
775
776 fp = fopen( str , "r" ) ;
777 if( fp == NULL ) ERREX( "Can't open input file!" ) ;
778
779 /*-- scan each line for a tag definition --*/
780
781 ntog = 0 ;
782 while( ntog < MAX_TAG_NUM ){
783 cpt = afni_fgets( buf , 256 , fp ); /* read line from disk */
784 if( cpt == NULL ) break ; /* nothing => exit */
785 if( buf[0] == '#' ) continue ; /* comment => skip this line */
786 if( buf[0] == '\n' ) continue ;
787 if( buf[0] == '\0' ) continue ;
788
789 lbuf = strlen(buf) ; /* skip whitespace at start */
790 jj = 0 ;
791 while( jj < lbuf && isspace(buf[jj]) ) jj++ ;
792 if( jj == lbuf ) continue ; /* was a blank line => skip */
793 if( buf[jj] == '#' ) continue ; /* comment */
794
795 /* scan for new label */
796
797 if( buf[jj] == '\'' || buf[jj] == '\"' ){ /* scan to matching quote */
798 quote = buf[jj] ; jj++ ; kk = jj ;
799 while( kk < lbuf && buf[kk] != quote ) kk++ ;
800 if( kk == lbuf ) kk-- ;
801 } else { /* scan to nonblank */
802 kk = jj+1 ;
803 while( kk < lbuf && !isspace(buf[kk]) ) kk++ ;
804 }
805 for( ii=0 ; ii < MAX_TAG_LABEL-1 && jj < kk ; ii++,jj++ )
806 new_label[ii] = buf[jj] ;
807 new_label[ii] = '\0' ;
808 if( strlen(new_label) == 0 ) continue ; /* error => skip */
809 jj = kk+1 ;
810
811 /* scan for x y z val ti */
812
813 new_set = new_ti = 0 ;
814 new_x = new_y = new_z = new_val = 0.0 ;
815 if( jj < lbuf-4 ){
816 kk = sscanf( buf+jj , "%f %f %f %f %d" ,
817 &new_x , &new_y , &new_z , &new_val , &new_ti ) ;
818 if( kk >= 3 ) new_set = 1 ; /* got x y z, at least */
819 }
820
821 /* set values */
822
823 strcpy( mytagset->tag[ntog].label , new_label ) ;
824 mytagset->tag[ntog].set = new_set ;
825 mytagset->tag[ntog].ti = new_ti ;
826 mytagset->tag[ntog].x = new_x ;
827 mytagset->tag[ntog].y = new_y ;
828 mytagset->tag[ntog].z = new_z ;
829 mytagset->tag[ntog].val = new_val ;
830 ntog++ ;
831 }
832
833 fclose(fp) ; /* done with file */
834
835 if( ntog == 0 ){ /* no tags ==> error */
836 sprintf(buf , "Couldn't read tagset from file %s" , str ) ;
837 ERREX(buf) ;
838 }
839
840 mytagset->num = ntog ;
841 }
842
843