1 #include "mrilib.h"
2
3 /*--------------------------------------------------------------------------*/
4
mri_topclip(MRI_IMAGE * im)5 float mri_topclip( MRI_IMAGE *im ) /* 28 Sep 2006 */
6 {
7 float cv , dv ;
8 ENTRY("mri_topclip") ;
9 cv = 3.11f * THD_cliplevel( im , 0.511f ) ;
10 dv = (float)mri_max( im ) ;
11 cv = MIN(cv,dv) ; RETURN(cv) ;
12 }
13
14 /*--------------------------------------------------------------------------
15 12 Aug 2001: compare with 3dClipLevel.c
16 - compute a clipping level for an image, to eliminate non-brain voxels
17 - negative voxels are ignored
18 - 05 Nov 2001: increased size of hist array to nhist+1 (from nhist), to
19 store properly elements [0..nhist] (d'oh).
20 ----------------------------------------------------------------------------*/
21
THD_cliplevel(MRI_IMAGE * im,float mfrac)22 float THD_cliplevel( MRI_IMAGE *im , float mfrac )
23 {
24 MRI_IMAGE *lim ;
25 double fac , sfac=1.0 , dsum ;
26 int nvox , *hist , ii,npos=0 , ncut,kk,ib , qq,nold ;
27 short *sar ;
28 byte *bar ;
29 int nhist , nneg=0 , nhalf ;
30
31 ENTRY("THD_cliplevel") ;
32 if( im == NULL ) RETURN(0.0f) ;
33
34 if( mfrac <= 0.0f || mfrac >= 0.99f ) mfrac = 0.50f ;
35
36 /*-- allocate histogram --*/
37
38 switch( im->kind ){
39 case MRI_short: nhist = 32767 ; lim = im ; break ;
40 case MRI_byte : nhist = 255 ; lim = im ; break ;
41 case MRI_float: nhist = 10000 ; lim = im ; break ; /* 20 Dec 2006 */
42
43 default:
44 if( im->kind == MRI_rgb ){
45 nhist = 255 ;
46 } else {
47 fac = mri_maxabs(im) ; if( fac < 1.0e-100 ) RETURN(0.0f) ;
48 sfac = 32767.0/fac ; nhist = 32767 ;
49 }
50 lim = mri_to_short( sfac , im ) ;
51 break ;
52 }
53
54 hist = (int *) calloc(sizeof(int),nhist+1) ; /* 05 Nov 2001: +1 */
55 nvox = lim->nvox ;
56
57 /*-- make histogram --*/
58
59 dsum = 0.0 ;
60 switch( lim->kind ){
61 default: break ;
62
63 case MRI_float:{ /* 20 Dec 2006: do the float->int conversion inline */
64 float *far = MRI_FLOAT_PTR(lim) ;
65 fac = mri_max(im) ; if( fac < 1.e-100 ){ free(hist); RETURN(0.0f); }
66 sfac = nhist / fac ;
67 for( ii=0 ; ii < nvox ; ii++ ){
68 if( far[ii] > 0.0f ){
69 kk = (int)(sfac*far[ii]+0.499) ;
70 if( kk <= nhist ){
71 hist[kk]++ ;
72 dsum += ((double)kk) * ((double)(kk)) ; npos++ ;
73 }
74 }
75 }
76 }
77 break ;
78
79 case MRI_short:
80 sar = MRI_SHORT_PTR(lim) ;
81 for( ii=0 ; ii < nvox ; ii++ ){
82 if( sar[ii] > 0 && sar[ii] <= nhist ){
83 hist[sar[ii]]++ ;
84 dsum += (double)(sar[ii])*(double)(sar[ii]); npos++;
85 } else if( sar[ii] < 0 )
86 nneg++ ;
87 }
88 break ;
89
90 case MRI_byte: /* there are no negative bytes */
91 bar = MRI_BYTE_PTR(lim) ;
92 for( ii=0 ; ii < nvox ; ii++ ){
93 if( bar[ii] > 0 ){
94 hist[bar[ii]]++ ;
95 dsum += (double)(bar[ii])*(double)(bar[ii]); npos++;
96 }
97 }
98 break ;
99 }
100
101 if( lim != im ) mri_free(lim) ;
102
103 if( npos <= 222 ){ free(hist); RETURN(0.0f); }
104
105 /*-- initialize cut position to include upper 65% of positive voxels --*/
106
107 qq = 0.65f * npos ; ib = rint(0.5*sqrt(dsum/npos)) ;
108 for( kk=0,ii=nhist-1 ; ii >= ib && kk < qq ; ii-- ) kk += hist[ii] ;
109
110 /*-- median adjustment algorithm:
111 we find a cut level so that it equals mfrac times
112 the median of all the values above the cut level. ---*/
113
114 ncut = ii ; qq = 0 ; /* qq is iteration count */
115 do{
116 for( npos=0,ii=ncut; ii < nhist; ii++ ) npos += hist[ii]; /* num >= cut */
117 nhalf = npos/2 ; /* half the number >= cut */
118 for( kk=0,ii=ncut ; ii < nhist && kk < nhalf ; ii++ ) /* find median */
119 kk += hist[ii] ; /* loop ends at ii=median bin */
120 nold = ncut ;
121 ncut = mfrac * ii ; /* new cut */
122 qq++ ;
123 } while( qq < 66 && ncut != nold ) ; /* usually only 2-3 iterations needed */
124
125 free(hist) ;
126
127 fac = ncut / sfac ;
128 if( fac > 1.e+38 ) fac = 1.e+38 ;
129 RETURN( (float)fac ) ;
130 }
131
132 /*-------------------------------------------------------------------------*/
133 /* (1) apply THD_cliplevel() algorithm above to absolute values.
134 (2) also, if mfrac < 0, then find 90% point on CDF of absolute values
135 and return smaller of this or the THD_cliplevel() result.
136 Purpose is to find some reasonable place to threshold an image for
137 visual effect only.
138 ---------------------------------------------------------------------------*/
139
THD_cliplevel_abs(MRI_IMAGE * im,float mfrac)140 float THD_cliplevel_abs( MRI_IMAGE *im , float mfrac )
141 {
142 MRI_IMAGE *fim ;
143 register float *far ;
144 register int ii ;
145 float val,tv ; int dotwo=0 ;
146
147 ENTRY("THD_cliplevel_abs") ;
148 if( im == NULL ) RETURN(0.0f) ;
149 fim = mri_to_float(im) ; if( fim == NULL ) RETURN(0.0f) ;
150 far = MRI_FLOAT_PTR(fim) ;
151 for( ii=0 ; ii < fim->nvox ; ii++ ) far[ii] = fabsf(far[ii]) ;
152 if( mfrac < 0.0f ){ dotwo = 1; mfrac = -mfrac; }
153 val = THD_cliplevel( fim , mfrac ) ;
154
155 if( dotwo ){
156 qsort_float( fim->nvox , far ) ;
157 ii = (int)(0.9*fim->nvox) ; tv = far[ii] ;
158 if( tv == 0.0f ){
159 for( ; ii < fim->nvox && far[ii] == 0.0f ; ii++ ) ; /* nada */
160 if( ii < fim->nvox ) tv = far[ii] ;
161 }
162 if( val > tv && tv > 0.0f ) val = tv ;
163 }
164
165 mri_free(fim) ; RETURN(val) ;
166 }
167
168 /*-------------------------------------------------------------------------*/
169 /*! Cliplevel for part of an image. Quick and easy to write!
170 Not very efficient, but this isn't a big CPU sink. [24 Oct 2006] */
171
THD_cliplevel_partial(MRI_IMAGE * im,float mfrac,int xa,int xb,int ya,int yb,int za,int zb)172 float THD_cliplevel_partial( MRI_IMAGE *im , float mfrac ,
173 int xa,int xb, int ya,int yb, int za,int zb )
174 {
175 MRI_IMAGE *qim ; float val ;
176
177 ENTRY("THD_cliplevel_partial") ;
178 qim = mri_cut_3D( im , xa,xb , ya,yb , za,zb ) ;
179 val = THD_cliplevel( qim , mfrac ) ;
180 mri_free(qim) ; RETURN(val) ;
181 }
182
183 /*-------------------------------------------------------------------------*/
184
185 typedef struct {
186 float clip_000, clip_100, clip_010, clip_110,
187 clip_001, clip_101, clip_011, clip_111 ;
188 float x0,x1,dxi , y0,y1,dyi , z0,z1,dzi ;
189 float clip_min , clip_max ;
190 } clipvec ;
191
192 /*-------------------------------------------------------------------------*/
193 /*! Get cliplevel for each octant about the center-of-mass. [24 Oct 2006] */
194
get_octant_clips(MRI_IMAGE * im,float mfrac)195 static clipvec get_octant_clips( MRI_IMAGE *im , float mfrac )
196 {
197 float xcm,ycm,zcm , sum,val , clip_min , clip_max ;
198 int ii,jj,kk , nx,ny,nz , ic,jc,kc , it,jt,kt , ijk ;
199 int icp,icm , jcp,jcm , kcp,kcm ;
200 clipvec cv ;
201
202 ENTRY("get_octant_clips") ;
203
204 cv.clip_000 = -1 ; /* flags error return */
205
206 if( im == NULL ) RETURN(cv) ;
207
208 nx = im->nx ; ny = im->ny ; nz = im->nz ;
209 it = nx-1 ; jt = ny-1 ; kt = nz-1 ;
210
211 /* compute CM of image */
212
213 mri_get_cmass_3D( im, &xcm, &ycm, &zcm ) ;
214
215 ic = (int)rint(xcm); jc = (int)rint(ycm); kc = (int)rint(zcm);
216
217 /* compute cliplevel in each octant about the CM */
218
219 val = 0.333f * THD_cliplevel( im,mfrac ) ;
220
221 ii = (int)rint(0.01*nx); if( ii < 1 ) ii = 1 ; /* 1% quadrant overlap */
222 jj = (int)rint(0.01*ny); if( jj < 1 ) jj = 1 ;
223 kk = (int)rint(0.01*nz); if( kk < 1 ) kk = 1 ;
224 icm = ic-ii ; icp = ic+ii ; if( icm < 0 ) icm = 0; if( icp > it ) icp = it ;
225 jcm = jc-jj ; jcp = jc+jj ; if( jcm < 0 ) jcm = 0; if( jcp > jt ) jcp = jt ;
226 kcm = kc-kk ; kcp = kc+kk ; if( kcm < 0 ) kcm = 0; if( kcp > kt ) kcp = kt ;
227
228 cv.clip_000 = THD_cliplevel_partial( im,mfrac, 0 ,icp, 0 ,jcp, 0 ,kcp );
229 cv.clip_100 = THD_cliplevel_partial( im,mfrac, icm,it , 0 ,jcp, 0 ,kcp );
230 cv.clip_010 = THD_cliplevel_partial( im,mfrac, 0 ,icp, jcm,jt , 0 ,kcp );
231 cv.clip_110 = THD_cliplevel_partial( im,mfrac, icm,it , jcm,jt , 0 ,kcp );
232 cv.clip_001 = THD_cliplevel_partial( im,mfrac, 0 ,icp, 0 ,jcp, kcm,kt );
233 cv.clip_101 = THD_cliplevel_partial( im,mfrac, icm,it , 0 ,jcp, kcm,kt );
234 cv.clip_011 = THD_cliplevel_partial( im,mfrac, 0 ,icp, jcm,jt , kcm,kt );
235 cv.clip_111 = THD_cliplevel_partial( im,mfrac, icm,it , jcm,jt , kcm,kt );
236
237 if( cv.clip_000 < val ) cv.clip_000 = val ; /* don't let */
238 if( cv.clip_100 < val ) cv.clip_100 = val ; /* clip levels */
239 if( cv.clip_010 < val ) cv.clip_010 = val ; /* get too */
240 if( cv.clip_110 < val ) cv.clip_110 = val ; /* small! */
241 if( cv.clip_001 < val ) cv.clip_001 = val ;
242 if( cv.clip_101 < val ) cv.clip_101 = val ;
243 if( cv.clip_011 < val ) cv.clip_011 = val ;
244 if( cv.clip_111 < val ) cv.clip_111 = val ;
245
246 clip_min = cv.clip_000 ;
247 clip_min = MIN(clip_min,cv.clip_100) ;
248 clip_min = MIN(clip_min,cv.clip_010) ;
249 clip_min = MIN(clip_min,cv.clip_110) ;
250 clip_min = MIN(clip_min,cv.clip_001) ;
251 clip_min = MIN(clip_min,cv.clip_101) ;
252 clip_min = MIN(clip_min,cv.clip_011) ;
253 clip_min = MIN(clip_min,cv.clip_111) ; cv.clip_min = clip_min ;
254
255 clip_max = cv.clip_000 ;
256 clip_max = MAX(clip_max,cv.clip_100) ;
257 clip_max = MAX(clip_max,cv.clip_010) ;
258 clip_max = MAX(clip_max,cv.clip_110) ;
259 clip_max = MAX(clip_max,cv.clip_001) ;
260 clip_max = MAX(clip_max,cv.clip_101) ;
261 clip_max = MAX(clip_max,cv.clip_011) ;
262 clip_max = MAX(clip_max,cv.clip_111) ; cv.clip_max = clip_max ;
263
264 /* (x0,y0,z0) = center of lowest octant
265 (x1,y1,z1) = center of highest octant */
266
267 cv.x0 = 0.5*ic ; cv.x1 = 0.5*(ic+it) ;
268 cv.y0 = 0.5*jc ; cv.y1 = 0.5*(jc+jt) ;
269 cv.z0 = 0.5*kc ; cv.z1 = 0.5*(kc+kt) ;
270 cv.dxi = (cv.x1 > cv.x0) ? 1.0/(cv.x1-cv.x0) : 0.0 ;
271 cv.dyi = (cv.y1 > cv.y0) ? 1.0/(cv.y1-cv.y0) : 0.0 ;
272 cv.dzi = (cv.z1 > cv.z0) ? 1.0/(cv.z1-cv.z0) : 0.0 ;
273
274 RETURN(cv) ;
275 }
276
277 /*--------------------------------------------------------------------------*/
278 /*! Return the cliplevel at any point,
279 tri-linearly interpolated between octant centers. [24 Oct 2006] */
280
pointclip(int ii,int jj,int kk,clipvec * cv)281 static INLINE float pointclip( int ii, int jj, int kk , clipvec *cv )
282 {
283 float x1,y1,z1 , x0,y0,z0 , val ;
284
285 /* get relative position in box defined by octant centers */
286
287 x1 = (ii-cv->x0)*cv->dxi; if(x1 < 0.0f) x1=0.0f; else if(x1 > 1.0f) x1=1.0f;
288 y1 = (jj-cv->y0)*cv->dyi; if(y1 < 0.0f) y1=0.0f; else if(y1 > 1.0f) y1=1.0f;
289 z1 = (kk-cv->z0)*cv->dzi; if(z1 < 0.0f) z1=0.0f; else if(z1 > 1.0f) z1=1.0f;
290
291 x0 = 1.0f-x1 ; y0 = 1.0f-y1 ; z0 = 1.0f-z1 ;
292
293 val = cv->clip_000 * x0*y0*z0 + cv->clip_100 * x1*y0*z0
294 + cv->clip_010 * x0*y1*z0 + cv->clip_110 * x1*y1*z0
295 + cv->clip_001 * x0*y0*z1 + cv->clip_101 * x1*y0*z1
296 + cv->clip_011 * x0*y1*z1 + cv->clip_111 * x1*y1*z1 ;
297 return val ;
298 }
299
300 /*--------------------------------------------------------------------------*/
301 /*! Return a cliplevel that varies gradually across the image.
302 At this time [24 Oct 2006], the clipping varies trilinearly
303 between octants. A fancier method would use higher order
304 polynomials, but I can't really be bothered with that now.
305 ----------------------------------------------------------------------------*/
306
THD_cliplevel_gradual(MRI_IMAGE * im,float mfrac)307 MRI_IMAGE * THD_cliplevel_gradual( MRI_IMAGE *im , float mfrac )
308 {
309 int ii,jj,kk,ijk , nx,ny,nz ;
310 MRI_IMAGE *cim ; float *car ;
311 clipvec bvec ;
312
313 ENTRY("THD_cliplevel_gradual") ;
314 if( im == NULL ) RETURN(NULL) ;
315
316 bvec = get_octant_clips( im , mfrac ) ;
317 if( bvec.clip_000 < 0.0 ) RETURN(NULL) ;
318
319 nx = im->nx ; ny = im->ny ; nz = im->nz ;
320
321 cim = mri_new_conforming( im , MRI_float ) ;
322 car = MRI_FLOAT_PTR(cim) ;
323
324 for( ijk=kk=0 ; kk < nz ; kk++ ){
325 for( jj=0 ; jj < ny ; jj++ ){
326 for( ii=0 ; ii < nx ; ii++,ijk++ ){
327 car[ijk] = pointclip( ii,jj,kk , &bvec ) ; /* cliplevel here */
328 }}}
329
330 RETURN(cim) ;
331 }
332
333 /*--------------------------------------------------------------------------*/
334
THD_cliplevel_search(MRI_IMAGE * im)335 float THD_cliplevel_search( MRI_IMAGE *im ) /* experimental */
336 {
337 #define NC 10
338 #define DC 0.05f
339 #define CB 0.10f
340 int qc , nmask[NC] ; float cc ; byte *mask ;
341
342 THD_automask_verbose(0) ;
343 INFO_message("\nTHD_cliplevel_search:") ;
344 for( qc=0 ; qc < NC ; qc++ ){
345 cc = DC * qc + CB ;
346 THD_automask_set_clipfrac(cc) ;
347 THD_automask_set_cheapo(1) ;
348 mask = mri_automask_image(im) ;
349 nmask[qc] = THD_countmask( im->nvox , mask ) ;
350 free(mask) ;
351 ININFO_message(" clfrac=%.2f nmask=%d (%.1f%%)",cc,nmask[qc],(100.0f*nmask[qc])/(float)(im->nvox)) ;
352 }
353 THD_automask_set_cheapo(0) ;
354 return 0.0f ;
355 }
356