00004 {
00005    float mfrac=0.50 , fac ;
00006    double dsum ;
00007    int nvox , iarg=1 , *hist , ii,npos=0 , ncut,nmed,kk,ib , qq,nold ;
00008    THD_3dim_dataset * dset ;
00009    short * sar ;
00010    byte * bar ;
00011    int nhist , nneg=0 , verb=0 , nhalf , iv,nvals ;
00012 
00013    if( argc < 2 || strcmp(argv[1],"-help") == 0 ){
00014       printf("Usage: 3dClipLevel [options] dataset\n"
00015              "Estimates the value at which to clip the anatomical dataset so\n"
00016              "that background regions are set to zero.\n"
00017              "Method:\n"
00018              "  Find the median of all positive values >= clip value.\n"
00019              "  Set the clip value to 0.50 of this median.\n"
00020              "  Repeat until the clip value doesn't change.\n"
00021              "Options:\n"
00022              "  -mfrac ff = Use the number ff instead of 0.50 in the algorithm.\n"
00023              "  -verb     = The clip value is always printed to stdout.  If\n"
00024              "                this option is used to select verbose output,\n"
00025              "                progress reports are printed to stderr as well.\n"
00026              "\n"
00027              "N.B.: This program only works with byte- and short-valued\n"
00028              "        datasets, and prints a warning message if any input\n"
00029              "        voxels are negative.  If the dataset has more than one\n"
00030              "        sub-brick, all sub-bricks are used to build the histogram.\n"
00031              "N.B.: Use at your own risk!  You might want to use the AFNI Histogram\n"
00032              "        plugin to see if the results are reasonable.  This program is\n"
00033              "        likely to produce bad results on images gathered with local\n"
00034              "        RF coils, or with pulse sequences with unusual contrasts.\n"
00035              "\n"
00036              "A csh command line for the truly adventurous:\n"
00037              "  afni -dset \"v1:time+orig<`3dClipLevel 'v1:time+orig[4]'` .. 10000>\"\n"
00038              "(the dataset is from the 'sample96.tgz' data samples).  Can you\n"
00039              "figure out what this does?\n"
00040              "(Hint: each type of quote \"'` means something different to csh.)\n"
00041             ) ;
00042       exit(0) ;
00043    }
00044 
00045    mainENTRY("3dClipLevel main"); machdep(); AFNI_logger("3dCliplevel",argc,argv);
00046 
00047    /*-- options --*/
00048 
00049    while( iarg < argc && argv[iarg][0] == '-' ){
00050 
00051       if( strcmp(argv[iarg],"-mfrac") == 0 ){
00052          mfrac = strtod( argv[++iarg] , NULL ) ;
00053          if( mfrac <= 0.0 ){ fprintf(stderr,"** ILLEGAL -mfrac\n");exit(1);}
00054          if( mfrac >= 1.0 ) mfrac *= 0.01 ;
00055          iarg++ ; continue ;
00056       }
00057 
00058       if( strncmp(argv[iarg],"-verb",5) == 0 ){
00059          verb++ ; iarg++ ; continue ;
00060       }
00061 
00062       fprintf(stderr,"** ILLEGAL option: %s\n",argv[iarg]) ; exit(1) ;
00063    }
00064 
00065    /*-- read data --*/
00066 
00067    dset = THD_open_dataset(argv[iarg]) ;
00068    if( !ISVALID_DSET(dset) ){ fprintf(stderr,"** CAN'T open dataset\n");exit(1); }
00069    if( DSET_BRICK_TYPE(dset,0) != MRI_short && DSET_BRICK_TYPE(dset,0) != MRI_byte ){
00070       fprintf(stderr,"** ILLEGAL dataset type\n");exit(1);
00071    }
00072    DSET_load(dset) ;
00073    if( !DSET_LOADED(dset) ){ fprintf(stderr,"** CAN'T load dataset\n");exit(1); }
00074 
00075    nvals = DSET_NVALS(dset) ;
00076 
00077    /*-- allocate histogram --*/
00078 
00079    switch( DSET_BRICK_TYPE(dset,0) ){
00080       case MRI_short: nhist = 32767 ; break ;
00081       case MRI_byte : nhist =   255 ; break ;
00082    }
00083 
00084    hist = (int *) calloc(sizeof(int),nhist+1) ;  /* 05 Nov 2001: +1 */
00085    nvox = DSET_NVOX(dset) ;
00086 
00087    /*-- make histogram --*/
00088 
00089    dsum = 0.0 ;
00090    for( iv=0 ; iv < nvals ; iv++ ){
00091       switch( DSET_BRICK_TYPE(dset,0) ){
00092          case MRI_short:
00093             sar =  DSET_ARRAY(dset,iv) ;
00094             for( ii=0 ; ii < nvox ; ii++ ){
00095                if( sar[ii] > 0 && sar[ii] <= nhist ){
00096                   hist[sar[ii]]++ ; dsum += (double)(sar[ii])*(double)(sar[ii]) ; npos++ ;
00097                } else if( sar[ii] < 0 )
00098                  nneg++ ;
00099             }
00100          break ;
00101    
00102          case MRI_byte:                       /* there are no negative bytes */
00103             bar =  DSET_ARRAY(dset,iv) ;
00104             for( ii=0 ; ii < nvox ; ii++ ){
00105                if( bar[ii] > 0 ){
00106                   hist[bar[ii]]++ ; dsum += (double)(bar[ii])*(double)(bar[ii]) ; npos++ ;
00107                }
00108             }
00109          break ;
00110       }
00111    }
00112 
00113    DSET_unload(dset) ;
00114 
00115    if( verb ) fprintf(stderr,"++ %d positive voxels\n",npos) ;
00116 
00117    if( npos <= 999 ){
00118       fprintf(stderr,"** TOO few positive voxel (%d)!\n",npos) ; exit(1) ;
00119    }
00120    if( nneg > 0 ){
00121       fprintf(stderr,"++ WARNING: ignored %d negative voxels!\n",nneg) ;
00122    }
00123 
00124    /*-- initialize cut position to include upper 65% of positive voxels --*/
00125 
00126    qq = 0.65 * npos ; ib = rint(0.5*sqrt(dsum/npos)) ;
00127    for( kk=0,ii=nhist-1 ; ii >= ib && kk < qq ; ii-- ) kk += hist[ii] ;
00128 
00129    /*-- algorithm --*/
00130 
00131    ncut = ii ; qq = 0 ;
00132    do{
00133       for( npos=0,ii=ncut ; ii < nhist ; ii++ ) npos += hist[ii] ; /* number >= cut */
00134       nhalf = npos/2 ;
00135       for( kk=0,ii=ncut ; ii < nhist && kk < nhalf ; ii++ )        /* find median */
00136          kk += hist[ii] ;
00137       if( verb )
00138          fprintf(stderr,"++ Clip=%d  Median=%d  Number>=Clip=%d\n",ncut,ii,npos) ;
00139       nold = ncut ;
00140       ncut = mfrac * ii ;                                          /* new cut */
00141       qq++ ;
00142    } while( qq < 20 && ncut != nold ) ;
00143 
00144    fac = DSET_BRICK_FACTOR(dset,0) ;
00145    if( fac > 0.0 )
00146       printf("%g\n",ncut*fac) ;
00147    else
00148       printf("%d\n",ncut) ;
00149 
00150    exit(0) ;
00151 }