1 /*****************************************************************************
2 Major portions of this software are copyrighted by the Medical College
3 of Wisconsin, 1994-2000, and are released under the Gnu General Public
4 License, Version 2. See the file README.Copyright for details.
5 ******************************************************************************/
6
7 #include "mrilib.h"
8
9 /*------- Adapted from plug_stats.c ---------*/
10
11 #define METH_MEAN 0
12 #define METH_SLOPE 1
13 #define METH_SIGMA 2
14 #define METH_CVAR 3
15
16 #define METH_MEDIAN 4 /* 14 Feb 2001 */
17 #define METH_MAD 5
18
19 #define METH_MAX 6
20 #define METH_MIN 7
21
22 #define METH_DW 8 /* KRH 3 Dec 2002 */
23
24 #define METH_SIGMA_NOD 9 /* KRH 27 Dec 2002 */
25 #define METH_CVAR_NOD 10 /* KRH 27 Dec 2002 */
26
27 #define METH_AUTOCORR 11 /* KRH 16 Jun 2004 */
28 #define METH_AUTOREGP 12 /* KRH 16 Jun 2004 */
29
30 #define METH_ABSMAX 13 /* KRH 15 Feb 2005 */
31
32 #define METH_ARGMAX 14 /* KRH 4 Aug 2005 */
33 #define METH_ARGMIN 15 /* KRH 4 Aug 2005 */
34 #define METH_ARGABSMAX 16 /* KRH 4 Aug 2005 */
35
36 #define METH_SUM 17 /* RWCox 24 Apr 2006 */
37 #define METH_DURATION 18 /* DRG 15 Jun 2007 */
38 #define METH_CENTROID 19
39 #define METH_CENTDURATION 20
40
41 #define METH_ABSSUM 21
42
43 #define METH_NZMEAN 22 /* DRG 03 Oct 2007 */
44 #define METH_ONSET 23 /* DRG 08 Jan 2008 */
45 #define METH_OFFSET 24
46
47 #define METH_ACCUMULATE 25 /* RCR 04 Mar 2008 */
48
49 #define METH_SUM_SQUARES 26 /* ZSS 17 Dec 2008 */
50
51 #define METH_BMV 27 /* RWC 16 Oct 2009 */
52
53 #define METH_ARGMIN1 28 /* ZSS Blizzard 2010 */
54 #define METH_ARGMAX1 29 /* ZSS Blizzard 2010 */
55 #define METH_ARGABSMAX1 30 /* ZSS Blizzard 2010 */
56
57 #define METH_CENTROMEAN 31 /* RWC 01 Nov 2010 */
58
59 #define METH_CVARINV 32 /* RWC 09 Aug 2011 */
60 #define METH_CVARINVNOD 33
61
62 #define METH_ZCOUNT 34
63 #define METH_NZMEDIAN 35 /* RCR 27 Jun 2012 */
64 #define METH_SIGNED_ABSMAX 36 /* RCR 31 Aug 2012 */
65 #define METH_L2_NORM 37 /* RCR 07 Jan 2013 */
66 #define METH_NZCOUNT 38
67
68 #define METH_NZSTDEV 39 /* RWC 29 Jul 2015 */
69
70 #define METH_PERCENTILE 40 /* RWC 05 May 2016 */
71 static int perc_val = -666;
72
73 #define METH_FIRSTVALUE 41 /* returns the 1st value - to avoid exiting on invalid 1-input-methods */
74 #define METH_TSNR 42 /* JKR 10 April 2017 */
75
76 #define METH_MEAN_SSD 43 /* 19 Mar 2018 */
77 #define METH_MEAN_SSDQ 44
78 #define METH_MEDIAN_ASD 45
79
80 #define METH_SKEWNESS 46 /* PDL 17 July 2020 */
81 #define METH_KURTOSIS 47 /* PDL 17 July 2020 */
82
83 #define MAX_NUM_OF_METHS 48
84
85 /* allow single inputs for some methods (test as we care to add) */
86 #define NUM_1_INPUT_METHODS 12
87 static int valid_1_input_methods[NUM_1_INPUT_METHODS]
88 = { METH_MEAN, METH_MAX, METH_MIN, METH_SUM,
89 METH_ARGMAX, METH_ARGMIN, METH_ARGABSMAX,METH_SUM,
90 METH_ABSSUM, METH_NZMEAN, METH_SUM_SQUARES, METH_FIRSTVALUE };
91
92
93 static int meth[MAX_NUM_OF_METHS] = {METH_MEAN};
94 static int nmeths = 0;
95 static char prefix[THD_MAX_PREFIX] = "stat" ;
96 static int datum = MRI_float ;
97 static int nscale = 0;
98 static float basepercent = 0.5; /* 50% assumed for duration unless user specified */
99
100 static int do_tdiff = 0 ; /* 25 May 2011 */
101
102 #if 1
103 # include "Tstat.h"
104 #else
105 static char *meth_names[] = {
106 "Mean" , "Slope" , "Std Dev" , "Coef of Var" ,
107 "Median" , "Med Abs Dev" , "Max" , "Min" ,
108 "Durbin-Watson" , "Std Dev(NOD)" , "Coef Var(NOD)" , "AutoCorr" ,
109 "AutoReg" , "Absolute Max" , "ArgMax" , "ArgMin" ,
110 "ArgAbsMax" , "Sum" , "Duration" , "Centroid" ,
111 "CentDuration" , "Absolute Sum" , "Non-zero Mean" , "Onset" ,
112 "Offset" , "Accumulate" , "SS" , "BiwtMidV" ,
113 "ArgMin+1" , "ArgMax+1" , "ArgAbsMax+1" , "CentroMean" ,
114 "CVarInv" , "CvarInv (NOD)", "ZeroCount" , "NZ Median" ,
115 "Signed Absmax" , "L2 Norm" , "NonZero Count" , "NZ Stdev" ,
116 "Percentile %d" , "FirstValue" , "TSNR" , "MSSD" ,
117 "MSSDsqrt" , "MASDx"
118 };
119 #endif
120
121 static void STATS_tsfunc( double tzero , double tdelta ,
122 int npts , float *ts , double ts_mean ,
123 double ts_slope , void *ud , int nbriks, float *val ) ;
124
125 static void autocorr( int npts, float ints[], int numVals, float outcoeff[] ) ;
126
127 static int Calc_duration(float *ts, int npts, float vmax, int max_index,
128 int *onset, int *offset);
129 static float Calc_centroid(float *ts, int npts);
130
usage_3dTstat(int detail)131 void usage_3dTstat(int detail)
132 {
133
134 printf(
135 "Usage: 3dTstat [options] dataset\n"
136 "Computes one or more voxel-wise statistics for a 3D+time dataset\n"
137 "and stores them in a bucket dataset. If no statistic option is\n"
138 "given, computes just the mean of each voxel time series.\n"
139 "Multiple statistics options may be given, and will result in\n"
140 "a multi-volume dataset.\n"
141 "\n"
142 "Statistics Options (note where detrending does/does not occur):\n"
143 "\n"
144 " -sum = compute sum of input voxels\n"
145 " -abssum = compute absolute sum of input voxels\n"
146 " -sos = compute sum of squares\n"
147 " -l2norm = compute L2 norm (sqrt(sum squares))\n"
148 " -mean = compute mean of input voxels\n"
149 " -slope = compute the slope of input voxels vs. time\n"
150 "\n"
151 " -stdev = compute standard deviation of input voxels\n"
152 " NB: input is detrended by first removing mean+slope\n"
153 " -stdevNOD = like '-stdev', but no initial detrending\n"
154 " -cvar = compute coefficient of variation of input:\n"
155 " voxels = stdev/fabs(mean)\n"
156 " NB: in stdev calc, input is detrended by removing mean+slope\n"
157 " -cvarNOD = like '-cvar', but no initial detrending in stdev calc\n"
158 " -cvarinv = 1.0/cvar = 'signal to noise ratio' [for Vinai]\n"
159 " NB: in stdev calc, input is detrended by removing mean+slope\n"
160 " -cvarinvNOD = like '-cvarinv', but no detrending in stdev calc\n"
161 "\n"
162 " -tsnr = compute temporal signal to noise ratio\n"
163 " fabs(mean)/stdev NOT DETRENDED (same as -cvarinvNOD)\n"
164 " -MAD = compute MAD (median absolute deviation) of\n"
165 " input voxels = median(|voxel-median(voxel)|)\n"
166 " [N.B.: the trend is NOT removed for this]\n"
167 " -DW = compute Durbin-Watson Statistic of input voxels\n"
168 " [N.B.: the trend IS removed for this]\n"
169 " -median = compute median of input voxels [undetrended]\n"
170 " -nzmedian = compute median of non-zero input voxels [undetrended]\n"
171 " -nzstdev = standard deviation of non-zero input voxel [undetrended]\n"
172 " -bmv = compute biweight midvariance of input voxels [undetrended]\n"
173 " [actually is 0.989*sqrt(biweight midvariance), to make]\n"
174 " [the value comparable to the standard deviation output]\n"
175 " -MSSD = Von Neumann's Mean of Successive Squared Differences\n"
176 " = average of sum of squares of first time difference\n"
177 " -MSSDsqrt = Sqrt(MSSD)\n"
178 " -MASDx = Median of absolute values of first time differences\n"
179 " times 1.4826 (to scale it like standard deviation)\n"
180 " = a robust alternative to MSSDsqrt\n"
181 " -min = compute minimum of input voxels [undetrended]\n"
182 " -max = compute maximum of input voxels [undetrended]\n"
183 " -absmax = compute absolute maximum of input voxels [undetrended]\n"
184 " -signed_absmax = (signed) value with absolute maximum [undetrended]\n"
185 " -percentile P = the P-th percentile point (0=min, 50=median, 100=max)\n"
186 " of the data in each voxel time series.\n"
187 " [this option can only be used once!]\n"
188 " -argmin = index of minimum of input voxels [undetrended]\n"
189 " -argmin1 = index + 1 of minimum of input voxels [undetrended]\n"
190 " -argmax = index of maximum of input voxels [undetrended]\n"
191 " -argmax1 = index + 1 of maximum of input voxels [undetrended]\n"
192 " -argabsmax = index of absolute maximum of input voxels [undetrended]\n"
193 " -argabsmax1= index +1 of absolute maximum of input voxels [undetrended]\n"
194 " -duration = compute number of points around max above a threshold\n"
195 " Use basepercent option to set limits\n"
196 " -onset = beginning of duration around max where value\n"
197 " exceeds basepercent\n"
198 " -offset = end of duration around max where value\n"
199 " exceeds basepercent\n"
200 " -centroid = compute centroid of data time curves\n"
201 " (sum(i*f(i)) / sum(f(i)))\n"
202 " -centduration = compute duration using centroid's index as center\n"
203 " -nzmean = compute mean of non-zero voxels\n"
204 " -zcount = count number of zero values at each voxel\n"
205 " -nzcount = count number of non zero values at each voxel\n"
206 "\n"
207 " -autocorr n = compute autocorrelation function and return\n"
208 " first n coefficients\n"
209 " -autoreg n = compute autoregression coefficients and return\n"
210 " first n coefficients\n"
211 " [N.B.: -autocorr 0 and/or -autoreg 0 will return number\n"
212 " coefficients equal to the length of the input data]\n"
213 "\n"
214 " -accumulate = accumulate time series values (partial sums)\n"
215 " val[i] = sum old_val[t] over t = 0..i\n"
216 " (output length = input length)\n"
217 "\n"
218 " -centromean = compute mean of middle 50%% of voxel values [undetrended]\n"
219 "\n"
220 " -skewness = measure of asymmetry in distribution - based on Pearson's\n"
221 " moment, coefficient of skewness.\n"
222 " -kurtosis = measure of the 'tailedness' of the probability distribution\n"
223 " - the fourth standardized moment. Never negative.\n"
224 " -firstvalue = first value in dataset - typically just placeholder\n\n"
225 " ** If no statistic option is given, then '-mean' is assumed **\n"
226 "\n"
227 "Other Options:\n"
228 " -tdiff = Means to take the first difference of each time\n"
229 " series before further processing.\n"
230 " -prefix p = Use string 'p' for the prefix of the\n"
231 " output dataset [DEFAULT = 'stat']\n"
232 " -datum d = use data type 'd' for the type of storage\n"
233 " of the output, where 'd' is one of\n"
234 " 'byte', 'short', or 'float' [DEFAULT=float]\n"
235 " -nscale = Do not scale output values when datum is byte or short.\n"
236 " Scaling is done by default.\n"
237 "\n"
238 " -basepercent nn = Percentage of maximum for duration calculation\n"
239 "\n"
240 " -mask mset Means to use the dataset 'mset' as a mask:\n"
241 " Only voxels with nonzero values in 'mset'\n"
242 " will be printed from 'dataset'. Note\n"
243 " that the mask dataset and the input dataset\n"
244 " must have the same number of voxels.\n"
245 "\n"
246 " -mrange a b Means to further restrict the voxels from\n"
247 " 'mset' so that only those mask values\n"
248 " between 'a' and 'b' (inclusive) will\n"
249 " be used. If this option is not given,\n"
250 " all nonzero values from 'mset' are used.\n"
251 " Note that if a voxel is zero in 'mset', then\n"
252 " it won't be included, even if a < 0 < b.\n"
253 "\n"
254 " -cmask 'opts' Means to execute the options enclosed in single\n"
255 " quotes as a 3dcalc-like program, and produce\n"
256 " produce a mask from the resulting 3D brick.\n"
257 " Examples:\n"
258 " -cmask '-a fred+orig[7] -b zork+orig[3] -expr step(a-b)'\n"
259 " produces a mask that is nonzero only where\n"
260 " the 7th sub-brick of fred+orig is larger than\n"
261 " the 3rd sub-brick of zork+orig.\n"
262 " -cmask '-a fred+orig -expr 1-bool(k-7)'\n"
263 " produces a mask that is nonzero only in the\n"
264 " 7th slice (k=7); combined with -mask, you\n"
265 " could use this to extract just selected voxels\n"
266 " from particular slice(s).\n"
267 " Notes: * You can use both -mask and -cmask in the same\n"
268 " run - in this case, only voxels present in\n"
269 " both masks will be dumped.\n"
270 " * Only single sub-brick calculations can be\n"
271 " used in the 3dcalc-like calculations -\n"
272 " if you input a multi-brick dataset here,\n"
273 " without using a sub-brick index, then only\n"
274 " its 0th sub-brick will be used.\n"
275 " * Do not use quotes inside the 'opts' string!\n"
276 "\n"
277 "\n"
278 "If you want statistics on a detrended dataset and the option\n"
279 "doesn't allow that, you can use program 3dDetrend first.\n"
280 "\n"
281 "The output is a bucket dataset. The input dataset may\n"
282 "use a sub-brick selection list, as in program 3dcalc.\n\n"
283 "*** If you are trying to compute the mean or std.dev. of multiple\n"
284 "datasets (not across time), use 3dMean or 3dmerge instead.\n"
285
286 ) ;
287
288 printf("\n"
289 "----------------- Processing 1D files with 3dTstat -----------------\n"
290 "To analyze a 1D file and get statistics on each of its columns,\n"
291 "you can do something like this:\n"
292 " 3dTstat -stdev -bmv -prefix stdout: file.1D\\'\n"
293 "where the \\' means to transpose the file on input, since 1D files\n"
294 "read into 3dXXX programs are interpreted as having the time direction\n"
295 "along the rows rather than down the columns. In this example, the\n"
296 "output is written to the screen, which could be captured with '>'\n"
297 "redirection. Note that if you don't give the '-prefix stdout:'\n"
298 "option, then the output will be written into a NIML-formatted 1D\n"
299 "dataset, which you might find slightly confusing (but still usable).\n"
300 ) ;
301
302 PRINT_COMPILE_DATE ;
303
304 return;
305 }
306
main(int argc,char * argv[])307 int main( int argc , char *argv[] )
308 {
309 THD_3dim_dataset *old_dset , *new_dset ; /* input and output datasets */
310 THD_3dim_dataset *mask_dset=NULL ;
311 float mask_bot=666.0 , mask_top=-666.0 ;
312 byte *cmask=NULL ; int ncmask=0 ;
313 byte *mmm = NULL ;
314 int mcount=0, verb=0;
315 int nopt, nbriks, ii ;
316 int addBriks = 0 ; /* n-1 sub-bricks out */
317 int fullBriks = 0 ; /* n sub-bricks out */
318 int tsout = 0 ; /* flag to output a time series (not a stat bucket) */
319 int numMultBriks,methIndex,brikIndex;
320
321 /*----- Help the pitiful user? -----*/
322
323
324 /* bureaucracy */
325 mainENTRY("3dTstat main"); machdep(); AFNI_logger("3dTstat",argc,argv);
326 PRINT_VERSION("3dTstat"); AUTHOR("KR Hammett & RW Cox");
327
328 /*--- scan command line for options ---*/
329
330 if (argc == 1) { usage_3dTstat(1); exit(0); } /* Bob's help shortcut */
331 nopt = 1 ;
332 nbriks = 0 ;
333 nmeths = 0 ;
334 verb = 0;
335 while( nopt < argc && argv[nopt][0] == '-' ){
336 if (strcmp(argv[nopt], "-h") == 0 || strcmp(argv[nopt], "-help") == 0) {
337 usage_3dTstat(strlen(argv[nopt]) > 3 ? 2:1);
338 exit(0);
339 }
340
341 if( strcmp(argv[nopt],"-verb") == 0 ){
342 verb++ ; nopt++ ; continue ;
343 }
344
345 /*-- methods --*/
346
347 if( strcasecmp(argv[nopt],"-methnum") == 0 ){ /* 20 Mar 2018 */
348 if( ++nopt == argc ) /* for plugin use */
349 ERROR_exit("no arg after '%s'",argv[nopt-1]);
350 meth[nmeths++] = (int)strtod(argv[nopt],NULL) ;
351 nbriks++ ;
352 nopt++ ; continue ;
353 }
354
355 if( strcasecmp(argv[nopt],"-centromean") == 0 ){ /* 01 Nov 2010 */
356 meth[nmeths++] = METH_CENTROMEAN ;
357 nbriks++ ;
358 nopt++ ; continue ;
359 }
360
361 if( strcasecmp(argv[nopt],"-bmv") == 0 ){
362 meth[nmeths++] = METH_BMV ;
363 nbriks++ ;
364 nopt++ ; continue ;
365 }
366
367 if( strcasecmp(argv[nopt],"-median") == 0 ){
368 meth[nmeths++] = METH_MEDIAN ;
369 nbriks++ ;
370 nopt++ ; continue ;
371 }
372
373 if( strcasecmp(argv[nopt],"-nzmedian") == 0 ){
374 meth[nmeths++] = METH_NZMEDIAN ;
375 nbriks++ ;
376 nopt++ ; continue ;
377 }
378
379 if( strcasecmp(argv[nopt],"-nzstdev") == 0 ){
380 meth[nmeths++] = METH_NZSTDEV ;
381 nbriks++ ;
382 nopt++ ; continue ;
383 }
384
385 if( strcasecmp(argv[nopt],"-DW") == 0 ){
386 meth[nmeths++] = METH_DW ;
387 nbriks++ ;
388 nopt++ ; continue ;
389 }
390
391 if( strcasecmp(argv[nopt],"-zcount") == 0 ){
392 meth[nmeths++] = METH_ZCOUNT ;
393 nbriks++ ;
394 nopt++ ; continue ;
395 }
396
397 if( strcasecmp(argv[nopt],"-nzcount") == 0 ){
398 meth[nmeths++] = METH_NZCOUNT ;
399 nbriks++ ;
400 nopt++ ; continue ;
401 }
402
403 if( strcasecmp(argv[nopt],"-MAD") == 0 ){
404 meth[nmeths++] = METH_MAD ;
405 nbriks++ ;
406 nopt++ ; continue ;
407 }
408
409 if( strcasecmp(argv[nopt],"-mean") == 0 ){
410 meth[nmeths++] = METH_MEAN ;
411 nbriks++ ;
412 nopt++ ; continue ;
413 }
414
415 if( strcasecmp(argv[nopt],"-percentile") == 0 ){ /* 05 May 2016 */
416 if( ++nopt >= argc )
417 ERROR_exit("need argument after option '%s'",argv[nopt-1]) ;
418 if( perc_val >= 0 )
419 ERROR_exit("you can't use option '%s' twice!",argv[nopt-1]) ;
420 perc_val = (int)strtod(argv[nopt],NULL) ;
421 if( perc_val < 0 || perc_val > 100 )
422 ERROR_exit("illegal number '%d' after option '%s'",
423 perc_val , argv[nopt-1]) ;
424 if( perc_val == 0 )
425 WARNING_message("option '%s 0' is the same as '-min'",argv[nopt-1]) ;
426 if( perc_val == 100 )
427 WARNING_message("option '%s 100' is the same as '-max'",argv[nopt-1]) ;
428 meth[nmeths++] = METH_PERCENTILE ;
429 nbriks++ ; nopt++ ; continue ;
430 }
431
432 if( strcasecmp(argv[nopt],"-MSSD") == 0 ){ /* 19 Mar 2018 */
433 meth[nmeths++] = METH_MEAN_SSD ;
434 nbriks++ ;
435 nopt++ ; continue ;
436 }
437 if( strcasecmp(argv[nopt],"-MSSDQ") == 0 ||
438 strcasecmp(argv[nopt],"-MSSDsqrt") == 0 ){ /* 20 Mar 2018 */
439 meth[nmeths++] = METH_MEAN_SSDQ ;
440 nbriks++ ;
441 nopt++ ; continue ;
442 }
443 if( strcasecmp(argv[nopt],"-MASDx") == 0 ){ /* 19 Mar 2018 */
444 meth[nmeths++] = METH_MEDIAN_ASD ;
445 nbriks++ ;
446 nopt++ ; continue ;
447 }
448 if( strcasecmp(argv[nopt],"-skewness") == 0 ){ /* 19 Mar 2018 */
449 meth[nmeths++] = METH_SKEWNESS ;
450 nbriks++ ;
451 nopt++ ; continue ;
452 }
453 if( strcasecmp(argv[nopt],"-kurtosis") == 0 ){ /* 19 Mar 2018 */
454 meth[nmeths++] = METH_KURTOSIS ;
455 nbriks++ ;
456 nopt++ ; continue ;
457 }
458
459 if( strcasecmp(argv[nopt],"-sum") == 0 ){
460 meth[nmeths++] = METH_SUM ;
461 nbriks++ ;
462 nopt++ ; continue ;
463 }
464 if( strcasecmp(argv[nopt],"-sos") == 0 ){
465 meth[nmeths++] = METH_SUM_SQUARES ;
466 nbriks++ ;
467 nopt++ ; continue ;
468 }
469
470 if( strcasecmp(argv[nopt],"-abssum") == 0 ){
471 meth[nmeths++] = METH_ABSSUM ;
472 nbriks++ ;
473 nopt++ ; continue ;
474 }
475
476 if( strcasecmp(argv[nopt],"-slope") == 0 ){
477 meth[nmeths++] = METH_SLOPE ;
478 nbriks++ ;
479 nopt++ ; continue ;
480 }
481
482 if( strcasecmp(argv[nopt],"-stdev") == 0 ||
483 strcasecmp(argv[nopt],"-sigma") == 0 ){
484
485 meth[nmeths++] = METH_SIGMA ;
486 nbriks++ ;
487 nopt++ ; continue ;
488 }
489
490 if( strcasecmp(argv[nopt],"-cvar") == 0 ){
491 meth[nmeths++] = METH_CVAR ;
492 nbriks++ ;
493 nopt++ ; continue ;
494 }
495
496 if( strcasecmp(argv[nopt],"-cvarinv") == 0 ){
497 meth[nmeths++] = METH_CVARINV ;
498 nbriks++ ;
499 nopt++ ; continue ;
500 }
501
502 if( strcasecmp(argv[nopt],"-stdevNOD") == 0 ||
503 strcasecmp(argv[nopt],"-sigmaNOD") == 0 ){ /* 07 Dec 2001 */
504
505 meth[nmeths++] = METH_SIGMA_NOD ;
506 nbriks++ ;
507 nopt++ ; continue ;
508 }
509
510 if( strcasecmp(argv[nopt],"-cvarNOD") == 0 ){ /* 07 Dec 2001 */
511 meth[nmeths++] = METH_CVAR_NOD ;
512 nbriks++ ;
513 nopt++ ; continue ;
514 }
515
516 if( strcasecmp(argv[nopt],"-cvarinvNOD") == 0 ){
517 meth[nmeths++] = METH_CVARINVNOD ;
518 nbriks++ ;
519 nopt++ ; continue ;
520 }
521
522 if( strcasecmp(argv[nopt],"-tsnr") == 0 ){ /* 10 April 2017 */
523 meth[nmeths++] = METH_TSNR ;
524 nbriks++ ;
525 nopt++ ; continue ;
526 }
527
528 if( strcasecmp(argv[nopt],"-min") == 0 ){
529 meth[nmeths++] = METH_MIN ;
530 nbriks++ ;
531 nopt++ ; continue ;
532 }
533
534 if( strcasecmp(argv[nopt],"-max") == 0 ){
535 meth[nmeths++] = METH_MAX ;
536 nbriks++ ;
537 nopt++ ; continue ;
538 }
539
540 if( strcasecmp(argv[nopt],"-absmax") == 0 ){
541 meth[nmeths++] = METH_ABSMAX ;
542 nbriks++ ;
543 nopt++ ; continue ;
544 }
545
546 if( strcasecmp(argv[nopt],"-signed_absmax") == 0 ){
547 meth[nmeths++] = METH_SIGNED_ABSMAX ;
548 nbriks++ ;
549 nopt++ ; continue ;
550 }
551
552 if( strcasecmp(argv[nopt],"-argmin") == 0 ){
553 meth[nmeths++] = METH_ARGMIN ;
554 nbriks++ ;
555 nopt++ ; continue ;
556 }
557
558 if( strcasecmp(argv[nopt],"-argmin1") == 0 ){
559 meth[nmeths++] = METH_ARGMIN1 ;
560 nbriks++ ;
561 nopt++ ; continue ;
562 }
563
564 if( strcasecmp(argv[nopt],"-argmax") == 0 ){
565 meth[nmeths++] = METH_ARGMAX ;
566 nbriks++ ;
567 nopt++ ; continue ;
568 }
569
570 if( strcasecmp(argv[nopt],"-argmax1") == 0 ){
571 meth[nmeths++] = METH_ARGMAX1 ;
572 nbriks++ ;
573 nopt++ ; continue ;
574 }
575
576 if( strcasecmp(argv[nopt],"-argabsmax") == 0 ){
577 meth[nmeths++] = METH_ARGABSMAX ;
578 nbriks++ ;
579 nopt++ ; continue ;
580 }
581
582 if( strcasecmp(argv[nopt],"-argabsmax1") == 0 ){
583 meth[nmeths++] = METH_ARGABSMAX1;
584 nbriks++ ;
585 nopt++ ; continue ;
586 }
587
588 if( strcasecmp(argv[nopt],"-duration") == 0 ){
589 meth[nmeths++] = METH_DURATION ;
590 nbriks++ ;
591 nopt++ ; continue ;
592 }
593
594 if( strcasecmp(argv[nopt],"-onset") == 0 ){
595 meth[nmeths++] = METH_ONSET ;
596 nbriks++ ;
597 nopt++ ; continue ;
598 }
599
600 if( strcasecmp(argv[nopt],"-offset") == 0 ){
601 meth[nmeths++] = METH_OFFSET ;
602 nbriks++ ;
603 nopt++ ; continue ;
604 }
605
606 if( strcasecmp(argv[nopt],"-centroid") == 0 ){
607 meth[nmeths++] = METH_CENTROID ;
608 nbriks++ ;
609 nopt++ ; continue ;
610 }
611 if( strcasecmp(argv[nopt],"-centduration") == 0 ){
612 meth[nmeths++] = METH_CENTDURATION ;
613 nbriks++ ;
614 nopt++ ; continue ;
615 }
616
617 if( strcasecmp(argv[nopt],"-nzmean") == 0 ){
618 meth[nmeths++] = METH_NZMEAN ;
619 nbriks++ ;
620 nopt++ ; continue ;
621 }
622
623 if( strncmp(argv[nopt],"-mask",5) == 0 ){
624 if( mask_dset != NULL )
625 ERROR_exit("Cannot have two -mask options!\n") ;
626 if( nopt+1 >= argc )
627 ERROR_exit("-mask option requires a following argument!\n");
628 mask_dset = THD_open_dataset( argv[++nopt] ) ;
629 if( mask_dset == NULL )
630 ERROR_exit("Cannot open mask dataset!\n") ;
631 if( DSET_BRICK_TYPE(mask_dset,0) == MRI_complex )
632 ERROR_exit("Cannot deal with complex-valued mask dataset!\n");
633 nopt++ ; continue ;
634 }
635
636 if( strncmp(argv[nopt],"-mrange",5) == 0 ){
637 if( nopt+2 >= argc )
638 ERROR_exit("-mrange option requires 2 following arguments!\n");
639 mask_bot = strtod( argv[++nopt] , NULL ) ;
640 mask_top = strtod( argv[++nopt] , NULL ) ;
641 if( mask_top < mask_top )
642 ERROR_exit("-mrange inputs are illegal!\n") ;
643 nopt++ ; continue ;
644 }
645
646 if( strcmp(argv[nopt],"-cmask") == 0 ){ /* 16 Mar 2000 */
647 if( nopt+1 >= argc )
648 ERROR_exit("-cmask option requires a following argument!\n");
649 cmask = EDT_calcmask( argv[++nopt] , &ncmask, 0 ) ;
650 if( cmask == NULL ) ERROR_exit("Can't compute -cmask!\n");
651 nopt++ ; continue ;
652 }
653
654 if( strcasecmp(argv[nopt],"-autocorr") == 0 ){
655 meth[nmeths++] = METH_AUTOCORR ;
656 if( ++nopt >= argc ) ERROR_exit("-autocorr needs an argument!\n");
657 meth[nmeths++] = atoi(argv[nopt++]);
658 if (meth[nmeths - 1] == 0) {
659 addBriks++;
660 } else {
661 nbriks+=meth[nmeths - 1] ;
662 }
663 continue ;
664 }
665
666 if( strcasecmp(argv[nopt],"-autoreg") == 0 ){
667 meth[nmeths++] = METH_AUTOREGP ;
668 if( ++nopt >= argc ) ERROR_exit("-autoreg needs an argument!\n");
669 meth[nmeths++] = atoi(argv[nopt++]);
670 if (meth[nmeths - 1] == 0) {
671 addBriks++;
672 } else {
673 nbriks+=meth[nmeths - 1] ;
674 }
675 continue ;
676 }
677
678 if( strcasecmp(argv[nopt],"-accumulate") == 0 ){ /* 4 Mar 2008 [rickr] */
679 meth[nmeths++] = METH_ACCUMULATE ;
680 meth[nmeths++] = -1; /* flag to add N (not N-1) output bricks */
681 fullBriks++;
682 tsout = 1; /* flag to output a timeseries */
683 nopt++ ; continue ;
684 }
685
686 if( strcasecmp(argv[nopt],"-l2norm") == 0 ){ /* 07 Jan 2013 [rickr] */
687 meth[nmeths++] = METH_L2_NORM ;
688 nbriks++ ;
689 nopt++ ; continue ;
690 }
691
692 /*-- prefix --*/
693
694 if( strcasecmp(argv[nopt],"-prefix") == 0 ){
695 if( ++nopt >= argc ) ERROR_exit("-prefix needs an argument!\n");
696 MCW_strncpy(prefix,argv[nopt],THD_MAX_PREFIX) ;
697 if( !THD_filename_ok(prefix) )
698 ERROR_exit("%s is not a valid prefix!\n",prefix);
699 nopt++ ; continue ;
700 }
701
702 /*-- tdiff --*/
703
704 if( strcasecmp(argv[nopt],"-tdiff") == 0 ){ /* 25 May 2011 */
705 do_tdiff = 1 ; nopt++ ; continue ;
706 }
707
708 /*-- nscale --*/
709
710 if( strcasecmp(argv[nopt],"-nscale") == 0 ){ /* 25 May 2011 */
711 nscale = 1 ; nopt++ ; continue ;
712 }
713
714 /*-- datum --*/
715
716 if( strcasecmp(argv[nopt],"-datum") == 0 ){
717 if( ++nopt >= argc ) ERROR_exit("-datum needs an argument!\n");
718 if( strcasecmp(argv[nopt],"short") == 0 ){
719 datum = MRI_short ;
720 } else if( strcasecmp(argv[nopt],"float") == 0 ){
721 datum = MRI_float ;
722 } else if( strcasecmp(argv[nopt],"byte") == 0 ){
723 datum = MRI_byte ;
724 } else {
725 ERROR_exit("-datum of type '%s' is not supported!\n",
726 argv[nopt] ) ;
727 }
728 nopt++ ; continue ;
729 }
730
731 /* base percentage for duration calcs */
732 if (strcasecmp (argv[nopt], "-basepercent") == 0) {
733 if( ++nopt >= argc ) ERROR_exit("-basepercent needs an argument!\n");
734 basepercent = strtod(argv[nopt], NULL);
735 if(basepercent>1) basepercent /= 100.0; /* assume integer percent if >1*/
736 nopt++; continue;
737 }
738
739 /*-- Quien sabe'? --*/
740
741 ERROR_message("Unknown option: %s\n",argv[nopt]) ;
742 suggest_best_prog_option(argv[0], argv[nopt]);
743 exit(1);
744 }
745
746 if (argc < 2) {
747 ERROR_message("Too few options, use -help for details");
748 exit(1);
749 }
750
751 /*--- If no options selected, default to single stat MEAN -- KRH ---*/
752
753 if (nmeths == 0) nmeths = 1;
754 if (nbriks == 0 && addBriks == 0 && fullBriks == 0) nbriks = 1;
755
756 /*----- read input dataset -----*/
757
758 if( nopt >= argc ) ERROR_exit(" No input dataset!?") ;
759
760 old_dset = THD_open_dataset( argv[nopt] ) ;
761 if( !ISVALID_DSET(old_dset) )
762 ERROR_exit("Can't open dataset %s\n",argv[nopt]);
763
764 nopt++ ;
765 if( nopt < argc )
766 WARNING_message("Trailing datasets on command line ignored: %s ...",argv[nopt]) ;
767
768 if( DSET_NVALS(old_dset) == 1 ){
769 WARNING_message("Input dataset has 1 sub-brick ==> -tdiff is turned off") ;
770 do_tdiff = 0 ;
771 }
772
773 /* no input volumes is bad, 1 volume applies to only certain methods */
774 /* 2 Nov 2010 [rickr] */
775 if( DSET_NVALS(old_dset) == 0 ) {
776 ERROR_exit("Time series is of length 0?\n") ;
777 }
778 else if( DSET_NVALS(old_dset) == 1 || (do_tdiff && DSET_NVALS(old_dset)==2) ) {
779 int methOK;
780 /* see if each method is valid for nvals == 1 */
781 for( methIndex = 0; methIndex < nmeths; methIndex++ ) {
782 methOK = 0;
783 for( ii = 0; ii < NUM_1_INPUT_METHODS; ii++ ) {
784 if( meth[methIndex] == valid_1_input_methods[ii] ) {
785 methOK = 1;
786 break;
787 }
788 else
789 meth[methIndex] = METH_FIRSTVALUE;
790 }
791 if( ! methOK )
792 WARNING_message("Using dataset with only %d values per voxel and giving back just the first value!" ,
793 DSET_NVALS(old_dset) ) ;
794 }
795 /* tell the library function that this case is okay */
796 g_thd_maker_allow_1brick = 1;
797 }
798
799 if( DSET_NUM_TIMES(old_dset) < 2 ){
800 WARNING_message("Input dataset is not 3D+time; assuming TR=1.0") ;
801 EDIT_dset_items( old_dset ,
802 ADN_ntt , DSET_NVALS(old_dset) ,
803 ADN_ttorg , 0.0 ,
804 ADN_ttdel , 1.0 ,
805 ADN_tunits , UNITS_SEC_TYPE ,
806 NULL ) ;
807 }
808
809 /* If one or more of the -autocorr/-autoreg options was called with */
810 /* an argument of 0, then I'll now add extra BRIKs for the N-1 data */
811 /* output points for each. */
812 nbriks += ((DSET_NVALS(old_dset)-1) * addBriks);
813 nbriks += ((DSET_NVALS(old_dset) ) * fullBriks);
814
815 /* ------------- Mask business -----------------*/
816 if( mask_dset == NULL ){
817 mmm = NULL ;
818 if( verb )
819 INFO_message("%d voxels in the entire dataset (no mask)\n",
820 DSET_NVOX(old_dset)) ;
821 } else {
822 if( DSET_NVOX(mask_dset) != DSET_NVOX(old_dset) )
823 ERROR_exit("Input and mask datasets are not same dimensions!\n");
824 mmm = THD_makemask( mask_dset , 0 , mask_bot, mask_top ) ;
825 mcount = THD_countmask( DSET_NVOX(old_dset) , mmm ) ;
826 if( mcount <= 0 ) ERROR_exit("No voxels in the mask!\n") ;
827 if( verb ) INFO_message("%d voxels in the mask\n",mcount) ;
828 DSET_delete(mask_dset) ;
829 }
830
831 if( cmask != NULL ){
832 if( ncmask != DSET_NVOX(old_dset) )
833 ERROR_exit("Input and cmask datasets are not same dimensions!\n");
834 if( mmm != NULL ){
835 for( ii=0 ; ii < DSET_NVOX(old_dset) ; ii++ )
836 mmm[ii] = (mmm[ii] && cmask[ii]) ;
837 free(cmask) ;
838 mcount = THD_countmask( DSET_NVOX(old_dset) , mmm ) ;
839 if( mcount <= 0 ) ERROR_exit("No voxels in the mask+cmask!\n") ;
840 if( verb ) INFO_message("%d voxels in the mask+cmask\n",mcount) ;
841 } else {
842 mmm = cmask ;
843 mcount = THD_countmask( DSET_NVOX(old_dset) , mmm ) ;
844 if( mcount <= 0 ) ERROR_exit("No voxels in the cmask!\n") ;
845 if( verb ) INFO_message("%d voxels in the cmask\n",mcount) ;
846 }
847 }
848
849 /*------------- ready to compute new dataset -----------*/
850
851 new_dset = MAKER_4D_to_typed_fbuc(
852 old_dset , /* input dataset */
853 prefix , /* output prefix */
854 datum , /* output datum */
855 0 , /* ignore count */
856 0 , /* can't detrend in maker function KRH 12/02*/
857 nbriks , /* number of briks */
858 STATS_tsfunc , /* timeseries processor */
859 NULL, /* data for tsfunc */
860 mmm,
861 nscale
862 ) ;
863
864 if( new_dset != NULL ){
865 tross_Copy_History( old_dset , new_dset ) ;
866 tross_Make_History( "3dTstat" , argc,argv , new_dset ) ;
867 for (methIndex = 0,brikIndex = 0; methIndex < nmeths;
868 methIndex++, brikIndex++) {
869 if ((meth[methIndex] == METH_AUTOCORR) ||
870 (meth[methIndex] == METH_ACCUMULATE) ||
871 (meth[methIndex] == METH_AUTOREGP)) {
872 numMultBriks = meth[methIndex+1];
873
874 /* note: this looks like it should be NV-1 4 Mar 2008 [rickr] */
875 if (numMultBriks == 0) numMultBriks = DSET_NVALS(old_dset)-1;
876 /* new flag for NVALS [rickr] */
877 if (numMultBriks == -1) numMultBriks = DSET_NVALS(old_dset);
878
879 for (ii = 1; ii <= numMultBriks; ii++) {
880 char tmpstr[25];
881 if (meth[methIndex] == METH_AUTOREGP) {
882 sprintf(tmpstr,"%s[%d](%d)",meth_names[meth[methIndex]],
883 numMultBriks,ii);
884 } else {
885 sprintf(tmpstr,"%s(%d)",meth_names[meth[methIndex]],ii);
886 }
887 EDIT_BRICK_LABEL(new_dset, (brikIndex + ii - 1), tmpstr) ;
888 }
889 methIndex++;
890 brikIndex += numMultBriks - 1;
891 } else if( meth[methIndex] == METH_PERCENTILE ){ /* 05 May 2016 */
892 char plabel[32] ;
893 sprintf(plabel,meth_names[METH_PERCENTILE],perc_val) ;
894 EDIT_BRICK_LABEL(new_dset, brikIndex, plabel ) ;
895 } else {
896 EDIT_BRICK_LABEL(new_dset, brikIndex, meth_names[meth[methIndex]]) ;
897 }
898 }
899
900 if( tsout ) /* then change output to a time series */
901 EDIT_dset_items( new_dset ,
902 ADN_ntt , brikIndex ,
903 ADN_ttorg , DSET_TIMEORIGIN(old_dset) ,
904 ADN_ttdel , DSET_TIMESTEP(old_dset) ,
905 ADN_tunits , DSET_TIMEUNITS(old_dset) ,
906 NULL ) ;
907
908 DSET_write( new_dset ) ;
909 WROTE_DSET( new_dset ) ;
910 } else {
911 ERROR_exit("Unable to compute output dataset!\n") ;
912 }
913
914 exit(0) ;
915 }
916
917 /**********************************************************************
918 Function that does the real work
919 ***********************************************************************/
920
STATS_tsfunc(double tzero,double tdelta,int npts,float * ts,double ts_mean,double ts_slope,void * ud,int nbriks,float * val)921 static void STATS_tsfunc( double tzero, double tdelta ,
922 int npts, float *ts ,
923 double ts_mean, double ts_slope,
924 void *ud, int nbriks, float *val )
925 {
926 static int ncall ;
927 int meth_index, ii , out_index, nzpts, onset, offset, duration;
928 float *ts_det, *ts_dif=NULL ;
929
930 /** is this a "notification"? **/
931
932 if( val == NULL ){
933
934 if( npts > 0 ){ /* the "start notification" */
935
936 ncall = 0 ; /* number of calls */
937
938 } else { /* the "end notification" */
939
940 /* nothing to do here */
941
942 }
943 return ;
944 }
945
946 /* RWC: first difference here [25 May 2011] */
947
948 if( do_tdiff ){
949 float tsm , tss ;
950 ts_dif = (float*)calloc(npts, sizeof(float)) ;
951 for( ii=1 ; ii < npts ; ii++ ) ts_dif[ii-1] = ts[ii]-ts[ii-1] ;
952 get_linear_trend( npts-1 , ts_dif , &tsm , &tss ) ;
953 ts_mean = (double)tsm ; ts_slope = (double)tss ;
954 npts-- ; ts = ts_dif ;
955 }
956
957 /* KRH make copy and detrend it right here.
958 Use ts_mean and ts_slope for detrend-ization. */
959
960 ts_det = (float*)calloc(npts, sizeof(float));
961 memcpy( ts_det, ts, npts * sizeof(float));
962
963 for( ii = 0; ii < npts; ii++)
964 ts_det[ii] -=
965 (ts_mean - (ts_slope * (npts - 1) * tdelta/2) + ts_slope * tdelta * ii) ;
966
967 /** OK, actually do some work **/
968
969 /* This main loop now uses meth_index AND out_index as loop variables, */
970 /* mainly to avoid rewriting code that worked. */
971
972 /* meth_index is an index into the static method array, which contains the */
973 /* list of methods to be executed (and will also contain an integer */
974 /* parameter specifying the number of return values if -autocorr n and/or */
975 /* -autoreg p have been requested). */
976 /* out_index is an index into the output array. */
977
978 for(meth_index=out_index=0 ; meth_index < nmeths; meth_index++,out_index++){
979 switch( meth[meth_index] ){
980
981 default:
982 case METH_FIRSTVALUE: val[out_index] = ts[0]; break ; /* placeholder drg 06/13/2016 */
983 case METH_MEAN: val[out_index] = ts_mean ; break ;
984
985 case METH_SUM: val[out_index] = ts_mean * npts; break; /* 24 Apr 2006 */
986
987 case METH_ABSSUM:{ /* 18 Jun 2006 */
988 register int ii ;
989 register float sum ;
990 sum = 0.0;
991 for(ii=0; ii< npts; ii++) sum += fabs(ts[ii]);
992 val[out_index] = sum;
993 }
994 break;
995
996 case METH_PERCENTILE:{ /* 05 May 2016 */
997 float* ts_copy; int ii , fi ;
998 ts_copy = (float*)calloc(npts, sizeof(float));
999 memcpy( ts_copy, ts, npts * sizeof(float));
1000 qsort_float(npts,ts_copy) ;
1001 fi = (perc_val * npts)*0.01f ;
1002 ii = (int)fi ; fi = fi - ii ;
1003 if( ii >= npts-1 ){ /* max */
1004 val[out_index] = ts_copy[npts-1] ;
1005 } else if( ii < 0 ){ /* should never happen */
1006 val[out_index] = ts_copy[0] ;
1007 } else { /* the usual suspects */
1008 val[out_index] = (1.0f-fi)*ts_copy[ii] + fi*ts_copy[ii+1] ;
1009 }
1010 free(ts_copy) ;
1011 }
1012 break ;
1013
1014 case METH_L2_NORM: /* 07 Jan 2013 [rickr] */
1015 case METH_SUM_SQUARES:{ /* 18 Dec 2008 */
1016 register int ii ;
1017 register float sum ;
1018 sum = 0.0;
1019 for(ii=0; ii< npts; ii++) sum += ts[ii]*ts[ii];
1020
1021 /* check multiple methods here */
1022 if ( meth[meth_index] == METH_SUM_SQUARES )
1023 val[out_index] = sum;
1024 else if ( meth[meth_index] == METH_L2_NORM ) {
1025 /* theory and practice do not always seem to agree, so... */
1026 if( sum >= 0.0 ) val[out_index] = sqrt(sum);
1027 else val[out_index] = 0.0;
1028 }
1029 }
1030 break;
1031
1032 case METH_SLOPE: val[out_index] = ts_slope ; break ;
1033
1034 case METH_CVAR_NOD: /* methods that depend on the mean and stdev */
1035 case METH_SIGMA_NOD:
1036 case METH_CVARINVNOD:
1037 case METH_CVAR:
1038 case METH_CVARINV:
1039 case METH_SIGMA:{
1040 register int ii ;
1041 register double sum ;
1042
1043 int mm = meth[meth_index] ;
1044 int nod = (mm == METH_CVAR_NOD) || /* no detrend flag */
1045 (mm == METH_SIGMA_NOD) ||
1046 (mm == METH_CVARINVNOD) ;
1047
1048 sum = 0.0 ;
1049 if( !nod ){ /* not no detrend ==> use detrended data */
1050 for( ii=0 ; ii < npts ; ii++ ) sum += ts_det[ii] * ts_det[ii] ;
1051 } else { /* use data as God gave it to us */
1052 for( ii=0 ; ii < npts ; ii++ ) sum += (ts[ii]-ts_mean)
1053 *(ts[ii]-ts_mean) ;
1054 }
1055
1056 sum = sqrt( sum/(npts-1.0) ) ; /* stdev */
1057
1058 if( mm == METH_SIGMA || mm == METH_SIGMA_NOD )
1059 val[out_index] = sum ;
1060 else if( mm == METH_CVAR || mm == METH_CVAR_NOD )
1061 val[out_index] = (ts_mean != 0.0) ? sum/fabs(ts_mean) : 0.0 ;
1062 else
1063 val[out_index] = (sum != 0.0) ? fabs(ts_mean)/sum : 0.0 ;
1064 }
1065 break ;
1066
1067 case METH_TSNR:{
1068 register int ii ;
1069 register double std = 0.0;
1070 register double sum = 0.0;
1071
1072 /* no detrending */
1073 for( ii=0 ; ii < npts ; ii++ ) sum += (ts[ii]-ts_mean)
1074 *(ts[ii]-ts_mean) ;
1075 std = sqrt( sum/(npts-1.0) ) ; /* stdev */
1076 val[out_index] = (std != 0.0) ? fabs(ts_mean)/std : 0.0 ;
1077 }
1078 break ;
1079
1080 case METH_MEAN_SSDQ:
1081 case METH_MEAN_SSD:{ /* 19 Mar 2018 */
1082 val[out_index] = cs_mean_square_sd( npts , ts ) ;
1083 if( meth[meth_index] == METH_MEAN_SSDQ )
1084 val[out_index] = sqrtf( MAX(val[out_index],0.0f) ) ;
1085 }
1086 break ;
1087
1088 case METH_MEDIAN_ASD:{ /* 19 Mar 2018 */
1089 val[out_index] = cs_median_abs_sd( npts , ts , NULL ) * 1.4826f ;
1090 }
1091 break ;
1092
1093 /* 14 Feb 2000: these 2 new methods disturb the array ts[] */
1094 /* 18 Dec 2002: these 2 methods no longer disturb the array ts[] */
1095
1096 case METH_MEDIAN:{
1097 float* ts_copy;
1098 ts_copy = (float*)calloc(npts, sizeof(float));
1099 memcpy( ts_copy, ts, npts * sizeof(float));
1100 val[out_index] = qmed_float( npts , ts_copy ) ;
1101 free(ts_copy);
1102 }
1103 break ;
1104
1105 case METH_NZMEDIAN:{ /* 27 Jun 2012 [rickr] */
1106 float* ts_copy;
1107 int lind, lnzcount;
1108 ts_copy = (float*)calloc(npts, sizeof(float));
1109 /* replace memcpy with non-zero copy */
1110 lnzcount=0;
1111 for (lind=0; lind < npts; lind++)
1112 if( ts[lind] ) ts_copy[lnzcount++] = ts[lind];
1113 /* and get the result from the possibly shortened array */
1114 if( lnzcount > 0 ) val[out_index] = qmed_float( lnzcount , ts_copy ) ;
1115 else val[out_index] = 0.0 ;
1116 free(ts_copy);
1117 }
1118 break ;
1119
1120 case METH_NZSTDEV:{ /* 29 Jul 2015 [RWC] */
1121 float* ts_copy;
1122 int lind, lnzcount;
1123 ts_copy = (float*)calloc(npts, sizeof(float));
1124 /* replace memcpy with non-zero copy */
1125 lnzcount=0;
1126 for (lind=0; lind < npts; lind++)
1127 if( ts[lind] ) ts_copy[lnzcount++] = ts[lind];
1128 /* and get the result from the possibly shortened array */
1129 if( lnzcount > 1 ) meansigma_float( lnzcount,ts_copy, NULL,val+out_index ) ;
1130 else val[out_index] = 0.0 ;
1131 free(ts_copy);
1132 }
1133 break ;
1134
1135 case METH_CENTROMEAN:{ /* 01 Nov 2010 */
1136 float* ts_copy;
1137 ts_copy = (float*)calloc(npts, sizeof(float));
1138 memcpy( ts_copy, ts, npts * sizeof(float));
1139 val[out_index] = centromean_float( npts , ts_copy ) ;
1140 free(ts_copy);
1141 }
1142 break ;
1143
1144 case METH_MAD:{
1145 float* ts_copy;
1146 register int ii ;
1147 register float vm ;
1148 ts_copy = (float*)calloc(npts, sizeof(float));
1149 memcpy( ts_copy, ts, npts * sizeof(float));
1150 vm = qmed_float( npts , ts_copy ) ;
1151 for( ii=0 ; ii < npts ; ii++ ) ts_copy[ii] = fabs(ts_copy[ii]-vm) ;
1152 val[out_index] = qmed_float( npts , ts_copy ) ;
1153 free(ts_copy);
1154 }
1155 break ;
1156
1157 case METH_BMV:{ /* 16 Oct 2009 */
1158 float bmv ;
1159 qmedmadbmv_float( npts , ts , NULL,NULL , &bmv ) ;
1160 val[out_index] = bmv ;
1161 }
1162 break ;
1163
1164 case METH_ARGMIN:
1165 case METH_ARGMIN1:
1166 case METH_MIN:{
1167 register int ii,outdex=0 ;
1168 register float vm=ts[0] ;
1169 for( ii=1 ; ii < npts ; ii++ ) if( ts[ii] < vm ) {
1170 vm = ts[ii] ;
1171 outdex = ii ;
1172 }
1173 if (meth[meth_index] == METH_MIN) {
1174 val[out_index] = vm ;
1175 } else if (meth[meth_index] == METH_ARGMIN) {
1176 val[out_index] = outdex ;
1177 } else {
1178 val[out_index] = outdex +1;
1179 }
1180 }
1181 break ;
1182
1183 case METH_ZCOUNT:{ /* 05 Apr 2012 */
1184 int ii , zc ;
1185 for( ii=zc=0 ; ii < npts ; ii++ ) if( ts[ii] == 0.0f ) zc++ ;
1186 val[out_index] = zc ;
1187 }
1188 break ;
1189
1190 case METH_NZCOUNT:{ /* Turkey 2014 */
1191 int ii , zc ;
1192 for( ii=zc=0 ; ii < npts ; ii++ ) if( ts[ii] == 0.0f ) zc++ ;
1193 val[out_index] = npts-zc ;
1194 }
1195 break ;
1196
1197 case METH_DW:{
1198 register int ii ;
1199 register float den=ts[0]*ts[0] ;
1200 register float num=0 ;
1201 for( ii=1 ; ii < npts ; ii++ ) {
1202 num = num + (ts_det[ii] - ts_det[ii-1])
1203 *(ts_det[ii] - ts_det[ii-1]);
1204 den = den + ts_det[ii] * ts_det[ii];
1205 }
1206 if (den == 0) {
1207 val[out_index] = 0 ;
1208 } else {
1209 val[out_index] = num/den ;
1210 }
1211 }
1212 break ;
1213
1214 case METH_ONSET:
1215 case METH_OFFSET:
1216 case METH_DURATION:
1217 case METH_ABSMAX:
1218 case METH_SIGNED_ABSMAX:
1219 case METH_ARGMAX:
1220 case METH_ARGMAX1:
1221 case METH_ARGABSMAX:
1222 case METH_ARGABSMAX1:
1223 case METH_MAX:
1224 case METH_CENTDURATION:
1225 case METH_CENTROID:{
1226 register int ii, outdex=0 ;
1227 register float vm=ts[0];
1228 if ( (meth[meth_index] == METH_ABSMAX) ||
1229 (meth[meth_index] == METH_ARGABSMAX) ||
1230 (meth[meth_index] == METH_ARGABSMAX1) ) {
1231 vm = fabs(vm) ;
1232 for( ii=1 ; ii < npts ; ii++ ) {
1233 if( fabs(ts[ii]) > vm ) {
1234 vm = fabs(ts[ii]) ;
1235 outdex = ii ;
1236 }
1237 }
1238 } else if ( meth[meth_index] == METH_SIGNED_ABSMAX ) {
1239 /* for P Hamilton 31 Aug 2012 [rickr] */
1240 register float avm=fabs(vm);
1241 for( ii=1 ; ii < npts ; ii++ ) {
1242 if( fabs(ts[ii]) > avm ) {
1243 vm = ts[ii] ;
1244 avm = fabs(vm) ;
1245 outdex = ii ;
1246 }
1247 }
1248 } else {
1249 for( ii=1 ; ii < npts ; ii++ ) {
1250 if( ts[ii] > vm ) {
1251 vm = ts[ii] ;
1252 outdex = ii ;
1253 }
1254 }
1255 }
1256 switch( meth[meth_index] ){
1257 default:
1258 case METH_ABSMAX:
1259 case METH_SIGNED_ABSMAX:
1260 case METH_MAX:
1261 val[out_index] = vm ;
1262 break;
1263
1264 case METH_ONSET:
1265 case METH_OFFSET:
1266 case METH_DURATION:
1267 duration = Calc_duration(ts, npts, vm, outdex,
1268 &onset,&offset);
1269 switch(meth[meth_index]) {
1270 case METH_ONSET: val[out_index] = onset; break;
1271 case METH_OFFSET: val[out_index] = offset; break;
1272 case METH_DURATION: val[out_index] = duration; break;
1273 }
1274
1275 break;
1276
1277 case METH_ARGMAX:
1278 case METH_ARGABSMAX:
1279 val[out_index] = outdex ;
1280 break;
1281
1282 case METH_ARGMAX1:
1283 case METH_ARGABSMAX1:
1284 val[out_index] = outdex +1;
1285 break;
1286
1287 case METH_CENTROID:
1288 case METH_CENTDURATION: {
1289 float cm;
1290 cm = Calc_centroid(ts, npts);
1291 if(meth[meth_index]== METH_CENTDURATION)
1292 val[out_index] = Calc_duration(ts, npts, vm, (int) cm,
1293 &onset,&offset);
1294 else
1295 val[out_index] = cm;
1296 }
1297 break;
1298 }
1299 }
1300 break ;
1301
1302 case METH_NZMEAN: {
1303 register int ii ;
1304 register float sum ;
1305 sum = 0.0;
1306 nzpts = 0;
1307 for( ii=0 ; ii < npts ; ii++ ) {
1308 if( ts[ii] != 0.0 ) {
1309 sum += ts[ii] ;
1310 nzpts++;
1311 }
1312 }
1313 if(npts>0)
1314 val[out_index] = sum / nzpts;
1315 else
1316 val[out_index] = 0.0;
1317 }
1318 break;
1319
1320 case METH_AUTOCORR:{
1321 int numVals;
1322 float *ts_corr;
1323 /* for these new methods, the extra, needed integer */
1324 /* parameter is stored in the static array "meth", */
1325 /* in the element right after the indicator for */
1326 /* this method. This parameter indicates the number*/
1327 /* of values to return, or 0 for the same length as */
1328 /* the input array. */
1329 numVals = meth[meth_index + 1];
1330 if (numVals == 0) numVals = npts - 1;
1331 ts_corr = (float*)calloc(numVals,sizeof(float));
1332 /* Call the autocorrelation function, in this file. */
1333 autocorr(npts,ts,numVals,ts_corr);
1334 /* Copy the values into the output array val, which */
1335 /* will be returned to the fbuc MAKER caller to */
1336 /* populate the appropriate BRIKs with the data. */
1337 for( ii = 0; ii < numVals; ii++) {
1338 val[out_index + ii] = ts_corr[ii];
1339 }
1340 /* Although meth_index will be incremented by the */
1341 /* for loop, it needs to be incremented an extra */
1342 /* time to account for the extra parameter we */
1343 /* pulled from the meth array. */
1344 meth_index++;
1345 /* Similarly, though the out_index will be incremented */
1346 /* by the for loop, we have added potentially several */
1347 /* values to the output array, and we need to account */
1348 /* for that here. */
1349 out_index+=(numVals - 1);
1350 free(ts_corr);
1351 }
1352 break ;
1353
1354 case METH_AUTOREGP:{
1355 int numVals,kk,mm;
1356 float *ts_corr, *y, *z;
1357 float alpha, beta, tmp;
1358
1359 /* For these new methods, the extra, needed integer */
1360 /* parameter is stored in the static array "meth", */
1361 /* in the element right after the indicator for */
1362 /* this method. This parameter indicates the number*/
1363 /* of values to return, or 0 for the same length as */
1364 /* the input array. */
1365 numVals = meth[meth_index + 1];
1366 if (numVals == 0) numVals = npts - 1;
1367
1368 /* Allocate space for the autocorrelation coefficients, */
1369 /* result, and a temp array for calculations. */
1370 /* Correlation coeff array is larger, because we must */
1371 /* disregard the r0, which is identically 1. */
1372 /* Might fix this in autocorr function */
1373 ts_corr = (float*)malloc((numVals) * sizeof(float));
1374 y = (float*)malloc(numVals * sizeof(float));
1375 z = (float*)malloc(numVals * sizeof(float));
1376
1377 /* Call autocorr function to generate the autocorrelation */
1378 /* coefficients. */
1379 autocorr(npts,ts,numVals,ts_corr);
1380
1381 /* Durbin algorithm for solving Yule-Walker equations. */
1382 /* The Yule-Walker equations specify the autoregressive */
1383 /* coefficients in terms of the autocorrelation coefficients. */
1384 /* R*Phi = r, where r is vector of correlation coefficients, */
1385 /* R is Toeplitz matrix formed from r, and Phi is a vector of */
1386 /* the autoregression coefficients. */
1387
1388 /* In this implementation, 'y' is 'Phi' above and 'ts_corr' is 'r' */
1389
1390 y[0] = -ts_corr[0];
1391 alpha = -ts_corr[0];
1392 beta = 1;
1393 for (kk = 0 ; kk < (numVals - 1) ; kk++ ) {
1394 beta = (1 - alpha * alpha ) * beta ;
1395 tmp = 0;
1396 for ( mm = 0 ; mm <= kk ; mm++) {
1397 tmp = tmp + ts_corr[kk - mm] * y[mm];
1398 }
1399 alpha = - ( ts_corr[kk+1] + tmp ) / beta ;
1400 for ( mm = 0 ; mm <= kk ; mm++ ) {
1401 z[mm] = y[mm] + alpha * y[kk - mm];
1402 }
1403 for ( mm = 0 ; mm <= kk ; mm++ ) {
1404 y[mm] = z[mm];
1405 }
1406 y[kk+1] = alpha ;
1407 }
1408
1409 /* Copy the values into the output array val, which */
1410 /* will be returned to the fbuc MAKER caller to */
1411 /* populate the appropriate BRIKs with the data. */
1412 for( ii = 0; ii < numVals; ii++) {
1413 val[out_index + ii] = y[ii];
1414 if (!isfinite(y[ii])){
1415 WARNING_message("BAD FLOAT y[%d]=%f; Call#%d\n",ii,y[ii],ncall);
1416 val[out_index + ii] = 0.0f ;
1417 }
1418 }
1419 /* Although meth_index will be incremented by the */
1420 /* for loop, it needs to be incremented an extra */
1421 /* time to account for the extra parameter we */
1422 /* pulled from the meth array. */
1423 meth_index++;
1424 /* Similarly, though the out_index will be incremented */
1425 /* by the for loop, we have added potentially several */
1426 /* values to the output array, and we need to account */
1427 /* for that here. */
1428 out_index+=(numVals - 1);
1429 free(ts_corr);
1430 free(y);
1431 free(z);
1432 }
1433 break ;
1434
1435 case METH_ACCUMULATE:{
1436 register double sum = 0.0 ;
1437 register int ii ;
1438
1439 meth_index++; /* go past dummy zero */
1440
1441 for( ii = 0; ii < npts; ii++) {
1442 sum += ts[ii]; /* keep track as double */
1443 val[out_index + ii] = sum;
1444 }
1445
1446 out_index+=(npts - 1);
1447 }
1448 break ;
1449
1450 /* New methods, added by PDL, 17 July 2020 */
1451
1452 case METH_SKEWNESS:{
1453
1454 val[out_index] = getSkewness(ts, npts);
1455 }
1456 break;
1457
1458 case METH_KURTOSIS:{
1459
1460 val[out_index] = getKurtosis(ts, npts);
1461 }
1462 break;
1463
1464 }
1465 }
1466
1467 free(ts_det); if( ts_dif == ts ) free(ts_dif) ;
1468 ncall++ ; return ;
1469 }
1470
1471
autocorr(int npts,float in_ts[],int numVals,float outcoeff[])1472 static void autocorr( int npts, float in_ts[], int numVals, float outcoeff[] )
1473 {
1474 /* autocorr is the inv fourier xform of the fourier xform abs squared */
1475
1476 int ii,nfft;
1477 double scaler;
1478 complex *cxar = NULL;
1479 float tsmean ;
1480 static int nfft_old = 0 ;
1481
1482 /* Calculate size for FFT, including padding for eliminating overlap */
1483 /* from circular convolution */
1484 nfft = csfft_nextup_even(npts * 2 - 1);
1485
1486 if( nfft != nfft_old ){
1487 nfft_old = nfft ;
1488 INFO_message("FFT length = %d",nfft) ;
1489 }
1490
1491 cxar = (complex *) calloc( sizeof(complex) , nfft ) ;
1492
1493 /* compute mean of input */
1494 for( tsmean=0.0f,ii=0 ; ii < npts ; ii++ ) tsmean += in_ts[ii] ;
1495 tsmean /= npts ;
1496
1497 /* Populate complex array with input (real-only) time series */
1498 for( ii=0 ; ii < npts ; ii++ ){
1499 cxar[ii].r = in_ts[ii] - tsmean; cxar[ii].i = 0.0;
1500 }
1501 /* Zero-pad input outside range of original time series */
1502 for( ii=npts ; ii < nfft ; ii++ ){ cxar[ii].r = cxar[ii].i = 0.0; }
1503
1504 /* Calculate nfft-point FFT of input time series */
1505 /* First argument is "mode." -1 value indicates */
1506 /* negative exponent in fourier sum, which means */
1507 /* this is an FFT and not an inverse FFT. */
1508 /* Output will be complex and symmetrical. */
1509 csfft_cox( -1 , nfft , cxar ) ;
1510
1511 /* Use macro to calculate absolute square of FFT */
1512 for( ii=0 ; ii < nfft ; ii++ ){
1513 cxar[ii].r = CSQR(cxar[ii]) ; cxar[ii].i = 0 ;
1514 }
1515
1516 /* Take inverse FFT of result. First function called */
1517 /* sets a static variable that the output should be */
1518 /* scaled by 1/N. This plus the mode argument of 1 */
1519 /* indicate that this is an inverse FFT. */
1520 csfft_scale_inverse( 1 ) ;
1521 csfft_cox( 1 , nfft, cxar ) ;
1522
1523 /* Copy the requested number of coefficients to the */
1524 /* output array outcoeff. These will be copied into */
1525 /* the output BRIKs or used for further calculations */
1526 /* in the calling function, STATS_tsfunc() above. */
1527 /* The output coefficients are scaled by 1/(M-p) to */
1528 /* provide an unbiased estimate, and also scaled by */
1529 /* the final value of the zeroth coefficient. */
1530 scaler = cxar[0].r/npts;
1531 if( scaler != 0.0 ) scaler = 1.0 / scaler ;
1532 for (ii = 0 ; ii < numVals ; ii++ ) {
1533 outcoeff[ii] = cxar[ii+1].r * ( scaler /(npts - (ii+1)) ) ;
1534 }
1535 free(cxar);
1536 cxar = NULL;
1537 }
1538
1539 /* calculate the duration of a peak in number of subbricks */
1540 /* duration is the number of points at or above the threshold*/
1541 static int
Calc_duration(float ts[],int npts,float vmax,int max_index,int * onset,int * offset)1542 Calc_duration(float ts[], int npts, float vmax, int max_index,
1543 int *onset, int *offset)
1544 {
1545 float minlimit;
1546 int i;
1547 ENTRY("Calc_duration");
1548 /* find beginning - onset - first point before max that falls below min */
1549 minlimit = basepercent * vmax;
1550 i = max_index; /* for centroid option we need to start at max_index*/
1551 *onset = -1;
1552 *offset = npts;
1553 while(i>=0) {
1554 if(ts[i]<minlimit) {
1555 *onset = i;
1556 break;
1557 }
1558 i--;
1559 }
1560
1561 /* find end of peak - offset - first point after max that falls below min */
1562 /* this starting index needs to be the same for centroid or peak*/
1563 i = max_index + 1;
1564 while(i<npts) {
1565 if(ts[i]<minlimit) {
1566 *offset = i;
1567 break;
1568 }
1569 i++;
1570 }
1571 RETURN(*offset - *onset -1);
1572 }
1573
1574 /* calculate the centroid of a time series*/
1575 /* centroid or center of mass is defined as Sum(i*f(i)) / Sum(f(i)) */
1576 static float
Calc_centroid(float ts[],int npts)1577 Calc_centroid(float ts[], int npts)
1578 {
1579 int i;
1580 double sum=0.0, tvsum = 0.0;
1581
1582 ENTRY("Calc_centroid");
1583 for(i=0;i<npts;i++) {
1584 sum += ts[i];
1585 tvsum += i*ts[i];
1586 }
1587
1588 RETURN(tvsum / sum);
1589 }
1590