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