#include "mrilib.h"
int main( int argc , char *argv[] )
{
THD_3dim_dataset *dsepi=NULL , *dsant=NULL ;
int nopt=1 ;
char *prefix=NULL , str[1024] ;
int nx=1 , ny=5 , nz=0 , ii ;
float step=1.0 ;
int verb=0 ;
THD_fvec3 fv ;
if( argc < 2 || strcmp(argv[1],"-help") == 0 ){
printf("Usage: 3dAnatNudge [options]\n"
"Moves the anat dataset around to best overlap the epi dataset.\n"
"\n"
"OPTIONS:\n"
" -anat aaa = aaa is an 'scalped' (3dIntracranial) high-resolution\n"
" anatomical dataset [a mandatory option]\n"
" -epi eee = eee is an EPI dataset [a mandatory option]\n"
" The first [0] sub-brick from each dataset is used,\n"
" unless otherwise specified on the command line.\n"
" -prefix ppp = ppp is the prefix of the output dataset;\n"
" this dataset will differ from the input only\n"
" in its name and its xyz-axes origin\n"
" [default=don't write new dataset]\n"
" -step sss = set the step size to be sss times the voxel size\n"
" in the anat dataset [default=1.0]\n"
" -x nx = search plus and minus nx steps along the EPI\n"
" -y ny dataset's x-axis; similarly for ny and the\n"
" -z nz y-axis, and for nz and the z-axis\n"
" [default: nx=1 ny=5 nz=0]\n"
" -verb = print progress reports (this is a slow program)\n"
"\n"
"NOTES\n"
"*Systematically moves the anat dataset around and find the shift\n"
" that maximizes overlap between the anat dataset and the EPI\n"
" dataset. No rotations are done.\n"
"*Note that if you use -prefix, a new dataset will be created that\n"
" is a copy of the anat, except that it's origin will be shifted\n"
" and it will have a different ID code than the anat. If you want\n"
" to use this new dataset as the anatomy parent for the EPI\n"
" datasets, you'll have to use\n"
" 3drefit -apar ppp+orig eee1+orig eee2+orig ...\n"
"*If no new dataset is written (no -prefix option), then you\n"
" can use the 3drefit command emitted at the end to modify\n"
" the origin of the anat dataset. (Assuming you trust the\n"
" results - visual inspection is recommended!)\n"
"*The reason the default search grid is mostly along the EPI y-axis\n"
" is that axis is usually the phase-encoding direction, which is\n"
" most subject to displacement due to off-resonance effects.\n"
"*Note that the time this program takes will be proportional to\n"
" (2*nx+1)*(2*ny+1)*(2*nz+1), so using a very large search grid\n"
" will result in a very large usage of CPU time.\n"
"*Recommended usage:\n"
" + Make a 1-brick function volume from a typical EPI dataset:\n"
" 3dbucket -fbuc -prefix epi_fb epi+orig\n"
" + Use 3dIntracranial to scalp a T1-weighted volume:\n"
" 3dIntracranial -anat spgr+orig -prefix spgr_st\n"
" + Use 3dAnatNudge to produce a shifted anat dataset\n"
" 3dAnatNudge -anat spgr_st+orig -epi epi_fb+orig -prefix spgr_nudge\n"
" + Start AFNI and look at epi_fb overlaid in color on the\n"
" anat datasets spgr_st+orig and spgr_nudge+orig, to see if the\n"
" nudged dataset seems like a better fit.\n"
" + Delete the nudged dataset spgr_nudge.\n"
" + If the nudged dataset DOES look better, then apply the\n"
" 3drefit command output by 3dAnatNudge to spgr+orig.\n"
"*Note that the x-, y-, and z-axes for the epi and anat datasets\n"
" may point in different directions (e.g., axial SPGR and\n"
" coronal EPI). The 3drefit command applies to the anat\n"
" dataset, NOT to the EPI dataset.\n"
"*If the program runs successfully, the only thing set to stdout\n"
" will be the 3drefit command string; all other messages go to\n"
" stderr. This can be useful if you want to capture the command\n"
" to a shell variable and then execute it, as in the following\n"
" csh fragment:\n"
" set cvar = `3dAnatNudge ...`\n"
" if( $cvar[1] == \"3drefit\" ) $cvar\n"
" The test on the first sub-string in cvar allows for the\n"
" possibility that the program fails, or that the optimal\n"
" nudge is zero.\n"
) ;
exit(0) ;
}
mainENTRY("3dAnatNudge"); machdep(); AFNI_logger("3dAnatNudge",argc,argv);
PRINT_VERSION("3dAnatNudge");
while( nopt < argc ){
/*-- -anat --*/
if( strcmp(argv[nopt],"-anat") == 0 ){
if( dsant != NULL ){
fprintf(stderr,"** Can't have 2 -anat options!\n"); exit(1);
}
dsant = THD_open_dataset( argv[++nopt] ) ;
CHECK_OPEN_ERROR(dsant,argv[nopt]) ;
nopt++ ; continue ;
}
/*-- -epi --*/
if( strcmp(argv[nopt],"-epi") == 0 ){
if( dsepi != NULL ){
fprintf(stderr,"** Can't have 2 -epi options!\n"); exit(1);
}
dsepi = THD_open_dataset( argv[++nopt] ) ;
CHECK_OPEN_ERROR(dsepi,argv[nopt]) ;
nopt++ ; continue ;
}
/*-- -prefix --*/
if( strcmp(argv[nopt],"-prefix") == 0 ){
prefix = argv[++nopt] ;
if( !THD_filename_ok(prefix) ){
fprintf(stderr,"** Illegal prefix\n"); exit(1);
}
nopt++ ; continue ;
}
/*-- -verb --*/
if( strcmp(argv[nopt],"-verb") == 0 ){
verb = 1 ; nopt++ ; continue ;
}
/*-- -step --*/
if( strcmp(argv[nopt],"-step") == 0 ){
step = strtod( argv[++nopt] , NULL ) ;
if( step <= 0.0 || step > 2.0 ){
fprintf(stderr,"** Illegal value of -step\n"); exit(1);
}
nopt++ ; continue ;
}
/*-- -x or -y or -z --*/
if( strcmp(argv[nopt],"-x") == 0 ){
nx = strtol( argv[++nopt] , NULL , 10 ) ;
if( nx < 0 || nx > 30 ){
fprintf(stderr,"** Illegal value of -x (range 0 30)\n"); exit(1);
}
nopt++ ; continue ;
}
if( strcmp(argv[nopt],"-y") == 0 ){
ny = strtol( argv[++nopt] , NULL , 10 ) ;
if( ny < 0 || ny > 30 ){
fprintf(stderr,"** Illegal value of -y (range 0 30)\n"); exit(1);
}
nopt++ ; continue ;
}
if( strcmp(argv[nopt],"-z") == 0 ){
nz = strtol( argv[++nopt] , NULL , 10 ) ;
if( nz < 0 || nz > 30 ){
fprintf(stderr,"** Illegal value of -z (range 0 30)\n"); exit(1);
}
nopt++ ; continue ;
}
/*-- bad news --*/
fprintf(stderr,"** Illegal option: %s\n",argv[nopt]); exit(1);
}
if( dsant == NULL || dsepi == NULL ){
fprintf(stderr,"** Must use both -anat and -epi options!\n"); exit(1);
}
/* do the actual work! */
fv = THD_autonudge( dsepi,0 , dsant,0 , step , nx,ny,nz , verb ) ;
/* see if we got anything back */
if( fv.xyz[0] == 0.0 && fv.xyz[1] == 0.0 && fv.xyz[2] == 0.0 ){
fprintf(stderr,"# Best nudge is nothing.\n") ; exit(0) ;
}
fprintf(stderr,"++ Best nudge is %.4f%c %.4f%c %.4f%c\n" ,
-fv.xyz[0] , ORIENT_typestr[dsant->daxes->xxorient][0] ,
-fv.xyz[1] , ORIENT_typestr[dsant->daxes->yyorient][0] ,
-fv.xyz[2] , ORIENT_typestr[dsant->daxes->zzorient][0] ) ;
/* print 3drefit command */
strcpy(str,"3drefit") ; ii = strlen(str) ;
if( fv.xyz[0] != 0.0 ){
sprintf(str+ii," -dxorigin %.4f",fv.xyz[0]); ii = strlen(str);
}
if( fv.xyz[1] != 0.0 ){
sprintf(str+ii," -dyorigin %.4f",fv.xyz[1]); ii = strlen(str);
}
if( fv.xyz[2] != 0.0 ){
sprintf(str+ii," -dzorigin %.4f",fv.xyz[2]); ii = strlen(str);
}
sprintf(str+ii," %s",DSET_HEADNAME(dsant)) ;
printf("%s\n",str) ;
/* produce new dataset */
if( prefix != NULL ){
THD_fvec3 ov , nv ;
LOAD_FVEC3(ov,DSET_XORG(dsant),DSET_YORG(dsant),DSET_ZORG(dsant)) ;
nv = ADD_FVEC3(ov,fv) ;
DSET_load(dsant) ;
dsant->dblk->diskptr->storage_mode = STORAGE_BY_BRICK ; /* 14 Jan 2004 */
dsant->idcode = MCW_new_idcode() ;
EDIT_dset_items( dsant ,
ADN_prefix , prefix ,
ADN_xyzorg , nv ,
ADN_none ) ;
if( THD_is_file(DSET_HEADNAME(dsant)) ){
fprintf(stderr,"** Can't overwrite existing dataset %s\n",
DSET_HEADNAME(dsant) ) ;
exit(1) ;
}
tross_Make_History( "3dAnatNudge" , argc,argv , dsant ) ;
DSET_write(dsant) ;
fprintf(stderr,"++ Wrote new dataset %s\n",DSET_BRIKNAME(dsant)) ;
}
exit(0) ;
}
syntax highlighted by Code2HTML, v. 0.9.1