forked from neurodebian/afni_removeme_eventually
-
Notifications
You must be signed in to change notification settings - Fork 0
/
3dFDR.c
1348 lines (1058 loc) · 44 KB
/
3dFDR.c
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
/*****************************************************************************
Major portions of this software are copyrighted by the Medical College
of Wisconsin, 2002, and are released under the Gnu General Public
License, Version 2. See the file README.Copyright for details.
******************************************************************************/
/*---------------------------------------------------------------------------*/
/*
This program implements the False Discovery Rate (FDR) algorithm for
thresholding of voxelwise statistics.
File: 3dFDR.c
Author: B. Douglas Ward
Date: 31 January 2002
*/
/*---------------------------------------------------------------------------*/
#define PROGRAM_NAME "3dFDR" /* name of this program */
#define PROGRAM_AUTHOR "B. Douglas Ward" /* program author */
#define PROGRAM_INITIAL "31 January 2002" /* date of initial program release */
#define PROGRAM_LATEST "18 January 2008" /* date of last program revision */
/*---------------------------------------------------------------------------*/
/*
Include header files and source code files.
*/
#include "mrilib.h"
/*---------------------------------------------------------------------------*/
/*
Structure declarations
*/
struct voxel;
typedef struct voxel
{
int ixyz; /* voxel index */
float pvalue; /* input p-value or output q-value */
struct voxel * next_voxel; /* pointer to next voxel in list */
} voxel;
/*-------------------------- global data ------------------------------------*/
#define FDR_MAX_LL 10000 /* maximum number of linked lists of voxels */
static int FDR_quiet = 0; /* flag for suppress screen output */
static int FDR_list = 0; /* flag for list voxel q-values */
static float FDR_mask_thr = 1.0; /* mask threshold */
static float FDR_cn = 1.0; /* statistical constant c(N) */
static int FDR_nxyz = 0; /* dataset dimensions in voxels */
static int FDR_nthr = 0; /* number of voxels in mask */
static char * FDR_input_filename = NULL; /* input 3d functional dataset */
static char * FDR_input1D_filename = NULL; /* input list of p-values */
static char * FDR_mask_filename = NULL; /* input 3d mask dataset */
static char * FDR_output_prefix = NULL; /* name for output 3d dataset */
static byte * FDR_mask = NULL; /* mask for voxels above thr. */
static float * FDR_input1D_data = NULL; /* input array of p-values */
static voxel * FDR_head_voxel[FDR_MAX_LL]; /* linked lists of voxels */
static char * commandline = NULL; /* command line for history notes */
static THD_3dim_dataset * FDR_dset = NULL; /* input dataset */
/*---------------------------------------------------------------------------*/
/*** 18 Jan 2008 changes [RWCox]:
* replace FDR calculation with mri_fdrize() unless -old is given
* add -force option to force conversion, assuming input is p-values
* add -pmask and -nopmask options
* add -float option
* add -qval option
-----------------------------------------------------------------------------*/
static int FDR_old = 0 ; /* new mode is on by default */
static int FDR_force = 0 ; /* only if the user asks */
static int FDR_pmask = 1 ; /* on by default in new mode */
static int FDR_float = 0 ; /* must be turned on by user */
static int FDR_qval = 0 ; /* must be turned on by user */
static int FDR_curve = 0 ; /* hidden option: -curve */
/*---------------------------------------------------------------------------*/
/** macro to open a dataset and make it ready for processing **/
#define DOPEN(ds,name) \
do{ int pv ; (ds) = THD_open_dataset((name)) ; \
CHECK_OPEN_ERROR((ds),(name)) ; \
DSET_load((ds)) ; \
pv = DSET_PRINCIPAL_VALUE((ds)) ; \
if( DSET_ARRAY((ds),pv) == NULL ){ \
fprintf(stderr,"*** Can't access data: %s\n",(name)) ; exit(1); } \
if( DSET_BRICK_TYPE((ds),pv) == MRI_complex ){ \
fprintf(stderr,"*** Can't use complex data: %s\n",(name)) ; exit(1);\
} \
break ; } while (0)
/*---------------------------------------------------------------------------*/
/** macro to return pointer to correct location in brick for current processing **/
#define SUB_POINTER(ds,vv,ind,ptr) \
do{ switch( DSET_BRICK_TYPE((ds),(vv)) ){ \
default: fprintf(stderr,"\n*** Illegal datum! ***\n");exit(1); \
case MRI_short:{ short * fim = (short *) DSET_ARRAY((ds),(vv)) ; \
(ptr) = (void *)( fim + (ind) ) ; \
} break ; \
case MRI_byte:{ byte * fim = (byte *) DSET_ARRAY((ds),(vv)) ; \
(ptr) = (void *)( fim + (ind) ) ; \
} break ; \
case MRI_float:{ float * fim = (float *) DSET_ARRAY((ds),(vv)) ; \
(ptr) = (void *)( fim + (ind) ) ; \
} break ; } break ; } while(0)
/*---------------------------------------------------------------------------*/
/*
Print error message and stop.
*/
#if 0
void FDR_error (char * message)
{
fprintf (stderr, "%s Error: %s \n", PROGRAM_NAME, message);
exit(1);
}
#else
# define FDR_error(s) ERROR_exit("3dFDR -- %s",s)
#endif
/*---------------------------------------------------------------------------*/
/** macro to test a malloc-ed pointer for validity **/
#define MTEST(ptr) \
if((ptr)==NULL) \
( FDR_error ("Cannot allocate memory") )
/*---------------------------------------------------------------------------*/
/*
Display help file.
*/
void FDR_Syntax(void)
{
printf(
"This program implements the False Discovery Rate (FDR) algorithm for \n"
"thresholding of voxelwise statistics. \n"
" \n"
"Program input consists of a functional dataset containing one (or more) \n"
"statistical sub-bricks. Output consists of a bucket dataset with one \n"
"sub-brick for each input sub-brick. For non-statistical input sub-bricks, \n"
"the output is a copy of the input. However, statistical input sub-bricks \n"
"are replaced by their corresponding FDR values, as follows: \n"
" \n"
"For each voxel, the minimum value of q is determined such that \n"
" E(FDR) <= q \n"
"leads to rejection of the null hypothesis in that voxel. Only voxels inside\n"
"the user specified mask will be considered. These q-values are then mapped\n"
"to z-scores for compatibility with the AFNI statistical threshold display: \n"
" \n"
" stat ==> p-value ==> FDR q-value ==> FDR z-score \n"
" \n"
);
printf(
"Usage: \n"
" 3dFDR \n"
" -input fname fname = filename of input 3d functional dataset \n"
" OR \n"
" -input1D dname dname = .1D file containing column of p-values \n"
" \n"
" -mask_file mname Use mask values from file mname. \n"
" *OR* Note: If file mname contains more than 1 sub-brick, \n"
" -mask mname the mask sub-brick must be specified! \n"
" Default: No mask \n"
" ** Generally speaking, you really should use a mask \n"
" to avoid counting non-brain voxels. However, with \n"
" the changes described below, the program will \n"
" automatically ignore voxels where the statistics \n"
" are set to 0, so if the program that created the \n"
" dataset used a mask, then you don't need one here. \n"
" \n"
" -mask_thr m Only voxels whose corresponding mask value is \n"
" greater than or equal to m in absolute value will \n"
" be considered. Default: m=1 \n"
" \n"
" Constant c(N) depends on assumption about p-values: \n"
" -cind c(N) = 1 p-values are independent across N voxels \n"
" -cdep c(N) = sum(1/i), i=1,...,N any joint distribution \n"
" Default: c(N) = 1 \n"
" \n"
" -quiet Flag to suppress screen output \n"
" \n"
" -list Write sorted list of voxel q-values to screen \n"
" \n"
" -prefix pname Use 'pname' for the output dataset prefix name. \n"
" OR \n"
" -output pname \n"
" \n"
"\n") ;
printf(
"===========================================================================\n"
"\n"
"January 2008: Changes to 3dFDR\n"
"------------------------------\n"
"The default mode of operation of 3dFDR has altered somewhat:\n"
"\n"
" * Voxel p-values of exactly 1 (e.g., from t=0 or F=0 or correlation=0)\n"
" are ignored by default; in the old mode of operation, they were\n"
" included in the count which goes into the FDR algorithm. The old\n"
" process tends to increase the q-values and so decrease the z-scores.\n"
"\n"
" * The array of voxel p-values are now sorted via Quicksort, rather than\n"
" by binning, as in the old mode. This (by itself) probably has no\n"
" discernible effect on the results, but should be faster.\n"
"\n"
"New Options:\n"
"------------\n"
" -old = Use the old mode of operation (for compatibility/nostalgia)\n"
" -new = Use the new mode of operation [now the default]\n"
" N.B.: '-list' does not work in the new mode!\n"
" -pmask = Instruct the program to ignore p=1 voxels\n"
" [the default in the new mode, but not in the old mode]\n"
" N.B.: voxels that were masked in 3dDeconvolve (etc.)\n"
" will have their statistics set to 0, which means p=1,\n"
" which means that such voxels are implicitly masked\n"
" with '-new', and so don't need to be explicitly\n"
" masked with the '-mask' option.\n"
" -nopmask = Instruct the program to count p=1 voxels\n"
" [the default in the old mode, but NOT in the new mode]\n"
" -force = Force the conversion of all sub-bricks, even if they\n"
" are not marked as with a statistical code; such\n"
" sub-bricks are treated as though they were p-values.\n"
" -float = Force the output of z-scores in floating point format.\n"
" -qval = Force the output of q-values rather than z-scores.\n"
" N.B.: A smaller q-value is more significant!\n"
" [-float is strongly recommended when -qval is used]\n"
"\n"
"* To be clear, you can use '-new -nopmask' to have the new mode of computing\n"
" carried out, but with p=1 voxels included (which should give results\n"
" nearly identical to '-old').\n"
"\n"
"* Or you can use '-old -pmask' to use the old mode of computing but where\n"
" p=1 voxels are not counted (which should give results virtually\n"
" identical to '-new').\n"
"\n"
"* However, the combination of '-new', '-nopmask' and '-mask_file' does not\n"
" work -- if you try it, '-pmask' will be turned back on and a warning\n"
" message printed to aid your path towards elucidation and enlightenment.\n"
"\n"
"Other Notes:\n"
"------------\n"
"* '3drefit -addFDR' can be used to add FDR curves of z(q) as a function\n"
" of threshold for all statistic sub-bricks in a dataset; in turn, these\n"
" curves let you see the (estimated) q-value as you move the threshold\n"
" slider in AFNI.\n"
" - Since 3drefit doesn't have a '-mask' option, you will have to mask\n"
" statistical sub-bricks yourself via 3dcalc (if desired):\n"
" 3dcalc -a stat+orig -b mask+orig -expr 'a*step(b)' -prefix statmm\n"
" - '-addFDR' runs as if '-new -pmask' were given to 3dFDR, so that\n"
" stat values == 0 are ignored in the FDR calculations.\n"
" - most AFNI statistical programs now automatically add FDR curves to\n"
" the output dataset header, so you can see the q-value as you adjust\n"
" the threshold slider.\n"
"\n"
"* q-values are estimates of the False Discovery Rate at a given threshold;\n"
" that is, about 5%% of all voxels with q <= 0.05 (z >= 1.96) are\n"
" (presumably) 'false positive' detections, and the other 95%% are\n"
" (presumably) 'true positives'. Of course, there is no way to tell\n"
" which above-threshold voxels are 'true' detections and which are 'false'.\n"
"\n"
"* Note the use of the words 'estimate' and 'about' in the above statement!\n"
" In particular, the accuracy of the q-value calculation depends on the\n"
" assumption that the p-values calculated from the input statistics are\n"
" correctly distributed (e.g., that the DOF parameters are correct).\n"
"\n"
"* The z-score is the conversion of the q-value to a double-sided tail\n"
" probability of the unit Gaussian N(0,1) distribution; that is, z(q)\n"
" is the value such that if x is a N(0,1) random variable, then\n"
" Prob[|x|>z] = q: for example, z(0.05) = 1.95996.\n"
" The reason for using z-scores here is simply that their range is\n"
" highly compressed relative to the range of q-values\n"
" (e.g., z(1e-9) = 6.10941), so z-scores are easily stored as shorts,\n"
" whereas q-values are much better stored as floats.\n"
"\n"
"* Changes above by RWCox -- 18 Jan 2008 == Cary Grant's Birthday!\n"
"\n"
"26 Mar 2009 -- Yet Another Change [RWCox]\n"
"-----------------------------------------\n"
"* FDR calculations in AFNI now 'adjust' the q-values downwards by\n"
" estimating the number of true negatives [m0 in the statistics\n"
" literature], and then reporting\n"
" q_new = q_old * m0 / m, where m = number of voxels being tested.\n"
" If you do NOT want this adjustment, then set environment variable\n"
" AFNI_DONT_ADJUST_FDR to YES. You can do this on the 3dFDR command\n"
" line with the option '-DAFNI_DONT_ADJUST_FDR=YES'\n"
"\n"
"For Further Reading and Amusement\n"
"---------------------------------\n"
"* cf. http://en.wikipedia.org/wiki/False_discovery_rate [Easy overview of FDR]\n"
"* cf. http://dx.doi.org/10.1093/bioinformatics/bti448 [False Negative Rate]\n"
"* cf. http://dx.doi.org/10.1093/biomet/93.3.491 [m0 adjustment idea]\n"
"* cf. C implementation in mri_fdrize.c [trust in the Source]\n"
"* cf. http://afni.nimh.nih.gov/pub/dist/doc/misc/FDR/FDR_Jan2008.pdf\n"
) ;
PRINT_COMPILE_DATE ;
exit(0) ;
}
/*---------------------------------------------------------------------------*/
/*
Read the arguments, load the global variables
*/
void read_options ( int argc , char * argv[] )
{
int nopt = 1 ; /* count of input arguments */
char message[80]; /* error message */
if( AFNI_yesenv("AFNI_FLOATIZE") ) FDR_float = 1 ;
/*----- main loop over input options -----*/
while( nopt < argc )
{
/*----- -input fname -----*/
if (strcmp(argv[nopt], "-input") == 0)
{
nopt++;
if (nopt >= argc) FDR_error ("need argument after -input ");
FDR_input_filename = (char *) malloc (sizeof(char)*THD_MAX_NAME);
MTEST (FDR_input_filename);
strcpy (FDR_input_filename, argv[nopt]);
nopt++;
continue;
}
/*----- -input1D dname -----*/
if (strcmp(argv[nopt], "-input1D") == 0)
{
nopt++;
if (nopt >= argc) FDR_error ("need argument after -input1D ");
FDR_input1D_filename = (char *) malloc (sizeof(char)*THD_MAX_NAME);
MTEST (FDR_input1D_filename);
strcpy (FDR_input1D_filename, argv[nopt]);
nopt++;
continue;
}
/*----- -mask_file mname -----*/
if (strcmp(argv[nopt], "-mask_file") == 0 || strcmp(argv[nopt],"-mask") == 0)
{
nopt++;
if (nopt >= argc) FDR_error ("need argument after -mask_file ");
FDR_mask_filename = (char *) malloc (sizeof(char)*THD_MAX_NAME);
MTEST (FDR_mask_filename);
strcpy (FDR_mask_filename, argv[nopt]);
nopt++;
continue;
}
/*----- -mask_thr m -----*/
if( strcmp(argv[nopt],"-mask_thr") == 0 ){
float fval;
nopt++ ;
if( nopt >= argc ){
FDR_error (" need 1 argument after -mask_thr");
}
sscanf (argv[nopt], "%f", &fval);
if (fval < 0.0){
FDR_error (" Require mask_thr >= 0.0 ");
}
FDR_mask_thr = fval;
nopt++; continue;
}
/*---- -force & -old & -pmask & -float etc. [18 Jan 2008] -----*/
if( strcmp(argv[nopt],"-force") == 0 ){
FDR_force = 1 ; nopt++ ; continue ;
}
if( strcmp(argv[nopt],"-qval") == 0 ){
FDR_qval = 1 ; nopt++ ; continue ;
}
if( strcmp(argv[nopt],"-old") == 0 ){
FDR_old = 1 ; FDR_pmask = 0 ; nopt++ ; continue ;
}
if( strcmp(argv[nopt],"-new") == 0 ){
FDR_old = -1 ; FDR_pmask = 1 ; nopt++ ; continue ;
}
if( strcmp(argv[nopt],"-pmask") == 0 ){
FDR_pmask = 1 ; nopt++ ; continue ;
}
if( strcmp(argv[nopt],"-nopmask") == 0 ){
FDR_pmask = 0 ; nopt++ ; continue ;
}
if( strcmp(argv[nopt],"-float") == 0 ){
FDR_float = 1 ; nopt++ ; continue ;
}
#if 0
if( strcmp(argv[nopt],"-curve") == 0 ){ /* hidden option */
FDR_curve = 1 ; nopt++ ; continue ;
}
#endif
/*----- -cind -----*/
if( strcmp(argv[nopt],"-cind") == 0 ){
FDR_cn = 1.0;
nopt++ ; continue ;
}
/*----- -cdep -----*/
if( strcmp(argv[nopt],"-cdep") == 0 ){
FDR_cn = -1.0;
nopt++ ; continue ;
}
/*----- -quiet -----*/
if( strcmp(argv[nopt],"-quiet") == 0 ){
FDR_quiet = 1;
nopt++ ; continue ;
}
if( strcmp(argv[nopt],"-verb") == 0 ){
FDR_quiet = 0 ; nopt++ ; continue ;
}
/*----- -list -----*/
if( strcmp(argv[nopt],"-list") == 0 ){
FDR_list = 1;
nopt++ ; continue ;
}
/*----- -prefix prefix -----*/
if( strcmp(argv[nopt],"-prefix") == 0 ||
strcmp(argv[nopt],"-output") == 0 ){
nopt++ ;
if( nopt >= argc ){
FDR_error (" need argument after -prefix!");
}
FDR_output_prefix = (char *) malloc (sizeof(char) * THD_MAX_PREFIX);
MCW_strncpy( FDR_output_prefix , argv[nopt++] , THD_MAX_PREFIX ) ;
continue ;
}
/*----- unknown command -----*/
sprintf(message,"Unrecognized command line option: %s\n", argv[nopt]);
FDR_error (message);
} /*----- end of loop over command line arguments -----*/
if( FDR_old == 0 && FDR_pmask ){ /* the new way is on by default */
fprintf(stderr,"\n"
"++ The 'new' method of FDR calculation is on by default; in particular:\n"
" + * Voxel p-values of exactly 1 (e.g., from t=0 or F=0 or correlation=0)\n"
" + are ignored by default; in the old mode of operation, they were\n"
" + included in the count which goes into the FDR algorithm. The old\n"
" + process tends to increase the q-values and so decrease the z-scores.\n"
" + * If you wish to do the FDR conversion using the old mode, use '-old'\n"
" + on the command line. For more information, use '3dFDR -help'.\n"
" + * If you don't want to see this message again, use the '-new' option.\n"
"++ RWCox - 18 Jan 2008\n\n"
) ;
}
if( FDR_old < 1 && FDR_pmask == 0 && FDR_mask != NULL ){ /* 29 Jan 2008 */
fprintf(stderr,"\n"
"++ In the 'new' method of FDR calculation, options '-nopmask' and\n"
" + -mask_file are incompatible. Am now turning '-pmask' back on\n"
" + so that the mask can be used.\n"
"++ RWCox - 29 Jan 2008\n\n"
) ;
FDR_pmask = 1 ;
}
return ;
}
/*---------------------------------------------------------------------------*/
/*
Read time series from specified file name. This file name may have
a column selector attached.
*/
float * read_time_series
(
char * ts_filename, /* time series file name (plus column index) */
int * ts_length /* output value for time series length */
)
{
char message[THD_MAX_NAME]; /* error message */
char * cpt; /* pointer to column suffix */
char filename[THD_MAX_NAME]; /* time series file name w/o column index */
char subv[THD_MAX_NAME]; /* string containing column index */
MRI_IMAGE * im, * flim; /* pointers to image structures
-- used to read 1D ASCII */
float * far; /* pointer to MRI_IMAGE floating point data */
int nx; /* number of time points in time series */
int ny; /* number of columns in time series file */
int iy; /* time series file column index */
int ipt; /* time point index */
float * ts_data = NULL; /* input time series data */
/*----- First, check for empty filename -----*/
if (ts_filename == NULL)
FDR_error ("Missing input time series file name");
/*----- Read the time series file -----*/
flim = mri_read_1D(ts_filename) ;
if (flim == NULL)
{
sprintf (message, "Unable to read time series file: %s", ts_filename);
FDR_error (message);
}
far = MRI_FLOAT_PTR(flim);
nx = flim->nx;
ny = flim->ny; iy = 0 ;
if( ny > 1 ){
fprintf(stderr,"WARNING: time series %s has more than 1 column\n",ts_filename);
}
/*----- Save the time series data -----*/
*ts_length = nx;
ts_data = (float *) malloc (sizeof(float) * nx);
MTEST (ts_data);
for (ipt = 0; ipt < nx; ipt++)
ts_data[ipt] = far[ipt + iy*nx];
mri_free (flim); flim = NULL;
return (ts_data);
}
/*---------------------------------------------------------------------------*/
/*
Routine to check whether one output file already exists.
*/
void check_one_output_file
(
THD_3dim_dataset * dset_time, /* input 3d+time data set */
char * filename /* name of output file */
)
{
char message[THD_MAX_NAME]; /* error message */
THD_3dim_dataset * new_dset=NULL; /* output afni data set pointer */
int ierror; /* number of errors in editing data */
ENTRY("check_one_output_file") ;
/*----- make an empty copy of input dataset -----*/
new_dset = EDIT_empty_copy( dset_time ) ;
ierror = EDIT_dset_items( new_dset ,
ADN_prefix , filename ,
ADN_label1 , filename ,
ADN_self_name , filename ,
ADN_type , ISHEAD(dset_time) ? HEAD_FUNC_TYPE :
GEN_FUNC_TYPE ,
ADN_none ) ;
if( ierror > 0 )
{
sprintf (message,
"*** %d errors in attempting to create output dataset!\n",
ierror);
FDR_error (message);
}
if( THD_is_file(new_dset->dblk->diskptr->header_name) )
{
sprintf (message,
"Output dataset file %s already exists "
" -- cannot continue! ",
new_dset->dblk->diskptr->header_name);
FDR_error (message);
}
/*----- deallocate memory -----*/
THD_delete_3dim_dataset( new_dset , False ) ; new_dset = NULL ;
EXRETURN ;
}
/*---------------------------------------------------------------------------*/
/*
Routine to initialize the program: get all operator inputs; create mask
for voxels above mask threshold.
*/
void initialize_program (int argc, char * argv[])
{
int iv; /* index number of sub-brick */
void * vfim = NULL; /* sub-brick data pointer */
float * ffim = NULL; /* sub-brick data in floating point format */
int ixyz; /* voxel index */
int nx, ny, nz, nxyz; /* numbers of voxels in input dataset */
int mx=0, my=0, mz=0, mxyz; /* numbers of voxels in mask dataset */
int nthr=0; /* number of voxels above mask threshold */
char message[80]; /* error message */
int ibin; /* p-value bin index */
/*-- 20 Apr 2001: addto the arglist, if user wants to [RWCox] --*/
machdep() ;
{ int new_argc ; char ** new_argv ;
addto_args( argc , argv , &new_argc , &new_argv ) ;
if( new_argv != NULL ){ argc = new_argc ; argv = new_argv ; }
}
/*----- Save command line for history notes -----*/
commandline = tross_commandline( PROGRAM_NAME , argc,argv ) ;
/*----- Does user request help menu? -----*/
if( argc < 2 || strcmp(argv[1],"-help") == 0 ) FDR_Syntax() ;
/*----- Add to program log -----*/
AFNI_logger (PROGRAM_NAME,argc,argv);
/*----- Read input options -----*/
read_options( argc , argv ) ;
/*----- Open the mask dataset -----*/
if (FDR_mask_filename != NULL)
{
if (!FDR_quiet)
printf ("Reading mask dataset: %s \n", FDR_mask_filename);
DOPEN (FDR_dset, FDR_mask_filename);
if (FDR_dset == NULL)
{
sprintf (message, "Cannot open mask dataset %s", FDR_mask_filename);
FDR_error (message);
}
if (DSET_NVALS(FDR_dset) != 1)
WARNING_message("Mask dataset: using sub-brick #0") ;
/*----- Get dimensions of mask dataset -----*/
mx = DSET_NX(FDR_dset);
my = DSET_NY(FDR_dset);
mz = DSET_NZ(FDR_dset);
mxyz = mx*my*mz;
/*----- Allocate memory for float data -----*/
ffim = (float *) malloc (sizeof(float) * mxyz); MTEST (ffim);
/*----- Convert mask dataset sub-brick to floats (in ffim) -----*/
iv = 0 ;
SUB_POINTER (FDR_dset, iv, 0, vfim);
EDIT_coerce_scale_type (mxyz, DSET_BRICK_FACTOR(FDR_dset,iv),
DSET_BRICK_TYPE(FDR_dset,iv), vfim, /* input */
MRI_float , ffim); /* output */
/*----- Allocate memory for mask volume -----*/
FDR_mask = (byte *) malloc (sizeof(byte) * mxyz);
MTEST (FDR_mask);
/*----- Create mask of voxels above mask threshold -----*/
nthr = 0;
for (ixyz = 0; ixyz < mxyz; ixyz++){
if (fabs(ffim[ixyz]) >= FDR_mask_thr){ FDR_mask[ixyz] = 1; nthr++; }
else FDR_mask[ixyz] = 0;
}
if (!FDR_quiet)
printf ("Number of voxels above mask threshold = %d \n", nthr);
if (nthr < 1)
FDR_error ("No voxels above mask threshold. Cannot continue.");
/*----- Delete floating point sub-brick -----*/
if (ffim != NULL) { free (ffim); ffim = NULL; }
/*----- Delete mask dataset -----*/
THD_delete_3dim_dataset (FDR_dset, False); FDR_dset = NULL ;
}
/*----- Get the input data -----*/
if (FDR_input1D_filename != NULL)
{
/*----- Read the input .1D file -----*/
if (!FDR_quiet) printf ("Reading input data: %s \n",
FDR_input1D_filename);
FDR_input1D_data = read_time_series (FDR_input1D_filename, &nxyz);
if (FDR_input1D_data == NULL)
{
sprintf (message, "Unable to read input .1D data file: %s",
FDR_input1D_filename);
FDR_error (message);
}
if (nxyz < 1)
{
sprintf (message, "No p-values in input .1D data file: %s",
FDR_input1D_filename);
FDR_error (message);
}
FDR_nxyz = nxyz;
FDR_nthr = nxyz;
}
else
{
/*----- Open the input 3D dataset -----*/
if (!FDR_quiet) printf ("Reading input dataset: %s \n",
FDR_input_filename);
FDR_dset = THD_open_dataset(FDR_input_filename);
CHECK_OPEN_ERROR(FDR_dset,FDR_input_filename);
/*----- Get dimensions of input dataset -----*/
nx = DSET_NX(FDR_dset);
ny = DSET_NY(FDR_dset);
nz = DSET_NZ(FDR_dset);
nxyz = nx*ny*nz;
/*----- Check for compatible dimensions -----*/
if (FDR_mask != NULL)
{
if ((nx != mx) || (ny != my) || (nz != mz))
FDR_error ("Mask and input dataset have incompatible dimensions");
FDR_nxyz = nxyz;
FDR_nthr = nthr;
}
else
{
FDR_nxyz = nxyz;
FDR_nthr = nxyz;
}
/*----- Check whether output dataset already exists -----*/
if( THD_deathcon() ) check_one_output_file (FDR_dset, FDR_output_prefix);
}
/*----- Initialize constant c(N) -----*/
if (FDR_cn < 0.0)
{
double cn;
cn = 0.0;
for (ixyz = 1; ixyz <= FDR_nthr; ixyz++)
cn += 1.0 / ixyz;
FDR_cn = cn;
if (!FDR_quiet)
printf ("c(N) = %f \n", FDR_cn);
}
/*----- Initialize voxel pointers -----*/
for (ibin = 0; ibin < FDR_MAX_LL; ibin++)
FDR_head_voxel[ibin] = NULL;
return ;
}
/*---------------------------------------------------------------------------*/
/*
Create an empty voxel.
*/
voxel * create_voxel ()
{
voxel * voxel_ptr = NULL;
voxel_ptr = (voxel *) malloc (sizeof(voxel));
MTEST (voxel_ptr);
voxel_ptr->ixyz = 0;
voxel_ptr->pvalue = 0.0;
voxel_ptr->next_voxel = NULL;
return (voxel_ptr);
}
/*---------------------------------------------------------------------------*/
/*
Add a new voxel to the linked list of voxels.
*/
voxel * add_voxel (voxel * new_voxel, voxel * head_voxel)
{
voxel * voxel_ptr = NULL;
if ((head_voxel == NULL) || (new_voxel->pvalue >= head_voxel->pvalue))
{
new_voxel->next_voxel = head_voxel;
head_voxel = new_voxel;
}
else
{
voxel_ptr = head_voxel;
while ((voxel_ptr->next_voxel != NULL) &&
(new_voxel->pvalue < voxel_ptr->next_voxel->pvalue))
voxel_ptr = voxel_ptr->next_voxel;
new_voxel->next_voxel = voxel_ptr->next_voxel;
voxel_ptr->next_voxel = new_voxel;
}
return (head_voxel);
}
/*---------------------------------------------------------------------------*/
/*
Create and initialize a new voxel, and add to list of voxels.
*/
voxel * new_voxel (int ixyz, float pvalue, voxel * head_voxel)
{
voxel * voxel_ptr = NULL;
voxel_ptr = create_voxel ();
voxel_ptr->ixyz = ixyz;
voxel_ptr->pvalue = pvalue;
head_voxel = add_voxel (voxel_ptr, head_voxel);
return (head_voxel);
}
/*---------------------------------------------------------------------------*/
/*
Deallocate memory for all voxel lists.
*/
void delete_all_voxels ()
{
int ibin;
voxel * voxel_ptr = NULL; /* pointer to current voxel */
voxel * next_voxel = NULL; /* pointer to next voxel */
for (ibin = 0; ibin < FDR_MAX_LL; ibin++)
{
voxel_ptr = FDR_head_voxel[ibin];
while (voxel_ptr != NULL)
{
next_voxel = voxel_ptr->next_voxel;
free (voxel_ptr);
voxel_ptr = next_voxel;
}
FDR_head_voxel[ibin] = NULL;
}
}
/*---------------------------------------------------------------------------*/
/*
Save voxel contents of all voxels into float array (sub-brick).
*/
void save_all_voxels (float * far)
{
int ixyz, ibin;
voxel * voxel_ptr = NULL; /* pointer to voxel */
/*----- Initialize all voxels to zero -----*/
for (ixyz = 0; ixyz < FDR_nxyz; ixyz++)
far[ixyz] = 0.0;
for (ibin = 0; ibin < FDR_MAX_LL; ibin++)
{
voxel_ptr = FDR_head_voxel[ibin];
while (voxel_ptr != NULL)
{
far[voxel_ptr->ixyz] = voxel_ptr->pvalue;
voxel_ptr = voxel_ptr->next_voxel;
}
}
}
/*---------------------------------------------------------------------------*/
/*
Calculate FDR z-scores for all voxels within one volume.
*/
void process_volume (float * ffim, int statcode, float * stataux)
{
int ixyz; /* voxel index */
int icount; /* count of sorted p-values */
float fval; /* voxel input statistical value */
float pval; /* voxel input stat. p-value */
float qval; /* voxel FDR q-value */
float zval; /* voxel FDR z-score */
float qval_min; /* smallest previous q-value */
voxel * head_voxel = NULL; /* linked list of voxels */
voxel * voxel_ptr = NULL; /* pointer to current voxel */
int ibin; /* p-value bin */
int * iarray = NULL; /* output array of voxel indices */
float * parray = NULL; /* output array of voxel p-values */
float * qarray = NULL; /* output array of voxel FDR q-values */
float * zarray = NULL; /* output array of voxel FDR z-scores */
float numer ;
/*------------ 18 Jan 2008: use the 'new' method? ------------*/
if( FDR_old < 1 ){
MRI_IMAGE *qim ; int flags=0 ;
qim = mri_new_vol_empty( FDR_nxyz,1,1 , MRI_float ) ;
mri_fix_data_pointer( ffim , qim ) ;
if( FDR_mask != NULL ){
float zz = (FUNC_IS_STAT(statcode)) ? 0.0f : 1.0f ;
for( ixyz=0 ; ixyz < FDR_nxyz ; ixyz++ )
if( !FDR_mask[ixyz] ) ffim[ixyz] = zz ;
}
if( FDR_curve ){ /* hidden option: produce t-z curve */
floatvec *fv = mri_fdr_curve( qim , statcode , stataux ) ;
if( fv == NULL ) ERROR_message("mri_fdr_curve fails!") ;
else {
printf("# FDR thresh-z curve\n") ;
for( ixyz=0 ; ixyz < fv->nar ; ixyz++ )
printf("%g %g\n", fv->x0+ixyz*fv->dx , fv->ar[ixyz] ) ;
}
exit(0) ;
} else { /* normal operation: convert to z(q) or q */
if( FDR_pmask == 0 ) flags |= 1 ; /* compatibility mode */
if( FDR_cn > 1.0f ) flags |= 2 ; /* dependency flag */
if( FDR_qval ) flags |= 4 ; /* qval flag */
(void)mri_fdrize( qim , statcode,stataux , flags ) ;
}
mri_clear_data_pointer(qim); mri_free(qim);
return ;
}
/*---------------- back to the 'old' method ------------------*/