00001
00002
00003
00004
00005
00006
00007 #include "afni.h"
00008
00009 #ifndef ALLOW_PLUGINS
00010 # error "Plugins not properly set up -- see machdep.h"
00011 #endif
00012
00013 static char helpstring[] =
00014 "The purpose of this plugin is to fit a Beta distribution to\n"
00015 "an empirical histogram. From this, you may be able to infer\n"
00016 "significance. The Beta distribution has two parameters (a,b).\n"
00017 "If the noise in an FMRI time series is Gaussian-white, and\n"
00018 "linear regression is used (as in FIM), then in theory you\n"
00019 "should have a=M/2 and b=(N-M-L)/2, where M=number of signal\n"
00020 "parameter in the regressors (1 or 2), L=number of ort parameters,\n"
00021 "and N=number of time points used. This theory does not usually\n"
00022 "reflect the truth very well. With this plugin, you can directly\n"
00023 "estimate a and b, and use this to select a threshold.\n"
00024 "\n"
00025 "Source : Dataset specifies the input dataset\n"
00026 " Brick specifies which sub-brick contains the\n"
00027 " correlation coefficient (FIM) or R**2 value\n"
00028 " (3dDeconvolve or 3dNLfim)\n"
00029 " Square = Yes if the plugin should square the\n"
00030 " brick values prior to using them\n"
00031 " (use this for correlations)\n"
00032 " No if the values from the brick should\n"
00033 " be used as is (use this for R**2)\n"
00034 "a Params : Range of values to search for the a parameter\n"
00035 "b Params : Range of values to search for the b parameter\n"
00036 " H last = top value to plot in the histogram\n"
00037 " (leave 0 to let the plugin choose)\n"
00038 "Misc : N ran = number of random (a,b) values to start\n"
00039 " the search with\n"
00040 " % cut = how much of the cumulative histogram to\n"
00041 " use in the search for (a,b)\n"
00042 " Hsqrt = No to plot a normal histogram\n"
00043 " Yes to plot the y-axis as the square-root\n"
00044 " of the normal histogram counts (this is\n"
00045 " used to emphasize the smaller bins)\n"
00046 "Mask : Used to select a sub-region from which the data\n"
00047 " source values will be taken\n"
00048 " Dataset = which dataset has the mask\n"
00049 " Brick = which sub-brick has the mask\n"
00050 "Range : Specifies the range of mask values to use\n"
00051 "Extra : Used to specify (a,b) for an extra plot on top\n"
00052 " of the empirical histogram and the best fit Beta\n"
00053 " distribution. When I use this, I pick (a,b) to\n"
00054 " be the theoretical white noise values.\n"
00055 "\n"
00056 "The output is a graph of the correlation (or R**2) histogram,\n"
00057 "with the best fit Beta(a,b) distribution overlaid. The graph\n"
00058 "is labeled with the (a,b) parameters. To pick significance\n"
00059 "for FMRI, you might choose a cutoff value at the very tail\n"
00060 "of the fitted Beta(a,b) distribution, assuming that it does\n"
00061 "look like a good fit to the central part of the data.\n"
00062 "\n"
00063 "-- RWCox - June 2000\n"
00064 ;
00065
00066 static char * BFIT_main( PLUGIN_interface * ) ;
00067
00068 #define NYESNO 2
00069 static char * YESNO_strings[NYESNO] = { "No" , "Yes" } ;
00070
00071
00072 DEFINE_PLUGIN_PROTOTYPE
00073
00074 PLUGIN_interface * PLUGIN_init(int ncall)
00075 {
00076 PLUGIN_interface * plint ;
00077
00078 if( ncall > 0 ) return NULL ;
00079
00080
00081
00082 plint = PLUTO_new_interface( "Histogram: BFit" ,
00083 "Betafit Histogram" ,
00084 helpstring ,
00085 PLUGIN_CALL_VIA_MENU , BFIT_main ) ;
00086
00087 PLUTO_add_hint( plint , "Histogram: Betafit" ) ;
00088
00089 PLUTO_set_sequence( plint , "A:afniinfo:dsethistog" ) ;
00090
00091
00092
00093 PLUTO_add_option( plint , "Source" , "Source" , TRUE ) ;
00094
00095 PLUTO_add_dataset( plint ,
00096 "Dataset" ,
00097 ANAT_ALL_MASK ,
00098 FUNC_ALL_MASK ,
00099 DIMEN_3D_MASK |
00100 BRICK_ALLREAL_MASK
00101 ) ;
00102 PLUTO_add_number( plint , "Brick" , 0,9999,0, 0,1 ) ;
00103 PLUTO_add_string( plint , "Square" , NYESNO , YESNO_strings , 1 ) ;
00104
00105
00106
00107 PLUTO_add_option( plint , "a Params" , "Params" , TRUE ) ;
00108 PLUTO_add_number( plint , "a bot" , 2,50 ,1 , 5 , 1 ) ;
00109 PLUTO_add_number( plint , "a top" , 2,500,1 , 20 , 1 ) ;
00110
00111 PLUTO_add_option( plint , "b Params" , "Params" , TRUE ) ;
00112 PLUTO_add_number( plint , "b bot" , 10,400 ,0 , 10 , 1 ) ;
00113 PLUTO_add_number( plint , "b top" , 10,9999,0 , 200 , 1 ) ;
00114 PLUTO_add_number( plint , "H last", 0,1000,-1 , 0,1 ) ;
00115
00116 PLUTO_add_option( plint , "Misc" , "Params" , TRUE ) ;
00117 PLUTO_add_number( plint , "N ran" , 10,1000,-2 , 100 , 1 ) ;
00118 PLUTO_add_number( plint , "% cut" , 20,90,0 , 70,1 ) ;
00119 PLUTO_add_string( plint , "HSqrt" , NYESNO , YESNO_strings , 0 ) ;
00120
00121
00122
00123 PLUTO_add_option( plint , "Mask" , "Mask" , FALSE ) ;
00124 PLUTO_add_dataset( plint , "Dataset" ,
00125 ANAT_ALL_MASK , FUNC_ALL_MASK ,
00126 DIMEN_ALL_MASK | BRICK_ALLREAL_MASK ) ;
00127 PLUTO_add_number( plint , "Brick" , 0,9999,0 , 0,1 ) ;
00128
00129
00130
00131 PLUTO_add_option( plint , "Range" , "Range" , FALSE ) ;
00132 PLUTO_add_number( plint , "Bottom" , -99999,99999, 1, 0,1 ) ;
00133 PLUTO_add_number( plint , "Top" , -99999,99999,-1, 0,1 ) ;
00134
00135
00136
00137 PLUTO_add_option( plint , "Extra" , "Extra" , FALSE ) ;
00138 PLUTO_add_number( plint , "a" , 2,50,1 , 5 , 1 ) ;
00139 PLUTO_add_number( plint , "b" , 10,999,0 , 200 , 1 ) ;
00140
00141 return plint ;
00142 }
00143
00144 #include "betafit.c"
00145
00146
00147
00148 static char * BFIT_main( PLUGIN_interface * plint )
00149 {
00150 MCW_idcode * idc ;
00151 THD_3dim_dataset * input_dset , * mask_dset = NULL ;
00152
00153 BFIT_data * bfd ;
00154 BFIT_result * bfr ;
00155
00156 int nvals,ival , nran,nvox , nbin , miv , sqr,sqt ;
00157 float abot,atop,bbot,btop,pcut , eps,eps1 , hlast ;
00158 float *bval , *cval ;
00159 double aa,bb,xc ;
00160 double chq,ccc,cdf ;
00161 int ihqbot,ihqtop ;
00162
00163 int mcount,mgood , ii , jj , ibot,itop ;
00164 float mask_bot=666.0 , mask_top=-666.0 , hbot,htop,dbin ;
00165 char buf[THD_MAX_NAME+128] , tbuf[THD_MAX_NAME+128] , * tag ;
00166 int * hbin , * jbin,*kbin=NULL , *jist[2] ;
00167 MRI_IMAGE * flim ;
00168
00169 double aext=-1.0,bext=-1.0 ;
00170
00171
00172
00173
00174 if( plint == NULL )
00175 return "************************\n"
00176 "BFIT_main: NULL input\n"
00177 "************************" ;
00178
00179
00180
00181 PLUTO_next_option(plint) ;
00182 idc = PLUTO_get_idcode(plint) ;
00183 input_dset = PLUTO_find_dset(idc) ;
00184 if( input_dset == NULL )
00185 return "****************************\n"
00186 "BFIT_main: bad input dataset\n"
00187 "****************************" ;
00188
00189 nvox = DSET_NVOX(input_dset) ;
00190 nvals = DSET_NVALS(input_dset) ;
00191 ival = (int) PLUTO_get_number(plint) ;
00192 if( ival < 0 || ival >= nvals )
00193 return "**************************\n"
00194 "BFIT_main: bad Brick index\n"
00195 "**************************" ;
00196
00197 DSET_load(input_dset) ;
00198 if( DSET_ARRAY(input_dset,0) == NULL )
00199 return "*****************************\n"
00200 "BFIT_main: can't load dataset\n"
00201 "*****************************" ;
00202
00203 tag = PLUTO_get_string(plint) ;
00204 sqr = PLUTO_string_index(tag,NYESNO,YESNO_strings) ;
00205
00206
00207
00208 PLUTO_next_option(plint) ;
00209 abot = PLUTO_get_number(plint) ;
00210 atop = PLUTO_get_number(plint) ;
00211 if( atop <= abot )
00212 return "*** atop <= abot! ***" ;
00213
00214 PLUTO_next_option(plint) ;
00215 bbot = PLUTO_get_number(plint) ;
00216 btop = PLUTO_get_number(plint) ;
00217 if( atop <= abot )
00218 return "*** btop <= bbot! ***" ;
00219 hlast = PLUTO_get_number(plint) ;
00220
00221 PLUTO_next_option(plint) ;
00222 nran = (int) PLUTO_get_number(plint) ;
00223 pcut = PLUTO_get_number(plint) ;
00224
00225 tag = PLUTO_get_string(plint) ;
00226 sqt = PLUTO_string_index(tag,NYESNO,YESNO_strings) ;
00227
00228
00229
00230 while( (tag=PLUTO_get_optiontag(plint)) != NULL ){
00231
00232
00233
00234 if( strcmp(tag,"Mask") == 0 ){
00235
00236 idc = PLUTO_get_idcode(plint) ;
00237 mask_dset = PLUTO_find_dset(idc) ;
00238
00239 if( mask_dset == NULL ){
00240 return "******************************\n"
00241 "BFIT_main: bad mask dataset\n"
00242 "******************************" ;
00243 }
00244
00245 if( DSET_NVOX(mask_dset) != nvox ){
00246 return "************************************************************\n"
00247 "BFIT_main: mask input dataset doesn't match source dataset\n"
00248 "************************************************************" ;
00249 }
00250
00251 miv = (int) PLUTO_get_number(plint) ;
00252 if( miv >= DSET_NVALS(mask_dset) || miv < 0 ){
00253 return "****************************************************\n"
00254 "BFIT_main: mask dataset sub-brick index is illegal\n"
00255 "****************************************************" ;
00256 }
00257
00258 DSET_load(mask_dset) ;
00259 if( DSET_ARRAY(mask_dset,miv) == NULL ){
00260 return "*************************************\n"
00261 "BFIT_main: can't load mask dataset\n"
00262 "*************************************" ;
00263 }
00264 continue ;
00265 }
00266
00267
00268
00269 if( strcmp(tag,"Range") == 0 ){
00270 if( mask_dset == NULL ){
00271 return "******************************************\n"
00272 "BFIT_main: Can't use Range without Mask\n"
00273 "******************************************" ;
00274 }
00275
00276 mask_bot = PLUTO_get_number(plint) ;
00277 mask_top = PLUTO_get_number(plint) ;
00278 continue ;
00279 }
00280
00281
00282
00283 if( strcmp(tag,"Extra") == 0 ){
00284 aext = PLUTO_get_number(plint) ;
00285 bext = PLUTO_get_number(plint) ;
00286 continue ;
00287 }
00288 }
00289
00290
00291
00292
00293 bfd = BFIT_prepare_dataset( input_dset , ival , sqr ,
00294 mask_dset , miv , mask_bot , mask_top ) ;
00295
00296 if( bfd == NULL ) return "*** BFIT_prepare_dataset fails ***" ;
00297
00298 bfr = BFIT_compute( bfd ,
00299 pcut , abot,atop , bbot,btop , nran,200 ) ;
00300
00301 if( bfr == NULL ){
00302 BFIT_free_data( bfd ) ;
00303 return "*** BFIT_compute fails! ***" ;
00304 }
00305
00306 itop = bfr->itop ;
00307 mgood = bfr->mgood ;
00308
00309 ibot = bfd->ibot ;
00310 bval = bfd->bval ;
00311 cval = bfd->cval ;
00312 mcount = bfd->mcount ;
00313
00314 xc = bfr->xcut ;
00315 aa = bfr->a ;
00316 bb = bfr->b ;
00317 eps = bfr->eps ;
00318 eps1 = 1.0 - eps ;
00319 if( eps1 > 1.0 ) eps1 = 1.0 ;
00320 eps1 = (mcount-ibot) * eps1 ;
00321
00322
00323
00324
00325
00326 if( !sqr ){
00327 hbot = 0.0 ; htop = 1.0 ; nbin = 200 ;
00328 if( bval[mcount-1] < 1.0 ) htop = bval[mcount-1] ;
00329 dbin = (htop-hbot)/nbin ;
00330
00331 hbin = (int *) calloc((nbin+1),sizeof(int)) ;
00332 jbin = (int *) calloc((nbin+1),sizeof(int)) ;
00333
00334 for( ii=0 ; ii < nbin ; ii++ ){
00335 jbin[ii] = (int)( eps1 * ( beta_t2p(hbot+ii*dbin,aa,bb)
00336 -beta_t2p(hbot+ii*dbin+dbin,aa,bb) ) ) ;
00337 }
00338
00339 jist[0] = jbin ;
00340
00341 flim = mri_new_vol_empty( mcount-ibot,1,1 , MRI_float ) ;
00342 mri_fix_data_pointer( bval+ibot , flim ) ;
00343 mri_histogram( flim , hbot,htop , TRUE , nbin,hbin ) ;
00344
00345
00346
00347 if( aext > 0.0 ){
00348 kbin = (int *) calloc((nbin+1),sizeof(int)) ;
00349 jist[1] = kbin ;
00350 for( ii=0 ; ii < nbin ; ii++ ){
00351 kbin[ii] = (int)( eps1 * ( beta_t2p(hbot+ii*dbin,aext,bext)
00352 -beta_t2p(hbot+ii*dbin+dbin,aext,bext) ) ) ;
00353 }
00354 }
00355
00356 } else {
00357
00358 double hb,ht ;
00359 htop = 1.0 ; nbin = 200 ;
00360 if( bval[mcount-1] < 1.0 ) htop = sqrt(bval[mcount-1]) ;
00361 hbot = -htop ;
00362 dbin = (htop-hbot)/nbin ;
00363
00364 hbin = (int *) calloc((nbin+1),sizeof(int)) ;
00365 jbin = (int *) calloc((nbin+1),sizeof(int)) ;
00366
00367 for( ii=0 ; ii < nbin ; ii++ ){
00368 hb = hbot+ii*dbin ; ht = hb+dbin ;
00369 hb = hb*hb ; ht = ht*ht ;
00370 if( hb > ht ){ double qq=hb ; hb=ht ; ht=qq ; }
00371 jbin[ii] = (int)( 0.5*eps1 * ( beta_t2p(hb,aa,bb)
00372 -beta_t2p(ht,aa,bb) ) ) ;
00373 }
00374
00375 jist[0] = jbin ;
00376
00377 flim = mri_new_vol_empty( mcount-ibot,1,1 , MRI_float ) ;
00378 mri_fix_data_pointer( cval+ibot , flim ) ;
00379 mri_histogram( flim , hbot,htop , TRUE , nbin,hbin ) ;
00380
00381
00382
00383 if( aext > 0.0 ){
00384 kbin = (int *) calloc((nbin+1),sizeof(int)) ;
00385 jist[1] = kbin ;
00386 for( ii=0 ; ii < nbin ; ii++ ){
00387 hb = hbot+ii*dbin ; ht = hb+dbin ;
00388 hb = hb*hb ; ht = ht*ht ;
00389 if( hb > ht ){ double qq=hb ; hb=ht ; ht=qq ; }
00390 kbin[ii] = (int)( 0.5*eps1 * ( beta_t2p(hb,aext,bext)
00391 -beta_t2p(ht,aext,bext) ) ) ;
00392 }
00393 }
00394 }
00395
00396 sprintf(buf,"%s[%d] a=%.2f b=%.2f \\epsilon=%.2f %%=%.0f",
00397 DSET_FILECODE(input_dset),ival,aa,bb,eps,pcut ) ;
00398
00399 ccc = bfr->q_chisq ;
00400
00401
00402
00403 if( sqt ){
00404 for( ii=0 ; ii < nbin ; ii++ ){
00405 hbin[ii] = (int) sqrt( (double)(100*hbin[ii]+0.5) ) ;
00406 jbin[ii] = (int) sqrt( (double)(100*jbin[ii]+0.5) ) ;
00407 if( kbin!=NULL )
00408 kbin[ii] = (int) sqrt( (double)(100*kbin[ii]+0.5) ) ;
00409 }
00410 }
00411
00412
00413
00414 sprintf(tbuf,"\\beta fit: cutoff=%.2f nvox=%d q(\\chi^2)=%8.2e",
00415 (sqr)?sqrt(xc):xc , mgood , ccc ) ;
00416 if( sqt ){
00417 ii = strlen(tbuf) ;
00418 sprintf( tbuf+ii , " \\surd ogram" ) ;
00419 }
00420
00421 if( hlast > 0.0 ){
00422 hbin[nbin-1] = jbin[nbin-1] = hlast ;
00423 if( kbin != NULL ) kbin[nbin-1] = hlast ;
00424 }
00425
00426 PLUTO_histoplot( nbin,hbot,htop,hbin ,
00427 tbuf,NULL,buf , (kbin==NULL)?1:2 , jist ) ;
00428
00429
00430
00431 mri_clear_data_pointer(flim) ; mri_free(flim) ;
00432 free(hbin) ; free(jbin) ; if( kbin != NULL ) free(kbin);
00433
00434 BFIT_free_data(bfd) ; BFIT_free_result(bfr) ;
00435 return NULL ;
00436 }
