Doxygen Source Code Documentation
Main Page Alphabetical List Data Structures File List Data Fields Globals Search
3dClipLevel.c
Go to the documentation of this file.00001 #include "mrilib.h"
00002
00003 int main( int argc , char * argv[] )
00004 {
00005 float mfrac=0.50 , fac ;
00006 double dsum ;
00007 int nvox , iarg=1 , *hist , ii,npos=0 , ncut,nmed,kk,ib , qq,nold ;
00008 THD_3dim_dataset * dset ;
00009 short * sar ;
00010 byte * bar ;
00011 int nhist , nneg=0 , verb=0 , nhalf , iv,nvals ;
00012
00013 if( argc < 2 || strcmp(argv[1],"-help") == 0 ){
00014 printf("Usage: 3dClipLevel [options] dataset\n"
00015 "Estimates the value at which to clip the anatomical dataset so\n"
00016 "that background regions are set to zero.\n"
00017 "Method:\n"
00018 " Find the median of all positive values >= clip value.\n"
00019 " Set the clip value to 0.50 of this median.\n"
00020 " Repeat until the clip value doesn't change.\n"
00021 "Options:\n"
00022 " -mfrac ff = Use the number ff instead of 0.50 in the algorithm.\n"
00023 " -verb = The clip value is always printed to stdout. If\n"
00024 " this option is used to select verbose output,\n"
00025 " progress reports are printed to stderr as well.\n"
00026 "\n"
00027 "N.B.: This program only works with byte- and short-valued\n"
00028 " datasets, and prints a warning message if any input\n"
00029 " voxels are negative. If the dataset has more than one\n"
00030 " sub-brick, all sub-bricks are used to build the histogram.\n"
00031 "N.B.: Use at your own risk! You might want to use the AFNI Histogram\n"
00032 " plugin to see if the results are reasonable. This program is\n"
00033 " likely to produce bad results on images gathered with local\n"
00034 " RF coils, or with pulse sequences with unusual contrasts.\n"
00035 "\n"
00036 "A csh command line for the truly adventurous:\n"
00037 " afni -dset \"v1:time+orig<`3dClipLevel 'v1:time+orig[4]'` .. 10000>\"\n"
00038 "(the dataset is from the 'sample96.tgz' data samples). Can you\n"
00039 "figure out what this does?\n"
00040 "(Hint: each type of quote \"'` means something different to csh.)\n"
00041 ) ;
00042 exit(0) ;
00043 }
00044
00045 mainENTRY("3dClipLevel main"); machdep(); AFNI_logger("3dCliplevel",argc,argv);
00046
00047
00048
00049 while( iarg < argc && argv[iarg][0] == '-' ){
00050
00051 if( strcmp(argv[iarg],"-mfrac") == 0 ){
00052 mfrac = strtod( argv[++iarg] , NULL ) ;
00053 if( mfrac <= 0.0 ){ fprintf(stderr,"** ILLEGAL -mfrac\n");exit(1);}
00054 if( mfrac >= 1.0 ) mfrac *= 0.01 ;
00055 iarg++ ; continue ;
00056 }
00057
00058 if( strncmp(argv[iarg],"-verb",5) == 0 ){
00059 verb++ ; iarg++ ; continue ;
00060 }
00061
00062 fprintf(stderr,"** ILLEGAL option: %s\n",argv[iarg]) ; exit(1) ;
00063 }
00064
00065
00066
00067 dset = THD_open_dataset(argv[iarg]) ;
00068 if( !ISVALID_DSET(dset) ){ fprintf(stderr,"** CAN'T open dataset\n");exit(1); }
00069 if( DSET_BRICK_TYPE(dset,0) != MRI_short && DSET_BRICK_TYPE(dset,0) != MRI_byte ){
00070 fprintf(stderr,"** ILLEGAL dataset type\n");exit(1);
00071 }
00072 DSET_load(dset) ;
00073 if( !DSET_LOADED(dset) ){ fprintf(stderr,"** CAN'T load dataset\n");exit(1); }
00074
00075 nvals = DSET_NVALS(dset) ;
00076
00077
00078
00079 switch( DSET_BRICK_TYPE(dset,0) ){
00080 case MRI_short: nhist = 32767 ; break ;
00081 case MRI_byte : nhist = 255 ; break ;
00082 }
00083
00084 hist = (int *) calloc(sizeof(int),nhist+1) ;
00085 nvox = DSET_NVOX(dset) ;
00086
00087
00088
00089 dsum = 0.0 ;
00090 for( iv=0 ; iv < nvals ; iv++ ){
00091 switch( DSET_BRICK_TYPE(dset,0) ){
00092 case MRI_short:
00093 sar = DSET_ARRAY(dset,iv) ;
00094 for( ii=0 ; ii < nvox ; ii++ ){
00095 if( sar[ii] > 0 && sar[ii] <= nhist ){
00096 hist[sar[ii]]++ ; dsum += (double)(sar[ii])*(double)(sar[ii]) ; npos++ ;
00097 } else if( sar[ii] < 0 )
00098 nneg++ ;
00099 }
00100 break ;
00101
00102 case MRI_byte:
00103 bar = DSET_ARRAY(dset,iv) ;
00104 for( ii=0 ; ii < nvox ; ii++ ){
00105 if( bar[ii] > 0 ){
00106 hist[bar[ii]]++ ; dsum += (double)(bar[ii])*(double)(bar[ii]) ; npos++ ;
00107 }
00108 }
00109 break ;
00110 }
00111 }
00112
00113 DSET_unload(dset) ;
00114
00115 if( verb ) fprintf(stderr,"++ %d positive voxels\n",npos) ;
00116
00117 if( npos <= 999 ){
00118 fprintf(stderr,"** TOO few positive voxel (%d)!\n",npos) ; exit(1) ;
00119 }
00120 if( nneg > 0 ){
00121 fprintf(stderr,"++ WARNING: ignored %d negative voxels!\n",nneg) ;
00122 }
00123
00124
00125
00126 qq = 0.65 * npos ; ib = rint(0.5*sqrt(dsum/npos)) ;
00127 for( kk=0,ii=nhist-1 ; ii >= ib && kk < qq ; ii-- ) kk += hist[ii] ;
00128
00129
00130
00131 ncut = ii ; qq = 0 ;
00132 do{
00133 for( npos=0,ii=ncut ; ii < nhist ; ii++ ) npos += hist[ii] ;
00134 nhalf = npos/2 ;
00135 for( kk=0,ii=ncut ; ii < nhist && kk < nhalf ; ii++ )
00136 kk += hist[ii] ;
00137 if( verb )
00138 fprintf(stderr,"++ Clip=%d Median=%d Number>=Clip=%d\n",ncut,ii,npos) ;
00139 nold = ncut ;
00140 ncut = mfrac * ii ;
00141 qq++ ;
00142 } while( qq < 20 && ncut != nold ) ;
00143
00144 fac = DSET_BRICK_FACTOR(dset,0) ;
00145 if( fac > 0.0 )
00146 printf("%g\n",ncut*fac) ;
00147 else
00148 printf("%d\n",ncut) ;
00149
00150 exit(0) ;
00151 }