Doxygen Source Code Documentation

Main Page Alphabetical List Data Structures File List Data Fields Globals Search
thd_cliplevel.c

00001 #include "mrilib.h"
00002 
00003 /*--------------------------------------------------------------------------
00004    12 Aug 2001: compare with 3dClipLevel.c
00005    - compute a clipping level for an image, to eliminate non-brain voxels
00006    05 Nov 2001: increased size of hist array to nhist+1 (from nhist), to
00007                 store properly elements [0..nhist] (d'oh).
00008 ----------------------------------------------------------------------------*/
00009 
00010 float THD_cliplevel( MRI_IMAGE *im , float mfrac )
00011 {
00012    MRI_IMAGE *lim ;
00013    float fac , sfac=1.0 ;
00014    double dsum ;
00015    int nvox , *hist , ii,npos=0 , ncut,kk,ib , qq,nold ;
00016    short *sar ;
00017    byte  *bar ;
00018    int nhist , nneg=0 , nhalf ;
00019 
00020 ENTRY("THD_cliplevel") ;
00021    if( im == NULL ) RETURN(0.0) ;
00022 
00023    if( mfrac <= 0.0 || mfrac >= 0.99 ) mfrac = 0.50 ;
00024 
00025    /*-- allocate histogram --*/
00026 
00027    switch( im->kind ){
00028       case MRI_short: nhist = 32767 ; lim = im ; break ;
00029       case MRI_byte : nhist =   255 ; lim = im ; break ;
00030 
00031       default:
00032         fac = mri_maxabs(im) ; if( fac == 0.0 ) RETURN(0.0) ;
00033         sfac = 32767.0/fac ; nhist = 32767 ;
00034         lim = mri_to_short( sfac , im ) ;
00035       break ;
00036    }
00037 
00038    hist = (int *) calloc(sizeof(int),nhist+1) ;  /* 05 Nov 2001: +1 */
00039    nvox = im->nvox ;
00040 
00041    /*-- make histogram --*/
00042 
00043    dsum = 0.0 ;
00044    switch( lim->kind ){
00045       default: break ;
00046 
00047       case MRI_short:
00048          sar =  MRI_SHORT_PTR(lim) ;
00049          for( ii=0 ; ii < nvox ; ii++ ){
00050             if( sar[ii] > 0 && sar[ii] <= nhist ){
00051                hist[sar[ii]]++ ;
00052                dsum += (double)(sar[ii])*(double)(sar[ii]); npos++;
00053             } else if( sar[ii] < 0 )
00054               nneg++ ;
00055          }
00056       break ;
00057 
00058       case MRI_byte:                       /* there are no negative bytes */
00059          bar =  MRI_BYTE_PTR(lim) ;
00060          for( ii=0 ; ii < nvox ; ii++ ){
00061             if( bar[ii] > 0 ){
00062                hist[bar[ii]]++ ;
00063                dsum += (double)(bar[ii])*(double)(bar[ii]); npos++;
00064             }
00065          }
00066       break ;
00067    }
00068 
00069    if( npos <= 999 ){ free(hist); if(lim!=im)mri_free(lim); RETURN(0.0); }
00070 
00071    /*-- initialize cut position to include upper 65% of positive voxels --*/
00072 
00073    qq = 0.65 * npos ; ib = rint(0.5*sqrt(dsum/npos)) ;
00074    for( kk=0,ii=nhist-1 ; ii >= ib && kk < qq ; ii-- ) kk += hist[ii] ;
00075 
00076    /*-- algorithm --*/
00077 
00078    ncut = ii ; qq = 0 ;
00079    do{
00080       for( npos=0,ii=ncut ; ii < nhist ; ii++ ) npos += hist[ii] ; /* number >= cut */
00081       nhalf = npos/2 ;
00082       for( kk=0,ii=ncut ; ii < nhist && kk < nhalf ; ii++ )        /* find median */
00083          kk += hist[ii] ;
00084       nold = ncut ;
00085       ncut = mfrac * ii ;                                          /* new cut */
00086       qq++ ;
00087    } while( qq < 20 && ncut != nold ) ;
00088 
00089    free(hist) ; if( lim != im ) mri_free(lim) ;
00090 
00091    RETURN( (ncut/sfac) ) ;
00092 }
AFNI/NIfTI Server

Sections

Personal tools

Navigation

Doxygen Source Code Documentation

thd_cliplevel.c
