00014 {
00015    int iarg ;
00016    THD_3dim_dataset * dset ;
00017    float bper = 60.0 , bmin = 1 ;
00018 
00019    int nxyz , ii , kk , nbin , sval , sum , nbot , a,b,c , nbase,npeak,ntop , nvox ;
00020    int * fbin ;
00021    short * sfim ;
00022    int   kmin[NPMAX] , kmax[NPMAX] ;
00023    int   nmin        , nmax        ;
00024 
00025    /*-- help? --*/
00026 
00027    if( argc < 2 || strncmp(argv[1],"-help",4) == 0 ){
00028       fprintf(stderr,
00029              "Find the clipping point on a histogram of a 3D brick\n"
00030              "Usage: 3dhclip [options] dataset\n"
00031              "\n"
00032              "Options:\n"
00033              "  -bper bb   Means to take the base percentage point\n"
00034              "              on the cumulative histogram as 'bb'\n"
00035              "              [default bb = 60]\n"
00036              "  -bmin cc   Means to take the minimum value to consider\n"
00037              "              as 'cc' [default cc = 1]\n"
00038          ) ;
00039 
00040       printf("\n" MASTER_SHORTHELP_STRING ) ;
00041 
00042       exit(0) ;
00043    }
00044 
00045    /*-- process options --*/
00046 
00047    iarg = 1 ;
00048    while( iarg < argc && argv[iarg][0] == '-' ){
00049 
00050       if( strcmp(argv[iarg],"-bper") == 0 ){
00051          bper = strtod( argv[++iarg] , NULL ) ;
00052          if(bper<1 || bper>99){fprintf(stderr,"** Illegal -bper value\n");exit(1);}
00053          iarg++ ; continue ;
00054       }
00055 
00056       if( strcmp(argv[iarg],"-bmin") == 0 ){
00057          bmin = strtod( argv[++iarg] , NULL ) ;
00058          if(bmin<0){fprintf(stderr,"** Illegal -bmin value\n");exit(1);}
00059          iarg++ ; continue ;
00060       }
00061 
00062       fprintf(stderr,"** Illegal option: %s\n",argv[iarg]) ; exit(1) ;
00063    }
00064 
00065    /*-- load dataset --*/
00066 
00067    if( iarg >=  argc ){fprintf(stderr,"** No dataset?!\n");exit(1);}
00068 
00069    dset = THD_open_dataset(argv[iarg]) ;
00070    if( dset == NULL ){
00071       fprintf(stderr,"** Can't open dataset %s\n",argv[iarg]) ;
00072       exit(1) ;
00073    }
00074 
00075    if( DSET_BRICK_TYPE(dset,0) != MRI_short || DSET_BRICK_FACTOR(dset,0) != 0.0 ){
00076      fprintf(stderr,"** Program only works with short-value datasets!\n") ;
00077      exit(1) ;
00078    }
00079 
00080    DSET_load(dset) ;
00081    sfim = DSET_ARRAY(dset,0) ;
00082    if( sfim == NULL ){fprintf(stderr,"** Can't load dataset brick\n");exit(1);}
00083 
00084    /*-- make histogram of shorts --*/
00085 
00086 
00087    fbin = (int *) malloc( sizeof(int) * 32768 ) ;
00088    for( kk=0 ; kk < 32768 ; kk++ ) fbin[kk] = 0 ;
00089 
00090    nvox = 0 ; nxyz = DSET_NVOX(dset) ;
00091 
00092    for( ii=0 ; ii < nxyz ; ii++ ){
00093       kk = sfim[ii] ; if( kk >= 0 ){ fbin[kk]++ ; nvox++ ; }
00094    }
00095 
00096    DSET_unload(dset) ;
00097 
00098    /*-- find largest value --*/
00099 
00100    for( kk=32767 ; kk > 0 ; kk-- ) if( fbin[kk] > 0 ) break ;
00101    if( kk == 0 ){fprintf(stderr,"** All voxels are zero!\n");exit(1);}
00102    nbin = kk+1 ;
00103 
00104    /*-- find bper point in cumulative distribution --*/
00105 
00106    sval = 0.01 * bper * nvox ;
00107    sum  = 0 ;
00108    for( kk=0 ; kk < nbin ; kk++ ){
00109       sum += fbin[kk] ; if( sum >= sval ) break ;
00110    }
00111    nbot = kk ; if( nbot == 0 ) nbot = 1 ; if( bmin > nbot ) nbot = bmin ;
00112    if( nbot >= nbin-9 ){fprintf(stderr,"** Base point on histogram too high\n");exit(1);}
00113 
00114    /*-- smooth histogram --*/
00115 
00116    b = fbin[nbot-1] ; c = fbin[nbot] ;
00117    for( kk=nbot ; kk < nbin ; kk++ ){
00118       a = b ; b = c ; c = fbin[kk+1] ; fbin[kk] = 0.25*(a+c+2*b) ;
00119    }
00120 
00121    /*-- find minima and maxima above bper point --*/
00122 
00123    nmin = nmax = 0 ;
00124    for( kk=nbot+1 ; kk < nbin ; kk++ ){
00125       if( fbin[kk] < fbin[kk-1] && fbin[kk] < fbin[kk+1] && nmin < NPMAX ){
00126          kmin[nmin++] = kk ;
00127       } else if( fbin[kk] > fbin[kk-1] && fbin[kk] > fbin[kk+1] && nmax < NPMAX ){
00128          kmax[nmax++] = kk ;
00129       }
00130    }
00131 
00132 #if 0
00133    for( kk=0 ; kk < nmin ; kk++ )
00134       printf("Min: %d has %d\n",kmin[kk],fbin[kmin[kk]]) ;
00135 
00136    for( kk=0 ; kk < nmax ; kk++ )
00137       printf("Max: %d has %d\n",kmax[kk],fbin[kmax[kk]]) ;
00138 #endif
00139 
00140    /*-- find the largest two maxima --*/
00141 
00142    if( nmax == 0 ){fprintf(stderr,"** No maxima above base point\n");exit(1);}
00143 
00144    if( nmax <= 2 ){
00145       npeak = kmax[0] ; ntop = 0 ;
00146    } else {
00147       int f1,f2 , k1,k2 , fk , klow,kup ;
00148 
00149       k1 = 0 ; f1 = fbin[kmax[0]] ;
00150       k2 = 1 ; f2 = fbin[kmax[1]] ;
00151       if( f1 < f2 ){
00152          k1 = 1 ; f1 = fbin[kmax[1]] ;
00153          k2 = 0 ; f2 = fbin[kmax[0]] ;
00154       }
00155 
00156       for( kk=2 ; kk < nmax ; kk++ ){
00157          fk = fbin[kmax[kk]] ;
00158          if( fk > f1 ){
00159             f2 = f1 ; k2 = k1 ;
00160             f1 = fk ; k1 = kk ;
00161          } else if( fk > f2 ){
00162             f2 = fk ; k2 = kk ;
00163          }
00164       }
00165       npeak = MIN( kmax[k1] , kmax[k2] ) ;  /* smaller bin of the 2 top peaks */
00166 
00167       /* find valley between 2 peaks */
00168 
00169       ntop  = MAX( kmax[k1] , kmax[k2] ) ;
00170 
00171       fk = fbin[ntop] ; klow = ntop ;
00172       for( kk=ntop-1 ; kk >= npeak ; kk-- ){
00173          if( fbin[kk] < fk ){ fk = fbin[kk] ; klow = kk ; }
00174       }
00175       fk = 0.16 * fk ; kup = MIN( nbin-1 , ntop+3*(ntop-klow+2) ) ;
00176       for( kk=ntop+1 ; kk <= kup ; kk++ ) if( fbin[kk] < fk ) break ;
00177       ntop = kk ;
00178    }
00179 
00180    for( kk=npeak-1 ; kk > 0 ; kk-- )
00181       if( fbin[kk] < fbin[kk-1] && fbin[kk] < fbin[kk+1] ) break ;
00182    nbase = kk ;
00183 
00184    printf("dataset: %s -- ",argv[iarg]) ;
00185    printf("base = %d  peak = %d  top = %d\n",nbase,npeak,ntop) ;
00186 
00187    exit(0) ;
00188 }