Doxygen Source Code Documentation
Main Page Alphabetical List Data Structures File List Data Fields Globals Search
3dnoise.c
Go to the documentation of this file.00001
00002
00003
00004
00005
00006
00007 #include "mrilib.h"
00008
00009 #define NBIN 32768
00010 static int histo[NBIN] ;
00011
00012 void init_histo(void)
00013 {
00014 int ii ;
00015 for( ii=0 ; ii < NBIN ; ii++ ) histo[ii] = 0 ;
00016 return ;
00017 }
00018
00019 int load_histo( THD_3dim_dataset * dset , int iv )
00020 {
00021 int ii , nvox ;
00022 MRI_IMAGE * bim ;
00023 short * bar ;
00024
00025 if( !ISVALID_DSET(dset) ) return 0 ;
00026 if( iv < 0 || iv >= DSET_NVALS(dset) ) return 0 ;
00027 if( DSET_BRICK_TYPE(dset,iv) != MRI_short ) return 0 ;
00028
00029 DSET_load(dset) ; bar = DSET_ARRAY(dset,iv) ;
00030 if( bar == NULL ) return 0 ;
00031
00032 nvox = DSET_NVOX(dset) ;
00033 for( ii=0 ; ii < nvox ; ii++ )
00034 if( bar[ii] < 0 ) return 0 ;
00035
00036 for( ii=0 ; ii < nvox ; ii++ ) histo[bar[ii]]++ ;
00037 return nvox ;
00038 }
00039
00040 double chfit( double mu )
00041 {
00042 int ntop = (2.5 * mu + 0.5) ;
00043 double ftop , ntot , ccc , eee ;
00044 int ii ;
00045
00046 if( ntop <= 1 || ntop >= NBIN/2 ) return -1.0 ;
00047
00048 ntot = 0 ;
00049 for( ii=0 ; ii <= ntop ; ii++ ) ntot += histo[ii] ;
00050
00051 ftop = 1.0 - exp(-(ntop/mu)*(ntop/mu)) ;
00052 ntot *= ftop ;
00053
00054 ccc = 0.0 ;
00055 for( ii=0 ; ii <= ntop ; ii++ ){
00056 eee = ntot * ( exp(-(ii/mu)*(ii/mu)) -
00057 exp(-((ii+1)/mu)*((ii+1)/mu)) ) ;
00058
00059 if( eee < 1.0 ) eee = 1.0 ;
00060
00061 ccc += (histo[ii] - eee) * (histo[ii] - eee) / eee ;
00062 }
00063
00064 ccc /= (ntop+1) ; return ccc ;
00065 }
00066
00067 int main( int argc , char * argv[] )
00068 {
00069 THD_3dim_dataset * dset ;
00070 double mu , ccc , mbest,cbest , perc , snr=2.5 , nlxx ;
00071 int ii , narg=1 , blast=0 , iv , ncut , nnn , nvox , nl=0 ;
00072 int vmax ;
00073 short * bar ;
00074
00075 if( argc < 2 || strcmp(argv[1],"-help") == 0 ){
00076 printf("Usage: 3dnoise [-blast] [-snr fac] [-nl x ] datasets ...\n"
00077 "Estimates noise level in 3D datasets, and optionally\n"
00078 "set voxels below the noise threshold to zero.\n"
00079 "This only works on datasets that are stored as shorts,\n"
00080 "and whose elements are all nonnegative.\n"
00081 " -blast = Set values at or below the cutoff to zero.\n"
00082 " In 3D+time datasets, a spatial location\n"
00083 " is set to zero only if a majority of time\n"
00084 " points fall below the cutoff; in that case\n"
00085 " all the values at that location are zeroed.\n"
00086 " -snr fac = Set cutoff to 'fac' times the estimated\n"
00087 " noise level. Default fac = 2.5. What to\n"
00088 " use for this depends strongly on your MRI\n"
00089 " system -- I often use 5, but our true SNR\n"
00090 " is about 100 for EPI.\n"
00091 " -nl x = Set the noise level to 'x', skipping the\n"
00092 " estimation procedure. Also sets fac=1.0.\n"
00093 " You can use program 3dClipLevel to get an\n"
00094 " estimate of a value for 'x'.\n"
00095 "Author -- RW Cox\n"
00096 ) ;
00097 exit(0) ;
00098 }
00099
00100 while( narg < argc && argv[narg][0] == '-' ){
00101
00102 if( strcmp(argv[narg],"-blast") == 0 ){
00103 blast = 1 ; narg++ ; continue ;
00104 }
00105
00106 if( strcmp(argv[narg],"-snr") == 0 ){
00107 narg++ ;
00108 snr = strtod( argv[narg] , NULL ) ;
00109 if( snr <= 0.0 ){fprintf(stderr,"Illegal snr value!\n");exit(1);}
00110 narg++ ; continue ;
00111 }
00112
00113 if( strcmp(argv[narg],"-nl") == 0 ){
00114 narg++ ; nl = 1 ;
00115 nlxx = strtod( argv[narg] , NULL ) ; snr = 1.0 ;
00116 if( nlxx <= 0.0 ){fprintf(stderr,"Illegal nl value!\n");exit(1);}
00117 narg++ ; continue ;
00118 }
00119
00120 fprintf(stderr,"Unknown option: %s\n",argv[narg]) ; exit(1) ;
00121 }
00122 if( narg >= argc ){fprintf(stderr,"No datasets?\n");exit(1);}
00123
00124 for( ; narg < argc ; narg++ ){
00125 dset = THD_open_one_dataset(argv[narg]) ;
00126 if( dset == NULL ){
00127 printf("*** Can't open %s\n",argv[narg]) ;
00128 continue ;
00129 }
00130
00131 fprintf(stderr,"%s",argv[narg]) ; fflush(stdout) ;
00132
00133 init_histo() ;
00134 if( blast )
00135 THD_force_malloc_type( dset->dblk , DATABLOCK_MEM_MALLOC ) ;
00136 nvox = DSET_NVOX(dset) ;
00137
00138 vmax = 0 ;
00139 for( iv=0 ; iv < DSET_NVALS(dset) ; iv++ ){
00140 ii = load_histo(dset,iv) ;
00141 if( ii <= 0 ){ fprintf(stderr,": Can't load data, or illegal data!\n"); break; }
00142 bar = DSET_ARRAY(dset,iv) ;
00143 for( ii=0 ; ii < nvox ; ii++ ) if( vmax < bar[ii] ) vmax = bar[ii] ;
00144 }
00145 if( iv < DSET_NVALS(dset) ) continue ;
00146 if( vmax < 40 ){ fprintf(stderr,": Didn't fit noise model!\n"); continue; }
00147
00148 printf(":") ; fflush(stdout) ;
00149
00150 #define DMU 0.5
00151 if( !nl ){
00152 mu = 1.0 ; mbest = mu ; cbest = chfit(mu) ; ii = 0 ;
00153 do {
00154 mu += DMU ; ccc = chfit(mu) ;
00155 if( ccc > cbest ) break ;
00156 cbest = ccc ; mbest = mu ; ii++ ;
00157 if( mu > 0.05 * vmax ){ ii=0 ; break; }
00158 } while( 1 ) ;
00159 if( ii <= 0 ){ fprintf(stderr," Didn't fit noise model!\n"); continue; }
00160 } else {
00161 mbest = nlxx ;
00162 }
00163
00164 ncut = (int) (snr * mbest) ;
00165 nnn = 0 ;
00166 for( ii=0 ; ii <= ncut ; ii++ ) nnn += histo[ii] ;
00167 perc = (100.0*nnn) / (double)(DSET_NVOX(dset)*DSET_NVALS(dset)) ;
00168 fprintf(stderr," Cutoff=%d Count=%d [%4.1f%%]",ncut,nnn,perc) ;
00169
00170 if( blast && nnn > 0 ){
00171 int nkilled ;
00172 if( DSET_NUM_TIMES(dset) <= 1 ){
00173 for( iv=0 ; iv < DSET_NVALS(dset) ; iv++ ){
00174 bar = DSET_ARRAY(dset,iv) ;
00175 for( ii=0 ; ii < nvox ; ii++ ) if( bar[ii] <= ncut ) bar[ii]=0;
00176 }
00177 nkilled = nnn ;
00178 } else {
00179 int nvh = DSET_NVALS(dset)/2 , qq ;
00180 for( ii=0,nkilled=0 ; ii < nvox ; ii++ ){
00181 for( iv=0,qq=0 ; iv < DSET_NVALS(dset) ; iv++ ){
00182 bar = DSET_ARRAY(dset,iv) ; if( bar[ii] <= ncut ) qq++ ;
00183 }
00184 if( qq >= nvh ){
00185 for( iv=0 ; iv < DSET_NVALS(dset) ; iv++ ){
00186 bar = DSET_ARRAY(dset,iv) ; bar[ii]=0;
00187 }
00188 nkilled += DSET_NVALS(dset) ;
00189 }
00190 }
00191 }
00192 fprintf(stderr,"--blasted %d voxels\n",nkilled) ;
00193 { char str[128] ;
00194 sprintf(str,"3dnoise -- blasted %d voxels",nkilled) ;
00195 tross_Append_History( dset , str ) ;
00196 }
00197 DSET_write(dset) ;
00198 }
00199
00200 THD_delete_3dim_dataset( dset , False ) ;
00201 }
00202
00203 exit(0) ;
00204 }