1 /*****************************************************************************
2 Major portions of this software are copyrighted by the Medical College
3 of Wisconsin, 1994-2001, and are released under the Gnu General Public
4 License, Version 2. See the file README.Copyright for details.
5 ******************************************************************************/
6
7 #ifdef USE_SUNPERF /** for Solaris **/
8 # include <sunperf.h>
9 #endif
10
11 /*
12 This program calculates a nonlinear regression for each voxel of the input
13 AFNI 3d+time data set. The nonlinear regression is calculated by means of
14 a least squares fit to the signal plus noise models which are specified
15 by the user.
16
17 File: 3dNLfim.c
18 Author: B. Douglas Ward
19 Date: 19 June 1997
20
21 Mod: Added external time reference capability (Rongyan Zhang)
22 Date: 10 August 1997
23
24 Mod: Added option for absolute noise parameter constraints.
25 Date: 22 August 1997
26
27 Mod: Added options for percent signal change above baseline, and
28 calculation of area under the signal above baseline.
29 Date: 26 November 1997
30
31 Mod: Print error message if user fails to specify the signal model or
32 the noise model.
33 Date: 26 December 1997
34
35 Mod: Extensive changes required to implement the 'bucket' dataset.
36 Date: 14 January 1998
37
38 Mod: Added the -inTR option.
39 22 July 1998 -- RWCox
40
41 Mod: Incorporated THD_extract_series routine.
42 Date: 19 April 1999
43
44 Mod: Added -sfit and -snfit options to write out the signal and
45 the signal+noise model time series fit for each voxel
46 to a 3d+time dataset.
47 Date: 08 July 1999
48
49 Mod: Added novar flag to eliminate unnecessary calculations.
50 Date: 13 July 1999
51
52 Mod: Added changes for incorporating History notes.
53 Date: 09 September 1999
54
55 Mod: Adjust F-statistics if parameter constraints force a parameter
56 to be a constant.
57 Date: 08 February 2000
58
59 Mod: Changes for output of R^2 (coefficient of multiple determination),
60 and standard deviation of residuals from full model fit.
61 Added global variable calc_tstats.
62 Also, added screen display of p-values.
63 Date: 10 May 2000
64
65 Mod: Corrected error in initializing of output data type (MRI_short)
66 in routine write_3dtime.
67 Date: 17 May 2000
68
69 Mod: Added -mask option. (Adapted from: 3dpc.c)
70 Date: 18 May 2000
71
72 Mod: Changed "return" at end of program to exit(0). Also, increased
73 maximum number of model parameters.
74 Date: 08 August 2001
75
76 Mod: Added call to AFNI_logger.
77 Date: 15 August 2001
78
79 Mod: Changes to allow -jobs option.
80 Date: 07 May 2003 - RWCox.
81
82 Mod: Added options -aux_name, -aux_fval and -voxel_count.
83 Date: 25 Jan 2006 [rickr]
84
85 Mod: Removed options -aux_name and -aux_fval, and the globals
86 require linking to afni, too.
87 Date: 30 Jan 2006 [rickr]
88
89 Mod: Added NEWUOA stuff (mostly to simplex.c, actually).
90 Date: 20 Jul 2006 [RWCox]
91
92 Mod: Copied memmap code from 3dDeconvolve.c
93 Date: 24 Oct 2006 [DRG]
94
95 Mod: Limit reports to nth voxels via progress option
96 Date: 25 Oct 2006 [DRG]
97
98 Mod: Limit g_voxel_count reports to every 10th voxel.
99 Date: 26 Oct 2006 [rickr]
100 */
101
102 /*---------------------------------------------------------------------------*/
103
104 #define PROGRAM_NAME "3dNLfim" /* name of this program */
105 #define PROGRAM_AUTHOR "B. Douglas Ward" /* program author */
106 #define PROGRAM_INITIAL "19 June 1997" /* date of initial program release */
107 #define PROGRAM_LATEST "24 Oct 2006 - DRG" /* date of latest program revision */
108
109 /*---------------------------------------------------------------------------*/
110
111 #define DEFAULT_NRAND 19999
112 #define DEFAULT_NBEST 9
113 #define DEFAULT_FDISP 999.0
114 #define DEFAULT_PROGRESS 10000
115
116 #include <stdio.h>
117 #include <math.h>
118 #include <stdlib.h>
119
120 #include "mrilib.h"
121
122 static int do_FDR = 1 ; /* 07 Oct 2008 */
123
124 /*---------------------------------------------------------------------------*/
125 /*--------- Global variables for multiple process execution - RWCox. --------*/
126 /*--------- All names start with "proc_", so search for that string. --------*/
127
128 #if !defined(DONT_USE_SHM) && !defined(DONT_USE_FORK) && !defined(CYGWIN)
129
130 # include "thd_iochan.h" /* prototypes for shm stuff */
131
132 # define PROC_MAX 32 /* max num processes */
133
134 static int proc_numjob = 1 ; /* num processes */
135 static pid_t proc_pid[PROC_MAX] ; /* IDs of processes */
136 static int proc_shmid = 0 ; /* shared memory ID */
137 static char *proc_shmptr = NULL; /* pointer to shared memory */
138 static long long proc_shmsize = 0 ; /* total size of shared memory */
139
140 static int proc_shm_arnum = 0 ; /* num arrays in shared memory */
141 static float ***proc_shm_ar = NULL; /* *proc_shm_ar[i] = ptr to array #i */
142 static int *proc_shm_arsiz = NULL; /* proc_shm_arsiz[i] = floats in #i */
143
144 static int proc_vox_bot[PROC_MAX] ; /* 1st voxel to use in each process */
145 static int proc_vox_top[PROC_MAX] ; /* last voxel (+1) in each process */
146
147 static int proc_ind = 0 ; /* index of THIS job */
148
149 #else /* can't use multiple processes */
150
151 # define proc_numjob 1 /* flag that only 1 process is in use */
152 # define proc_ind 0 /* index of THIS job */
153
154 #endif
155
156 /*---------------------------------------------------------------------------*/
157
158 #include "matrix.h"
159 #include "simplex.h"
160 #include "NLfit.h"
161
162 #include "simplex.c" /* 20 Jul 2006 - now includes NEWUOA variables [N_*] */
163 #include "NLfit.c"
164
165 typedef struct NL_options
166 {
167 char * bucket_filename; /* file name for bucket dataset */
168 int numbricks; /* number of sub-bricks in bucket dataset */
169 int * brick_type; /* indicates type of sub-brick */
170 int * brick_coef; /* regression coefficient number for sub-brick*/
171 char ** brick_label; /* character string label for sub-brick */
172
173 } NL_options;
174
175
176 /*---------------------------------------------------------------------------*/
177 /*
178 Global data
179 */
180
181 /***** 22 July 1998 -- RWCox:
182 Modified to allow DELT to be set from the TR of the input file *****/
183
184 static float DELT = 1.0; /* default */
185 static int inTR = 0 ; /* set to 1 if -inTR option is used */
186 static float dsTR = 0.0 ; /* TR of the input file */
187 static float uuTR = 0.0 ;
188 static int g_voxel_count = 0; /* display current voxel counter */
189 /* 25 Jan 2006 [rickr] */
190
191 static char *commandline = NULL ; /* command line for history notes */
192
193 static byte *mask_vol = NULL; /* mask volume */
194 static int mask_nvox = 0; /* number of voxels in mask volume */
195 static int output_datum = ILLEGAL_TYPE ;
196 static byte *gmask = NULL ;
197
198 /*---------------------------------------------------------------------------*/
199 /*
200 Routine to display 3dNLfim help menu.
201 */
202
display_help_menu()203 void display_help_menu()
204 {
205 printf
206 (
207 "This program calculates a nonlinear regression for each voxel of the \n"
208 "input AFNI 3d+time data set. The nonlinear regression is calculated \n"
209 "by means of a least squares fit to the signal plus noise models which \n"
210 "are specified by the user. \n"
211 "\n"
212 "Usage with terminal options:\n"
213 "\n"
214 "3dNLfim\n"
215 " -help show this help\n"
216 " -help_models show model help from any that have it\n"
217 " (can come via AFNI_MODEL_HELP_ALL)\n"
218 "\n"
219 " One can get help for an individual model, *if* it exists, by\n"
220 " setting a similar environment variable, and providing some\n"
221 " non-trivial function option (like -load_models), e.g.,\n"
222 "\n"
223 " 3dNLfim -DAFNI_MODEL_HELP_CONV_PRF_6=Y -load_models\n"
224 "\n"
225 " Indifidual help should be available for any model with help\n"
226 " via -help_models.\n"
227 "\n"
228 " -load_models simply load all models and exit\n"
229 " (this is for testing or getting model help)\n"
230 "\n"
231 "General usage:\n"
232 "\n"
233 "3dNLfim\n"
234 "-input fname fname = filename of 3d + time data file for input \n"
235 "[-mask mset] Use the 0 sub-brick of dataset 'mset' as a mask \n"
236 " to indicate which voxels to analyze (a sub-brick \n"
237 " selector is allowed) [default = use all voxels] \n"
238 "[-ignore num] num = skip this number of initial images in the \n"
239 " time series for regresion analysis; default = 0 \n"
240 " ****N.B.: default ignore value changed from 3 to 0, \n"
241 " on 04 Nov 2008 (BHO day). \n"
242 "[-inTR] set delt = TR of the input 3d+time dataset \n"
243 " [The default is to compute with delt = 1.0 ] \n"
244 " [The model functions are calculated using a \n"
245 " time grid of: 0, delt, 2*delt, 3*delt, ... ] \n"
246 "[-TR delt] directly set the TR of the time series model; \n"
247 " can be useful if the input file is a .1D file \n"
248 " (transposed with the \\' operator) \n"
249 "[-time fname] fname = ASCII file containing each time point \n"
250 " in the time series. Defaults to even spacing \n"
251 " given by TR (this option overrides -inTR). \n"
252 "-signal slabel slabel = name of (non-linear) signal model \n"
253 "-noise nlabel nlabel = name of (linear) noise model \n"
254 "-sconstr k c d constraints for kth signal parameter: \n"
255 " c <= gs[k] <= d \n"
256 " **N.B.: It is important to set the parameter \n"
257 " constraints with care! \n"
258 " **N.B.: -sconstr and -nconstr options must appear \n"
259 " AFTER -signal and -noise on the command line \n"
260 "-nconstr k c d constraints for kth noise parameter: \n"
261 " c+b[k] <= gn[k] <= d+b[k] \n"
262 "[-nabs] use absolute constraints for noise parameters: \n"
263 " c <= gn[k] <= d [default=relative, as above] \n"
264 "[-nrand n] n = number of random test points [default=%d] \n"
265 "[-nbest b] b = use b best test points to start [default=%d] \n"
266 "[-rmsmin r] r = minimum rms error to reject reduced model \n"
267 "[-fdisp fval] display (to screen) results for those voxels \n"
268 " whose f-statistic is > fval [default=%.1f] \n"
269 "[-progress ival] display (to screen) results for those voxels \n"
270 " every ival number of voxels \n"
271 "[-voxel_count] display (to screen) the current voxel index \n"
272 " \n"
273 "--- These options choose the least-square minimization algorithm --- \n"
274 " \n"
275 "[-SIMPLEX] use Nelder-Mead simplex method [default] \n"
276 "[-POWELL] use Powell's NEWUOA method instead of the \n"
277 " Nelder-Mead simplex method to find the \n"
278 " nonlinear least-squares solution \n"
279 " [slower; usually more accurate, but not always!] \n"
280 "[-BOTH] use both Powell's and Nelder-Mead method \n"
281 " [slowest, but should be most accurate] \n"
282 " \n"
283 "--- These options generate individual AFNI 2 sub-brick datasets --- \n"
284 "--- [All these options must be AFTER options -signal and -noise]--- \n"
285 " \n"
286 "[-freg fname] perform f-test for significance of the regression; \n"
287 " output 'fift' is written to prefix filename fname\n"
288 "[-frsqr fname] calculate R^2 (coef. of multiple determination); \n"
289 " store along with f-test for regression; \n"
290 " output 'fift' is written to prefix filename fname\n"
291 "[-fsmax fname] estimate signed maximum of signal; store along \n"
292 " with f-test for regression; output 'fift' is \n"
293 " written to prefix filename fname \n"
294 "[-ftmax fname] estimate time of signed maximum; store along \n"
295 " with f-test for regression; output 'fift' is \n"
296 " written to prefix filename fname \n"
297 "[-fpsmax fname] calculate (signed) maximum percentage change of \n"
298 " signal from baseline; output 'fift' is \n"
299 " written to prefix filename fname \n"
300 "[-farea fname] calculate area between signal and baseline; store \n"
301 " with f-test for regression; output 'fift' is \n"
302 " written to prefix filename fname \n"
303 "[-fparea fname] percentage area of signal relative to baseline; \n"
304 " store with f-test for regression; output 'fift' \n"
305 " is written to prefix filename fname \n"
306 "[-fscoef k fname] estimate kth signal parameter gs[k]; store along \n"
307 " with f-test for regression; output 'fift' is \n"
308 " written to prefix filename fname \n"
309 "[-fncoef k fname] estimate kth noise parameter gn[k]; store along \n"
310 " with f-test for regression; output 'fift' is \n"
311 " written to prefix filename fname \n"
312 "[-tscoef k fname] perform t-test for significance of the kth signal \n"
313 " parameter gs[k]; output 'fitt' is written \n"
314 " to prefix filename fname \n"
315 "[-tncoef k fname] perform t-test for significance of the kth noise \n"
316 " parameter gn[k]; output 'fitt' is written \n"
317 " to prefix filename fname \n"
318 " \n"
319 "--- These options generate one AFNI 'bucket' type dataset --- \n"
320 " \n"
321 "[-bucket n prefixname] create one AFNI 'bucket' dataset containing \n"
322 " n sub-bricks; n=0 creates default output; \n"
323 " output 'bucket' is written to prefixname \n"
324 "The mth sub-brick will contain: \n"
325 "[-brick m scoef k label] kth signal parameter regression coefficient\n"
326 "[-brick m ncoef k label] kth noise parameter regression coefficient \n"
327 "[-brick m tmax label] time at max. abs. value of signal \n"
328 "[-brick m smax label] signed max. value of signal \n"
329 "[-brick m psmax label] signed max. value of signal as percent \n"
330 " above baseline level \n"
331 "[-brick m area label] area between signal and baseline \n"
332 "[-brick m parea label] signed area between signal and baseline \n"
333 " as percent of baseline area \n"
334 "[-brick m tscoef k label] t-stat for kth signal parameter coefficient\n"
335 "[-brick m tncoef k label] t-stat for kth noise parameter coefficient \n"
336 "[-brick m resid label] std. dev. of the full model fit residuals \n"
337 "[-brick m rsqr label] R^2 (coefficient of multiple determination)\n"
338 "[-brick m fstat label] F-stat for significance of the regression \n"
339 "\n"
340 "[-noFDR] Don't write the FDR (q vs. threshold)\n"
341 " curves into the output dataset.\n"
342 " (Same as 'setenv AFNI_AUTOMATIC_FDR NO')\n"
343 " \n"
344 " --- These options write time series fit for --- \n"
345 " --- each voxel to an AFNI 3d+time dataset --- \n"
346 " \n"
347 "[-sfit fname] fname = prefix for output 3d+time signal model fit \n"
348 "[-snfit fname] fname = prefix for output 3d+time signal+noise fit \n"
349 " \n"
350 , DEFAULT_NRAND , DEFAULT_NBEST , DEFAULT_FDISP
351 );
352
353
354 #ifdef PROC_MAX
355 printf( "\n"
356 " -jobs J Run the program with 'J' jobs (sub-processes).\n"
357 " On a multi-CPU machine, this can speed the\n"
358 " program up considerably. On a single CPU\n"
359 " machine, using this option is silly.\n"
360 " J should be a number from 1 up to the\n"
361 " number of CPU sharing memory on the system.\n"
362 " J=1 is normal (single process) operation.\n"
363 " The maximum allowed value of J is %d.\n"
364 " * For more information on parallelizing, see\n"
365 " https://sscc.nimh.nih.gov/afni/doc/misc/afni_parallelize/index_html/view\n"
366 " * Use -mask to get more speed; cf. 3dAutomask.\n"
367 , PROC_MAX ) ;
368 #endif
369
370 printf(
371 "\n"
372 "----------------------------------------------------------------------\n"
373 "Signal Models (see the appropriate model_*.c file for exact details) :\n"
374 "\n"
375 " Null : No Signal\n"
376 " (no parameters)\n"
377 " see model_null.c\n"
378 "\n"
379 " SineWave_AP : Sinusoidal Response\n"
380 " (amplitude, phase)\n"
381 " see model_sinewave_ap.c\n"
382 "\n"
383 " SquareWave_AP : Square Wave Response\n"
384 " (amplitude, phase)\n"
385 " see model_squarewave_ap.c\n"
386 "\n"
387 " TrnglWave_AP : Triangular Wave Response\n"
388 " (amplitude, phase)\n"
389 " see model_trnglwave_ap.c\n"
390 "\n"
391 " SineWave_APF : Sinusoidal Wave Response\n"
392 " (amplitude, phase, frequency)\n"
393 " see model_sinewave_apf.c\n"
394 "\n"
395 " SquareWave_APF : Sinusoidal Wave Response\n"
396 " (amplitude, phase, frequency)\n"
397 " see model_squarewave_apf.c\n"
398 "\n"
399 " TrnglWave_APF : Sinusoidal Wave Response\n"
400 " (amplitude, phase, frequency)\n"
401 " see model_trnglwave_apf.c\n"
402 "\n"
403 " Exp : Exponential Function\n"
404 " (a,b): a * exp(b * t)\n"
405 " see model_exp.c\n"
406 "\n"
407 " DiffExp : Differential-Exponential Drug Response\n"
408 " (t0, k, alpha1, alpha2)\n"
409 " see model_diffexp.c\n"
410 "\n"
411 " GammaVar : Gamma-Variate Function Drug Response\n"
412 " (t0, k, r, b)\n"
413 " see model_gammavar.c\n"
414 "\n"
415 " Beta : Beta Distribution Model\n"
416 " (t0, tf, k, alpha, beta)\n"
417 " see model_beta.c\n"
418 "\n"
419 "\n"
420 " * The following convolved functions are generally convolved with\n"
421 " the time series in AFNI_CONVMODEL_REF, allowing one to specify\n"
422 " multiple event onsets, varying durations and varying resopnse\n"
423 " magnitudes.\n"
424 "\n"
425 " ConvGamma : Gamma Vairate Response Model\n"
426 " (t0, amp, r, b)\n"
427 " see model_convgamma.c\n"
428 "\n"
429 " ConvGamma2a : Gamma Convolution with 2 Input Time Series\n"
430 " (t0, r, b)\n"
431 " see model_convgamma2a.c\n"
432 "\n"
433 " ConvDiffGam : Difference of 2 Gamma Variates\n"
434 " (A0, T0, E0, D0, A1, T1, E1, D1)\n"
435 " see model_conv_diffgamma.c\n"
436 " for help : setenv AFNI_MODEL_HELP_CONVDIFFGAM YES\n"
437 " 3dNLfim -signal ConvDiffGam\n"
438 "\n"
439 " demri_3 : Dynamic (contrast) Enhanced MRI\n"
440 " (K_trans, Ve, k_ep)\n"
441 " see model_demri_3.c\n"
442 " for help : setenv AFNI_MODEL_HELP_DEMRI_3 YES\n"
443 " 3dNLfim -signal demri_3\n"
444 "\n"
445 " ADC : Diffusion Signal Model\n"
446 " (So, D)\n"
447 " see model_diffusion.c\n"
448 "\n"
449 " michaelis_menton : Michaelis/Menten Concentration Model\n"
450 " (v, vmax, k12, k21, mag)\n"
451 " see model_michaelis_menton.c\n"
452 "\n"
453 " Expr2 : generic (3dcalc-like) expression with\n"
454 " exactly 2 'free' parameters and using\n"
455 " symbol 't' as the time variable;\n"
456 " see model_expr2.c for details.\n"
457 "\n"
458 " ConvCosine4 : 4-piece Cosine Convolution Model\n"
459 " (A, C1, C2, M1, M2, M3, M4)\n"
460 " see model_conv_cosine4.c\n"
461 " for help : setenv AFNI_MODEL_HELP_CONV_COSINE4 YES\n"
462 " 3dNLfim -signal ConvCosine4\n"
463 "\n"
464 " Conv_PRF : 4-param Population Receptive Field Model\n"
465 " (A, X, Y, sigma)\n"
466 " see model_conv_PRF.c\n"
467 " for help : setenv AFNI_MODEL_HELP_CONV_PRF YES\n"
468 " 3dNLfim -signal bunnies\n"
469 "\n"
470 " Conv_PRF_6 : 6-param Population Receptive Field Model\n"
471 " (A, X, Y, sigma, sigrat, theta)\n"
472 " see model_conv_PRF_6.c\n"
473 " for help : setenv AFNI_MODEL_HELP_CONV_PRF_6 YES\n"
474 " 3dNLfim -signal bunnies\n"
475 "\n"
476 " Conv_PRF_DOG : 6-param 'Difference of Gaussians' PRF Model\n"
477 " (as Conv_PRF, but with second A and sigma)\n"
478 " (A, X, Y, sig, A2, sig2)\n"
479 " see model_conv_PRF_DOG.c\n"
480 " for help : setenv AFNI_MODEL_HELP_CONV_PRF_DOG YES\n"
481 " 3dNLfim -signal bunnies\n"
482 "\n"
483 "----------------------------------------\n"
484 "Noise Models (see the appropriate model_*.c file for exact details) :\n"
485 "\n"
486 " Zero : Zero Noise Model\n"
487 " (no parameters)\n"
488 " see model_zero.c\n"
489 "\n"
490 " Constant : Constant Noise Model\n"
491 " (constant)\n"
492 " see model_constant.c\n"
493 "\n"
494 " Linear : Linear Noise Model\n"
495 " (constant, linear)\n"
496 " see model_linear.c\n"
497 "\n"
498 " Linear+Ort : Linear+Ort Noise Model\n"
499 " (constant, linear, Ort)\n"
500 " see model_linplusort.c\n"
501 "\n"
502 " Quadratic : Quadratic Noise Model\n"
503 " (constant, linear, quadratic)\n"
504 " see model_quadratic.c\n"
505 ) ;
506 PRINT_COMPILE_DATE ; exit(0);
507 }
508
509
510 /*---------------------------------------------------------------------------*/
511
512 /** macro to test a malloc-ed pointer for validity **/
513
514 #define MTEST(ptr) \
515 if((ptr)==NULL) \
516 ( NLfit_error ("Cannot allocate memory") )
517
518
519 /*---------------------------------------------------------------------------*/
520 /*
521 Routine to initialize the input options.
522 */
523
initialize_options(int * ignore,vfp * nmodel,vfp * smodel,int * r,int * p,char *** npname,char *** spname,float ** min_nconstr,float ** max_nconstr,float ** min_sconstr,float ** max_sconstr,int * nabs,int * nrand,int * nbest,float * rms_min,float * fdisp,int * progress,char ** input_filename,char ** tfilename,char ** freg_filename,char ** frsqr_filename,char *** fncoef_filename,char *** fscoef_filename,char *** tncoef_filename,char *** tscoef_filename,char ** sfit_filename,char ** snfit_filename,NL_options * option_data)524 void initialize_options
525 (
526 int * ignore, /* delete this number of points from time series */
527 vfp * nmodel, /* pointer to noise model */
528 vfp * smodel, /* pointer to signal model */
529 int * r, /* number of parameters in the noise model */
530 int * p, /* number of parameters in the signal model */
531 char *** npname, /* noise parameter names */
532 char *** spname, /* signal parameter names */
533 float ** min_nconstr, /* minimum parameter constraints for noise model */
534 float ** max_nconstr, /* maximum parameter constraints for noise model */
535 float ** min_sconstr, /* minimum parameter constraints for signal model */
536 float ** max_sconstr, /* maximum parameter constraints for signal model */
537 int * nabs, /* use absolute constraints for noise parameters */
538 int * nrand, /* number of random vectors to generate */
539 int * nbest, /* number of random vectors to keep */
540 float * rms_min, /* minimum rms error to reject reduced model */
541 float * fdisp, /* minimum f-statistic for display */
542 int *progress, /* nth voxel to show report */
543 char ** input_filename, /* file name of input 3d+time dataset */
544 char ** tfilename, /* file name for time point series */
545 char ** freg_filename, /* file name for regression f-statistics */
546 char ** frsqr_filename, /* file name for R^2 statistics */
547 char *** fncoef_filename, /* file name for noise model parameters */
548 char *** fscoef_filename, /* file name for signal model parameters */
549 char *** tncoef_filename, /* file name for noise model t-statistics */
550 char *** tscoef_filename, /* file name for signal model t-statistics */
551 char ** sfit_filename, /* file name for 3d+time fitted signal model */
552 char ** snfit_filename, /* file name for 3d+time fitted signal+noise */
553 NL_options * option_data /* bucket dataset options */
554 )
555
556 {
557 int ip; /* parameter index */
558
559
560 /*----- initialize default values -----*/
561 *ignore = 0;
562 *nabs = 0;
563 *nrand = DEFAULT_NRAND;
564 *nbest = DEFAULT_NBEST;
565 *rms_min = 0.0;
566 *fdisp = DEFAULT_FDISP;
567 *progress = DEFAULT_PROGRESS;
568 *smodel = NULL;
569 *nmodel = NULL;
570 *r = -1;
571 *p = -1;
572
573
574 /*----- allocate memory for noise parameter names -----*/
575 *npname = (char **) malloc (sizeof(char *) * MAX_PARAMETERS);
576 MTEST (*npname);
577 for (ip = 0; ip < MAX_PARAMETERS; ip++)
578 {
579 (*npname)[ip] = (char *) malloc (sizeof(char) * MAX_NAME_LENGTH);
580 MTEST ((*npname)[ip]);
581 }
582
583
584 /*----- allocate memory for signal parameter names -----*/
585 *spname = (char **) malloc (sizeof(char *) * MAX_PARAMETERS);
586 MTEST (*spname);
587 for (ip = 0; ip < MAX_PARAMETERS; ip++)
588 {
589 (*spname)[ip] = (char *) malloc (sizeof(char) * MAX_NAME_LENGTH);
590 MTEST ((*spname)[ip]);
591 }
592
593
594 /*----- allocate memory for parameter constraints -----*/
595 *min_nconstr = (float *) malloc (sizeof(float) * MAX_PARAMETERS);
596 MTEST (*min_nconstr);
597 *max_nconstr = (float *) malloc (sizeof(float) * MAX_PARAMETERS);
598 MTEST (*max_nconstr);
599 *min_sconstr = (float *) malloc (sizeof(float) * MAX_PARAMETERS);
600 MTEST (*min_sconstr);
601 *max_sconstr = (float *) malloc (sizeof(float) * MAX_PARAMETERS);
602 MTEST (*max_sconstr);
603
604
605 /*----- allocate memory space and initialize pointers for filenames -----*/
606 *input_filename = NULL;
607 *tfilename = NULL;
608 *freg_filename = NULL;
609 *frsqr_filename = NULL;
610 *fncoef_filename = (char **) malloc (sizeof(char *) * MAX_PARAMETERS);
611 MTEST (*fncoef_filename);
612 *fscoef_filename = (char **) malloc (sizeof(char *) * MAX_PARAMETERS);
613 MTEST (*fscoef_filename);
614 *tncoef_filename = (char **) malloc (sizeof(char *) * MAX_PARAMETERS);
615 MTEST (*tncoef_filename);
616 *tscoef_filename = (char **) malloc (sizeof(char *) * MAX_PARAMETERS);
617 MTEST (*tscoef_filename);
618 for (ip = 0; ip < MAX_PARAMETERS; ip++)
619 {
620 (*fncoef_filename)[ip] = NULL;
621 (*fscoef_filename)[ip] = NULL;
622 (*tncoef_filename)[ip] = NULL;
623 (*tscoef_filename)[ip] = NULL;
624 }
625 *sfit_filename = NULL;
626 *snfit_filename = NULL;
627
628
629 /*----- initialize bucket dataset options -----*/
630 option_data->bucket_filename = NULL;
631 option_data->numbricks = -1;
632 option_data->brick_type = NULL;
633 option_data->brick_coef = NULL;
634 option_data->brick_label = NULL;
635
636 }
637
638
639 /*---------------------------------------------------------------------------*/
640 /*
641 Routine to get user specified input options.
642 */
643
get_options(int argc,char ** argv,int * ignore,char ** nname,char ** sname,vfp * nmodel,vfp * smodel,int * r,int * p,char *** npname,char *** spname,float ** min_nconstr,float ** max_nconstr,float ** min_sconstr,float ** max_sconstr,int * nabs,int * nrand,int * nbest,float * rms_min,float * fdisp,int * progress,char ** input_filename,char ** tfilename,char ** freg_filename,char ** frsqr_filename,char ** fsmax_filename,char ** ftmax_filename,char ** fpmax_filename,char ** farea_filename,char ** fparea_filename,char *** fncoef_filename,char *** fscoef_filename,char *** tncoef_filename,char *** tscoef_filename,char ** sfit_filename,char ** snfit_filename,THD_3dim_dataset ** dset_time,int * nxyz,int * ts_length,NL_options * option_data)644 void get_options
645 (
646 int argc, /* number of input arguments */
647 char ** argv, /* array of input arguments */
648 int * ignore, /* delete this number of points from time series */
649 char ** nname, /* name of noise model */
650 char ** sname, /* name of signal model */
651 vfp * nmodel, /* pointer to noise model */
652 vfp * smodel, /* pointer to signal model */
653 int * r, /* number of parameters in the noise model */
654 int * p, /* number of parameters in the signal model */
655 char *** npname, /* noise parameter names */
656 char *** spname, /* signal parameter names */
657 float ** min_nconstr, /* minimum parameter constraints for noise model */
658 float ** max_nconstr, /* maximum parameter constraints for noise model */
659 float ** min_sconstr, /* minimum parameter constraints for signal model */
660 float ** max_sconstr, /* maximum parameter constraints for signal model */
661 int * nabs, /* use absolute constraints for noise parameters */
662 int * nrand, /* number of random vectors to generate */
663 int * nbest, /* number of random vectors to keep */
664 float * rms_min, /* minimum rms error to reject reduced model */
665 float * fdisp, /* minimum f-statistic for display */
666 int * progress, /* progress report every nth voxel */
667 char ** input_filename, /* file name of input 3d+time dataset */
668 char ** tfilename, /* file name for time point series */
669 char ** freg_filename, /* file name for regression f-statistics */
670 char ** frsqr_filename, /* file name for R^2 statistics */
671 char ** fsmax_filename, /* file name for signal signed maximum */
672 char ** ftmax_filename, /* file name for time of signed maximum */
673 char ** fpmax_filename, /* file name for max. percentage change */
674 char ** farea_filename, /* file name for area under the signal */
675 char ** fparea_filename, /* file name for percent area under the signal */
676 char *** fncoef_filename, /* file name for noise model parameters */
677 char *** fscoef_filename, /* file name for signal model parameters */
678 char *** tncoef_filename, /* file name for noise model t-statistics */
679 char *** tscoef_filename, /* file name for signal model t-statistics */
680 char ** sfit_filename, /* file name for 3d+time fitted signal model */
681 char ** snfit_filename, /* file name for 3d+time fitted signal+noise */
682
683 THD_3dim_dataset ** dset_time, /* input 3d+time data set */
684 int * nxyz, /* number of voxels in image */
685 int * ts_length, /* length of time series data */
686 NL_options * option_data /* bucket dataset options */
687 )
688
689 #define MAX_BRICKS 100
690 {
691 int nopt = 1; /* input option argument counter */
692 int ival, index; /* integer input */
693 float fval; /* float input */
694 char message[MAX_NAME_LENGTH]; /* error message */
695 int ok; /* boolean for specified model exists */
696
697 NLFIT_MODEL_array * model_array = NULL; /* array of SO models */
698 int im; /* model index */
699 int ibrick; /* sub-brick index */
700 int nbricks; /* number of bricks in the bucket */
701
702
703 /*----- does user request help menu? -----*/
704 if (argc < 2 || strcmp(argv[1], "-help") == 0) display_help_menu();
705
706
707 /*----- add to program log -----*/
708 AFNI_logger (PROGRAM_NAME,argc,argv);
709
710 /*----- to get help on all models, set general env var, load models -----*/
711 /* and bail 17 May 2018 [rickr] */
712 if ( !strcmp(argv[1], "-help_models") || !strcmp(argv[1], "-load_models") ) {
713 if ( !strcmp(argv[1], "-help_models") )
714 AFNI_setenv("AFNI_MODEL_HELP_ALL YES");
715 model_array = NLFIT_get_many_MODELs ();
716 exit(0);
717 }
718
719 /*----- initialize the model array -----*/
720 model_array = NLFIT_get_many_MODELs ();
721 if ((model_array == NULL) || (model_array->num == 0))
722 NLfit_error ("Unable to locate any models");
723
724 /*----- initialize the input options -----*/
725 initialize_options (ignore, nmodel, smodel, r, p, npname, spname,
726 min_nconstr, max_nconstr, min_sconstr, max_sconstr, nabs,
727 nrand, nbest, rms_min, fdisp, progress,
728 input_filename, tfilename, freg_filename,
729 frsqr_filename, fncoef_filename, fscoef_filename,
730 tncoef_filename, tscoef_filename,
731 sfit_filename, snfit_filename, option_data);
732
733 if( *ignore == 0 )
734 INFO_message("NOTICE: Default value of -ignore is now 0;\n"
735 " Before 04 Nov 2008, default value was 3" ) ;
736
737 /*----- main loop over input options -----*/
738 *input_filename = NULL;
739 while (nopt < argc )
740 {
741 /*----- -noFDR -----*/
742
743 if( strcasecmp(argv[nopt],"-noFDR") == 0 ){
744 do_FDR = 0 ; nopt++ ; continue ;
745 }
746
747 /*----- -input filename -----*/
748 if (strcmp(argv[nopt], "-input") == 0)
749 {
750 nopt++;
751 if (nopt >= argc) NLfit_error ("need argument after -input ");
752 /* *input_filename = malloc (sizeof(char) * MAX_NAME_LENGTH); */
753 /* do not depend on input name lengths [2 Oct 2019 rickr] */
754 *input_filename = nifti_strdup(argv[nopt]);
755 MTEST (*input_filename);
756 /* strcpy (*input_filename, argv[nopt]); */
757 /* open input dataset - was THD_open_one_dataset, allow sub-bricks */
758 *dset_time = THD_open_dataset (*input_filename);
759 if ((*dset_time) == NULL)
760 {
761 sprintf (message,
762 "Unable to open data file: %s", *input_filename);
763 NLfit_error (message);
764 }
765 DSET_UNMSEC( *dset_time ) ; /* RWCox */
766
767 DSET_load(*dset_time) ; CHECK_LOAD_ERROR(*dset_time) ;
768
769 *nxyz = (*dset_time)->dblk->diskptr->dimsizes[0]
770 * (*dset_time)->dblk->diskptr->dimsizes[1]
771 * (*dset_time)->dblk->diskptr->dimsizes[2] ;
772 #if 0
773 *ts_length = DSET_NUM_TIMES(*dset_time);
774
775 /* verify that we seem to have a time series */
776 if( DSET_NVALS(*dset_time) != *ts_length )
777 WARNING_message("dataset num_times (%d) != nvals (%d)\n"
778 " --> no time axis could be a problem!\n",
779 *ts_length, DSET_NVALS(*dset_time));
780 #else
781 *ts_length = DSET_NVALS(*dset_time);
782 if( *ts_length > 1 )
783 INFO_message("time series length = %d points",*ts_length) ;
784 else
785 ERROR_exit ("time series length = %d points",*ts_length) ;
786 #endif
787
788 dsTR = DSET_TIMESTEP(*dset_time) ;
789
790 if(output_datum==ILLEGAL_TYPE) { /* if output_datum type is not specified by user*/
791 output_datum = DSET_BRICK_TYPE(*dset_time,0); /* get datum type from dataset */
792 }
793 nopt++;
794 continue;
795 }
796 /**** -datum type ****/
797
798 if( strcmp(argv[nopt],"-datum") == 0 ){
799 if( ++nopt >= argc )
800 NLfit_error("need an argument after -datum!\n") ;
801 if( strcmp(argv[nopt],"short") == 0 ){
802 output_datum = MRI_short ;
803 } else if( strcmp(argv[nopt],"float") == 0 ){
804 output_datum = MRI_float ;
805 } else {
806 ERROR_exit("-datum of type '%s' not supported in 3dNLfim!\n",argv[nopt]) ;
807 }
808 nopt++ ; continue ; /* go to next arg */
809 }
810
811 if( strncmp(argv[nopt],"-float",6) == 0 ){ /* 04 Dec 2008 */
812 output_datum = MRI_float ; nopt++ ; continue ;
813 }
814
815
816 /**** -mask mset [18 May 2000] ****/
817
818 if( strcmp(argv[nopt],"-mask") == 0 ){
819 THD_3dim_dataset * mset ; int ii,mc ;
820 nopt++ ;
821 if (nopt >= argc) NLfit_error ("need argument after -mask!") ;
822 mset = THD_open_dataset( argv[nopt] ) ;
823 if (mset == NULL) NLfit_error ("can't open -mask dataset!") ;
824
825 mask_vol = THD_makemask( mset , 0 , 1.0,0.0 ) ;
826 mask_nvox = DSET_NVOX(mset) ;
827 DSET_delete(mset) ;
828
829 if (mask_vol == NULL ) NLfit_error ("can't use -mask dataset!") ;
830 for( ii=mc=0 ; ii < mask_nvox ; ii++ ) if (mask_vol[ii]) mc++ ;
831 if (mc == 0) NLfit_error ("mask is all zeros!") ;
832 printf("++ %d voxels in mask %s\n",mc,argv[nopt]) ;
833 nopt++ ; continue ;
834 }
835
836
837 /*----- 22 July 1998: the -inTR option -----*/
838
839 if( strcmp(argv[nopt],"-inTR") == 0 ){
840 inTR = 1 ;
841 nopt++ ; continue ;
842 }
843
844 if( strcmp(argv[nopt],"-TR") == 0 ){
845 uuTR = (float)strtod(argv[++nopt],NULL) ;
846 if( uuTR <= 0.0f )
847 ERROR_message("value after -TR is not positive!") ;
848 nopt++ ; continue ;
849 }
850
851 /*----- -ignore num -----*/
852 if (strcmp(argv[nopt], "-ignore") == 0)
853 {
854 nopt++;
855 if (nopt >= argc) NLfit_error ("need argument after -ignore ");
856 sscanf (argv[nopt], "%d", &ival);
857 if (ival < 0)
858 NLfit_error ("illegal argument after -ignore ");
859 *ignore = ival;
860 nopt++;
861 continue;
862 }
863
864 /*----- -time filename -----*/
865 if (strcmp(argv[nopt], "-time") == 0)
866 {
867 nopt++;
868 if (nopt >= argc) NLfit_error ("need argument after -time ");
869 /* do not depend on input name lengths */
870 *tfilename = nifti_strdup(argv[nopt]);
871 MTEST (*tfilename);
872 nopt++;
873 continue;
874 }
875
876 /*----- -nabs -----*/
877 if (strcmp(argv[nopt], "-nabs") == 0)
878 {
879 *nabs = 1;
880 nopt++;
881 continue;
882 }
883
884 /*----- -POWELL -----*/
885
886 if( strcmp(argv[nopt],"-POWELL") == 0 ){ /* 20 Jul 2006 */
887 N_newuoa = 1 ; nopt++ ; continue ;
888 }
889 if( strcmp(argv[nopt],"-SIMPLEX") == 0 ){
890 N_newuoa = 0 ; nopt++ ; continue ;
891 }
892 if( strcmp(argv[nopt],"-BOTH") == 0 ){
893 N_newuoa = 2 ; nopt++ ; continue ;
894 }
895
896
897 /*----- -signal slabel -----*/
898 if (strcmp(argv[nopt], "-signal") == 0)
899 {
900 nopt++;
901 if (nopt >= argc) NLfit_error ("need argument after -signal ");
902 *sname = malloc (sizeof(char) * MAX_NAME_LENGTH);
903 MTEST (*sname);
904 strcpy (*sname, argv[nopt]);
905 initialize_signal_model (model_array, *sname,
906 smodel, p, *spname,
907 *min_sconstr, *max_sconstr);
908 nopt++;
909 continue;
910 }
911
912
913 /*----- -noise nlabel -----*/
914 if (strcmp(argv[nopt], "-noise") == 0)
915 {
916 nopt++;
917 if (nopt >= argc) NLfit_error ("need argument after -noise ");
918 *nname = malloc (sizeof(char) * MAX_NAME_LENGTH);
919 MTEST (*nname);
920 strcpy (*nname, argv[nopt]);
921 initialize_noise_model (model_array, *nname,
922 nmodel, r, *npname,
923 *min_nconstr, *max_nconstr);
924 nopt++;
925 continue;
926 }
927
928
929 /*----- -nrand n -----*/
930 if (strcmp(argv[nopt], "-nrand") == 0)
931 {
932 nopt++;
933 if (nopt >= argc) NLfit_error ("need argument after -nrand ");
934 sscanf (argv[nopt], "%d", &ival);
935 if (ival <= 0)
936 NLfit_error ("illegal argument after -nrand ");
937 *nrand = ival;
938 nopt++;
939 continue;
940 }
941
942
943 /*----- -nbest n -----*/
944 if (strcmp(argv[nopt], "-nbest") == 0)
945 {
946 nopt++;
947 if (nopt >= argc) NLfit_error ("need argument after -nbest ");
948 sscanf (argv[nopt], "%d", &ival);
949 if (ival <= 0)
950 NLfit_error ("illegal argument after -nbest ");
951 *nbest = ival;
952 nopt++;
953 continue;
954 }
955
956
957 /*----- -rmsmin r -----*/
958 if (strcmp(argv[nopt], "-rmsmin") == 0)
959 {
960 nopt++;
961 if (nopt >= argc) NLfit_error ("need argument after -rmsmin ");
962 sscanf (argv[nopt], "%f", &fval);
963 if (fval < 0.0)
964 NLfit_error ("illegal argument after -rmsmin ");
965 *rms_min = fval;
966 nopt++;
967 continue;
968 }
969
970
971 /*----- -fdisp fval -----*/
972 if (strcmp(argv[nopt], "-fdisp") == 0)
973 {
974 nopt++;
975 if (nopt >= argc) NLfit_error ("need argument after -fdisp ");
976 sscanf (argv[nopt], "%f", &fval);
977 *fdisp = fval;
978 nopt++;
979 continue;
980 }
981
982 /*----- -progress ival -----*/
983 if (strcmp(argv[nopt], "-progress") == 0)
984 {
985 nopt++;
986 if (nopt >= argc) NLfit_error ("need argument after -progress ");
987 sscanf (argv[nopt], "%d", &ival);
988 if (ival < 1)
989 NLfit_error("illegal argument after -progress ");
990 *progress = ival;
991 nopt++;
992 continue;
993 }
994
995 /*----- -voxel_count -----*/
996 if (strcmp(argv[nopt], "-voxel_count") == 0)
997 {
998 nopt++;
999 g_voxel_count = 1;
1000 continue;
1001 }
1002
1003 /*----- -jobs J -----*/
1004 if( strcmp(argv[nopt],"-jobs") == 0 ){ /* RWCox */
1005 nopt++ ;
1006 if (nopt >= argc) NLfit_error ("need J parameter after -jobs ");
1007 #ifdef PROC_MAX
1008 proc_numjob = strtol(argv[nopt],NULL,10) ;
1009 if( proc_numjob < 1 ){
1010 fprintf(stderr,"** setting number of processes to 1!\n") ;
1011 proc_numjob = 1 ;
1012 } else if( proc_numjob > PROC_MAX ){
1013 fprintf(stderr,"** setting number of processes to %d!\n",PROC_MAX);
1014 proc_numjob = PROC_MAX ;
1015 }
1016 #else
1017 fprintf(stderr,"** -jobs not supported in this version\n") ;
1018 #endif
1019 nopt++; continue;
1020 }
1021
1022
1023 /*----- check that user has specified the signal and noise models -----*/
1024 if ((*smodel) == NULL) ERROR_exit("Must specify signal model before '%s'",argv[nopt]);
1025 if ((*nmodel) == NULL) ERROR_exit("Must specify noise model before '%s'" ,argv[nopt]);
1026
1027 /*----- -sconstr k min max -----*/
1028 if (strcmp(argv[nopt], "-sconstr") == 0)
1029 {
1030 nopt++;
1031 if (nopt+2 >= argc) NLfit_error("need 3 arguments after -sconstr ");
1032
1033 sscanf (argv[nopt], "%d", &ival);
1034 if ((ival < 0) || (ival >= *p)) {
1035 fprintf(stderr,"*p = %d, ival = %d\n", *p, ival);
1036 ERROR_exit("Illegal 'k' index after -sconstr: legal range is 0..%d",*p-1) ;
1037 }
1038 index = ival;
1039 nopt++;
1040
1041 sscanf (argv[nopt], "%f", &fval);
1042 (*min_sconstr)[index] = fval;
1043 nopt++;
1044
1045 sscanf (argv[nopt], "%f", &fval);
1046 (*max_sconstr)[index] = fval;
1047 nopt++;
1048 continue;
1049 }
1050
1051
1052 /*----- -nconstr k min max -----*/
1053 if (strcmp(argv[nopt], "-nconstr") == 0)
1054 {
1055 nopt++;
1056 if (nopt+2 >= argc) NLfit_error("need 3 arguments after -nconstr ");
1057
1058 sscanf (argv[nopt], "%d", &ival);
1059 if ((ival < 0) || (ival >= *r))
1060 NLfit_error ("illegal argument after -nconstr ");
1061 index = ival;
1062 nopt++;
1063
1064 sscanf (argv[nopt], "%f", &fval);
1065 (*min_nconstr)[index] = fval;
1066 nopt++;
1067
1068 sscanf (argv[nopt], "%f", &fval);
1069 (*max_nconstr)[index] = fval;
1070 nopt++;
1071 continue;
1072 }
1073
1074
1075 /*----- -freg filename -----*/
1076 if (strcmp(argv[nopt], "-freg") == 0)
1077 {
1078 nopt++;
1079 if (nopt >= argc) NLfit_error ("need argument after -freg ");
1080 /* do not depend on input name lengths */
1081 *freg_filename = nifti_strdup(argv[nopt]);
1082 MTEST (*freg_filename);
1083 nopt++;
1084 continue;
1085 }
1086
1087
1088 /*----- -frsqr filename -----*/
1089 if (strcmp(argv[nopt], "-frsqr") == 0)
1090 {
1091 nopt++;
1092 if (nopt >= argc) NLfit_error ("need argument after -frsqr ");
1093 /* do not depend on input name lengths */
1094 *frsqr_filename = nifti_strdup(argv[nopt]);
1095 MTEST (*frsqr_filename);
1096 nopt++;
1097 continue;
1098 }
1099
1100
1101 /*----- -fsmax filename -----*/
1102 if (strcmp(argv[nopt], "-fsmax") == 0)
1103 {
1104 nopt++;
1105 if (nopt >= argc) NLfit_error ("need argument after -fsmax ");
1106 /* do not depend on input name lengths */
1107 *fsmax_filename = nifti_strdup(argv[nopt]);
1108 MTEST (*fsmax_filename);
1109 nopt++;
1110 continue;
1111 }
1112
1113
1114 /*----- -ftmax filename -----*/
1115 if (strcmp(argv[nopt], "-ftmax") == 0)
1116 {
1117 nopt++;
1118 if (nopt >= argc) NLfit_error ("need argument after -ftmax ");
1119 /* do not depend on input name lengths */
1120 *ftmax_filename = nifti_strdup(argv[nopt]);
1121 MTEST (*ftmax_filename);
1122 nopt++;
1123 continue;
1124 }
1125
1126
1127 /*----- -fpmax filename -----*/
1128 if (strcmp(argv[nopt], "-fpmax") == 0)
1129 {
1130 nopt++;
1131 if (nopt >= argc) NLfit_error ("need argument after -fpmax ");
1132 /* do not depend on input name lengths */
1133 *fpmax_filename = nifti_strdup(argv[nopt]);
1134 MTEST (*fpmax_filename);
1135 nopt++;
1136 continue;
1137 }
1138
1139
1140 /*----- -farea filename -----*/
1141 if (strcmp(argv[nopt], "-farea") == 0)
1142 {
1143 nopt++;
1144 if (nopt >= argc) NLfit_error ("need argument after -farea ");
1145 /* do not depend on input name lengths [2 Oct 2019 rickr] */
1146 *farea_filename = nifti_strdup(argv[nopt]);
1147 MTEST (*farea_filename);
1148 nopt++;
1149 continue;
1150 }
1151
1152
1153 /*----- -fparea filename -----*/
1154 if (strcmp(argv[nopt], "-fparea") == 0)
1155 {
1156 nopt++;
1157 if (nopt >= argc) NLfit_error ("need argument after -fparea ");
1158 /* do not depend on input name lengths [2 Oct 2019 rickr] */
1159 *fparea_filename = nifti_strdup(argv[nopt]);
1160 MTEST (*fparea_filename);
1161 nopt++;
1162 continue;
1163 }
1164
1165
1166 /*----- -fscoef k filename -----*/
1167 if (strcmp(argv[nopt], "-fscoef") == 0)
1168 {
1169 nopt++;
1170 if (nopt+1 >= argc) NLfit_error ("need 2 arguments after -fscoef ");
1171 sscanf (argv[nopt], "%d", &ival);
1172 if ((ival < 0) || (ival >= *p))
1173 NLfit_error ("illegal argument after -fscoef ");
1174 index = ival;
1175 nopt++;
1176
1177 /* do not depend on input name lengths [2 Oct 2019 rickr] */
1178 (*fscoef_filename)[index] = nifti_strdup(argv[nopt]);
1179 MTEST ((*fscoef_filename)[index]);
1180
1181 nopt++;
1182 continue;
1183 }
1184
1185
1186 /*----- -fncoef k filename -----*/
1187 if (strcmp(argv[nopt], "-fncoef") == 0)
1188 {
1189 nopt++;
1190 if (nopt+1 >= argc) NLfit_error ("need 2 arguments after -fncoef ");
1191 sscanf (argv[nopt], "%d", &ival);
1192 if ((ival < 0) || (ival >= *r))
1193 NLfit_error ("illegal argument after -fncoef ");
1194 index = ival;
1195 nopt++;
1196
1197 /* do not depend on input name lengths [2 Oct 2019 rickr] */
1198 (*fncoef_filename)[index] = nifti_strdup(argv[nopt]);
1199 MTEST ((*fncoef_filename)[index]);
1200
1201 nopt++;
1202 continue;
1203 }
1204
1205
1206 /*----- -tscoef k filename -----*/
1207 if (strcmp(argv[nopt], "-tscoef") == 0)
1208 {
1209 nopt++;
1210 if (nopt+1 >= argc) NLfit_error ("need 2 arguments after -tscoef ");
1211 sscanf (argv[nopt], "%d", &ival);
1212 if ((ival < 0) || (ival >= *p))
1213 NLfit_error ("illegal argument after -tscoef ");
1214 index = ival;
1215 nopt++;
1216
1217 /* do not depend on input name lengths [2 Oct 2019 rickr] */
1218 (*tscoef_filename)[index] = nifti_strdup(argv[nopt]);
1219 MTEST ((*tscoef_filename)[index]);
1220
1221 calc_tstats = 1;
1222
1223 nopt++;
1224 continue;
1225 }
1226
1227
1228 /*----- -tncoef k filename -----*/
1229 if (strcmp(argv[nopt], "-tncoef") == 0)
1230 {
1231 nopt++;
1232 if (nopt+1 >= argc) NLfit_error ("need 2 arguments after -tncoef ");
1233 sscanf (argv[nopt], "%d", &ival);
1234 if ((ival < 0) || (ival >= *r))
1235 NLfit_error ("illegal argument after -tncoef ");
1236 index = ival;
1237 nopt++;
1238
1239 /* do not depend on input name lengths [2 Oct 2019 rickr] */
1240 (*tncoef_filename)[index] = nifti_strdup(argv[nopt]);
1241 MTEST ((*tncoef_filename)[index]);
1242
1243 calc_tstats = 1;
1244
1245 nopt++;
1246 continue;
1247 }
1248
1249
1250 /*----- -bucket n prefixname -----*/
1251 if (strcmp(argv[nopt], "-bucket") == 0)
1252 {
1253 nopt++;
1254 if (nopt+1 >= argc) NLfit_error ("need 2 arguments after -bucket ");
1255 sscanf (argv[nopt], "%d", &ival);
1256 if ((ival < 0) || (ival > MAX_BRICKS))
1257 NLfit_error ("illegal argument after -bucket ");
1258 nopt++;
1259
1260 /* do not depend on input name lengths [2 Oct 2019 rickr] */
1261 option_data->bucket_filename = nifti_strdup(argv[nopt]);
1262 MTEST (option_data->bucket_filename);
1263
1264 /*----- set number of sub-bricks in the bucket -----*/
1265 if (ival == 0)
1266 nbricks = (*p) + (*r) + 8;
1267 else
1268 nbricks = ival;
1269 option_data->numbricks = nbricks;
1270
1271 /*----- allocate memory and initialize bucket dataset options -----*/
1272 option_data->brick_type = malloc (sizeof(int) * nbricks);
1273 MTEST (option_data->brick_type);
1274 option_data->brick_coef = malloc (sizeof(int) * nbricks);
1275 MTEST (option_data->brick_coef);
1276 option_data->brick_label = malloc (sizeof(char *) * nbricks);
1277 MTEST (option_data->brick_label);
1278 for (ibrick = 0; ibrick < nbricks; ibrick++)
1279 {
1280 option_data->brick_type[ibrick] = -1;
1281 option_data->brick_coef[ibrick] = -1;
1282 option_data->brick_label[ibrick] =
1283 malloc (sizeof(char) * MAX_NAME_LENGTH);
1284 MTEST (option_data->brick_label[ibrick]);
1285 }
1286
1287
1288 if (ival == 0)
1289 /*----- throw (almost) everything into the bucket -----*/
1290 {
1291 for (ibrick = 0; ibrick < (*r); ibrick++)
1292 {
1293 option_data->brick_type[ibrick] = FUNC_FIM_TYPE;
1294 option_data->brick_coef[ibrick] = ibrick;
1295 strcpy (option_data->brick_label[ibrick], (*npname)[ibrick]);
1296 }
1297
1298 for (ibrick = (*r); ibrick < (*p) + (*r); ibrick++)
1299 {
1300 option_data->brick_type[ibrick] = FUNC_FIM_TYPE;
1301 option_data->brick_coef[ibrick] = ibrick;
1302 strcpy (option_data->brick_label[ibrick],
1303 (*spname)[ibrick-(*r)]);
1304 }
1305
1306 ibrick = (*p) + (*r);
1307 option_data->brick_type[ibrick] = FUNC_FIM_TYPE;
1308 option_data->brick_coef[ibrick] = ibrick;
1309 strcpy (option_data->brick_label[ibrick], "Signal TMax");
1310
1311 ibrick++;
1312 option_data->brick_type[ibrick] = FUNC_FIM_TYPE;
1313 option_data->brick_coef[ibrick] = ibrick;
1314 strcpy (option_data->brick_label[ibrick], "Signal SMax");
1315
1316 ibrick++;
1317 option_data->brick_type[ibrick] = FUNC_FIM_TYPE;
1318 option_data->brick_coef[ibrick] = ibrick;
1319 strcpy (option_data->brick_label[ibrick], "Signal % SMax");
1320
1321 ibrick++;
1322 option_data->brick_type[ibrick] = FUNC_FIM_TYPE;
1323 option_data->brick_coef[ibrick] = ibrick;
1324 strcpy (option_data->brick_label[ibrick], "Signal Area");
1325
1326 ibrick++;
1327 option_data->brick_type[ibrick] = FUNC_FIM_TYPE;
1328 option_data->brick_coef[ibrick] = ibrick;
1329 strcpy (option_data->brick_label[ibrick], "Signal % Area");
1330
1331 ibrick++;
1332 option_data->brick_type[ibrick] = FUNC_FIM_TYPE;
1333 option_data->brick_coef[ibrick] = ibrick;
1334 strcpy (option_data->brick_label[ibrick], "Sigma Resid");
1335
1336 ibrick++;
1337 option_data->brick_type[ibrick] = FUNC_FIM_TYPE;
1338 option_data->brick_coef[ibrick] = ibrick;
1339 strcpy (option_data->brick_label[ibrick], "R^2");
1340
1341 ibrick++;
1342 option_data->brick_type[ibrick] = FUNC_FT_TYPE;
1343 option_data->brick_coef[ibrick] = ibrick;
1344 strcpy (option_data->brick_label[ibrick], "F-stat Regression");
1345
1346 }
1347
1348 nopt++;
1349 continue;
1350 }
1351
1352
1353 /*----- -brick m type k label -----*/
1354 if (strcmp(argv[nopt], "-brick") == 0)
1355 {
1356 nopt++;
1357 if (nopt+2 >= argc)
1358 NLfit_error ("need more arguments after -brick ");
1359 sscanf (argv[nopt], "%d", &ibrick);
1360 if ((ibrick < 0) || (ibrick >= option_data->numbricks))
1361 NLfit_error ("illegal argument after -brick ");
1362 nopt++;
1363
1364 if (strcmp(argv[nopt], "scoef") == 0)
1365 {
1366 option_data->brick_type[ibrick] = FUNC_FIM_TYPE;
1367
1368 nopt++;
1369 sscanf (argv[nopt], "%d", &ival);
1370 if ((ival < 0) || (ival > (*p)))
1371 NLfit_error ("illegal argument after scoef ");
1372 option_data->brick_coef[ibrick] = ival + (*r);
1373
1374 nopt++;
1375 if (nopt >= argc)
1376 NLfit_error ("need more arguments after -brick ");
1377 strcpy (option_data->brick_label[ibrick], argv[nopt]);
1378 }
1379
1380 else if (strcmp(argv[nopt], "ncoef") == 0)
1381 {
1382 option_data->brick_type[ibrick] = FUNC_FIM_TYPE;
1383
1384 nopt++;
1385 sscanf (argv[nopt], "%d", &ival);
1386 if ((ival < 0) || (ival > (*r)))
1387 NLfit_error ("illegal argument after ncoef ");
1388 option_data->brick_coef[ibrick] = ival;
1389
1390 nopt++;
1391 if (nopt >= argc)
1392 NLfit_error ("need more arguments after -brick ");
1393 strcpy (option_data->brick_label[ibrick], argv[nopt]);
1394 }
1395
1396 else if (strcmp(argv[nopt], "tmax") == 0)
1397 {
1398 option_data->brick_type[ibrick] = FUNC_FIM_TYPE;
1399 option_data->brick_coef[ibrick] = (*r) + (*p);
1400 nopt++;
1401 if (nopt >= argc)
1402 NLfit_error ("need more arguments after -brick ");
1403 strcpy (option_data->brick_label[ibrick], argv[nopt]);
1404 }
1405
1406 else if (strcmp(argv[nopt], "smax") == 0)
1407 {
1408 option_data->brick_type[ibrick] = FUNC_FIM_TYPE;
1409 option_data->brick_coef[ibrick] = (*r) + (*p) + 1;
1410 nopt++;
1411 if (nopt >= argc)
1412 NLfit_error ("need more arguments after -brick ");
1413 strcpy (option_data->brick_label[ibrick], argv[nopt]);
1414 }
1415
1416 else if (strcmp(argv[nopt], "psmax") == 0)
1417 {
1418 option_data->brick_type[ibrick] = FUNC_FIM_TYPE;
1419 option_data->brick_coef[ibrick] = (*r) + (*p) + 2;
1420 nopt++;
1421 if (nopt >= argc)
1422 NLfit_error ("need more arguments after -brick ");
1423 strcpy (option_data->brick_label[ibrick], argv[nopt]);
1424 }
1425
1426 else if (strcmp(argv[nopt], "area") == 0)
1427 {
1428 option_data->brick_type[ibrick] = FUNC_FIM_TYPE;
1429 option_data->brick_coef[ibrick] = (*r) + (*p) + 3;
1430 nopt++;
1431 if (nopt >= argc)
1432 NLfit_error ("need more arguments after -brick ");
1433 strcpy (option_data->brick_label[ibrick], argv[nopt]);
1434 }
1435
1436 else if (strcmp(argv[nopt], "parea") == 0)
1437 {
1438 option_data->brick_type[ibrick] = FUNC_FIM_TYPE;
1439 option_data->brick_coef[ibrick] = (*r) + (*p) + 4;
1440 nopt++;
1441 if (nopt >= argc)
1442 NLfit_error ("need more arguments after -brick ");
1443 strcpy (option_data->brick_label[ibrick], argv[nopt]);
1444 }
1445
1446 else if (strcmp(argv[nopt], "resid") == 0)
1447 {
1448 option_data->brick_type[ibrick] = FUNC_FIM_TYPE;
1449 option_data->brick_coef[ibrick] = (*r) + (*p) + 5;
1450 nopt++;
1451 if (nopt >= argc)
1452 NLfit_error ("need more arguments after -brick ");
1453 strcpy (option_data->brick_label[ibrick], argv[nopt]);
1454 }
1455
1456 else if (strcmp(argv[nopt], "rsqr") == 0)
1457 {
1458 option_data->brick_type[ibrick] = FUNC_FIM_TYPE;
1459 option_data->brick_coef[ibrick] = (*r) + (*p) + 6;
1460 nopt++;
1461 if (nopt >= argc)
1462 NLfit_error ("need more arguments after -brick ");
1463 strcpy (option_data->brick_label[ibrick], argv[nopt]);
1464 }
1465
1466 else if (strcmp(argv[nopt], "fstat") == 0)
1467 {
1468 option_data->brick_type[ibrick] = FUNC_FT_TYPE;
1469 option_data->brick_coef[ibrick] = (*r) + (*p) + 7;
1470 nopt++;
1471 if (nopt >= argc)
1472 NLfit_error ("need more arguments after -brick ");
1473 strcpy (option_data->brick_label[ibrick], argv[nopt]);
1474 }
1475
1476 else if (strcmp(argv[nopt], "tscoef") == 0)
1477 {
1478 option_data->brick_type[ibrick] = FUNC_TT_TYPE;
1479
1480 nopt++;
1481 sscanf (argv[nopt], "%d", &ival);
1482 if ((ival < 0) || (ival > (*p)))
1483 NLfit_error ("illegal argument after tscoef ");
1484 option_data->brick_coef[ibrick] = ival + (*r);
1485
1486 nopt++;
1487 if (nopt >= argc)
1488 NLfit_error ("need more arguments after -brick ");
1489 strcpy (option_data->brick_label[ibrick], argv[nopt]);
1490
1491 calc_tstats = 1;
1492 }
1493
1494 else if (strcmp(argv[nopt], "tncoef") == 0)
1495 {
1496 option_data->brick_type[ibrick] = FUNC_TT_TYPE;
1497
1498 nopt++;
1499 sscanf (argv[nopt], "%d", &ival);
1500 if ((ival < 0) || (ival > (*r)))
1501 NLfit_error ("illegal argument after tncoef ");
1502 option_data->brick_coef[ibrick] = ival;
1503
1504 nopt++;
1505 if (nopt >= argc)
1506 NLfit_error ("need more arguments after -brick ");
1507 strcpy (option_data->brick_label[ibrick], argv[nopt]);
1508
1509 calc_tstats = 1;
1510 }
1511
1512 else NLfit_error ("unable to interpret options after -brick ");
1513
1514 nopt++;
1515 continue;
1516 }
1517
1518
1519 /*----- -sfit filename -----*/
1520 if (strcmp(argv[nopt], "-sfit") == 0)
1521 {
1522 nopt++;
1523 if (nopt >= argc) NLfit_error ("need argument after -sfit ");
1524 /* do not depend on input name lengths [2 Oct 2019 rickr] */
1525 *sfit_filename = nifti_strdup(argv[nopt]);
1526 MTEST (*sfit_filename);
1527 nopt++;
1528 continue;
1529 }
1530
1531
1532 /*----- -snfit filename -----*/
1533 if (strcmp(argv[nopt], "-snfit") == 0)
1534 {
1535 nopt++;
1536 if (nopt >= argc) NLfit_error ("need argument after -snfit ");
1537 /* do not depend on input name lengths [2 Oct 2019 rickr] */
1538 *snfit_filename = nifti_strdup(argv[nopt]);
1539 MTEST (*snfit_filename);
1540 nopt++;
1541 continue;
1542 }
1543
1544
1545
1546 /*----- unknown command -----*/
1547 sprintf(message,"Unrecognized command line option: %s\n", argv[nopt]);
1548 NLfit_error (message);
1549
1550 }
1551
1552 /* maybe there was no input 13 Sep 2013 [rickr] */
1553 if ( *input_filename == NULL ) ERROR_exit("missing -input option");
1554
1555 if (*ignore < 0) {
1556 fprintf(stderr,"\n** NOTICE: "
1557 "ignore option now defaults to 0 instead of 3.\n\n");
1558 *ignore = 0;
1559 }
1560
1561 if( uuTR > 0.0f ) dsTR = uuTR ;
1562
1563 /*----- discard the model array -----*/
1564 DESTROY_MODEL_ARRAY (model_array) ;
1565
1566 }
1567
1568
1569 /*---------------------------------------------------------------------------*/
1570 /*
1571 Routine to check whether one output file already exists.
1572 */
1573
check_one_output_file(THD_3dim_dataset * dset_time,char * filename)1574 void check_one_output_file
1575 (
1576 THD_3dim_dataset * dset_time, /* input 3d+time data set */
1577 char * filename /* name of output file */
1578 )
1579
1580 {
1581 THD_3dim_dataset * new_dset=NULL; /* output afni data set pointer */
1582 int ierror; /* number of errors in editing data */
1583 char message[MAX_NAME_LENGTH]; /* error message */
1584
1585
1586 /*----- make an empty copy of input dataset -----*/
1587 new_dset = EDIT_empty_copy( dset_time ) ;
1588
1589
1590 ierror = EDIT_dset_items( new_dset ,
1591 ADN_prefix , filename ,
1592 ADN_label1 , filename ,
1593 ADN_self_name , filename ,
1594 ADN_type , ISHEAD(dset_time) ? HEAD_FUNC_TYPE :
1595 GEN_FUNC_TYPE ,
1596 ADN_none ) ;
1597
1598 if( ierror > 0 )
1599 {
1600 fprintf(stderr,
1601 "*** %d errors in attempting to create output dataset!\n",
1602 ierror ) ;
1603 exit(1) ;
1604 }
1605
1606 if( THD_is_file(new_dset->dblk->diskptr->header_name) )
1607 {
1608 fprintf(stderr, "** file already exists: %s\n",
1609 new_dset->dblk->diskptr->header_name);
1610 NLfit_error ("Output dataset file already exists");
1611 }
1612
1613 /*----- deallocate memory -----*/
1614 THD_delete_3dim_dataset( new_dset , False ) ; new_dset = NULL ;
1615
1616 }
1617
1618
1619 /*---------------------------------------------------------------------------*/
1620 /*
1621 Routine to check whether output files already exist.
1622 */
1623
check_output_files(char * freg_filename,char * frsqr_filename,char * fsmax_filename,char * ftmax_filename,char * fpmax_filename,char * farea_filename,char * fparea_filename,char ** fncoef_filename,char ** fscoef_filename,char ** tncoef_filename,char ** tscoef_filename,char * bucket_filename,char * sfit_filename,char * snfit_filename,THD_3dim_dataset * dset_time)1624 void check_output_files
1625 (
1626 char * freg_filename, /* file name for regression f-statistics */
1627 char * frsqr_filename, /* file name for R^2 statistics */
1628 char * fsmax_filename, /* file name for signal signed maximum */
1629 char * ftmax_filename, /* file name for epoch of signed maximum */
1630 char * fpmax_filename, /* file name for max. percentage change */
1631 char * farea_filename, /* file name for area under the signal */
1632 char * fparea_filename, /* file name for % area under the signal */
1633 char ** fncoef_filename, /* file name for noise model parameters */
1634 char ** fscoef_filename, /* file name for signal model parameters */
1635 char ** tncoef_filename, /* file name for noise model t-statistics */
1636 char ** tscoef_filename, /* file name for signal model t-statistics */
1637 char * bucket_filename, /* file name for bucket dataset */
1638 char * sfit_filename, /* file name for 3d+time fitted signal model */
1639 char * snfit_filename, /* file name for 3d+time fitted signal+noise */
1640 THD_3dim_dataset * dset_time /* input 3d+time data set */
1641 )
1642
1643 {
1644 int ip; /* parameter index */
1645
1646 if (freg_filename != NULL)
1647 check_one_output_file (dset_time, freg_filename);
1648
1649 if (frsqr_filename != NULL)
1650 check_one_output_file (dset_time, frsqr_filename);
1651
1652 if (fsmax_filename != NULL)
1653 check_one_output_file (dset_time, fsmax_filename);
1654
1655 if (ftmax_filename != NULL)
1656 check_one_output_file (dset_time, ftmax_filename);
1657
1658 if (fpmax_filename != NULL)
1659 check_one_output_file (dset_time, fpmax_filename);
1660
1661 if (farea_filename != NULL)
1662 check_one_output_file (dset_time, farea_filename);
1663
1664 if (fparea_filename != NULL)
1665 check_one_output_file (dset_time, fparea_filename);
1666
1667 for (ip = 0; ip < MAX_PARAMETERS; ip++)
1668 {
1669 if (fncoef_filename[ip] != NULL)
1670 check_one_output_file (dset_time, fncoef_filename[ip]);
1671 if (fscoef_filename[ip] != NULL)
1672 check_one_output_file (dset_time, fscoef_filename[ip]);
1673 if (tncoef_filename[ip] != NULL)
1674 check_one_output_file (dset_time, tncoef_filename[ip]);
1675 if (tscoef_filename[ip] != NULL)
1676 check_one_output_file (dset_time, tscoef_filename[ip]);
1677 }
1678
1679 if (bucket_filename != NULL)
1680 check_one_output_file (dset_time, bucket_filename);
1681
1682
1683 if (sfit_filename != NULL)
1684 check_one_output_file (dset_time, sfit_filename);
1685 if (snfit_filename != NULL)
1686 check_one_output_file (dset_time, snfit_filename);
1687
1688
1689 }
1690
1691
1692 /*---------------------------------------------------------------------------*/
1693 /*
1694 Routine to check for valid inputs.
1695 */
1696
check_for_valid_inputs(int nxyz,int r,int p,float * min_nconstr,float * max_nconstr,float * min_sconstr,float * max_sconstr,int nrand,int nbest,char * freg_filename,char * frsqr_filename,char * fsmax_filename,char * ftmax_filename,char * fpmax_filename,char * farea_filename,char * fparea_filename,char ** fncoef_filename,char ** fscoef_filename,char ** tncoef_filename,char ** tscoef_filename,char * bucket_filename,char * sfit_filename,char * snfit_filename,THD_3dim_dataset * dset_time)1697 void check_for_valid_inputs
1698 (
1699 int nxyz, /* number of voxels in image */
1700 int r, /* number of parameters in the noise model */
1701 int p, /* number of parameters in the signal model */
1702 float * min_nconstr, /* minimum parameter constraints for noise model */
1703 float * max_nconstr, /* maximum parameter constraints for noise model */
1704 float * min_sconstr, /* minimum parameter constraints for signal model */
1705 float * max_sconstr, /* maximum parameter constraints for signal model */
1706 int nrand, /* number of random vectors to generate */
1707 int nbest, /* number of random vectors to keep */
1708
1709 char * freg_filename, /* file name for regression f-statistics */
1710 char * frsqr_filename, /* file name for R^2 statistics */
1711 char * fsmax_filename, /* file name for signal signed maximum */
1712 char * ftmax_filename, /* file name for epoch of signed maximum */
1713 char * fpmax_filename, /* file name for max. percentage change */
1714 char * farea_filename, /* file name for area under the signal */
1715 char * fparea_filename, /* file name for % area under the signal */
1716 char ** fncoef_filename, /* file name for noise model parameters */
1717 char ** fscoef_filename, /* file name for signal model parameters */
1718 char ** tncoef_filename, /* file name for noise model t-statistics */
1719 char ** tscoef_filename, /* file name for signal model t-statistics */
1720 char * bucket_filename, /* file name for bucket dataset */
1721 char * sfit_filename, /* file name for 3d+time fitted signal model */
1722 char * snfit_filename, /* file name for 3d+time fitted signal+noise */
1723 THD_3dim_dataset * dset_time /* input 3d+time data set */
1724 )
1725
1726 {
1727 int ip; /* parameter index */
1728
1729
1730 /*----- check if mask is right size -----*/
1731 if (mask_vol != NULL)
1732 if (mask_nvox != nxyz)
1733 NLfit_error ("mask and input dataset bricks don't match in size");
1734
1735 /*----- 03 Feb 2009 - global mask for misfit reporting -----*/
1736 gmask = mask_vol ;
1737 if( gmask == NULL && nxyz > 666 ){
1738 MRI_IMAGE *qim ; int mc ;
1739 qim = THD_rms_brick( dset_time ) ;
1740 gmask = mri_automask_image( qim ) ;
1741 mri_free( qim ) ;
1742 mc = THD_countmask( DSET_NVOX(dset_time) , gmask ) ;
1743 if( mc <= 1 ){ free(gmask) ; gmask = NULL ; }
1744 }
1745 EDIT_set_misfit_mask(gmask) ;
1746
1747
1748 /*----- check for valid constraints -----*/
1749 for (ip = 0; ip < r; ip++)
1750 if (min_nconstr[ip] > max_nconstr[ip])
1751 NLfit_error ("Must have minimum constraints <= maximum constraints");
1752 for (ip = 0; ip < p; ip++)
1753 if (min_sconstr[ip] > max_sconstr[ip])
1754 NLfit_error("Must have mininum constraints <= maximum constraints");
1755
1756
1757 /*----- must have nbest <= nrand -----*/
1758 if (nrand < nbest)
1759 NLfit_error ("Must have nbest <= nrand");
1760
1761
1762 /*----- check whether any of the output files already exist -----*/
1763 if( THD_deathcon() )
1764 check_output_files (freg_filename, frsqr_filename,
1765 fsmax_filename, ftmax_filename,
1766 fpmax_filename, farea_filename, fparea_filename,
1767 fncoef_filename, fscoef_filename,
1768 tncoef_filename, tscoef_filename, bucket_filename,
1769 sfit_filename, snfit_filename, dset_time);
1770
1771 }
1772
1773 /*---------------------------------------------------------------------------*/
1774 /*
1775 Allocate volume memory and fill with zeros.
1776 */
1777
zero_fill_volume(float ** fvol,int nxyz)1778 void zero_fill_volume (float ** fvol, int nxyz)
1779 {
1780 int ixyz;
1781
1782 if( proc_numjob == 1 ){ /* 1 process ==> allocate locally */
1783
1784 *fvol = (float *) malloc (sizeof(float) * nxyz); MTEST(*fvol);
1785 for (ixyz = 0; ixyz < nxyz; ixyz++)
1786 (*fvol)[ixyz] = 0.0;
1787
1788 }
1789 #ifdef PROC_MAX
1790 else { /* multiple processes ==> prepare for shared memory */
1791 /* by remembering what to do */
1792
1793 proc_shm_arnum++ ;
1794 proc_shm_arsiz = (int *) realloc( proc_shm_arsiz ,
1795 sizeof(int) *proc_shm_arnum ) ;
1796 proc_shm_ar = (float ***) realloc( proc_shm_ar ,
1797 sizeof(float **)*proc_shm_arnum ) ;
1798 proc_shm_arsiz[proc_shm_arnum-1] = nxyz ;
1799 proc_shm_ar[proc_shm_arnum-1] = fvol ;
1800
1801 /* actual allocation and pointer assignment (to *fvol)
1802 will take place in function proc_finalize_shm_volumes() */
1803 }
1804 #endif
1805 }
1806
1807
1808 /*---------------------------------------------------------------------------*/
1809
1810 #ifdef PROC_MAX
1811
1812 /*** signal handler for fatal errors -- make sure shared memory is deleted ***/
1813
1814 #include <signal.h>
1815
proc_sigfunc(int sig)1816 void proc_sigfunc(int sig)
1817 {
1818 char * sname ; int ii ;
1819 static volatile int fff=0 ;
1820 if( fff ) _exit(1); else fff=1 ;
1821 switch(sig){
1822 default: sname = "unknown" ; break ;
1823 case SIGHUP: sname = "SIGHUP" ; break ;
1824 case SIGTERM: sname = "SIGTERM" ; break ;
1825 case SIGILL: sname = "SIGILL" ; break ;
1826 case SIGKILL: sname = "SIGKILL" ; break ;
1827 case SIGPIPE: sname = "SIGPIPE" ; break ;
1828 case SIGSEGV: sname = "SIGSEGV" ; break ;
1829 case SIGBUS: sname = "SIGBUS" ; break ;
1830 case SIGINT: sname = "SIGINT" ; break ;
1831 case SIGFPE: sname = "SIGFPE" ; break ;
1832 }
1833 if( proc_shmid > 0 ){
1834 shmctl( proc_shmid , IPC_RMID , NULL ) ; proc_shmid = 0 ;
1835 }
1836 fprintf(stderr,"Fatal Signal %d (%s) received in job #%d\n",
1837 sig,sname,proc_ind) ;
1838 exit(1) ;
1839 }
1840
proc_atexit(void)1841 void proc_atexit( void ) /*** similarly - atexit handler ***/
1842 {
1843 if( proc_shmid > 0 )
1844 shmctl( proc_shmid , IPC_RMID , NULL ) ;
1845 }
1846
1847 /*---------------------------------------------------------------*/
1848 /*** This function is called to allocate all output
1849 volumes at once in shared memory, and set their pointers ***/
1850 /** 24 Oct 2006: use mmap() instead of shared memory, if possible **/
1851
1852 #include <sys/mman.h>
1853 #if !defined(MAP_ANON) && defined(MAP_ANONYMOUS)
1854 # define MAP_ANON MAP_ANONYMOUS
1855 #endif
1856
1857
proc_finalize_shm_volumes(void)1858 void proc_finalize_shm_volumes(void)
1859 {
1860 char kstr[32] ; int ii ;
1861 long long psum , twogig = 2ll * 1024ll * 1024ll * 1024ll ;
1862
1863 if( proc_shm_arnum == 0 ) return ; /* should never happen */
1864
1865 /** 21 Oct 2005: check in big arithmetic how much mem we need **/
1866
1867 psum = 0 ;
1868 for( ii=0 ; ii < proc_shm_arnum ; ii++ )
1869 psum += (long long)proc_shm_arsiz[ii] ;
1870 psum *= sizeof(float) ;
1871 #ifdef MAP_ANON
1872 if( psum >= twogig &&
1873 ( sizeof(void *) < 8 || sizeof(size_t) < 8 ) ) /* too much for 32-bit */
1874 #else
1875 if( psum >= twogig ) /* too much for shmem */
1876 #endif
1877 ERROR_exit(
1878 "Total shared memory needed = %lld >= %lld (2 GB)\n"
1879 "** SUGGESTION 1: Use 3dAutobox to automatically eliminate non-brain\n"
1880 " areas from the 3d+time input dataset and reduce memory \n"
1881 " requirements, e.g.\n"
1882 " 3dAutobox -prefix Subj1AllRuns_Smaller -input Subj1AllRuns\n"
1883 " Then run 3dNLfim again with the smaller 3d+time input dataset\n"
1884 "\n"
1885 "** SUGGESTION 2: Use 3dZcutup to slice dataset into smaller pieces\n"
1886 "** and then 3dZcat to glue results back together.\n"
1887 "\n"
1888 "** SUGGESTION 3: Run on a 64-bit computer system, instead of 32-bit.\n"
1889 , psum,twogig) ;
1890 else
1891 INFO_message("total shared memory needed = %s bytes (about %s)" ,
1892 commaized_integer_string(psum) ,
1893 approximate_number_string((double)psum) ) ;
1894
1895 proc_shmsize = psum ; /* global variable */
1896 #if 0
1897 /* old code */
1898 if( proc_shm_arnum == 0 ) return ; /* should never happen */
1899
1900 proc_shmsize = 0 ; /* add up sizes of */
1901 for( ii=0 ; ii < proc_shm_arnum ; ii++ ) /* all arrays for */
1902 proc_shmsize += proc_shm_arsiz[ii] ; /* shared memory */
1903
1904 proc_shmsize *= sizeof(float) ; /* convert to byte count */
1905 #endif
1906
1907 /* create shared memory segment */
1908 #ifdef MAP_ANON /** 24 Oct 2005: use mmap() instead of shmem **/
1909
1910 proc_shmptr = mmap( (void *)0 , (size_t)psum ,
1911 PROT_READ | PROT_WRITE , MAP_ANON | MAP_SHARED ,
1912 -1 , 0 ) ;
1913 if( proc_shmptr == NULL || proc_shmptr == (void *)(-1) ){
1914 perror("** FATAL ERROR: Can't create shared mmap() segment\n"
1915 "** Unix message") ;
1916 INFO_message("\n\n"
1917 "** SUGGESTION: Use 3dZcutup to slice dataset into smaller pieces\n"
1918 "** and then 3dZcat to glue results back together.\n"
1919 "** SUGGESTION: Run on a 64-bit computer system, instead of 32-bit.\n ");
1920 exit(1) ;
1921 }
1922
1923 #else /** old code: shared memory segment **/
1924
1925 UNIQ_idcode_fill( kstr ) ; /* unique string "key" */
1926 proc_shmid = shm_create( kstr , proc_shmsize ) ; /* thd_iochan.c */
1927 if( proc_shmid < 0 ){
1928 fprintf(stderr,"\n** Can't create shared memory of size %d!\n"
1929 "** Try re-running without -jobs option!\n" ,
1930 proc_shmsize ) ;
1931
1932 /** if failed, print out some advice on how to tune SHMMAX **/
1933
1934 { char *cmd=NULL ;
1935 #if defined(LINUX)
1936 cmd = "cat /proc/sys/kernel/shmmax" ;
1937 #elif defined(SOLARIS)
1938 cmd = "/usr/sbin/sysdef | grep SHMMAX" ;
1939 #elif defined(DARWIN)
1940 cmd = "sysctl -n kern.sysv.shmmax" ;
1941 #endif
1942 if( cmd != NULL ){
1943 FILE *fp = popen( cmd , "r" ) ;
1944 if( fp != NULL ){
1945 long long smax=0 ;
1946 fscanf(fp,"%lld",&smax) ; pclose(fp) ;
1947 if( smax > 0 )
1948 fprintf(stderr ,
1949 "\n"
1950 "** POSSIBLY USEFUL ADVICE:\n"
1951 "** Current max shared memory size = %u bytes.\n"
1952 "** For information on how to change this, see\n"
1953 "** https://afni.nimh.nih.gov/afni/parallize.htm\n"
1954 "** and also contact your system administrator.\n"
1955 , smax ) ;
1956 }
1957 }
1958 }
1959 exit(1) ;
1960 }
1961 #endif /* MAP_ANON */
1962
1963 /* set a signal handler to catch most fatal errors and
1964 delete the shared memory segment if program crashes */
1965
1966 signal(SIGPIPE,proc_sigfunc) ; signal(SIGSEGV,proc_sigfunc) ;
1967 signal(SIGINT ,proc_sigfunc) ; signal(SIGFPE ,proc_sigfunc) ;
1968 signal(SIGBUS ,proc_sigfunc) ; signal(SIGHUP ,proc_sigfunc) ;
1969 signal(SIGTERM,proc_sigfunc) ; signal(SIGILL ,proc_sigfunc) ;
1970 signal(SIGKILL,proc_sigfunc) ; signal(SIGPIPE,proc_sigfunc) ;
1971 atexit(proc_atexit) ; /* or when the program exits */
1972
1973 #ifdef MAP_ANON
1974 fprintf(stderr , "++ mmap() memory allocated: %lld bytes\n" ,
1975 proc_shmsize ) ;
1976 #else
1977 fprintf(stderr , "++ Shared memory allocated: %lld bytes at id=%d\n" ,
1978 proc_shmsize , proc_shmid ) ;
1979 #endif
1980
1981 /* get pointer to shared memory segment we just created */
1982 #ifndef MAP_ANON
1983 if(proc_shmid > 0) {
1984 proc_shmptr = shm_attach( proc_shmid ) ; /* thd_iochan.c */
1985 if( proc_shmptr == NULL ){
1986 fprintf(stderr,"\n** Can't attach to shared memory!?\n"
1987 "** This is bizarre.\n" ) ;
1988 shmctl( proc_shmid , IPC_RMID , NULL ) ;
1989 exit(1) ;
1990 }
1991 }
1992 #endif
1993
1994 /* clear all shared memory to zero*/
1995
1996 memset( proc_shmptr , 0 , proc_shmsize ) ;
1997
1998 /* fix the local pointers to arrays in shared memory */
1999
2000 *proc_shm_ar[0] = (float *) proc_shmptr ;
2001 for( ii=1 ; ii < proc_shm_arnum ; ii++ )
2002 *proc_shm_ar[ii] = *proc_shm_ar[ii-1] + proc_shm_arsiz[ii] ;
2003
2004 return;
2005 }
2006
2007 /*-------------------------------------------------------------*/
2008 /*** This function replaces free();
2009 it won't try to free things stored in the shared memory ***/
2010
proc_free(void * ptr)2011 void proc_free( void *ptr )
2012 {
2013 int ii ;
2014
2015 if( ptr == NULL ) return ;
2016 /* from if ( proc_shmid == 0 ) 25 Oct 2006 DRG */
2017 if( proc_shmptr == NULL ){ free(ptr); return; } /* no shm */
2018 for( ii=0 ; ii < proc_shm_arnum ; ii++ )
2019 if( ((float *)ptr) == *proc_shm_ar[ii] ) return;
2020 free(ptr); return;
2021 }
2022
2023 #undef free /* replace use of library free() */
2024 #define free proc_free /* with proc_free() function */
2025
2026 #endif /* PROC_MAX */
2027
2028 /*---------------------------------------------------------------------------*/
2029 /*
2030 Perform all program initialization.
2031 */
2032
initialize_program(int argc,char ** argv,int * ignore,char ** nname,char ** sname,vfp * nmodel,vfp * smodel,int * r,int * p,char *** npname,char *** spname,float ** min_nconstr,float ** max_nconstr,float ** min_sconstr,float ** max_sconstr,int * nabs,int * nrand,int * nbest,float * rms_min,float * fdisp,int * progress,char ** input_filename,char ** tfilename,char ** freg_filename,char ** frsqr_filename,char ** fsmax_filename,char ** ftmax_filename,char ** fpmax_filename,char ** farea_filename,char ** fparea_filename,char *** fncoef_filename,char *** fscoef_filename,char *** tncoef_filename,char *** tscoef_filename,char ** sfit_filename,char ** snfit_filename,THD_3dim_dataset ** dset_time,int * nxyz,int * ts_length,float *** x_array,float ** ts_array,float ** par_rdcd,float ** par_full,float ** tpar_full,float ** rmsreg_vol,float ** freg_vol,float ** rsqr_vol,float ** smax_vol,float ** tmax_vol,float ** pmax_vol,float ** area_vol,float ** parea_vol,float *** ncoef_vol,float *** scoef_vol,float *** tncoef_vol,float *** tscoef_vol,float *** sfit_vol,float *** snfit_vol,NL_options * option_data)2033 void initialize_program
2034 (
2035 int argc, /* number of input arguments */
2036 char ** argv, /* array of input arguments */
2037 int * ignore, /* delete this number of points from time series */
2038 char ** nname, /* name of noise model */
2039 char ** sname, /* name of signal model */
2040 vfp * nmodel, /* pointer to noise model */
2041 vfp * smodel, /* pointer to signal model */
2042 int *r, /* number of parameters in the noise model */
2043 int *p, /* number of parameters in the signal model */
2044 char *** npname, /* noise parameter names */
2045 char *** spname, /* signal parameter names */
2046 float ** min_nconstr, /* minimum parameter constraints for noise model */
2047 float ** max_nconstr, /* maximum parameter constraints for noise model */
2048 float ** min_sconstr, /* minimum parameter constraints for signal model */
2049 float ** max_sconstr, /* maximum parameter constraints for signal model */
2050 int * nabs, /* use absolute constraints for noise parameters */
2051 int * nrand, /* number of random vectors to generate */
2052 int * nbest, /* number of random vectors to keep */
2053 float * rms_min, /* minimum rms error to reject reduced model */
2054 float * fdisp, /* minimum f-statistic for display */
2055 int *progress, /* progress report every nth voxel */
2056 char ** input_filename, /* file name of input 3d+time dataset */
2057 char ** tfilename, /* file name for time point series */
2058 char ** freg_filename, /* file name for regression f-statistics */
2059 char ** frsqr_filename, /* file name for R^2 statistics */
2060 char ** fsmax_filename, /* file name for signal signed maximum */
2061 char ** ftmax_filename, /* file name for time of signed maximum */
2062 char ** fpmax_filename, /* file name for max. percentage change */
2063 char ** farea_filename, /* file name for area under the signal */
2064 char ** fparea_filename, /* file name for % area under the signal */
2065 char *** fncoef_filename, /* file name for noise model parameters */
2066 char *** fscoef_filename, /* file name for signal model parameters */
2067 char *** tncoef_filename, /* file name for noise model t-statistics */
2068 char *** tscoef_filename, /* file name for signal model t-statistics */
2069 char ** sfit_filename, /* file name for 3d+time fitted signal model */
2070 char ** snfit_filename, /* file name for 3d+time fitted signal+noise */
2071
2072 THD_3dim_dataset ** dset_time, /* input 3d+time data set */
2073 int * nxyz, /* number of voxels in image */
2074 int * ts_length, /* length of time series data */
2075
2076 float *** x_array, /* independent variable matrix */
2077 float ** ts_array, /* input time series array */
2078
2079 float ** par_rdcd, /* estimated parameters for the reduced model */
2080 float ** par_full, /* estimated parameters for the full model */
2081 float ** tpar_full, /* t-statistic of the parameters in full model */
2082
2083 float ** rmsreg_vol, /* root-mean-square error for the full model */
2084 float ** freg_vol, /* f-statistic for the full regression model */
2085 float ** rsqr_vol, /* R^2 volume data */
2086 float ** smax_vol, /* volume of signed maximum of signal */
2087 float ** tmax_vol, /* volume of epoch of signed maximum of signal */
2088 float ** pmax_vol, /* max. percentage change due to signal */
2089 float ** area_vol, /* area between signal and baseline */
2090 float ** parea_vol, /* percentage area between signal and baseline */
2091 float *** ncoef_vol, /* volume of noise model parameters */
2092 float *** scoef_vol, /* volume of signal model parameters */
2093 float *** tncoef_vol, /* volume of noise model t-statistics */
2094 float *** tscoef_vol, /* volume of signal model t-statistics */
2095 float *** sfit_vol, /* voxelwise 3d+time fitted signal model */
2096 float *** snfit_vol, /* voxelwise 3d+time fitted signal+noise model */
2097
2098 NL_options * option_data /* bucket dataset options */
2099
2100 )
2101
2102 {
2103 int dimension; /* dimension of full model */
2104 int ip; /* parameter index */
2105 int it; /* time index */
2106 MRI_IMAGE * im, * flim; /* pointers to image structures
2107 -- used to read 1D ASCII */
2108 int nt; /* number of points in 1D x data file */
2109 float * tar;
2110 int ixyz; /* voxel index */
2111
2112
2113 /*----- save command line for history notes -----*/
2114 commandline = tross_commandline( PROGRAM_NAME , argc,argv ) ;
2115
2116
2117 /*----- get command line inputs -----*/
2118 get_options(argc, argv, ignore, nname, sname, nmodel, smodel,
2119 r, p, npname, spname,
2120 min_nconstr, max_nconstr, min_sconstr, max_sconstr, nabs,
2121 nrand, nbest, rms_min, fdisp, progress, input_filename, tfilename,
2122 freg_filename, frsqr_filename, fsmax_filename,
2123 ftmax_filename, fpmax_filename, farea_filename, fparea_filename,
2124 fncoef_filename, fscoef_filename,
2125 tncoef_filename, tscoef_filename, sfit_filename, snfit_filename,
2126 dset_time, nxyz, ts_length, option_data);
2127
2128
2129 /*----- check for valid inputs -----*/
2130 check_for_valid_inputs (*nxyz, *r, *p, *min_nconstr, *max_nconstr,
2131 *min_sconstr, *max_sconstr, *nrand, *nbest,
2132 *freg_filename, *frsqr_filename,
2133 *fsmax_filename, *ftmax_filename,
2134 *fpmax_filename, *farea_filename, *fparea_filename,
2135 *fncoef_filename, *fscoef_filename,
2136 *tncoef_filename, *tscoef_filename,
2137 *sfit_filename, *snfit_filename,
2138 option_data->bucket_filename, *dset_time);
2139
2140
2141 /*----- allocate space for input time series -----*/
2142 *ts_length = *ts_length - *ignore;
2143 *ts_array = (float *) malloc (sizeof(float) * (*ts_length));
2144 MTEST (*ts_array);
2145
2146
2147 /*----- allocate space for independent variable matrix -----*/
2148 *x_array = (float **) malloc (sizeof(float *) * (*ts_length));
2149 MTEST (*x_array);
2150 for (it = 0; it < (*ts_length); it++)
2151 {
2152 (*x_array)[it] = (float *) malloc (sizeof(float) * 3);
2153 MTEST ((*x_array)[it]);
2154 }
2155
2156
2157 /*----- initialize independent variable matrix -----*/
2158 if (*tfilename == NULL)
2159 {
2160 if( (inTR || uuTR > 0.0f) && dsTR > 0.0 ){ /* 22 July 1998 */
2161 DELT = dsTR ;
2162 fprintf(stderr,"--- computing with TR = %g\n",DELT) ;
2163 }
2164 for (it = 0; it < (*ts_length); it++)
2165 {
2166 (*x_array)[it][0] = 1.0;
2167 (*x_array)[it][1] = it * DELT;
2168 (*x_array)[it][2] = (it * DELT) * (it * DELT);
2169 }
2170 }
2171 else
2172 {
2173 flim = mri_read_1D(*tfilename) ;
2174 if (flim == NULL)
2175 NLfit_error ("Unable to read time reference file \n");
2176 nt = flim -> nx;
2177 if (nt < (*ts_length))
2178 NLfit_error ("Time reference array is too short");
2179 tar = MRI_FLOAT_PTR(flim) ;
2180 for (it = 0; it < (*ts_length); it++)
2181 {
2182 (*x_array)[it][0] = 1.0;
2183 (*x_array)[it][1] = tar[it] ;
2184 (*x_array)[it][2] = tar[it] * tar[it];
2185 }
2186 mri_free (flim);
2187 }
2188
2189 /*--- 24 Jul 2006: special change to x_array[][2] for Linear+Ort [RWCox] ---*/
2190
2191 if( strcmp(*nname,"Linear+Ort") == 0 ){
2192 char *fname ; MRI_IMAGE *fim ; int nx ; float *far ;
2193 fname = my_getenv("AFNI_ORTMODEL_REF") ;
2194 if( fname == NULL )
2195 ERROR_exit("Linear+Ort model: 'AFNI_ORTMODEL_REF' not set") ;
2196
2197 fim = mri_read_1D(fname) ;
2198 if( fim == NULL || fim->nx < 2 )
2199 ERROR_exit(
2200 "Linear+Ort model: can't read file AFNI_ORTMODEL_REF='%s'",fname) ;
2201
2202 nx = fim->nx ; far = MRI_FLOAT_PTR(fim) ;
2203 if( nx != (*ts_length) || fim->ny > 1 )
2204 WARNING_message(
2205 "Linear+Ort model: AFNI_ORTMODEL_REF='%s' has nx=%d ny=%d; data length=%d",
2206 fname , nx , fim->ny , *ts_length ) ;
2207 else
2208 INFO_message(
2209 "Linear+Ort model: AFNI_ORTMODEL_REF='%s' loaded; length=%d" ,
2210 fname , nx ) ;
2211
2212 for( it=0 ; it < (*ts_length); it++)
2213 (*x_array)[it][2] = (it < nx) ? far[it] : 0.0f ;
2214 }
2215
2216 /*----- allocate memory space for parameters -----*/
2217 dimension = (*r) + (*p);
2218 *par_rdcd = (float *) malloc (sizeof(float) * dimension);
2219 MTEST (*par_rdcd);
2220 *par_full = (float *) malloc (sizeof(float) * dimension);
2221 MTEST (*par_full);
2222 *tpar_full = (float *) malloc (sizeof(float) * dimension);
2223 MTEST (*tpar_full);
2224
2225
2226 /*----- allocate memory space for volume data -----*/
2227 zero_fill_volume( freg_vol , *nxyz ) ;
2228
2229 if ((*freg_filename != NULL) || (option_data->bucket_filename != NULL))
2230 zero_fill_volume( rmsreg_vol , *nxyz ) ;
2231
2232 if ((*frsqr_filename != NULL) || (option_data->bucket_filename != NULL))
2233 zero_fill_volume( rsqr_vol , *nxyz ) ;
2234
2235 if ((*fsmax_filename != NULL) || (option_data->bucket_filename != NULL))
2236 zero_fill_volume( smax_vol , *nxyz ) ;
2237
2238 if ((*ftmax_filename != NULL) || (option_data->bucket_filename != NULL))
2239 zero_fill_volume( tmax_vol , *nxyz ) ;
2240
2241
2242 if ((*fpmax_filename != NULL) || (option_data->bucket_filename != NULL))
2243 zero_fill_volume( pmax_vol , *nxyz ) ;
2244
2245 if ((*farea_filename != NULL) || (option_data->bucket_filename != NULL))
2246 zero_fill_volume( area_vol , *nxyz ) ;
2247
2248 if ((*fparea_filename != NULL) || (option_data->bucket_filename != NULL))
2249 zero_fill_volume( parea_vol , *nxyz ) ;
2250
2251
2252 *ncoef_vol = (float **) malloc (sizeof(float *) * (*r));
2253 MTEST (*ncoef_vol);
2254 *tncoef_vol = (float **) malloc (sizeof(float *) * (*r));
2255 MTEST (*tncoef_vol);
2256
2257 for (ip = 0; ip < (*r); ip++)
2258 {
2259 if (((*fncoef_filename)[ip] != NULL) || ((*tncoef_filename)[ip] != NULL)
2260 || (option_data->bucket_filename != NULL))
2261 zero_fill_volume( &((*ncoef_vol)[ip]) , *nxyz ) ;
2262 else
2263 (*ncoef_vol)[ip] = NULL;
2264
2265 if (((*tncoef_filename)[ip] != NULL)
2266 || (option_data->bucket_filename != NULL))
2267 zero_fill_volume( &((*tncoef_vol)[ip]) , *nxyz ) ;
2268 else
2269 (*tncoef_vol)[ip] = NULL;
2270 }
2271
2272
2273 *scoef_vol = (float **) malloc (sizeof(float *) * (*p));
2274 MTEST (*scoef_vol);
2275 *tscoef_vol = (float **) malloc (sizeof(float *) * (*p));
2276 MTEST (*tscoef_vol);
2277
2278 for (ip = 0; ip < (*p); ip++)
2279 {
2280 if (((*fscoef_filename)[ip] != NULL) || ((*tscoef_filename)[ip] != NULL)
2281 || (option_data->bucket_filename != NULL))
2282 zero_fill_volume( &((*scoef_vol)[ip]) , *nxyz ) ;
2283 else
2284 (*scoef_vol)[ip] = NULL;
2285
2286 if (((*tscoef_filename)[ip] != NULL)
2287 || (option_data->bucket_filename != NULL))
2288 zero_fill_volume( &((*tscoef_vol)[ip]) , *nxyz ) ;
2289 else
2290 (*tscoef_vol)[ip] = NULL;
2291 }
2292
2293
2294 /*----- Allocate memory space for 3d+time fitted signal model -----*/
2295 if (*sfit_filename != NULL)
2296 {
2297 *sfit_vol = (float **) malloc (sizeof(float *) * (*ts_length));
2298 MTEST(*sfit_vol);
2299
2300 for (it = 0; it < *ts_length; it++)
2301 zero_fill_volume( &((*sfit_vol)[it]) , *nxyz ) ;
2302 }
2303
2304 /*----- Allocate memory space for 3d+time fitted signal+noise model -----*/
2305 if (*snfit_filename != NULL)
2306 {
2307 *snfit_vol = (float **) malloc (sizeof(float *) * (*ts_length));
2308 MTEST(*snfit_vol);
2309
2310 for (it = 0; it < *ts_length; it++)
2311 zero_fill_volume( &((*snfit_vol)[it]) , *nxyz ) ;
2312 }
2313
2314 #ifdef PROC_MAX
2315 if( proc_numjob > 1 ) proc_finalize_shm_volumes() ; /* RWCox */
2316 #endif
2317
2318 }
2319
2320
2321 /*---------------------------------------------------------------------------*/
2322 /*
2323 Get the time series for one voxel from the AFNI 3d+time data set.
2324 */
2325
read_ts_array(THD_3dim_dataset * dset_time,int iv,int ts_length,int ignore,float * ts_array)2326 void read_ts_array
2327 (
2328 THD_3dim_dataset * dset_time, /* input 3d+time data set */
2329 int iv, /* get time series for this voxel */
2330 int ts_length, /* length of time series array */
2331 int ignore, /* delete this number of points from time series */
2332 float * ts_array /* input time series data for voxel #iv */
2333 )
2334
2335 {
2336 MRI_IMAGE * im; /* intermediate float data */
2337 float * ar; /* pointer to float data */
2338 int it; /* time index */
2339
2340 ENTRY("read_ts_array") ;
2341
2342
2343 /*----- Extract time series from 3d+time data set into MRI_IMAGE -----*/
2344 STATUS("extracting time series") ;
2345 im = THD_extract_series (iv, dset_time, 0);
2346
2347
2348 /*----- Verify extraction -----*/
2349 if (im == NULL) NLfit_error ("Unable to extract data from 3d+time dataset");
2350
2351
2352 /*----- Now extract time series from MRI_IMAGE -----*/
2353 STATUS("loading into ts_array") ;
2354 ar = MRI_FLOAT_PTR (im);
2355 for (it = 0; it < ts_length; it++)
2356 {
2357 ts_array[it] = ar[it+ignore];
2358 }
2359
2360
2361 /*----- Release memory -----*/
2362 mri_free (im); im = NULL;
2363 EXRETURN ;
2364 }
2365
2366
2367 /*---------------------------------------------------------------------------*/
2368 /*
2369 Print out the given time series.
2370 */
2371
write_ts_array(int ts_length,float * ts_array)2372 void write_ts_array
2373 (
2374 int ts_length, /* length of time series array */
2375 float * ts_array /* time series data to be printed */
2376 )
2377
2378 {
2379 int it; /* time index */
2380
2381 for (it = 0; it < ts_length; it++)
2382 printf ("%d %f \n", it, ts_array[it]);
2383 }
2384
2385
2386 /*---------------------------------------------------------------------------*/
2387 /*
2388 Save results for output later.
2389 */
2390
save_results(int iv,vfp nmodel,vfp smodel,int r,int p,int novar,int ts_length,float ** x_array,float * par_full,float * tpar_full,float rmsreg,float freg,float rsqr,float smax,float tmax,float pmax,float area,float parea,float * rmsreg_vol,float * freg_vol,float * rsqr_vol,float * smax_vol,float * tmax_vol,float * pmax_vol,float * area_vol,float * parea_vol,float ** ncoef_vol,float ** scoef_vol,float ** tncoef_vol,float ** tscoef_vol,float ** sfit_vol,float ** snfit_vol)2391 void save_results
2392 (
2393 int iv, /* current voxel number */
2394 vfp nmodel, /* pointer to noise model */
2395 vfp smodel, /* pointer to signal model */
2396 int r, /* number of parameters in the noise model */
2397 int p, /* number of parameters in the signal model */
2398 int novar, /* flag for insufficient variation in the data */
2399 int ts_length, /* length of time series data */
2400 float ** x_array, /* independent variable matrix */
2401 float * par_full, /* estimated parameters for the full model */
2402 float * tpar_full, /* t-statistic of the parameters in full model */
2403 float rmsreg, /* root-mean-square error for the full model */
2404 float freg, /* f-statistic for the full regression model */
2405 float rsqr, /* R^2 (coef. of multiple determination) */
2406 float smax, /* signed maximum of signal */
2407 float tmax, /* epoch of signed maximum of signal */
2408 float pmax, /* percentage change due to signal */
2409 float area, /* area between signal and baseline */
2410 float parea, /* percentage area between signal and baseline */
2411
2412 float * rmsreg_vol, /* root-mean-square for the full regression model */
2413 float * freg_vol, /* f-statistic for the full regression model */
2414 float * rsqr_vol, /* R^2 volume data */
2415 float * smax_vol, /* signed maximum of signal volume data */
2416 float * tmax_vol, /* epoch of signed maximum volume data */
2417 float * pmax_vol, /* percentage change in signal volume data */
2418 float * area_vol, /* area underneath signal volume data */
2419 float * parea_vol, /* percentage area underneath signal volume data */
2420 float ** ncoef_vol, /* volume of noise model parameters */
2421 float ** scoef_vol, /* volume of signal model parameters */
2422 float ** tncoef_vol, /* volume of noise model t-statistics */
2423 float ** tscoef_vol, /* volume of signal model t-statistics */
2424 float ** sfit_vol, /* voxelwise 3d+time fitted signal model */
2425 float ** snfit_vol /* voxelwise 3d+time fitted signal+noise model */
2426 )
2427
2428 {
2429 int ip; /* parameter index */
2430 int it; /* time index */
2431 float * s_array; /* fitted signal model time series */
2432 float * n_array; /* fitted noise model time series */
2433
2434
2435 /*----- save regression results into volume data -----*/
2436 if (freg_vol != NULL) freg_vol[iv] = freg;
2437 if (rmsreg_vol != NULL) rmsreg_vol[iv] = rmsreg;
2438 if (rsqr_vol != NULL) rsqr_vol[iv] = rsqr;
2439
2440 /*----- save signed maximum and epoch of signed maximum of signal -----*/
2441 if (smax_vol != NULL) smax_vol[iv] = smax;
2442 if (tmax_vol != NULL) tmax_vol[iv] = tmax;
2443
2444 /*----- save percentage change and area beneath the signal -----*/
2445 if (pmax_vol != NULL) pmax_vol[iv] = pmax;
2446 if (area_vol != NULL) area_vol[iv] = area;
2447 if (parea_vol != NULL) parea_vol[iv] = parea;
2448
2449 /*----- save noise parameter estimates -----*/
2450 for (ip = 0; ip < r; ip++)
2451 {
2452 if (ncoef_vol[ip] != NULL) ncoef_vol[ip][iv] = par_full[ip];
2453 if (tncoef_vol[ip] != NULL) tncoef_vol[ip][iv] = tpar_full[ip];
2454 }
2455
2456 /*----- save signal parameter estimates -----*/
2457 for (ip = 0; ip < p; ip++)
2458 {
2459 if (scoef_vol[ip] != NULL) scoef_vol[ip][iv] = par_full[ip+r];
2460 if (tscoef_vol[ip] != NULL) tscoef_vol[ip][iv] = tpar_full[ip+r];
2461 }
2462
2463
2464 /*----- save fitted signal model time series -----*/
2465 if (sfit_vol != NULL)
2466 {
2467 if (novar)
2468 {
2469 for (it = 0; it < ts_length; it++)
2470 sfit_vol[it][iv] = 0.0;
2471 }
2472 else
2473 {
2474 s_array = (float *) malloc (sizeof(float) * (ts_length));
2475 MTEST (s_array);
2476
2477 smodel (par_full + r, ts_length, x_array, s_array);
2478
2479 for (it = 0; it < ts_length; it++)
2480 sfit_vol[it][iv] = s_array[it];
2481
2482 free (s_array); s_array = NULL;
2483 }
2484 }
2485
2486
2487 /*----- save fitted signal+noise model time series -----*/
2488 if (snfit_vol != NULL)
2489 {
2490 n_array = (float *) malloc (sizeof(float) * (ts_length));
2491 MTEST (n_array);
2492 nmodel (par_full, ts_length, x_array, n_array);
2493
2494 for (it = 0; it < ts_length; it++)
2495 {
2496 snfit_vol[it][iv] = n_array[it];
2497 }
2498
2499 free (n_array); n_array = NULL;
2500
2501
2502 if (! novar)
2503 {
2504 s_array = (float *) malloc (sizeof(float) * (ts_length));
2505 MTEST (s_array);
2506 smodel (par_full + r, ts_length, x_array, s_array);
2507
2508 for (it = 0; it < ts_length; it++)
2509 {
2510 snfit_vol[it][iv] += s_array[it];
2511 }
2512
2513 free (s_array); s_array = NULL;
2514 }
2515 }
2516
2517 }
2518
2519 /*---------------------------------------------------------------------------*/
2520 /*
2521 Routine to write one AFNI data set.
2522
2523 This data set may be either 'fitt' type (intensity + t-statistic)
2524 or 'fift' type (intensity + F-statistic).
2525
2526 The intensity data is in ffim, and the corresponding statistic is in ftr.
2527
2528 numdof = numerator degrees of freedom
2529 dendof = denominator degrees of freedom
2530
2531 Note: if dendof = 0, then write 'fitt' type data set;
2532 if dendof > 0, then write 'fift' type data set.
2533 */
2534
write_afni_data(char * input_filename,int nxyz,char * filename,float * ffim,float * ftr,int numdof,int dendof)2535 void write_afni_data (char * input_filename, int nxyz, char * filename,
2536 float * ffim, float * ftr, int numdof, int dendof)
2537 {
2538 int ii; /* voxel index */
2539 THD_3dim_dataset * dset=NULL; /* input afni data set pointer */
2540 THD_3dim_dataset * new_dset=NULL; /* output afni data set pointer */
2541 int iv; /* sub-brick index */
2542 int ierror; /* number of errors in editing data */
2543 int ibuf[32]; /* integer buffer */
2544 float fbuf[MAX_STAT_AUX]; /* float buffer */
2545 float fimfac; /* scale factor for short data */
2546 /* int output_datum; */ /* data type for output data */
2547 short * tsp; /* 2nd sub-brick data pointer */
2548 void * vdif; /* 1st sub-brick data pointer */
2549 int func_type; /* afni data set type */
2550 float top, func_scale_short; /* parameters for scaling data */
2551 char label[THD_MAX_NAME]; /* label for output file history */
2552 int nbad; /* number of bad voxels in volume */
2553
2554 /*----- read input dataset -----*/
2555 dset = THD_open_dataset (input_filename) ; /* was THD_open_one...*/
2556 if( ! ISVALID_3DIM_DATASET(dset) ){
2557 fprintf(stderr,"*** Unable to open dataset file %s\n",
2558 input_filename);
2559 exit(1) ;
2560 }
2561
2562 /*-- make an empty copy of this dataset, for eventual output --*/
2563 new_dset = EDIT_empty_copy( dset ) ;
2564
2565
2566 /*----- Record history of dataset -----*/
2567 tross_Copy_History( dset , new_dset ) ;
2568 snprintf (label, THD_MAX_NAME, "Output prefix: %s", filename);
2569 if( commandline != NULL )
2570 tross_multi_Append_History( new_dset , commandline,label,NULL ) ;
2571 else
2572 tross_Append_History ( new_dset, label);
2573
2574
2575 iv = DSET_PRINCIPAL_VALUE(dset) ;
2576
2577 fprintf(stderr,"--output datum is %d\n", output_datum);
2578 /* output_datum = DSET_BRICK_TYPE(dset,iv) ;
2579 if( output_datum == MRI_byte ) output_datum = MRI_short ;
2580 */
2581
2582 ibuf[0] = output_datum ; ibuf[1] = MRI_short ;
2583
2584 if (dendof == 0) func_type = FUNC_TT_TYPE;
2585 else func_type = FUNC_FT_TYPE;
2586
2587 ierror = EDIT_dset_items( new_dset ,
2588 ADN_prefix , filename ,
2589 ADN_label1 , filename ,
2590 ADN_self_name , filename ,
2591 ADN_type , ISHEAD(dset) ? HEAD_FUNC_TYPE :
2592 GEN_FUNC_TYPE ,
2593 ADN_func_type , func_type ,
2594 ADN_nvals , FUNC_nvals[func_type] ,
2595 ADN_datum_array , ibuf ,
2596 ADN_ntt , 0 ,
2597 ADN_malloc_type, DATABLOCK_MEM_MALLOC ,
2598 ADN_none ) ;
2599
2600 if( ierror > 0 ){
2601 fprintf(stderr,
2602 "*** %d errors in attempting to create output dataset!\n", ierror ) ;
2603 exit(1) ;
2604 }
2605
2606 if( THD_is_file(new_dset->dblk->diskptr->header_name) ){
2607 fprintf(stderr,
2608 "*** Output dataset file %s already exists--cannot continue!\a\n",
2609 new_dset->dblk->diskptr->header_name ) ;
2610 exit(1) ;
2611 }
2612
2613 /*----- deleting exemplar dataset -----*/
2614 THD_delete_3dim_dataset( dset , False ) ; dset = NULL ;
2615
2616 nbad = thd_floatscan( nxyz , ffim ) ;
2617 if(nbad)
2618 fprintf(stderr,"++ %d bad floating point values in combined dataset\n",nbad);
2619
2620 /*----- allocate memory for output data -----*/
2621 vdif = (void *) malloc( mri_datum_size(output_datum) * nxyz );
2622 if (vdif == NULL) NLfit_error ("Unable to allocate memory for vdif");
2623 tsp = (short *) malloc( sizeof(short) * nxyz );
2624 if (tsp == NULL) NLfit_error ("Unable to allocate memory for tsp");
2625
2626 /*----- attach bricks to new data set -----*/
2627 mri_fix_data_pointer (vdif, DSET_BRICK(new_dset,0));
2628 mri_fix_data_pointer (tsp, DSET_BRICK(new_dset,1));
2629
2630
2631 /*----- convert data type to output specification -----*/
2632 fimfac = EDIT_coerce_autoscale_new (nxyz,
2633 MRI_float, ffim,
2634 output_datum, vdif);
2635 if (fimfac != 0.0) fimfac = 1.0 / fimfac;
2636
2637 #define TOP_SS 32700
2638
2639 if (dendof == 0) /* t-statistic */
2640 {
2641 top = TOP_SS/FUNC_TT_SCALE_SHORT;
2642 func_scale_short = FUNC_TT_SCALE_SHORT;
2643 }
2644 else /* F-statistic */
2645 {
2646 top = TOP_SS/FUNC_FT_SCALE_SHORT;
2647 func_scale_short = FUNC_FT_SCALE_SHORT;
2648 }
2649
2650 for (ii = 0; ii < nxyz; ii++)
2651 {
2652 if (ftr[ii] > top)
2653 tsp[ii] = TOP_SS;
2654 else if (ftr[ii] < -top)
2655 tsp[ii] = -TOP_SS;
2656 else if (ftr[ii] >= 0.0)
2657 tsp[ii] = (short) (func_scale_short * ftr[ii] + 0.5);
2658 else
2659 tsp[ii] = (short) (func_scale_short * ftr[ii] - 0.5);
2660
2661 }
2662
2663
2664 /*----- write afni data set -----*/
2665 INFO_message("Writing combined dataset into %s\n",DSET_BRIKNAME(new_dset)) ;
2666
2667 fbuf[0] = numdof;
2668 fbuf[1] = dendof;
2669 for( ii=2 ; ii < MAX_STAT_AUX ; ii++ ) fbuf[ii] = 0.0 ;
2670 (void) EDIT_dset_items( new_dset , ADN_stat_aux , fbuf , ADN_none ) ;
2671
2672 fbuf[0] = (output_datum == MRI_short && fimfac != 1.0 ) ? fimfac : 0.0 ;
2673 fbuf[1] = 1.0 / func_scale_short ;
2674 (void) EDIT_dset_items( new_dset , ADN_brick_fac , fbuf , ADN_none ) ;
2675
2676 if( do_FDR && !AFNI_noenv("AFNI_AUTOMATIC_FDR") )
2677 { int ii = THD_create_all_fdrcurves( new_dset ) ;
2678 if( ii > 0 ) ININFO_message("created %d FDR curves in header",ii) ;
2679 }
2680
2681 THD_load_statistics( new_dset ) ;
2682 THD_write_3dim_dataset( NULL,NULL , new_dset , True ) ;
2683
2684 /*----- deallocate memory -----*/
2685 THD_delete_3dim_dataset( new_dset , False ) ; new_dset = NULL ;
2686
2687 }
2688
2689
2690 /*---------------------------------------------------------------------------*/
2691 /*
2692 Routine to write one bucket data set.
2693 */
2694
write_bucket_data(int q,int p,int numdof,int dendof,int nxyz,int n,float * rmsreg_vol,float * freg_vol,float * rsqr_vol,float * smax_vol,float * tmax_vol,float * pmax_vol,float * area_vol,float * parea_vol,float ** ncoef_vol,float ** scoef_vol,float ** tncoef_vol,float ** tscoef_vol,char * input_filename,NL_options * option_data)2695 void write_bucket_data
2696 (
2697 int q, /* number of parameters in the noise model */
2698 int p, /* number of parameters in the signal model */
2699 int numdof, /* numerator dof for F-statistic */
2700 int dendof, /* denominator dof for F-statistic */
2701 int nxyz, /* number of voxels in image */
2702 int n, /* length of time series data */
2703
2704 float * rmsreg_vol, /* root-mean-square error for the full model */
2705 float * freg_vol, /* f-statistic for the full regression model */
2706 float * rsqr_vol, /* R^2 volume data */
2707 float * smax_vol, /* signed maximum of signal volume data */
2708 float * tmax_vol, /* epoch of signed maximum volume data */
2709 float * pmax_vol, /* percentage change in signal volume data */
2710 float * area_vol, /* area underneath signal volume data */
2711 float * parea_vol, /* percentage area underneath signal volume data */
2712 float ** ncoef_vol, /* volume of noise model parameters */
2713 float ** scoef_vol, /* volume of signal model parameters */
2714 float ** tncoef_vol, /* volume of noise model t-statistics */
2715 float ** tscoef_vol, /* volume of signal model t-statistics */
2716
2717 char * input_filename,
2718
2719 NL_options * option_data /* user input options */
2720 )
2721
2722 {
2723 const float EPSILON = 1.0e-10;
2724
2725 THD_3dim_dataset * old_dset = NULL; /* prototype dataset */
2726 THD_3dim_dataset * new_dset = NULL; /* output bucket dataset */
2727 char * output_prefix; /* prefix name for bucket dataset */
2728 char * output_session; /* directory for bucket dataset */
2729 int nbricks, ib; /* number of sub-bricks in bucket dataset */
2730 short ** bar = NULL; /* bar[ib] points to data for sub-brick #ib */
2731 float ** far = NULL; /* far[ib] points to data for sub-brick #ib */
2732 float factor; /* factor is new scale factor for sub-brick #ib */
2733 int brick_type; /* indicates statistical type of sub-brick */
2734 int brick_coef; /* regression coefficient index for sub-brick */
2735 char * brick_label; /* character string label for sub-brick */
2736 int ierror; /* number of errors in editing data */
2737 float *volume=NULL; /* volume of floating point data */
2738 int dimension; /* dimension of full model = p + q */
2739 char label[THD_MAX_NAME]; /* label for output file history */
2740 void * imptr; /* pointer to volume in correct datum type to actually write out*/
2741 int nbad; /* number of bad floating point values in volume */
2742
2743 /*----- initialize local variables -----*/
2744 nbricks = option_data->numbricks;
2745 output_prefix = option_data->bucket_filename;
2746 output_session = (char *) malloc (sizeof(char) * THD_MAX_NAME);
2747 strcpy (output_session, "./");
2748 dimension = p + q;
2749
2750
2751 /*----- allocate memory -----*/
2752 if(output_datum==MRI_short) {
2753 bar = (short **) malloc (sizeof(short *) * nbricks);
2754 MTEST (bar);
2755 }
2756 else {
2757 far = (float **) malloc (sizeof(float *) * nbricks);
2758 MTEST (far);
2759 }
2760
2761 /*----- read first dataset -----*/
2762 old_dset = THD_open_dataset (input_filename); /* was THD_open_one...*/
2763
2764
2765 /*-- make an empty copy of this dataset, for eventual output --*/
2766 new_dset = EDIT_empty_copy (old_dset);
2767
2768
2769 /*----- Record history of dataset -----*/
2770 tross_Copy_History( old_dset , new_dset ) ;
2771 if( commandline != NULL )
2772 tross_Append_History( new_dset , commandline ) ;
2773 snprintf (label, THD_MAX_NAME, "Output prefix: %s", output_prefix);
2774 tross_Append_History ( new_dset, label);
2775
2776
2777 /*----- Modify some structural properties. Note that the nbricks
2778 just make empty sub-bricks, without any data attached. -----*/
2779 ierror = EDIT_dset_items (new_dset,
2780 ADN_prefix, output_prefix,
2781 ADN_directory_name, output_session,
2782 ADN_type, HEAD_FUNC_TYPE,
2783 ADN_func_type, FUNC_BUCK_TYPE,
2784 ADN_ntt, 0, /* no time */
2785 ADN_nvals, nbricks,
2786 ADN_malloc_type, DATABLOCK_MEM_MALLOC ,
2787 ADN_none ) ;
2788
2789 if( ierror > 0 )
2790 {
2791 fprintf(stderr,
2792 "*** %d errors in attempting to create output dataset!\n",
2793 ierror);
2794 exit(1);
2795 }
2796
2797 if (THD_is_file(DSET_HEADNAME(new_dset)))
2798 {
2799 fprintf(stderr,
2800 "*** Output dataset file %s already exists--cannot continue!\n",
2801 DSET_HEADNAME(new_dset));
2802 exit(1);
2803 }
2804
2805
2806 /*----- deleting exemplar dataset -----*/
2807 THD_delete_3dim_dataset( old_dset , False ); old_dset = NULL ;
2808
2809 nbad = 0;
2810
2811 /*----- loop over new sub-brick index, attach data array with
2812 EDIT_substitute_brick then put some strings into the labels and
2813 keywords, and modify the sub-brick scaling factor -----*/
2814 for (ib = 0; ib < nbricks; ib++)
2815 {
2816 /*----- get information about this sub-brick -----*/
2817 brick_type = option_data->brick_type[ib];
2818 brick_coef = option_data->brick_coef[ib];
2819 brick_label = option_data->brick_label[ib];
2820
2821 if (brick_type == FUNC_FIM_TYPE)
2822 {
2823 if (brick_coef < q)
2824 volume = ncoef_vol[brick_coef];
2825 else if (brick_coef < p+q)
2826 volume = scoef_vol[brick_coef-q];
2827 else if (brick_coef == p+q)
2828 volume = tmax_vol;
2829 else if (brick_coef == p+q+1)
2830 volume = smax_vol;
2831 else if (brick_coef == p+q+2)
2832 volume = pmax_vol;
2833 else if (brick_coef == p+q+3)
2834 volume = area_vol;
2835 else if (brick_coef == p+q+4)
2836 volume = parea_vol;
2837 else if (brick_coef == p+q+5)
2838 volume = rmsreg_vol;
2839 else if (brick_coef == p+q+6)
2840 volume = rsqr_vol;
2841 }
2842 else if (brick_type == FUNC_TT_TYPE)
2843 {
2844 if (brick_coef < q)
2845 volume = tncoef_vol[brick_coef];
2846 else if (brick_coef < p+q)
2847 volume = tscoef_vol[brick_coef-q];
2848 EDIT_BRICK_TO_FITT (new_dset, ib, dendof);
2849 }
2850 else if (brick_type == FUNC_FT_TYPE)
2851 {
2852 volume = freg_vol;
2853 EDIT_BRICK_TO_FIFT (new_dset, ib, numdof, dendof);
2854 }
2855
2856 nbad += thd_floatscan( nxyz , volume ) ;
2857
2858 /*----- allocate memory for output sub-brick -----*/
2859 if(output_datum==MRI_short) {
2860 bar[ib] = (short *) malloc (sizeof(short) * nxyz);
2861 MTEST (bar[ib]);
2862 factor = EDIT_coerce_autoscale_new (nxyz, MRI_float, volume,
2863 MRI_short, bar[ib]);
2864 imptr = bar[ib];
2865 }
2866 else {
2867 far[ib] = (float *) malloc (sizeof(float) * nxyz);
2868 MTEST (far[ib]);
2869 factor = EDIT_coerce_autoscale_new (nxyz, MRI_float, volume,
2870 output_datum, far[ib]);
2871 imptr = far[ib];
2872 }
2873 if (factor < EPSILON) factor = 0.0;
2874 else factor = 1.0 / factor;
2875
2876 if( output_datum == MRI_short )
2877 EDIT_misfit_report( DSET_FILECODE(new_dset) , ib ,
2878 nxyz , factor , imptr , volume ) ;
2879
2880
2881 /*----- edit the sub-brick -----*/
2882 EDIT_BRICK_LABEL (new_dset, ib, brick_label);
2883 EDIT_BRICK_FACTOR (new_dset, ib, factor);
2884
2885 /*----- attach image pointer to be sub-brick #ib -----*/
2886 EDIT_substitute_brick (new_dset, ib, output_datum, imptr);
2887
2888 }
2889
2890 if(nbad)
2891 fprintf(stderr,"++ %d bad floating point values in dataset\n",nbad);
2892
2893 /*----- write bucket data set -----*/
2894
2895 INFO_message("Writing bucket dataset: %s",DSET_BRIKNAME(new_dset)) ;
2896
2897 if( do_FDR && !AFNI_noenv("AFNI_AUTOMATIC_FDR") )
2898 { int ii = THD_create_all_fdrcurves( new_dset ) ;
2899 if( ii > 0 ) ININFO_message("created %d FDR curves in header",ii) ;
2900 }
2901
2902 THD_load_statistics (new_dset);
2903 THD_write_3dim_dataset( NULL,NULL , new_dset , True ) ;
2904
2905 /*----- deallocate memory -----*/
2906 /* if(output_datum==MRI_short) {*/
2907 THD_delete_3dim_dataset( new_dset , False ) ; new_dset = NULL ;
2908 /*} */
2909
2910 }
2911
2912
2913 /*---------------------------------------------------------------------------*/
2914 /*
2915 Routine to write one AFNI 3d+time data set.
2916 */
2917
2918
write_3dtime(char * input_filename,int ts_length,float ** vol_array,char * output_filename)2919 void write_3dtime
2920 (
2921 char * input_filename, /* file name of input 3d+time dataset */
2922 int ts_length, /* length of time series data */
2923 float ** vol_array, /* output time series volume data */
2924 char * output_filename /* output afni data set file name */
2925 )
2926
2927 {
2928 const float EPSILON = 1.0e-10;
2929
2930 THD_3dim_dataset * dset = NULL; /* input afni data set pointer */
2931 THD_3dim_dataset * new_dset = NULL; /* output afni data set pointer */
2932 int ib; /* sub-brick index */
2933 int ierror; /* number of errors in editing data */
2934 int nxyz; /* total number of voxels */
2935 float factor; /* factor is new scale factor for sub-brick #ib */
2936 short ** bar = NULL; /* bar[ib] points to data for sub-brick #ib */
2937 float ** far = NULL; /* far[ib] points to data for sub-brick #ib */
2938 float * fbuf; /* float buffer */
2939 float * volume; /* pointer to volume of data */
2940 char label[THD_MAX_NAME]; /* label for output file history */
2941 int nbad; /* number of voxels in volume with bad floating point values */
2942
2943 /*----- Initialize local variables -----*/
2944 dset = THD_open_dataset (input_filename); /* was THD_open_one...*/
2945 nxyz = dset->daxes->nxx * dset->daxes->nyy * dset->daxes->nzz;
2946
2947
2948 /*----- allocate memory -----*/
2949 if(output_datum==MRI_short) {
2950 bar = (short **) malloc (sizeof(short *) * ts_length); MTEST (bar);
2951 }
2952 else {
2953 far = (float **) malloc (sizeof(float *) * ts_length); MTEST (far);
2954 }
2955
2956 fbuf = (float *) malloc (sizeof(float) * ts_length); MTEST (fbuf);
2957 for (ib = 0; ib < ts_length; ib++) fbuf[ib] = 0.0;
2958
2959
2960 /*-- make an empty copy of the prototype dataset, for eventual output --*/
2961 new_dset = EDIT_empty_copy (dset);
2962
2963
2964 /*----- Record history of dataset -----*/
2965 tross_Copy_History( dset , new_dset ) ;
2966 if( commandline != NULL )
2967 tross_Append_History( new_dset , commandline ) ;
2968 snprintf (label, THD_MAX_NAME, "Output prefix: %s", output_filename);
2969 tross_Append_History ( new_dset, label);
2970
2971
2972 /*----- delete prototype dataset -----*/
2973 THD_delete_3dim_dataset (dset, False); dset = NULL ;
2974
2975
2976 ierror = EDIT_dset_items (new_dset,
2977 ADN_prefix, output_filename,
2978 ADN_label1, output_filename,
2979 ADN_self_name, output_filename,
2980 ADN_malloc_type, DATABLOCK_MEM_MALLOC,
2981 ADN_datum_all, output_datum,
2982 ADN_nvals, ts_length,
2983 ADN_ntt, ts_length,
2984 ADN_none);
2985
2986 if( ierror > 0 ){
2987 fprintf(stderr,
2988 "*** %d errors in attempting to create output dataset!\n", ierror ) ;
2989 exit(1) ;
2990 }
2991
2992 if( THD_is_file(new_dset->dblk->diskptr->header_name) ){
2993 fprintf(stderr,
2994 "*** Output dataset file %s already exists--cannot continue!\a\n",
2995 new_dset->dblk->diskptr->header_name ) ;
2996 exit(1) ;
2997 }
2998
2999 nbad = 0;
3000
3001 /*----- attach bricks to new data set -----*/
3002 for (ib = 0; ib < ts_length; ib++)
3003 {
3004
3005 /*----- Set pointer to appropriate volume -----*/
3006 volume = vol_array[ib];
3007 nbad += thd_floatscan( nxyz , volume ) ;
3008 if(output_datum==MRI_short) {
3009 /*----- Allocate memory for output sub-brick -----*/
3010 bar[ib] = (short *) malloc (sizeof(short) * nxyz);
3011 MTEST (bar[ib]);
3012
3013 /*----- Convert data type to short for this sub-brick -----*/
3014 factor = EDIT_coerce_autoscale_new (nxyz, MRI_float, volume,
3015 output_datum, bar[ib]);
3016 if (factor < EPSILON) factor = 0.0;
3017 else factor = 1.0 / factor;
3018 fbuf[ib] = factor;
3019 /*----- attach bar[ib] to be sub-brick #ib -----*/
3020 mri_fix_data_pointer (bar[ib], DSET_BRICK(new_dset,ib));
3021
3022 EDIT_misfit_report( DSET_FILECODE(new_dset) , ib ,
3023 nxyz , factor , bar[ib] , volume ) ;
3024 }
3025 else {
3026 /*----- Allocate memory for output sub-brick -----*/
3027 far[ib] = (float *) malloc (sizeof(float) * nxyz);
3028 MTEST (far[ib]);
3029
3030 /*----- Convert data type to float for this sub-brick -----*/
3031 factor = EDIT_coerce_autoscale_new (nxyz, MRI_float, volume,
3032 output_datum, far[ib]);
3033 if (factor < EPSILON) factor = 0.0;
3034 else factor = 1.0 / factor;
3035 fbuf[ib] = factor;
3036 /*----- attach far[ib] to be sub-brick #ib -----*/
3037 mri_fix_data_pointer (far[ib], DSET_BRICK(new_dset,ib));
3038 }
3039 }
3040 if(nbad)
3041 WARNING_message("%d bad floating point values in dataset\n",nbad);
3042
3043 /*----- write afni data set -----*/
3044
3045 (void) EDIT_dset_items( new_dset , ADN_brick_fac , fbuf , ADN_none ) ;
3046
3047 INFO_message("Writing 3D+time dataset %s\n",DSET_BRIKNAME(new_dset)) ;
3048
3049 THD_load_statistics (new_dset);
3050 THD_write_3dim_dataset (NULL, NULL, new_dset, True);
3051
3052
3053 /*----- deallocate memory -----*/
3054 /* if(output_datum==MRI_short) {*/
3055 THD_delete_3dim_dataset (new_dset, False); new_dset = NULL ;
3056 free (fbuf); fbuf = NULL;
3057 /* }*/
3058 }
3059
3060
3061 /*---------------------------------------------------------------------------*/
3062 /*
3063 Write out the user requested output files.
3064 */
3065
output_results(int r,int p,float * min_nconstr,float * max_nconstr,float * min_sconstr,float * max_sconstr,int nxyz,int ts_length,float * rmsreg_vol,float * freg_vol,float * rsqr_vol,float * smax_vol,float * tmax_vol,float * pmax_vol,float * area_vol,float * parea_vol,float ** ncoef_vol,float ** scoef_vol,float ** tncoef_vol,float ** tscoef_vol,float ** sfit_vol,float ** snfit_vol,char * input_filename,char * freg_filename,char * frsqr_filename,char * fsmax_filename,char * ftmax_filename,char * fpmax_filename,char * farea_filename,char * fparea_filename,char ** fncoef_filename,char ** fscoef_filename,char ** tncoef_filename,char ** tscoef_filename,char * sfit_filename,char * snfit_filename,NL_options * option_data)3066 void output_results
3067 (
3068 int r, /* number of parameters in the noise model */
3069 int p, /* number of parameters in the signal model */
3070 float * min_nconstr, /* minimum parameter constraints for noise model */
3071 float * max_nconstr, /* maximum parameter constraints for noise model */
3072 float * min_sconstr, /* minimum parameter constraints for signal model */
3073 float * max_sconstr, /* maximum parameter constraints for signal model */
3074 int nxyz, /* number of voxels in image */
3075 int ts_length, /* length of time series data */
3076 float * rmsreg_vol, /* root-mean-square error for the full model */
3077 float * freg_vol, /* f-statistic for the full regression model */
3078 float * rsqr_vol, /* R^2 volume data */
3079 float * smax_vol, /* signed maximum of signal volume data */
3080 float * tmax_vol, /* epoch of signed maximum volume data */
3081 float * pmax_vol, /* percentage change in signal volume data */
3082 float * area_vol, /* area underneath signal volume data */
3083 float * parea_vol, /* percentage area underneath signal volume data */
3084 float ** ncoef_vol, /* volume of noise model parameters */
3085 float ** scoef_vol, /* volume of signal model parameters */
3086 float ** tncoef_vol, /* volume of noise model t-statistics */
3087 float ** tscoef_vol, /* volume of signal model t-statistics */
3088 float ** sfit_vol, /* voxelwise 3d+time fitted signal model */
3089 float ** snfit_vol, /* voxelwise 3d+time fitted signal+noise model */
3090
3091 char * input_filename, /* file name of input 3d+time dataset */
3092 char * freg_filename, /* file name for f-statistics */
3093 char * frsqr_filename, /* file name for R^2 statistics */
3094 char * fsmax_filename, /* file name for signal signed maximum */
3095 char * ftmax_filename, /* file name for time of signed maximum */
3096 char * fpmax_filename, /* file name for percentage signal change */
3097 char * farea_filename, /* file name for area underneath signal */
3098 char * fparea_filename, /* file name for % area underneath signal */
3099 char ** fncoef_filename, /* file name for noise model parameters */
3100 char ** fscoef_filename, /* file name for signal model parameters */
3101 char ** tncoef_filename, /* file name for noise model t-statistics */
3102 char ** tscoef_filename, /* file name for signal model t-statistics */
3103 char * sfit_filename, /* file name for 3d+time fitted signal model */
3104 char * snfit_filename, /* file name for 3d+time fitted signal+noise */
3105
3106 NL_options * option_data /* user input options */
3107
3108 )
3109
3110 {
3111 int ip; /* parameter index */
3112 int dimension; /* dimension of full model = r + p */
3113 int numdof, dendof; /* numerator and denominator degrees of freedom */
3114
3115
3116 /*----- Initialization -----*/
3117 dimension = r + p;
3118 numdof = p;
3119 dendof = ts_length - dimension;
3120
3121
3122 /*----- Adjust dof if constraints force a parameter to be a constant -----*/
3123 for (ip = 0; ip < r; ip++)
3124 if (min_nconstr[ip] == max_nconstr[ip])
3125 dendof++;
3126
3127 for (ip = 0; ip < p; ip++)
3128 if (min_sconstr[ip] == max_sconstr[ip])
3129 {
3130 numdof--;
3131 dendof++;
3132 }
3133
3134
3135 /*----- write the bucket dataset -----*/
3136 if (option_data->numbricks > 0)
3137 write_bucket_data (r, p, numdof, dendof, nxyz, ts_length, rmsreg_vol,
3138 freg_vol, rsqr_vol, smax_vol, tmax_vol, pmax_vol,
3139 area_vol, parea_vol, ncoef_vol, scoef_vol,
3140 tncoef_vol, tscoef_vol, input_filename, option_data);
3141
3142
3143 /*----- write out the f-statistics for the regression -----*/
3144 if (freg_filename != NULL)
3145 write_afni_data (input_filename, nxyz, freg_filename,
3146 rmsreg_vol, freg_vol, numdof, dendof);
3147
3148
3149 /*----- write out the R^2 (coefficient of multiple determination) -----*/
3150 if (frsqr_filename != NULL)
3151 write_afni_data (input_filename, nxyz, frsqr_filename,
3152 rsqr_vol, freg_vol, numdof, dendof);
3153
3154
3155 /*----- write out the signed maximum of signal estimates -----*/
3156 if (fsmax_filename != NULL)
3157 write_afni_data (input_filename, nxyz, fsmax_filename,
3158 smax_vol, freg_vol, numdof, dendof);
3159
3160
3161 /*----- write out the epoch of signed maximum of signal estimates -----*/
3162 if (ftmax_filename != NULL)
3163 write_afni_data (input_filename, nxyz, ftmax_filename,
3164 tmax_vol, freg_vol, numdof, dendof);
3165
3166
3167 /*----- write out the max. percentage change due to signal -----*/
3168 if (fpmax_filename != NULL)
3169 write_afni_data (input_filename, nxyz, fpmax_filename,
3170 pmax_vol, freg_vol, numdof, dendof);
3171
3172
3173 /*----- write out the area between the signal and baseline -----*/
3174 if (farea_filename != NULL)
3175 write_afni_data (input_filename, nxyz, farea_filename,
3176 area_vol, freg_vol, numdof, dendof);
3177
3178
3179 /*----- write out the percentage area between the signal and baseline -----*/
3180 if (fparea_filename != NULL)
3181 write_afni_data (input_filename, nxyz, fparea_filename,
3182 parea_vol, freg_vol, numdof, dendof);
3183
3184
3185 /*----- write noise model parameters -----*/
3186 for (ip = 0; ip < r; ip++)
3187 {
3188 if (tncoef_filename[ip] != NULL)
3189 write_afni_data (input_filename, nxyz, tncoef_filename[ip],
3190 ncoef_vol[ip], tncoef_vol[ip], dendof, 0);
3191
3192 if (fncoef_filename[ip] != NULL)
3193 write_afni_data (input_filename, nxyz, fncoef_filename[ip],
3194 ncoef_vol[ip], freg_vol, numdof, dendof);
3195 }
3196
3197
3198 /*----- write signal model parameters -----*/
3199 for (ip = 0; ip < p; ip++)
3200 {
3201 if (tscoef_filename[ip] != NULL)
3202 write_afni_data (input_filename, nxyz, tscoef_filename[ip],
3203 scoef_vol[ip], tscoef_vol[ip], dendof, 0);
3204
3205 if (fscoef_filename[ip] != NULL)
3206 write_afni_data (input_filename, nxyz, fscoef_filename[ip],
3207 scoef_vol[ip], freg_vol, numdof, dendof);
3208 }
3209
3210
3211 /*----- write the fitted 3d+time signal model -----*/
3212 if (sfit_filename != NULL)
3213 {
3214 write_3dtime (input_filename, ts_length, sfit_vol, sfit_filename);
3215 }
3216
3217
3218 /*----- write the fitted 3d+time signal+noise model -----*/
3219 if (snfit_filename != NULL)
3220 {
3221 write_3dtime (input_filename, ts_length, snfit_vol, snfit_filename);
3222 }
3223
3224 }
3225
3226
3227 /*---------------------------------------------------------------------------*/
3228 /*
3229 Release all allocated memory space.
3230 */
3231
terminate_program(int r,int p,int ts_length,float *** x_array,float ** ts_array,char ** nname,char *** npname,float ** par_rdcd,float ** min_nconstr,float ** max_nconstr,char ** sname,char *** spname,float ** par_full,float ** tpar_full,float ** min_sconstr,float ** max_sconstr,float ** rmsreg_vol,float ** freg_vol,float ** rsqr_vol,float ** smax_vol,float ** tmax_vol,float ** pmax_vol,float ** area_vol,float ** parea_vol,float *** ncoef_vol,float *** scoef_vol,float *** tncoef_vol,float *** tscoef_vol,float *** sfit_vol,float *** snfit_vol,char ** input_filename,char ** freg_filename,char ** frsqr_filename,char ** fsmax_filename,char ** ftmax_filename,char ** fpmax_filename,char ** farea_filename,char ** fparea_filename,char *** fncoef_filename,char *** fscoef_filename,char *** tncoef_filename,char *** tscoef_filename,char ** sfit_filename,char ** snfit_filename)3232 void terminate_program
3233 (
3234 int r, /* number of parameters in the noise model */
3235 int p, /* number of parameters in the signal model */
3236 int ts_length, /* length of time series data */
3237 float *** x_array, /* independent variable matrix */
3238 float ** ts_array, /* input time series array */
3239 char ** nname, /* noise model name */
3240 char *** npname, /* noise parameter labels */
3241 float ** par_rdcd, /* estimated parameters for the reduced model */
3242 float ** min_nconstr, /* min parameter constraints for noise model */
3243 float ** max_nconstr, /* max parameter constraints for noise model */
3244 char ** sname, /* signal model name */
3245 char *** spname, /* signal parameter labels */
3246 float ** par_full, /* estimated parameters for the full model */
3247 float ** tpar_full, /* t-statistic of parameters in full model */
3248 float ** min_sconstr, /* min parameter constraints for signal model */
3249 float ** max_sconstr, /* max parameter constraints for signal model */
3250
3251 float ** rmsreg_vol, /* rms error for the full regression model */
3252 float ** freg_vol, /* f-statistic for the full regression model */
3253 float ** rsqr_vol, /* R^2 volume data */
3254 float ** smax_vol, /* signed max. of signal volume data */
3255 float ** tmax_vol, /* epoch of signed max. volume data */
3256 float ** pmax_vol, /* percentage change in signal volume data */
3257 float ** area_vol, /* area underneath signal volume data */
3258 float ** parea_vol, /* percent area underneath signal volume data */
3259 float *** ncoef_vol, /* noise model parameters volume data */
3260 float *** scoef_vol, /* signal model parameters volume data */
3261 float *** tncoef_vol, /* noise model t-statistics volume data */
3262 float *** tscoef_vol, /* signal model t-statistics volume data */
3263 float *** sfit_vol, /* voxelwise 3d+time fitted signal model */
3264 float *** snfit_vol, /* voxelwise 3d+time fitted signal+noise */
3265
3266 char ** input_filename, /* file name of input 3d+time dataset */
3267 char ** freg_filename, /* file name for regression f-statistics */
3268 char ** frsqr_filename, /* file name for R^2 statistics */
3269 char ** fsmax_filename, /* file name for signal signed maximum */
3270 char ** ftmax_filename, /* file name for epoch of signed maximum */
3271 char ** fpmax_filename, /* file name for percentage signal change */
3272 char ** farea_filename, /* file name for area underneath signal */
3273 char ** fparea_filename, /* file name for % area underneath signal */
3274 char *** fncoef_filename, /* file name for noise model parameters */
3275 char *** fscoef_filename, /* file name for signal model parameters */
3276 char *** tncoef_filename, /* file name for noise model t-statistics */
3277 char *** tscoef_filename, /* file name for signal model t-statistics */
3278 char ** sfit_filename, /* file name for 3d+time fitted signal model */
3279 char ** snfit_filename /* file name for 3d+time fitted signal+noise */
3280 )
3281
3282 {
3283 int ip; /* parameter index */
3284 int it; /* time index */
3285
3286
3287 /*----- deallocate space for model names and parameters labels -----*/
3288 if (*nname != NULL)
3289 { free (*nname); *nname = NULL; }
3290 if (*sname != NULL)
3291 { free (*sname); *sname = NULL; }
3292 for (ip = 0; ip < MAX_PARAMETERS; ip++)
3293 {
3294 if ((*npname)[ip] != NULL)
3295 { free ((*npname)[ip]); (*npname)[ip] = NULL; }
3296 if ((*spname)[ip] != NULL)
3297 { free ((*spname)[ip]); (*spname)[ip] = NULL; }
3298 }
3299
3300 if (*npname != NULL)
3301 { free (*npname); *npname = NULL; }
3302 if (*spname != NULL)
3303 { free (*spname); *spname = NULL; }
3304
3305
3306 /*----- deallocate memory for parameter constraints -----*/
3307 if (*min_nconstr != NULL) { free (*min_nconstr); *min_nconstr = NULL; }
3308 if (*max_nconstr != NULL) { free (*max_nconstr); *max_nconstr = NULL; }
3309 if (*min_sconstr != NULL) { free (*min_sconstr); *min_sconstr = NULL; }
3310 if (*max_sconstr != NULL) { free (*max_sconstr); *max_sconstr = NULL; }
3311
3312
3313 /*----- deallocate memory space for filenames -----*/
3314 if (*input_filename != NULL)
3315 { free (*input_filename); *input_filename = NULL; }
3316 if (*freg_filename != NULL)
3317 { free (*freg_filename); *freg_filename = NULL; }
3318 if (*frsqr_filename != NULL)
3319 { free (*frsqr_filename); *frsqr_filename = NULL; }
3320 if (*fsmax_filename != NULL)
3321 { free (*fsmax_filename); *fsmax_filename = NULL; }
3322 if (*ftmax_filename != NULL)
3323 { free (*ftmax_filename); *ftmax_filename = NULL; }
3324 if (*fpmax_filename != NULL)
3325 { free (*fpmax_filename); *fpmax_filename = NULL; }
3326 if (*farea_filename != NULL)
3327 { free (*farea_filename); *farea_filename = NULL; }
3328 if (*fparea_filename != NULL)
3329 { free (*fparea_filename); *fparea_filename = NULL; }
3330
3331 for (ip = 0; ip < MAX_PARAMETERS; ip++)
3332 {
3333 if ((*fncoef_filename)[ip] != NULL)
3334 { free ((*fncoef_filename)[ip]); (*fncoef_filename)[ip] = NULL; }
3335 if ((*fscoef_filename)[ip] != NULL)
3336 { free ((*fscoef_filename)[ip]); (*fscoef_filename)[ip] = NULL; }
3337 if ((*tncoef_filename)[ip] != NULL)
3338 { free ((*tncoef_filename)[ip]); (*tncoef_filename)[ip] = NULL; }
3339 if ((*tscoef_filename)[ip] != NULL)
3340 { free ((*tscoef_filename)[ip]); (*tscoef_filename)[ip] = NULL; }
3341 }
3342
3343 if (*fncoef_filename != NULL)
3344 { free (*fncoef_filename); *fncoef_filename = NULL; }
3345 if (*fscoef_filename != NULL)
3346 { free (*fscoef_filename); *fscoef_filename = NULL; }
3347 if (*tncoef_filename != NULL)
3348 { free (*tncoef_filename); *tncoef_filename = NULL; }
3349 if (*tscoef_filename != NULL)
3350 { free (*tscoef_filename); *tscoef_filename = NULL; }
3351
3352 if (*sfit_filename != NULL)
3353 { free (*sfit_filename); *sfit_filename = NULL; }
3354 if (*snfit_filename != NULL)
3355 { free (*snfit_filename); *snfit_filename = NULL; }
3356
3357
3358 /*----- deallocate space for input time series -----*/
3359 if (*ts_array != NULL) { free (*ts_array); *ts_array = NULL; }
3360
3361
3362 /*----- deallocate space for independent variable matrix -----*/
3363 for (it = 0; it < ts_length; it++)
3364 if ((*x_array)[it] != NULL)
3365 { free ((*x_array)[it]); (*x_array)[it] = NULL; }
3366 if (*x_array != NULL) { free (*x_array); *x_array = NULL; }
3367
3368
3369 /*----- deallocate space for parameters -----*/
3370 if (*par_rdcd != NULL) { free (*par_rdcd); *par_rdcd = NULL; }
3371 if (*par_full != NULL) { free (*par_full); *par_full = NULL; }
3372 if (*tpar_full != NULL) { free (*tpar_full); *tpar_full = NULL; }
3373
3374
3375 /*----- deallocate space for volume data -----*/
3376 if (*rmsreg_vol != NULL) { free (*rmsreg_vol); *rmsreg_vol = NULL; }
3377 if (*freg_vol != NULL) { free (*freg_vol); *freg_vol = NULL; }
3378 if (*rsqr_vol != NULL) { free (*rsqr_vol); *rsqr_vol = NULL; }
3379 if (*smax_vol != NULL) { free (*smax_vol); *smax_vol = NULL; }
3380 if (*tmax_vol != NULL) { free (*tmax_vol); *tmax_vol = NULL; }
3381 if (*pmax_vol != NULL) { free (*pmax_vol); *pmax_vol = NULL; }
3382 if (*area_vol != NULL) { free (*area_vol); *area_vol = NULL; }
3383 if (*parea_vol != NULL) { free (*parea_vol); *parea_vol = NULL; }
3384
3385 if (*ncoef_vol != NULL)
3386 {
3387 for (ip = 0; ip < r; ip++)
3388 if ((*ncoef_vol)[ip] != NULL)
3389 { free ((*ncoef_vol)[ip]); (*ncoef_vol)[ip] = NULL; }
3390 free (*ncoef_vol); *ncoef_vol = NULL;
3391 }
3392
3393 if (*tncoef_vol != NULL)
3394 {
3395 for (ip = 0; ip < r; ip++)
3396 if ((*tncoef_vol)[ip] != NULL)
3397 { free ((*tncoef_vol)[ip]); (*tncoef_vol)[ip] = NULL; }
3398 free (*tncoef_vol); *tncoef_vol = NULL;
3399 }
3400
3401 if (*scoef_vol != NULL)
3402 {
3403 for (ip = 0; ip < p; ip++)
3404 if ((*scoef_vol)[ip] != NULL)
3405 { free ((*scoef_vol)[ip]); (*scoef_vol)[ip] = NULL; }
3406 free (*scoef_vol); *scoef_vol = NULL;
3407 }
3408
3409 if (*tscoef_vol != NULL)
3410 {
3411 for (ip = 0; ip < p; ip++)
3412 if ((*tscoef_vol)[ip] != NULL)
3413 { free ((*tscoef_vol)[ip]); (*tscoef_vol)[ip] = NULL; }
3414 free (*tscoef_vol); *tscoef_vol = NULL;
3415 }
3416
3417 if (*sfit_vol != NULL)
3418 {
3419 for (it = 0; it < ts_length; it++)
3420 if ((*sfit_vol)[it] != NULL)
3421 { free ((*sfit_vol)[it]); (*sfit_vol)[it] = NULL; }
3422 free (*sfit_vol); *sfit_vol = NULL;
3423 }
3424
3425 if (*snfit_vol != NULL)
3426 {
3427 for (it = 0; it < ts_length; it++)
3428 if ((*snfit_vol)[it] != NULL)
3429 { free ((*snfit_vol)[it]); (*snfit_vol)[it] = NULL; }
3430 free (*snfit_vol); *snfit_vol = NULL;
3431 }
3432
3433 }
3434
3435
3436 /*---------------------------------------------------------------------------*/
3437
main(int argc,char ** argv)3438 int main
3439 (
3440 int argc, /* number of input arguments */
3441 char ** argv /* array of input arguments */
3442 )
3443
3444 {
3445 /*----- declare input option variables -----*/
3446 NL_options option_data; /* bucket dataset options */
3447 int nabs; /* use absolute constraints for noise parameters */
3448 int nrand; /* number of random vectors to generate */
3449 int nbest; /* number of random vectors to keep */
3450 float rms_min; /* minimum rms error to reject reduced model */
3451 float fdisp; /* minimum f-statistic for display */
3452 int progress; /* show progress report every n voxels */
3453
3454 /*----- declare time series variables -----*/
3455 THD_3dim_dataset * dset_time = NULL; /* input 3d+time data set */
3456 int ts_length; /* length of time series data */
3457 int ignore; /* delete this number of points from time series */
3458 float ** x_array = NULL; /* independent variable matrix */
3459 float * ts_array = NULL; /* input time series array */
3460 int nxyz; /* number of voxels in image */
3461 int iv; /* voxel counter */
3462
3463 /*----- declare reduced (noise) model variables -----*/
3464 char * nname = NULL; /* noise model name */
3465 vfp nmodel; /* pointer to noise model */
3466 int r; /* number of parameters in the noise model */
3467 char ** npname = NULL; /* noise parameter labels */
3468 float * par_rdcd = NULL; /* estimated parameters for the reduced model */
3469 float sse_rdcd; /* error sum of squares for the reduced model */
3470 float * min_nconstr = NULL; /* min parameter constraints for noise model */
3471 float * max_nconstr = NULL; /* max parameter constraints for noise model */
3472
3473 /*----- declare full (signal+noise) model variables -----*/
3474 char * sname = NULL; /* signal model name */
3475 vfp smodel; /* pointer to signal model */
3476 int p; /* number of parameters in the signal model */
3477 char ** spname = NULL; /* signal parameter labels */
3478 float * par_full = NULL; /* estimated parameters for the full model */
3479 float sse_full; /* error sum of squares for the full model */
3480 float * tpar_full = NULL; /* t-statistic of parameters in full model */
3481 float freg; /* f-statistic for the full regression model */
3482 float rmsreg; /* rms error for the full regression model */
3483 float rsqr; /* R^2 (coef. of multiple determination) */
3484 float smax; /* signed maximum of signal */
3485 float tmax; /* epoch of signed maximum of signal */
3486 float pmax; /* percentage change due to signal */
3487 float area; /* area between signal and baseline */
3488 float parea; /* percent area between signal and baseline */
3489 float * min_sconstr = NULL; /* min parameter constraints for signal model */
3490 float * max_sconstr = NULL; /* max parameter constraints for signal model */
3491
3492 /*----- declare output volume data -----*/
3493 float * rmsreg_vol = NULL; /* rms error for the full regression model */
3494 float * freg_vol = NULL; /* f-statistic for the full regression model */
3495 float * rsqr_vol = NULL; /* R^2 volume data */
3496 float * smax_vol = NULL; /* signed max. of signal volume data */
3497 float * tmax_vol = NULL; /* epoch of signed max. volume data */
3498 float * pmax_vol = NULL; /* max. percentage change due to signal */
3499 float * area_vol = NULL; /* area between signal and baseline */
3500 float * parea_vol = NULL; /* percent area between signal and baseline */
3501 float ** ncoef_vol = NULL; /* noise model parameters volume data */
3502 float ** scoef_vol = NULL; /* signal model parameters volume data */
3503 float ** tncoef_vol = NULL; /* noise model t-statistics volume data */
3504 float ** tscoef_vol = NULL; /* signal model t-statistics volume data */
3505 float ** sfit_vol = NULL; /* voxelwise 3d+time fitted signal model */
3506 float ** snfit_vol = NULL; /* voxelwise 3d+time fitted signal+noise */
3507
3508 /*----- declare file name variables -----*/
3509 char * input_filename = NULL; /* file name of input 3d+time dataset */
3510 char * tfilename = NULL; /* file name of time points */
3511 char * freg_filename = NULL; /* file name for regression f-statistics */
3512 char * frsqr_filename= NULL; /* file name for R^2 statistics */
3513 char * fsmax_filename = NULL; /* file name for signal signed maximum */
3514 char * ftmax_filename = NULL; /* file name for time of signed maximum */
3515 char * fpmax_filename = NULL; /* file name for max. percentage change */
3516 char * farea_filename = NULL; /* file name for area under the signal */
3517 char * fparea_filename = NULL; /* file name for % area under the signal */
3518 char ** fncoef_filename = NULL; /* file name for noise model parameters */
3519 char ** fscoef_filename = NULL; /* file name for signal model parameters */
3520 char ** tncoef_filename = NULL; /* file name for noise model t-statistics */
3521 char ** tscoef_filename = NULL; /* file name for signal model t-statistics */
3522 char * sfit_filename = NULL; /* file name for fitted signal model */
3523 char * snfit_filename = NULL; /* file name for fitted signal+noise model */
3524
3525 char * label; /* report results for one voxel */
3526 int novar; /* flag for insufficient variation in the data */
3527
3528 int ixyz_bot , ixyz_top ;
3529 int voxel_count_index = 0 ; /* for g_voxel_count output 26 Oct 2006 */
3530
3531 /*----- start the elapsed time counter -----*/
3532 (void) COX_clock_time() ;
3533
3534 /*----- Identify software -----*/
3535 #if 0
3536 printf ("\n\n");
3537 printf ("Program: %s \n", PROGRAM_NAME);
3538 printf ("Author: %s \n", PROGRAM_AUTHOR);
3539 printf ("Initial Release: %s \n", PROGRAM_INITIAL);
3540 printf ("Latest Revision: %s \n", PROGRAM_LATEST);
3541 printf ("\n");
3542 #endif
3543
3544 /*-- 20 Apr 2001: addto the arglist, if user wants to [RWCox] --*/
3545 PRINT_VERSION("3dNLfim") ; AUTHOR(PROGRAM_AUTHOR);
3546 mainENTRY("3dNLfim main") ; machdep() ;
3547 THD_check_AFNI_version("3dNLfim") ;
3548
3549 {
3550 int new_argc ; char ** new_argv ;
3551 addto_args( argc , argv , &new_argc , &new_argv ) ;
3552 if( new_argv != NULL ){ argc = new_argc ; argv = new_argv ; }
3553 }
3554
3555
3556 /*----- program initialization -----*/
3557 initialize_program (argc, argv, &ignore,
3558 &nname, &sname, &nmodel, &smodel,
3559 &r, &p, &npname, &spname,
3560 &min_nconstr, &max_nconstr, &min_sconstr, &max_sconstr,
3561 &nabs, &nrand, &nbest, &rms_min, &fdisp,&progress,
3562 &input_filename, &tfilename,
3563 &freg_filename, &frsqr_filename,
3564 &fsmax_filename, &ftmax_filename, &fpmax_filename,
3565 &farea_filename, &fparea_filename, &fncoef_filename,
3566 &fscoef_filename, &tncoef_filename, &tscoef_filename,
3567 &sfit_filename, &snfit_filename,
3568 &dset_time, &nxyz, &ts_length, &x_array, &ts_array,
3569 &par_rdcd, &par_full, &tpar_full,
3570 &rmsreg_vol, &freg_vol, &rsqr_vol,
3571 &smax_vol, &tmax_vol, &pmax_vol, &area_vol, &parea_vol,
3572 &ncoef_vol, &scoef_vol, &tncoef_vol, &tscoef_vol,
3573 &sfit_vol, &snfit_vol, &option_data);
3574
3575 #if 0
3576 #ifdef SAVE_RAN
3577 RAN_setup( nmodel , smodel , r , p , nabs ,
3578 min_nconstr, max_nconstr,
3579 min_sconstr, max_sconstr,
3580 ts_length, x_array, nrand ) ;
3581 #endif
3582 #endif
3583
3584 INFO_message(
3585 "At each voxel, will use %d best of %d random parameter sets",
3586 nbest , nrand ) ;
3587
3588 ixyz_bot = 0 ; ixyz_top = nxyz ; /* RWCox */
3589
3590 #ifdef PROC_MAX
3591 if( proc_numjob > 1 ){ /*---- set up multiple processes ----*/
3592 int vv , nvox=nxyz , nper , pp , nv ;
3593 pid_t newpid ;
3594
3595 /* count number of voxels to compute with into nvox */
3596
3597 if( mask_vol != NULL ){
3598 for( vv=nvox=0 ; vv < nxyz ; vv++ )
3599 if( mask_vol[vv] != 0 ) nvox++ ;
3600 }
3601
3602 if( nvox < proc_numjob ){ /* too few voxels for multiple jobs? */
3603
3604 proc_numjob = 1 ;
3605
3606 } else { /* prepare jobs */
3607
3608 /* split voxels between jobs evenly */
3609
3610 #if 0
3611 if( g_voxel_count ){
3612 g_voxel_count = 0 ;
3613 WARNING_message("-voxel_count disabled by -jobs") ;
3614 }
3615 #endif
3616
3617 nper = nvox / proc_numjob ; /* # voxels per job */
3618 if( mask_vol == NULL ){
3619 proc_vox_bot[0] = 0 ;
3620 for( pp=0 ; pp < proc_numjob ; pp++ ){
3621 proc_vox_top[pp] = proc_vox_bot[pp] + nper ;
3622 if( pp < proc_numjob-1 ) proc_vox_bot[pp+1] = proc_vox_top[pp] ;
3623 }
3624 proc_vox_top[proc_numjob-1] = nxyz ;
3625 } else {
3626 proc_vox_bot[0] = 0 ;
3627 for( pp=0 ; pp < proc_numjob ; pp++ ){
3628 for( nv=0,vv=proc_vox_bot[pp] ; /* count ahead until */
3629 nv < nper && vv < nxyz ; vv++ ){ /* find nper voxels */
3630 if( mask_vol[vv] != 0 ) nv++ ; /* inside the mask */
3631 }
3632 proc_vox_top[pp] = vv ;
3633 if( pp < proc_numjob-1 ) proc_vox_bot[pp+1] = proc_vox_top[pp] ;
3634 }
3635 proc_vox_top[proc_numjob-1] = nxyz ;
3636 }
3637
3638 /* make sure dataset is in memory before forks */
3639
3640 DSET_load(dset_time) ; /* so dataset will be common */
3641
3642 /* start processes */
3643
3644 fprintf(stderr,"++ Voxels in dataset: %d\n",nxyz) ;
3645 if( nvox < nxyz )
3646 fprintf(stderr,"++ Voxels in mask: %d\n",nvox) ;
3647 fprintf(stderr,"++ Voxels per job: %d\n",nper) ;
3648
3649 for( pp=1 ; pp < proc_numjob ; pp++ ){
3650 ixyz_bot = proc_vox_bot[pp] ; /* these 3 variables */
3651 ixyz_top = proc_vox_top[pp] ; /* are for the process */
3652 proc_ind = pp ; /* we're about to fork */
3653 newpid = fork() ;
3654 if( newpid == -1 ){
3655 fprintf(stderr,"** Can't fork job #%d! Error exit!\n",pp);
3656 exit(1) ;
3657 }
3658 if( newpid == 0 ) break ; /* I'm the child */
3659 proc_pid[pp] = newpid ; /* I'm the parent */
3660 iochan_sleep(10) ;
3661 }
3662 if( pp == proc_numjob ){ /* only in the parent */
3663 ixyz_bot = proc_vox_bot[0] ; /* set the 3 control */
3664 ixyz_top = proc_vox_top[0] ; /* variables needed */
3665 proc_ind = 0 ; /* below */
3666 }
3667 fprintf(stderr,"++ Job #%d: processing voxels %d to %d; elapsed time=%.3f\n",
3668 proc_ind,ixyz_bot,ixyz_top-1,COX_clock_time()) ;
3669 }
3670 }
3671 #endif /* PROC_MAX */
3672
3673 if( proc_numjob == 1 )
3674 fprintf(stderr,"++ Calculations starting; elapsed time=%.3f\n",
3675 COX_clock_time()) ;
3676
3677
3678 /*----- loop over voxels in the data set -----*/
3679 for (iv = ixyz_bot; iv < ixyz_top; iv++)
3680 {
3681 /*----- check for mask -----*/
3682 if (mask_vol != NULL)
3683 if (mask_vol[iv] == 0) continue;
3684
3685 /*----- display progress for user (1-based) -----*/
3686 if (g_voxel_count && proc_ind == 0 )
3687 {
3688 /* only print every 100th 26 Oct 2006 [rickr] */
3689 if( voxel_count_index % 100 == 0 )
3690 fprintf(stderr,"\r++ voxel count: %8d (of %d)", iv+1, ixyz_top);
3691 voxel_count_index++ ;
3692 }
3693
3694 AFNI_store_dset_index( iv, -1); /* set voxel index 03 Nov 2006 */
3695 /*----- read the time series for voxel iv -----*/
3696 read_ts_array (dset_time, iv, ts_length, ignore, ts_array);
3697
3698
3699 /*----- calculate the reduced (noise) model -----*/
3700 STATUS("call calc_reduced_model") ;
3701 calc_reduced_model (ts_length, r, x_array, ts_array,
3702 par_rdcd, &sse_rdcd);
3703
3704
3705 /*----- calculate the full (signal+noise) model -----*/
3706 STATUS("call calc_full_model") ;
3707 calc_full_model (nmodel, smodel, r, p,
3708 min_nconstr, max_nconstr, min_sconstr, max_sconstr,
3709 ts_length, x_array, ts_array, par_rdcd, sse_rdcd, nabs,
3710 nrand, nbest, rms_min, par_full, &sse_full, &novar);
3711
3712
3713 /*----- calculate statistics for the full model -----*/
3714 STATUS("call analyze_results") ;
3715 analyze_results (nmodel, smodel, r, p, novar,
3716 min_nconstr, max_nconstr, min_sconstr, max_sconstr,
3717 ts_length, x_array,
3718 par_rdcd, sse_rdcd, par_full, sse_full,
3719 &rmsreg, &freg, &rsqr, &smax, &tmax, &pmax,
3720 &area, &parea, tpar_full);
3721
3722
3723 /*----- report results for this voxel -----*/
3724 if ((freg >= fdisp && proc_ind == 0 )
3725 && (iv % progress == 0))
3726 {
3727 report_results (nname, sname, r, p, npname, spname, ts_length,
3728 par_rdcd, sse_rdcd, par_full, sse_full, tpar_full,
3729 rmsreg, freg, rsqr, smax, tmax, pmax,
3730 area, parea, &label);
3731 printf ("\n\nVoxel #%d\n", iv);
3732 printf ("%s \n", label);
3733 }
3734
3735 /*----- save results for this voxel into volume data -----*/
3736 save_results (iv, nmodel, smodel, r, p, novar, ts_length, x_array,
3737 par_full, tpar_full, rmsreg, freg, rsqr, smax,
3738 tmax, pmax, area, parea, rmsreg_vol,
3739 freg_vol, rsqr_vol, smax_vol,
3740 tmax_vol, pmax_vol, area_vol, parea_vol,ncoef_vol,
3741 scoef_vol, tncoef_vol, tscoef_vol, sfit_vol, snfit_vol);
3742 }
3743 if(g_voxel_count) fputc('\n',stderr);
3744
3745
3746 /*-- if this is a child process, we're done.
3747 if this is the parent process, wait for the children --*/
3748
3749 #ifdef PROC_MAX
3750 if( proc_numjob > 1 ){
3751 if( proc_ind > 0 ){ /* death of child */
3752 fprintf(stderr,"++ Job #%d finished; elapsed time=%.3f\n",
3753 proc_ind,COX_clock_time()) ;
3754 _exit(0) ;
3755
3756 } else { /* parent waits for children */
3757 int pp ;
3758 fprintf(stderr,"++ Job #0 waiting for children to finish; elapsed time=%.3f\n",COX_clock_time()) ;
3759 for( pp=1 ; pp < proc_numjob ; pp++ )
3760 waitpid( proc_pid[pp] , NULL , 0 ) ;
3761 fprintf(stderr,"++ Job #0 now finishing up; elapsed time=%.3f\n",
3762 COX_clock_time()) ;
3763 }
3764
3765 /* when get to here, only parent process is left alive,
3766 and all the results are in the shared memory segment arrays */
3767 }
3768 #endif
3769 if( proc_numjob == 1 )
3770 fprintf(stderr,"++ Calculations finished; elapsed time=%.3f\n",
3771 COX_clock_time()) ;
3772
3773
3774 /*----- delete input dataset -----*/
3775 THD_delete_3dim_dataset( dset_time , False ) ; dset_time = NULL ;
3776
3777
3778 /*----- write requested output files -----*/
3779 output_results (r, p, min_nconstr, max_nconstr, min_sconstr, max_sconstr,
3780 nxyz, ts_length, rmsreg_vol, freg_vol, rsqr_vol,
3781 smax_vol, tmax_vol, pmax_vol, area_vol, parea_vol,
3782 ncoef_vol, scoef_vol, tncoef_vol, tscoef_vol,
3783 sfit_vol, snfit_vol,
3784 input_filename, freg_filename, frsqr_filename,
3785 fsmax_filename, ftmax_filename,
3786 fpmax_filename, farea_filename, fparea_filename,
3787 fncoef_filename, fscoef_filename,
3788 tncoef_filename, tscoef_filename,
3789 sfit_filename, snfit_filename, &option_data);
3790
3791
3792 /*----- end of program -----*/
3793 terminate_program (r, p, ts_length, &x_array, &ts_array,
3794 &nname, &npname, &par_rdcd, &min_nconstr, &max_nconstr,
3795 &sname, &spname, &par_full, &tpar_full,
3796 &min_sconstr, &max_sconstr,
3797 &rmsreg_vol, &freg_vol, &rsqr_vol,
3798 &smax_vol, &tmax_vol, &pmax_vol, &area_vol, &parea_vol,
3799 &ncoef_vol, &scoef_vol, &tncoef_vol, &tscoef_vol,
3800 &sfit_vol, &snfit_vol, &input_filename,
3801 &freg_filename, &frsqr_filename,
3802 &fsmax_filename, &ftmax_filename,
3803 &fpmax_filename, &farea_filename, &fparea_filename,
3804 &fncoef_filename, &fscoef_filename,
3805 &tncoef_filename, &tscoef_filename,
3806 &sfit_filename, &snfit_filename);
3807
3808 if( nwin_pow > 0 && (nwin_sim > 0 || nwin_stp > 0) )
3809 INFO_message(
3810 "# best via POWELL=%d; via SIMPLEX=%d; SIMPLEX then POWELL=%d",
3811 nwin_pow , nwin_sim , nwin_stp ) ;
3812
3813 INFO_message("Program finished; elapsed time=%.3f\n",COX_clock_time()) ;
3814 exit (0);
3815 }
3816