Doxygen Source Code Documentation
thd_cliplevel.c File Reference
#include "mrilib.h"Go to the source code of this file.
Functions | |
| float | THD_cliplevel (MRI_IMAGE *im, float mfrac) |
Function Documentation
|
||||||||||||
|
Definition at line 10 of file thd_cliplevel.c. References calloc, ENTRY, free, MRI_IMAGE::kind, MRI_BYTE_PTR, mri_free(), mri_maxabs(), MRI_SHORT_PTR, mri_to_short(), MRI_IMAGE::nvox, and RETURN. Referenced by main(), mri_automask_image(), mri_warp3D_align_setup(), THD_autobbox(), THD_autonudge(), THD_orient_guess(), and THD_outlier_count().
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 }
|
