00001
00002
00003
00004
00005
00006
00007 #include "mrilib.h"
00008 #include "betafit.c"
00009
00010
00011
00012 int main( int argc , char * argv[] )
00013 {
00014 BFIT_data * bfd , * nfd ;
00015 BFIT_result * bfr , * nfr ;
00016
00017 int nvals,ival , nvox , nbin , miv , sqr=0 ;
00018 float pcut , eps,eps1 ;
00019 float *bval , *cval ;
00020 double aa,bb,xc,xth ;
00021 double chq,ccc,cdf ;
00022 int ihqbot,ihqtop ;
00023
00024 int mcount,mgood , ii , jj , ibot,itop ;
00025
00026 int narg=1 , nran=1000 ;
00027 float abot= 0.5 , atop= 4.0 ;
00028 float bbot=10.0 , btop=200.0 ;
00029 float pbot=50.0 , ptop= 80.0 ;
00030
00031 int nboot=0 ;
00032 double aboot,bboot,tboot , pthr=1.e-4 ;
00033 float asig , bsig , tsig , abcor ;
00034
00035 THD_3dim_dataset * input_dset , * mask_dset=NULL ;
00036 float mask_bot=666.0 , mask_top=-666.0 ;
00037 byte * mmm=NULL ;
00038
00039 if( argc < 2 || strcmp(argv[1],"-help") == 0 ){
00040 printf("Usage: 3dbetafit [options] dataset\n"
00041 "Fits a beta distribution to the values in a brick.\n"
00042 "\n"
00043 "Options:\n"
00044 " -arange abot atop = Sets the search range for parameter\n"
00045 " 'a' to abot..atop.\n"
00046 " [default is 0.5 .. 4.0]\n"
00047 "\n"
00048 " -brange bbot btop = Sets the search range for parameter\n"
00049 " 'b' to bbot..btop\n"
00050 " [default is 10 .. 200]\n"
00051 "\n"
00052 " -prange pbot ptop = Will evaluate for percent cutoffs\n"
00053 " from pbot to ptop (steps of 1%%)\n"
00054 " [default is 50 .. 80]\n"
00055 "\n"
00056 " -bootstrap N = Does N bootstrap evaluations to\n"
00057 " compute variance of 'a' and 'b'\n"
00058 " estimates [default is no bootstrap]\n"
00059 "\n"
00060 " -mask mset = A mask dataset to indicate which\n"
00061 " voxels are to be used\n"
00062 " -mrange b t = Use only mask values in range from\n"
00063 " 'b' to 't' (inclusive)\n"
00064 "\n"
00065 " -sqr = Flag to square the data from the dataset\n"
00066 " -pthr p = Sets p-value of cutoff for threshold evaluation\n"
00067 " [default = 1.e-4]\n"
00068 ) ;
00069 exit(0) ;
00070 }
00071
00072
00073
00074 while( narg < argc && argv[narg][0] == '-' ){
00075
00076 if( strcmp(argv[narg],"-pthr") == 0 ){
00077 pthr = strtod(argv[++narg],NULL) ;
00078 if( pthr <= 0.0 || pthr >= 1.0 ){
00079 fprintf(stderr,"*** Illegal value after -pthr!\n");exit(1);
00080 }
00081 narg++ ; continue;
00082 }
00083
00084 if( strcmp(argv[narg],"-sqr") == 0 ){
00085 sqr = 1 ; narg++ ; continue;
00086 }
00087
00088 if( strcmp(argv[narg],"-arange") == 0 ){
00089 abot = strtod(argv[++narg],NULL) ;
00090 atop = strtod(argv[++narg],NULL) ;
00091 if( abot < 0.1 || abot > atop ){
00092 fprintf(stderr,"*** Illegal value after -arange!\n");exit(1);
00093 }
00094 narg++ ; continue;
00095 }
00096
00097 if( strcmp(argv[narg],"-brange") == 0 ){
00098 bbot = strtod(argv[++narg],NULL) ;
00099 btop = strtod(argv[++narg],NULL) ;
00100 if( bbot < 0.1 || bbot > btop ){
00101 fprintf(stderr,"*** Illegal value after -brange!\n");exit(1);
00102 }
00103 narg++ ; continue;
00104 }
00105
00106 if( strcmp(argv[narg],"-prange") == 0 ){
00107 pbot = strtod(argv[++narg],NULL) ;
00108 ptop = strtod(argv[++narg],NULL) ;
00109 if( pbot < 30.0 || pbot > ptop ){
00110 fprintf(stderr,"*** Illegal value after -prange!\n");exit(1);
00111 }
00112 narg++ ; continue;
00113 }
00114
00115 if( strcmp(argv[narg],"-bootstrap") == 0 ){
00116 nboot = strtod(argv[++narg],NULL) ;
00117 if( nboot < 10 ){
00118 fprintf(stderr,"*** Illegal value after -bootstrap!\n");exit(1);
00119 }
00120 narg++ ; continue;
00121 }
00122
00123 if( strncmp(argv[narg],"-mask",5) == 0 ){
00124 if( mask_dset != NULL ){
00125 fprintf(stderr,"*** Cannot have two -mask options!\n") ; exit(1) ;
00126 }
00127 if( narg+1 >= argc ){
00128 fprintf(stderr,"*** -mask option requires a following argument!\n");
00129 exit(1) ;
00130 }
00131 mask_dset = THD_open_dataset( argv[++narg] ) ;
00132 if( mask_dset == NULL ){
00133 fprintf(stderr,"*** Cannot open mask dataset!\n") ; exit(1) ;
00134 }
00135 if( DSET_BRICK_TYPE(mask_dset,0) == MRI_complex ){
00136 fprintf(stderr,"*** Cannot deal with complex-valued mask dataset!\n");
00137 exit(1) ;
00138 }
00139 narg++ ; continue ;
00140 }
00141
00142 if( strncmp(argv[narg],"-mrange",5) == 0 ){
00143 if( narg+2 >= argc ){
00144 fprintf(stderr,"*** -mrange option requires 2 following arguments!\n")
00145 ;
00146 exit(1) ;
00147 }
00148 mask_bot = strtod( argv[++narg] , NULL ) ;
00149 mask_top = strtod( argv[++narg] , NULL ) ;
00150 if( mask_top < mask_top ){
00151 fprintf(stderr,"*** -mrange inputs are illegal!\n") ; exit(1) ;
00152 }
00153 narg++ ; continue ;
00154 }
00155
00156 fprintf(stderr,"*** Illegal option: %s\n",argv[narg]) ; exit(1) ;
00157 }
00158
00159 if( narg >= argc ){
00160 fprintf(stderr,"*** No dataset argument on command line!?\n");exit(1);
00161 }
00162
00163 input_dset = THD_open_dataset( argv[narg] ) ;
00164 if( input_dset == NULL ){
00165 fprintf(stderr,"*** Can't open dataset %s\n",argv[narg]); exit(1);
00166 }
00167
00168 nvox = DSET_NVOX(input_dset) ;
00169
00170
00171
00172 DSET_load(input_dset) ;
00173 if( !DSET_LOADED(input_dset) ){
00174 fprintf(stderr,"*** Couldn't load dataset brick!\n");exit(1);
00175 }
00176
00177 if( DSET_BRICK_STATCODE(input_dset,0) == FUNC_COR_TYPE ) sqr = 1 ;
00178
00179 bfd = BFIT_prepare_dataset( input_dset , 0 , sqr ,
00180 mask_dset , 0 , mask_bot , mask_top ) ;
00181
00182 if( bfd == NULL ){
00183 fprintf(stderr,"*** Couldn't prepare data from input dataset!\n");
00184 exit(1) ;
00185 }
00186
00187 DSET_delete(mask_dset) ; DSET_delete(input_dset) ;
00188
00189 for( pcut=pbot ; pcut <= ptop ; pcut += 1.0 ){
00190 bfr = BFIT_compute( bfd ,
00191 pcut , abot,atop , bbot,btop , nran,200 ) ;
00192 if( bfr == NULL ){
00193 fprintf(stderr,"*** Can't compute betafit at pcut=%f\n",pcut) ;
00194 exit(1) ;
00195 }
00196
00197 itop = bfr->itop ;
00198 mgood = bfr->mgood ;
00199
00200 ibot = bfd->ibot ;
00201 bval = bfd->bval ;
00202 cval = bfd->cval ;
00203 mcount = bfd->mcount ;
00204
00205 xc = bfr->xcut ;
00206 aa = bfr->a ;
00207 bb = bfr->b ;
00208 eps = bfr->eps ;
00209
00210 xth = beta_p2t( pthr , aa,bb ) ;
00211 if( sqr ) xth = sqrt(xth) ;
00212 if( sqr ) xc = sqrt(xc ) ;
00213
00214 if( nboot > 0 ){
00215 asig = bsig = tsig = abcor = 0.0 ;
00216 for( ii=0 ; ii < nboot ; ii++ ){
00217 nfd = BFIT_bootstrap_sample( bfd ) ;
00218 nfr = BFIT_compute( nfd ,
00219 pcut , abot,atop , bbot,btop , nran,200 ) ;
00220 aboot = nfr->a ;
00221 bboot = nfr->b ;
00222 tboot = beta_p2t( pthr , aboot,bboot ) ;
00223 if( sqr ) tboot = sqrt(tboot) ;
00224 BFIT_free_result(nfr) ;
00225 BFIT_free_data (nfd) ;
00226 asig += SQR( aboot-aa ) ;
00227 bsig += SQR( bboot-bb ) ;
00228 tsig += SQR( tboot-xth) ;
00229 abcor+= (aboot-aa)*(bboot-bb) ;
00230 }
00231 asig = sqrt(asig/nboot) ;
00232 bsig = sqrt(bsig/nboot) ;
00233 tsig = sqrt(tsig/nboot) ;
00234 abcor = abcor/(nboot*asig*bsig) ;
00235 }
00236
00237 if( nboot <= 0 ){
00238 printf("%3.0f%%: a=%.3f b=%.3f eps=%.3f cutoff=%.3f qfit=%8.2e thr=%.3f\n",
00239 pcut , aa,bb,eps , xc , bfr->q_chisq , xth ) ;
00240 } else {
00241 printf("%3.0f%%: a=%.3f[%.3f] b=%.3f[%.3f] [%.3f] eps=%.3f cutoff=%.3f qfit=%8.2e thr=%.3f[%.3f]\n",
00242 pcut , aa,asig,bb,bsig,abcor,eps , xc , bfr->q_chisq , xth,tsig ) ;
00243 }
00244 fflush(stdout) ;
00245
00246 BFIT_free_result(bfr) ;
00247 }
00248 exit(0) ;
00249 }
