00001 #include "mrilib.h"
00002 
00003 #undef ALLOW_FILLIN  /* 28 May 2002 */
00004 
00005 int main( int argc , char * argv[] )
00006 {
00007    THD_3dim_dataset *dset , *mset ;
00008    char *prefix = "automask" ;
00009    byte *mask ;
00010    int iarg=1 , fillin=0 , nmask,nfill , dilate=0 , dd ;
00011    float SIhh=0.0 ;        /* 06 Mar 2003 */
00012    int   SIax=0 , SIbot,SItop ;
00013    int   verb=1 ;
00014 
00015    if( argc < 2 || strcmp(argv[1],"-help") == 0 ){
00016       printf("Usage: 3dAutomask [options] dataset\n"
00017              "Input dataset is EPI 3D+time.\n"
00018              "Output dataset is a brain-only mask dataset.\n"
00019              "Method:\n"
00020              " + Uses 3dClipLevel algorithm to find clipping level.\n"
00021              " + Keeps only the largest connected component of the\n"
00022              "   supra-threshold voxels, after an erosion/dilation step.\n"
00023              " + Writes result as a 'fim' type of functional dataset.\n"
00024              "Options:\n"
00025              "  -prefix ppp = Write mask into dataset with prefix 'ppp'.\n"
00026              "                 [default='automask']\n"
00027              "  -q          = Don't write progress messages (i.e., be quiet).\n"
00028              "  -eclip      = After creating the mask, remove exterior\n"
00029              "                 voxels below the clip threshold.\n"
00030              "  -dilate nd  = Dilate the mask outwards 'nd' times.\n"
00031 #ifdef ALLOW_FILLIN
00032              "  -fillin nnn = Fill in holes inside the mask of width up\n"
00033              "                 to 'nnn' voxels. [default=0=no fillin]\n"
00034 #endif
00035              "  -SI hh      = After creating the mask, find the most superior\n"
00036              "                 voxel, then zero out everything more than 'hh'\n"
00037              "                 millimeters inferior to that.  hh=130 seems to\n"
00038              "                 be decent (for human brains).\n"
00039             ) ;
00040       exit(0) ;
00041    }
00042 
00043    mainENTRY("3dAutomask main"); machdep(); AFNI_logger("3dAutomask",argc,argv);
00044    PRINT_VERSION("3dAutomask") ;
00045 
00046    /*-- options --*/
00047 
00048    while( iarg < argc && argv[iarg][0] == '-' ){
00049 
00050       if( strcmp(argv[iarg],"-SI") == 0 ){        /* 06 Mar 2003 */
00051         SIhh = strtod( argv[++iarg] , NULL ) ;
00052         if( SIhh <= 0.0 )
00053           ERROR_exit("-SI value %f is illegal!\n",SIhh) ;
00054         iarg++ ; continue ;
00055       }
00056 
00057       if( strcmp(argv[iarg],"-eclip") == 0 ){     /* 28 Oct 2003 */
00058         THD_automask_extclip(1) ;
00059         iarg++ ; continue ;
00060       }
00061 
00062       if( strcmp(argv[iarg],"-q") == 0 ){     /* 28 Oct 2003 */
00063         verb = 0 ; iarg++ ; continue ;
00064       }
00065 
00066       if( strcmp(argv[iarg],"-prefix") == 0 ){
00067          prefix = argv[++iarg] ;
00068          if( !THD_filename_ok(prefix) )
00069            ERROR_exit("-prefix %s is illegal!\n",prefix) ;
00070          iarg++ ; continue ;
00071       }
00072 
00073 #ifdef ALLOW_FILLIN
00074       if( strcmp(argv[iarg],"-fillin") == 0 ){
00075          fillin = strtol( argv[++iarg] , NULL , 10 ) ;
00076          if( fillin <  0 )
00077            ERROR_exit("-fillin %s is illegal!\n",argv[iarg]) ;
00078          else if( fillin > 0 )
00079            fillin = (fillin+2) / 2 ;
00080          iarg++ ; continue ;
00081       }
00082 #endif
00083 
00084       if( strcmp(argv[iarg],"-dilate") == 0 ){
00085          dilate = strtol( argv[++iarg] , NULL , 10 ) ;
00086          if( dilate < 0 )
00087            ERROR_exit("-dilate %s is illegal!\n",argv[iarg]);
00088          iarg++ ; continue ;
00089       }
00090 
00091       ERROR_exit("ILLEGAL option: %s\n",argv[iarg]) ;
00092    }
00093 
00094    /*-- read data --*/
00095 
00096    dset = THD_open_dataset(argv[iarg]) ;
00097    if( !ISVALID_DSET(dset) ) ERROR_exit("Can't open dataset %s\n",argv[iarg]);
00098    if( DSET_BRICK_TYPE(dset,0) != MRI_short &&
00099        DSET_BRICK_TYPE(dset,0) != MRI_byte  &&
00100        DSET_BRICK_TYPE(dset,0) != MRI_float   ){
00101       ERROR_exit("Illegal dataset datum type\n") ;
00102    }
00103    if( verb ) INFO_message("Loading dataset %s\n",argv[iarg]) ;
00104    DSET_load(dset) ;
00105    if( !DSET_LOADED(dset) ) ERROR_exit("Can't load dataset %s\n",argv[iarg]);
00106 
00107    /*** do all the real work now ***/
00108 
00109    if( verb ) INFO_message("Forming automask\n") ;
00110    if( verb ) THD_automask_verbose(1) ;
00111    mask = THD_automask( dset ) ;
00112    if( mask == NULL )
00113      ERROR_exit("Mask creation fails for unknown reasons!\n");
00114 
00115    /* 30 Aug 2002 (modified 05 Mar 2003 to do fillin, etc, after dilation) */
00116 
00117    if( dilate > 0 ){
00118      int ii,nx,ny,nz , nmm ;
00119      if( verb ) INFO_message("Dilating automask\n") ;
00120      nx = DSET_NX(dset) ; ny = DSET_NY(dset) ; nz = DSET_NZ(dset) ;
00121      nmm = 1 ;
00122      ii  = rint(0.032*nx) ; nmm = MAX(nmm,ii) ;
00123      ii  = rint(0.032*ny) ; nmm = MAX(nmm,ii) ;
00124      ii  = rint(0.032*nz) ; nmm = MAX(nmm,ii) ;
00125      for( dd=0 ; dd < dilate ; dd++ ){
00126        THD_mask_dilate           ( nx,ny,nz , mask, 3   ) ;
00127        THD_mask_fillin_completely( nx,ny,nz , mask, nmm ) ;
00128      }
00129      nmm = nx*ny*nz ;
00130      for( ii=0 ; ii < nmm ; ii++ ) mask[ii] = !mask[ii] ;
00131      THD_mask_clust( nx,ny,nz, mask ) ;
00132      for( ii=0 ; ii < nmm ; ii++ ) mask[ii] = !mask[ii] ;
00133    }
00134 
00135    /* 18 Apr 2002: print voxel count */
00136 
00137    nmask = THD_countmask( DSET_NVOX(dset) , mask ) ;
00138    if( verb ) INFO_message("%d voxels in the mask [out of %d: %.2f%%]\n",
00139                  nmask,DSET_NVOX(dset), (100.0*nmask)/DSET_NVOX(dset) ) ;
00140    if( nmask == 0 )
00141       ERROR_exit("No voxels? Quitting without saving mask\n");
00142 
00143    /* 18 Apr 2002: maybe fill in voxels */
00144 
00145 #ifdef ALLOW_FILLIN
00146    if( fillin > 0 ){
00147      nfill = THD_mask_fillin_completely(
00148                  DSET_NX(dset),DSET_NY(dset),DSET_NZ(dset), mask, fillin ) ;
00149      if( verb ) INFO_message("%d voxels filled in; %d voxels total\n",
00150                              nfill,nfill+nmask ) ;
00151    }
00152 #endif
00153 
00154    /** 04 Jun 2002: print cut plane report **/
00155 
00156    { int nx=DSET_NX(dset), ny=DSET_NY(dset), nz=DSET_NZ(dset), nxy=nx*ny ;
00157      int ii,jj,kk ;
00158 
00159 #if 0
00160      { int xm=-1,xp=-1,ym=-1,yp=-1,zm=-1,zp=-1 ;
00161        THD_autobbox( dset , &xm,&xp , &ym,&yp , &zm,&zp ) ;
00162        INFO_message("Auto bbox: x=%d..%d  y=%d..%d  z=%d..%d\n",
00163                      xm,xp,ym,yp,zm,zp ) ;
00164      }
00165 #endif
00166 
00167      for( ii=0 ; ii < nx ; ii++ )
00168        for( kk=0 ; kk < nz ; kk++ )
00169          for( jj=0 ; jj < ny ; jj++ )
00170            if( mask[ii+jj*nx+kk*nxy] ) goto CP5 ;
00171      CP5: if( verb )
00172            INFO_message("first %3d x-planes are zero [from %c]\n",
00173                         ii,ORIENT_tinystr[dset->daxes->xxorient][0]) ;
00174      if( SIhh > 0.0 && ORIENT_tinystr[dset->daxes->xxorient][0] == 'S' ){
00175        SIax = 1 ; SIbot = ii + (int)(SIhh/fabs(DSET_DX(dset))+0.5) ; SItop = nx-1 ;
00176      }
00177 
00178      for( ii=nx-1 ; ii >= 0 ; ii-- )
00179        for( kk=0 ; kk < nz ; kk++ )
00180          for( jj=0 ; jj < ny ; jj++ )
00181            if( mask[ii+jj*nx+kk*nxy] ) goto CP6 ;
00182      CP6: if( verb )
00183            INFO_message("last  %3d x-planes are zero [from %c]\n",
00184                         nx-1-ii,ORIENT_tinystr[dset->daxes->xxorient][1]) ;
00185      if( SIhh > 0.0 && ORIENT_tinystr[dset->daxes->xxorient][1] == 'S' ){
00186        SIax = 1 ; SIbot = 0 ; SItop = ii - (int)(SIhh/fabs(DSET_DX(dset))+0.5) ;
00187      }
00188 
00189      for( jj=0 ; jj < ny ; jj++ )
00190        for( kk=0 ; kk < nz ; kk++ )
00191          for( ii=0 ; ii < nx ; ii++ )
00192            if( mask[ii+jj*nx+kk*nxy] ) goto CP3 ;
00193      CP3: if( verb )
00194            INFO_message("first %3d y-planes are zero [from %c]\n",
00195                         jj,ORIENT_tinystr[dset->daxes->yyorient][0]) ;
00196      if( SIhh > 0.0 && ORIENT_tinystr[dset->daxes->yyorient][0] == 'S' ){
00197        SIax = 2 ; SIbot = jj + (int)(SIhh/fabs(DSET_DY(dset))+0.5) ; SItop = ny-1 ;
00198      }
00199 
00200      for( jj=ny-1 ; jj >= 0 ; jj-- )
00201        for( kk=0 ; kk < nz ; kk++ )
00202          for( ii=0 ; ii < nx ; ii++ )
00203            if( mask[ii+jj*nx+kk*nxy] ) goto CP4 ;
00204      CP4: if( verb )
00205            INFO_message("last  %3d y-planes are zero [from %c]\n",
00206                         ny-1-jj,ORIENT_tinystr[dset->daxes->yyorient][1]) ;
00207      if( SIhh > 0.0 && ORIENT_tinystr[dset->daxes->yyorient][1] == 'S' ){
00208        SIax = 2 ; SIbot = 0 ; SItop = jj - (int)(SIhh/fabs(DSET_DY(dset))+0.5) ;
00209      }
00210 
00211      for( kk=0 ; kk < nz ; kk++ )
00212        for( jj=0 ; jj < ny ; jj++ )
00213          for( ii=0 ; ii < nx ; ii++ )
00214            if( mask[ii+jj*nx+kk*nxy] ) goto CP1 ;
00215      CP1: if( verb )
00216            INFO_message("first %3d z-planes are zero [from %c]\n",
00217                         kk,ORIENT_tinystr[dset->daxes->zzorient][0]) ;
00218      if( SIhh > 0.0 && ORIENT_tinystr[dset->daxes->zzorient][0] == 'S' ){
00219        SIax = 3 ; SIbot = kk + (int)(SIhh/fabs(DSET_DZ(dset))+0.5) ; SItop = nz-1 ;
00220      }
00221 
00222      for( kk=nz-1 ; kk >= 0 ; kk-- )
00223        for( jj=0 ; jj < ny ; jj++ )
00224          for( ii=0 ; ii < nx ; ii++ )
00225            if( mask[ii+jj*nx+kk*nxy] ) goto CP2 ;
00226      CP2: if( verb )
00227            INFO_message("last  %3d z-planes are zero [from %c]\n",
00228                         nz-1-kk,ORIENT_tinystr[dset->daxes->zzorient][1]) ;
00229      if( SIhh > 0.0 && ORIENT_tinystr[dset->daxes->zzorient][1] == 'S' ){
00230        SIax = 3 ; SIbot = 0 ; SItop = kk - (int)(SIhh/fabs(DSET_DZ(dset))+0.5) ;
00231      }
00232 
00233      /* 06 Mar 2003: cut off stuff below SIhh mm from most Superior point */
00234 
00235      if( SIax > 0 && SIbot <= SItop ){
00236        char *cax="xyz" ;
00237        if( verb )
00238          INFO_message("SI clipping mask along axis %c from %d..%d\n" ,
00239                       cax[SIax-1] , SIbot,SItop ) ;
00240        switch( SIax ){
00241          case 1:
00242            for( ii=SIbot ; ii <= SItop ; ii++ )
00243              for( kk=0 ; kk < nz ; kk++ )
00244                for( jj=0 ; jj < ny ; jj++ ) mask[ii+jj*nx+kk*nxy] = 0 ;
00245          break ;
00246          case 2:
00247            for( jj=SIbot ; jj <= SItop ; jj++ )
00248              for( kk=0 ; kk < nz ; kk++ )
00249                for( ii=0 ; ii < nx ; ii++ ) mask[ii+jj*nx+kk*nxy] = 0 ;
00250          break ;
00251          case 3:
00252            for( kk=SIbot ; kk <= SItop ; kk++ )
00253              for( jj=0 ; jj < ny ; jj++ )
00254                for( ii=0 ; ii < nx ; ii++ ) mask[ii+jj*nx+kk*nxy] = 0 ;
00255          break ;
00256        }
00257        nmask = THD_countmask( DSET_NVOX(dset) , mask ) ;
00258        if( verb )
00259          INFO_message("%d voxels left [out of %d]\n",nmask,DSET_NVOX(dset)) ;
00260      }
00261    }
00262 
00263    DSET_unload( dset ) ;  /* don't need data any more */
00264 
00265    /* create output dataset */
00266 
00267    mset = EDIT_empty_copy( dset ) ;
00268    EDIT_dset_items( mset ,
00269                       ADN_prefix     , prefix   ,
00270                       ADN_datum_all  , MRI_byte ,
00271                       ADN_nvals      , 1        ,
00272                       ADN_ntt        , 0        ,
00273                       ADN_type       , HEAD_FUNC_TYPE ,
00274                       ADN_func_type  , FUNC_FIM_TYPE ,
00275                     ADN_none ) ;
00276    EDIT_substitute_brick( mset , 0 , MRI_byte , mask ) ;
00277 
00278    /* 16 Apr 2002: make history */
00279 
00280    tross_Copy_History( dset , mset ) ;
00281    tross_Make_History( "3dAutomask", argc,argv, mset ) ;
00282 
00283    DSET_write( mset ) ;
00284    if( verb ) WROTE_DSET(mset) ;
00285    if( verb ) INFO_message("CPU time = %f sec\n",COX_cpu_time()) ;
00286    exit(0) ;
00287 }