1 /*****************************************************************************
2 Major portions of this software are copyrighted by the Medical College
3 of Wisconsin, 1994-2000, and are released under the Gnu General Public
4 License, Version 2. See the file README.Copyright for details.
5 ******************************************************************************/
6
7 #include "mrilib.h"
8
9 /*---------------------------------------------------------------------------------*/
10
11 #define NPMAX 128
12
main(int argc,char * argv[])13 int main( int argc , char * argv[] )
14 {
15 int iarg ;
16 THD_3dim_dataset * dset ;
17 float bper = 60.0 , bmin = 1 ;
18
19 int nxyz , ii , kk , nbin , sval , sum , nbot , a,b,c , nbase,npeak,ntop , nvox ;
20 int * fbin ;
21 short * sfim ;
22 int kmin[NPMAX] , kmax[NPMAX] ;
23 int nmin , nmax ;
24
25 /*-- help? --*/
26
27 if( argc < 2 || strncmp(argv[1],"-help",4) == 0 ){
28 fprintf(stderr,
29 "Find the clipping point on a histogram of a 3D brick\n"
30 "Usage: 3dhclip [options] dataset\n"
31 "\n"
32 "Options:\n"
33 " -bper bb Means to take the base percentage point\n"
34 " on the cumulative histogram as 'bb'\n"
35 " [default bb = 60]\n"
36 " -bmin cc Means to take the minimum value to consider\n"
37 " as 'cc' [default cc = 1]\n"
38 ) ;
39
40 printf("\n" MASTER_SHORTHELP_STRING ) ;
41
42 PRINT_COMPILE_DATE ; exit(0) ;
43 }
44
45 /*-- process options --*/
46
47 iarg = 1 ;
48 while( iarg < argc && argv[iarg][0] == '-' ){
49
50 if( strcmp(argv[iarg],"-bper") == 0 ){
51 bper = strtod( argv[++iarg] , NULL ) ;
52 if(bper<1 || bper>99){fprintf(stderr,"** Illegal -bper value\n");exit(1);}
53 iarg++ ; continue ;
54 }
55
56 if( strcmp(argv[iarg],"-bmin") == 0 ){
57 bmin = strtod( argv[++iarg] , NULL ) ;
58 if(bmin<0){fprintf(stderr,"** Illegal -bmin value\n");exit(1);}
59 iarg++ ; continue ;
60 }
61
62 fprintf(stderr,"** Illegal option: %s\n",argv[iarg]) ; exit(1) ;
63 }
64
65 /*-- load dataset --*/
66
67 if( iarg >= argc ){fprintf(stderr,"** No dataset?!\n");exit(1);}
68
69 dset = THD_open_dataset(argv[iarg]) ;
70 if( dset == NULL ){
71 fprintf(stderr,"** Can't open dataset %s\n",argv[iarg]) ;
72 exit(1) ;
73 }
74
75 if( DSET_BRICK_TYPE(dset,0) != MRI_short || DSET_BRICK_FACTOR(dset,0) != 0.0 ){
76 fprintf(stderr,"** Program only works with short-value datasets!\n") ;
77 exit(1) ;
78 }
79
80 DSET_load(dset) ;
81 sfim = DSET_ARRAY(dset,0) ;
82 if( sfim == NULL ){fprintf(stderr,"** Can't load dataset brick\n");exit(1);}
83
84 /*-- make histogram of shorts --*/
85
86
87 fbin = (int *) malloc( sizeof(int) * 32768 ) ;
88 for( kk=0 ; kk < 32768 ; kk++ ) fbin[kk] = 0 ;
89
90 nvox = 0 ; nxyz = DSET_NVOX(dset) ;
91
92 for( ii=0 ; ii < nxyz ; ii++ ){
93 kk = sfim[ii] ; if( kk >= 0 ){ fbin[kk]++ ; nvox++ ; }
94 }
95
96 DSET_unload(dset) ;
97
98 /*-- find largest value --*/
99
100 for( kk=32767 ; kk > 0 ; kk-- ) if( fbin[kk] > 0 ) break ;
101 if( kk == 0 ){fprintf(stderr,"** All voxels are zero!\n");exit(1);}
102 nbin = kk+1 ;
103
104 /*-- find bper point in cumulative distribution --*/
105
106 sval = 0.01 * bper * nvox ;
107 sum = 0 ;
108 for( kk=0 ; kk < nbin ; kk++ ){
109 sum += fbin[kk] ; if( sum >= sval ) break ;
110 }
111 nbot = kk ; if( nbot == 0 ) nbot = 1 ; if( bmin > nbot ) nbot = bmin ;
112 if( nbot >= nbin-9 ){fprintf(stderr,"** Base point on histogram too high\n");exit(1);}
113
114 /*-- smooth histogram --*/
115
116 b = fbin[nbot-1] ; c = fbin[nbot] ;
117 for( kk=nbot ; kk < nbin ; kk++ ){
118 a = b ; b = c ; c = fbin[kk+1] ; fbin[kk] = 0.25*(a+c+2*b) ;
119 }
120
121 /*-- find minima and maxima above bper point --*/
122
123 nmin = nmax = 0 ;
124 for( kk=nbot+1 ; kk < nbin ; kk++ ){
125 if( fbin[kk] < fbin[kk-1] && fbin[kk] < fbin[kk+1] && nmin < NPMAX ){
126 kmin[nmin++] = kk ;
127 } else if( fbin[kk] > fbin[kk-1] && fbin[kk] > fbin[kk+1] && nmax < NPMAX ){
128 kmax[nmax++] = kk ;
129 }
130 }
131
132 #if 0
133 for( kk=0 ; kk < nmin ; kk++ )
134 printf("Min: %d has %d\n",kmin[kk],fbin[kmin[kk]]) ;
135
136 for( kk=0 ; kk < nmax ; kk++ )
137 printf("Max: %d has %d\n",kmax[kk],fbin[kmax[kk]]) ;
138 #endif
139
140 /*-- find the largest two maxima --*/
141
142 if( nmax == 0 ){fprintf(stderr,"** No maxima above base point\n");exit(1);}
143
144 if( nmax <= 2 ){
145 npeak = kmax[0] ; ntop = 0 ;
146 } else {
147 int f1,f2 , k1,k2 , fk , klow,kup ;
148
149 k1 = 0 ; f1 = fbin[kmax[0]] ;
150 k2 = 1 ; f2 = fbin[kmax[1]] ;
151 if( f1 < f2 ){
152 k1 = 1 ; f1 = fbin[kmax[1]] ;
153 k2 = 0 ; f2 = fbin[kmax[0]] ;
154 }
155
156 for( kk=2 ; kk < nmax ; kk++ ){
157 fk = fbin[kmax[kk]] ;
158 if( fk > f1 ){
159 f2 = f1 ; k2 = k1 ;
160 f1 = fk ; k1 = kk ;
161 } else if( fk > f2 ){
162 f2 = fk ; k2 = kk ;
163 }
164 }
165 npeak = MIN( kmax[k1] , kmax[k2] ) ; /* smaller bin of the 2 top peaks */
166
167 /* find valley between 2 peaks */
168
169 ntop = MAX( kmax[k1] , kmax[k2] ) ;
170
171 fk = fbin[ntop] ; klow = ntop ;
172 for( kk=ntop-1 ; kk >= npeak ; kk-- ){
173 if( fbin[kk] < fk ){ fk = fbin[kk] ; klow = kk ; }
174 }
175 fk = 0.16 * fk ; kup = MIN( nbin-1 , ntop+3*(ntop-klow+2) ) ;
176 for( kk=ntop+1 ; kk <= kup ; kk++ ) if( fbin[kk] < fk ) break ;
177 ntop = kk ;
178 }
179
180 for( kk=npeak-1 ; kk > 0 ; kk-- )
181 if( fbin[kk] < fbin[kk-1] && fbin[kk] < fbin[kk+1] ) break ;
182 nbase = kk ;
183
184 printf("dataset: %s -- ",argv[iarg]) ;
185 printf("base = %d peak = %d top = %d\n",nbase,npeak,ntop) ;
186
187 exit(0) ;
188 }
189