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    /*-- scan command line --*/
00057 
00058    iarg = 1 ;
00059    while( iarg < argc && argv[iarg][0] == '-' ){
00060 
00061       if( strcmp(argv[iarg],"-mask") == 0 ){         /* 06 Mar 2003 */
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    /*-- read input --*/
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    /*-- compute output --*/
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 }