1 /*****************************************************************************
2 Major portions of this software are copyrighted by the Medical College
3 of Wisconsin, 1994-2000, and are released under the Gnu General Public
4 License, Version 2. See the file README.Copyright for details.
5 ******************************************************************************/
6
7 #include "mrilib.h"
8 #include "betafit.c"
9
10 /*-----------------------------------------------------------------------*/
11
main(int argc,char * argv[])12 int main( int argc , char * argv[] )
13 {
14 BFIT_data * bfd , * nfd ;
15 BFIT_result * bfr , * nfr ;
16
17 int nvals,ival , nvox , nbin , miv , sqr=0 ;
18 float pcut , eps,eps1 ;
19 float *bval , *cval ;
20 double aa,bb,xc,xth ;
21 double chq,ccc,cdf ;
22 int ihqbot,ihqtop ;
23
24 int mcount,mgood , ii , jj , ibot,itop ;
25
26 int narg=1 , nran=1000 ;
27 float abot= 0.5 , atop= 4.0 ;
28 float bbot=10.0 , btop=200.0 ;
29 float pbot=50.0 , ptop= 80.0 ;
30
31 int nboot=0 ;
32 double aboot,bboot,tboot , pthr=1.e-4 ;
33 float asig , bsig , tsig , abcor ;
34
35 THD_3dim_dataset * input_dset , * mask_dset=NULL ;
36 float mask_bot=666.0 , mask_top=-666.0 ;
37 byte * mmm=NULL ;
38
39 if( argc < 2 || strcmp(argv[1],"-help") == 0 ){
40 printf("Usage: 3dbetafit [options] dataset\n"
41 "Fits a beta distribution to the values in a brick.\n"
42 "\n"
43 "Options:\n"
44 " -arange abot atop = Sets the search range for parameter\n"
45 " 'a' to abot..atop.\n"
46 " [default is 0.5 .. 4.0]\n"
47 "\n"
48 " -brange bbot btop = Sets the search range for parameter\n"
49 " 'b' to bbot..btop\n"
50 " [default is 10 .. 200]\n"
51 "\n"
52 " -prange pbot ptop = Will evaluate for percent cutoffs\n"
53 " from pbot to ptop (steps of 1%%)\n"
54 " [default is 50 .. 80]\n"
55 "\n"
56 " -bootstrap N = Does N bootstrap evaluations to\n"
57 " compute variance of 'a' and 'b'\n"
58 " estimates [default is no bootstrap]\n"
59 "\n"
60 " -mask mset = A mask dataset to indicate which\n"
61 " voxels are to be used\n"
62 " -mrange b t = Use only mask values in range from\n"
63 " 'b' to 't' (inclusive)\n"
64 "\n"
65 " -sqr = Flag to square the data from the dataset\n"
66 " -pthr p = Sets p-value of cutoff for threshold evaluation\n"
67 " [default = 1.e-4]\n"
68 ) ;
69 PRINT_COMPILE_DATE ; exit(0) ;
70 }
71
72 /* scan command-line args */
73
74 while( narg < argc && argv[narg][0] == '-' ){
75
76 if( strcmp(argv[narg],"-pthr") == 0 ){
77 pthr = strtod(argv[++narg],NULL) ;
78 if( pthr <= 0.0 || pthr >= 1.0 ){
79 fprintf(stderr,"*** Illegal value after -pthr!\n");exit(1);
80 }
81 narg++ ; continue;
82 }
83
84 if( strcmp(argv[narg],"-sqr") == 0 ){
85 sqr = 1 ; narg++ ; continue;
86 }
87
88 if( strcmp(argv[narg],"-arange") == 0 ){
89 abot = strtod(argv[++narg],NULL) ;
90 atop = strtod(argv[++narg],NULL) ;
91 if( abot < 0.1 || abot > atop ){
92 fprintf(stderr,"*** Illegal value after -arange!\n");exit(1);
93 }
94 narg++ ; continue;
95 }
96
97 if( strcmp(argv[narg],"-brange") == 0 ){
98 bbot = strtod(argv[++narg],NULL) ;
99 btop = strtod(argv[++narg],NULL) ;
100 if( bbot < 0.1 || bbot > btop ){
101 fprintf(stderr,"*** Illegal value after -brange!\n");exit(1);
102 }
103 narg++ ; continue;
104 }
105
106 if( strcmp(argv[narg],"-prange") == 0 ){
107 pbot = strtod(argv[++narg],NULL) ;
108 ptop = strtod(argv[++narg],NULL) ;
109 if( pbot < 30.0 || pbot > ptop ){
110 fprintf(stderr,"*** Illegal value after -prange!\n");exit(1);
111 }
112 narg++ ; continue;
113 }
114
115 if( strcmp(argv[narg],"-bootstrap") == 0 ){
116 nboot = strtod(argv[++narg],NULL) ;
117 if( nboot < 10 ){
118 fprintf(stderr,"*** Illegal value after -bootstrap!\n");exit(1);
119 }
120 narg++ ; continue;
121 }
122
123 if( strncmp(argv[narg],"-mask",5) == 0 ){
124 if( mask_dset != NULL ){
125 fprintf(stderr,"*** Cannot have two -mask options!\n") ; exit(1) ;
126 }
127 if( narg+1 >= argc ){
128 fprintf(stderr,"*** -mask option requires a following argument!\n");
129 exit(1) ;
130 }
131 mask_dset = THD_open_dataset( argv[++narg] ) ;
132 if( mask_dset == NULL ){
133 fprintf(stderr,"*** Cannot open mask dataset!\n") ; exit(1) ;
134 }
135 if( DSET_BRICK_TYPE(mask_dset,0) == MRI_complex ){
136 fprintf(stderr,"*** Cannot deal with complex-valued mask dataset!\n");
137 exit(1) ;
138 }
139 narg++ ; continue ;
140 }
141
142 if( strncmp(argv[narg],"-mrange",5) == 0 ){
143 if( narg+2 >= argc ){
144 fprintf(stderr,"*** -mrange option requires 2 following arguments!\n")
145 ;
146 exit(1) ;
147 }
148 mask_bot = strtod( argv[++narg] , NULL ) ;
149 mask_top = strtod( argv[++narg] , NULL ) ;
150 if( mask_top < mask_top ){
151 fprintf(stderr,"*** -mrange inputs are illegal!\n") ; exit(1) ;
152 }
153 narg++ ; continue ;
154 }
155
156 fprintf(stderr,"*** Illegal option: %s\n",argv[narg]) ; exit(1) ;
157 }
158
159 if( narg >= argc ){
160 fprintf(stderr,"*** No dataset argument on command line!?\n");exit(1);
161 }
162
163 input_dset = THD_open_dataset( argv[narg] ) ;
164 if( input_dset == NULL ){
165 fprintf(stderr,"*** Can't open dataset %s\n",argv[narg]); exit(1);
166 }
167
168 nvox = DSET_NVOX(input_dset) ;
169
170 /* load data from dataset */
171
172 DSET_load(input_dset) ; CHECK_LOAD_ERROR(input_dset) ;
173
174 if( DSET_BRICK_STATCODE(input_dset,0) == FUNC_COR_TYPE ) sqr = 1 ;
175
176 bfd = BFIT_prepare_dataset( input_dset , 0 , sqr ,
177 mask_dset , 0 , mask_bot , mask_top ) ;
178
179 if( bfd == NULL ){
180 fprintf(stderr,"*** Couldn't prepare data from input dataset!\n");
181 exit(1) ;
182 }
183
184 DSET_delete(mask_dset) ; DSET_delete(input_dset) ;
185
186 for( pcut=pbot ; pcut <= ptop ; pcut += 1.0 ){
187 bfr = BFIT_compute( bfd ,
188 pcut , abot,atop , bbot,btop , nran,200 ) ;
189 if( bfr == NULL ){
190 fprintf(stderr,"*** Can't compute betafit at pcut=%f\n",pcut) ;
191 exit(1) ;
192 }
193
194 itop = bfr->itop ;
195 mgood = bfr->mgood ;
196
197 ibot = bfd->ibot ;
198 bval = bfd->bval ;
199 cval = bfd->cval ;
200 mcount = bfd->mcount ;
201
202 xc = bfr->xcut ;
203 aa = bfr->a ;
204 bb = bfr->b ;
205 eps = bfr->eps ;
206
207 xth = beta_p2t( pthr , aa,bb ) ;
208 if( sqr ) xth = sqrt(xth) ;
209 if( sqr ) xc = sqrt(xc ) ;
210
211 if( nboot > 0 ){
212 asig = bsig = tsig = abcor = 0.0 ;
213 for( ii=0 ; ii < nboot ; ii++ ){
214 nfd = BFIT_bootstrap_sample( bfd ) ;
215 nfr = BFIT_compute( nfd ,
216 pcut , abot,atop , bbot,btop , nran,200 ) ;
217 aboot = nfr->a ;
218 bboot = nfr->b ;
219 tboot = beta_p2t( pthr , aboot,bboot ) ;
220 if( sqr ) tboot = sqrt(tboot) ;
221 BFIT_free_result(nfr) ;
222 BFIT_free_data (nfd) ;
223 asig += SQR( aboot-aa ) ;
224 bsig += SQR( bboot-bb ) ;
225 tsig += SQR( tboot-xth) ;
226 abcor+= (aboot-aa)*(bboot-bb) ;
227 }
228 asig = sqrt(asig/nboot) ;
229 bsig = sqrt(bsig/nboot) ;
230 tsig = sqrt(tsig/nboot) ;
231 abcor = abcor/(nboot*asig*bsig) ;
232 }
233
234 if( nboot <= 0 ){
235 printf("%3.0f%%: a=%.3f b=%.3f eps=%.3f cutoff=%.3f qfit=%8.2e thr=%.3f\n",
236 pcut , aa,bb,eps , xc , bfr->q_chisq , xth ) ;
237 } else {
238 printf("%3.0f%%: a=%.3f[%.3f] b=%.3f[%.3f] [%.3f] eps=%.3f cutoff=%.3f qfit=%8.2e thr=%.3f[%.3f]\n",
239 pcut , aa,asig,bb,bsig,abcor,eps , xc , bfr->q_chisq , xth,tsig ) ;
240 }
241 fflush(stdout) ;
242
243 BFIT_free_result(bfr) ;
244 }
245 exit(0) ;
246 }
247