1 /* ----------------------------------------------------------------------
2 * Do stuff.
3 *
4 *
5 * Author: R Reynolds 18 July 2017
6 */
7
8 static char * g_history[] =
9 {
10 "----------------------------------------------------------------------\n"
11 "history (of 3dTto1D):\n"
12 "\n",
13 "0.0 18 Jul 2017\n"
14 " (Rick Reynolds of the National Institutes of Health, SSCC/DIRP/NIMH)\n"
15 " -initial version\n"
16 "1.0 19 Jul 2017: initial release\n"
17 "1.1 18 Aug 2017: modify help\n"
18 "1.2 01 Feb 2018: added -4095_count/frac/warn\n"
19 };
20
21 static char g_version[] = "3dTto1D version 1.2, 1 February 2018";
22
23 #include "mrilib.h"
24
25
26 /*--------------- method enumeration ---------------*/
27 #define T21_METH_UNDEF 0
28 #define T21_METH_ENORM 1 /* enorm */
29 #define T21_METH_RMS 2 /* rms */
30 #define T21_METH_SRMS 3 /* srms */
31 #define T21_METH_S_SRMS 4 /* shift_srms */
32 #define T21_METH_MDIFF 5 /* mdiff */
33 #define T21_METH_SMDIFF 6 /* smdiff */
34 #define T21_METH_4095_COUNT 7 /* 4095_count */
35 #define T21_METH_4095_FRAC 8 /* 4095_frac */
36 #define T21_METH_4095_WARN 9 /* 4095_warn */
37
38 /*--------------- global options struct ---------------*/
39 typedef struct
40 {
41 THD_3dim_dataset * inset;
42 char * mask_name;
43 char * prefix;
44 int method;
45 int automask;
46 int verb;
47
48 /* put parameters here, too */
49 int nt;
50 float * result; /* of length nt */
51 byte * mask;
52 } options_t;
53
54 options_t g_opts;
55
56 /*--------------- prototypes ---------------*/
57
58 int show_help (void);
59 int check_dims (options_t *);
60 int compute_results (options_t *);
61 int compute_4095 (options_t *, int);
62 int compute_enorm (options_t *, int);
63 int compute_meandiff (options_t *, int);
64 int fill_mask (options_t *);
65 int meth_name_to_index (char *);
66 char * meth_index_to_name (int);
67 int process_opts (options_t *, int, char *[] );
68 int write_results (options_t *);
69
70 /*--------------- main routine ---------------*/
main(int argc,char * argv[])71 int main( int argc, char *argv[] )
72 {
73 options_t * opts = &g_opts;
74 int rv;
75
76 if( argc < 2 ) { show_help(); return 0; }
77
78 mainENTRY("3dTto1D main"); machdep(); AFNI_logger("3dTto1D",argc,argv);
79
80 /* process command line arguments (and read dataset and mask) */
81 rv = process_opts(opts, argc, argv);
82 if( rv ) RETURN(rv < 0); /* only a negative return is considered failure */
83
84 if( check_dims(opts) ) RETURN(1);
85
86 /* evaluation of rv now depends on the method, usually non-zero is bad */
87 rv = compute_results(opts);
88
89 /* for 4095_warn, return now, regardless */
90 if( opts->method == T21_METH_4095_WARN ) RETURN(rv);
91
92 /* otherwise, any non-zero return is a failure */
93 if ( rv ) RETURN(1);
94
95 if( write_results(opts) ) RETURN(1);
96
97 RETURN(0);
98 }
99
100
101 /* just make sure we have sufficient data for computations */
check_dims(options_t * opts)102 int check_dims(options_t * opts)
103 {
104 int nt, nvox, nmask;
105
106 ENTRY("check_dims");
107
108 nt = DSET_NVALS(opts->inset);
109 nvox = DSET_NVOX(opts->inset);
110 if( opts->mask ) nmask = THD_countmask( nvox, opts->mask );
111 else nmask = nvox;
112
113 /* make sure we have something to compute */
114 if( nvox < 1 ) {
115 ERROR_message("input dataset must have at least 1 voxel");
116 RETURN(1);
117 } else if( nmask < 1 ) {
118 ERROR_message("input mask must have at least 1 voxel");
119 RETURN(1);
120 } else if( nt < 2 ) {
121 ERROR_message("input dataset must have at least 2 time points");
122 RETURN(1);
123 }
124
125 RETURN(0);
126 }
127
128
write_results(options_t * opts)129 int write_results(options_t * opts)
130 {
131 FILE * fp;
132
133 int c;
134
135 ENTRY("write_results");
136
137 if( ! opts->result ) {
138 ERROR_message("no results to write!");
139 RETURN(1);
140 }
141
142 if( ! opts->prefix ) fp = stdout;
143 else if ( !strcmp(opts->prefix, "-") ) fp = stdout;
144 else if ( !strcmp(opts->prefix, "stdout") ) fp = stdout;
145 else if ( !strcmp(opts->prefix, "stderr") ) fp = stderr;
146 else {
147 fp = fopen(opts->prefix, "w");
148 if( ! fp ) {
149 ERROR_message("failed to open '%s' for writing", opts->prefix);
150 RETURN(1);
151 }
152 }
153
154 if( opts->verb ) {
155 if ( fp == stdout ) INFO_message("writing to stdout...");
156 else if( fp == stdout ) INFO_message("writing to stderr...");
157 else if( opts->prefix ) INFO_message("writing to '%s'...", opts->prefix);
158 else INFO_message("writing to unknown?!?");
159 }
160
161 /* actually write results */
162 for( c=0; c < opts->nt; c++ )
163 fprintf(fp, "%f\n", opts->result[c]);
164
165 if( fp != stdout && fp != stderr )
166 fclose(fp);
167
168 RETURN(0);
169 }
170
171
compute_results(options_t * opts)172 int compute_results(options_t * opts)
173 {
174 int method = opts->method;
175
176 ENTRY("compute_results");
177
178 /* the first 3 methods are varieties of enorm/dvars */
179 if( method == T21_METH_ENORM ||
180 method == T21_METH_RMS ||
181 method == T21_METH_SRMS ||
182 method == T21_METH_S_SRMS) RETURN(compute_enorm(opts, method));
183
184 if( method == T21_METH_MDIFF ||
185 method == T21_METH_SMDIFF ) RETURN(compute_meandiff(opts, method));
186
187 if( method == T21_METH_4095_COUNT ||
188 method == T21_METH_4095_FRAC ||
189 method == T21_METH_4095_WARN ) RETURN(compute_4095(opts, method));
190
191 ERROR_message("unknown method index %d", method);
192
193 RETURN(1);
194 }
195
196
197 /* this computes the mean abs first (backwards) diff */
compute_meandiff(options_t * opts,int method)198 int compute_meandiff(options_t * opts, int method)
199 {
200 double * dwork, dscale=1.0, gmean, mdiff;
201 float * fdata, fdiff;
202 byte * mask = opts->mask; /* save 6 characters... */
203 int nt, nvox, vind, tind, nmask;
204
205 ENTRY("compute_meandiff");
206
207 nt = DSET_NVALS(opts->inset);
208 nvox = DSET_NVOX(opts->inset);
209
210 /* note how many voxels this is over */
211 if( opts->mask ) nmask = THD_countmask( nvox, opts->mask );
212 else nmask = nvox;
213
214 if( opts->verb )
215 INFO_message("computing %s, nvox = %d, nmask = %d, nt = %d",
216 meth_index_to_name(method), nvox, nmask, nt);
217
218 /* check_dims() has been called, so we have sufficient data */
219
220 dwork = calloc(nt, sizeof(double));
221 fdata = calloc(nt, sizeof(float));
222
223 /* could steal fdata, but that might be unethical (plus, garbage on err) */
224 opts->result = calloc(nt, sizeof(float));
225
226 /* use gmean to get mean across all masked voxels and time */
227 gmean = 0.0;
228 mdiff = 0.0;
229 for( vind=0; vind < nvox; vind++ ) {
230 if( opts->mask && ! opts->mask[vind] ) continue;
231
232 if( THD_extract_array(vind, opts->inset, 0, fdata) ) {
233 ERROR_message("failed to exract data at index %d\n", vind);
234 free(dwork); free(fdata); RETURN(1);
235 }
236
237 /* accumuate squared differences; dwork[0] is already 0 */
238 gmean += fdata[0]; /* and accumlate for mean */
239 for( tind=1; tind < nt; tind++ ) {
240 gmean += fdata[tind]; /* accumlate for mean */
241 fdiff = fabs(fdata[tind]-fdata[tind-1]);
242 mdiff += fdiff; /* accumulate for mean diff */
243 dwork[tind] += fdiff;
244 }
245 }
246
247 /* convert dsum to dmean and then a scalar:
248 * method = 1 : enorm sqrt(ss)
249 * 2 : rms sqrt(ss/nmask) = stdev (biased) of first diffs
250 * 3 : srms sqrt(ss/nmask) / abs(grand mean)
251 */
252 gmean = fabs(gmean/nmask/nt); /* global masked mean */
253 mdiff /= nmask*nt; /* global masked mean first diff */
254
255 /* set the enorm scalar, based on the method */
256 if( method == T21_METH_MDIFF ) {
257 if( opts->verb ) INFO_message("mean diff uses no scaling");
258 dscale = nmask;
259 } else if ( method == T21_METH_SMDIFF ) {
260 if( opts->verb ) INFO_message("scaled mean diff = mdiff/gmean");
261 if( gmean == 0.0 ) {
262 ERROR_message("values have zero mean, failing (use mdiff, instead)");
263 free(dwork); free(fdata); RETURN(1);
264 } else if( gmean < 1.0 ) {
265 WARNING_message("values might be de-meaned, consider mdiff");
266 }
267 dscale = nmask*gmean;
268 }
269
270 if( opts->verb )
271 INFO_message("global mean = %f, mean diff = %f, mdiff/gmean = %f",
272 gmean, mdiff, mdiff/gmean);
273
274 /* scale and store the result */
275 for( tind=1; tind < nt; tind++ )
276 opts->result[tind] = dwork[tind] / dscale;
277 opts->nt = nt;
278
279 free(dwork);
280 free(fdata);
281
282 if( opts->verb > 1 ) INFO_message("successfully computed mean diff");
283
284 RETURN(0);
285 }
286
287
288 /* count voxels or mask fraction that hits a global maximum of 4095
289 *
290 * method = 4095_count : return masked counts
291 * 4095_frac : return masked fractions
292 * 4095_warn : return limit warning
293 *
294 * return 1 on error, 0 otherwise
295 */
compute_4095(options_t * opts,int method)296 int compute_4095(options_t * opts, int method)
297 {
298 double * dwork;
299 float * fdata, gmax, fval;
300 byte * mask = opts->mask; /* save 6 characters... */
301 int nt, nvox, vind, tind, nmask, found, nlocal;
302
303 ENTRY("compute_4095");
304
305 nt = DSET_NVALS(opts->inset);
306 nvox = DSET_NVOX(opts->inset);
307
308 /* note how many voxels this is over */
309 if( opts->mask ) nmask = THD_countmask( nvox, opts->mask );
310 else nmask = nvox;
311
312 if( opts->verb )
313 INFO_message("computing %s, nvox = %d, nmask = %d, nt = %d",
314 meth_index_to_name(method), nvox, nmask, nt);
315
316 /* allocate result early, in case there is nothing to do */
317 opts->result = calloc(nt, sizeof(float));
318
319 if( nmask == 0 ) {
320 ERROR_message("empty mask");
321 RETURN(0);
322 }
323
324 /* check_dims() has been called, so we have sufficient data */
325 fdata = calloc(nt, sizeof(float));
326
327 /* do the work */
328 gmax = 0.0;
329 found = 0;
330 for( vind=0; vind < nvox; vind++ ) {
331 if( opts->mask && ! opts->mask[vind] ) continue;
332
333 if( THD_extract_array(vind, opts->inset, 0, fdata) ) {
334 ERROR_message("failed to exract data at index %d\n", vind);
335 free(fdata); RETURN(1);
336 }
337
338 /* start counting */
339 for( tind=0; tind < nt; tind++ ) {
340 fval = fdata[tind];
341
342 /* track global max */
343 if( fval > gmax ) gmax = fval;
344
345 /* if we pass the limit, quit working */
346 if( gmax > 4095.0 ) break;
347
348 /* track hits */
349 if( fval == 4095.0 ) opts->result[tind]++;
350 }
351 }
352
353 /* no longer needed */
354 free(fdata);
355
356 /* babble to user */
357 if( opts->verb )
358 INFO_message("global max = %f", gmax);
359
360 /* if warning, decide and return */
361 if( method == T21_METH_4095_WARN ) {
362 if( gmax != 4095 ) {
363 if( opts->verb ) INFO_message("max of %g is okay", gmax);
364 else { printf("0\n"); fflush(stdout); }
365 RETURN(0);
366 }
367
368 if( opts->verb ) WARNING_message("suspicious max of exactly 4095");
369 else { printf("1\n"); fflush(stdout); }
370
371 RETURN(0);
372 }
373
374 /* if frac, scale */
375 if( method == T21_METH_4095_FRAC ) {
376 for( tind=0; tind < nt; tind++ ) opts->result[tind] /= nmask;
377 }
378
379 if( opts->verb > 1 ) INFO_message("successfully ran 4095 test");
380
381 RETURN(0);
382 }
383
384
385 /* this is basically enorm, with scaling variants
386 *
387 * convert dsum to dmean and then a scalar:
388 * method = 1 : enorm sqrt(ss)
389 * 2 : rms sqrt(ss/nvox) = stdev (biased)
390 * 3 : srms sqrt(ss/nvox) / grand mean
391 */
compute_enorm(options_t * opts,int method)392 int compute_enorm(options_t * opts, int method)
393 {
394 double * dwork, dscale=1.0, gmean, mdiff;
395 float * fdata, fdiff;
396 byte * mask = opts->mask; /* save 6 characters... */
397 int nt, nvox, vind, tind, nmask;
398
399 ENTRY("compute_enorm");
400
401 nt = DSET_NVALS(opts->inset);
402 nvox = DSET_NVOX(opts->inset);
403
404 /* note how many voxels this is over */
405 if( opts->mask ) nmask = THD_countmask( nvox, opts->mask );
406 else nmask = nvox;
407
408 if( opts->verb )
409 INFO_message("computing %s, nvox = %d, nmask = %d, nt = %d",
410 meth_index_to_name(method), nvox, nmask, nt);
411
412 /* check_dims() has been called, so we have sufficient data */
413
414 dwork = calloc(nt, sizeof(double));
415 fdata = calloc(nt, sizeof(float));
416
417 /* could steal fdata, but that might be unethical (plus, garbage on err) */
418 opts->result = calloc(nt, sizeof(float));
419
420 /* use gmean to get mean across all masked voxels and time */
421 gmean = 0.0;
422 mdiff = 0.0;
423 for( vind=0; vind < nvox; vind++ ) {
424 if( opts->mask && ! opts->mask[vind] ) continue;
425
426 if( THD_extract_array(vind, opts->inset, 0, fdata) ) {
427 ERROR_message("failed to exract data at index %d\n", vind);
428 free(dwork); free(fdata); RETURN(1);
429 }
430
431 /* accumuate squared differences; dwork[0] is already 0 */
432 gmean += fdata[0]; /* and accumlate for mean */
433 for( tind=1; tind < nt; tind++ ) {
434 gmean += fdata[tind]; /* accumlate for mean */
435 fdiff = fdata[tind]-fdata[tind-1];
436 mdiff += fabs(fdiff); /* accumulate for mean diff */
437 dwork[tind] += fdiff*fdiff;
438 }
439 }
440
441 /* convert dsum to dmean and then a scalar:
442 * method = 1 : enorm sqrt(ss)
443 * 2 : rms sqrt(ss/nmask) = stdev (biased) of first diffs
444 * 3 : srms sqrt(ss/nmask) / abs(grand mean)
445 */
446 gmean = fabs(gmean/nmask/nt); /* global masked mean */
447 mdiff /= nmask*nt; /* global masked mean first diff */
448
449 /* set the enorm scalar, based on the method */
450 if( method == T21_METH_ENORM ) {
451 if( opts->verb ) INFO_message("writing enorm : uses no scaling");
452 dscale = 1.0;
453 } else if ( method == T21_METH_RMS ) {
454 if( opts->verb ) INFO_message("writing rms = dvars = enorm/sqrt(nvox)");
455 dscale = sqrt(nmask);
456 } else if ( method == T21_METH_SRMS || method == T21_METH_S_SRMS ) {
457 if( opts->verb ) INFO_message("scaled dvars = srms = rms/gmean");
458 if( gmean == 0.0 ) {
459 ERROR_message("values have zero mean, failing (use rms, instead)");
460 free(dwork); free(fdata); RETURN(1);
461 } else if( gmean < 1.0 ) {
462 WARNING_message("values might be de-meaned, consider rms");
463 }
464 dscale = sqrt(nmask)*gmean;
465 }
466
467 /* babble to user */
468 if( opts->verb )
469 INFO_message("global mean = %f, mean diff = %f, mdiff/gmean = %f",
470 gmean, mdiff, mdiff/gmean);
471 if( opts->verb > 2 ) INFO_message("scaling enorm down by %f", dscale);
472
473 /* scale and store the result */
474 for( tind=1; tind < nt; tind++ )
475 opts->result[tind] = sqrt(dwork[tind]) / dscale;
476 opts->nt = nt;
477
478 /* if computing a mean diff, apply the shift */
479 if( method == T21_METH_S_SRMS ) {
480 mdiff /= gmean;
481 if( opts->verb ) INFO_message("shift by scaled mean diff %f", mdiff);
482 for( tind=1; tind < nt; tind++ )
483 opts->result[tind] -= mdiff;
484 }
485
486 free(dwork);
487 free(fdata);
488
489 if( opts->verb > 1 ) INFO_message("successfully computed enorm");
490
491 RETURN(0);
492 }
493
494
fill_mask(options_t * opts)495 int fill_mask(options_t * opts)
496 {
497 THD_3dim_dataset * mset;
498 int nvox;
499
500 ENTRY("fill_mask");
501
502 if( opts->automask ) {
503 if( opts->verb ) INFO_message("creating automask...");
504
505 opts->mask = THD_automask(opts->inset);
506 if( ! opts->mask ) {
507 ERROR_message("failed to apply -automask");
508 RETURN(1);
509 }
510
511 RETURN(0);
512 }
513
514 if( opts->mask_name ) {
515 if( opts->verb )
516 INFO_message("reading mask dset from %s...", opts->mask_name);
517
518 mset = THD_open_dataset( opts->mask_name );
519 if( ! mset ) ERROR_exit("cannot open mask dset '%s'", opts->mask_name);
520 nvox = DSET_NVOX(opts->inset);
521 if( DSET_NVOX(mset) != nvox ) {
522 ERROR_message("mask does not have the same voxel count as input");
523 RETURN(1);
524 }
525
526 /* fill mask array and mask_nxyz, remove mask dset */
527 DSET_load(mset); CHECK_LOAD_ERROR(mset);
528
529 opts->mask = THD_makemask(mset, 0, 1, 0);
530 DSET_delete(mset);
531
532 if( ! opts->mask ) {
533 ERROR_message("cannot make mask from '%s'", opts->mask_name);
534 RETURN(1);
535 }
536
537 if( opts->verb > 1 )
538 INFO_message("have mask with %d voxels", nvox);
539 }
540
541 RETURN(0);
542 }
543
544
545 /* ----------------------------------------------------------------------*/
546 /* method name to and from index */
meth_name_to_index(char * name)547 int meth_name_to_index(char * name)
548 {
549 if( ! name ) return T21_METH_UNDEF;
550
551 if( ! strcasecmp(name, "enorm") ) return T21_METH_ENORM;
552
553 if( ! strcasecmp(name, "rms") ||
554 ! strcasecmp(name, "dvars") ) return T21_METH_RMS;
555
556 if( ! strcasecmp(name, "srms") ||
557 ! strcasecmp(name, "cvar") ) return T21_METH_SRMS;
558
559 if( ! strcasecmp(name, "shift_srms") ||
560 ! strcasecmp(name, "s_srms") ) return T21_METH_S_SRMS;
561
562 if( ! strcasecmp(name, "mdiff") ) return T21_METH_MDIFF;
563
564 if( ! strcasecmp(name, "smdiff") ) return T21_METH_SMDIFF;
565
566 if( ! strcasecmp(name, "4095_count") ) return T21_METH_4095_COUNT;
567 if( ! strcasecmp(name, "4095_frac") ) return T21_METH_4095_FRAC;
568 if( ! strcasecmp(name, "4095_warn") ) return T21_METH_4095_WARN;
569
570 /* be explicit, since we have a case for this */
571 if( ! strcmp(name, "undefined") ) return T21_METH_UNDEF;
572
573 return T21_METH_UNDEF;
574 }
575
meth_index_to_name(int method)576 char * meth_index_to_name(int method)
577 {
578 if( method == T21_METH_ENORM ) return "enorm";
579 if( method == T21_METH_RMS ) return "rms";
580 if( method == T21_METH_SRMS ) return "srms";
581 if( method == T21_METH_S_SRMS ) return "shift_srms";
582 if( method == T21_METH_MDIFF ) return "mdiff";
583 if( method == T21_METH_SMDIFF ) return "smdiff";
584 if( method == T21_METH_4095_COUNT ) return "4095_count";
585 if( method == T21_METH_4095_FRAC ) return "4095_frac";
586 if( method == T21_METH_4095_WARN ) return "4095_warn";
587
588 return "undefined";
589 }
590
591 /* ----------------------------------------------------------------------*/
show_help(void)592 int show_help(void)
593 {
594 ENTRY("show_help");
595
596 printf(
597 "-------------------------------------------------------------------------\n"
598 "3dTto1D - collapse a 4D time series to a 1D time series\n"
599 "\n"
600 "The program takes as input a 4D times series (a list of 1D time series)\n"
601 "and optionally a mask, and computes from in a 1D time series using some\n"
602 "method applied to the first (backward) differences. Methods include:\n"
603 "\n"
604 " enorm : the Euclidean norm\n"
605 " rms/dvars : root mean square (DVARS)\n"
606 " srms : rms scaled down by global mean \n"
607 " shift_srms : srms shifted by the global mean\n"
608 " mdiff : mean abs(diff)\n"
609 " smdiff : mdiff scaled down by global mean\n"
610 " 4095_count : count voxels with max of exactly 4095\n"
611 " 4095_frac : fraction of masked voxels with max of exactly 4095\n"
612 " 4095_warn : warn if short datum and max of exactly 4095\n"
613 "\n"
614 "More details are provided after the examples.\n"
615 "\n"
616 "--------------------------------------------------\n"
617 "examples:\n"
618 "\n"
619 "E1. compute SRMS of EPI time series within a brain mask\n"
620 " (might be good for censoring, and is comparable across subjects)\n"
621 "\n"
622 " 3dTto1D -input epi_r1+orig -mask mask.auto.nii.gz -method srms \\\n"
623 " -prefix epi_srms.1D\n"
624 "\n"
625 "E2. compute DVARS of EPI time series within a brain mask\n"
626 " (similarly good for censoring, but not comparable across subjects)\n"
627 "\n"
628 " 3dTto1D -input epi_r1+orig -mask mask.auto.nii.gz -method dvars \\\n"
629 " -prefix epi_dvars.1D\n"
630 "\n"
631 "E3. compute ENORM of motion parameters\n"
632 " (as is done by afni_proc.py via 1d_tool.py)\n"
633 "\n"
634 " Note that 1D inputs will generally need the transpose operator,\n"
635 " applied by appending an escaped ' to the -input dataset.\n"
636 "\n"
637 " 3dTto1D -input dfile.r01.1D\\' -method enorm -prefix enorm.r01.1D\n"
638 "\n"
639 "E4. warn if short data and max is 4095\n"
640 "\n"
641 " 3dTto1D -input epi+orig -method 4095_warn\n"
642 "\n"
643 "--------------------------------------------------\n"
644 "methods:\n"
645 "\n"
646 " Since the initial step is generally to compute the first (backwards)\n"
647 " difference, call that dataset TDIFF. The value at any voxel of TDIFF\n"
648 " is the same as the input, minus the value at the prior time point.\n"
649 " TDIFF is defined as 0 at time point 0.\n"
650 "\n"
651 " method enorm\n"
652 "\n"
653 " This is the Euclidean norm.\n"
654 "\n"
655 " Starting with the TDIFF dataset, the value at each time point is\n"
656 " the Euclidean norm for that volume (or list of values). This is\n"
657 " the same as the L2-norm, and which is often applied to the motion\n"
658 " parameters for censoring.\n"
659 "\n"
660 " enorm = sqrt(sum squares)\n"
661 "\n"
662 "\n"
663 " method rms (or dvars)\n"
664 "\n"
665 " RMS = DVARS = enorm/sqrt(nvox).\n"
666 "\n"
667 " The RMS (root mean square) is the same as the enorm divided by\n"
668 " sqrt(nvox). It is like a standard deviation, but without removal\n"
669 " of the mean (per time point).\n"
670 "\n"
671 " rms = dvars = enorm/sqrt(nvox) = sqrt(sum squares/nvox)\n"
672 "\n"
673 " This is the RMS of backward differences first described by Smyser\n"
674 " et. al., 2010, for motion detection, and later renamed to DVARS by\n"
675 " Power et. al., 2012.\n"
676 "\n"
677 " * DVARS survives a resampling, where it would be unchanged if every\n"
678 " voxel were listed multiple times, for example.\n"
679 "\n"
680 " * DVARS does not survive a scaling, it scales with the data.\n"
681 " This is why the SRMS method was introduced.\n"
682 "\n"
683 "\n"
684 " method srms (or cvar) (= scaled rms = dvars/mean)\n"
685 "\n"
686 " This result is basically the coefficient of variation, but without\n"
687 " removal of each volume mean.\n"
688 "\n"
689 " This is the same as dvars divided by the global mean, gmean.\n"
690 "\n"
691 " srms = dvars/gmean = enorm/sqrt(nvox)/gmean\n"
692 "\n"
693 " * SRMS survives both a resampling and scaling of the data. Since it\n"
694 " is unchanged with any data scaling (unlike DVARS), values are\n"
695 " comparable across subjects and studies.\n"
696 "\n"
697 " *** The above 3 curves will look identical, subject to scaling.\n"
698 "\n"
699 "\n"
700 " method shift_srms (= srms - meandiff)\n"
701 "\n"
702 " This is simply the SRMS curve shifted down by the global mean of\n"
703 " (the absolute values of) the first differences. This is probably\n"
704 " useless.\n"
705 "\n"
706 "\n"
707 " method mdiff (mean diff = mean abs(first diff))\n"
708 "\n"
709 " Again, starting with the first backward difference, TDIFF, this\n"
710 " is just the mean absolute value, per time point.\n"
711 "\n"
712 "\n"
713 " method smdiff (scaled mean diff = mdiff/mean)\n"
714 "\n"
715 " This is the mean diff scaled by the global mean.\n"
716 "\n"
717 " method 4095_count\n"
718 "\n"
719 " At each time point, output the number of (masked) voxels that are\n"
720 " exactly 4095.\n"
721 "\n"
722 " method 4095_frac\n"
723 "\n"
724 " At each time point, output the fraction of (masked) voxels that\n"
725 " are exactly 4095.\n"
726 "\n"
727 " method 4095_warn\n"
728 "\n"
729 " Simply warn whether the maximum is exactly 4095 (so no -prefix).\n"
730 "\n"
731 "--------------------------------------------------\n"
732 "informational command arguments:\n"
733 "\n"
734 " -help : show this help\n"
735 " -hist : show program history\n"
736 " -ver : show program version\n"
737 "\n"
738 "--------------------------------------------------\n"
739 "required command arguments:\n"
740 "\n"
741 " -input DSET : specify input dataset\n"
742 "\n"
743 " e.g. -input epi_r1+orig\n"
744 " e.g. -input dfile.r01.1D\\'\n"
745 "\n"
746 " Specify the input dataset to be processed. This should be a set\n"
747 " of 3D time series. If the input is in 1D format, a transpose\n"
748 " operator will typically be required.\n"
749 "\n"
750 " -method METHOD : specify 4D to 1D conversion method\n"
751 "\n"
752 " e.g. -method srms\n"
753 " e.g. -method DVARS\n"
754 " e.g. -method dvars\n"
755 " e.g. -method enorm\n"
756 "\n"
757 " Details of the computational methods are at the top of the help.\n"
758 " The methods (which are case insensitive) include:\n"
759 "\n"
760 " enorm : Euclidean norm of first differences\n"
761 " = sqrt(sum squares(first diffs))\n"
762 "\n"
763 " rms : RMS (root mean square) of first differences\n"
764 " = DVARS = enorm/sqrt(nvox)\n"
765 "\n"
766 " srms : scaled (by grand mean) RMS of first differences\n"
767 " = DVARS/mean\n"
768 "\n"
769 " * seems like the most useful method for censoring\n"
770 "\n"
771 " s_srms : SRMS shifted by grand mean abs of first diffs\n"
772 " = SRMS - mean(abs(first diffs))\n"
773 "\n"
774 " mdiff : mean absolute first differences\n"
775 " = mean(abs(first diff))\n"
776 "\n"
777 " smdiff : mdiff scaled by grand mean\n"
778 " = mdiff/mean\n"
779 "\n"
780 " 4095_count : count of voxels that are exactly 4095\n"
781 "\n"
782 " 4095_frac : fraction of voxels that are exactly 4095\n"
783 " = 4095_count/(mask size)\n"
784 "\n"
785 " 4095_warn : state whether global max is exactly 4095\n"
786 " (no 1D output)\n"
787 "\n"
788 "--------------------------------------------------\n"
789 "optional command arguments:\n"
790 "\n"
791 " -automask : restrict computation to automask\n"
792 " -mask MSET : restrict computation to given mask\n"
793 " -prefix PREFIX : specify output file\n"
794 " e.g. -prefix SVAR_run1.1D\n"
795 " default: -prefix stdout\n"
796 " -verb LEVEL : specify verbose level\n"
797 " e.g. -verb 2\n"
798 " default: -verb 1\n"
799 "\n"
800 "--------------------------------------------------\n"
801 "R Reynolds July, 2017\n"
802 "-------------------------------------------------------------------------\n"
803 );
804
805 printf("%s\ncompiled: %s\n\n", g_version, __DATE__);
806
807 RETURN(0);
808 }
809
810
811 /* ----------------------------------------------------------------------
812 * fill the options_t struct
813 * return 1 on (acceptable) termination
814 * 0 on continue
815 * -1 on error
816 */
process_opts(options_t * opts,int argc,char * argv[])817 int process_opts(options_t * opts, int argc, char * argv[] )
818 {
819 int ac;
820
821 ENTRY("process_opts");
822
823 memset(opts, 0, sizeof(options_t)); /* init everything to 0 */
824 opts->method = T21_METH_UNDEF; /* be explicit, required option */
825 opts->verb = 1;
826
827 ac = 1;
828 while( ac < argc ) {
829
830 /* check for terminal options */
831 if( ! strcmp(argv[ac],"-help") || ! strcmp(argv[ac],"-h") ) {
832 show_help();
833 RETURN(1);
834 } else if( strcmp(argv[ac],"-hist") == 0 ) {
835 int c, len = sizeof(g_history)/sizeof(char *);
836 for( c = 0; c < len; c++) fputs(g_history[c], stdout);
837 putchar('\n');
838 RETURN(1);
839 } else if( strcmp(argv[ac],"-ver") == 0 ) {
840 puts(g_version);
841 RETURN(1);
842 }
843
844 /* the remaining options are alphabetical */
845
846 else if( strcmp(argv[ac],"-automask") == 0 ) {
847 opts->automask = 1;
848 ac++; continue;
849 }
850 else if( ! strcmp(argv[ac],"-input") || ! strcmp(argv[ac], "-inset") ) {
851 if( opts->inset ) ERROR_exit("-input already applied");
852
853 if( ++ac >= argc ) ERROR_exit("need argument after '-input'");
854
855 opts->inset = THD_open_dataset( argv[ac] );
856 opts->nt = DSET_NVALS(opts->inset);
857 if( ! opts->inset ) ERROR_exit("cannot open dset '%s'", argv[ac]);
858
859 DSET_load(opts->inset); CHECK_LOAD_ERROR(opts->inset);
860 ac++; continue;
861 }
862 else if( ! strcmp(argv[ac],"-mask") ) {
863 if( opts->mask_name ) ERROR_exit("-mask already applied");
864 if( ++ac >= argc ) ERROR_exit("need argument after '-mask'");
865
866 opts->mask_name = argv[ac];
867 ac++; continue;
868 }
869
870 else if( ! strcmp(argv[ac],"-method") ) {
871 if( ++ac >= argc ) ERROR_exit("need argument after '-method'");
872
873 opts->method = meth_name_to_index(argv[ac]);
874 if( opts->method == T21_METH_UNDEF )
875 ERROR_exit("illegal -method '%s'", argv[ac]);
876 ac++; continue;
877 }
878
879 else if( strcmp(argv[ac],"-prefix") == 0 ) {
880 if( ++ac >= argc ) ERROR_exit("need argument after '-prefix'");
881
882 opts->prefix = argv[ac];
883
884 if( !THD_filename_ok(opts->prefix) )
885 ERROR_exit("Illegal name after '-prefix'");
886 ac++; continue;
887 }
888
889 else if( strcmp(argv[ac],"-verb") == 0 ) {
890 if( ++ac >= argc ) ERROR_exit("need argument after '-verb'");
891
892 opts->verb = atoi(argv[ac]);
893 ac++; continue;
894 }
895
896 ERROR_message("** unknown option '%s'\n",argv[ac]);
897 RETURN(-1);
898 }
899
900 /* require data and method */
901 if( !opts->inset ) ERROR_exit("missing -input dataset");
902 if( opts->method == T21_METH_UNDEF ) ERROR_exit("missing -method option");
903
904 if( opts->verb > 1 ) {
905 INFO_message("input dataset loaded\n");
906 INFO_message("will apply method %s", meth_index_to_name(opts->method));
907 }
908
909 if( opts->automask && opts->mask_name )
910 ERROR_exit("cannot apply both -mask and -mset");
911
912 if( opts->automask || opts->mask_name ) {
913 if( fill_mask(opts) ) RETURN(1);
914 }
915
916 RETURN(0);
917 }
918
919