Doxygen Source Code Documentation
Main Page Alphabetical List Data Structures File List Data Fields Globals Search
3dhclip.c
Go to the documentation of this file.00001
00002
00003
00004
00005
00006
00007 #include "mrilib.h"
00008
00009
00010
00011 #define NPMAX 128
00012
00013 int main( int argc , char * argv[] )
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
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
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
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
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
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
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
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
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
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] ) ;
00166
00167
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 }