Doxygen Source Code Documentation
3dAutomask.c File Reference
#include "mrilib.h"
Go to the source code of this file.
Functions | |
int | main (int argc, char *argv[]) |
Function Documentation
|
compute the overall minimum and maximum voxel values for a dataset Definition at line 5 of file 3dAutomask.c. References ADN_datum_all, ADN_func_type, ADN_none, ADN_ntt, ADN_nvals, ADN_prefix, ADN_type, AFNI_logger(), argc, COX_cpu_time(), THD_3dim_dataset::daxes, DSET_BRICK_TYPE, DSET_DX, DSET_DY, DSET_DZ, DSET_load, DSET_LOADED, DSET_NVOX, DSET_NX, DSET_NY, DSET_NZ, DSET_unload, DSET_write, EDIT_dset_items(), EDIT_empty_copy(), EDIT_substitute_brick(), ERROR_exit(), FUNC_FIM_TYPE, HEAD_FUNC_TYPE, INFO_message(), ISVALID_DSET, machdep(), mainENTRY, MAX, nz, PRINT_VERSION, strtod(), THD_autobbox(), THD_automask(), THD_automask_extclip(), THD_automask_verbose(), THD_countmask(), THD_filename_ok(), THD_mask_clust(), THD_mask_dilate(), THD_mask_fillin_completely(), THD_open_dataset(), tross_Copy_History(), tross_Make_History(), WROTE_DSET, THD_dataxes::xxorient, THD_dataxes::yyorient, and THD_dataxes::zzorient.
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 } |