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       /*-- -anat --*/
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       /*-- -epi --*/
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       /*-- -prefix --*/
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       /*-- -verb --*/
00128 
00129       if( strcmp(argv[nopt],"-verb") == 0 ){
00130         verb = 1 ; nopt++ ; continue ;
00131       }
00132 
00133       /*-- -step --*/
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       /*-- -x or -y or -z --*/
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       /*-- bad news --*/
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    /* do the actual work! */
00179 
00180    fv = THD_autonudge( dsepi,0 , dsant,0 , step , nx,ny,nz , verb ) ;
00181 
00182    /* see if we got anything back */
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    /* print 3drefit command */
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    /* produce new dataset */
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 ; /* 14 Jan 2004 */
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 }