00001 #include "mrilib.h"
00002 
00003 THD_fvec3 THD_cmass( THD_3dim_dataset *xset , int iv , byte *mmm ) ;
00004 
00005 int main( int argc , char * argv[] )
00006 {
00007    int narg=1, do_automask=0 , iv , nxyz , do_set=0 ;
00008    THD_3dim_dataset *xset ;
00009    byte *mmm=NULL ; int nmask=0 , nvox_mask ;
00010    THD_fvec3 cmv , setv ;
00011 
00012    /*-- read command line arguments --*/
00013 
00014    if( argc < 2 || strncmp(argv[1],"-help",5) == 0 ){
00015       printf("Usage: 3dCM [options] dset\n"
00016              "Output = center of mass of dataset, to stdout.\n"
00017              "  -mask mset   Means to use the dataset 'mset' as a mask:\n"
00018              "                 Only voxels with nonzero values in 'mset'\n"
00019              "                 will be averaged from 'dataset'.  Note\n"
00020              "                 that the mask dataset and the input dataset\n"
00021              "                 must have the same number of voxels.\n"
00022              "  -automask    Generate the mask automatically.\n"
00023              "  -set x y z   After computing the CM of the dataset, set the\n"
00024              "                 origin fields in the header so that the CM\n"
00025              "                 will be at (x,y,z) in DICOM coords.\n"
00026             ) ;
00027       exit(0) ;
00028    }
00029 
00030    narg = 1 ;
00031    while( narg < argc && argv[narg][0] == '-' ){
00032 
00033       if( strcmp(argv[narg],"-set") == 0 ){
00034         float xset,yset,zset ;
00035         if( narg+3 >= argc ){
00036           fprintf(stderr,"*** -set need 3 args following!\n") ; exit(1) ;
00037         }
00038         xset = strtod( argv[++narg] , NULL ) ;
00039         yset = strtod( argv[++narg] , NULL ) ;
00040         zset = strtod( argv[++narg] , NULL ) ;
00041         LOAD_FVEC3(setv,xset,yset,zset) ; do_set = 1 ;
00042         narg++ ; continue ;
00043       }
00044 
00045       if( strncmp(argv[narg],"-mask",5) == 0 ){
00046         THD_3dim_dataset *mask_dset ;
00047         if( mmm != NULL ){
00048           fprintf(stderr,"*** Cannot have two -mask options!\n") ; exit(1) ;
00049         }
00050         if( do_automask ){
00051           fprintf(stderr,"*** Can't have -mask and -automask!\n") ; exit(1) ;
00052         }
00053         if( narg+1 >= argc ){
00054           fprintf(stderr,"*** -mask option requires a following argument!\n");
00055           exit(1) ;
00056         }
00057         mask_dset = THD_open_dataset( argv[++narg] ) ;
00058         if( mask_dset == NULL ){
00059           fprintf(stderr,"*** Cannot open mask dataset!\n") ; exit(1) ;
00060         }
00061         if( DSET_BRICK_TYPE(mask_dset,0) == MRI_complex ){
00062           fprintf(stderr,"*** Cannot deal with complex-valued mask dataset!\n");
00063           exit(1) ;
00064         }
00065         mmm = THD_makemask( mask_dset , 0 , 1.0,0.0 ) ;
00066         nvox_mask = DSET_NVOX(mask_dset) ;
00067         nmask = THD_countmask( nvox_mask , mmm ) ;
00068         if( mmm == NULL || nmask <= 0 ){
00069           fprintf(stderr,"*** Can't make mask from dataset %s\n",argv[narg-1]);
00070           exit(1) ;
00071         }
00072         DSET_delete( mask_dset ) ;
00073         narg++ ; continue ;
00074       }
00075 
00076       if( strcmp(argv[narg],"-automask") == 0 ){
00077         if( mmm != NULL ){
00078           fprintf(stderr,"*** Can't have -mask and -automask!\n") ; exit(1) ;
00079         }
00080         do_automask = 1 ; narg++ ; continue ;
00081       }
00082 
00083       fprintf(stderr,"*** Unknown option: %s\n",argv[narg]) ; exit(1) ;
00084    }
00085 
00086    /* should have at least 1 more argument */
00087 
00088    if( argc <= narg ){
00089      fprintf(stderr,"*** No input dataset!?\n") ; exit(1) ;
00090    }
00091 
00092    for( ; narg < argc ; narg++ ){
00093      xset = THD_open_dataset( argv[narg] ) ;
00094      if( xset == NULL ){
00095        fprintf(stderr,"+++ Can't open dataset %s -- skipping\n",argv[narg]);
00096        continue ;
00097      }
00098      DSET_load(xset) ;
00099      if( !DSET_LOADED(xset) ){
00100        fprintf(stderr,"+++ Can't load dataset %s -- skipping\n",argv[narg]);
00101        DSET_delete(xset) ; continue ;
00102      }
00103      if( do_automask ){
00104        if( mmm != NULL ){ free(mmm); mmm = NULL; }
00105        mmm = THD_automask(xset) ;
00106        nvox_mask = DSET_NVOX(xset) ;
00107        nmask = THD_countmask( nvox_mask , mmm ) ;
00108        if( mmm == NULL || nmask <= 0 ){
00109          fprintf(stderr,"+++ Can't make automask from dataset %s -- skipping\n",argv[narg]) ;
00110          DSET_delete(xset) ; continue ;
00111        }
00112      }
00113      nxyz = DSET_NVOX(xset) ;
00114      if( mmm != NULL && nxyz != nvox_mask ){
00115        fprintf(stderr,"+++ Mask/Dataset grid size mismatch at %s\n -- skipping\n",argv[narg]) ;
00116        DSET_delete(xset) ; continue ;
00117      }
00118 
00119      cmv = THD_cmass( xset , 0 , mmm ) ;
00120      printf("%g  %g  %g\n",cmv.xyz[0],cmv.xyz[1],cmv.xyz[2]) ;
00121      DSET_unload(xset) ;
00122 
00123      if( do_set ){
00124        THD_fvec3 dv , ov ;
00125        if( !DSET_IS_BRIK(xset) || DSET_IS_MASTERED(xset) ){
00126          fprintf(stderr,"+++ Can't modify CM of dataset %s\n",argv[narg]) ;
00127        } else {
00128          LOAD_FVEC3(ov,DSET_XORG(xset),DSET_YORG(xset),DSET_ZORG(xset)) ;
00129          ov = THD_3dmm_to_dicomm( xset , ov ) ;
00130          dv = SUB_FVEC3(setv,cmv) ;
00131          ov = ADD_FVEC3(dv,ov) ;
00132          ov = THD_dicomm_to_3dmm( xset , ov ) ;
00133          xset->daxes->xxorg = ov.xyz[0] ;
00134          xset->daxes->yyorg = ov.xyz[1] ;
00135          xset->daxes->zzorg = ov.xyz[2] ;
00136          fprintf(stderr,"+++ Rewriting header %s\n",DSET_HEADNAME(xset)) ;
00137          DSET_write_header( xset ) ;
00138        }
00139      }
00140      DSET_delete(xset) ;
00141   }
00142 }
00143 
00144 /*----------------------------------------------------------------------*/
00145 /*! Get the center of mass of this volume, in DICOM coords.
00146 ------------------------------------------------------------------------*/
00147 
00148 THD_fvec3 THD_cmass( THD_3dim_dataset *xset , int iv , byte *mmm )
00149 {
00150    THD_fvec3 cmv ;
00151    MRI_IMAGE *im ;
00152    float *far , icm,jcm,kcm ;
00153    int ii , nvox ;
00154 
00155    LOAD_FVEC3(cmv,0,0,0) ;
00156 
00157    nvox = DSET_NVOX(xset) ;
00158    im   = mri_to_float( DSET_BRICK(xset,iv) ) ;
00159                              if( im  == NULL ) return cmv ;
00160    far = MRI_FLOAT_PTR(im) ; if( far == NULL ) return cmv ;
00161 
00162    if( mmm != NULL ){
00163      for( ii=0 ; ii < nvox ; ii++ )
00164        if( mmm[ii] ) far[ii] = 0.0 ;
00165    }
00166 
00167    mri_get_cmass_3D( im , &icm,&jcm,&kcm ) ; mri_free(im) ;
00168    LOAD_FVEC3(cmv,icm,jcm,kcm) ;
00169    cmv = THD_3dfind_to_3dmm( xset , cmv ) ;
00170    cmv = THD_3dmm_to_dicomm( xset , cmv ) ;
00171    return cmv ;
00172 }