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 "afni.h"
8
9 #ifndef ALLOW_PLUGINS
10 # error "Plugins not properly set up -- see machdep.h"
11 #endif
12
13 static char helpstring[] =
14 "The purpose of this plugin is to fit a Beta distribution to\n"
15 "an empirical histogram. From this, you may be able to infer\n"
16 "significance. The Beta distribution has two parameters (a,b).\n"
17 "If the noise in an FMRI time series is Gaussian-white, and\n"
18 "linear regression is used (as in FIM), then in theory you\n"
19 "should have a=M/2 and b=(N-M-L)/2, where M=number of signal\n"
20 "parameter in the regressors (1 or 2), L=number of ort parameters,\n"
21 "and N=number of time points used. This theory does not usually\n"
22 "reflect the truth very well. With this plugin, you can directly\n"
23 "estimate a and b, and use this to select a threshold.\n"
24 "\n"
25 "Source : Dataset specifies the input dataset\n"
26 " Brick specifies which sub-brick contains the\n"
27 " correlation coefficient (FIM) or R**2 value\n"
28 " (3dDeconvolve or 3dNLfim)\n"
29 " Square = Yes if the plugin should square the\n"
30 " brick values prior to using them\n"
31 " (use this for correlations)\n"
32 " No if the values from the brick should\n"
33 " be used as is (use this for R**2)\n"
34 "a Params : Range of values to search for the a parameter\n"
35 "b Params : Range of values to search for the b parameter\n"
36 " H last = top value to plot in the histogram\n"
37 " (leave 0 to let the plugin choose)\n"
38 "Misc : N ran = number of random (a,b) values to start\n"
39 " the search with\n"
40 " % cut = how much of the cumulative histogram to\n"
41 " use in the search for (a,b)\n"
42 " Hsqrt = No to plot a normal histogram\n"
43 " Yes to plot the y-axis as the square-root\n"
44 " of the normal histogram counts (this is\n"
45 " used to emphasize the smaller bins)\n"
46 "Mask : Used to select a sub-region from which the data\n"
47 " source values will be taken\n"
48 " Dataset = which dataset has the mask\n"
49 " Brick = which sub-brick has the mask\n"
50 "Range : Specifies the range of mask values to use\n"
51 "Extra : Used to specify (a,b) for an extra plot on top\n"
52 " of the empirical histogram and the best fit Beta\n"
53 " distribution. When I use this, I pick (a,b) to\n"
54 " be the theoretical white noise values.\n"
55 "\n"
56 "The output is a graph of the correlation (or R**2) histogram,\n"
57 "with the best fit Beta(a,b) distribution overlaid. The graph\n"
58 "is labeled with the (a,b) parameters. To pick significance\n"
59 "for FMRI, you might choose a cutoff value at the very tail\n"
60 "of the fitted Beta(a,b) distribution, assuming that it does\n"
61 "look like a good fit to the central part of the data.\n"
62 "\n"
63 "-- RWCox - June 2000\n"
64 ;
65
66 static char * BFIT_main( PLUGIN_interface * ) ;
67
68 #define NYESNO 2
69 static char * YESNO_strings[NYESNO] = { "No" , "Yes" } ;
70
71
72 DEFINE_PLUGIN_PROTOTYPE
73
PLUGIN_init(int ncall)74 PLUGIN_interface * PLUGIN_init(int ncall)
75 {
76 PLUGIN_interface * plint ;
77
78 if( ncall > 0 ) return NULL ; /* only one interface */
79 CHECK_IF_ALLOWED("HISTOGRAMBFIT","Histogram: BFit") ; /* 30 Sep 2016 */
80
81 /*-- set titles and call point --*/
82
83 plint = PLUTO_new_interface( "Histogram: BFit" ,
84 "Betafit Histogram" ,
85 helpstring ,
86 PLUGIN_CALL_VIA_MENU , BFIT_main ) ;
87
88 PLUTO_add_hint( plint , "Histogram: Betafit" ) ;
89
90 PLUTO_set_sequence( plint , "A:afniinfo:dsethistog" ) ;
91
92 /*-- first line of input --*/
93
94 PLUTO_add_option( plint , "Source" , "Source" , TRUE ) ;
95
96 PLUTO_add_dataset( plint ,
97 "Dataset" , /* label next to button */
98 ANAT_ALL_MASK , /* take any anat datasets */
99 FUNC_ALL_MASK , /* only allow fim funcs */
100 DIMEN_3D_MASK | /* need 3D+time datasets */
101 BRICK_ALLREAL_MASK /* need real-valued datasets */
102 ) ;
103 PLUTO_add_number( plint , "Brick" , 0,9999,0, 0,1 ) ;
104 PLUTO_add_string( plint , "Square" , NYESNO , YESNO_strings , 1 ) ;
105
106 /*-- second line of input --*/
107
108 PLUTO_add_option( plint , "a Params" , "Params" , TRUE ) ;
109 PLUTO_add_number( plint , "a bot" , 2,50 ,1 , 5 , 1 ) ;
110 PLUTO_add_number( plint , "a top" , 2,500,1 , 20 , 1 ) ;
111
112 PLUTO_add_option( plint , "b Params" , "Params" , TRUE ) ;
113 PLUTO_add_number( plint , "b bot" , 10,400 ,0 , 10 , 1 ) ;
114 PLUTO_add_number( plint , "b top" , 10,9999,0 , 200 , 1 ) ;
115 PLUTO_add_number( plint , "H last", 0,1000,-1 , 0,1 ) ;
116
117 PLUTO_add_option( plint , "Misc" , "Params" , TRUE ) ;
118 PLUTO_add_number( plint , "N ran" , 10,1000,-2 , 100 , 1 ) ;
119 PLUTO_add_number( plint , "% cut" , 20,90,0 , 70,1 ) ;
120 PLUTO_add_string( plint , "HSqrt" , NYESNO , YESNO_strings , 0 ) ;
121
122 /*-- (optional) line of input --*/
123
124 PLUTO_add_option( plint , "Mask" , "Mask" , FALSE ) ;
125 PLUTO_add_dataset( plint , "Dataset" ,
126 ANAT_ALL_MASK , FUNC_ALL_MASK ,
127 DIMEN_ALL_MASK | BRICK_ALLREAL_MASK ) ;
128 PLUTO_add_number( plint , "Brick" , 0,9999,0 , 0,1 ) ;
129
130 /*-- (optional) line of input --*/
131
132 PLUTO_add_option( plint , "Range" , "Range" , FALSE ) ;
133 PLUTO_add_number( plint , "Bottom" , -99999,99999, 1, 0,1 ) ;
134 PLUTO_add_number( plint , "Top" , -99999,99999,-1, 0,1 ) ;
135
136 /*-- (optional) line of input --*/
137
138 PLUTO_add_option( plint , "Extra" , "Extra" , FALSE ) ;
139 PLUTO_add_number( plint , "a" , 2,50,1 , 5 , 1 ) ;
140 PLUTO_add_number( plint , "b" , 10,999,0 , 200 , 1 ) ;
141
142 return plint ;
143 }
144
145 #include "betafit.c"
146
147 /*--------------------------------------------------------------------*/
148
BFIT_main(PLUGIN_interface * plint)149 static char * BFIT_main( PLUGIN_interface * plint )
150 {
151 MCW_idcode * idc ;
152 THD_3dim_dataset * input_dset , * mask_dset = NULL ;
153
154 BFIT_data * bfd ;
155 BFIT_result * bfr ;
156
157 int nvals,ival , nran,nvox , nbin , miv=0 , sqr,sqt ;
158 float abot,atop,bbot,btop,pcut , eps,eps1 , hlast ;
159 float *bval , *cval ;
160 double aa,bb,xc ;
161 double chq,ccc,cdf ;
162 int ihqbot,ihqtop ;
163
164 int mcount,mgood , ii , jj , ibot,itop ;
165 float mask_bot=666.0 , mask_top=-666.0 , hbot,htop,dbin ;
166 char buf[THD_MAX_NAME+128] , tbuf[THD_MAX_NAME+128] , * tag ;
167 int * hbin , * jbin,*kbin=NULL , *jist[2] ;
168 MRI_IMAGE * flim ;
169
170 double aext=-1.0,bext=-1.0 ;
171
172 /*--------------------------------------------------------------------*/
173 /*----- Check inputs from AFNI to see if they are reasonable-ish -----*/
174
175 if( plint == NULL )
176 return "************************\n"
177 "BFIT_main: NULL input\n"
178 "************************" ;
179
180 /*-- read 1st line --*/
181
182 PLUTO_next_option(plint) ;
183 idc = PLUTO_get_idcode(plint) ;
184 input_dset = PLUTO_find_dset(idc) ;
185 if( input_dset == NULL )
186 return "****************************\n"
187 "BFIT_main: bad input dataset\n"
188 "****************************" ;
189
190 nvox = DSET_NVOX(input_dset) ;
191 nvals = DSET_NVALS(input_dset) ;
192 ival = (int) PLUTO_get_number(plint) ;
193 if( ival < 0 || ival >= nvals )
194 return "**************************\n"
195 "BFIT_main: bad Brick index\n"
196 "**************************" ;
197
198 DSET_load(input_dset) ;
199 if( DSET_ARRAY(input_dset,0) == NULL )
200 return "*****************************\n"
201 "BFIT_main: can't load dataset\n"
202 "*****************************" ;
203
204 tag = PLUTO_get_string(plint) ;
205 sqr = PLUTO_string_index(tag,NYESNO,YESNO_strings) ;
206
207 /*-- read 2nd line --*/
208
209 PLUTO_next_option(plint) ;
210 abot = PLUTO_get_number(plint) ;
211 atop = PLUTO_get_number(plint) ;
212 if( atop <= abot )
213 return "*** atop <= abot! ***" ;
214
215 PLUTO_next_option(plint) ;
216 bbot = PLUTO_get_number(plint) ;
217 btop = PLUTO_get_number(plint) ;
218 if( atop <= abot )
219 return "*** btop <= bbot! ***" ;
220 hlast = PLUTO_get_number(plint) ;
221
222 PLUTO_next_option(plint) ;
223 nran = (int) PLUTO_get_number(plint) ;
224 pcut = PLUTO_get_number(plint) ;
225
226 tag = PLUTO_get_string(plint) ;
227 sqt = PLUTO_string_index(tag,NYESNO,YESNO_strings) ;
228
229 /*-- read optional lines --*/
230
231 while( (tag=PLUTO_get_optiontag(plint)) != NULL ){
232
233 /*-- Mask itself --*/
234
235 if( strcmp(tag,"Mask") == 0 ){
236
237 idc = PLUTO_get_idcode(plint) ;
238 mask_dset = PLUTO_find_dset(idc) ;
239
240 if( mask_dset == NULL ){
241 return "******************************\n"
242 "BFIT_main: bad mask dataset\n"
243 "******************************" ;
244 }
245
246 if( DSET_NVOX(mask_dset) != nvox ){
247 return "************************************************************\n"
248 "BFIT_main: mask input dataset doesn't match source dataset\n"
249 "************************************************************" ;
250 }
251
252 miv = (int) PLUTO_get_number(plint) ;
253 if( miv >= DSET_NVALS(mask_dset) || miv < 0 ){
254 return "****************************************************\n"
255 "BFIT_main: mask dataset sub-brick index is illegal\n"
256 "****************************************************" ;
257 }
258
259 DSET_load(mask_dset) ;
260 if( DSET_ARRAY(mask_dset,miv) == NULL ){
261 return "*************************************\n"
262 "BFIT_main: can't load mask dataset\n"
263 "*************************************" ;
264 }
265 continue ;
266 }
267
268 /*-- Mask range of values --*/
269
270 if( strcmp(tag,"Range") == 0 ){
271 if( mask_dset == NULL ){
272 return "******************************************\n"
273 "BFIT_main: Can't use Range without Mask\n"
274 "******************************************" ;
275 }
276
277 mask_bot = PLUTO_get_number(plint) ;
278 mask_top = PLUTO_get_number(plint) ;
279 continue ;
280 }
281
282 /*-- Extra plot --*/
283
284 if( strcmp(tag,"Extra") == 0 ){
285 aext = PLUTO_get_number(plint) ;
286 bext = PLUTO_get_number(plint) ;
287 continue ;
288 }
289 }
290
291 /*------------------------------------------------------*/
292 /*---------- At this point, the inputs are OK ----------*/
293
294 bfd = BFIT_prepare_dataset( input_dset , ival , sqr ,
295 mask_dset , miv , mask_bot , mask_top ) ;
296
297 if( bfd == NULL ) return "*** BFIT_prepare_dataset fails ***" ;
298
299 bfr = BFIT_compute( bfd ,
300 pcut , abot,atop , bbot,btop , nran,200 ) ;
301
302 if( bfr == NULL ){
303 BFIT_free_data( bfd ) ;
304 return "*** BFIT_compute fails! ***" ;
305 }
306
307 itop = bfr->itop ;
308 mgood = bfr->mgood ;
309
310 ibot = bfd->ibot ;
311 bval = bfd->bval ;
312 cval = bfd->cval ;
313 mcount = bfd->mcount ;
314
315 xc = bfr->xcut ;
316 aa = bfr->a ;
317 bb = bfr->b ;
318 eps = bfr->eps ;
319 eps1 = 1.0 - eps ;
320 if( eps1 > 1.0 ) eps1 = 1.0 ;
321 eps1 = (mcount-ibot) * eps1 ;
322
323 /*-- compute and plot histogram --*/
324
325 /* original data was already squared (e.g., R**2 values) */
326
327 if( !sqr ){
328 hbot = 0.0 ; htop = 1.0 ; nbin = 200 ;
329 if( bval[mcount-1] < 1.0 ) htop = bval[mcount-1] ;
330 dbin = (htop-hbot)/nbin ;
331
332 hbin = (int *) calloc((nbin+1),sizeof(int)) ; /* actual histogram */
333 jbin = (int *) calloc((nbin+1),sizeof(int)) ; /* theoretical fit */
334
335 for( ii=0 ; ii < nbin ; ii++ ){ /* beta fit */
336 jbin[ii] = (int)( eps1 * ( beta_t2p(hbot+ii*dbin,aa,bb)
337 -beta_t2p(hbot+ii*dbin+dbin,aa,bb) ) ) ;
338 }
339
340 jist[0] = jbin ;
341
342 flim = mri_new_vol_empty( mcount-ibot,1,1 , MRI_float ) ;
343 mri_fix_data_pointer( bval+ibot , flim ) ;
344 mri_histogram( flim , hbot,htop , TRUE , nbin,hbin ) ;
345
346 /* "extra" histogram (nominal values?) */
347
348 if( aext > 0.0 ){
349 kbin = (int *) calloc((nbin+1),sizeof(int)) ;
350 jist[1] = kbin ;
351 for( ii=0 ; ii < nbin ; ii++ ){ /* beta fit */
352 kbin[ii] = (int)( eps1 * ( beta_t2p(hbot+ii*dbin,aext,bext)
353 -beta_t2p(hbot+ii*dbin+dbin,aext,bext) ) ) ;
354 }
355 }
356
357 } else { /* original data was not squared (e.g., correlations) */
358
359 double hb,ht ;
360 htop = 1.0 ; nbin = 200 ;
361 if( bval[mcount-1] < 1.0 ) htop = sqrt(bval[mcount-1]) ;
362 hbot = -htop ;
363 dbin = (htop-hbot)/nbin ;
364
365 hbin = (int *) calloc((nbin+1),sizeof(int)) ; /* actual histogram */
366 jbin = (int *) calloc((nbin+1),sizeof(int)) ; /* theoretical fit */
367
368 for( ii=0 ; ii < nbin ; ii++ ){ /* beta fit */
369 hb = hbot+ii*dbin ; ht = hb+dbin ;
370 hb = hb*hb ; ht = ht*ht ;
371 if( hb > ht ){ double qq=hb ; hb=ht ; ht=qq ; }
372 jbin[ii] = (int)( 0.5*eps1 * ( beta_t2p(hb,aa,bb)
373 -beta_t2p(ht,aa,bb) ) ) ;
374 }
375
376 jist[0] = jbin ;
377
378 flim = mri_new_vol_empty( mcount-ibot,1,1 , MRI_float ) ;
379 mri_fix_data_pointer( cval+ibot , flim ) ;
380 mri_histogram( flim , hbot,htop , TRUE , nbin,hbin ) ;
381
382 /* nominal fit */
383
384 if( aext > 0.0 ){
385 kbin = (int *) calloc((nbin+1),sizeof(int)) ;
386 jist[1] = kbin ;
387 for( ii=0 ; ii < nbin ; ii++ ){ /* beta fit */
388 hb = hbot+ii*dbin ; ht = hb+dbin ;
389 hb = hb*hb ; ht = ht*ht ;
390 if( hb > ht ){ double qq=hb ; hb=ht ; ht=qq ; }
391 kbin[ii] = (int)( 0.5*eps1 * ( beta_t2p(hb,aext,bext)
392 -beta_t2p(ht,aext,bext) ) ) ;
393 }
394 }
395 }
396
397 sprintf(buf,"%s[%d] a=%.2f b=%.2f \\epsilon=%.2f %%=%.0f",
398 DSET_FILECODE(input_dset),ival,aa,bb,eps,pcut ) ;
399
400 ccc = bfr->q_chisq ;
401
402 /* blow up histogram details by sqrt-ing, if ordered */
403
404 if( sqt ){
405 for( ii=0 ; ii < nbin ; ii++ ){
406 hbin[ii] = (int) sqrt( (double)(100*hbin[ii]+0.5) ) ;
407 jbin[ii] = (int) sqrt( (double)(100*jbin[ii]+0.5) ) ;
408 if( kbin!=NULL )
409 kbin[ii] = (int) sqrt( (double)(100*kbin[ii]+0.5) ) ;
410 }
411 }
412
413 /* and plot */
414
415 sprintf(tbuf,"\\beta fit: cutoff=%.2f nvox=%d q(\\chi^2)=%8.2e",
416 (sqr)?sqrt(xc):xc , mgood , ccc ) ;
417 if( sqt ){
418 ii = strlen(tbuf) ;
419 sprintf( tbuf+ii , " \\surd ogram" ) ;
420 }
421
422 if( hlast > 0.0 ){
423 hbin[nbin-1] = jbin[nbin-1] = hlast ;
424 if( kbin != NULL ) kbin[nbin-1] = hlast ;
425 }
426
427 PLUTO_histoplot( nbin,hbot,htop,hbin ,
428 tbuf,NULL,buf , (kbin==NULL)?1:2 , jist ) ;
429
430 /* cleanup */
431
432 mri_clear_data_pointer(flim) ; mri_free(flim) ;
433 free(hbin) ; free(jbin) ; if( kbin != NULL ) free(kbin);
434
435 BFIT_free_data(bfd) ; BFIT_free_result(bfr) ;
436 return NULL ;
437 }
438