#include "mrilib.h"

int main( int argc , char * argv[] )
{
   float mfrac=0.50f , val ;
   MRI_IMAGE *medim ;
   THD_3dim_dataset *dset ;
   char *gprefix=NULL ;
   int iarg=1 ;

   if( argc < 2 || strcmp(argv[1],"-help") == 0 ){
     printf(
       "Usage: 3dClipLevel [options] dataset\n"
       "Estimates the value at which to clip the anatomical dataset so\n"
       "  that background regions are set to zero.\n"
       "\n"
       "The program's output is a single number sent to stdout.  This\n"
       "  value can be 'captured' to a shell variable using the backward\n"
       "  single quote operator; a trivial csh/tcsh example is\n"
       "\n"
       "    set ccc = `3dClipLevel -mfrac 0.333 Elvis+orig`\n"
       "    3dcalc -a Elvis+orig -expr \"step(a-$ccc)\" -prefix Presley\n"
       "\n"
       "Algorithm:\n"
       "  (a) Set some initial clip value using wizardry (AKA 'variance').\n"
       "  (b) Find the median of all positive values >= clip value.\n"
       "  (c) Set the clip value to 0.50 of this median.\n"
       "  (d) Loop back to (b) until the clip value doesn't change.\n"
       "This method was made up out of nothing, based on histogram gazing.\n"
       "\n"
       "Options:\n"
       "--------\n"
       "  -mfrac ff = Use the number ff instead of 0.50 in the algorithm.\n"
       "\n"
       "  -grad ppp = In addition to using the 'one size fits all routine',\n"
       "              also compute a 'gradual' clip level as a function\n"
       "              of voxel position, and output that to a dataset with\n"
       "              prefix 'ppp'.\n"
       "             [This is the same 'gradual' clip level that is now the\n"
       "              default in 3dAutomask - as of 24 Oct 2006.\n"
       "              You can use this option to see how 3dAutomask clips\n"
       "              the dataset as its first step.  The algorithm above is\n"
       "              is used in each octant of the dataset, and then these\n"
       "              8 values are interpolated to cover the whole volume.]\n"
       "Notes:\n"
       "------\n"
       "* Use at your own risk!  You might want to use the AFNI Histogram\n"
       "    plugin to see if the results are reasonable.  This program is\n"
       "    likely to produce bad results on images gathered with local\n"
       "    RF coils, or with pulse sequences with unusual contrasts.\n"
       "\n"
       "* For brain images, most brain voxels seem to be in the range from\n"
       "    the clip level (mfrac=0.5) to about 3-3.5 times the clip level.\n"
       "    - In T1-weighted images, voxels above that level are usually\n"
       "      blood vessels (e.g., inflow artifact brightens them).\n"
       "\n"
       "* If the input dataset has more than 1 sub-brick, the data is\n"
       "    analyzed on the median volume -- at each voxel, the median\n"
       "    of all sub-bricks at that voxel is computed, and then this\n"
       "    median volume is used in the histogram algorithm.\n"
       "\n"
       "* If the input dataset is short- or byte-valued, the output will\n"
       "    be an integer; otherwise, the output is a float value.\n"
       "------\n"
       "Author: Emperor Zhark -- Sadistic Galactic Domination since 1994!\n\n"
     ) ;
     exit(0) ;
   }

   mainENTRY("3dClipLevel main"); machdep(); AFNI_logger("3dCliplevel",argc,argv);

   /*-- options --*/

   while( iarg < argc && argv[iarg][0] == '-' ){

     if( strcmp(argv[iarg],"-mfrac") == 0 || strcmp(argv[iarg],"-clfrac") == 0 ){
       mfrac = strtod( argv[++iarg] , NULL ) ;
       if( mfrac <= 0.0f ) ERROR_exit("Illegal -mfrac '%s'",argv[iarg]) ;
       while( mfrac >= 1.0f ) mfrac *= 0.01f ;
       iarg++ ; continue ;
     }

     if( strcmp(argv[iarg],"-grad") == 0 ){
       gprefix = argv[++iarg] ;
       if( !THD_filename_ok(gprefix) ) ERROR_exit("Illegal -grad '%s'",argv[iarg]) ;
       iarg++ ; continue ;
     }

     ERROR_exit("Unknown option: %s\n",argv[iarg]) ;
   }

   /*-- read input dataset --*/

   dset = THD_open_dataset(argv[iarg]) ;
   CHECK_OPEN_ERROR(dset,argv[iarg]) ;

   /*-- get median at each voxel --*/

   medim = THD_median_brick( dset ) ;
   if( medim == NULL )       ERROR_exit("Can't load dataset '%s'",argv[iarg]);

   DSET_unload(dset) ;  /* no longer needed */

   val = THD_cliplevel( medim , mfrac ) ;  /* floating point clip level */

   if( !THD_need_brick_factor(dset) &&     /* convert to integer? */
       DSET_datum_constant(dset)    &&
      (DSET_BRICK_TYPE(dset,0)==MRI_short || DSET_BRICK_TYPE(dset,0)==MRI_byte) )
     val = (float)rint((double)val) ;

   printf("%g\n",val) ;        /***** write the output value! *****/

   /*--- 25 Oct 2006: create the gradual clip dataset? ---*/

   if( gprefix != NULL && val > 0.0f ){
     THD_3dim_dataset *gset ;
     MRI_IMAGE *gim ;
     gim = THD_cliplevel_gradual( medim , mfrac ) ;
     if( gim == NULL ) ERROR_exit("Can't compute gradual clip?!") ;
     gset = EDIT_empty_copy(dset) ;
     EDIT_dset_items( gset ,
                        ADN_nvals , 1 ,
                        ADN_ntt   , 0 ,
                        ADN_prefix, gprefix ,
                      ADN_none ) ;
     EDIT_substitute_brick( gset , 0 , MRI_float , MRI_FLOAT_PTR(gim) ) ;
     tross_Copy_History( dset , gset ) ;
     tross_Make_History( "3dClipLevel" , argc,argv , gset ) ;
     DSET_write(gset) ; DSET_delete(gset) ;
   }

   exit(0) ;
}