00001
00002
00003
00004
00005
00006 #include "mrilib.h"
00007
00008
00009
00010 #define MATVEC_FOR 1
00011 #define MATVEC_BAC 2
00012
00013 static THD_vecmat dicom_in2out , dicom_out2in ;
00014
00015
00016
00017 void warp_dicom_in2out( float xin , float yin , float zin ,
00018 float *xout, float *yout, float *zout )
00019 {
00020 THD_fvec3 xxx ;
00021 LOAD_FVEC3( xxx , xin,yin,zin ) ;
00022 xxx = VECMAT_VEC( dicom_in2out , xxx ) ;
00023 UNLOAD_FVEC3( xxx , *xout , *yout , *zout ) ;
00024 }
00025
00026
00027
00028 void warp_dicom_out2in( float xin , float yin , float zin ,
00029 float *xout, float *yout, float *zout )
00030 {
00031 THD_fvec3 xxx ;
00032 LOAD_FVEC3( xxx , xin,yin,zin ) ;
00033 xxx = VECMAT_VEC( dicom_out2in , xxx ) ;
00034 UNLOAD_FVEC3( xxx , *xout , *yout , *zout ) ;
00035 }
00036
00037
00038
00039 int main( int argc , char * argv[] )
00040 {
00041 THD_3dim_dataset *inset , *outset=NULL , *newgset=NULL ;
00042 int nxin,nyin,nzin,nxyzin,nvals , ival ;
00043 int nxout,nyout,nzout,nxyzout ;
00044 char * prefix = "warped" ;
00045 int nopt=1 , verb=0 , zpad=0 , fsl=0 ;
00046 int use_matvec=0 ;
00047 float xbot,xtop , ybot,ytop , zbot,ztop ;
00048 int use_newgrid=0 ;
00049 float ddd_newgrid=0.0 , fac ;
00050 MRI_IMAGE *inim , *outim , *wim ;
00051 void *newggg=NULL ; int gflag=0 ;
00052
00053 int tta2mni=0 , mni2tta=0 ;
00054 int matdefined=0 ;
00055
00056 THD_coorder cord ;
00057 ATR_float *atr_matfor=NULL , *atr_matinv=NULL ;
00058
00059 #if 0
00060 MRI_IMAGE *matflim=NULL ;
00061 float *matflar=NULL ;
00062 #endif
00063
00064
00065
00066 if( argc < 2 || strcmp(argv[1],"-help") == 0 ){
00067 printf("\n"
00068 "Usage: 3dWarp [options] dataset\n"
00069 "Warp (spatially transform) a 3D dataset.\n"
00070 "--------------------------\n"
00071 "Transform Defining Options: [exactly one of these must be used]\n"
00072 "--------------------------\n"
00073 " -matvec_in2out mmm = Read a 3x4 affine transform matrix+vector\n"
00074 " from file 'mmm':\n"
00075 " x_out = Matrix x_in + Vector\n"
00076 "\n"
00077 " -matvec_out2in mmm = Read a 3x4 affine transform matrix+vector\n"
00078 " from file 'mmm':\n"
00079 " x_in = Matrix x_out + Vector\n"
00080 "\n"
00081 " ** N.B.: The coordinate vectors described above are\n"
00082 " defined in DICOM ('RAI') coordinate order.\n"
00083 " (Also see the '-fsl_matvec option, below.)\n"
00084 " ** N.B.: Using the special name 'IDENTITY' for 'mmm'\n"
00085 " means to use the identity matrix.\n"
00086 " ** N.B.: You can put the matrix on the command line\n"
00087 " directly by using an argument of the form\n"
00088 " 'MATRIX(a11,a12,a13,a14,a21,a22,a23,a24,a31,a32,a33,a34)'\n"
00089 " in place of 'mmm', where the aij values are the\n"
00090 " matrix entries (aij = i-th row, j-th column),\n"
00091 " separated by commas.\n"
00092 " * You will need the 'forward single quotes' around\n"
00093 " the argument.\n"
00094 "\n"
00095 " -tta2mni = Transform a dataset in Talairach-Tournoux Atlas\n"
00096 " coordinates to MNI-152 coordinates.\n"
00097 " -mni2tta = Transform a dataset in MNI-152 coordinates to\n"
00098 " Talairach-Tournoux Atlas coordinates.\n"
00099 "\n"
00100 " -matparent mset = Read in the matrix from WARPDRIVE_MATVEC_*\n"
00101 " attributes in the header of dataset 'mset',\n"
00102 " which must have been created by program\n"
00103 " 3dWarpDrive. In this way, you can apply\n"
00104 " a transformation matrix computed from\n"
00105 " in 3dWarpDrive to another dataset.\n"
00106 "\n"
00107 " ** N.B.: The above option is analogous to the -rotparent\n"
00108 " option in program 3drotate. Use of -matparent\n"
00109 " should be limited to datasets whose spatial\n"
00110 " coordinate system corresponds to that which\n"
00111 " was used for input to 3dWarpDrive (i.e., the\n"
00112 " input to 3dWarp should overlay properly with\n"
00113 " the input to 3dWarpDrive that generated the\n"
00114 " -matparent dataset).\n"
00115 " Sample usages:\n"
00116 " 3dWarpDrive -affine_general -base d1+orig -prefix d2WW -twopass -input d2+orig\n"
00117 " 3dWarp -matparent d2WW+orig -prefix epi2WW epi2+orig\n"
00118 #if 0
00119 "\n"
00120 " -matfile mname = Read in the file 'mname', which consists\n"
00121 " of 12 ASCII-formatted numbers per line.\n"
00122 " The i-th input line defines a 3x4 matrix\n"
00123 " which is used to transform the i-th\n"
00124 " sub-brick of the input dataset.\n"
00125 #endif
00126 "\n"
00127 "-----------------------\n"
00128 "Other Transform Options:\n"
00129 "-----------------------\n"
00130 " -linear }\n"
00131 " -cubic } = Chooses spatial interpolation method.\n"
00132 " -NN } = [default = linear]\n"
00133 " -quintic }\n"
00134 "\n"
00135 " -fsl_matvec = Indicates that the matrix file 'mmm' uses FSL\n"
00136 " ordered coordinates ('LPI'). For use with\n"
00137 " matrix files from FSL and SPM.\n"
00138 "\n"
00139 " -newgrid ddd = Tells program to compute new dataset on a\n"
00140 " new 3D grid, with spacing of 'ddd' mmm.\n"
00141 " * If this option is given, then the new\n"
00142 " 3D region of space covered by the grid\n"
00143 " is computed by warping the 8 corners of\n"
00144 " the input dataset, then laying down a\n"
00145 " regular grid with spacing 'ddd'.\n"
00146 " * If this option is NOT given, then the\n"
00147 " new dataset is computed on the old\n"
00148 " dataset's grid.\n"
00149 "\n"
00150 " -gridset ggg = Tells program to compute new dataset on the\n"
00151 " same grid as dataset 'ggg'.\n"
00152 "\n"
00153 " -zpad N = Tells program to pad input dataset with 'N'\n"
00154 " planes of zeros on all sides before doing\n"
00155 " transformation.\n"
00156 "---------------------\n"
00157 "Miscellaneous Options:\n"
00158 "---------------------\n"
00159 #if 0
00160 " -verb = Print out some information along the way.\n"
00161 #endif
00162 " -prefix ppp = Sets the prefix of the output dataset.\n"
00163 "\n"
00164 ) ;
00165 exit(0) ;
00166 }
00167
00168
00169
00170 mainENTRY("3dWarp main"); machdep(); AFNI_logger("3dWarp",argc,argv);
00171 PRINT_VERSION("3dWarp") ;
00172
00173 THD_coorder_fill( my_getenv("AFNI_ORIENT") , &cord ) ;
00174
00175
00176
00177 while( nopt < argc && argv[nopt][0] == '-' ){
00178
00179
00180
00181 #if 0
00182 if( strncmp(argv[nopt],"-matfile",8) == 0 ){
00183 if( matdefined )
00184 ERROR_exit("-matfile: Matrix already defined!\n");
00185 if( ++nopt >= argc )
00186 ERROR_exit("Need argument after -matfile!\n");
00187 matflim = mri_read_ascii( argv[nopt] ) ;
00188 if( matflim == NULL )
00189 ERROR_exit("Can't read from -matfile %s\n",argv[nopt]);
00190 if( matflim->nx < 12 )
00191 ERROR_exit("Need 12 numbers per line in -matfile %s\n",argv[nopt]) ;
00192
00193 matflar = MRI_FLOAT_PTR(matflim) ;
00194 matdefined = 1 ; nopt++ ; continue ;
00195 }
00196 #endif
00197
00198
00199
00200 if( strncmp(argv[nopt],"-matparent",10) == 0 ){
00201 ATR_float *atr ; float *matar ; THD_3dim_dataset *matparset=NULL ;
00202
00203 if( matdefined )
00204 ERROR_exit("-matparent: Matrix already defined!\n");
00205 if( ++nopt >= argc )
00206 ERROR_exit("Need argument after -matparent!\n");
00207 matparset = THD_open_dataset( argv[nopt] ) ;
00208 if( matparset == NULL )
00209 ERROR_exit("Can't open -matparent %s\n",argv[nopt]);
00210 atr = THD_find_float_atr( matparset->dblk, "WARPDRIVE_MATVEC_INV_000000" );
00211 if( atr != NULL ) atr_matinv = (ATR_float *)THD_copy_atr( (ATR_any *)atr ) ;
00212 atr = THD_find_float_atr( matparset->dblk, "WARPDRIVE_MATVEC_FOR_000000" );
00213 if( atr != NULL ) atr_matfor = (ATR_float *)THD_copy_atr( (ATR_any *)atr ) ;
00214 if( atr_matinv == NULL || atr_matinv->nfl < 12 )
00215 ERROR_exit("-matparent %s doesn't have WARPDRIVE attributes!?\n",argv[nopt]) ;
00216
00217 matar = atr_matinv->fl ;
00218 if( strstr(argv[nopt-1],"INV") == NULL ){
00219 LOAD_MAT ( dicom_in2out.mm, matar[0],matar[1],matar[2],
00220 matar[4],matar[5],matar[6],
00221 matar[8],matar[9],matar[10] ) ;
00222 LOAD_FVEC3( dicom_in2out.vv, matar[3],matar[7],matar[11] ) ;
00223 dicom_out2in = INV_VECMAT( dicom_in2out ) ;
00224 } else {
00225 LOAD_MAT ( dicom_out2in.mm, matar[0],matar[1],matar[2],
00226 matar[4],matar[5],matar[6],
00227 matar[8],matar[9],matar[10] ) ;
00228 LOAD_FVEC3( dicom_out2in.vv, matar[3],matar[7],matar[11] ) ;
00229 dicom_in2out = INV_VECMAT( dicom_out2in ) ;
00230 }
00231
00232 DSET_delete(matparset) ;
00233
00234 use_matvec = MATVEC_FOR ; matdefined = 1 ; nopt++ ; continue ;
00235 }
00236
00237
00238
00239 if( strcmp(argv[nopt],"-tta2mni") == 0 ){
00240 if( matdefined )
00241 ERROR_exit("-tta2mni: Matrix already defined!\n");
00242 matdefined = 1 ; tta2mni = 1 ; nopt++ ; continue ;
00243 }
00244
00245 if( strcmp(argv[nopt],"-mni2tta") == 0 ){
00246 if( matdefined )
00247 ERROR_exit("-mni2tta: Matrix already defined!\n");
00248 matdefined = 1 ; mni2tta = 1 ; nopt++ ; continue ;
00249 }
00250
00251
00252
00253 if( strcmp(argv[nopt],"-fsl_matvec") == 0 ){
00254 fsl = 1 ; nopt++ ; continue ;
00255 }
00256
00257
00258
00259 if( strcmp(argv[nopt],"-gridset") == 0 ){
00260 if( use_newgrid )
00261 ERROR_exit("Can't use -gridset with -newgrid!\n");
00262 if( newgset != NULL )
00263 ERROR_exit("Can't use -gridset twice!\n");
00264 if( ++nopt >= argc )
00265 ERROR_exit("Need argument after -gridset!\n");
00266 newgset = THD_open_dataset( argv[nopt] ) ;
00267 if( newgset == NULL )
00268 ERROR_exit("Can't open -gridset %s\n",argv[nopt]);
00269 nopt++ ; continue ;
00270 }
00271
00272
00273
00274 if( strcmp(argv[nopt],"-zpad") == 0 ){
00275 if( ++nopt >= argc )
00276 ERROR_exit("Need argument after -zpad!\n");
00277 zpad = (int) strtod( argv[nopt] , NULL ) ;
00278 if( zpad < 0 )
00279 ERROR_exit("Illegal negative argument after -zpad!\n");
00280 else if( zpad == 0 )
00281 INFO_message("NOTA BENE: -zpad is 0 ==> ignored\n") ;
00282 nopt++ ; continue ;
00283 }
00284
00285
00286
00287 if( strcmp(argv[nopt],"-newgrid") == 0 ){
00288 if( newgset != NULL )
00289 ERROR_exit("Can't use -newgrid with -gridset!\n");
00290 if( use_newgrid )
00291 ERROR_exit("Can't use -newgrid twice!\n");
00292 if( ++nopt >= argc )
00293 ERROR_exit("Need argument after -newgrid!\n");
00294
00295 ddd_newgrid = strtod( argv[nopt] , NULL ) ;
00296 if( ddd_newgrid <= 0.0 )
00297 ERROR_exit("Illegal argument after -newgrid!\n");
00298
00299 use_newgrid = 1 ; nopt++ ; continue ;
00300 }
00301
00302
00303
00304 if( strcmp(argv[nopt],"-NN") == 0 ){
00305 mri_warp3D_method( MRI_NN ) ; nopt++ ; continue ;
00306 }
00307 if( strcmp(argv[nopt],"-linear") == 0 ){
00308 mri_warp3D_method( MRI_LINEAR ) ; nopt++ ; continue ;
00309 }
00310 if( strcmp(argv[nopt],"-cubic") == 0 ){
00311 mri_warp3D_method( MRI_CUBIC ) ; nopt++ ; continue ;
00312 }
00313 if( strcmp(argv[nopt],"-quintic") == 0 ){
00314 mri_warp3D_method( MRI_QUINTIC ); nopt++ ; continue ;
00315 }
00316
00317
00318
00319 if( strcmp(argv[nopt],"-prefix") == 0 ){
00320 if( ++nopt >= argc )
00321 ERROR_exit("Need an argument after -prefix!\n");
00322 if( !THD_filename_ok(argv[nopt]) )
00323 ERROR_exit("-prefix argument is invalid!\n");
00324
00325 prefix = argv[nopt] ; nopt++ ; continue ;
00326 }
00327
00328
00329
00330 if( strncmp(argv[nopt],"-verbose",5) == 0 ){
00331 verb++ ; nopt++ ; continue ;
00332 }
00333
00334
00335
00336 if( strncmp(argv[nopt],"-matvec_",8) == 0 ){
00337 MRI_IMAGE *matim ; float *matar , dm ;
00338
00339 if( matdefined )
00340 ERROR_exit("%s: matrix already defined!\n",argv[nopt]);
00341 if( ++nopt >= argc )
00342 ERROR_exit("Need an argument after -matvec!\n");
00343 if( strcmp(argv[nopt],"IDENTITY") == 0 ){
00344 matim = mri_new( 4,3 , MRI_float ) ;
00345 matar = MRI_FLOAT_PTR(matim) ;
00346 matar[0] = matar[5] = matar[10] = 1.0 ;
00347 } else if( strncmp(argv[nopt],"MATRIX(",7) == 0 ){
00348 int nn ;
00349 matim = mri_new( 4,3 , MRI_float ) ;
00350 matar = MRI_FLOAT_PTR(matim) ;
00351 nn = sscanf(argv[nopt],"MATRIX(%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,%f",
00352 matar+0 , matar+1 , matar+2 , matar+3 ,
00353 matar+4 , matar+5 , matar+6 , matar+7 ,
00354 matar+8 , matar+9 , matar+10, matar+11 ) ;
00355 if( nn < 12 )
00356 ERROR_exit("Could only decode %d numbers from %s\n",nn,argv[nopt]) ;
00357 } else {
00358 matim = mri_read_ascii( argv[nopt] ) ;
00359 if( matim == NULL )
00360 ERROR_exit("Can't read -matvec file %s\n",argv[nopt]);
00361 if( matim->nx != 4 || matim->ny < 3 )
00362 ERROR_exit("-matvec file not a 3x4 matrix!\n");
00363 matar = MRI_FLOAT_PTR(matim) ;
00364 }
00365
00366 use_matvec = (strstr(argv[nopt-1],"_in2out") != NULL) ? MATVEC_FOR : MATVEC_BAC ;
00367
00368 switch( use_matvec ){
00369
00370 case MATVEC_FOR:
00371 LOAD_MAT ( dicom_in2out.mm, matar[0],matar[1],matar[2],
00372 matar[4],matar[5],matar[6],
00373 matar[8],matar[9],matar[10] ) ;
00374 LOAD_FVEC3( dicom_in2out.vv, matar[3],matar[7],matar[11] ) ;
00375
00376 dm = MAT_DET( dicom_in2out.mm ) ;
00377 if( dm == 0.0 )
00378 ERROR_exit("Determinant of matrix is 0\n");
00379
00380 dicom_out2in = INV_VECMAT( dicom_in2out ) ;
00381 break ;
00382
00383 case MATVEC_BAC:
00384 LOAD_MAT ( dicom_out2in.mm, matar[0],matar[1],matar[2],
00385 matar[4],matar[5],matar[6],
00386 matar[8],matar[9],matar[10] ) ;
00387 LOAD_FVEC3( dicom_out2in.vv, matar[3],matar[7],matar[11] ) ;
00388
00389 dm = MAT_DET( dicom_out2in.mm ) ;
00390 if( dm == 0.0 )
00391 ERROR_exit("Determinant of matrix is 0\n");
00392
00393 dicom_in2out = INV_VECMAT( dicom_out2in ) ;
00394 break ;
00395 }
00396
00397 matdefined = 1 ; nopt++ ; continue ;
00398 }
00399
00400
00401
00402 ERROR_exit("Unknown option %s\n",argv[nopt]) ;
00403
00404 }
00405
00406
00407
00408 if( !matdefined )
00409 ERROR_exit("No transformation on command line!?\n");
00410
00411 if( tta2mni || mni2tta ) fsl = 0 ;
00412
00413
00414
00415 if( use_matvec && fsl ){
00416 dicom_in2out.mm.mat[0][2] = -dicom_in2out.mm.mat[0][2] ;
00417 dicom_in2out.mm.mat[1][2] = -dicom_in2out.mm.mat[1][2] ;
00418 dicom_in2out.mm.mat[2][0] = -dicom_in2out.mm.mat[2][0] ;
00419 dicom_in2out.mm.mat[2][1] = -dicom_in2out.mm.mat[2][1] ;
00420 dicom_in2out.vv.xyz[0] = -dicom_in2out.vv.xyz[0] ;
00421 dicom_in2out.vv.xyz[1] = -dicom_in2out.vv.xyz[1] ;
00422
00423 dicom_out2in.mm.mat[0][2] = -dicom_out2in.mm.mat[0][2] ;
00424 dicom_out2in.mm.mat[1][2] = -dicom_out2in.mm.mat[1][2] ;
00425 dicom_out2in.mm.mat[2][0] = -dicom_out2in.mm.mat[2][0] ;
00426 dicom_out2in.mm.mat[2][1] = -dicom_out2in.mm.mat[2][1] ;
00427 dicom_out2in.vv.xyz[0] = -dicom_out2in.vv.xyz[0] ;
00428 dicom_out2in.vv.xyz[1] = -dicom_out2in.vv.xyz[1] ;
00429 }
00430
00431
00432
00433 if( nopt != argc-1 ){
00434 ERROR_message("Command line should have exactly 1 dataset!") ;
00435 ERROR_exit ("Whereas there seems to be %d of them!\n", argc-nopt ) ;
00436 }
00437
00438
00439
00440 inset = THD_open_dataset( argv[nopt] ) ;
00441 if( !ISVALID_DSET(inset) )
00442 ERROR_exit("Can't open dataset %s\n",argv[nopt]);
00443
00444
00445
00446 if( tta2mni || mni2tta ){
00447 if( inset->view_type != VIEW_TALAIRACH_TYPE ){
00448 WARNING_message("Input dataset not in +tlrc view!\n") ;
00449 } else {
00450 ATR_string *ast ;
00451 ast = THD_find_string_atr( inset->dblk , "TLRC_SPACE" ) ;
00452 if( ast != NULL ){
00453 if( strcmp(ast->ch,"MNI-152") == 0 && tta2mni ){
00454 WARNING_message("Input dataset appears to be MNI-152 already!\n");
00455 }
00456 }
00457 }
00458 }
00459
00460
00461
00462 if( use_newgrid ){
00463 newggg = &ddd_newgrid ; gflag = WARP3D_NEWGRID ;
00464 } else if( newgset != NULL ){
00465 newggg = newgset ; gflag = WARP3D_NEWDSET ;
00466 }
00467
00468 if( use_matvec ){
00469 outset = THD_warp3D_affine( inset , dicom_out2in ,
00470 newggg , prefix , zpad , gflag ) ;
00471 } else if( tta2mni ){
00472 outset = THD_warp3D_tta2mni( inset , newggg , prefix,zpad,gflag ) ;
00473 } else if( mni2tta ){
00474 outset = THD_warp3D_mni2tta( inset , newggg , prefix,zpad,gflag ) ;
00475 }
00476
00477 if( outset == NULL )
00478 ERROR_exit("THD_warp3D() fails for some reason!\n") ;
00479
00480
00481
00482 tross_Copy_History( inset , outset ) ;
00483 tross_Make_History( "3dWarp" , argc,argv , outset ) ;
00484
00485
00486
00487 if( tta2mni ){
00488 THD_set_string_atr( outset->dblk , "TLRC_SPACE" , "MNI-152" ) ;
00489 }
00490
00491
00492
00493 if( atr_matfor != NULL ) THD_insert_atr( outset->dblk, (ATR_any *)atr_matfor );
00494 if( atr_matinv != NULL ) THD_insert_atr( outset->dblk, (ATR_any *)atr_matinv );
00495
00496
00497
00498 DSET_delete( inset ) ;
00499 DSET_write( outset ) ; if( verb ) WROTE_DSET(outset) ;
00500 exit(0) ;
00501 }