Doxygen Source Code Documentation
3dnoise.c File Reference
#include "mrilib.h"Go to the source code of this file.
Defines | |
| #define | NBIN 32768 |
| #define | DMU 0.5 |
Functions | |
| void | init_histo (void) |
| int | load_histo (THD_3dim_dataset *dset, int iv) |
| double | chfit (double mu) |
| int | main (int argc, char *argv[]) |
Variables | |
| int | histo [NBIN] |
Define Documentation
|
|
|
|
|
Definition at line 9 of file 3dnoise.c. Referenced by chfit(), and init_histo(). |
Function Documentation
|
|
Definition at line 40 of file 3dnoise.c. Referenced by main().
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 }
|
|
|
Definition at line 12 of file 3dnoise.c. Referenced by main().
|
|
||||||||||||
|
Definition at line 19 of file 3dnoise.c. References DSET_ARRAY, DSET_BRICK_TYPE, DSET_load, DSET_NVALS, DSET_NVOX, histo, and ISVALID_DSET. Referenced by main().
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 }
|
|
||||||||||||
|
---------- Adapted from 3dZeropad.c by RWCox - 08 Aug 2001 ----------* Definition at line 67 of file 3dnoise.c. References argc, chfit(), DATABLOCK_MEM_MALLOC, THD_3dim_dataset::dblk, DSET_ARRAY, DSET_NUM_TIMES, DSET_NVALS, DSET_NVOX, DSET_write, histo, init_histo(), load_histo(), snr, strtod(), THD_delete_3dim_dataset(), THD_force_malloc_type(), THD_open_one_dataset(), and tross_Append_History().
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]) ; /** NOT THD_open_dataset()! **/
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 { /* 24 Mar 1998 */
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 }
|
Variable Documentation
|
|
Definition at line 10 of file 3dnoise.c. Referenced by chfit(), init_histo(), load_histo(), and main(). |
