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    /* scan command-line args */
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    /* load data from dataset */
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 }