Doxygen Source Code Documentation
Main Page Alphabetical List Data Structures File List Data Fields Globals Search
3dWinsor.c
Go to the documentation of this file.00001
00002
00003
00004
00005
00006
00007 #include "mrilib.h"
00008
00009
00010
00011 int main( int argc , char * argv[] )
00012 {
00013 float irad=1.5 ;
00014 int nrep=1 , cbot=-1,ctop=-1 ;
00015 char * prefix = "winsor" ;
00016 int keepzero = 0 , clipval = 0;
00017 int iarg ;
00018 THD_3dim_dataset *inset , *outset , *mset=NULL ;
00019 byte *mask=NULL ;
00020
00021 if( argc < 2 || strcmp(argv[1],"-help") == 0 ){
00022 printf("Usage: 3dWinsor [options] dataset\n"
00023 "Apply a 3D 'Winsorizing' filter to a short-valued dataset.\n"
00024 "\n"
00025 "Options:\n"
00026 " -irad rr = include all points within 'distance'\n"
00027 " rr in the operation, where distance\n"
00028 " is defined as sqrt(i*i+j*j+k*k), and\n"
00029 " (i,j,k) are voxel index offsets\n"
00030 " [default rr=1.5]\n"
00031 "\n"
00032 " -cbot bb = set bottom clip index to bb\n"
00033 " [default = 20%% of the number of points]\n"
00034 " -ctop tt = set top clip index to tt\n"
00035 " [default = 80%% of the number of points]\n"
00036 "\n"
00037 " -nrep nn = repeat filter nn times [default nn=1]\n"
00038 " if nn < 0, means to repeat filter until\n"
00039 " less than abs(n) voxels change\n"
00040 "\n"
00041 " -keepzero = don't filter voxels that are zero\n"
00042 " -clip xx = set voxels at or below 'xx' to zero\n"
00043 "\n"
00044 " -prefix pp = use 'pp' as the prefix for the output\n"
00045 " dataset [default pp='winsor']\n"
00046 "\n"
00047 " -mask mmm = use 'mmm' as a mask dataset - voxels NOT\n"
00048 " in the mask won't be filtered\n"
00049 ) ;
00050 exit(0) ;
00051 }
00052
00053 mainENTRY("3dWinsor main"); machdep(); AFNI_logger("3dWinsor",argc,argv);
00054 PRINT_VERSION("3dWinsor") ;
00055
00056
00057
00058 iarg = 1 ;
00059 while( iarg < argc && argv[iarg][0] == '-' ){
00060
00061 if( strcmp(argv[iarg],"-mask") == 0 ){
00062 mset = THD_open_dataset( argv[++iarg] ) ;
00063 DSET_load(mset) ;
00064 if( mset == NULL || !DSET_LOADED(mset) ){
00065 fprintf(stderr,"** Can't load mask dataset %s\n",argv[iarg]); exit(1);
00066 }
00067 iarg++ ; continue ;
00068 }
00069
00070 if( strcmp(argv[iarg],"-clip") == 0 ){
00071 clipval = strtol(argv[++iarg],NULL,10) ;
00072 iarg++ ; continue ;
00073 }
00074
00075 if( strcmp(argv[iarg],"-keepzero") == 0 ){
00076 keepzero = 1 ;
00077 iarg++ ; continue ;
00078 }
00079
00080 if( strcmp(argv[iarg],"-irad") == 0 ){
00081 irad = strtod( argv[++iarg] , NULL ) ;
00082 if( irad < 1.0 || irad > 4.0 ){
00083 fprintf(stderr,"*** Illegal value after -irad!\n"); exit(1);
00084 }
00085 iarg++ ; continue ;
00086 }
00087
00088 if( strcmp(argv[iarg],"-cbot") == 0 ){
00089 cbot = strtod( argv[++iarg] , NULL ) ;
00090 if( cbot < 1 ){
00091 fprintf(stderr,"*** Illegal value after -cbot!\n"); exit(1);
00092 }
00093 iarg++ ; continue ;
00094 }
00095
00096 if( strcmp(argv[iarg],"-ctop") == 0 ){
00097 ctop = strtod( argv[++iarg] , NULL ) ;
00098 if( ctop < 1 ){
00099 fprintf(stderr,"*** Illegal value after -ctop!\n"); exit(1);
00100 }
00101 iarg++ ; continue ;
00102 }
00103
00104 if( strcmp(argv[iarg],"-nrep") == 0 ){
00105 nrep = strtod( argv[++iarg] , NULL ) ;
00106 if( nrep == 0 ){
00107 fprintf(stderr,"*** Illegal value after -nrep!\n"); exit(1);
00108 }
00109 iarg++ ; continue ;
00110 }
00111
00112 if( strcmp(argv[iarg],"-prefix") == 0 ){
00113 prefix = argv[++iarg] ;
00114 if( !THD_filename_ok(prefix) ){
00115 fprintf(stderr,"*** Illegal value after -prefix!\n"); exit(1);
00116 }
00117 iarg++ ; continue ;
00118 }
00119
00120 fprintf(stderr,"*** Unknown option: %s\n",argv[iarg]); exit(1) ;
00121 }
00122
00123 if( iarg >= argc ){
00124 fprintf(stderr,"*** No dataset name on command line?\n"); exit(1);
00125 }
00126
00127
00128
00129 inset = THD_open_dataset( argv[iarg] ) ;
00130 if( inset == NULL ){
00131 fprintf(stderr,"*** Can't open dataset: %s\n",argv[iarg]); exit(1);
00132 }
00133
00134 if( DSET_BRICK_TYPE(inset,0) != MRI_short ){
00135 fprintf(stderr,"*** Dataset not stored as shorts!\n"); exit(1);
00136 }
00137
00138 if( mset != NULL ){
00139 if( DSET_NVOX(mset) != DSET_NVOX(inset) ){
00140 fprintf(stderr,"** Mask and input datasets don't match in size!\n"); exit(1);
00141 }
00142 mask = THD_makemask( mset , 0 , 1.0,0.0 ) ;
00143 if( mask == NULL ){
00144 fprintf(stderr,"** Can't make mask from mask dataset?!\n"); exit(1);
00145 }
00146 fprintf(stderr,"++ %d voxels in the mask\n",THD_countmask(DSET_NVOX(mset),mask));
00147 DSET_delete(mset) ;
00148 }
00149
00150 DSET_load(inset) ;
00151 if( !DSET_LOADED(inset) ){
00152 fprintf(stderr,"*** Can't read dataset into memory!\n"); exit(1);
00153 }
00154
00155 if( DSET_NVALS(inset) > 1 ){
00156 fprintf(stderr,"+++ WARNING: only processing sub-brick #0\n") ;
00157 }
00158
00159
00160
00161 outset = WINsorize( inset, nrep, cbot, ctop, irad, prefix, keepzero,clipval , mask ) ;
00162
00163 if( outset == NULL ){
00164 fprintf(stderr,"*** Can't compute Winsor filter!\n"); exit(1);
00165 }
00166
00167 tross_Copy_History( inset , outset ) ;
00168 tross_Make_History( "3dWinsor" , argc,argv , outset ) ;
00169 DSET_write(outset) ;
00170 fprintf(stderr,"++ output dataset: %s\n",DSET_BRIKNAME(outset)) ;
00171 exit(0) ;
00172 }