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 /***********************************************************************
14 Plugin to plot a histogram of a sub-brick.
15 [Interface is adapted from plug_maskave.c - RWCox - 30 Sep 1999]
16 01 Mar 2001 -- Modified to include Max Count -- RWCox
17 ************************************************************************/
18
19 static char * HISTO_main( PLUGIN_interface * ) ;
20
21 static PLUGIN_interface * CORREL_init(void) ;
22
23 static char helpstring[] =
24 " Purpose: Plot a histogram of data from a dataset brick.\n"
25 "\n"
26 " Source: Dataset = data to be processed\n"
27 " Sub-brick = which one to use\n\n"
28 " Values: Bottom = minimum value from dataset to include\n"
29 " Top = maximum value from dataset to include\n\n"
30 " Bins: Number = number of bins to use\n"
31 " Max Count = maximum count per bin\n\n"
32 " Smooth = +/- bin range to smooth histogram\n"
33 " Mask: Dataset = masking dataset\n"
34 " Sub-brick = which one to use\n\n"
35 " Range: Bottom = min value from mask dataset to use\n"
36 " Top = max value from mask dataset to use\n"
37 " [if Bottom > Top, then all nonzero mask voxels will be used; ]\n"
38 " [if Bottom <= Top, then only nonzero mask voxels in this range]\n"
39 " [ will be used in computing the statistics. ]\n"
40 " Aboot: If activated, then only voxels within a distance of Radius mm\n"
41 " of the current crosshairs will be used in the histogram.\n"
42 " Output: Name of the ascii file to which histogram values are written.\n"
43 "\n"
44 " Author -- RW Cox - 30 September 1999\n"
45 " Output feature added by V Roopchansingh, Feb 2002\n"
46 ;
47
48 /***********************************************************************
49 Set up the interface to the user
50 ************************************************************************/
51
52
53 DEFINE_PLUGIN_PROTOTYPE
54
PLUGIN_init(int ncall)55 PLUGIN_interface * PLUGIN_init( int ncall )
56 {
57 PLUGIN_interface * plint ;
58
59 #if 0
60 if( ncall > 0 ) return NULL ; /* only one interface */
61 #else
62 if( ncall == 1 ) return CORREL_init() ;
63 if( ncall > 1 ) return NULL ;
64 #endif
65
66 /*-- set titles and call point --*/
67
68 plint = PLUTO_new_interface( "Histogram" ,
69 "Histogram of Dataset Brick" ,
70 helpstring ,
71 PLUGIN_CALL_VIA_MENU , HISTO_main ) ;
72
73 PLUTO_add_hint( plint , "Histogram of Dataset Brick" ) ;
74
75 PLUTO_set_sequence( plint , "A:afniinfo:dsethistog" ) ;
76
77 PLUTO_set_runlabels( plint , "Plot+Keep" , "Plot+Close" ) ; /* 04 Nov 2003 */
78
79 /*-- first line of input --*/
80
81 PLUTO_add_option( plint , "Source" , "Source" , TRUE ) ;
82 PLUTO_add_dataset(plint , "Dataset" ,
83 ANAT_ALL_MASK , FUNC_ALL_MASK ,
84 DIMEN_ALL_MASK | BRICK_ALLREAL_MASK ) ;
85
86 PLUTO_add_number( plint , "Sub-brick" , 0,9999,0 , 0,1 ) ;
87
88 /*-- second line of input --*/
89
90 PLUTO_add_option( plint , "Values" , "Values" , FALSE ) ;
91 PLUTO_add_number( plint , "Bottom" , -99999,99999, 0, 1,1 ) ;
92 PLUTO_add_number( plint , "Top" , -99999,99999, 0, -1,1 ) ;
93
94 /*-- third line of input --*/
95
96 PLUTO_add_option( plint , "Bins" , "Bins" , FALSE ) ;
97 PLUTO_add_number( plint , "Number" , 10,1000,0, 100,1 ) ;
98 PLUTO_add_number( plint , "Max Count" , 0,999999999,0 , 0,1 ) ;
99 PLUTO_add_number( plint , "Smooth" , 0,99,0 , 0,1 ) ; /* 03 Dec 2004 */
100
101 /*-- fourth line of input --*/
102
103 PLUTO_add_option( plint , "Mask" , "Mask" , FALSE ) ;
104 PLUTO_add_dataset( plint , "Dataset" ,
105 ANAT_ALL_MASK , FUNC_ALL_MASK ,
106 DIMEN_ALL_MASK | BRICK_ALLREAL_MASK ) ;
107 PLUTO_add_number( plint , "Sub-brick" , 0,9999,0 , 0,1 ) ;
108
109 /*-- fifth line of input --*/
110
111 PLUTO_add_option( plint , "Range" , "Range" , FALSE ) ;
112 PLUTO_add_number( plint , "Bottom" , -99999,99999, 0, 1,1 ) ;
113 PLUTO_add_number( plint , "Top" , -99999,99999, 0, -1,1 ) ;
114
115 /*-- sixth line of input [20 Mar 2001] --*/
116
117 PLUTO_add_option( plint , "Aboot" , "Aboot" , FALSE ) ;
118 PLUTO_add_number( plint , "Radius" , 2,100,0,10,1 ) ;
119
120 /*-- seventh line of input [05 Feb 2002 - VR] --*/
121
122 PLUTO_add_option( plint , "Output", "Output", FALSE ) ;
123 PLUTO_add_string( plint , "Filename", 0 , NULL , 20 ) ;
124
125 return plint ;
126 }
127
128 /***************************************************************************
129 Main routine for this plugin (will be called from AFNI).
130 ****************************************************************************/
131
HISTO_main(PLUGIN_interface * plint)132 static char * HISTO_main( PLUGIN_interface * plint )
133 {
134 MCW_idcode * idc ;
135 THD_3dim_dataset * input_dset , * mask_dset=NULL ;
136 int iv , mcount , nvox , ii,jj , nbin=-1 , do_mval=0,mval ;
137 float mask_bot=666.0 , mask_top=-666.0 , hbot,htop ;
138 float val_bot=666.0 , val_top=-666.0 ;
139 char * tag , * str , buf[THD_MAX_NAME+16] ;
140 byte * mmm ;
141 MRI_IMAGE * flim ;
142 float * flar ;
143 int * hbin ;
144 int smooth=0 ; /* 03 Dec 2004 */
145
146 int miv=0 ;
147
148 int maxcount=0 ; /* 01 Mar 2001 */
149 float hrad=0.0 ; /* 20 Mar 2001 */
150
151 char * histout=NULL ; /* 05 Feb 2002 - VR */
152 FILE * HISTOUT=NULL ; /* 05 Feb 2002 - VR */
153 int writehist=0 ; /* 05 Feb 2002 - VR */
154 float dx ; /* 05 Feb 2002 - VR */
155
156 /*--------------------------------------------------------------------*/
157 /*----- Check inputs from AFNI to see if they are reasonable-ish -----*/
158
159 if( plint == NULL )
160 return "***********************\n"
161 "HISTO_main: NULL input\n"
162 "***********************" ;
163
164 /*-- read 1st line --*/
165
166 PLUTO_next_option(plint) ;
167 idc = PLUTO_get_idcode(plint) ;
168 input_dset = PLUTO_find_dset(idc) ;
169 if( input_dset == NULL )
170 return "******************************\n"
171 "HISTO_main: bad input dataset\n"
172 "******************************" ;
173
174 iv = (int) PLUTO_get_number(plint) ;
175 if( iv >= DSET_NVALS(input_dset) || iv < 0 )
176 return "********************************\n"
177 "HISTO_main: bad input sub-brick\n"
178 "********************************" ;
179
180 DSET_load(input_dset) ;
181 if( DSET_ARRAY(input_dset,iv) == NULL )
182 return "*******************************\n"
183 "HISTO_main: can't load dataset\n"
184 "*******************************" ;
185 nvox = DSET_NVOX(input_dset) ;
186
187 /*-- read optional lines --*/
188
189 while( (tag=PLUTO_get_optiontag(plint)) != NULL ){
190
191 /*-- Dataset range of values --*/
192
193 if( strcmp(tag,"Values") == 0 ){
194 val_bot = PLUTO_get_number(plint) ;
195 val_top = PLUTO_get_number(plint) ;
196 do_mval = (val_bot < val_top) ;
197 continue ;
198 }
199
200 /*-- Number of bins --*/
201
202 if( strcmp(tag,"Bins") == 0 ){
203 nbin = PLUTO_get_number(plint) ;
204 maxcount = PLUTO_get_number(plint) ;
205 smooth = PLUTO_get_number(plint) ; /* 03 Dec 2004 */
206 continue ;
207 }
208
209 /*-- Mask range of values --*/
210
211 if( strcmp(tag,"Range") == 0 ){
212 if( mask_dset == NULL )
213 return "*****************************************\n"
214 "HISTO_main: Can't use Range without Mask\n"
215 "*****************************************" ;
216
217 mask_bot = PLUTO_get_number(plint) ;
218 mask_top = PLUTO_get_number(plint) ;
219 continue ;
220 }
221
222 /*-- Mask itself --*/
223
224 if( strcmp(tag,"Mask") == 0 ){
225
226 idc = PLUTO_get_idcode(plint) ;
227 mask_dset = PLUTO_find_dset(idc) ;
228
229 if( mask_dset == NULL )
230 return "*****************************\n"
231 "HISTO_main: bad mask dataset\n"
232 "*****************************" ;
233
234 if( DSET_NVOX(mask_dset) != nvox )
235 return "***********************************************************\n"
236 "HISTO_main: mask input dataset doesn't match source dataset\n"
237 "***********************************************************" ;
238
239 miv = (int) PLUTO_get_number(plint) ; /* 06 Aug 1998 */
240 if( miv >= DSET_NVALS(mask_dset) || miv < 0 )
241 return "***************************************************\n"
242 "HISTO_main: mask dataset sub-brick index is illegal\n"
243 "***************************************************" ;
244
245 DSET_load(mask_dset) ;
246 if( DSET_ARRAY(mask_dset,miv) == NULL )
247 return "************************************\n"
248 "HISTO_main: can't load mask dataset\n"
249 "************************************" ;
250 continue ;
251 }
252
253 /*-- 20 Mar 2001: Aboot --*/
254
255 if( strcmp(tag,"Aboot") == 0 ){
256 hrad = PLUTO_get_number(plint) ;
257 continue ;
258 }
259
260 /*-- 05 Feb 2002: Output - VR --*/
261
262 if( strcmp(tag,"Output") == 0 ){
263 histout = PLUTO_get_string(plint) ;
264 writehist = 1 ;
265 continue ;
266 }
267 }
268
269 /*------------------------------------------------------*/
270 /*---------- At this point, the inputs are OK ----------*/
271
272 /*-- build a byte mask array --*/
273
274 if( mask_dset == NULL ){
275 mmm = (byte *) malloc( sizeof(byte) * nvox ) ;
276 if( mmm == NULL )
277 return " \n*** Can't malloc workspace! ***\n" ;
278 memset( mmm , 1, nvox ) ; mcount = nvox ;
279 } else {
280
281 mmm = THD_makemask( mask_dset , miv , mask_bot , mask_top ) ;
282 if( mmm == NULL )
283 return " \n*** Can't make mask for some reason! ***\n" ;
284 mcount = THD_countmask( nvox , mmm ) ;
285
286 if( !EQUIV_DSETS(mask_dset,input_dset) ) DSET_unload(mask_dset) ;
287 if( mcount < 3 ){
288 free(mmm) ;
289 return " \n*** Less than 3 voxels survive the mask! ***\n" ;
290 }
291 sprintf(buf," \n"
292 " %d voxels in the mask\n"
293 " out of %d dataset voxels\n ",mcount,nvox) ;
294 PLUTO_popup_transient(plint,buf) ;
295 }
296
297 /*-- 20 Mar 2001: modify mask via Aboot Radius, if present --*/
298
299 if( hrad > 0.0 ){
300 MCW_cluster * cl ;
301 short *di,*dj,*dk ;
302 int nd , xx,yy,zz , dd , nx,ny,nz,nxy, nx1,ny1,nz1 , ip,jp,kp ;
303
304 cl = MCW_build_mask( fabs(DSET_DX(input_dset)) ,
305 fabs(DSET_DY(input_dset)) ,
306 fabs(DSET_DZ(input_dset)) , hrad ) ;
307
308 if( cl == NULL || cl->num_pt < 6 ){
309 KILL_CLUSTER(cl);
310 PLUTO_popup_transient(plint, " \n"
311 " Aboot Radius too small\n"
312 " for this dataset!\n" ) ;
313 } else {
314 ADDTO_CLUSTER(cl,0,0,0,0) ;
315 di = cl->i ; dj = cl->j ; dk = cl->k ; nd = cl->num_pt ;
316 nx = DSET_NX(input_dset) ; nx1 = nx-1 ;
317 ny = DSET_NY(input_dset) ; ny1 = ny-1 ; nxy = nx*ny ;
318 nz = DSET_NZ(input_dset) ; nz1 = nz-1 ;
319 xx = plint->im3d->vinfo->i1 ;
320 yy = plint->im3d->vinfo->j2 ;
321 zz = plint->im3d->vinfo->k3 ;
322 for( dd=0 ; dd < nd ; dd++ ){
323 ip = xx+di[dd] ; if( ip < 0 || ip > nx1 ) continue ;
324 jp = yy+dj[dd] ; if( jp < 0 || jp > ny1 ) continue ;
325 kp = zz+dk[dd] ; if( kp < 0 || kp > nz1 ) continue ;
326 mmm[ip+jp*nx+kp*nxy]++ ;
327 }
328 KILL_CLUSTER(cl) ;
329 for( dd=0 ; dd < nvox ; dd++ ) if( mmm[dd] == 1 ) mmm[dd] = 0 ;
330 mcount = THD_countmask( nvox , mmm ) ;
331
332 if( mcount < 3 ){
333 free(mmm) ;
334 return " \n*** Less than 3 voxels survive the mask+radius! ***\n" ;
335 }
336 sprintf(buf," \n"
337 " %d voxels in the mask+radius\n"
338 " out of %d dataset voxels\n ",mcount,nvox) ;
339 PLUTO_popup_transient(plint,buf) ;
340 }
341 }
342
343 /*-- check for text output of histogram - 05 Feb 2002 - VR --*/
344
345 if ( writehist )
346 {
347 static char hbuf[1024] ;
348 if ( ( histout == NULL ) || ( strlen (histout) == 0 ) ){
349 sprintf( hbuf , "%s.histog" , DSET_PREFIX(input_dset) ) ;
350 } else {
351 strcpy( hbuf , histout ) ;
352 if( strstr(hbuf,".hist") == NULL ) strcat( hbuf , ".histog" ) ;
353 }
354 histout = hbuf ;
355
356 if (THD_is_file(histout))
357 {
358 free(mmm) ;
359
360 return "*******************************\n"
361 "Outfile exists, won't overwrite\n"
362 "*******************************\n" ;
363 }
364 else {
365 HISTOUT = fopen (histout, "w") ;
366 if( HISTOUT == NULL ){
367 free(mmm) ;
368 return "**********************************\n"
369 "Can't open Outfile for some reason\n"
370 "**********************************\n" ;
371 }
372 }
373 }
374
375 /*-- allocate an array to histogrammatize --*/
376
377 flim = mri_new( mcount , 1 , MRI_float ) ;
378 flar = MRI_FLOAT_PTR(flim) ;
379
380 /*-- load values into this array --*/
381
382 switch( DSET_BRICK_TYPE(input_dset,iv) ){
383 default:
384 free(mmm) ; mri_free(flim) ;
385 return "*** Can't use source dataset -- illegal data type! ***" ;
386
387 case MRI_short:{
388 short *bar = (short *) DSET_ARRAY(input_dset,iv) ;
389 float mfac = DSET_BRICK_FACTOR(input_dset,iv) ;
390 if( mfac == 0.0 ) mfac = 1.0 ;
391 if( do_mval ){
392 float val ;
393 for( ii=jj=0 ; ii < nvox ; ii++ ){
394 if( mmm[ii] ){
395 val = mfac*bar[ii] ;
396 if( val >= val_bot && val <= val_top ) flar[jj++] = val ;
397 }
398 }
399 mval = jj ;
400 } else {
401 for( ii=jj=0 ; ii < nvox ; ii++ )
402 if( mmm[ii] ) flar[jj++] = mfac*bar[ii] ;
403 }
404 }
405 break ;
406
407 case MRI_byte:{
408 byte *bar = (byte *) DSET_ARRAY(input_dset,iv) ;
409 float mfac = DSET_BRICK_FACTOR(input_dset,iv) ;
410 if( mfac == 0.0 ) mfac = 1.0 ;
411 if( do_mval ){
412 float val ;
413 for( ii=jj=0 ; ii < nvox ; ii++ ){
414 if( mmm[ii] ){
415 val = mfac*bar[ii] ;
416 if( val >= val_bot && val <= val_top ) flar[jj++] = val ;
417 }
418 }
419 mval = jj ;
420 } else {
421 for( ii=jj=0 ; ii < nvox ; ii++ )
422 if( mmm[ii] ) flar[jj++] = mfac*bar[ii] ;
423 }
424 }
425 break ;
426
427 case MRI_float:{
428 float *bar = (float *) DSET_ARRAY(input_dset,iv) ;
429 float mfac = DSET_BRICK_FACTOR(input_dset,iv) ;
430 if( mfac == 0.0 ) mfac = 1.0 ;
431 if( do_mval ){
432 float val ;
433 for( ii=jj=0 ; ii < nvox ; ii++ ){
434 if( mmm[ii] ){
435 val = mfac*bar[ii] ;
436 if( val >= val_bot && val <= val_top ) flar[jj++] = val ;
437 }
438 }
439 mval = jj ;
440 } else {
441 for( ii=jj=0 ; ii < nvox ; ii++ )
442 if( mmm[ii] ) flar[jj++] = mfac*bar[ii] ;
443 }
444 }
445 break ;
446 }
447
448 if( do_mval ){
449 if( mval == 0 ){
450 free(mmm) ; mri_free(flim) ;
451 return "*** Can't use source dataset -- no data in Values range ***" ;
452 }
453 flim->nx = flim->nvox = mcount = mval ;
454 }
455
456 /*-- set range and size of histogram --*/
457
458 if( val_bot > val_top ){
459 hbot = mri_min(flim) ; htop = mri_max(flim) ;
460 if( hbot >= htop ){
461 free(mmm) ; mri_free(flim) ;
462 return "***********************************\n"
463 "Selected voxels have no data range!\n"
464 "***********************************" ;
465 }
466 } else {
467 hbot = val_bot ; htop = val_top ;
468 }
469
470 if( nbin < 10 || nbin > 1000 ){
471 switch( DSET_BRICK_TYPE(input_dset,iv) ){
472 case MRI_float:
473 nbin = (int) sqrt((double)mcount) ;
474 break ;
475
476 case MRI_short:
477 case MRI_byte:{
478 float mfac = DSET_BRICK_FACTOR(input_dset,iv) ;
479 if( mfac == 0.0 || mfac == 1.0 )
480 nbin = (int)( htop - hbot ) ;
481 else
482 nbin = (int) sqrt((double)mcount) ;
483 }
484 break ;
485
486 default: break ;
487 }
488 if( nbin < 10 ) nbin = 10 ; else if( nbin > 1000 ) nbin = 1000 ;
489 }
490
491 /*-- actually compute and plot histogram --*/
492
493 hbin = (int *) calloc((nbin+1),sizeof(int)) ;
494
495 mri_histogram( flim , hbot,htop , TRUE , nbin,hbin ) ;
496
497 if( smooth > 0 ){ /* 03 Dec 2004 */
498 int nwid=smooth , *gbin=(int *)calloc((nbin+1),sizeof(int)) , ibot,itop ;
499 float ws,wss , *wt ;
500
501 ws = 0.0 ;
502 wt = (float *)malloc(sizeof(float)*(2*nwid+1)) ;
503 for( ii=0 ; ii <= 2*nwid ; ii++ ){
504 wt[ii] = nwid-abs(nwid-ii) + 0.5f ;
505 ws += wt[ii] ;
506 }
507 for( ii=0 ; ii <= 2*nwid ; ii++ ) wt[ii] /= ws ;
508
509 for( jj=0 ; jj <= nbin ; jj++ ){
510 ibot = jj-nwid ; if( ibot < 0 ) ibot = 0 ;
511 itop = jj+nwid ; if( itop > nbin ) itop = nbin ;
512 ws = wss = 0.0 ;
513 for( ii=ibot ; ii <= itop ; ii++ ){
514 ws += wt[nwid-jj+ii] * hbin[ii] ; wss += wt[nwid-jj+ii] ;
515 }
516 gbin[jj] = rint(ws/wss) ;
517 }
518 memcpy(hbin,gbin,sizeof(int)*(nbin+1)) ;
519 free((void *)wt) ; free((void *)gbin) ;
520 }
521
522 if( maxcount > 0 ){
523 for( ii=0 ; ii <= nbin ; ii++ ) hbin[ii] = MIN( hbin[ii] , maxcount ) ;
524 }
525 hrad = AFNI_numenv("AFNI_1DPLOT_THIK") ;
526 if( hrad <= 0.0f || hrad >= 0.02f ) hrad = 0.004f ;
527 plot_ts_setTHIK(hrad) ; plot_ts_setthik(0.0015f) ;
528 plot_ts_xypush(0,-1) ;
529 sprintf(buf,"\\noesc %s[%d] %s voxels",DSET_FILECODE(input_dset),iv,commaized_integer_string(mcount));
530 PLUTO_histoplot( nbin,hbot,htop,hbin , NULL , NULL , buf , 0,NULL ) ;
531
532 /*-- 05 Feb 2002: Output - VR --*/
533
534 if ( HISTOUT != NULL )
535 {
536 if( hbot >= htop ){ hbot = 0.0 ; htop = nbin ;}
537
538 dx = (htop-hbot)/nbin ;
539
540 for( ii=0 ; ii <= nbin ; ii++ )
541 fprintf (HISTOUT, "%12.6f %13d \n", hbot+ii*dx, hbin[ii]) ;
542
543 fclose (HISTOUT) ;
544
545 fprintf (stderr, "%s written to disk \n", histout) ;
546 }
547
548 /*-- go home to mama --*/
549
550 free(hbin) ; free(mmm) ; mri_free(flim) ; return NULL ;
551 }
552
553 /******************************************************************************/
554
555 #define FIT_FISHER
556 #undef DO_GREEN
557
558 static char c_helpstring[] =
559 "Purpose: Plot a histogram of correlation coefficient of a 3D+time\n"
560 " dataset with an input time series file.\n"
561 "\n"
562 "Source: Dataset = 3D+time dataset to use\n"
563 " Ignore = number of points at beginning to ignore\n"
564 "\n"
565 "Vector: 1D File = time series to correlate the dataset with\n"
566 " Polort = degree of polynomial detrending to use\n"
567 "\n"
568 "Mask: Dataset = optional mask dataset to use, to restrict the voxels\n"
569 " which will be processed\n"
570 " Sub-brick = which sub-brick of the mask dataset to use\n"
571 "\n"
572 "Range: Bottom = minimum value of mask dataset to use\n"
573 " Top = maximum value of mask dataset to use\n"
574 " [default is to use all nonzero values in the mask]\n"
575 "\n"
576 "All selected voxel time series from the source are correlated with the\n"
577 "input 1D vector, using the same algorithms as the AFNI internal FIM.\n"
578 "The array of correlation coefficients is then put into 100 bins, ranging\n"
579 "from -1.0 to 1.0, and the histogram graph is popped up to the display.\n"
580 "\n"
581 "Overlaid on the histogram are other graphs:\n"
582 " * The first [red] is a normal fit to the Fisher z-transform of the\n"
583 " correlation coefficient.\n"
584 #ifdef DO_GREEN
585 " * The second [green] is the nominal fit to the number of degrees of\n"
586 " freedom, assuming the data time series are normal white noise.\n"
587 #endif
588 "\n"
589 "-- Bob Cox - October 1999\n"
590 ;
591
592 static char * CORREL_main( PLUGIN_interface * ) ;
593
CORREL_init(void)594 static PLUGIN_interface * CORREL_init(void)
595 {
596 PLUGIN_interface * plint ;
597
598 /*-- set titles and call point --*/
599
600 plint = PLUTO_new_interface( "Histogram: CC" ,
601 "Histogram of Correlation Coefficient" ,
602 c_helpstring ,
603 PLUGIN_CALL_VIA_MENU , CORREL_main ) ;
604
605 PLUTO_add_hint( plint , "Histogram: Correlation Coefficient" ) ;
606
607 PLUTO_set_sequence( plint , "A:afniinfo:dset" ) ;
608
609 /*-- first line of input --*/
610
611 PLUTO_add_option( plint , "Source" , "Source" , TRUE ) ;
612
613 PLUTO_add_dataset( plint ,
614 "Dataset" , /* label next to button */
615 ANAT_ALL_MASK , /* take any anat datasets */
616 FUNC_FIM_MASK , /* only allow fim funcs */
617 DIMEN_4D_MASK | /* need 3D+time datasets */
618 BRICK_ALLREAL_MASK /* need real-valued datasets */
619 ) ;
620
621 PLUTO_add_number( plint ,
622 "Ignore" , /* label next to chooser */
623 0 , /* smallest possible value */
624 999 , /* largest possible value */
625 0 , /* decimal shift (none in this case) */
626 4 , /* default value */
627 TRUE /* allow user to edit value? */
628 ) ;
629
630 /*-- second line of input --*/
631
632 PLUTO_add_option( plint , "Vector" , "Vector" , TRUE ) ;
633
634 PLUTO_add_timeseries( plint , "1D File" ) ;
635
636 PLUTO_add_number( plint , "Polort" , 0,MAX_POLORT,0,1,FALSE ) ;
637
638 /*-- (optional) third line of input --*/
639
640 PLUTO_add_option( plint , "Mask" , "Mask" , FALSE ) ;
641 PLUTO_add_dataset( plint , "Dataset" ,
642 ANAT_ALL_MASK , FUNC_ALL_MASK ,
643 DIMEN_ALL_MASK | BRICK_ALLREAL_MASK ) ;
644 PLUTO_add_number( plint , "Sub-brick" , 0,9999,0 , 0,1 ) ;
645
646 /*-- (optional) fourth line of input --*/
647
648 PLUTO_add_option( plint , "Range" , "Range" , FALSE ) ;
649 PLUTO_add_number( plint , "Bottom" , -99999,99999, 1, 0,1 ) ;
650 PLUTO_add_number( plint , "Top" , -99999,99999,-1, 0,1 ) ;
651
652
653 return plint ;
654 }
655
656 /*-----------------------------------------------------------------------------*/
657
658 #include "uuu.c"
659
CORREL_main(PLUGIN_interface * plint)660 static char * CORREL_main( PLUGIN_interface * plint )
661 {
662 MCW_idcode * idc ;
663 THD_3dim_dataset * input_dset , * mask_dset = NULL ;
664 MRI_IMAGE * tsim , * flim ;
665 int ignore , nvox , ntim , polort , miv=0 , it , ip , nupdt , nbin ;
666 int mcount , ii , jj ;
667 float * tsar ;
668 float mask_bot=666.0 , mask_top=-666.0 , hbot,htop ;
669 byte * mmm ;
670 char buf[THD_MAX_NAME+16] , * tag ;
671
672 PCOR_references * pc_ref ;
673 PCOR_voxel_corr * pc_vc ;
674 int fim_nref ;
675 float * ref_vec , * vval , * zval , * aval ;
676 int * hbin , * jbin , * kbin , *jist[2] ;
677 float sum , sumq , dbin , gval,rval,gg , sqp , zmid,zmed,zsig ;
678 float pstar,zstar,zplus,zminus,psum,msum ;
679
680 /*--------------------------------------------------------------------*/
681 /*----- Check inputs from AFNI to see if they are reasonable-ish -----*/
682
683 if( plint == NULL )
684 return "************************\n"
685 "CORREL_main: NULL input\n"
686 "************************" ;
687
688 /*-- read 1st line --*/
689
690 PLUTO_next_option(plint) ;
691 idc = PLUTO_get_idcode(plint) ;
692 input_dset = PLUTO_find_dset(idc) ;
693 if( input_dset == NULL )
694 return "*******************************\n"
695 "CORREL_main: bad input dataset\n"
696 "*******************************" ;
697
698 nvox = DSET_NVOX(input_dset) ;
699 ntim = DSET_NUM_TIMES(input_dset) ;
700
701 ignore = (int) PLUTO_get_number(plint) ;
702 if( ignore >= ntim-5 || ignore < 0 )
703 return "******************************\n"
704 "CORREL_main: bad ignore count\n"
705 "******************************" ;
706
707 DSET_load(input_dset) ;
708 if( DSET_ARRAY(input_dset,0) == NULL )
709 return "********************************\n"
710 "CORREL_main: can't load dataset\n"
711 "********************************" ;
712
713 /*-- read 2nd line --*/
714
715 PLUTO_next_option(plint) ;
716 tsim = PLUTO_get_timeseries(plint) ;
717 if( tsim == NULL || tsim->nx < ntim )
718 return "*****************************\n"
719 "CORREL_main: bad input vector\n"
720 "*****************************" ;
721
722 flim = mri_to_float(tsim) ; tsar = MRI_FLOAT_PTR(flim) ;
723
724 polort = (int) PLUTO_get_number(plint) ;
725
726 /*-- read optional lines --*/
727
728 while( (tag=PLUTO_get_optiontag(plint)) != NULL ){
729
730 /*-- Mask range of values --*/
731
732 if( strcmp(tag,"Range") == 0 ){
733 if( mask_dset == NULL ){
734 mri_free(flim) ;
735 return "******************************************\n"
736 "CORREL_main: Can't use Range without Mask\n"
737 "******************************************" ;
738 }
739
740 mask_bot = PLUTO_get_number(plint) ;
741 mask_top = PLUTO_get_number(plint) ;
742 continue ;
743 }
744
745 /*-- Mask itself --*/
746
747 if( strcmp(tag,"Mask") == 0 ){
748
749 idc = PLUTO_get_idcode(plint) ;
750 mask_dset = PLUTO_find_dset(idc) ;
751
752 if( mask_dset == NULL ){
753 mri_free(flim) ;
754 return "******************************\n"
755 "CORREL_main: bad mask dataset\n"
756 "******************************" ;
757 }
758
759 if( DSET_NVOX(mask_dset) != nvox ){
760 mri_free(flim) ;
761 return "************************************************************\n"
762 "CORREL_main: mask input dataset doesn't match source dataset\n"
763 "************************************************************" ;
764 }
765
766 miv = (int) PLUTO_get_number(plint) ;
767 if( miv >= DSET_NVALS(mask_dset) || miv < 0 ){
768 mri_free(flim) ;
769 return "****************************************************\n"
770 "CORREL_main: mask dataset sub-brick index is illegal\n"
771 "****************************************************" ;
772 }
773
774 DSET_load(mask_dset) ;
775 if( DSET_ARRAY(mask_dset,miv) == NULL ){
776 mri_free(flim) ;
777 return "*************************************\n"
778 "CORREL_main: can't load mask dataset\n"
779 "*************************************" ;
780 }
781 continue ;
782 }
783 }
784
785 /*------------------------------------------------------*/
786 /*---------- At this point, the inputs are OK ----------*/
787
788 /*-- build a byte mask array --*/
789
790 if( mask_dset == NULL ){
791 mmm = (byte *) malloc( sizeof(byte) * nvox ) ;
792 if( mmm == NULL )
793 return " \n*** Can't malloc workspace! ***\n" ;
794 memset( mmm , 1, nvox ) ; mcount = nvox ;
795 } else {
796
797 mmm = THD_makemask( mask_dset , miv , mask_bot , mask_top ) ;
798 if( mmm == NULL )
799 return " \n*** Can't make mask for some reason! ***\n" ;
800 mcount = THD_countmask( nvox , mmm ) ;
801
802 if( !EQUIV_DSETS(mask_dset,input_dset) ) DSET_unload(mask_dset) ;
803 if( mcount < 3 ){
804 free(mmm) ;
805 return " \n*** Less than 3 voxels survive the mask! ***\n" ;
806 }
807 sprintf(buf," \n"
808 " %d voxels in the mask\n"
809 " out of %d dataset voxels\n ",mcount,nvox) ;
810 PLUTO_popup_transient(plint,buf) ;
811 }
812
813 /*-- setup to do the FIM calculation --*/
814
815 fim_nref = polort+2 ;
816 ref_vec = (float *) malloc( sizeof(float) * fim_nref ) ;
817 vval = (float *) malloc( sizeof(float) * mcount ) ;
818
819 pc_ref = new_PCOR_references( fim_nref ) ;
820 pc_vc = new_PCOR_voxel_corr( mcount , fim_nref ) ;
821
822 /*-- do FIM --*/
823
824 for( nupdt=0,it=ignore ; it < ntim ; it++ ){
825 if( tsar[it] >= WAY_BIG ) continue ; /* skip this */
826
827 ref_vec[0] = 1.0 ;
828 for( ip=1 ; ip <= polort ; ip++ )
829 ref_vec[ip] = ref_vec[ip-1] * ((float)it) ;
830
831 ref_vec[ip] = tsar[it] ; /* vector value */
832
833 update_PCOR_references( ref_vec , pc_ref ) ;
834
835 /*-- load values into vval --*/
836
837 switch( DSET_BRICK_TYPE(input_dset,it) ){
838
839 case MRI_short:{
840 short * bar = (short *) DSET_ARRAY(input_dset,it) ;
841 float mfac = DSET_BRICK_FACTOR(input_dset,it) ;
842 if( mfac == 0.0 ) mfac = 1.0 ;
843 for( ii=jj=0 ; ii < nvox ; ii++ )
844 if( mmm[ii] ) vval[jj++] = mfac*bar[ii] ;
845 }
846 break ;
847
848 case MRI_byte:{
849 byte * bar = (byte *) DSET_ARRAY(input_dset,it) ;
850 float mfac = DSET_BRICK_FACTOR(input_dset,it) ;
851 if( mfac == 0.0 ) mfac = 1.0 ;
852 for( ii=jj=0 ; ii < nvox ; ii++ )
853 if( mmm[ii] ) vval[jj++] = mfac*bar[ii] ;
854 }
855 break ;
856
857 case MRI_float:{
858 float * bar = (float *) DSET_ARRAY(input_dset,it) ;
859 float mfac = DSET_BRICK_FACTOR(input_dset,it) ;
860 if( mfac == 0.0 ) mfac = 1.0 ;
861 for( ii=jj=0 ; ii < nvox ; ii++ )
862 if( mmm[ii] ) vval[jj++] = mfac*bar[ii] ;
863 }
864 break ;
865
866 default: break ;
867 }
868
869 PCOR_update_float( vval , pc_ref , pc_vc ) ;
870 nupdt++ ;
871 }
872
873 free(ref_vec) ; mri_free(flim) ; free(mmm) ;
874
875 /*-- get correlation coefficient --*/
876
877 PCOR_get_pcor( pc_ref , pc_vc , vval ) ;
878
879 free_PCOR_references(pc_ref) ; free_PCOR_voxel_corr(pc_vc) ;
880
881 /*-- compute statistics --*/
882
883 sum = 0.0 ;
884 for( ii=0 ; ii < mcount ; ii++ ) sum += vval[ii] ;
885 sum /= mcount ; sumq = 0.0 ;
886 for( ii=0 ; ii < mcount ; ii++ )
887 sumq += (vval[ii]-sum)*(vval[ii]-sum) ;
888 sumq = sqrt(sumq/mcount) ;
889
890 /*-- get robust statistics of Fisher z-transform --*/
891
892 zval = (float *) malloc( sizeof(float) * mcount ) ;
893 aval = (float *) malloc( sizeof(float) * mcount ) ;
894 for( ii=0 ; ii < mcount ; ii++ ) zval[ii] = atanh(vval[ii]) ;
895 qsort_float( mcount , zval ) ;
896 if( mcount%2 == 1 ) /* median */
897 zmid = zval[mcount/2] ;
898 else
899 zmid = 0.5 * ( zval[mcount/2] + zval[mcount/2-1] ) ;
900
901 for( ii=0 ; ii < mcount ; ii++ ) aval[ii] = fabs(zval[ii]-zmid) ;
902 qsort_float( mcount , aval ) ;
903 if( mcount%2 == 1 ) /* MAD = median absolute deviation */
904 zmed = aval[mcount/2] ;
905 else
906 zmed = 0.5 * ( aval[mcount/2] + aval[mcount/2-1] ) ;
907 zsig = 1.4826 * zmed ; /* estimate st.dev. */
908 /* 1/1.4826 = sqrt(2)*erfinv(0.5) */
909 free(aval) ;
910 free(zval) ;
911
912 /*-- do histogram --*/
913
914 hbot = -1.0 ; htop = 1.0 ; nbin = 100 ; dbin = (htop-hbot)/nbin ;
915
916 hbin = (int *) calloc((nbin+1),sizeof(int)) ;
917 jbin = (int *) calloc((nbin+1),sizeof(int)) ; /* 04 Oct 1999 */
918
919 #ifndef FIT_FISHER
920 sqp = 1.0/(sqrt(2.0*PI)*sumq) ; /* Gaussian fit */
921 for( ii=0 ; ii < nbin ; ii++ ){ /* to rho data */
922 gval = hbot + (ii+0.5)*dbin - sum ;
923 gval = sqp * exp( -0.5*gval*gval/(sumq*sumq) ) ;
924 jbin[ii] = (int)( mcount * dbin * gval + 0.5 ) ;
925 }
926 #else
927 sqp = 1.0/(sqrt(2.0*PI)*zsig) ; /* Gaussian fit */
928 for( ii=0 ; ii < nbin ; ii++ ){ /* to z(rho) */
929 rval = hbot + (ii+0.5)*dbin ;
930 gval = atanh(rval) - zmid ;
931 gval = sqp * exp( -0.5*gval*gval/(zsig*zsig) )/sqrt(1.0-rval*rval) ;
932 jbin[ii] = (int)( mcount * dbin * gval + 0.5 ) ;
933 }
934 sum = tanh(zmid) ;
935 sumq = 0.5*( tanh(zmid+zmed) - tanh(zmid-zmed) ) ;
936 #endif
937
938 kbin = (int *) calloc((nbin+1),sizeof(int)) ;
939 for( ii=0 ; ii < nbin ; ii++ ){
940 gval = correl_t2p( fabs(hbot+ii*dbin) ,
941 (double)nupdt , (double)1 , (double)(polort+1) ) ;
942 sqp = correl_t2p( fabs(hbot+(ii+1)*dbin) ,
943 (double)nupdt , (double)1 , (double)(polort+1) ) ;
944 kbin[ii] = (int)( 0.5*mcount*fabs(gval-sqp) ) ;
945 }
946 jist[0] = jbin ; jist[1] = kbin ;
947
948 flim = mri_new_vol_empty( mcount,1,1 , MRI_float ) ;
949 mri_fix_data_pointer( vval , flim ) ;
950 mri_histogram( flim , hbot,htop , TRUE , nbin,hbin ) ;
951
952 { char * ps = my_getenv("PTAIL") ;
953 float pp=0.0 ;
954 if( ps != NULL ) pp = strtod(ps,NULL) ;
955 set_unusuality_tail(pp) ;
956 }
957
958 for( ii=0 ; ii < mcount ; ii++ ) vval[ii] = -vval[ii] ;
959 msum = unusuality( mcount , vval ) ;
960 for( ii=0 ; ii < mcount ; ii++ ) vval[ii] = -vval[ii] ;
961 psum = unusuality( mcount , vval ) ;
962
963 sprintf(buf,"%s^.%s:\\rho_{mid}=%.2f\\pm%.2f,u=%.1f,%.1f",
964 DSET_FILECODE(input_dset),
965 (tsim->name != NULL) ? THD_trailname(tsim->name,0) : " " ,
966 sum,sumq , psum,msum ) ;
967 plot_ts_xypush(0,-1) ;
968 #ifdef DO_GREEN
969 PLUTO_histoplot( nbin,hbot,htop,hbin ,
970 "Correlation Coefficient",NULL,buf , 2,jist ) ;
971 #else
972 PLUTO_histoplot( nbin,hbot,htop,hbin ,
973 "Correlation Coefficient",NULL,buf , 1,jist ) ;
974 #endif
975
976 mri_clear_data_pointer(flim) ; mri_free(flim) ;
977 free(vval) ; free(hbin) ; free(jbin) ; free(kbin) ;
978 return NULL ;
979 }
980