#include "mrilib.h"
/*--------------------------------------------------------------------------*/
float mri_topclip( MRI_IMAGE *im ) /* 28 Sep 2006 */
{
float cv , dv ;
ENTRY("mri_topclip") ;
cv = 3.11f * THD_cliplevel( im , 0.511f ) ;
dv = (float)mri_max( im ) ;
cv = MIN(cv,dv) ; RETURN(cv) ;
}
/*--------------------------------------------------------------------------
12 Aug 2001: compare with 3dClipLevel.c
- compute a clipping level for an image, to eliminate non-brain voxels
- negative voxels are ignored
- 05 Nov 2001: increased size of hist array to nhist+1 (from nhist), to
store properly elements [0..nhist] (d'oh).
----------------------------------------------------------------------------*/
float THD_cliplevel( MRI_IMAGE *im , float mfrac )
{
MRI_IMAGE *lim ;
float fac , sfac=1.0 ;
double dsum ;
int nvox , *hist , ii,npos=0 , ncut,kk,ib , qq,nold ;
short *sar ;
byte *bar ;
int nhist , nneg=0 , nhalf ;
ENTRY("THD_cliplevel") ;
if( im == NULL ) RETURN(0.0f) ;
if( mfrac <= 0.0f || mfrac >= 0.99f ) mfrac = 0.50f ;
/*-- allocate histogram --*/
switch( im->kind ){
case MRI_short: nhist = 32767 ; lim = im ; break ;
case MRI_byte : nhist = 255 ; lim = im ; break ;
case MRI_float: nhist = 10000 ; lim = im ; break ; /* 20 Dec 2006 */
default:
if( im->kind == MRI_rgb ){
nhist = 255 ;
} else {
fac = (float)mri_maxabs(im) ; if( fac == 0.0f ) RETURN(0.0f) ;
sfac = 32767.0f/fac ; nhist = 32767 ;
}
lim = mri_to_short( sfac , im ) ;
break ;
}
hist = (int *) calloc(sizeof(int),nhist+1) ; /* 05 Nov 2001: +1 */
nvox = lim->nvox ;
/*-- make histogram --*/
dsum = 0.0 ;
switch( lim->kind ){
default: break ;
case MRI_float:{ /* 20 Dec 2006: do the float->int conversion inline */
float *far = MRI_FLOAT_PTR(lim) ;
fac = (float)mri_max(im) ;
if( fac <= 0.0f ){ free(hist); RETURN(0.0f); }
sfac = nhist / fac ;
for( ii=0 ; ii < nvox ; ii++ ){
if( far[ii] > 0.0f ){
kk = (int)(sfac * far[ii]+0.499f) ;
if( kk <= nhist ){
hist[kk]++ ;
dsum += ((double)kk) * ((double)(kk)) ; npos++ ;
}
}
}
}
break ;
case MRI_short:
sar = MRI_SHORT_PTR(lim) ;
for( ii=0 ; ii < nvox ; ii++ ){
if( sar[ii] > 0 && sar[ii] <= nhist ){
hist[sar[ii]]++ ;
dsum += (double)(sar[ii])*(double)(sar[ii]); npos++;
} else if( sar[ii] < 0 )
nneg++ ;
}
break ;
case MRI_byte: /* there are no negative bytes */
bar = MRI_BYTE_PTR(lim) ;
for( ii=0 ; ii < nvox ; ii++ ){
if( bar[ii] > 0 ){
hist[bar[ii]]++ ;
dsum += (double)(bar[ii])*(double)(bar[ii]); npos++;
}
}
break ;
}
if( lim != im ) mri_free(lim) ;
if( npos <= 222 ){ free(hist); RETURN(0.0f); }
/*-- initialize cut position to include upper 65% of positive voxels --*/
qq = 0.65f * npos ; ib = rint(0.5*sqrt(dsum/npos)) ;
for( kk=0,ii=nhist-1 ; ii >= ib && kk < qq ; ii-- ) kk += hist[ii] ;
/*-- median adjustment algorithm:
we find a cut level so that it equals mfrac times
the median of all the values above the cut level. ---*/
ncut = ii ; qq = 0 ;
do{
for( npos=0,ii=ncut; ii < nhist; ii++ ) npos += hist[ii]; /* num >= cut */
nhalf = npos/2 ; /* half the number >= cut */
for( kk=0,ii=ncut ; ii < nhist && kk < nhalf ; ii++ ) /* find median */
kk += hist[ii] ; /* loop ends at ii=median bin */
nold = ncut ;
ncut = mfrac * ii ; /* new cut */
qq++ ;
} while( qq < 66 && ncut != nold ) ;
free(hist) ;
RETURN( (ncut/sfac) ) ;
}
/*-------------------------------------------------------------------------*/
float THD_cliplevel_abs( MRI_IMAGE *im , float mfrac )
{
MRI_IMAGE *fim ;
register float *far ;
register int ii ;
float val,tv ; int dotwo ;
ENTRY("THD_cliplevel_abs") ;
if( im == NULL ) RETURN(0.0f) ;
fim = mri_to_float(im) ; if( fim == NULL ) RETURN(0.0f) ;
far = MRI_FLOAT_PTR(fim) ;
for( ii=0 ; ii < fim->nvox ; ii++ ) far[ii] = fabsf(far[ii]) ;
if( mfrac < 0.0f ){ dotwo = 1; mfrac = -mfrac; }
val = THD_cliplevel( fim , mfrac ) ;
if( dotwo ){
qsort_float( fim->nvox , far ) ;
ii = (int)(0.9*fim->nvox) ; tv = far[ii] ;
if( tv == 0.0f ){
for( ; ii < fim->nvox && far[ii] == 0.0f ; ii++ ) ; /* nada */
if( ii < fim->nvox ) tv = far[ii] ;
}
if( val > tv && tv > 0.0f ) val = tv ;
}
mri_free(fim) ; RETURN(val) ;
}
/*-------------------------------------------------------------------------*/
/*! Cliplevel for part of an image. Quick and easy to write!
Not very efficient, but this isn't a big CPU sink. [24 Oct 2006] */
float THD_cliplevel_partial( MRI_IMAGE *im , float mfrac ,
int xa,int xb, int ya,int yb, int za,int zb )
{
MRI_IMAGE *qim ; float val ;
ENTRY("THD_cliplevel_partial") ;
qim = mri_cut_3D( im , xa,xb , ya,yb , za,zb ) ;
val = THD_cliplevel( qim , mfrac ) ;
mri_free(qim) ; RETURN(val) ;
}
/*-------------------------------------------------------------------------*/
typedef struct {
float clip_000, clip_100, clip_010, clip_110,
clip_001, clip_101, clip_011, clip_111 ;
float x0,x1,dxi , y0,y1,dyi , z0,z1,dzi ;
float clip_min , clip_max ;
} clipvec ;
/*-------------------------------------------------------------------------*/
/*! Get cliplevel for each octant about the center-of-mass. [24 Oct 2006] */
static clipvec get_octant_clips( MRI_IMAGE *im , float mfrac )
{
float xcm,ycm,zcm , sum,val , clip_min , clip_max ;
int ii,jj,kk , nx,ny,nz , ic,jc,kc , it,jt,kt , ijk ;
int icp,icm , jcp,jcm , kcp,kcm ;
clipvec cv ;
ENTRY("get_octant_clips") ;
cv.clip_000 = -1 ; /* flags error return */
if( im == NULL ) RETURN(cv) ;
nx = im->nx ; ny = im->ny ; nz = im->nz ;
it = nx-1 ; jt = ny-1 ; kt = nz-1 ;
/* compute CM of image */
mri_get_cmass_3D( im, &xcm, &ycm, &zcm ) ;
ic = (int)rint(xcm); jc = (int)rint(ycm); kc = (int)rint(zcm);
/* compute cliplevel in each octant about the CM */
val = 0.5f * THD_cliplevel( im,mfrac ) ;
ii = (int)rint(0.01*nx); if( ii < 1 ) ii = 1 ; /* 1% quadrant overlap */
jj = (int)rint(0.01*ny); if( jj < 1 ) jj = 1 ;
kk = (int)rint(0.01*nz); if( kk < 1 ) kk = 1 ;
icm = ic-ii ; icp = ic+ii ; if( icm < 0 ) icm = 0; if( icp > it ) icp = it ;
jcm = jc-jj ; jcp = jc+jj ; if( jcm < 0 ) jcm = 0; if( jcp > jt ) jcp = jt ;
kcm = kc-kk ; kcp = kc+kk ; if( kcm < 0 ) kcm = 0; if( kcp > kt ) kcp = kt ;
cv.clip_000 = THD_cliplevel_partial( im,mfrac, 0 ,icp, 0 ,jcp, 0 ,kcp );
cv.clip_100 = THD_cliplevel_partial( im,mfrac, icm,it , 0 ,jcp, 0 ,kcp );
cv.clip_010 = THD_cliplevel_partial( im,mfrac, 0 ,icp, jcm,jt , 0 ,kcp );
cv.clip_110 = THD_cliplevel_partial( im,mfrac, icm,it , jcm,jt , 0 ,kcp );
cv.clip_001 = THD_cliplevel_partial( im,mfrac, 0 ,icp, 0 ,jcp, kcm,kt );
cv.clip_101 = THD_cliplevel_partial( im,mfrac, icm,it , 0 ,jcp, kcm,kt );
cv.clip_011 = THD_cliplevel_partial( im,mfrac, 0 ,icp, jcm,jt , kcm,kt );
cv.clip_111 = THD_cliplevel_partial( im,mfrac, icm,it , jcm,jt , kcm,kt );
if( cv.clip_000 < val ) cv.clip_000 = val ; /* don't let */
if( cv.clip_100 < val ) cv.clip_100 = val ; /* clip levels */
if( cv.clip_010 < val ) cv.clip_010 = val ; /* get too */
if( cv.clip_110 < val ) cv.clip_110 = val ; /* small! */
if( cv.clip_001 < val ) cv.clip_001 = val ;
if( cv.clip_101 < val ) cv.clip_101 = val ;
if( cv.clip_011 < val ) cv.clip_011 = val ;
if( cv.clip_111 < val ) cv.clip_111 = val ;
clip_min = cv.clip_000 ;
clip_min = MIN(clip_min,cv.clip_100) ;
clip_min = MIN(clip_min,cv.clip_010) ;
clip_min = MIN(clip_min,cv.clip_110) ;
clip_min = MIN(clip_min,cv.clip_001) ;
clip_min = MIN(clip_min,cv.clip_101) ;
clip_min = MIN(clip_min,cv.clip_011) ;
clip_min = MIN(clip_min,cv.clip_111) ; cv.clip_min = clip_min ;
clip_max = cv.clip_000 ;
clip_max = MAX(clip_max,cv.clip_100) ;
clip_max = MAX(clip_max,cv.clip_010) ;
clip_max = MAX(clip_max,cv.clip_110) ;
clip_max = MAX(clip_max,cv.clip_001) ;
clip_max = MAX(clip_max,cv.clip_101) ;
clip_max = MAX(clip_max,cv.clip_011) ;
clip_max = MAX(clip_max,cv.clip_111) ; cv.clip_max = clip_max ;
/* (x0,y0,z0) = center of lowest octant
(x1,y1,z1) = center of highest octant */
cv.x0 = 0.5*ic ; cv.x1 = 0.5*(ic+it) ;
cv.y0 = 0.5*jc ; cv.y1 = 0.5*(jc+jt) ;
cv.z0 = 0.5*kc ; cv.z1 = 0.5*(kc+kt) ;
cv.dxi = (cv.x1 > cv.x0) ? 1.0/(cv.x1-cv.x0) : 0.0 ;
cv.dyi = (cv.y1 > cv.y0) ? 1.0/(cv.y1-cv.y0) : 0.0 ;
cv.dzi = (cv.z1 > cv.z0) ? 1.0/(cv.z1-cv.z0) : 0.0 ;
RETURN(cv) ;
}
/*--------------------------------------------------------------------------*/
/*! Return the cliplevel at any point,
tri-linearly interpolated between octant centers. [24 Oct 2006] */
static INLINE float pointclip( int ii, int jj, int kk , clipvec *cv )
{
float x1,y1,z1 , x0,y0,z0 , val ;
/* get relative position in box defined by octant centers */
x1 = (ii-cv->x0)*cv->dxi; if(x1 < 0.0f) x1=0.0f; else if(x1 > 1.0f) x1=1.0f;
y1 = (jj-cv->y0)*cv->dyi; if(y1 < 0.0f) y1=0.0f; else if(y1 > 1.0f) y1=1.0f;
z1 = (kk-cv->z0)*cv->dzi; if(z1 < 0.0f) z1=0.0f; else if(z1 > 1.0f) z1=1.0f;
x0 = 1.0f-x1 ; y0 = 1.0f-y1 ; z0 = 1.0f-z1 ;
val = cv->clip_000 * x0*y0*z0 + cv->clip_100 * x1*y0*z0
+ cv->clip_010 * x0*y1*z0 + cv->clip_110 * x1*y1*z0
+ cv->clip_001 * x0*y0*z1 + cv->clip_101 * x1*y0*z1
+ cv->clip_011 * x0*y1*z1 + cv->clip_111 * x1*y1*z1 ;
return val ;
}
/*--------------------------------------------------------------------------*/
/*! Return a cliplevel that varies gradually across the image.
At this time [24 Oct 2006], the clipping varies trilinearly
between octants. A fancier method would use higher order
polynomials, but I can't really be bothered with that now.
----------------------------------------------------------------------------*/
MRI_IMAGE * THD_cliplevel_gradual( MRI_IMAGE *im , float mfrac )
{
int ii,jj,kk,ijk , nx,ny,nz ;
MRI_IMAGE *cim ; float *car ;
clipvec bvec ;
ENTRY("THD_cliplevel_gradual") ;
if( im == NULL ) RETURN(NULL) ;
bvec = get_octant_clips( im , mfrac ) ;
if( bvec.clip_000 < 0.0 ) RETURN(NULL) ;
nx = im->nx ; ny = im->ny ; nz = im->nz ;
cim = mri_new_conforming( im , MRI_float ) ;
car = MRI_FLOAT_PTR(cim) ;
for( ijk=kk=0 ; kk < nz ; kk++ ){
for( jj=0 ; jj < ny ; jj++ ){
for( ii=0 ; ii < nx ; ii++,ijk++ ){
car[ijk] = pointclip( ii,jj,kk , &bvec ) ; /* cliplevel here */
}}}
RETURN(cim) ;
}
syntax highlighted by Code2HTML, v. 0.9.1