00001 #include "mrilib.h"
00002
00003 int main( int argc , char *argv[] )
00004 {
00005 THD_3dim_dataset *dsepi=NULL , *dsant=NULL ;
00006 int nopt=1 ;
00007 char *prefix=NULL , str[1024] ;
00008 int nx=1 , ny=5 , nz=0 , ii ;
00009 float step=1.0 ;
00010 int verb=0 ;
00011 THD_fvec3 fv ;
00012
00013 if( argc < 2 || strcmp(argv[1],"-help") == 0 ){
00014 printf("Usage: 3dAnatNudge [options]\n"
00015 "Moves the anat dataset around to best overlap the epi dataset.\n"
00016 "\n"
00017 "OPTIONS:\n"
00018 " -anat aaa = aaa is an 'scalped' (3dIntracranial) high-resolution\n"
00019 " anatomical dataset [a mandatory option]\n"
00020 " -epi eee = eee is an EPI dataset [a mandatory option]\n"
00021 " The first [0] sub-brick from each dataset is used,\n"
00022 " unless otherwise specified on the command line.\n"
00023 " -prefix ppp = ppp is the prefix of the output dataset;\n"
00024 " this dataset will differ from the input only\n"
00025 " in its name and its xyz-axes origin\n"
00026 " [default=don't write new dataset]\n"
00027 " -step sss = set the step size to be sss times the voxel size\n"
00028 " in the anat dataset [default=1.0]\n"
00029 " -x nx = search plus and minus nx steps along the EPI\n"
00030 " -y ny dataset's x-axis; similarly for ny and the\n"
00031 " -z nz y-axis, and for nz and the z-axis\n"
00032 " [default: nx=1 ny=5 nz=0]\n"
00033 " -verb = print progress reports (this is a slow program)\n"
00034 "\n"
00035 "NOTES\n"
00036 "*Systematically moves the anat dataset around and find the shift\n"
00037 " that maximizes overlap between the anat dataset and the EPI\n"
00038 " dataset. No rotations are done.\n"
00039 "*Note that if you use -prefix, a new dataset will be created that\n"
00040 " is a copy of the anat, except that it's origin will be shifted\n"
00041 " and it will have a different ID code than the anat. If you want\n"
00042 " to use this new dataset as the anatomy parent for the EPI\n"
00043 " datasets, you'll have to use\n"
00044 " 3drefit -apar ppp+orig eee1+orig eee2+orig ...\n"
00045 "*If no new dataset is written (no -prefix option), then you\n"
00046 " can use the 3drefit command emitted at the end to modify\n"
00047 " the origin of the anat dataset. (Assuming you trust the\n"
00048 " results - visual inspection is recommended!)\n"
00049 "*The reason the default search grid is mostly along the EPI y-axis\n"
00050 " is that axis is usually the phase-encoding direction, which is\n"
00051 " most subject to displacement due to off-resonance effects.\n"
00052 "*Note that the time this program takes will be proportional to\n"
00053 " (2*nx+1)*(2*ny+1)*(2*nz+1), so using a very large search grid\n"
00054 " will result in a very large usage of CPU time.\n"
00055 "*Recommended usage:\n"
00056 " + Make a 1-brick function volume from a typical EPI dataset:\n"
00057 " 3dbucket -fbuc -prefix epi_fb epi+orig\n"
00058 " + Use 3dIntracranial to scalp a T1-weighted volume:\n"
00059 " 3dIntracranial -anat spgr+orig -prefix spgr_st\n"
00060 " + Use 3dAnatNudge to produce a shifted anat dataset\n"
00061 " 3dAnatNudge -anat spgr_st+orig -epi epi_fb+orig -prefix spgr_nudge\n"
00062 " + Start AFNI and look at epi_fb overlaid in color on the\n"
00063 " anat datasets spgr_st+orig and spgr_nudge+orig, to see if the\n"
00064 " nudged dataset seems like a better fit.\n"
00065 " + Delete the nudged dataset spgr_nudge.\n"
00066 " + If the nudged dataset DOES look better, then apply the\n"
00067 " 3drefit command output by 3dAnatNudge to spgr+orig.\n"
00068 "*Note that the x-, y-, and z-axes for the epi and anat datasets\n"
00069 " may point in different directions (e.g., axial SPGR and\n"
00070 " coronal EPI). The 3drefit command applies to the anat\n"
00071 " dataset, NOT to the EPI dataset.\n"
00072 "*If the program runs successfully, the only thing set to stdout\n"
00073 " will be the 3drefit command string; all other messages go to\n"
00074 " stderr. This can be useful if you want to capture the command\n"
00075 " to a shell variable and then execute it, as in the following\n"
00076 " csh fragment:\n"
00077 " set cvar = `3dAnatNudge ...`\n"
00078 " if( $cvar[1] == \"3drefit\" ) $cvar\n"
00079 " The test on the first sub-string in cvar allows for the\n"
00080 " possibility that the program fails, or that the optimal\n"
00081 " nudge is zero.\n"
00082 ) ;
00083 exit(0) ;
00084 }
00085
00086 mainENTRY("3dAnatNudge"); machdep(); AFNI_logger("3dAnatNudge",argc,argv);
00087 PRINT_VERSION("3dAnatNudge");
00088
00089 while( nopt < argc ){
00090
00091
00092
00093 if( strcmp(argv[nopt],"-anat") == 0 ){
00094 if( dsant != NULL ){
00095 fprintf(stderr,"** Can't have 2 -anat options!\n"); exit(1);
00096 }
00097 dsant = THD_open_dataset( argv[++nopt] ) ;
00098 if( !ISVALID_DSET(dsant) ){
00099 fprintf(stderr,"** Can't open dataset %s\n",argv[nopt]); exit(1);
00100 }
00101 nopt++ ; continue ;
00102 }
00103
00104
00105
00106 if( strcmp(argv[nopt],"-epi") == 0 ){
00107 if( dsepi != NULL ){
00108 fprintf(stderr,"** Can't have 2 -epi options!\n"); exit(1);
00109 }
00110 dsepi = THD_open_dataset( argv[++nopt] ) ;
00111 if( !ISVALID_DSET(dsepi) ){
00112 fprintf(stderr,"** Can't open dataset %s\n",argv[nopt]); exit(1);
00113 }
00114 nopt++ ; continue ;
00115 }
00116
00117
00118
00119 if( strcmp(argv[nopt],"-prefix") == 0 ){
00120 prefix = argv[++nopt] ;
00121 if( !THD_filename_ok(prefix) ){
00122 fprintf(stderr,"** Illegal prefix\n"); exit(1);
00123 }
00124 nopt++ ; continue ;
00125 }
00126
00127
00128
00129 if( strcmp(argv[nopt],"-verb") == 0 ){
00130 verb = 1 ; nopt++ ; continue ;
00131 }
00132
00133
00134
00135 if( strcmp(argv[nopt],"-step") == 0 ){
00136 step = strtod( argv[++nopt] , NULL ) ;
00137 if( step <= 0.0 || step > 2.0 ){
00138 fprintf(stderr,"** Illegal value of -step\n"); exit(1);
00139 }
00140 nopt++ ; continue ;
00141 }
00142
00143
00144
00145 if( strcmp(argv[nopt],"-x") == 0 ){
00146 nx = strtol( argv[++nopt] , NULL , 10 ) ;
00147 if( nx < 0 || nx > 10 ){
00148 fprintf(stderr,"** Illegal value of -x\n"); exit(1);
00149 }
00150 nopt++ ; continue ;
00151 }
00152
00153 if( strcmp(argv[nopt],"-y") == 0 ){
00154 ny = strtol( argv[++nopt] , NULL , 10 ) ;
00155 if( ny < 0 || ny > 10 ){
00156 fprintf(stderr,"** Illegal value of -y\n"); exit(1);
00157 }
00158 nopt++ ; continue ;
00159 }
00160
00161 if( strcmp(argv[nopt],"-z") == 0 ){
00162 nz = strtol( argv[++nopt] , NULL , 10 ) ;
00163 if( nz < 0 || nz > 10 ){
00164 fprintf(stderr,"** Illegal value of -z\n"); exit(1);
00165 }
00166 nopt++ ; continue ;
00167 }
00168
00169
00170
00171 fprintf(stderr,"** Illegal option: %s\n",argv[nopt]); exit(1);
00172 }
00173
00174 if( dsant == NULL || dsepi == NULL ){
00175 fprintf(stderr,"** Must use both -anat and -epi options!\n"); exit(1);
00176 }
00177
00178
00179
00180 fv = THD_autonudge( dsepi,0 , dsant,0 , step , nx,ny,nz , verb ) ;
00181
00182
00183
00184 if( fv.xyz[0] == 0.0 && fv.xyz[1] == 0.0 && fv.xyz[2] == 0.0 ){
00185 fprintf(stderr,"# Best nudge is nothing.\n") ; exit(0) ;
00186 }
00187
00188 fprintf(stderr,"++ Best nudge is %.4f%c %.4f%c %.4f%c\n" ,
00189 -fv.xyz[0] , ORIENT_typestr[dsant->daxes->xxorient][0] ,
00190 -fv.xyz[1] , ORIENT_typestr[dsant->daxes->yyorient][0] ,
00191 -fv.xyz[2] , ORIENT_typestr[dsant->daxes->zzorient][0] ) ;
00192
00193
00194
00195 strcpy(str,"3drefit") ; ii = strlen(str) ;
00196 if( fv.xyz[0] != 0.0 ){
00197 sprintf(str+ii," -dxorigin %.4f",fv.xyz[0]); ii = strlen(str);
00198 }
00199 if( fv.xyz[1] != 0.0 ){
00200 sprintf(str+ii," -dyorigin %.4f",fv.xyz[1]); ii = strlen(str);
00201 }
00202 if( fv.xyz[2] != 0.0 ){
00203 sprintf(str+ii," -dzorigin %.4f",fv.xyz[2]); ii = strlen(str);
00204 }
00205 sprintf(str+ii," %s",DSET_HEADNAME(dsant)) ;
00206 printf("%s\n",str) ;
00207
00208
00209
00210 if( prefix != NULL ){
00211 THD_fvec3 ov , nv ;
00212 LOAD_FVEC3(ov,DSET_XORG(dsant),DSET_YORG(dsant),DSET_ZORG(dsant)) ;
00213 nv = ADD_FVEC3(ov,fv) ;
00214 DSET_load(dsant) ;
00215 dsant->dblk->diskptr->storage_mode = STORAGE_BY_BRICK ;
00216 dsant->idcode = MCW_new_idcode() ;
00217 EDIT_dset_items( dsant ,
00218 ADN_prefix , prefix ,
00219 ADN_xyzorg , nv ,
00220 ADN_none ) ;
00221 if( THD_is_file(DSET_HEADNAME(dsant)) ){
00222 fprintf(stderr,"** Can't overwrite existing dataset %s\n",
00223 DSET_HEADNAME(dsant) ) ;
00224 exit(1) ;
00225 }
00226 tross_Make_History( "3dAnatNudge" , argc,argv , dsant ) ;
00227 DSET_write(dsant) ;
00228 fprintf(stderr,"++ Wrote new dataset %s\n",DSET_BRIKNAME(dsant)) ;
00229 }
00230
00231 exit(0) ;
00232 }