00001 #include "mrilib.h"
00002 #include "rickr/r_new_resam_dset.h"
00003 
00004 int main( int argc , char *argv[] )
00005 {
00006    MRI_IMAGE *imin, *imout , *imout_orig;
00007    THD_3dim_dataset *iset, *oset , *ooset;
00008    char *prefix = "SpatNorm" ;
00009    int iarg , verb=0, OrigSpace = 0 ;
00010    THD_ivec3 orixyz , nxyz ;
00011    THD_fvec3 dxyz , orgxyz, originRAIfv, fv2;
00012 
00013 
00014    /*--- get help here or get help somewhere ---*/
00015 
00016    if( argc < 2 ){
00017      printf("Usage: 3dSpatNorm [options] dataset\n"
00018             "\n"
00019             "Options:\n"
00020             "  -prefix ppp = Write output dataset using 'ppp' for the prefix.\n"
00021             "  -orig_space = Write output dataset using the same grid as dataset.\n"
00022             "  -verb       = Write out progress reports\n"
00023            ) ;
00024      exit(0) ;
00025    }
00026 
00027    /*--- options ---*/
00028 
00029    iarg = 1 ;
00030    OrigSpace = 0;
00031    while( iarg < argc && argv[iarg][0] == '-' ){
00032 
00033      /* -prefix */
00034 
00035      if( strcmp(argv[iarg],"-prefix") == 0 ){
00036        if( ++iarg >= argc ){
00037          fprintf(stderr,"**ERROR: -prefix requires another argument!\n") ;
00038          exit(1) ;
00039        }
00040        prefix = strdup(argv[iarg]) ;
00041        if( !THD_filename_ok(prefix) ){
00042          fprintf(stderr,"**ERROR: -prefix value contains forbidden characters!\n") ;
00043          exit(1) ;
00044        }
00045        iarg++ ; continue ;
00046      }
00047 
00048      /* -verb */
00049 
00050      if( strncmp(argv[iarg],"-verb",5) == 0 ){
00051        verb++ ; iarg++ ; continue ;
00052      }
00053      
00054      if( strncmp(argv[iarg],"-orig_space",10) == 0 ){
00055        OrigSpace = 1 ; iarg++ ; continue ;
00056      }
00057      
00058      fprintf(stderr,"**ERROR: %s is unknown option!\n",argv[iarg]) ;
00059      exit(1) ;
00060    }
00061 
00062    if( iarg >= argc ){
00063      fprintf(stderr,"**ERROR: no input dataset name on command line?!\n") ;
00064      exit(1) ;
00065    }
00066 
00067    /*--- read dataset ---*/
00068 
00069    iset = THD_open_dataset( argv[iarg] ) ;
00070    if( !ISVALID_DSET(iset) ){
00071      fprintf(stderr,"**ERROR: can't open dataset %s\n",argv[iarg]) ;
00072      exit(1) ;
00073    }
00074 
00075    /*--- get median brick --*/
00076 
00077    if( verb ) fprintf(stderr,"++3dSpatNorm: loading dataset\n") ;
00078 
00079    imin = THD_median_brick( iset ) ;
00080    if( imin == NULL ){
00081      fprintf(stderr,"**ERROR: can't load dataset %s\n",argv[iarg]) ;
00082      exit(1) ;
00083    }
00084    imin->dx = fabs(iset->daxes->xxdel) ;
00085    imin->dy = fabs(iset->daxes->yydel) ;
00086    imin->dz = fabs(iset->daxes->zzdel) ;
00087    
00088    mri_brainormalize_initialize(imin->dz, imin->dy, imin->dz); /* To get around the #define for voxel counts and dimensions */
00089    
00090    /* me needs the origin of this dset in RAI world */
00091    LOAD_FVEC3(originRAIfv , iset->daxes->xxorg , iset->daxes->yyorg , iset->daxes->zzorg) ;
00092    originRAIfv = THD_3dmm_to_dicomm( iset , originRAIfv ) ;
00093 
00094    LOAD_FVEC3(fv2 , iset->daxes->xxorg + (iset->daxes->nxx-1)*iset->daxes->xxdel ,
00095                     iset->daxes->yyorg + (iset->daxes->nyy-1)*iset->daxes->yydel ,
00096                     iset->daxes->zzorg + (iset->daxes->nzz-1)*iset->daxes->zzdel  ) ;
00097    fv2 = THD_3dmm_to_dicomm( iset , fv2 ) ;
00098 
00099    if( originRAIfv.xyz[0] > fv2.xyz[0] ) { float tf; tf = originRAIfv.xyz[0]; originRAIfv.xyz[0] = fv2.xyz[0]; fv2.xyz[0] = tf; } 
00100    if( originRAIfv.xyz[1] > fv2.xyz[1] ) { float tf; tf = originRAIfv.xyz[1]; originRAIfv.xyz[1] = fv2.xyz[1]; fv2.xyz[1] = tf; }
00101    if( originRAIfv.xyz[2] > fv2.xyz[2] ) { float tf; tf = originRAIfv.xyz[2]; originRAIfv.xyz[2] = fv2.xyz[2]; fv2.xyz[2] = tf; }
00102    
00103    if (1 && verb) {
00104       fprintf(stderr,"++3dSpatNorm (ZSS): RAI origin info: %f %f %f\n", originRAIfv.xyz[0], originRAIfv.xyz[1], originRAIfv.xyz[2]);
00105    }
00106    
00107    
00108    DSET_unload( iset ) ;  /* don't need this data no more */
00109 
00110    /*-- convert image to shorts, if appropriate --*/
00111 
00112    if( DSET_BRICK_TYPE(iset,0) == MRI_short ||
00113        DSET_BRICK_TYPE(iset,0) == MRI_byte    ){
00114 
00115      imout = mri_to_short(1.0,imin) ;
00116      mri_free(imin) ; imin = imout ;
00117    }
00118 
00119    /*--- normalize image spatially ---*/
00120 
00121    mri_brainormalize_verbose( verb ) ;
00122    if (OrigSpace) {
00123       imout = mri_brainormalize( imin , iset->daxes->xxorient,
00124                                      iset->daxes->yyorient,
00125                                      iset->daxes->zzorient , &imout_orig, NULL) ;
00126    } else {
00127       imout = mri_brainormalize( imin , iset->daxes->xxorient,
00128                                      iset->daxes->yyorient,
00129                                      iset->daxes->zzorient , NULL, NULL) ;
00130    }
00131    mri_free( imin ) ;
00132 
00133    if( imout == NULL ){
00134      fprintf(stderr,"**ERROR: normalization fails!?\n"); exit(1);
00135    }
00136    
00137    if (OrigSpace) {
00138       if( verb ) fprintf(stderr,"++3dSpatNorm: Output in Orignal space\n") ;
00139       mri_free( imout ) ;
00140       imout = imout_orig; imout->xo = originRAIfv.xyz[0]; imout->yo = originRAIfv.xyz[1]; imout->zo = originRAIfv.xyz[2]; 
00141       imout_orig = NULL;
00142    } else {
00143       if( verb ) fprintf(stderr,"++3dSpatNorm: Output in SpatNorm space\n") ;
00144    }
00145    
00146 #if 0
00147    if( AFNI_yesenv("WATERSHED") ){
00148      imin = mri_watershedize( imout , 0.10 ) ;
00149      if( imin != NULL ){ mri_free(imout); imout = imin; }
00150    }
00151 #endif
00152 
00153    /*--- create output dataset ---*/
00154    if( verb )
00155      fprintf(stderr,"++3dSpatNorm: Creating output dset\n") ;
00156 
00157    oset = EDIT_empty_copy( NULL ) ;
00158 
00159    tross_Copy_History( iset , oset ) ;
00160    tross_Make_History( "3dSpatNorm" , argc,argv , oset ) ;
00161 
00162    LOAD_IVEC3( nxyz   , imout->nx    , imout->ny    , imout->nz    ) ;
00163    LOAD_FVEC3( dxyz   , imout->dx    , imout->dy    , imout->dz    ) ;
00164    LOAD_FVEC3( orgxyz , imout->xo    , imout->yo    , imout->zo    ) ;
00165    LOAD_IVEC3( orixyz , ORI_R2L_TYPE , ORI_A2P_TYPE , ORI_I2S_TYPE ) ;
00166 
00167    if( verb )
00168      fprintf(stderr,"++3dSpatNorm: EDIT_dset_items\n") ;
00169    EDIT_dset_items( oset ,
00170                       ADN_prefix      , prefix ,
00171                       ADN_datum_all   , imout->kind ,
00172                       ADN_nxyz        , nxyz ,
00173                       ADN_xyzdel      , dxyz ,
00174                       ADN_xyzorg      , orgxyz ,
00175                       ADN_xyzorient   , orixyz ,
00176                       ADN_malloc_type , DATABLOCK_MEM_MALLOC ,
00177                       ADN_view_type   , VIEW_ORIGINAL_TYPE ,
00178                       ADN_type        , HEAD_ANAT_TYPE ,
00179                       ADN_func_type   , ANAT_BUCK_TYPE ,
00180                     ADN_none ) ;
00181 
00182    if( verb )
00183      fprintf(stderr,"++3dSpatNorm: EDIT_substitute_brick\n") ;
00184    EDIT_substitute_brick( oset , 0 , imout->kind , mri_data_pointer(imout) ) ;
00185 
00186    if (OrigSpace) {
00187       if( verb )
00188          fprintf(stderr,"++3dSpatNorm: Changing orientation from RAI\n") ;
00189       ooset = r_new_resam_dset ( oset, iset,    0,      0,      0,      NULL, MRI_NN, NULL);
00190       if (!ooset) {
00191          fprintf(stderr,"**ERROR: Failed to reslice!?\n"); exit(1);
00192       }
00193       DSET_delete(oset); oset = ooset; ooset = NULL;
00194    }
00195 
00196    DSET_write(oset) ;
00197    if( verb )
00198      fprintf(stderr,"++3dSpatNorm: wrote dataset %s\n",DSET_BRIKNAME(oset)) ;
00199    
00200    exit(0) ;
00201 }