00001
00002
00003
00004
00005
00006
00007 #include "vecmat.h"
00008 #include "mrilib.h"
00009
00010
00011
00012 #define ERREX(str) (fprintf(stderr,"** %s\n",str),exit(1))
00013
00014 THD_dmat33 DBLE_mat_to_dicomm ( THD_3dim_dataset * ) ;
00015
00016 #define ROTATION 1
00017 #define AFFINE 2
00018 #define ROTSCL 3
00019
00020
00021
00022 int main( int argc , char * argv[] )
00023 {
00024 THD_dfvec3 *xx , *yy , dv ;
00025 int nvec , ii,jj, iarg ;
00026 THD_dvecmat rt , rtinv ;
00027 THD_dmat33 pp,ppt , rr ;
00028 THD_dfvec3 tt ;
00029
00030 THD_3dim_dataset *mset=NULL , *dset=NULL ;
00031 double *ww=NULL ;
00032 int nww=0 ;
00033 int keeptags=1 , wtval=0 , verb=0 , dummy=0 ;
00034 char * prefix = "tagalign" , *mfile=NULL ;
00035
00036 float *fvol , cbot,ctop , dsum ;
00037 int nval , nvox , clipit , ival ;
00038
00039 float matar[12] ;
00040
00041 int use_3dWarp=1 , matrix_type=ROTATION ;
00042
00043
00044
00045 if( argc < 2 || strcmp(argv[1],"-help") == 0 ){
00046 printf("Usage: 3dTagalign [options] dset\n"
00047 "Rotates/translates dataset 'dset' to be aligned with the master,\n"
00048 "using the tagsets embedded in their .HEAD files.\n"
00049 "\n"
00050 "Options:\n"
00051 " -master mset = Use dataset 'mset' as the master dataset\n"
00052 " [this is a nonoptional option]\n"
00053 #if 0
00054 "\n"
00055 " -wtval = Use the numerical value attached to each tag\n"
00056 " in the master as weighting factor for that tag\n"
00057 " -wt1D ff = Read the *.1D file to use as a vector of weights\n"
00058 " for each tag\n"
00059 " N.B.: The default is to weight each tag the same.\n"
00060 " If you use weights, a larger weight means to\n"
00061 " count aligning that tag as more important.\n"
00062 #endif
00063 "\n"
00064 " -nokeeptags = Don't put transformed locations of dset's tags\n"
00065 " into the output dataset [default = keep tags]\n"
00066 "\n"
00067 " -matvec mfile = Write the matrix+vector of the transformation to\n"
00068 " file 'mfile'. This can be used as input to the\n"
00069 " '-matvec_out2in' option of 3dWarp, if you want\n"
00070 " to align other datasets in the same way (e.g.,\n"
00071 " functional datasets).\n"
00072 #if 0
00073 " N.B.: The matrix+vector of the transformation is also\n"
00074 " saved in the .HEAD file of the output dataset\n"
00075 " (unless -dummy is used). You can use the\n"
00076 " output dataset from 3dTagalign as the dataset\n"
00077 " for 3drotates's '-matvec_dset' option; this\n"
00078 " lets you avoid using an intermediate mfile.\n"
00079 "\n"
00080 " -3dWarp = Use the '3dWarp' function to transform the dataset,\n"
00081 " rather than the '3drotate' function. The output\n"
00082 " dataset will be computed on the same grid as the\n"
00083 " master dataset.\n"
00084 " N.B.: This option is implied by '-affine' and '-rotscl'.\n"
00085 #endif
00086 "\n"
00087 " -rotate = Compute the best transformation as a rotation + shift.\n"
00088 " This is the default.\n"
00089 "\n"
00090 " -affine = Compute the best transformation as a general affine\n"
00091 " map rather than just a rotation + shift. In all\n"
00092 " cases, the transformation from input to output\n"
00093 " coordinates is of the form\n"
00094 " [out] = [R] [in] + [V]\n"
00095 " where [R] is a 3x3 matrix and [V] is a 3-vector.\n"
00096 " By default, [R] is computed as a proper (det=1)\n"
00097 " rotation matrix (3 parameters). The '-affine'\n"
00098 " option says to fit [R] as a general matrix\n"
00099 " (9 parameters).\n"
00100 " N.B.: An affine transformation can rotate, rescale, and\n"
00101 " shear the volume. Be sure to look at the dataset\n"
00102 " before and after to make sure things are OK.\n"
00103 "\n"
00104 " -rotscl = Compute transformation as a rotation times an isotropic\n"
00105 " scaling; that is, [R] is an orthogonal matrix times\n"
00106 " a scalar.\n"
00107 " N.B.: '-affine' and '-rotscl' do unweighted least squares.\n"
00108 "\n"
00109 " -prefix pp = Use 'pp' as the prefix for the output dataset.\n"
00110 " [default = 'tagalign']\n"
00111 " -verb = Print progress reports\n"
00112 " -dummy = Don't actually rotate the dataset, just compute\n"
00113 " the transformation matrix and vector. If\n"
00114 " '-matvec' is used, the mfile will be written.\n"
00115 "\n"
00116 "Nota Bene:\n"
00117 "* Cubic interpolation is used. The transformation is carried out\n"
00118 " using the same methods as program 3dWarp.\n"
00119 "\n"
00120 "Author: RWCox - 16 Jul 2000, etc.\n"
00121 ) ;
00122 exit(0) ;
00123 }
00124
00125
00126
00127 iarg = 1 ;
00128 while( iarg < argc && argv[iarg][0] == '-' ){
00129
00130
00131
00132 if( strcmp(argv[iarg],"-rotate") == 0 ){
00133 matrix_type = ROTATION ; use_3dWarp = 1 ;
00134 iarg++ ; continue ;
00135 }
00136
00137
00138
00139 if( strcmp(argv[iarg],"-affine") == 0 ){
00140 matrix_type = AFFINE ; use_3dWarp = 1 ;
00141 iarg++ ; continue ;
00142 }
00143
00144
00145
00146 if( strcmp(argv[iarg],"-rotscl") == 0 ){
00147 matrix_type = ROTSCL ; use_3dWarp = 1 ;
00148 iarg++ ; continue ;
00149 }
00150
00151 #if 0
00152
00153
00154 if( strcmp(argv[iarg],"-3dWarp") == 0 ){
00155 use_3dWarp = 1 ;
00156 iarg++ ; continue ;
00157 }
00158 #endif
00159
00160
00161
00162 if( strcmp(argv[iarg],"-master") == 0 ){
00163 if( mset != NULL ) ERREX("Can only have one -master option") ;
00164 if( ++iarg >= argc ) ERREX("Need an argument after -master") ;
00165
00166 mset = THD_open_dataset( argv[iarg] ) ;
00167
00168 if( mset == NULL ) ERREX("Can't open -master dataset") ;
00169 if( mset->tagset == NULL ) ERREX("No tags in -master dataset") ;
00170 if( TAGLIST_COUNT(mset->tagset) < 3 ) ERREX("Not enough tags in -master dataset") ;
00171
00172 for( nvec=ii=0 ; ii < TAGLIST_COUNT(mset->tagset) ; ii++ )
00173 if( TAG_SET(TAGLIST_SUBTAG(mset->tagset,ii)) ) nvec++ ;
00174
00175 if( nvec < 3 ) ERREX("Not enough tags set in -master dataset") ;
00176
00177 if( nvec < TAGLIST_COUNT(mset->tagset) )
00178 fprintf(stderr,"++ WARNING: not all tags are set in -master dataset\n") ;
00179
00180 if( verb ) fprintf(stderr,"++ Found %d tags in -master dataset\n",nvec) ;
00181
00182 iarg++ ; continue ;
00183 }
00184
00185 #if 0
00186
00187
00188 if( strcmp(argv[iarg],"-wtval") == 0 ){
00189 if( ww != NULL ) ERREX("Can't have -wtval after -wt1D") ;
00190 wtval++ ;
00191 iarg++ ; continue ;
00192 }
00193
00194
00195
00196 if( strcmp(argv[iarg],"-wt1D") == 0 ){
00197 MRI_IMAGE * wtim ; float * wtar ;
00198
00199 if( wtval ) ERREX("Can't have -wt1D after -wtval") ;
00200 if( ww != NULL ) ERREX("Can't have two -wt1D options!") ;
00201 if( ++iarg >= argc ) ERREX("Need an argument after -wt1D") ;
00202
00203 wtim = mri_read_1D( argv[iarg] ) ;
00204
00205 if( wtim == NULL ) ERREX("Can't read -wtim file") ;
00206 if( wtim->ny > 1 ) ERREX("-wtim file has more than one columm") ;
00207
00208 wtar = MRI_FLOAT_PTR(wtim) ;
00209 ww = (double *) malloc(sizeof(double)*wtim->nx) ; nww = wtim->nx ;
00210 for( ii=0 ; ii < nww ; ii++ ){
00211 ww[ii] = (double) wtar[ii] ;
00212 if( ww[ii] < 0.0 ) ERREX("Negative value found in -wt1D file") ;
00213 }
00214
00215 mri_free(wtim) ;
00216 iarg++ ; continue ;
00217 }
00218 #endif
00219
00220
00221
00222 if( strcmp(argv[iarg],"-nokeeptags") == 0 ){
00223 keeptags = 0 ;
00224 iarg++ ; continue ;
00225 }
00226
00227
00228
00229 if( strncmp(argv[iarg],"-verb",5) == 0 ){
00230 verb++ ;
00231 iarg++ ; continue ;
00232 }
00233
00234
00235
00236 if( strcmp(argv[iarg],"-dummy") == 0 ){
00237 dummy++ ;
00238 iarg++ ; continue ;
00239 }
00240
00241
00242
00243 if( strcmp(argv[iarg],"-prefix") == 0 ){
00244 if( ++iarg >= argc ) ERREX("Need an argument after -prefix") ;
00245 prefix = argv[iarg] ;
00246 if( !THD_filename_ok(prefix) ) ERREX("-prefix string is illegal") ;
00247 iarg++ ; continue ;
00248 }
00249
00250
00251
00252 if( strcmp(argv[iarg],"-matvec") == 0 ){
00253 if( ++iarg >= argc ) ERREX("Need an argument after -matvec") ;
00254 mfile = argv[iarg] ;
00255 if( !THD_filename_ok(mfile) ) ERREX("-matvec string is illegal") ;
00256 iarg++ ; continue ;
00257 }
00258
00259
00260
00261
00262 fprintf(stderr,"** Unknown option: %s\n",argv[iarg]) ; exit(1) ;
00263
00264 }
00265
00266 if( mset == NULL ) ERREX("No -master option found on command line") ;
00267
00268 #if 0
00269 if( ww != NULL && nww < nvec ) ERREX("Not enough weights found in -wt1D file") ;
00270
00271
00272
00273 if( wtval ){
00274 ww = (double *) malloc(sizeof(double)*nvec) ; nww = nvec ;
00275 for( ii=jj=0 ; ii < TAGLIST_COUNT(mset->tagset) ; ii++ ){
00276 if( TAG_SET(TAGLIST_SUBTAG(mset->tagset,ii)) ){
00277 ww[jj] = (double) TAG_VAL(TAGLIST_SUBTAG(mset->tagset,ii)) ;
00278
00279 if( ww[jj] < 0.0 ) ERREX("Negative value found in -master tag values") ;
00280 jj++ ;
00281 }
00282 }
00283 }
00284 #endif
00285
00286
00287
00288 if( iarg >= argc ) ERREX("No input dataset?") ;
00289
00290 dset = THD_open_dataset( argv[iarg] ) ;
00291
00292 if( dset == NULL ) ERREX("Can't open input dataset") ;
00293 if( dset->tagset == NULL ) ERREX("No tags in input dataset") ;
00294 if( TAGLIST_COUNT(dset->tagset) !=
00295 TAGLIST_COUNT(mset->tagset) ) ERREX("Tag counts don't match in -master and input") ;
00296
00297
00298
00299 for( ii=0 ; ii < TAGLIST_COUNT(mset->tagset) ; ii++ ){
00300 if( TAG_SET(TAGLIST_SUBTAG(mset->tagset,ii)) !=
00301 TAG_SET(TAGLIST_SUBTAG(dset->tagset,ii)) )
00302 ERREX("Set tags don't match in -master and input") ;
00303 }
00304
00305
00306
00307 xx = (THD_dfvec3 *) malloc( sizeof(THD_dfvec3) * nvec ) ;
00308 yy = (THD_dfvec3 *) malloc( sizeof(THD_dfvec3) * nvec ) ;
00309 dsum = 0.0 ;
00310 for( ii=jj=0 ; ii < nvec ; ii++ ){
00311 if( TAG_SET(TAGLIST_SUBTAG(mset->tagset,ii)) ){
00312
00313 LOAD_DFVEC3( xx[jj] ,
00314 TAG_X( TAGLIST_SUBTAG(mset->tagset,ii) ) ,
00315 TAG_Y( TAGLIST_SUBTAG(mset->tagset,ii) ) ,
00316 TAG_Z( TAGLIST_SUBTAG(mset->tagset,ii) ) ) ;
00317
00318 LOAD_DFVEC3( yy[jj] ,
00319 TAG_X( TAGLIST_SUBTAG(dset->tagset,ii) ) ,
00320 TAG_Y( TAGLIST_SUBTAG(dset->tagset,ii) ) ,
00321 TAG_Z( TAGLIST_SUBTAG(dset->tagset,ii) ) ) ;
00322
00323 dv = SUB_DFVEC3( xx[jj] , yy[jj] ) ;
00324 dsum += dv.xyz[0]*dv.xyz[0] + dv.xyz[1]*dv.xyz[1] + dv.xyz[2]*dv.xyz[2] ;
00325
00326 jj++ ;
00327 }
00328 }
00329
00330 dsum = sqrt(dsum/nvec) ;
00331 fprintf(stderr,"++ RMS distance between tags before = %.2f mm\n" , dsum ) ;
00332
00333
00334
00335 switch( matrix_type ){
00336 default:
00337 case ROTATION:
00338 rt = DLSQ_rot_trans( nvec , yy , xx , ww ) ;
00339 break ;
00340
00341 case AFFINE:
00342 rt = DLSQ_affine ( nvec , yy , xx ) ;
00343 break ;
00344
00345 case ROTSCL:
00346 rt = DLSQ_rotscl ( nvec , yy , xx , (DSET_NZ(dset)==1) ? 2 : 3 ) ;
00347 break ;
00348 }
00349 rtinv = INV_DVECMAT(rt) ;
00350
00351
00352
00353 nval = 0 ;
00354 for( ii=0 ; ii < 3 ; ii++ ){
00355 dsum = rt.vv.xyz[ii] ; nval += thd_floatscan(1,&dsum) ;
00356 for( jj=0 ; jj < 3 ; jj++ ){
00357 dsum = rt.mm.mat[ii][jj] ; nval += thd_floatscan(1,&dsum) ;
00358 }
00359 }
00360 if( nval > 0 ){
00361 fprintf(stderr,"** Floating point errors during calculation\n"
00362 "** of transform matrix and translation vector\n" ) ;
00363 exit(1) ;
00364 }
00365
00366
00367
00368 dsum = DMAT_DET(rt.mm) ;
00369
00370 if( dsum == 0.0 || (matrix_type == ROTATION && fabs(dsum-1.0) > 0.01) ){
00371 fprintf(stderr,"** Invalid transform matrix computed: tags dependent?\n"
00372 "** computed [matrix] and [vector] follow:\n" ) ;
00373
00374 for( ii=0 ; ii < 3 ; ii++ )
00375 fprintf(stderr," [ %10.5f %10.5f %10.5f ] [ %10.5f ] \n",
00376 rt.mm.mat[ii][0],rt.mm.mat[ii][1],rt.mm.mat[ii][2],rt.vv.xyz[ii] );
00377
00378 exit(1) ;
00379 }
00380
00381
00382
00383 if( verb ){
00384 fprintf(stderr,"++ Matrix & Vector [Dicom: x=R-L; y=A-P; z=I-S]\n") ;
00385 for( ii=0 ; ii < 3 ; ii++ )
00386 fprintf(stderr," %10.5f %10.5f %10.5f %10.5f\n",
00387 rt.mm.mat[ii][0],rt.mm.mat[ii][1],rt.mm.mat[ii][2],rt.vv.xyz[ii] );
00388 }
00389
00390 if( matrix_type == ROTATION || matrix_type == ROTSCL ){
00391 double theta, costheta , dist , fac=1.0 ;
00392
00393 if( matrix_type == ROTSCL ){
00394 fac = DMAT_DET(rt.mm); fac = fabs(fac);
00395 if( DSET_NZ(dset) == 1 ) fac = sqrt(fac) ;
00396 else fac = pow(fac,0.33333333);
00397 }
00398
00399 costheta = 0.5 * sqrt(1.0 + DMAT_TRACE(rt.mm)/fac ) ;
00400 theta = 2.0 * acos(costheta) * 180/3.14159265 ;
00401 dist = SIZE_DFVEC3(rt.vv) ;
00402
00403 fprintf(stderr,"++ Total rotation=%.2f degrees; translation=%.2f mm; scaling=%.2f\n",
00404 theta,dist,fac) ;
00405 }
00406
00407 if( mfile ){
00408 FILE * mp ;
00409
00410 if( THD_is_file(mfile) )
00411 fprintf(stderr,"++ Warning: -matvec will overwrite file %s\n",mfile) ;
00412
00413 mp = fopen(mfile,"w") ;
00414 if( mp == NULL ){
00415 fprintf(stderr,"** Can't write to -matvec %s\n",mfile) ;
00416 } else {
00417 for( ii=0 ; ii < 3 ; ii++ )
00418 fprintf(mp," %10.5f %10.5f %10.5f %10.5f\n",
00419 rt.mm.mat[ii][0],rt.mm.mat[ii][1],rt.mm.mat[ii][2],rt.vv.xyz[ii] );
00420 fclose(mp) ;
00421 if( verb ) fprintf(stderr,"++ Wrote matrix+vector to %s\n",mfile) ;
00422 }
00423 }
00424
00425 if( dummy ){
00426 fprintf(stderr,"++ This was a -dummy run: no output dataset\n") ; exit(0) ;
00427 }
00428
00429
00430
00431
00432 #if 0
00433 if( !use_3dWarp ){
00434
00435
00436
00437
00438 pp = DBLE_mat_to_dicomm( dset ) ;
00439 ppt = TRANSPOSE_DMAT(pp) ;
00440 rr = DMAT_MUL(ppt,rt.mm) ; rr = DMAT_MUL(rr,pp) ; tt = DMATVEC(ppt,rt.vv) ;
00441
00442
00443
00444
00445 DSET_mallocize(dset) ;
00446 DSET_load( dset ) ;
00447 if( !DSET_LOADED(dset) ) ERREX("Can't load input dataset bricks") ;
00448 dset->idcode = MCW_new_idcode() ;
00449 dset->dblk->diskptr->storage_mode = STORAGE_BY_BRICK ;
00450 EDIT_dset_items( dset ,
00451 ADN_prefix , prefix ,
00452 ADN_label1 , prefix ,
00453 ADN_none ) ;
00454
00455 if( THD_is_file(dset->dblk->diskptr->header_name) ){
00456 fprintf(stderr,
00457 "** Output file %s already exists -- cannot continue!\n",
00458 dset->dblk->diskptr->header_name ) ;
00459 exit(1) ;
00460 }
00461
00462 tross_Make_History( "3dTagalign" , argc,argv , dset ) ;
00463
00464
00465
00466 if( keeptags ){
00467 THD_dfvec3 rv ;
00468
00469 dsum = 0.0 ;
00470 for( jj=ii=0 ; ii < TAGLIST_COUNT(dset->tagset) ; ii++ ){
00471 if( TAG_SET(TAGLIST_SUBTAG(dset->tagset,ii)) ){
00472 rv = DMATVEC( rt.mm , yy[jj] ) ;
00473 rv = ADD_DFVEC3( rt.vv , rv ) ;
00474
00475 dv = SUB_DFVEC3( xx[jj] , rv ) ;
00476 dsum += dv.xyz[0]*dv.xyz[0] + dv.xyz[1]*dv.xyz[1]
00477 + dv.xyz[2]*dv.xyz[2] ;
00478
00479 UNLOAD_DFVEC3( rv , TAG_X( TAGLIST_SUBTAG(dset->tagset,ii) ) ,
00480 TAG_Y( TAGLIST_SUBTAG(dset->tagset,ii) ) ,
00481 TAG_Z( TAGLIST_SUBTAG(dset->tagset,ii) ) ) ;
00482
00483 jj++ ;
00484 }
00485 }
00486 dsum = sqrt(dsum/nvec) ;
00487 fprintf(stderr,"++ RMS distance between tags after = %.2f mm\n" , dsum ) ;
00488
00489 } else {
00490 myXtFree(dset->tagset) ;
00491 }
00492
00493
00494
00495 if( verb ) fprintf(stderr,"++ computing output BRIK") ;
00496
00497 nvox = DSET_NVOX(dset) ;
00498 nval = DSET_NVALS(dset) ;
00499 fvol = (float *) malloc( sizeof(float) * nvox ) ;
00500
00501 THD_rota_method( MRI_HEPTIC ) ;
00502 clipit = 1 ;
00503
00504 for( ival=0 ; ival < nval ; ival++ ){
00505
00506
00507
00508 EDIT_coerce_type( nvox ,
00509 DSET_BRICK_TYPE(dset,ival),DSET_ARRAY(dset,ival) ,
00510 MRI_float,fvol ) ;
00511
00512 if( clipit ){
00513 register int ii ; register float bb,tt ;
00514 bb = tt = fvol[0] ;
00515 for( ii=1 ; ii < nvox ; ii++ ){
00516 if( fvol[ii] < bb ) bb = fvol[ii] ;
00517 else if( fvol[ii] > tt ) tt = fvol[ii] ;
00518 }
00519 cbot = bb ; ctop = tt ;
00520 }
00521
00522 if( verb && nval < 5 ) fprintf(stderr,".") ;
00523
00524
00525
00526 THD_rota_vol_matvec( DSET_NX(dset) , DSET_NY(dset) , DSET_NZ(dset) ,
00527 fabs(DSET_DX(dset)) , fabs(DSET_DY(dset)) ,
00528 fabs(DSET_DZ(dset)) ,
00529 fvol , rr , tt ) ;
00530
00531 if( verb ) fprintf(stderr,".") ;
00532
00533 if( clipit ){
00534 register int ii ; register float bb,tt ;
00535 bb = cbot ; tt = ctop ;
00536 for( ii=0 ; ii < nvox ; ii++ ){
00537 if( fvol[ii] < bb ) fvol[ii] = bb ;
00538 else if( fvol[ii] > tt ) fvol[ii] = tt ;
00539 }
00540 }
00541
00542 if( verb && nval < 5 ) fprintf(stderr,".") ;
00543
00544
00545
00546 EDIT_coerce_type( nvox, MRI_float,fvol ,
00547 DSET_BRICK_TYPE(dset,ival),DSET_ARRAY(dset,ival) );
00548
00549 }
00550
00551 if( verb ) fprintf(stderr,":") ;
00552
00553
00554
00555 UNLOAD_DMAT(rt.mm,matar[0],matar[1],matar[2],
00556 matar[4],matar[5],matar[6],
00557 matar[8],matar[9],matar[10] ) ;
00558 UNLOAD_DFVEC3(rt.vv,matar[3],matar[7],matar[11]) ;
00559 THD_set_atr( dset->dblk, "TAGALIGN_MATVEC", ATR_FLOAT_TYPE, 12, matar ) ;
00560
00561
00562
00563 dset->dblk->master_nvals = 0 ;
00564 DSET_write(dset) ;
00565
00566 if( verb ) fprintf(stderr,"\n") ;
00567
00568 } else
00569 #endif
00570 {
00571
00572 THD_3dim_dataset *oset ;
00573 THD_vecmat tran ;
00574
00575 #if 0
00576 DFVEC3_TO_FVEC3( rt.vv , tran.vv ) ;
00577 DMAT_TO_MAT ( rt.mm , tran.mm ) ;
00578 #else
00579 DFVEC3_TO_FVEC3( rtinv.vv , tran.vv ) ;
00580 DMAT_TO_MAT ( rtinv.mm , tran.mm ) ;
00581 #endif
00582
00583 mri_warp3D_method( MRI_CUBIC ) ;
00584 oset = THD_warp3D_affine( dset, tran, mset, prefix, 0, WARP3D_NEWDSET ) ;
00585 if( oset == NULL ){
00586 fprintf(stderr,"** ERROR: THD_warp3D() fails!\n"); exit(1);
00587 }
00588
00589 tross_Copy_History( dset , oset ) ;
00590 tross_Make_History( "3dTagalign" , argc,argv , oset ) ;
00591
00592 UNLOAD_DMAT(rt.mm,matar[0],matar[1],matar[2],
00593 matar[4],matar[5],matar[6],
00594 matar[8],matar[9],matar[10] ) ;
00595 UNLOAD_DFVEC3(rt.vv,matar[3],matar[7],matar[11]) ;
00596 THD_set_atr( oset->dblk, "TAGALIGN_MATVEC", ATR_FLOAT_TYPE, 12, matar ) ;
00597
00598
00599
00600 if( keeptags ){
00601 THD_dfvec3 rv ;
00602
00603 oset->tagset = myXtNew(THD_usertaglist) ;
00604 *(oset->tagset) = *(dset->tagset) ;
00605
00606 dsum = 0.0 ;
00607 for( jj=ii=0 ; ii < TAGLIST_COUNT(oset->tagset) ; ii++ ){
00608 if( TAG_SET(TAGLIST_SUBTAG(oset->tagset,ii)) ){
00609 rv = DMATVEC( rt.mm , yy[jj] ) ;
00610 rv = ADD_DFVEC3( rt.vv , rv ) ;
00611
00612 dv = SUB_DFVEC3( xx[jj] , rv ) ;
00613 dsum += dv.xyz[0]*dv.xyz[0] + dv.xyz[1]*dv.xyz[1]
00614 + dv.xyz[2]*dv.xyz[2] ;
00615
00616 UNLOAD_DFVEC3( rv , TAG_X( TAGLIST_SUBTAG(oset->tagset,ii) ) ,
00617 TAG_Y( TAGLIST_SUBTAG(oset->tagset,ii) ) ,
00618 TAG_Z( TAGLIST_SUBTAG(oset->tagset,ii) ) ) ;
00619
00620 jj++ ;
00621 }
00622 }
00623 dsum = sqrt(dsum/nvec) ;
00624 fprintf(stderr,"++ RMS distance between tags after = %.2f mm\n" , dsum ) ;
00625 }
00626
00627 DSET_write(oset) ;
00628
00629 }
00630
00631 exit(0) ;
00632 }
00633
00634
00635
00636 #include "thd.h"
00637
00638
00639
00640
00641
00642
00643 THD_dmat33 DBLE_mat_to_dicomm( THD_3dim_dataset * dset )
00644 {
00645 THD_dmat33 tod ;
00646
00647 LOAD_ZERO_DMAT(tod) ;
00648
00649 switch( dset->daxes->xxorient ){
00650 case ORI_R2L_TYPE: tod.mat[0][0] = 1.0 ; break ;
00651 case ORI_L2R_TYPE: tod.mat[0][0] = -1.0 ; break ;
00652 case ORI_P2A_TYPE: tod.mat[1][0] = -1.0 ; break ;
00653 case ORI_A2P_TYPE: tod.mat[1][0] = 1.0 ; break ;
00654 case ORI_I2S_TYPE: tod.mat[2][0] = 1.0 ; break ;
00655 case ORI_S2I_TYPE: tod.mat[2][0] = -1.0 ; break ;
00656
00657 default: THD_FATAL_ERROR("illegal xxorient code") ;
00658 }
00659
00660 switch( dset->daxes->yyorient ){
00661 case ORI_R2L_TYPE: tod.mat[0][1] = 1.0 ; break ;
00662 case ORI_L2R_TYPE: tod.mat[0][1] = -1.0 ; break ;
00663 case ORI_P2A_TYPE: tod.mat[1][1] = -1.0 ; break ;
00664 case ORI_A2P_TYPE: tod.mat[1][1] = 1.0 ; break ;
00665 case ORI_I2S_TYPE: tod.mat[2][1] = 1.0 ; break ;
00666 case ORI_S2I_TYPE: tod.mat[2][1] = -1.0 ; break ;
00667
00668 default: THD_FATAL_ERROR("illegal yyorient code") ;
00669 }
00670
00671 switch( dset->daxes->zzorient ){
00672 case ORI_R2L_TYPE: tod.mat[0][2] = 1.0 ; break ;
00673 case ORI_L2R_TYPE: tod.mat[0][2] = -1.0 ; break ;
00674 case ORI_P2A_TYPE: tod.mat[1][2] = -1.0 ; break ;
00675 case ORI_A2P_TYPE: tod.mat[1][2] = 1.0 ; break ;
00676 case ORI_I2S_TYPE: tod.mat[2][2] = 1.0 ; break ;
00677 case ORI_S2I_TYPE: tod.mat[2][2] = -1.0 ; break ;
00678
00679 default: THD_FATAL_ERROR("illegal zxorient code") ;
00680 }
00681
00682 return tod ;
00683 }