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 /*---------------------------------------------------------------------------
8 This program is revised for 3D+time data calculation,
9 [Raoqiong Tong, August 1997]
10
11 Added ability to use a 1D time series file as a "dataset" -- see TS variables.
12 [RW Cox, April 1998]
13
14 Added ability to operate on 3D bucket datasets.
15 [RW Cox, April 1998]
16
17 Added ability to use sub-brick selectors on input datasets.
18 [RW Cox, Jan 1999]
19
20 Modified output to scale each sub-brick to shorts/bytes separately
21 [RW Cox, Mar 1999]
22
23 Modifed sub-brick selection of type "-b3 name+view" to mangle dataset
24 into form "name+view[3]", since that code works better on 3D+time.
25 Modified TS_reader to use new mri_read_1D() function, instead of
26 mri_read_ascii().
27 Added -histpar option.
28 Added the _dshift stuff.
29 [RW Cox, Nov 1999]
30
31 Modified help menu
32 [P Christidis, July 2005]
33
34 Removed the '-b3' type of input from the help menu
35 [RW Cox, May 2010]
36 ----------------------------------------------------------------------------*/
37
38 #include "mrilib.h"
39 #include "parser.h"
40
41 #undef SHOW_B3
42
43 /*-------------------------- global data --------------------------*/
44
45 static int CALC_datum = ILLEGAL_TYPE ;
46 static int CALC_nvox = -1 ;
47 static PARSER_code * CALC_code = NULL ;
48 static int ntime[26] ;
49 static int ntime_max = 0 ;
50 static int CALC_fscale = 0 ; /* 16 Mar 1998 */
51 static int CALC_gscale = 0 ; /* 01 Apr 1999 */
52 static int CALC_nscale = 0 ; /* 15 Jun 2000 */
53
54 static int CALC_histpar = -1; /* 22 Nov 1999 */
55
56 static int CALC_usetemp = 0 ; /* 18 Oct 2005 */
57
58 #undef ALLOW_FDR /* this was purely experimental */
59 #ifdef ALLOW_FDR
60 static int CALC_fdrize = 0 ; /* 17 Jan 2008 */
61 #endif
62
63 #define ALLOW_SORT /* this is not experimental, but is pretty useless */
64 #ifdef ALLOW_SORT
65 static int CALC_sort = 0 ; /* 22 Jan 2008 */
66 #endif
67
68 static int CALC_do_isolas = 0 ; /* 27 Nov 2019 */
69
70 /*---------- dshift stuff [22 Nov 1999] ----------*/
71
72 #define DSHIFT_MODE_STOP 0
73 #define DSHIFT_MODE_WRAP 1
74 #define DSHIFT_MODE_ZERO 2
75
76 static int CALC_dshift [26] ; /* 22 Nov 1999 */
77 static int CALC_dshift_i [26] ;
78 static int CALC_dshift_j [26] ;
79 static int CALC_dshift_k [26] ;
80 static int CALC_dshift_l [26] ;
81 static int CALC_dshift_mode[26] ;
82
83 static int CALC_dshift_mode_current = DSHIFT_MODE_STOP ;
84 static int CALC_has_timeshift = 0 ;
85
86 /*------------------------------------------------*/
87
88 static int CALC_has_sym[26] ; /* 15 Sep 1999 */
89 static char abet[] = "abcdefghijklmnopqrstuvwxyz" ;
90
91 #define HAS_I CALC_has_sym[ 8]
92 #define HAS_J CALC_has_sym[ 9]
93 #define HAS_K CALC_has_sym[10]
94 #define HAS_X CALC_has_sym[23]
95 #define HAS_Y CALC_has_sym[24]
96 #define HAS_Z CALC_has_sym[25]
97 #define HAS_T CALC_has_sym[19] /* 19 Nov 1999 */
98 #define HAS_L CALC_has_sym[11] /* 19 Nov 1999 */
99 #define HAS_N CALC_has_sym[13] /* 29 Jul 2010 */
100
101 #define PREDEFINED_MASK ((1<< 8)|(1<< 9)|(1<<10)|(1<<11)| \
102 (1<<19)|(1<<23)|(1<<24)|(1<<25)|(1<<13) )
103
104 static int CALC_has_predefined = 0 ; /* 19 Nov 1999 */
105 static int CALC_has_xyz = 0 ; /* 17 May 2005 */
106 static int CALC_mangle_xyz = 0 ; /* 17 May 2005 */
107
108 #define MANGLE_NONE 0
109 #define MANGLE_RAI 1
110 #define MANGLE_LPI 2
111
112 static THD_3dim_dataset * CALC_dset[26] ;
113 static int CALC_type[26] ;
114 static byte ** CALC_byte[26] ;
115 static short ** CALC_short[26] ;
116 static float ** CALC_float[26] ;
117 static complex ** CALC_complex[26] ; /* 10 Mar 2006 */
118 static int CALC_cxcode[26] ;
119 static float * CALC_ffac[26] ;
120 static int CALC_noffac[26] ; /* 14 Nov 2003 */
121
122 static int CALC_verbose = 0 ; /* 30 April 1998 */
123
124 static char CALC_output_prefix[THD_MAX_NAME] = "calc" ;
125
126 static char CALC_session[THD_MAX_NAME] = "./" ;
127
128 static MRI_IMAGE * TS_flim[26] ; /* 17 Apr 1998 */
129 static float * TS_flar[26] ;
130 static int TS_nmax = 0 ;
131 static int TS_make = 0 ;
132 static float TS_dt = 1.0 ; /* 13 Aug 2001 */
133
134 static MRI_IMAGE * IJKAR_flim[26] ; /* 22 Feb 2005 */
135 static float * IJKAR_flar[26] ;
136 static int IJKAR_dcod[26] ;
137
138 /* this macro tells if a variable (index 0..25) is defined,
139 either by a time series file or an input dataset - 16 Nov 1999 */
140
141 #define VAR_DEFINED(kv) \
142 (TS_flim[kv] != NULL || IJKAR_flim[kv] != NULL || \
143 CALC_dset[kv] != NULL || CALC_dshift[kv] >= 0 )
144
145 static float Rfac = 0.299 ; /* 10 Feb 2002: for RGB inputs */
146 static float Gfac = 0.587 ;
147 static float Bfac = 0.114 ;
148
149 static int CALC_taxis_num = 0 ; /* 28 Apr 2003 */
150
151 #define CX_REALPART 0 /* 10 Mar 2006: complex to real methods */
152 #define CX_IMAGPART 1
153 #define CX_MAGNITUDE 2
154 #define CX_PHASE 3
155 static int CUR_cxcode = CX_MAGNITUDE ; /* default conversion method */
156
157 /*----------------------------------------------------------------------------*/
158 static char *tempfnam = NULL ; /* 18 Oct 2005: -usetemp stuff */
159 static FILE *tempfile = NULL ;
160
calc_atexit(void)161 static void calc_atexit(void) /*-- called by exit(): delete tempfile --*/
162 {
163 if( tempfile != NULL ){ fclose(tempfile) ; tempfile = NULL ; }
164 if( tempfnam != NULL ){
165 INFO_message("Deleting -usetemp file %s",tempfnam) ;
166 remove(tempfnam) ; tempfnam = NULL ;
167 }
168 return ;
169 }
170
171 /*--------------------------- prototypes ---------------------------*/
172 void CALC_read_opts( int , char ** ) ;
173 void CALC_Syntax(void) ;
174 int TS_reader( int , char * ) ;
175 int IJKAR_reader( int , char * ) ;
176
177 int remove_isolated_stuff( int nnx, int nny, int nnz, float *far, int maxite ) ;
178
179 /*--------------------------------------------------------------------
180 Read a time series file into TS variable number ival.
181 Returns -1 if an error occured, 0 otherwise.
182 ----------------------------------------------------------------------*/
183
TS_reader(int ival,char * fname)184 int TS_reader( int ival , char *fname )
185 {
186 MRI_IMAGE *tsim ;
187
188 if( ival < 0 || ival >= 26 ) return -1 ;
189
190 tsim = mri_read_1D( fname ) ; /* 16 Nov 1999: replaces mri_read_ascii */
191 if( tsim == NULL ) return -1 ;
192 if( tsim->nx < 2 ){ mri_free(tsim) ; return -1 ; }
193
194 TS_flim[ival] = tsim ;
195 TS_nmax = MAX( TS_nmax , TS_flim[ival]->nx ) ;
196 TS_flar[ival] = MRI_FLOAT_PTR( TS_flim[ival] ) ;
197 return 0 ;
198 }
199
200 /*--------------------------------------------------------------------
201 Read a time series file into IJK variable number ival.
202 Returns -1 if an error occured, 0 otherwise.
203 ----------------------------------------------------------------------*/
204
IJKAR_reader(int ival,char * fname)205 int IJKAR_reader( int ival , char *fname ) /* 22 Feb 2005 */
206 {
207 MRI_IMAGE *tsim ;
208
209 if( ival < 0 || ival >= 26 ) return -1 ;
210
211 tsim = mri_read_1D( fname ) ; /* 16 Nov 1999: replaces mri_read_ascii */
212 if( tsim == NULL ) return -1 ;
213 if( tsim->nx < 2 ){ mri_free(tsim) ; return -1 ; }
214
215 IJKAR_flim[ival] = tsim ;
216 IJKAR_flar[ival] = MRI_FLOAT_PTR( IJKAR_flim[ival] ) ;
217 return 0 ;
218 }
219
220 /*--------------------------------------------------------------------
221 read the arguments, load the global variables
222 ----------------------------------------------------------------------*/
223
CALC_read_opts(int argc,char * argv[])224 void CALC_read_opts( int argc , char * argv[] )
225 {
226 int nopt = 1 ;
227 int ids ;
228 int ii, kk;
229 char *srcspace=NULL, *eptr;
230
231 /* try to init CALC_mangle_xyz from env, first 27 Aug 2014 [rickr] */
232 eptr = my_getenv("AFNI_ORIENT");
233 if( eptr ) {
234 if ( ! strcasecmp(eptr, "rai") ) CALC_mangle_xyz = MANGLE_RAI;
235 else if( ! strcasecmp(eptr, "lpi") ) CALC_mangle_xyz = MANGLE_LPI;
236 }
237
238 for( ids=0 ; ids < 26 ; ids++ ){
239 CALC_dset[ids] = NULL ;
240 CALC_type[ids] = -1 ;
241 TS_flim[ids] = NULL ;
242 IJKAR_flim[ids] = NULL ; /* 22 Feb 2005 */
243
244 CALC_dshift[ids] = -1 ; /* 22 Nov 1999 */
245 CALC_dshift_mode[ids] = CALC_dshift_mode_current ;
246
247 CALC_noffac[ids] = 1 ; /* 14 Nov 2003 */
248 }
249
250 while( nopt < argc && argv[nopt][0] == '-' ){
251 #ifdef USE_TRACING
252 if( strncmp(argv[nopt],"-trace",5) == 0 ){ DBG_trace=1; nopt++; continue; }
253 if( strncmp(argv[nopt],"-TRACE",5) == 0 ){ DBG_trace=2; nopt++; continue; }
254 #endif
255
256 /**** -help ****/
257 if (!strcmp(argv[nopt], "-help")) {
258 CALC_Syntax();
259 exit(0);
260 }
261
262 /**** -isola [27 Nov 2019] ****/
263
264 if( strncasecmp(argv[nopt],"-isola",6) == 0 ){
265 if( strstr(argv[nopt],"-ISOLA") == 0 ) CALC_do_isolas = 99 ;
266 else CALC_do_isolas++ ;
267 nopt++ ; continue ;
268 }
269
270 /**** -dicom, -RAI, -LPI, -SPM [18 May 2005] ****/
271
272 if( strcasecmp(argv[nopt],"-dicom") == 0 || strcasecmp(argv[nopt],"-rai") == 0 ){
273 CALC_mangle_xyz = MANGLE_RAI ;
274 nopt++ ; continue ;
275 }
276
277 if( strcasecmp(argv[nopt],"-spm") == 0 || strcasecmp(argv[nopt],"-lpi") == 0 ){
278 CALC_mangle_xyz = MANGLE_LPI ;
279 nopt++ ; continue ;
280 }
281
282 /**** -rgbfac r g b [10 Feb 2003] ****/
283
284 if( strncasecmp(argv[nopt],"-rgbfac",7) == 0 ){
285 if( ++nopt >= argc )
286 ERROR_exit("need an argument after -rgbfac!\n") ;
287 Rfac = strtod( argv[nopt++] , NULL ) ;
288 Gfac = strtod( argv[nopt++] , NULL ) ;
289 Bfac = strtod( argv[nopt++] , NULL ) ;
290 if( Rfac == 0.0 && Gfac == 0.0 && Bfac == 0.0 )
291 ERROR_exit("All 3 factors after -rgbfac are zero!?\n");
292 continue ;
293 }
294
295 #ifdef ALLOW_FDR
296 /**** -fdrize ****/ /* [[ not in -help at this time !! ]] */
297
298 if( strcasecmp(argv[nopt],"-fdrize") == 0 ){ /* 17 Jan 2008 */
299 CALC_fdrize++ ; CALC_datum = MRI_float ; nopt++ ; continue ;
300 }
301 #endif
302
303 #ifdef ALLOW_SORT
304 /**** -sort ****/
305
306 if( strcmp(argv[nopt],"-sort") == 0 ){ /* 22 Jan 2008 */
307 CALC_sort = 1 ; nopt++ ; continue ;
308 }
309 if( strcmp(argv[nopt],"-SORT") == 0 ){
310 CALC_sort = -1 ; nopt++ ; continue ;
311 }
312 #endif
313
314 /**** -cx2r code [10 Mar 2006] ***/
315
316 if( strncasecmp(argv[nopt],"-cx2r",5) == 0 ){
317 if( ++nopt >= argc )
318 ERROR_exit("need an argument after -cx2r!\n") ;
319 if( strncasecmp(argv[nopt],"real",4) == 0 ) CUR_cxcode=CX_REALPART ;
320 else if( strncasecmp(argv[nopt],"imag",4) == 0 ) CUR_cxcode=CX_IMAGPART ;
321 else if( strncasecmp(argv[nopt],"mag",3) == 0 ||
322 strncasecmp(argv[nopt],"abs",3) == 0 ) CUR_cxcode=CX_MAGNITUDE;
323 else if( strncasecmp(argv[nopt],"pha",3) == 0 ||
324 strncasecmp(argv[nopt],"arc",3) == 0 ) CUR_cxcode=CX_PHASE ;
325 else {
326 CUR_cxcode=CX_MAGNITUDE;
327 WARNING_message("Don't understand '-cx2r %s' - using ABS",argv[nopt]);
328 }
329 nopt++ ; continue ;
330 }
331
332 /**** -taxis N:dt [28 Apr 2003] ****/
333
334 if( strncasecmp(argv[nopt],"-taxis",6) == 0 ){
335 char *cpt ;
336 if( ++nopt >= argc )
337 ERROR_exit("need an argument after -taxis!\n") ;
338 CALC_taxis_num = strtod( argv[nopt] , &cpt ) ;
339 if( CALC_taxis_num < 2 )
340 ERROR_exit("N value after -taxis must be bigger than 1!\n");
341 if( *cpt == ':' ){
342 float dt = strtod( cpt+1 , &cpt ) ;
343 if( dt > 0.0 ){
344 TS_dt = dt ;
345 if( *cpt == 'm' && *(cpt+1) == 's' ) TS_dt *= 0.001 ; /* 09 Mar 2004 */
346 } else {
347 WARNING_message("time step value in '-taxis %s' not legal!\n",argv[nopt]);
348 }
349 }
350 nopt++ ; continue ; /* go to next arg */
351 }
352
353 /**** -usetemp [18 Oct 2005] ****/
354
355 if( strncmp(argv[nopt],"-usetemp",6) == 0 ){
356 CALC_usetemp = 1 ;
357 nopt++ ; continue ; /* go to next arg */
358 }
359
360 /**** -dt val [13 Aug 2001] ****/
361
362 if( strncasecmp(argv[nopt],"-dt",3) == 0 || strncmp(argv[nopt],"-TR",3) == 0 ){
363 char *cpt ;
364 if( ++nopt >= argc )
365 ERROR_exit("need an argument after -dt!\n") ;
366 TS_dt = strtod( argv[nopt] , &cpt ) ;
367 if( TS_dt <= 0.0 )
368 ERROR_exit("Illegal time step value after -dt!\n");
369 if( *cpt == 'm' && *(cpt+1) == 's' ) TS_dt *= 0.001 ; /* 09 Mar 2004 */
370 nopt++ ; continue ; /* go to next arg */
371 }
372
373 /**** -histpar letter [22 Nov 1999] ****/
374
375 if( strncasecmp(argv[nopt],"-histpar",5) == 0 ){
376 if( ++nopt >= argc )
377 ERROR_exit("need an argument after -histpar!\n") ;
378 if( argv[nopt][0] < 'a' || argv[nopt][0] > 'z')
379 ERROR_exit("argument after -histpar is illegal!\n");
380 CALC_histpar = (int) (argv[nopt][0] - 'a') ;
381
382 nopt++ ; continue ; /* go to next arg */
383 }
384
385 /**** -float and -short and -byte [22 Jan 2008] ****/
386
387 if( strcasecmp(argv[nopt],"-float") == 0 ){
388 CALC_datum = MRI_float ; nopt++ ; continue ;
389 }
390 if( strcasecmp(argv[nopt],"-short") == 0 ){
391 CALC_datum = MRI_short ; nopt++ ; continue ;
392 }
393 if( strcasecmp(argv[nopt],"-byte") == 0 ){
394 CALC_datum = MRI_byte ; nopt++ ; continue ;
395 }
396
397 /**** -datum type ****/
398
399 if( strncasecmp(argv[nopt],"-datum",6) == 0 ){
400 if( ++nopt >= argc )
401 ERROR_exit("need an argument after -datum!\n") ;
402 if( strcasecmp(argv[nopt],"short") == 0 ){
403 CALC_datum = MRI_short ;
404 } else if( strcasecmp(argv[nopt],"float") == 0 ){
405 CALC_datum = MRI_float ;
406 } else if( strcasecmp(argv[nopt],"byte") == 0 ){
407 CALC_datum = MRI_byte ;
408 #if 0
409 } else if( strcasecmp(argv[nopt],"complex") == 0 ){ /* not listed in help */
410 CALC_datum = MRI_complex ;
411 #endif
412 } else {
413 ERROR_exit("-datum of type '%s' not supported in 3dcalc!\n",argv[nopt]) ;
414 }
415 nopt++ ; continue ; /* go to next arg */
416 }
417
418 /**** -verbose [30 April 1998] ****/
419
420 if( strncasecmp(argv[nopt],"-verbose",5) == 0 ){
421 CALC_verbose = 1 ;
422 nopt++ ; continue ;
423 }
424
425 /**** -nscale [15 Jun 2000] ****/
426
427 if( strncasecmp(argv[nopt],"-nscale",6) == 0 ){
428 CALC_gscale = CALC_fscale = 0 ;
429 CALC_nscale = 1 ;
430 nopt++ ; continue ;
431 }
432
433 /**** -fscale [16 Mar 1998] ****/
434
435 if( strncasecmp(argv[nopt],"-fscale",6) == 0 ){
436 CALC_fscale = 1 ;
437 CALC_nscale = 0 ;
438 nopt++ ; continue ;
439 }
440
441 /**** -gscale [01 Apr 1999] ****/
442
443 if( strncasecmp(argv[nopt],"-gscale",6) == 0 ){
444 CALC_gscale = CALC_fscale = 1 ;
445 CALC_nscale = 0 ;
446 nopt++ ; continue ;
447 }
448
449 /**** -prefix prefix ****/
450
451 if( strncasecmp(argv[nopt],"-prefix",6) == 0 ){
452 nopt++ ;
453 if( nopt >= argc )
454 ERROR_exit("need argument after -prefix!\n") ;
455 if( !THD_filename_ok(argv[nopt]) )
456 ERROR_exit("argument after -prefix is not a good name: '%s'",argv[nopt]) ;
457 MCW_strncpy( CALC_output_prefix , argv[nopt++] , THD_MAX_NAME ) ;
458 continue ;
459 }
460 /**** -session directory ****/
461
462 if( strncasecmp(argv[nopt],"-session",6) == 0 ){
463 nopt++ ;
464 if( nopt >= argc )
465 ERROR_exit("need argument after -session!\n") ;
466 MCW_strncpy( CALC_session , argv[nopt++] , THD_MAX_NAME ) ;
467 continue ;
468 }
469
470 /**** -expr expression ****/
471
472 if( strncasecmp(argv[nopt],"-expr",4) == 0 ){
473 if( CALC_code != NULL )
474 ERROR_exit("you cannot have 2 -expr options!\n") ;
475 nopt++ ;
476 if( nopt >= argc )
477 ERROR_exit("need argument after -expr!\n") ;
478 PARSER_set_printout(1) ; /* 21 Jul 2003 */
479 CALC_code = PARSER_generate_code( argv[nopt++] ) ;
480 if( CALC_code == NULL )
481 ERROR_exit("illegal expression -- see the help for details") ;
482 PARSER_mark_symbols( CALC_code , CALC_has_sym ) ; /* 15 Sep 1999 */
483 continue ;
484 }
485
486 /**** -dsSTOP [22 Nov 1999] ****/
487
488 if( strncasecmp(argv[nopt],"-dsSTOP",6) == 0 ){
489 CALC_dshift_mode_current = DSHIFT_MODE_STOP ;
490 nopt++ ; continue ;
491 }
492
493 /**** -dsWRAP [22 Nov 1999] ****/
494
495 if( strncasecmp(argv[nopt],"-dsWRAP",6) == 0 ){
496 CALC_dshift_mode_current = DSHIFT_MODE_WRAP ;
497 nopt++ ; continue ;
498 }
499
500 /**** -dsZERO [22 Nov 1999] ****/
501
502 if( strncasecmp(argv[nopt],"-dsZERO",6) == 0 ){
503 CALC_dshift_mode_current = DSHIFT_MODE_ZERO ;
504 nopt++ ; continue ;
505 }
506
507 /**** -<letter>[number] dataset ****/
508
509 ids = strlen( argv[nopt] ) ;
510
511 if( (argv[nopt][1] >= 'a' && argv[nopt][1] <= 'z') &&
512 (ids == 2 ||
513 (ids > 2 && argv[nopt][2] >= '0' && argv[nopt][2] <= '9')) ){
514
515 int ival , nxyz , isub , ll ;
516 THD_3dim_dataset *dset ;
517
518 ival = argv[nopt][1] - 'a' ;
519 if( VAR_DEFINED(ival) )
520 ERROR_exit("Can't define %c symbol twice",argv[nopt][1]);
521
522 isub = (ids == 2) ? 0 : strtol(argv[nopt]+2,NULL,10) ;
523 if( isub < 0 )
524 ERROR_exit("Illegal negative sub-brick value: %s",argv[nopt]) ;
525
526 nopt++ ;
527 if( nopt >= argc )
528 ERROR_exit("need argument after %s",argv[nopt-1]);
529
530 /*-- 22 Feb 2005: allow for I:, J:, K: prefix --*/
531
532 ll = strlen(argv[nopt]) ;
533 if( ll >= 4 &&
534 strstr(argv[nopt],"1D") != NULL &&
535 argv[nopt][1] == ':' &&
536 (argv[nopt][0] == 'I' || argv[nopt][0] == 'i' ||
537 argv[nopt][0] == 'J' || argv[nopt][0] == 'j' ||
538 argv[nopt][0] == 'K' || argv[nopt][0] == 'k' ) ){
539
540 ll = IJKAR_reader( ival , argv[nopt]+2 ) ;
541 if( ll == 0 ){
542 switch( argv[nopt][0] ){
543 case 'I': case 'i': IJKAR_dcod[ival] = 8 ; break ;
544 case 'J': case 'j': IJKAR_dcod[ival] = 9 ; break ;
545 case 'K': case 'k': IJKAR_dcod[ival] = 10 ; break ;
546 }
547 nopt++ ; goto DSET_DONE ;
548 }
549 }
550
551 /*-- 17 Apr 1998: allow for a *.1D filename --*/
552
553 ll = strlen(argv[nopt]) ;
554 if( ll >= 4 && ( strstr(argv[nopt],".1D") != NULL ||
555 strstr(argv[nopt],"1D:") != NULL )
556 && strstr(argv[nopt],"'") == NULL ){
557
558 ll = TS_reader( ival , argv[nopt] ) ;
559 if( ll == 0 ){ nopt++ ; goto DSET_DONE ; }
560
561 /* get to here => something bad happened, so try it as a dataset */
562 }
563
564 /*-- 22 Nov 1999: allow for a differentially
565 subscripted name, as in "-b a[1,0,0,0]" --*/
566
567 if( (argv[nopt][0] >= 'a' && argv[nopt][0] <= 'z') && /* legal name */
568 ( (ll >= 3 && argv[nopt][1] == '[') || /* subscript */
569 (ll == 3 && /* OR */
570 (argv[nopt][1] == '+' || argv[nopt][1] == '-')) /* +- ijkl */
571 ) ){
572
573 int jds = argv[nopt][0] - 'a' ; /* actual dataset index */
574 int * ijkl ; /* array of subscripts */
575
576 /*- sanity checks -*/
577
578 if( ids > 2 )
579 ERROR_exit("Can't combine %s with differential subscripting %s",
580 argv[nopt-1],argv[nopt]) ;
581 if( CALC_dset[jds] == NULL )
582 ERROR_exit("Must define dataset %c before using it in %s",
583 argv[nopt][0] , argv[nopt] ) ;
584
585 /*- get subscripts -*/
586
587 if( argv[nopt][1] == '[' ){ /* format is [i,j,k,l] */
588 MCW_intlist_allow_negative(1) ;
589 ijkl = MCW_get_intlist( 9999 , argv[nopt]+1 ) ;
590 MCW_intlist_allow_negative(0) ;
591 if( ijkl == NULL || ijkl[0] != 4 )
592 ERROR_exit("Illegal differential subscripting %s\n",
593 argv[nopt] ) ;
594 } else { /* format is +i, -j, etc */
595 ijkl = (int *) malloc( sizeof(int) * 5 ) ;
596 ijkl[1] = ijkl[2] = ijkl[3] = ijkl[4] = 0 ; /* initialize */
597 switch( argv[nopt][2] ){
598 default:
599 ERROR_exit("Bad differential subscripting %s\n",argv[nopt]);
600
601 case 'i': ijkl[1] = (argv[nopt][1]=='+') ? 1 : -1 ; break ;
602 case 'j': ijkl[2] = (argv[nopt][1]=='+') ? 1 : -1 ; break ;
603 case 'k': ijkl[3] = (argv[nopt][1]=='+') ? 1 : -1 ; break ;
604 case 'l': ijkl[4] = (argv[nopt][1]=='+') ? 1 : -1 ; break ;
605 }
606 }
607
608 /*- more sanity checks -*/
609
610 if( ijkl[1]==0 && ijkl[2]==0 && ijkl[3]==0 && ijkl[4]==0 )
611 WARNING_message("differential subscript %s is all zero\n",argv[nopt]);
612
613 if( ntime[jds] == 1 && ijkl[4] != 0 ){
614 WARNING_message(
615 "differential subscript %s has nonzero time\n"
616 " + shift on base dataset with 1 sub-brick!\n"
617 " + Setting time shift to 0.\n" ,
618 argv[nopt] ) ;
619 ijkl[4] = 0 ;
620 }
621
622 /*- set values for later use -*/
623
624 CALC_dshift [ival] = jds ;
625 CALC_dshift_i[ival] = ijkl[1] ;
626 CALC_dshift_j[ival] = ijkl[2] ;
627 CALC_dshift_k[ival] = ijkl[3] ;
628 CALC_dshift_l[ival] = ijkl[4] ;
629
630 CALC_dshift_mode[ival] = CALC_dshift_mode_current ;
631 CALC_cxcode[ival] = CUR_cxcode ; /* 10 Mar 2006 */
632
633 CALC_has_timeshift = CALC_has_timeshift || (ijkl[4] != 0) ;
634
635 /*- time to trot, Bwana -*/
636
637 free(ijkl) ; nopt++ ; goto DSET_DONE ;
638
639 } /* end of _dshift */
640
641 /*-- meanwhile, back at the "normal" dataset opening ranch --*/
642
643 { char dname[512] ; /* 02 Nov 1999 */
644 char *fname = argv[nopt]; /* 8 May 2007 [rickr,dglen] */
645
646 if( ids > 2 ){ /* mangle name */
647 if( strstr(argv[nopt],"[") != NULL ){
648 ERROR_exit(
649 "Illegal combination of sub-brick specifiers: "
650 "%s %s\n" ,
651 argv[nopt-1] , argv[nopt] ) ;
652 }
653 sprintf(dname,"%s[%d]",argv[nopt++],isub) ; /* use sub-brick */
654 fname = dname ;
655 isub = 0 ; /* 0 of dname */
656 } else {
657 nopt++ ; /* don't mangle */
658 }
659 dset = THD_open_dataset( fname ) ; /* open it */
660 if( dset == NULL )
661 ERROR_exit("can't open dataset %s\n",fname) ;
662 }
663
664 /* set some parameters based on the dataset */
665
666 ntime[ival] = DSET_NVALS(dset) ;
667 if ( ids > 2 ) ntime[ival] = 1 ;
668 ntime_max = MAX( ntime_max, ntime[ival] );
669
670 nxyz = dset->daxes->nxx * dset->daxes->nyy * dset->daxes->nzz ;
671 if( CALC_nvox < 0 ){
672 CALC_nvox = nxyz ;
673 /* set space for any atlas datasets */
674 srcspace = THD_get_space(dset); /* update default space */
675 set_out_space(srcspace); /* make output space for mask dset */
676 }
677 else{
678 if(!equivalent_space(THD_get_space(dset))){
679 WARNING_message(
680 "Template space of dataset %s does not match space of first dataset\n"
681 "Continuing anyway, but be sure the datasets really do match\n",
682 argv[nopt-1]);
683 /* exit(1);*/
684 }
685 if( nxyz != CALC_nvox ){
686 ERROR_exit(
687 "dataset %s differs in size [%d voxels] from others [%d]\n"
688 " -->> Since calculations are done voxel-by-voxel and\n"
689 " volume-by-volume, datasets must match both in\n"
690 " number of voxels per volume, and number of volumes.\n"
691 " -->> Depending on what you are doing, program 3dresample\n"
692 " might be useful to you.\n" ,
693 argv[nopt-1] , nxyz , CALC_nvox ) ;
694 }
695 }
696
697 if( !DSET_datum_constant(dset) ){ /* 29 May 2003 */
698 float *far ;
699
700 WARNING_message("dataset %s has sub-bricks with different types\n"
701 " + ==> converting all sub-bricks to floats\n",
702 argv[nopt-1]);
703
704 DSET_mallocize(dset) ; DSET_load(dset) ; CHECK_LOAD_ERROR(dset) ;
705
706 for( ii=0 ; ii < ntime[ival] ; ii++ ){
707 if( DSET_BRICK_TYPE(dset,ii) != MRI_float ){
708 far = calloc( sizeof(float) , nxyz ) ;
709 if( far == NULL )
710 ERROR_exit("can't malloc space for conversion -- ran out of memory :(");
711 EDIT_coerce_scale_type( nxyz , DSET_BRICK_FACTOR(dset,ii) ,
712 DSET_BRICK_TYPE(dset,ii), DSET_ARRAY(dset,ii),
713 MRI_float , far ) ;
714 EDIT_substitute_brick( dset , ii , MRI_float , far ) ;
715 DSET_BRICK_FACTOR(dset,ii) = 0.0f ;
716 }
717 }
718 }
719
720 CALC_type[ival] = DSET_BRICK_TYPE(dset,isub) ;
721 CALC_dset[ival] = dset ;
722
723 /* load floating scale factors */
724 /* 14 Nov 2003: CALC_noffac[ival] signals there is no scale factor
725 (so can avoid the multiplication when loading values) */
726
727 CALC_ffac[ival] = (float *) malloc( sizeof(float) * ntime[ival] ) ;
728 if ( ntime[ival] == 1 ) {
729 CALC_ffac[ival][0] = DSET_BRICK_FACTOR( dset , isub) ;
730 if (CALC_ffac[ival][0] == 0.0 ) CALC_ffac[ival][0] = 1.0 ;
731 if( CALC_ffac[ival][0] != 1.0 ) CALC_noffac[ival] = 0 ; /* 14 Nov 2003 */
732 } else {
733 for (ii = 0 ; ii < ntime[ival] ; ii ++ ) {
734 CALC_ffac[ival][ii] = DSET_BRICK_FACTOR(dset, ii) ;
735 if (CALC_ffac[ival][ii] == 0.0 ) CALC_ffac[ival][ii] = 1.0;
736 if( CALC_ffac[ival][ii] != 1.0 ) CALC_noffac[ival] = 0 ; /* 14 Nov 2003 */
737 }
738 }
739
740 /* read data from disk */
741
742 if( CALC_verbose ){
743 int iv , nb ;
744 for( iv=nb=0 ; iv < DSET_NVALS(dset) ; iv++ )
745 nb += DSET_BRICK_BYTES(dset,iv) ;
746 INFO_message("Reading dataset %s (%d bytes)\n",argv[nopt-1],nb);
747 }
748
749 if( ! DSET_LOADED(dset) ){
750 DSET_load(dset) ; CHECK_LOAD_ERROR(dset) ;
751 }
752
753 /* set pointers for actual dataset arrays */
754
755 CALC_cxcode[ival] = CUR_cxcode ; /* 10 Mar 2006 */
756
757 switch (CALC_type[ival]) {
758 case MRI_short:
759 CALC_short[ival] = (short **) malloc( sizeof(short *) * ntime[ival] ) ;
760 if (ntime[ival] == 1 )
761 CALC_short[ival][0] = (short *) DSET_ARRAY(dset, isub) ;
762 else
763 for (ii=0; ii < ntime[ival]; ii++)
764 CALC_short[ival][ii] = (short *) DSET_ARRAY(dset, ii);
765 break;
766
767 case MRI_float:
768 CALC_float[ival] = (float **) malloc( sizeof(float *) * ntime[ival] ) ;
769 if (ntime[ival] == 1 )
770 CALC_float[ival][0] = (float *) DSET_ARRAY(dset, isub) ;
771 else
772 for (ii=0; ii < ntime[ival]; ii++)
773 CALC_float[ival][ii] = (float *) DSET_ARRAY(dset, ii);
774 break;
775
776 case MRI_byte:
777 CALC_byte[ival] = (byte **) malloc( sizeof(byte *) * ntime[ival] ) ;
778 if (ntime[ival] == 1 )
779 CALC_byte[ival][0] = (byte *) DSET_ARRAY(dset, isub) ;
780 else
781 for (ii=0; ii < ntime[ival]; ii++)
782 CALC_byte[ival][ii] = (byte *) DSET_ARRAY(dset, ii);
783 break;
784
785 case MRI_rgb: /* 10 Feb 2003 */
786 CALC_byte[ival] = (byte **) malloc( sizeof(byte *) * ntime[ival] ) ;
787 if (ntime[ival] == 1 )
788 CALC_byte[ival][0] = (byte *) DSET_ARRAY(dset, isub) ;
789 else
790 for (ii=0; ii < ntime[ival]; ii++)
791 CALC_byte[ival][ii] = (byte *) DSET_ARRAY(dset, ii);
792 break ;
793
794 case MRI_complex: /* 10 Mar 2006 */
795 CALC_complex[ival] = (complex **)malloc( sizeof(complex *) * ntime[ival] );
796 if( ntime[ival] == 1 )
797 CALC_complex[ival][0] = (complex *) DSET_ARRAY(dset, isub) ;
798 else
799 for (ii=0; ii < ntime[ival]; ii++)
800 CALC_complex[ival][ii] = (complex *) DSET_ARRAY(dset, ii);
801 break ;
802
803 default:
804 ERROR_exit("Dataset %s has illegal data type for 3dcalc: %s\n" ,
805 argv[nopt-1] , MRI_type_name[CALC_type[ival]] ) ;
806
807 } /* end of switch over type switch */
808
809 /* if -datum not given or implied yet, set the output datum now */
810
811 if( CALC_datum < 0 && CALC_type[ival] != MRI_rgb ){
812 if( CALC_type[ival] == MRI_complex ) CALC_datum = MRI_float ;
813 else CALC_datum = CALC_type[ival] ;
814 }
815
816 DSET_DONE: continue; /*** target for various goto statements above ***/
817
818 } /* end of dataset input */
819
820 ERROR_message("Unknown option: %s\n",argv[nopt]) ;
821 suggest_best_prog_option(argv[0], argv[nopt]);
822 exit(1);
823 } /* end of loop over options */
824
825 /*---------------------------------------*/
826 /*** cleanup: check for various errors ***/
827
828 if( nopt < argc )
829 ERROR_exit(
830 "Extra command line arguments puzzle me! argv[%d]=%s ...\n",nopt,argv[nopt]) ;
831
832 if( CALC_gscale && CALC_usetemp ) /* 18 Oct 2005 */
833 ERROR_exit("-gscale and -usetemp are incompatible!") ;
834
835 for( ids=0 ; ids < 26 ; ids++ ) if( CALC_dset[ids] != NULL ) break ;
836 if( ids == 26 ){
837 for( ids=0 ; ids < 26 ; ids++ ) if( TS_flim[ids] != NULL ) break ;
838 if( ids < 26 )
839 ERROR_exit("No actual input datasets given! "
840 "Use '1deval' for .1D file calculations.");
841 else
842 ERROR_exit("No actual input datasets given!") ;
843 }
844
845 /* 22 Feb 2005: check IJKAR inputs against 1st dataset found */
846
847 for( ii=0 ; ii < 26 ; ii++ ){
848 if( IJKAR_flim[ii] != NULL ){
849 int siz=0 ;
850 switch( IJKAR_dcod[ii] ){
851 case 8: siz = DSET_NX(CALC_dset[ids]) ; break ;
852 case 9: siz = DSET_NY(CALC_dset[ids]) ; break ;
853 case 10: siz = DSET_NZ(CALC_dset[ids]) ; break ;
854 }
855 if( IJKAR_flim[ii]->nx != siz )
856 WARNING_message("dimension mismatch between '-%c' and '%-c'\n",
857 'a'+ii , 'a'+ids ) ;
858 }
859 }
860
861 if( CALC_code == NULL ) ERROR_exit("No expression given!\n") ;
862
863 if( CALC_histpar >= 0 && CALC_dset[CALC_histpar] == NULL ){
864 WARNING_message("-histpar dataset not defined!\n") ;
865 CALC_histpar = -1 ;
866 }
867
868 for (ids=0; ids < 26; ids ++)
869 if (ntime[ids] > 1 && ntime[ids] != ntime_max ) {
870 ERROR_exit("Multi-brick datasets don't have same number of volumes [%d]!\n"
871 " -->> Since calculations are done voxel-by-voxel and\n"
872 " volume-by-volume, datasets must match both in\n"
873 " number of voxels per volume, and number of volumes.",ntime_max);
874 }
875
876 /* 17 Apr 1998: if all input datasets are 3D only (no time),
877 and if there are any input time series,
878 then the output must become 3D+time itself */
879
880 if( ntime_max == 1 && TS_nmax > 0 ){
881 ntime_max = TS_nmax ;
882 TS_make = 1 ; /* flag to force manufacture of a 3D+time dataset */
883 INFO_message(
884 "Calculating 3D+time[%d]"
885 " dataset from 3D datasets and time series, with dt=%g s\n" ,
886 ntime_max , TS_dt ) ;
887 }
888
889 if( CALC_taxis_num > 0 ){ /* 28 Apr 2003 */
890 if( ntime_max > 1 ){
891 WARNING_message("-taxis %d overriden by dataset input(s)\n",
892 CALC_taxis_num) ;
893 } else {
894 ntime_max = CALC_taxis_num ;
895 TS_make = 1 ;
896 INFO_message("Calculating 3D+time[%d]"
897 " dataset from 3D datasets and -taxis with dt=%g s\n" ,
898 ntime_max , TS_dt ) ;
899 }
900 }
901
902 /* 15 Apr 1999: check if each input dataset is used,
903 or if an undefined symbol is used. */
904
905 for (ids=0; ids < 26; ids ++){
906 if( VAR_DEFINED(ids) && !CALC_has_sym[ids] )
907 WARNING_message("input '%c' is not used in the expression\n" ,
908 abet[ids] ) ;
909
910 else if( !VAR_DEFINED(ids) && CALC_has_sym[ids] ){
911
912 if( ((1<<ids) & PREDEFINED_MASK) == 0 ){
913 WARNING_message( "symbol %c is used but not defined\n" , abet[ids] ) ;
914 } else {
915 CALC_has_predefined++ ;
916 INFO_message("Symbol %c using predefined value\n",abet[ids] ) ;
917 if( ids >= 23 ) CALC_has_xyz = 1 ;
918 }
919 }
920 }
921
922 return ;
923 }
924
925 /*------------------------------------------------------------------*/
926
CALC_Syntax(void)927 void CALC_Syntax(void)
928 {
929 printf(
930 "Program: 3dcalc \n"
931 "Author: RW Cox et al \n"
932 "\n"
933 "3dcalc - AFNI's calculator program ~1~ \n"
934 "\n"
935 " This program does voxel-by-voxel arithmetic on 3D datasets \n"
936 " (only limited inter-voxel computations are possible). \n"
937 "\n"
938 " The program assumes that the voxel-by-voxel computations are being \n"
939 " performed on datasets that occupy the same space and have the same \n"
940 " orientations. \n"
941 "\n"
942 " 3dcalc has a lot of input options, as its capabilities have grown \n"
943 " over the years. So this 'help' output has gotten kind of long. \n"
944 "\n"
945 " For simple voxel-wise averaging of datasets: cf. 3dMean \n"
946 " For averaging along the time axis: cf. 3dTstat \n"
947 " For smoothing in time: cf. 3dTsmooth \n"
948 " For statistics from a region around each voxel: cf. 3dLocalstat \n"
949 "\n"
950 "------------------------------------------------------------------------\n"
951 "Usage: ~1~ \n"
952 "----- \n"
953 " 3dcalc -a dsetA [-b dsetB...] \\ \n"
954 " -expr EXPRESSION \\ \n"
955 " [options] \n"
956 "\n"
957 "Examples: ~1~ \n"
958 "-------- \n"
959 "1. Average datasets together, on a voxel-by-voxel basis: \n"
960 "\n"
961 " 3dcalc -a fred+tlrc -b ethel+tlrc -c lucy+tlrc \\\n"
962 " -expr '(a+b+c)/3' -prefix subjects_mean \n"
963 "\n"
964 " Averaging datasets can also be done by programs 3dMean and 3dmerge. \n"
965 " Use 3dTstat to averaging across sub-bricks in a single dataset. \n"
966 "\n"
967 "2. Perform arithmetic calculations between the sub-bricks of a single \n"
968 " dataset by noting the sub-brick number on the command line: \n"
969 "\n"
970 " 3dcalc -a 'func+orig[2]' -b 'func+orig[4]' -expr 'sqrt(a*b)' \n"
971 "\n"
972 "3. Create a simple mask that consists only of values in sub-brick #0 \n"
973 " that are greater than 3.14159: \n"
974 "\n"
975 " 3dcalc -a 'func+orig[0]' -expr 'ispositive(a-3.14159)' \\\n"
976 " -prefix mask \n"
977 "\n"
978 "4. Normalize subjects' time series datasets to percent change values in \n"
979 " preparation for group analysis: \n"
980 "\n"
981 " Voxel-by-voxel, the example below divides each intensity value in \n"
982 " the time series (epi_r1+orig) with the voxel's mean value (mean+orig)\n"
983 " to get a percent change value. The 'ispositive' command will ignore \n"
984 " voxels with mean values less than 167 (i.e., they are labeled as \n"
985 " 'zero' in the output file 'percent_change+orig') and are most likely \n"
986 " background/noncortical voxels. \n"
987 "\n"
988 " 3dcalc -a epi_run1+orig -b mean+orig \\\n"
989 " -expr '100 * a/b * ispositive(b-167)' -prefix percent_chng \n"
990 "\n"
991 "5. Create a compound mask from a statistical dataset, where 3 stimuli \n"
992 " show activation. \n"
993 " NOTE: 'step' and 'ispositive' are identical expressions that can \n"
994 " be used interchangeably: \n"
995 "\n"
996 " 3dcalc -a 'func+orig[12]' -b 'func+orig[15]' -c 'func+orig[18]' \\\n"
997 " -expr 'step(a-4.2)*step(b-2.9)*step(c-3.1)' \\\n"
998 " -prefix compound_mask \n"
999 "\n"
1000 " In this example, all 3 statistical criteria must be met at once for \n"
1001 " a voxel to be selected (value of 1) in this mask. \n"
1002 "\n"
1003 "6. Same as example #5, but this time create a mask of 8 different values\n"
1004 " showing all combinations of activations (i.e., not only where \n"
1005 " everything is active, but also each stimulus individually, and all \n"
1006 " combinations). The output mask dataset labels voxel values as such: \n"
1007 "\n"
1008 " 0 = none active 1 = A only active 2 = B only active \n"
1009 " 3 = A and B only 4 = C only active 5 = A and C only \n"
1010 " 6 = B and C only 7 = all A, B, and C active \n"
1011 "\n"
1012 " 3dcalc -a 'func+orig[12]' -b 'func+orig[15]' -c 'func+orig[18]' \\\n"
1013 " -expr 'step(a-4.2)+2*step(b-2.9)+4*step(c-3.1)' \\\n"
1014 " -prefix mask_8 \n"
1015 "\n"
1016 " In displaying such a binary-encoded mask in AFNI, you would probably \n"
1017 " set the color display to have 8 discrete levels (the '#' menu). \n"
1018 "\n"
1019 "7. Create a region-of-interest mask comprised of a 3-dimensional sphere.\n"
1020 " Values within the ROI sphere will be labeled as '1' while values \n"
1021 " outside the mask will be labeled as '0'. Statistical analyses can \n"
1022 " then be done on the voxels within the ROI sphere. \n"
1023 "\n"
1024 " The example below puts a solid ball (sphere) of radius 3=sqrt(9) \n"
1025 " about the point with coordinates (x,y,z)=(20,30,70): \n"
1026 "\n"
1027 " 3dcalc -a anat+tlrc \\\n"
1028 " -expr 'step(9-(x-20)*(x-20)-(y-30)*(y-30)-(z-70)*(z-70))' \\\n"
1029 " -prefix ball \n"
1030 "\n"
1031 " The spatial meaning of (x,y,z) is discussed in the 'COORDINATES' \n"
1032 " section of this help listing (far below). \n"
1033 "\n"
1034 "8. Some datsets are 'short' (16 bit) integers with a scalar attached, \n"
1035 " which allow them to be smaller than float datasets and to contain \n"
1036 " fractional values. \n"
1037 "\n"
1038 " Dataset 'a' is always used as a template for the output dataset. For \n"
1039 " the examples below, assume that datasets d1+orig and d2+orig consist \n"
1040 " of small integers. \n"
1041 "\n"
1042 " a) When dividing 'a' by 'b', the result should be scaled, so that a \n"
1043 " value of 2.4 is not truncated to '2'. To avoid this truncation, \n"
1044 " force scaling with the -fscale option: \n"
1045 "\n"
1046 " 3dcalc -a d1+orig -b d2+orig -expr 'a/b' -prefix quot -fscale \n"
1047 "\n"
1048 " b) If it is preferable that the result is of type 'float', then set \n"
1049 " the output data type (datum) to float: \n"
1050 "\n"
1051 " 3dcalc -a d1+orig -b d2+orig -expr 'a/b' -prefix quot \\\n"
1052 " -datum float \n"
1053 "\n"
1054 " c) Perhaps an integral division is desired, so that 9/4=2, not 2.24. \n"
1055 " Force the results not to be scaled (opposite of example 8a) using \n"
1056 " the -nscale option: \n"
1057 "\n"
1058 " 3dcalc -a d1+orig -b d2+orig -expr 'a/b' -prefix quot -nscale \n"
1059 "\n"
1060 "9. Compare the left and right amygdala between the Talairach atlas, \n"
1061 " and the CA_N27_ML atlas. The result will be 1 if TT only, 2 if CA \n"
1062 " only, and 3 where they overlap. \n"
1063 "\n"
1064 " 3dcalc -a 'TT_Daemon::amygdala' -b 'CA_N27_ML::amygdala' \\\n"
1065 " -expr 'step(a)+2*step(b)' -prefix compare.maps \n"
1066 "\n"
1067 " (see 'whereami -help' for more information on atlases) \n"
1068 "\n"
1069 "10. Convert a dataset from AFNI short format storage to NIfTI-1 floating\n"
1070 " point (perhaps for input to an non-AFNI program that requires this):\n"
1071 "\n"
1072 " 3dcalc -a zork+orig -prefix zfloat.nii -datum float -expr 'a' \n"
1073 "\n"
1074 " This operation could also be performed with program 3dAFNItoNIFTI. \n"
1075 "\n"
1076 "11. Compute the edge voxels of a mask dataset. An edge voxel is one \n"
1077 " that shares some face with a non-masked voxel. This computation \n"
1078 " assumes 'a' is a binary mask (particularly for 'amongst'). \n"
1079 "\n"
1080 " 3dcalc -a mask+orig -prefix edge \\\n"
1081 " -b a+i -c a-i -d a+j -e a-j -f a+k -g a-k \\\n"
1082 " -expr 'a*amongst(0,b,c,d,e,f,g)' \n"
1083 "\n"
1084 " consider similar erode or dilate operations: \n"
1085 " erosion: -expr 'a*(1-amongst(0,b,c,d,e,f,g))' \n"
1086 " dilation: -expr 'amongst(1,a,b,c,d,e,f,g)' \n"
1087 "\n"
1088 "------------------------------------------------------------------------\n"
1089 "ARGUMENTS for 3dcalc (must be included on command line): ~1~ \n"
1090 "--------- \n"
1091 "\n"
1092 " -a dname = Read dataset 'dname' and call the voxel values 'a' in the\n"
1093 " expression (-expr) that is input below. Up to 26 dnames \n"
1094 " (-a, -b, -c, ... -z) can be included in a single 3dcalc \n"
1095 " calculation/expression. \n"
1096 " ** If some letter name is used in the expression, but \n"
1097 " not present in one of the dataset options here, then \n"
1098 " that variable is set to 0. \n"
1099 #ifdef SHOW_B3
1100 " ** If the letter is followed by a number, then that \n"
1101 " number is used to select the sub-brick of the dataset \n"
1102 " which will be used in the calculations. \n"
1103 " E.g., '-b3 dname' specifies that the variable 'b' \n"
1104 " refers to sub-brick '3' of that dataset \n"
1105 " (indexes in AFNI start at 0). \n"
1106 #endif
1107 " ** You can use the subscript '[]' method \n"
1108 " to select sub-bricks of datasets, as in \n"
1109 " -b dname+orig'[3]' \n"
1110 #ifdef SHOW_B3
1111 " rather than the older notation \n"
1112 " -b3 dname+orig \n"
1113 " The subscript notation is more flexible, as it can \n"
1114 " be used to select a collection of sub-bricks. \n"
1115 #endif
1116 " ** If you just want to test some 3dcalc expression, \n"
1117 " you can supply a dataset 'name' of the form \n"
1118 " jRandomDataset:64,64,16,40 \n"
1119 " to have the program create and use a dataset \n"
1120 " with a 3D 64x64x16 grid, with 40 time points, \n"
1121 " filled with random numbers (uniform on [-1,1]). \n"
1122 "\n"
1123 " -expr = Apply the expression - within quotes - to the input \n"
1124 " datasets (dnames), one voxel at time, to produce the \n"
1125 " output dataset. \n"
1126 " ** You must use 1 and only 1 '-expr' option! \n"
1127 "\n"
1128 " NOTE: If you want to average or sum up a lot of datasets, programs \n"
1129 " 3dTstat and/or 3dMean and/or 3dmerge are better suited for these \n"
1130 " purposes. A common request is to increase the number of input \n"
1131 " datasets beyond 26, but in almost all cases such users simply \n"
1132 " want to do simple addition! \n"
1133 "\n"
1134 " NOTE: If you want to include shell variables in the expression (or in \n"
1135 " the dataset sub-brick selection), then you should use double \n"
1136 " \"quotes\" and the '$' notation for the shell variables; this \n"
1137 " example uses csh notation to set the shell variable 'z': \n"
1138 "\n"
1139 " set z = 3.5 \n"
1140 " 3dcalc -a moose.nii -prefix goose.nii -expr \"a*$z\"\n"
1141 "\n"
1142 " The shell will not expand variables inside single 'quotes', \n"
1143 " and 3dcalc's parser will not understand the '$' character. \n"
1144 "\n"
1145 " NOTE: You can use the ccalc program to play with the expression \n"
1146 " evaluator, in order to get a feel for how it works and \n"
1147 " what it accepts. \n"
1148 "\n"
1149 "------------------------------------------------------------------------\n"
1150 ) ;
1151 printf(
1152 " OPTIONS for 3dcalc: ~1~ \n"
1153 " ------- \n"
1154 "\n"
1155 " -help = Show this help.\n"
1156 "\n"
1157 " -verbose = Makes the program print out various information as it \n"
1158 " progresses. \n"
1159 "\n"
1160 " -datum type= Coerce the output data to be stored as the given type, \n"
1161 " which may be byte, short, or float. \n"
1162 " [default = datum of first input dataset] \n"
1163 " -float } \n"
1164 " -short } = Alternative options to specify output data format. \n"
1165 " -byte } \n"
1166 "\n"
1167 " -fscale = Force scaling of the output to the maximum integer \n"
1168 " range. This only has effect if the output datum is byte \n"
1169 " or short (either forced or defaulted). This option is \n"
1170 " often necessary to eliminate unpleasant truncation \n"
1171 " artifacts. \n"
1172 " [The default is to scale only if the computed values \n"
1173 " seem to need it -- are all <= 1.0 or there is at \n"
1174 " least one value beyond the integer upper limit.] \n"
1175 "\n"
1176 " ** In earlier versions of 3dcalc, scaling (if used) was \n"
1177 " applied to all sub-bricks equally -- a common scale \n"
1178 " factor was used. This would cause trouble if the \n"
1179 " values in different sub-bricks were in vastly \n"
1180 " different scales. In this version, each sub-brick \n"
1181 " gets its own scale factor. To override this behavior,\n"
1182 " use the '-gscale' option. \n"
1183 "\n"
1184 " -gscale = Same as '-fscale', but also forces each output sub-brick \n"
1185 " to get the same scaling factor. This may be desirable \n"
1186 " for 3D+time datasets, for example. \n"
1187 " ** N.B.: -usetemp and -gscale are incompatible!! \n"
1188 "\n"
1189 " -nscale = Don't do any scaling on output to byte or short datasets.\n"
1190 " This may be especially useful when operating on mask \n"
1191 " datasets whose output values are only 0's and 1's. \n"
1192 " ** Only use this option if you are sure you \n"
1193 " want the output dataset to be integer-valued! \n"
1194 "\n"
1195 " -prefix pname = Use 'pname' for the output dataset prefix name. \n"
1196 " [default='calc'] \n"
1197 "\n"
1198 " -session dir = Use 'dir' for the output dataset session directory. \n"
1199 " [default='./'=current working directory] \n"
1200 " You can also include the output directory in the \n"
1201 " 'pname' parameter to the -prefix option. \n"
1202 "\n"
1203 " -usetemp = With this option, a temporary file will be created to \n"
1204 " hold intermediate results. This will make the program\n"
1205 " run slower, but can be useful when creating huge \n"
1206 " datasets that won't all fit in memory at once. \n"
1207 " * The program prints out the name of the temporary \n"
1208 " file; if 3dcalc crashes, you might have to delete \n"
1209 " this file manually. \n"
1210 " ** N.B.: -usetemp and -gscale are incompatible!! \n"
1211 "\n"
1212 " -dt tstep = Use 'tstep' as the TR for \"manufactured\" 3D+time \n"
1213 " *OR* datasets. \n"
1214 " -TR tstep = If not given, defaults to 1 second. \n"
1215 "\n"
1216 " -taxis N = If only 3D datasets are input (no 3D+time or .1D files),\n"
1217 " *OR* then normally only a 3D dataset is calculated. With \n"
1218 " -taxis N:tstep: this option, you can force the creation of a time axis\n"
1219 " of length 'N', optionally using time step 'tstep'. In\n"
1220 " such a case, you will probably want to use the pre- \n"
1221 " defined time variables 't' and/or 'k' in your \n"
1222 " expression, or each resulting sub-brick will be \n"
1223 " identical. For example: \n"
1224 " '-taxis 121:0.1' will produce 121 points in time, \n"
1225 " spaced with TR 0.1. \n"
1226 "\n"
1227 " N.B.: You can also specify the TR using the -dt option. \n"
1228 " N.B.: You can specify 1D input datasets using the \n"
1229 " '1D:n@val,n@val' notation to get a similar effect. \n"
1230 " For example: \n"
1231 " -dt 0.1 -w '1D:121@0' \n"
1232 " will have pretty much the same effect as \n"
1233 " -taxis 121:0.1\n"
1234 " N.B.: For both '-dt' and '-taxis', the 'tstep' value is in \n"
1235 " seconds. \n"
1236 "\n"
1237 " -rgbfac A B C = For RGB input datasets, the 3 channels (r,g,b) are \n"
1238 " collapsed to one for the purposes of 3dcalc, using the\n"
1239 " formula value = A*r + B*g + C*b \n"
1240 "\n"
1241 " The default values are A=0.299 B=0.587 C=0.114, which \n"
1242 " gives the grayscale intensity. To pick out the Green \n"
1243 " channel only, use '-rgbfac 0 1 0', for example. Note \n"
1244 " that each channel in an RGB dataset is a byte in the \n"
1245 " range 0..255. Thus, '-rgbfac 0.001173 0.002302 0.000447'\n"
1246 " will compute the intensity rescaled to the range 0..1.0\n"
1247 " (i.e., 0.001173=0.299/255, etc.) \n"
1248 "\n"
1249 " -cx2r METHOD = For complex input datasets, the 2 channels must be \n"
1250 " converted to 1 real number for calculation. The \n"
1251 " methods available are: REAL IMAG ABS PHASE \n"
1252 " * The default method is ABS = sqrt(REAL^2+IMAG^2) \n"
1253 " * PHASE = atan2(IMAG,REAL) \n"
1254 " * Multiple '-cx2r' options can be given: \n"
1255 " when a complex dataset is given on the command line,\n"
1256 " the most recent previous method will govern. \n"
1257 " This also means that for -cx2r to affect a variable \n"
1258 " it must precede it. For example, to compute the \n"
1259 " phase of data in 'a' you should use \n"
1260 " 3dcalc -cx2r PHASE -a dft.lh.TS.niml.dset -expr 'a' \n"
1261 " However, the -cx2r option will have no effect in \n"
1262 " 3dcalc -a dft.lh.TS.niml.dset -cx2r PHASE -expr 'a' \n"
1263 " which will produce the default ABS of 'a' \n"
1264 " The -cx2r option in the latter example only applies \n"
1265 " to variables that will be defined after it. \n"
1266 " When in doubt, check your output. \n"
1267 " * If a complex dataset is used in a differential \n"
1268 " subscript, then the most recent previous -cx2r \n"
1269 " method applies to the extraction; for example \n"
1270 " -cx2r REAL -a cx+orig -cx2r IMAG -b 'a[0,0,0,0]' \n"
1271 " means that variable 'a' refers to the real part \n"
1272 " of the input dataset and variable 'b' to the \n"
1273 " imaginary part of the input dataset. \n"
1274 " * 3dcalc cannot be used to CREATE a complex dataset! \n"
1275 " [See program 3dTwotoComplex for that purpose.] \n"
1276 #ifdef ALLOW_SORT
1277 "\n"
1278 " -sort = Sort each output brick separately, before output: \n"
1279 " -SORT 'sort' ==> increasing order, 'SORT' ==> decreasing. \n"
1280 " [This is useful only under unusual circumstances!] \n"
1281 " [Sorting is done in spatial indexes, not in time.] \n"
1282 " [Program 3dTsort will sort voxels along time axis] \n"
1283 #endif
1284 "\n"
1285 " -isola = After computation, remove isolated non-zero voxels. \n"
1286 " This option can be repeated to iterate the process; \n"
1287 " each copy of '-isola' will cause the isola removal \n"
1288 " process to be repeated one more time.\n"
1289 "\n"
1290 "------------------------------------------------------------------------\n"
1291 "DATASET TYPES: ~1~ \n"
1292 "------------- \n"
1293 "\n"
1294 " The most common AFNI dataset types are 'byte', 'short', and 'float'. \n"
1295 "\n"
1296 " A byte value is an 8-bit signed integer (0..255), a short value ia a \n"
1297 " 16-bit signed integer (-32768..32767), and a float value is a 32-bit \n"
1298 " real number. A byte value has almost 3 decimals of accuracy, a short \n"
1299 " has almost 5, and a float has approximately 7 (from a 23+1 bit \n"
1300 " mantissa). \n"
1301 "\n"
1302 " Datasets can also have a scalar attached to each sub-brick. The main \n"
1303 " use of this is allowing a short type dataset to take on non-integral \n"
1304 " values, while being half the size of a float dataset. \n"
1305 "\n"
1306 " As an example, consider a short dataset with a scalar of 0.0001. This \n"
1307 " could represent values between -32.768 and +32.767, at a resolution of \n"
1308 " 0.001. One could represnt the difference between 4.916 and 4.917, for \n"
1309 " instance, but not 4.9165. Each number has 15 bits of accuracy, plus a \n"
1310 " sign bit, which gives 4-5 decimal places of accuracy. If this is not \n"
1311 " enough, then it makes sense to use the larger type, float. \n"
1312 "\n"
1313 "------------------------------------------------------------------------\n"
1314 "3D+TIME DATASETS: ~1~ \n"
1315 "---------------- \n"
1316 "\n"
1317 " This version of 3dcalc can operate on 3D+time datasets. Each input \n"
1318 " dataset will be in one of these conditions: \n"
1319 "\n"
1320 " (A) Is a regular 3D (no time) dataset; or \n"
1321 " (B) Is a 3D+time dataset with a sub-brick index specified ('[3]'); or\n"
1322 " (C) Is a 3D+time dataset with no sub-brick index specified ('-b'). \n"
1323 "\n"
1324 " If there is at least one case (C) dataset, then the output dataset will\n"
1325 " also be 3D+time; otherwise it will be a 3D dataset with one sub-brick. \n"
1326 " When producing a 3D+time dataset, datasets in case (A) or (B) will be \n"
1327 " treated as if the particular brick being used has the same value at each\n"
1328 " point in time. \n"
1329 "\n"
1330 " Multi-brick 'bucket' datasets may also be used. Note that if multi-brick\n"
1331 " (bucket or 3D+time) datasets are used, the lowest letter dataset will \n"
1332 " serve as the template for the output; that is, '-b fred+tlrc' takes \n"
1333 " precedence over '-c wilma+tlrc'. (The program 3drefit can be used to \n"
1334 " alter the .HEAD parameters of the output dataset, if desired.) \n"
1335 "\n"
1336 "------------------------------------------------------------------------\n"
1337 MASTER_HELP_STRING
1338 "\n"
1339 CATENATE_HELP_STRING
1340
1341 #ifdef SHOW_B3
1342 "\n"
1343 "** WARNING: you cannot combine sub-brick selection of the form \n"
1344 " -b3 bambam+orig (the old method) \n"
1345 " with sub-brick selection of the form \n"
1346 " -b 'bambam+orig[3]' (the new method) \n"
1347 " If you try, the Doom of Mandos will fall upon you! \n"
1348 #endif
1349 "\n"
1350 "------------------------------------------------------------------------\n"
1351 "1D TIME SERIES: ~1~ \n"
1352 "-------------- \n"
1353 "\n"
1354 " You can also input a '*.1D' time series file in place of a dataset. \n"
1355 " In this case, the value at each spatial voxel at time index n will be \n"
1356 " the same, and will be the n-th value from the time series file. \n"
1357 " At least one true dataset must be input. If all the input datasets \n"
1358 " are 3D (single sub-brick) or are single sub-bricks from multi-brick \n"
1359 " datasets, then the output will be a 'manufactured' 3D+time dataset. \n"
1360 "\n"
1361 " For example, suppose that 'a3D+orig' is a 3D dataset: \n"
1362 "\n"
1363 " 3dcalc -a a3D+orig -b b.1D -expr \"a*b\"\n"
1364 "\n"
1365 " The output dataset will 3D+time with the value at (x,y,z,t) being \n"
1366 " computed by a3D(x,y,z)*b(t). The TR for this dataset will be set \n"
1367 " to 'tstep' seconds -- this could be altered later with program 3drefit.\n"
1368 " Another method to set up the correct timing would be to input an \n"
1369 " unused 3D+time dataset -- 3dcalc will then copy that dataset's time \n"
1370 " information, but simply do not use that dataset's letter in -expr. \n"
1371 "\n"
1372 " If the *.1D file has multiple columns, only the first read will be \n"
1373 " used in this program. You can select a column to be the first by \n"
1374 " using a sub-vector selection of the form 'b.1D[3]', which will \n"
1375 " choose the 4th column (since counting starts at 0). \n"
1376 "\n"
1377 " '{...}' row selectors can also be used - see the output of '1dcat -help'\n"
1378 " for more details on these. Note that if multiple timeseries or 3D+time\n"
1379 " or 3D bucket datasets are input, they must all have the same number of \n"
1380 " points along the 'time' dimension. \n"
1381 "\n"
1382 " N.B.: To perform calculations ONLY on .1D files, use program 1deval. \n"
1383 " 3dcalc takes .1D files for use in combination with 3D datasets! \n"
1384 "\n"
1385 " N.B.: If you auto-transpose a .1D file on the command line, (by ending \n"
1386 " the filename with \\'), then 3dcalc will NOT treat it as the \n"
1387 " special case described above, but instead will treat it as \n"
1388 " a normal dataset, where each row in the transposed input is a \n"
1389 " 'voxel' time series. This would allow you to do differential \n"
1390 " subscripts on 1D time series, which program 1deval does not \n"
1391 " implement. For example: \n"
1392 "\n"
1393 " 3dcalc -a '1D: 3 4 5 6'\\' -b a+l -expr 'sqrt(a+b)' -prefix - \n"
1394 "\n"
1395 " This technique allows expression evaluation on multi-column \n"
1396 " .1D files, which 1deval also does not implement. For example: \n"
1397 "\n"
1398 " 3dcalc -a '1D: 3 4 5 | 1 2 3'\\' -expr 'cbrt(a)' -prefix - \n"
1399 "\n"
1400 "------------------------------------------------------------------------\n"
1401 "'1D:' INPUT: ~1~ \n"
1402 "----------- \n"
1403 "\n"
1404 " You can input a 1D time series 'dataset' directly on the command line, \n"
1405 " without an external file. The 'filename for such input takes the \n"
1406 " general format \n"
1407 "\n"
1408 " '1D:n_1@val_1,n_2@val_2,n_3@val_3,...' \n"
1409 "\n"
1410 " where each 'n_i' is an integer and each 'val_i' is a float. For \n"
1411 " example \n"
1412 "\n"
1413 " -a '1D:5@0,10@1,5@0,10@1,5@0' \n"
1414 "\n"
1415 " specifies that variable 'a' be assigned to a 1D time series of 35, \n"
1416 " alternating in blocks between values 0 and value 1. \n"
1417 "\n"
1418 " You can combine 3dUndump with 3dcalc to create an all zero 3D+time \n"
1419 " dataset from 'thin air', as in the commands \n"
1420 " 3dUndump -dimen 128 128 32 -prefix AllZero_A -datum float \n"
1421 " 3dcalc -a AllZero_A+orig -b '1D: 100@' -expr 0 -prefix AllZero_B \n"
1422 " If you replace the '0' expression with 'gran(0,1)', you'd get a \n"
1423 " random 3D+time dataset, which might be useful for testing purposes. \n"
1424 "\n"
1425 "------------------------------------------------------------------------\n"
1426 "'I:*.1D' and 'J:*.1D' and 'K:*.1D' INPUT: ~1~ \n"
1427 "---------------------------------------- \n"
1428 "\n"
1429 " You can input a 1D time series 'dataset' to be defined as spatially \n"
1430 " dependent instead of time dependent using a syntax like: \n"
1431 "\n"
1432 " -c I:fred.1D \n"
1433 "\n"
1434 " This indicates that the n-th value from file fred.1D is to be associated\n"
1435 " with the spatial voxel index i=n (respectively j=n and k=n for 'J: and \n"
1436 " K: input dataset names). This technique can be useful if you want to \n"
1437 " scale each slice by a fixed constant; for example: \n"
1438 "\n"
1439 " -a dset+orig -b K:slicefactor.1D -expr 'a*b' \n"
1440 "\n"
1441 " In this example, the '-b' value only varies in the k-index spatial \n"
1442 " direction. \n"
1443 "\n"
1444 "------------------------------------------------------------------------\n"
1445 "COORDINATES and PREDEFINED VALUES: ~1~ \n"
1446 "--------------------------------- \n"
1447 "\n"
1448 " If you don't use '-x', '-y', or '-z' for a dataset, then the voxel \n"
1449 " spatial coordinates will be loaded into those variables. For example, \n"
1450 " the expression 'a*step(x*x+y*y+z*z-100)' will zero out all the voxels \n"
1451 " inside a 10 mm radius of the origin x=y=z=0. \n"
1452 "\n"
1453 " Similarly, the '-t' value, if not otherwise used by a dataset or *.1D \n"
1454 " input, will be loaded with the voxel time coordinate, as determined \n"
1455 " from the header file created for the OUTPUT. Please note that the units\n"
1456 " of this are variable; they might be in milliseconds, seconds, or Hertz.\n"
1457 " In addition, slices of the dataset might be offset in time from one \n"
1458 " another, and this is allowed for in the computation of 't'. Use program\n"
1459 " 3dinfo to find out the structure of your datasets, if you are not sure.\n"
1460 " If no input datasets are 3D+time, then the effective value of TR is \n"
1461 " tstep in the output dataset, with t=0 at the first sub-brick. \n"
1462 "\n"
1463 " Similarly, the '-i', '-j', and '-k' values, if not otherwise used, \n"
1464 " will be loaded with the voxel spatial index coordinates. The '-l' \n"
1465 " (letter 'ell') value will be loaded with the temporal index coordinate.\n"
1466 "\n"
1467 " The '-n' value, if not otherwise used, will be loaded with the overall \n"
1468 " voxel 1D index. For a 3D dataset, n = i + j*NX + k*NX*NY, where \n"
1469 " NX, NY, NZ are the array dimensions of the 3D grid. [29 Jul 2010] \n"
1470 "\n"
1471 " Otherwise undefined letters will be set to zero. In the future, new \n"
1472 " default values for other letters may be added. \n"
1473 "\n"
1474 " NOTE WELL: By default, the coordinate order of (x,y,z) is the order in \n"
1475 " ********* which the data array is stored on disk; this order is output\n"
1476 " by 3dinfo. The options below control can change this order:\n"
1477 "\n"
1478 " -dicom }= Sets the coordinates to appear in DICOM standard (RAI) order,\n"
1479 " -RAI }= (the AFNI standard), so that -x=Right, -y=Anterior , -z=Inferior,\n"
1480 " +x=Left , +y=Posterior, +z=Superior.\n"
1481 "\n"
1482 " -SPM }= Sets the coordinates to appear in SPM (LPI) order, \n"
1483 " -LPI }= so that -x=Left , -y=Posterior, -z=Inferior,\n"
1484 " +x=Right, +y=Anterior , +z=Superior.\n"
1485 "\n"
1486 " The -LPI/-RAI behavior can also be achieved via the AFNI_ORIENT \n"
1487 " environment variable (27 Aug, 2014). \n"
1488 "------------------------------------------------------------------------\n"
1489 "DIFFERENTIAL SUBSCRIPTS [22 Nov 1999]: ~1~ \n"
1490 "----------------------- \n"
1491 "\n"
1492 " Normal calculations with 3dcalc are strictly on a per-voxel basis:\n"
1493 " there is no 'cross-talk' between spatial or temporal locations.\n"
1494 " The differential subscript feature allows you to specify variables\n"
1495 " that refer to different locations, relative to the base voxel.\n"
1496 " For example,\n"
1497 " -a fred+orig -b 'a[1,0,0,0]' -c 'a[0,-1,0,0]' -d 'a[0,0,2,0]'\n"
1498 " means: symbol 'a' refers to a voxel in dataset fred+orig,\n"
1499 " symbol 'b' refers to the following voxel in the x-direction,\n"
1500 " symbol 'c' refers to the previous voxel in the y-direction\n"
1501 " symbol 'd' refers to the 2nd following voxel in the z-direction\n"
1502 "\n"
1503 " To use this feature, you must define the base dataset (e.g., 'a')\n"
1504 " first. Then the differentially subscripted symbols are defined\n"
1505 " using the base dataset symbol followed by 4 integer subscripts,\n"
1506 " which are the shifts in the x-, y-, z-, and t- (or sub-brick index)\n"
1507 " directions. For example,\n"
1508 "\n"
1509 " -a fred+orig -b 'a[0,0,0,1]' -c 'a[0,0,0,-1]' -expr 'median(a,b,c)'\n"
1510 "\n"
1511 " will produce a temporal median smoothing of a 3D+time dataset (this\n"
1512 " can be done more efficiently with program 3dTsmooth).\n"
1513 "\n"
1514 " Note that the physical directions of the x-, y-, and z-axes depend\n"
1515 " on how the dataset was acquired or constructed. See the output of\n"
1516 " program 3dinfo to determine what direction corresponds to what axis.\n"
1517 "\n"
1518 " For convenience, the following abbreviations may be used in place of\n"
1519 " some common subscript combinations:\n"
1520 "\n"
1521 " [1,0,0,0] == +i [-1, 0, 0, 0] == -i\n"
1522 " [0,1,0,0] == +j [ 0,-1, 0, 0] == -j\n"
1523 " [0,0,1,0] == +k [ 0, 0,-1, 0] == -k\n"
1524 " [0,0,0,1] == +l [ 0, 0, 0,-1] == -l\n"
1525 "\n"
1526 " The median smoothing example can thus be abbreviated as\n"
1527 "\n"
1528 " -a fred+orig -b a+l -c a-l -expr 'median(a,b,c)'\n"
1529 "\n"
1530 " When a shift calls for a voxel that is outside of the dataset range,\n"
1531 " one of three things can happen:\n"
1532 "\n"
1533 " STOP => shifting stops at the edge of the dataset\n"
1534 " WRAP => shifting wraps back to the opposite edge of the dataset\n"
1535 " ZERO => the voxel value is returned as zero\n"
1536 "\n"
1537 " Which one applies depends on the setting of the shifting mode at the\n"
1538 " time the symbol using differential subscripting is defined. The mode\n"
1539 " is set by one of the switches '-dsSTOP', '-dsWRAP', or '-dsZERO'. The\n"
1540 " default mode is STOP. Suppose that a dataset has range 0..99 in the\n"
1541 " x-direction. Then when voxel 101 is called for, the value returned is\n"
1542 "\n"
1543 " STOP => value from voxel 99 [didn't shift past edge of dataset]\n"
1544 " WRAP => value from voxel 1 [wrapped back through opposite edge]\n"
1545 " ZERO => the number 0.0 \n"
1546 "\n"
1547 " You can set the shifting mode more than once - the most recent setting\n"
1548 " on the command line applies when a differential subscript symbol is\n"
1549 " encountered.\n"
1550 "\n"
1551 "N.B.: You can also use program 3dLocalstat to process data from a\n"
1552 " spatial neighborhood of each voxel; for example, to compute\n"
1553 " the maximum over a sphere of radius 9 mm placed around\n"
1554 " each voxel:\n"
1555 " 3dLocalstat -nbhd 'SPHERE(9)' -stat max -prefix Amax9 A+orig\n"
1556 "\n"
1557 "------------------------------------------------------------------------\n"
1558 "ISSUES: ~1~\n"
1559 "------ \n"
1560 "\n"
1561 " * Complex-valued datasets cannot be processed, except via '-cx2r'.\n"
1562 " * This program is not very efficient (but is faster than it once was).\n"
1563 " * Differential subscripts slow the program down even more.\n"
1564 "\n"
1565 "------------------------------------------------------------------------\n"
1566 ) ;
1567
1568 printf(
1569 "------------------------------------------------------------------------\n"
1570 "EXPRESSIONS: ~1~\n"
1571 "----------- \n"
1572 "\n"
1573 " As noted above, datasets are referred to by single letter variable names.\n"
1574 PARSER_HELP_STRING
1575 "\n"
1576 "** If you modify a statistical sub-brick, you may want to use program\n"
1577 " '3drefit' to modify the dataset statistical auxiliary parameters.\n"
1578 "\n"
1579 "** Computations are carried out in double precision before being\n"
1580 " truncated to the final output 'datum'.\n"
1581 "\n"
1582 "** Note that the quotes around the expression are needed so the shell\n"
1583 " doesn't try to expand * characters, or interpret parentheses.\n"
1584 "\n"
1585 "** Try the 'ccalc' program to see how the expression evaluator works.\n"
1586 " The arithmetic parser and evaluator is written in Fortran-77 and\n"
1587 " is derived from a program written long ago by RW Cox to facilitate\n"
1588 " compiling on an array processor hooked up to a VAX. (It's a mess, but\n"
1589 " it works - somewhat slowly - but hey, computers are fast these days.)\n"
1590 ) ;
1591
1592 PRINT_COMPILE_DATE ;
1593 }
1594
1595 /*------------------------------------------------------------------*/
1596
main(int argc,char * argv[])1597 int main( int argc , char *argv[] )
1598 {
1599 #define VSIZE 1024
1600
1601 double * atoz[26] ;
1602 int ii , ids , jj, kk, kt, ll, jbot, jtop ;
1603 THD_3dim_dataset * new_dset=NULL ;
1604 float ** buf;
1605 double temp[VSIZE];
1606 int nbad ; /* 09 Aug 2000: check for bad results */
1607
1608 THD_ivec3 iv ;
1609 THD_fvec3 fv ;
1610 float xxx[VSIZE], yyy[VSIZE], zzz[VSIZE] ;
1611 int iii,jjj,kkk , nx,ny,nz,nxy ;
1612 THD_dataxes * daxes ;
1613
1614 size_t tempnum , tempsiz ;
1615
1616 /*** read input options ***/
1617
1618 /*-- 20 Apr 2001: addto the arglist, if user wants to [RWCox] --*/
1619
1620 mainENTRY("3dcalc main"); machdep() ;
1621 PRINT_VERSION("3dcalc") ; AUTHOR("A cast of thousands") ;
1622 THD_check_AFNI_version("3dcalc") ;
1623
1624 if(argc == 1){ CALC_Syntax(); exit(0); } /* Bob's help shortcut */
1625
1626 set_obliquity_report(0); /* silence obliquity */
1627
1628 { int new_argc ; char ** new_argv ;
1629 addto_args( argc , argv , &new_argc , &new_argv ) ;
1630 if( new_argv != NULL ){ argc = new_argc ; argv = new_argv ; }
1631 }
1632
1633 AFNI_logger("3dcalc",argc,argv) ;
1634
1635 for (ii=0; ii<26; ii++) ntime[ii] = 0 ;
1636
1637 if( AFNI_yesenv("AFNI_FLOATIZE") ) CALC_datum = MRI_float ;
1638
1639 CALC_read_opts( argc , argv ) ;
1640
1641 if( argc < 2) ERROR_message("Too few options, use -help for details");
1642
1643 /*** make output dataset ***/
1644
1645 if( ntime_max == 1 || TS_make == 1 ){
1646 for( ids=0 ; ids < 26 ; ids++ ) if( CALC_dset[ids] != NULL ) break ;
1647 } else {
1648 for( ids=0 ; ids < 26 ; ids++ ) if( CALC_dset[ids] != NULL &&
1649 ntime[ids] > 1 ) break ;
1650 }
1651 if( ids == 26 ) ERROR_exit("Can't find template dataset?! [this should never happen]") ;
1652
1653 new_dset = EDIT_empty_copy( CALC_dset[ids] ) ;
1654
1655 /* 23 May 2005: check input datasets for axis consistency */
1656
1657 for( iii=0 ; iii < 26 ; iii++ ){
1658 double angle;
1659 if( iii != ids &&
1660 CALC_dset[iii] != NULL) {
1661 if (!EQUIV_DATAXES(new_dset->daxes,CALC_dset[iii]->daxes) )
1662 WARNING_message("dataset '%c'=%s grid mismatch with %s\n",
1663 'a'+iii , DSET_BRIKNAME(CALC_dset[iii]) ,
1664 DSET_BRIKNAME(CALC_dset[ids]) ) ;
1665
1666 }
1667 /* allow for small angle differences 22 May 2015 [rickr] */
1668 angle = dset_obliquity_angle_diff(new_dset, CALC_dset[iii],
1669 OBLIQ_ANGLE_THRESH);
1670 if (angle > 0.0) {
1671 WARNING_message(
1672 "dataset '%c'=%s has an obliquity difference of %f degrees with %s\n",
1673 'a'+iii , DSET_BRIKNAME(CALC_dset[iii]) ,
1674 angle, DSET_BRIKNAME(CALC_dset[ids]) );
1675 }
1676 }
1677
1678 /** make history for new dataset */
1679
1680 if( CALC_histpar < 0 ){
1681 for( iii=jjj=0 ; iii < 26 ; iii++ ) /* count number of input datasets */
1682 if( CALC_dset[iii] != NULL ) jjj++ ;
1683 } else {
1684 ids = CALC_histpar ;
1685 jjj = 1 ;
1686 }
1687
1688 if( jjj == 1 || AFNI_yesenv("AFNI_SIMPLE_HISTORY") ){
1689 tross_Copy_History( CALC_dset[ids] , new_dset ) ;
1690 } else { /* 27 Feb 2003 */
1691 char hbuf[64] ;
1692 tross_Append_History( new_dset ,
1693 "===================================" ) ;
1694 tross_Append_History( new_dset ,
1695 "=== History of inputs to 3dcalc ===" ) ;
1696 for( iii=0 ; iii < 26 ; iii++ ){
1697 if( CALC_dset[iii] != NULL ){
1698 sprintf(hbuf,"=== Input %c:", 'a'+iii ) ;
1699 tross_Append_History( new_dset , hbuf ) ;
1700 tross_Addto_History( CALC_dset[iii] , new_dset ) ;
1701 }
1702 }
1703 tross_Append_History( new_dset ,
1704 "===================================" ) ;
1705 }
1706 tross_Make_History( "3dcalc" , argc,argv , new_dset ) ;
1707
1708 if( CALC_datum < 0 ) CALC_datum = MRI_float ; /* 10 Feb 2003 */
1709
1710 EDIT_dset_items( new_dset ,
1711 ADN_prefix , CALC_output_prefix ,
1712 ADN_directory_name , CALC_session ,
1713 ADN_datum_all , CALC_datum ,
1714 ADN_none ) ;
1715
1716 if( DSET_NVALS(new_dset) != ntime_max )
1717 EDIT_dset_items( new_dset , ADN_nvals , ntime_max , ADN_none ) ;
1718
1719 /* 17 Apr 1998: if we are making up a 3D+time dataset,
1720 we need to attach some time axis info to it */
1721
1722 if( TS_make ){
1723 EDIT_dset_items( new_dset ,
1724 ADN_ntt , ntime_max ,
1725 ADN_ttdel , TS_dt ,
1726 ADN_ttorg , 0.0 ,
1727 ADN_ttdur , 0.0 ,
1728 ADN_tunits , UNITS_SEC_TYPE ,
1729 ADN_none ) ;
1730 }
1731
1732 if( ISFUNC(new_dset) && ! ISFUNCBUCKET(new_dset) && new_dset->taxis != NULL )
1733 EDIT_dset_items( new_dset , ADN_func_type , FUNC_FIM_TYPE , ADN_none ) ;
1734 else if( ISANATBUCKET(new_dset) ) /* 30 Nov 1997 */
1735 EDIT_dset_items( new_dset , ADN_func_type , ANAT_EPI_TYPE , ADN_none ) ;
1736
1737 if( THD_deathcon() && THD_is_file(new_dset->dblk->diskptr->header_name) )
1738 ERROR_exit("Output file %s already exists -- cannot continue!\n",
1739 new_dset->dblk->diskptr->header_name ) ;
1740
1741 for (ids=0; ids<26; ids++)
1742 atoz[ids] = (double *) malloc(sizeof(double) * VSIZE ) ;
1743
1744 for( ids=0 ; ids < 26 ; ids++ ) /* initialize to all zeros */
1745 for (ii=0; ii<VSIZE; ii++)
1746 atoz[ids][ii] = 0.0 ;
1747
1748 /*** loop over time steps ***/
1749
1750 nx = DSET_NX(new_dset) ;
1751 ny = DSET_NY(new_dset) ;
1752 nz = DSET_NZ(new_dset) ;
1753 nxy = nx * ny ; daxes = new_dset->daxes ;
1754
1755 buf = (float **) malloc(sizeof(float *) * ntime_max);
1756
1757 if( CALC_usetemp ){ /* 18 Oct 2005: -usetemp? */
1758 tempfnam = UNIQ_idcode() ;
1759 tempfnam[0] = 'C'; tempfnam[1] = 'A';
1760 tempfnam[2] = 'L'; tempfnam[3] = 'C'; tempfnam[4] = '_' ;
1761 tempfile = fopen( tempfnam , "w+b" ) ;
1762 INFO_message("Creating -usetemp file %s",tempfnam) ;
1763 atexit(calc_atexit) ;
1764 }
1765 tempsiz = ((size_t)CALC_nvox) * sizeof(float) ;
1766
1767 for( kt=0 ; kt < ntime_max ; kt++ ){
1768
1769 if( CALC_verbose )
1770 INFO_message("Computing sub-brick %d\n",kt) ;
1771
1772 /* 30 April 1998: only malloc output space as it is needed */
1773
1774 buf[kt] = (float *)calloc(1,tempsiz);
1775 if( buf[kt] == NULL )
1776 ERROR_exit("Can't malloc output dataset sub-brick %d -- out of memory :(",kt) ;
1777
1778 /*** loop over voxels ***/
1779
1780 for ( ii = 0 ; ii < CALC_nvox ; ii += VSIZE ) {
1781
1782 jbot = ii ;
1783 jtop = MIN( ii + VSIZE , CALC_nvox ) ;
1784
1785 /* load (x,y,z) coords of these voxels into arrays, if needed */
1786
1787 if( CALC_has_xyz ){ /* 17 May 2005 */
1788 for( jj=jbot ; jj < jtop ; jj++ ){
1789 LOAD_IVEC3( iv , jj%nx , (jj%nxy)/nx , jj/nxy ) ;
1790 fv = THD_3dind_to_3dmm( new_dset , iv ) ;
1791 if( CALC_mangle_xyz )
1792 fv = THD_3dmm_to_dicomm(new_dset,fv) ;
1793 UNLOAD_FVEC3(fv,xxx[jj-jbot],yyy[jj-jbot],zzz[jj-jbot]) ;
1794 if( CALC_mangle_xyz == MANGLE_LPI ){
1795 xxx[jj-jbot] = -xxx[jj-jbot] ; yyy[jj-jbot] = -yyy[jj-jbot] ;
1796 }
1797 }
1798 }
1799
1800 /* loop over datasets or other symbol definitions */
1801
1802 for (ids = 0 ; ids < 26 ; ids ++ ) { /* the whole alphabet */
1803
1804 /* 17 Apr 1998: if a time series is used here instead of a dataset,
1805 just copy the single value (or zero) to all voxels. */
1806
1807 if( TS_flim[ids] != NULL ){
1808 if( jbot == 0 ){ /* only must do this on first vector at each time */
1809 double tval ;
1810 if( kt < TS_flim[ids]->nx ) tval = TS_flar[ids][kt] ;
1811 else tval = 0.0 ;
1812
1813 for (jj =jbot ; jj < jtop ; jj ++ )
1814 atoz[ids][jj-ii] = tval ;
1815 }
1816 }
1817
1818 /** 22 Feb 2005: IJKAR 1D arrays **/
1819
1820 else if( IJKAR_flim[ids] != NULL ){
1821 int ss , ix=IJKAR_flim[ids]->nx ;
1822
1823 switch( IJKAR_dcod[ids] ){
1824 case 8:
1825 for( jj=jbot ; jj < jtop ; jj++ ){
1826 ss = (jj%nx) ;
1827 atoz[ids][jj-jbot] = (ss < ix) ? IJKAR_flar[ids][ss] : 0.0 ;
1828 }
1829 break ;
1830 case 9:
1831 for( jj=jbot ; jj < jtop ; jj++ ){
1832 ss = ((jj%nxy)/nx) ;
1833 atoz[ids][jj-jbot] = (ss < ix) ? IJKAR_flar[ids][ss] : 0.0 ;
1834 }
1835 break ;
1836 case 10:
1837 for( jj=jbot ; jj < jtop ; jj++ ){
1838 ss = (jj/nxy) ;
1839 atoz[ids][jj-jbot] = (ss < ix) ? IJKAR_flar[ids][ss] : 0.0 ;
1840 }
1841 break ;
1842 }
1843 }
1844
1845 /** 22 Nov 1999: if a differentially subscripted dataset is here **/
1846
1847 else if( CALC_dshift[ids] >= 0 ){
1848 int jds = CALC_dshift[ids] ; /* actual dataset index */
1849 int kts , jjs , ix,jy,kz ;
1850 int id=CALC_dshift_i[ids] , jd=CALC_dshift_j[ids] ,
1851 kd=CALC_dshift_k[ids] , ld=CALC_dshift_l[ids] ;
1852 int ijkd = ((id!=0) || (jd!=0) || (kd!=0)) ;
1853 int dsx = DSET_NX(CALC_dset[jds]) - 1 ;
1854 int dsy = DSET_NY(CALC_dset[jds]) - 1 ;
1855 int dsz = DSET_NZ(CALC_dset[jds]) - 1 ;
1856 int dst = ntime[jds] - 1 ;
1857 int mode = CALC_dshift_mode[ids] , dun=0 ;
1858
1859 kts = kt + ld ; /* t shift */
1860 if( kts < 0 || kts > dst ){
1861 switch( mode ){
1862 case DSHIFT_MODE_ZERO:
1863 for( jj=jbot ; jj < jtop ; jj++ ) atoz[ids][jj-ii] = 0.0 ;
1864 dun = 1 ;
1865 break ;
1866 default:
1867 case DSHIFT_MODE_STOP:
1868 if( kts < 0 ) kts = 0 ;
1869 else if( kts > dst ) kts = dst ;
1870 break ;
1871 case DSHIFT_MODE_WRAP:
1872 while( kts < 0 ) kts += (dst+1) ;
1873 while( kts > dst ) kts -= (dst+1) ;
1874 break ;
1875 }
1876 }
1877
1878 if( !dun ){ /* must get some actual data */
1879 for( dun=0,jj=jbot ; jj < jtop ; jj++ ){ /* loop over voxels */
1880 jjs = jj ; /* nominal voxel spatial index */
1881 if( ijkd ){ /* if spatial shift is ordered */
1882 ix = DSET_index_to_ix(CALC_dset[jds],jj) ;
1883 jy = DSET_index_to_jy(CALC_dset[jds],jj) ;
1884 kz = DSET_index_to_kz(CALC_dset[jds],jj) ;
1885
1886 ix += id ; /* x shift */
1887 if( ix < 0 || ix > dsx ){
1888 switch( mode ){
1889 case DSHIFT_MODE_ZERO:
1890 atoz[ids][jj-ii] = 0.0 ; dun = 1 ;
1891 break ;
1892 default:
1893 case DSHIFT_MODE_STOP:
1894 if( ix < 0 ) ix = 0 ;
1895 else if( ix > dsx ) ix = dsx ;
1896 break ;
1897 case DSHIFT_MODE_WRAP:
1898 while( ix < 0 ) ix += (dsx+1) ;
1899 while( ix > dsx ) ix -= (dsx+1) ;
1900 break ;
1901 }
1902 }
1903 if( dun ){ dun=0; continue; } /* go to next jj */
1904
1905 jy += jd ; /* y shift */
1906 if( jy < 0 || jy > dsy ){
1907 switch( mode ){
1908 case DSHIFT_MODE_ZERO:
1909 atoz[ids][jj-ii] = 0.0 ; dun = 1 ;
1910 break ;
1911 default:
1912 case DSHIFT_MODE_STOP:
1913 if( jy < 0 ) jy = 0 ;
1914 else if( jy > dsy ) jy = dsy ;
1915 break ;
1916 case DSHIFT_MODE_WRAP:
1917 while( jy < 0 ) jy += (dsy+1) ;
1918 while( jy > dsy ) jy -= (dsy+1) ;
1919 break ;
1920 }
1921 }
1922 if( dun ){ dun=0; continue; } /* go to next jj */
1923
1924 kz += kd ; /* z shift */
1925 if( kz < 0 || kz > dsz ){
1926 switch( mode ){
1927 case DSHIFT_MODE_ZERO:
1928 atoz[ids][jj-ii] = 0.0 ; dun = 1 ;
1929 break ;
1930 default:
1931 case DSHIFT_MODE_STOP:
1932 if( kz < 0 ) kz = 0 ;
1933 else if( kz > dsz ) kz = dsz ;
1934 break ;
1935 case DSHIFT_MODE_WRAP:
1936 while( kz < 0 ) kz += (dsz+1) ;
1937 while( kz > dsz ) kz -= (dsz+1) ;
1938 break ;
1939 }
1940 }
1941 if( dun ){ dun=0; continue; } /* go to next jj */
1942
1943 jjs = DSET_ixyz_to_index(CALC_dset[jds],ix,jy,kz) ;
1944 } /* end of spatial shift index calculation */
1945
1946 switch( CALC_type[jds] ) { /* extract data */
1947 case MRI_short:
1948 atoz[ids][jj-ii] = CALC_short[jds][kts][jjs]
1949 * CALC_ffac[jds][kts];
1950 break ;
1951 case MRI_float:
1952 atoz[ids][jj-ii] = CALC_float[jds][kts][jjs]
1953 * CALC_ffac[jds][kts];
1954 break ;
1955 case MRI_byte:
1956 atoz[ids][jj-ii] = CALC_byte[jds][kts][jjs]
1957 * CALC_ffac[jds][kts];
1958 break ;
1959 case MRI_rgb:
1960 atoz[ids][jj-ii] = Rfac*CALC_byte[jds][kts][3*jjs ]
1961 +Gfac*CALC_byte[jds][kts][3*jjs+1]
1962 +Bfac*CALC_byte[jds][kts][3*jjs+2] ;
1963 break ;
1964 case MRI_complex:{ /* 10 Mar 2006 */
1965 complex cv=CALC_complex[jds][kts][jjs] ;
1966 float xx=cv.r, yy=cv.i , vv ;
1967 switch( CALC_cxcode[ids] ){ /* ids, NOT jds! */
1968 case CX_REALPART: vv = xx ; break ;
1969 case CX_IMAGPART: vv = yy ; break ;
1970 case CX_PHASE: vv = (xx==0.0f && yy==0.0f)
1971 ? 0.0f : atan2(yy,xx) ; break ;
1972 default:
1973 case CX_MAGNITUDE: vv = sqrt(xx*xx+yy*yy) ; break ;
1974 }
1975 atoz[ids][jj-ii] = vv ;
1976 }
1977 } /* end of data extraction switch */
1978 } /* end of loop over voxels */
1979 } /* end of getting actual data */
1980 } /* end of differential subscripted input */
1981
1982 /** the case of a 3D dataset (i.e., only 1 sub-brick) **/
1983
1984 else if ( ntime[ids] == 1 && CALC_type[ids] >= 0 ) {
1985 switch( CALC_type[ids] ) {
1986 case MRI_short:
1987 if( CALC_noffac[ids] ) /* 14 Nov 2003 */
1988 for (jj =jbot ; jj < jtop ; jj ++ )
1989 atoz[ids][jj-ii] = CALC_short[ids][0][jj] ;
1990 else
1991 for (jj =jbot ; jj < jtop ; jj ++ )
1992 atoz[ids][jj-ii] = CALC_short[ids][0][jj] * CALC_ffac[ids][0] ;
1993 break;
1994
1995 case MRI_float:
1996 if( CALC_noffac[ids] ) /* 14 Nov 2003 */
1997 for (jj =jbot ; jj < jtop ; jj ++ )
1998 atoz[ids][jj-ii] = CALC_float[ids][0][jj] ;
1999 else
2000 for (jj =jbot ; jj < jtop ; jj ++ )
2001 atoz[ids][jj-ii] = CALC_float[ids][0][jj] * CALC_ffac[ids][0] ;
2002 break;
2003
2004 case MRI_byte:
2005 if( CALC_noffac[ids] ) /* 14 Nov 2003 */
2006 for (jj =jbot ; jj < jtop ; jj ++ )
2007 atoz[ids][jj-ii] = CALC_byte[ids][0][jj] ;
2008 else
2009 for (jj =jbot ; jj < jtop ; jj ++ )
2010 atoz[ids][jj-ii] = CALC_byte[ids][0][jj] * CALC_ffac[ids][0] ;
2011 break;
2012
2013 case MRI_rgb:
2014 for (jj =jbot ; jj < jtop ; jj ++ )
2015 atoz[ids][jj-ii] = Rfac*CALC_byte[ids][0][3*jj ]
2016 +Gfac*CALC_byte[ids][0][3*jj+1]
2017 +Bfac*CALC_byte[ids][0][3*jj+2] ;
2018 break;
2019
2020 case MRI_complex:{ /* 10 Mar 2006 */
2021 complex cv ; float xx,yy,vv ;
2022 for( jj=jbot ; jj < jtop ; jj++ ){
2023 cv=CALC_complex[ids][0][jj] ; xx = cv.r ; yy = cv.i ;
2024 switch( CALC_cxcode[ids] ){
2025 case CX_REALPART: vv = xx ; break ;
2026 case CX_IMAGPART: vv = yy ; break ;
2027 case CX_PHASE: vv = (xx==0.0f && yy==0.0f)
2028 ? 0.0f : atan2(yy,xx) ; break ;
2029 default:
2030 case CX_MAGNITUDE: vv = sqrt(xx*xx+yy*yy) ; break ;
2031 }
2032 atoz[ids][jj-ii] = vv ;
2033 }
2034 }
2035 }
2036 } /** end of 3D dataset **/
2037
2038 /** the case of a 3D+time dataset (or a bucket, etc.) **/
2039
2040 else if( ntime[ids] > 1 && CALC_type[ids] >= 0 ) {
2041 switch ( CALC_type[ids] ) {
2042 case MRI_short:
2043 if( CALC_noffac[ids] )
2044 for (jj = jbot ; jj < jtop ; jj ++ )
2045 atoz[ids][jj-ii] = CALC_short[ids][kt][jj] ;
2046 else
2047 for (jj = jbot ; jj < jtop ; jj ++ )
2048 atoz[ids][jj-ii] = CALC_short[ids][kt][jj] * CALC_ffac[ids][kt];
2049 break;
2050
2051 case MRI_float:
2052 if( CALC_noffac[ids] )
2053 for (jj = jbot ; jj < jtop ; jj ++ )
2054 atoz[ids][jj-ii] = CALC_float[ids][kt][jj] ;
2055 else
2056 for (jj = jbot ; jj < jtop ; jj ++ )
2057 atoz[ids][jj-ii] = CALC_float[ids][kt][jj] * CALC_ffac[ids][kt];
2058 break;
2059
2060 case MRI_byte:
2061 if( CALC_noffac[ids] )
2062 for (jj = jbot ; jj < jtop ; jj ++ )
2063 atoz[ids][jj-ii] = CALC_byte[ids][kt][jj] ;
2064 else
2065 for (jj = jbot ; jj < jtop ; jj ++ )
2066 atoz[ids][jj-ii] = CALC_byte[ids][kt][jj] * CALC_ffac[ids][kt];
2067 break;
2068
2069 case MRI_rgb:
2070 for (jj =jbot ; jj < jtop ; jj ++ )
2071 atoz[ids][jj-ii] = Rfac*CALC_byte[ids][kt][3*jj ]
2072 +Gfac*CALC_byte[ids][kt][3*jj+1]
2073 +Bfac*CALC_byte[ids][kt][3*jj+2] ;
2074 break;
2075
2076 case MRI_complex:{ /* 10 Mar 2006 */
2077 complex cv ; float xx,yy,vv ;
2078 for( jj=jbot ; jj < jtop ; jj++ ){
2079 cv=CALC_complex[ids][kt][jj] ; xx = cv.r ; yy = cv.i ;
2080 switch( CALC_cxcode[ids] ){
2081 case CX_REALPART: vv = xx ; break ;
2082 case CX_IMAGPART: vv = yy ; break ;
2083 case CX_PHASE: vv = (xx==0.0f && yy==0.0f)
2084 ? 0.0f : atan2(yy,xx) ; break ;
2085 default:
2086 case CX_MAGNITUDE: vv = sqrt(xx*xx+yy*yy) ; break ;
2087 }
2088 atoz[ids][jj-ii] = vv ;
2089 }
2090 }
2091 }
2092 }
2093
2094 /* the case of a voxel (x,y,z) or (i,j,k) coordinate */
2095
2096 else if( CALC_has_predefined ) {
2097
2098 switch( ids ){
2099 case 23: /* x */
2100 if( HAS_X )
2101 for( jj=jbot ; jj < jtop ; jj++ )
2102 atoz[ids][jj-ii] = xxx[jj-ii] ;
2103 break ;
2104
2105 case 24: /* y */
2106 if( HAS_Y )
2107 for( jj=jbot ; jj < jtop ; jj++ )
2108 atoz[ids][jj-ii] = yyy[jj-ii] ;
2109 break ;
2110
2111 case 25: /* z */
2112 if( HAS_Z )
2113 for( jj=jbot ; jj < jtop ; jj++ )
2114 atoz[ids][jj-ii] = zzz[jj-ii] ;
2115 break ;
2116
2117 case 8: /* i */
2118 if( HAS_I )
2119 for( jj=jbot ; jj < jtop ; jj++ )
2120 atoz[ids][jj-ii] = (jj%nx) ;
2121 break ;
2122
2123 case 9: /* j */
2124 if( HAS_J )
2125 for( jj=jbot ; jj < jtop ; jj++ )
2126 atoz[ids][jj-ii] = ((jj%nxy)/nx) ;
2127 break ;
2128
2129 case 10: /* k */
2130 if( HAS_K )
2131 for( jj=jbot ; jj < jtop ; jj++ )
2132 atoz[ids][jj-ii] = (jj/nxy) ;
2133 break ;
2134
2135 case 19: /* t */
2136 if( HAS_T )
2137 for( jj=jbot ; jj < jtop ; jj++ )
2138 atoz[ids][jj-ii] = THD_timeof_vox(kt,jj,new_dset) ;
2139 break ;
2140
2141 case 11: /* l */
2142 if( HAS_L )
2143 for( jj=jbot ; jj < jtop ; jj++ )
2144 atoz[ids][jj-ii] = kt ;
2145 break ;
2146
2147 case 13: /* n */
2148 if( HAS_N )
2149 for( jj=jbot ; jj < jtop ; jj++ )
2150 atoz[ids][jj-ii] = jj ;
2151 break ;
2152
2153 } /* end of switch on symbol subscript */
2154
2155 } /* end of choice over data type (if-else cascade) */
2156 } /* end of loop over datasets/symbols */
2157
2158 /**** actually do the calculation work! ****/
2159
2160 PARSER_evaluate_vector(CALC_code, atoz, jtop-jbot, temp);
2161 for ( jj = jbot ; jj < jtop ; jj ++ )
2162 buf[kt][jj] = temp[jj-ii];
2163
2164 } /*---------- end of loop over space (voxels) ----------*/
2165
2166 /* 09 Aug 2000: check results for validity */
2167
2168 nbad = thd_floatscan( CALC_nvox , buf[kt] ) ;
2169 if( nbad > 0 )
2170 WARNING_message("%d bad floats replaced by 0 in sub-brick %d\n\a",
2171 nbad , kt ) ;
2172
2173 /* 30 April 1998: purge 3D+time sub-bricks if possible */
2174
2175 if( ! CALC_has_timeshift ){
2176 for( ids=0 ; ids < 26 ; ids ++ ){
2177 if( CALC_dset[ids] != NULL && ntime[ids] > 1 &&
2178 CALC_dset[ids]->dblk->malloc_type == DATABLOCK_MEM_MALLOC ){
2179
2180 void * ptr = DSET_ARRAY(CALC_dset[ids],kt) ;
2181 if( ptr != NULL ) free(ptr) ;
2182 mri_clear_data_pointer( DSET_BRICK(CALC_dset[ids],kt) ) ;
2183 }
2184 }
2185 }
2186
2187 /* 27 Nov 2019: remove isolas? */
2188
2189 if( CALC_do_isolas ){
2190 int ntot = remove_isolated_stuff( nx,ny,nz , buf[kt] , CALC_do_isolas ) ;
2191 if( ntot > 0 )
2192 ININFO_message("removed %d isolas from volume %d",ntot,kt) ;
2193 }
2194
2195 /* 18 Oct 2005: write to a temp file? */
2196
2197 if( tempfile != NULL ){
2198 tempnum = fwrite( buf[kt] , 1 , tempsiz , tempfile ) ;
2199 free( buf[kt] ) ; buf[kt] = NULL ;
2200 if( tempnum < tempsiz ){
2201 ERROR_message("-usetemp #%d: only %u bytes written, out of %u\n",
2202 kt , (unsigned)tempnum , (unsigned)tempsiz ) ;
2203 perror("** Unix error message") ;
2204 } else if( CALC_verbose )
2205 INFO_message("-usetemp #%d output %u bytes",kt,(unsigned)tempsiz) ;
2206 }
2207
2208 } /*-------------- end of loop over time steps -------------*/
2209
2210 for( ids=0 ; ids < 26 ; ids++ ){
2211 if( CALC_dset[ids] != NULL ) PURGE_DSET( CALC_dset[ids] ) ;
2212 if( TS_flim[ids] != NULL ) mri_free( TS_flim[ids] ) ;
2213 if( IJKAR_flim[ids]!= NULL ) mri_free( IJKAR_flim[ids] ) ;
2214 }
2215
2216 /*** attach new data to output brick ***/
2217
2218 if( tempfile != NULL ) rewind(tempfile) ; /* 18 Oct 2005 */
2219
2220 #undef TGET
2221 #define TGET(q) \
2222 if( tempfile != NULL ){ \
2223 buf[q] = (float *)calloc(1,tempsiz) ; \
2224 tempnum = fread( buf[q] , 1 , tempsiz , tempfile ) ; \
2225 if( tempnum < tempsiz ){ \
2226 ERROR_message("-usetemp #%d: only %u bytes read, out of %u\n", \
2227 (q) , (unsigned)tempnum , (unsigned)tempsiz ) ; \
2228 perror("** Unix error message") ; \
2229 } \
2230 }
2231
2232 switch( CALC_datum ){
2233
2234 default:
2235 ERROR_exit("Somehow ended up with CALC_datum = %d [this should never happen]",CALC_datum) ;
2236 exit(1) ;
2237
2238 /* the easy case! */
2239
2240 case MRI_float:{
2241 for( ii=0 ; ii < ntime_max ; ii++ ){
2242 TGET(ii) ; /* 18 Oct 2005: load from temp file? */
2243
2244 #ifdef ALLOW_SORT
2245 if( CALC_sort > 0 ) qsort_float ( CALC_nvox, buf[ii] ) ;
2246 else if( CALC_sort < 0 ) qsort_float_rev( CALC_nvox, buf[ii] ) ;
2247 #endif
2248
2249 #ifdef ALLOW_FDR /* only experimental, not useful -- cf. 3dFDR instead */
2250 if( CALC_fdrize && DSET_BRICK_STATCODE(new_dset,ii) > 0 ){
2251 MRI_IMAGE *qim=mri_new_vol_empty(CALC_nvox,1,1,MRI_float) ;
2252 mri_fix_data_pointer( buf[ii] , qim ) ;
2253 mri_fdrize(qim,DSET_BRICK_STATCODE(new_dset,ii),
2254 DSET_BRICK_STATAUX (new_dset,ii) , CALC_fdrize==2 ) ;
2255 mri_clear_data_pointer(qim); mri_free(qim);
2256 if( !ISFUNCBUCKET(new_dset) )
2257 EDIT_dset_items( new_dset ,
2258 ADN_type , HEAD_FUNC_TYPE,
2259 ADN_func_type, FUNC_BUCK_TYPE, ADN_none ) ;
2260 EDIT_BRICK_TO_FIZT(new_dset,ii) ;
2261 }
2262 #endif
2263
2264 EDIT_substitute_brick(new_dset, ii, MRI_float, buf[ii]);
2265 DSET_BRICK_FACTOR(new_dset, ii) = 0.0;
2266 }
2267 }
2268 break ;
2269
2270 /* the harder cases */
2271
2272 case MRI_byte: /* modified 31 Mar 1999 to scale each sub-brick */
2273 case MRI_short:{ /* with its own factor, rather than use the same */
2274 void **dfim ; /* factor for each sub-brick -- RWCox */
2275 float gtop=0.0f , fimfac , gtemp ;
2276
2277 if( CALC_verbose )
2278 INFO_message("Scaling output to type %s brick(s)\n",
2279 MRI_TYPE_name[CALC_datum] ) ;
2280
2281 dfim = (void **) malloc( sizeof(void *) * ntime_max ) ;
2282
2283 if( CALC_gscale ){ /* 01 Apr 1999: global scaling */
2284 gtop = 0.0 ;
2285 for( ii=0 ; ii < ntime_max ; ii++ ){
2286 gtemp = MCW_vol_amax( CALC_nvox , 1 , 1 , MRI_float, buf[ii] ) ;
2287 gtop = MAX( gtop , gtemp ) ;
2288 if( gtemp == 0.0 )
2289 WARNING_message("output sub-brick %d is all zeros!\n",ii) ;
2290 }
2291 }
2292
2293 for( ii=0 ; ii < ntime_max ; ii++ ) {
2294
2295 TGET(ii) ; /* 18 Oct 2005: temp load */
2296
2297 #ifdef ALLOW_SORT
2298 if( CALC_sort > 0 ) qsort_float ( CALC_nvox, buf[ii] ) ;
2299 else if( CALC_sort < 0 ) qsort_float_rev( CALC_nvox, buf[ii] ) ;
2300 #endif
2301
2302 /* get max of this sub-brick, if not doing global scaling */
2303
2304 if( ! CALC_gscale ){
2305 gtop = MCW_vol_amax( CALC_nvox , 1 , 1 , MRI_float, buf[ii] ) ;
2306 if( gtop == 0.0 )
2307 WARNING_message("output sub-brick %d is all zeros!\n",ii) ;
2308 }
2309
2310 /* compute scaling factor for this brick into fimfac */
2311
2312 if( CALC_fscale ){ /* 16 Mar 1998: forcibly scale */
2313 fimfac = (gtop > 0.0) ? MRI_TYPE_maxval[CALC_datum] / gtop : 0.0 ;
2314
2315 } else if( !CALC_nscale ){ /* maybe scale */
2316
2317 /* (gtop <= 1.0) to (gtop < 1.0) 23 Mar 2006 [rickr/dglen] */
2318 fimfac = (gtop > MRI_TYPE_maxval[CALC_datum]
2319 || (gtop > 0.0 && gtop < 1.0) )
2320 ? MRI_TYPE_maxval[CALC_datum] / gtop : 0.0 ;
2321
2322 if( fimfac == 0.0 && gtop > 0.0 ){ /* 28 Jul 2003: check for non-integers */
2323 float fv,iv ; int kk ;
2324 for( kk=0 ; kk < CALC_nvox ; kk++ ){
2325 fv = buf[ii][kk] ; iv = rint(fv) ;
2326 if( fabsf(fv-iv) >= 0.01 ){
2327 fimfac = MRI_TYPE_maxval[CALC_datum]/ gtop ; break ;
2328 }
2329 }
2330 }
2331
2332 } else { /* user says "don't scale" */
2333 fimfac = 0.0 ;
2334 }
2335
2336 if( CALC_verbose ){
2337 if( fimfac != 0.0 )
2338 INFO_message("Sub-brick %d scale factor = %f\n",ii,fimfac) ;
2339 else
2340 INFO_message("Sub-brick %d: no scale factor\n" ,ii) ;
2341 }
2342
2343 /* make space for output brick and scale into it */
2344
2345 dfim[ii] = (void *) malloc( mri_datum_size(CALC_datum) * CALC_nvox ) ;
2346 if( dfim[ii] == NULL ) ERROR_exit("malloc fails at output[%d] -- out of memory :(",ii);
2347
2348 if( CALC_datum == MRI_byte ){ /* 29 Nov 2004: check for bad byte-ization */
2349 int nneg ;
2350 for( nneg=jj=0 ; jj < CALC_nvox ; jj++ ) nneg += (buf[ii][jj] < 0.0f) ;
2351 if( nneg > 0 )
2352 WARNING_message(
2353 "sub-brick #%d has %d negative values set=0 in conversion to bytes\n",
2354 ii , nneg ) ;
2355 }
2356
2357 EDIT_coerce_scale_type( CALC_nvox , fimfac ,
2358 MRI_float, buf[ii] , CALC_datum,dfim[ii] ) ;
2359
2360 if( CALC_datum == MRI_short )
2361 EDIT_misfit_report( DSET_FILECODE(new_dset) , ii ,
2362 CALC_nvox , (fimfac != 0.0f) ? 1.0f/fimfac : 0.0f ,
2363 dfim[ii] , buf[ii] ) ;
2364
2365 free( buf[ii] ) ; buf[ii] = NULL ;
2366
2367 /* put result into output dataset */
2368
2369 EDIT_substitute_brick(new_dset, ii, CALC_datum, dfim[ii] );
2370
2371 DSET_BRICK_FACTOR(new_dset,ii) = (fimfac != 0.0) ? 1.0/fimfac : 0.0 ;
2372 }
2373 }
2374 break ;
2375 }
2376
2377 if( tempfile != NULL ){ /* 18 Oct 2005 */
2378 INFO_message("Deleting -usetemp file %s",tempfnam) ;
2379 fclose(tempfile) ; tempfile = NULL ;
2380 remove(tempfnam) ; tempfnam = NULL ;
2381 }
2382
2383 if( CALC_verbose ) INFO_message("Computing output statistics\n") ;
2384 THD_load_statistics( new_dset ) ;
2385
2386 DSET_BRICK_FDRCURVE_ALLKILL(new_dset) ; /* 24 Jan 2008 */
2387
2388 THD_write_3dim_dataset( NULL,NULL , new_dset , True ) ;
2389 WROTE_DSET(new_dset) ;
2390
2391 exit(0) ;
2392 }
2393
2394 /*----------------------------------------------------------------------------*/
2395 /* Remove isolated stuff in-place from a 2D or 3D image. [Nov 2019 - RWCox] */
2396 /*----------------------------------------------------------------------------*/
2397
2398 #define FF(i,j,k) far[(i)+(j)*nx+(k)*nxy]
2399 #define QQ(i,j,k) qar[(i)+(j)*nx+(k)*nxy]
2400
remove_isolated_stuff(int nnx,int nny,int nnz,float * far,int maxite)2401 int remove_isolated_stuff( int nnx, int nny, int nnz, float *far, int maxite )
2402 {
2403 int nx,ny,nz ,nxy,nxyz , ii,jj,kk,ll,cc , nblast,nite,ntot=0 ;
2404 float nb[27] ; /* temp data */
2405 float *qar ;
2406
2407 ENTRY("remove_isolated_stuff") ;
2408
2409 if( maxite < 1 || far == NULL ) RETURN(0) ;
2410
2411 nx = nnx ;
2412 ny = nny ;
2413 nz = nnz ;
2414 nxy = nx*ny ; nxyz = nx*ny*nz ;
2415
2416 ii = (nx < 5) + (ny < 5) + (nz < 5) ; /* how many small dimensions */
2417 if( ii > 1 ) RETURN(0) ;
2418
2419 qar = (float *)malloc(sizeof(float)*nxyz) ;
2420
2421 nblast = 666 ;
2422 for( nite=0 ; nite < maxite && nblast > 0 ; nite++ ){
2423 nblast = 0 ;
2424 memcpy( qar , far , sizeof(float)*nxyz ) ;
2425
2426 /* edit 3x3x3 volumes */
2427
2428 for( kk=1 ; kk < nz-1 ; kk++ ){
2429 for( jj=1 ; jj < ny-1 ; jj++ ){
2430 for( ii=1 ; ii < nx-1 ; ii++ ){
2431 if( FF(ii,jj,kk) != 0.0f ){ /* must have at least 3 nonzero */
2432 nb[ 0] = FF(ii-1,jj-1,kk-1) ; /* voxels in a 3x3x3 neighborhood */
2433 nb[ 1] = FF(ii ,jj-1,kk-1) ;
2434 nb[ 2] = FF(ii+1,jj-1,kk-1) ;
2435 nb[ 3] = FF(ii-1,jj ,kk-1) ;
2436 nb[ 4] = FF(ii ,jj ,kk-1) ;
2437 nb[ 5] = FF(ii+1,jj ,kk-1) ;
2438 nb[ 6] = FF(ii-1,jj+1,kk-1) ;
2439 nb[ 7] = FF(ii ,jj+1,kk-1) ;
2440 nb[ 8] = FF(ii+1,jj+1,kk-1) ;
2441 nb[ 9] = FF(ii-1,jj-1,kk ) ;
2442 nb[10] = FF(ii ,jj-1,kk ) ;
2443 nb[11] = FF(ii+1,jj-1,kk ) ;
2444 nb[12] = FF(ii-1,jj ,kk ) ;
2445 nb[13] = FF(ii ,jj ,kk ) ;
2446 nb[14] = FF(ii+1,jj ,kk ) ;
2447 nb[15] = FF(ii-1,jj-1,kk ) ;
2448 nb[16] = FF(ii ,jj-1,kk ) ;
2449 nb[17] = FF(ii+1,jj-1,kk ) ;
2450 nb[18] = FF(ii-1,jj+1,kk+1) ;
2451 nb[19] = FF(ii ,jj+1,kk+1) ;
2452 nb[20] = FF(ii+1,jj+1,kk+1) ;
2453 nb[21] = FF(ii-1,jj ,kk+1) ;
2454 nb[22] = FF(ii ,jj ,kk+1) ;
2455 nb[23] = FF(ii+1,jj ,kk+1) ;
2456 nb[24] = FF(ii-1,jj+1,kk+1) ;
2457 nb[25] = FF(ii ,jj+1,kk+1) ;
2458 nb[26] = FF(ii+1,jj+1,kk+1) ;
2459
2460 for( ll=cc=0 ; ll < 27 ; ll++ ) if( nb[ll] != 0.0f ) cc++ ;
2461 if( cc < 4 ){ QQ(ii,jj,kk) = 0.0f ; nblast++ ; }
2462 }
2463 } } }
2464
2465 /* edit 3x3 areas in each 2D orientation */
2466
2467 for( kk=0 ; kk < nz ; kk++ ){
2468 for( jj=1 ; jj < ny-1 ; jj++ ){
2469 for( ii=1 ; ii < nx-1 ; ii++ ){
2470 if( QQ(ii,jj,kk) != 0.0f ){ /* must have at least 2 nonzero */
2471 nb[ 0] = FF(ii-1,jj-1,kk) ; /* voxels in a 3x3 neighborhood */
2472 nb[ 1] = FF(ii ,jj-1,kk) ;
2473 nb[ 2] = FF(ii+1,jj-1,kk) ;
2474 nb[ 3] = FF(ii-1,jj ,kk) ;
2475 nb[ 4] = FF(ii ,jj ,kk) ;
2476 nb[ 5] = FF(ii+1,jj ,kk) ;
2477 nb[ 6] = FF(ii-1,jj+1,kk) ;
2478 nb[ 7] = FF(ii ,jj+1,kk) ;
2479 nb[ 8] = FF(ii+1,jj+1,kk) ;
2480 for( ll=cc=0 ; ll < 9 ; ll++ ) if( nb[ll] != 0.0f ) cc++ ;
2481 if( cc < 2 ){ QQ(ii,jj,kk) = 0.0f ; nblast++ ; }
2482 }
2483 } } }
2484
2485 for( jj=0 ; jj < ny ; jj++ ){
2486 for( kk=1 ; kk < nz-1 ; kk++ ){
2487 for( ii=1 ; ii < nx-1 ; ii++ ){
2488 if( QQ(ii,jj,kk) != 0.0f ){ /* must have at least 2 nonzero */
2489 nb[ 0] = FF(ii-1,jj,kk-1) ; /* voxels in a 3x3 neighborhood */
2490 nb[ 1] = FF(ii ,jj,kk-1) ;
2491 nb[ 2] = FF(ii+1,jj,kk-1) ;
2492 nb[ 3] = FF(ii-1,jj,kk ) ;
2493 nb[ 4] = FF(ii ,jj,kk ) ;
2494 nb[ 5] = FF(ii+1,jj,kk ) ;
2495 nb[ 6] = FF(ii-1,jj,kk+1) ;
2496 nb[ 7] = FF(ii ,jj,kk+1) ;
2497 nb[ 8] = FF(ii+1,jj,kk+1) ;
2498 for( ll=cc=0 ; ll < 9 ; ll++ ) if( nb[ll] != 0.0f ) cc++ ;
2499 if( cc < 2 ){ QQ(ii,jj,kk) = 0.0f ; nblast++ ; }
2500 }
2501 } } }
2502
2503 for( ii=0 ; ii < nx ; ii++ ){
2504 for( jj=1 ; jj < ny-1 ; jj++ ){
2505 for( kk=1 ; kk < nz-1 ; kk++ ){
2506 if( QQ(ii,jj,kk) != 0.0f ){ /* must have at least 2 nonzero */
2507 nb[ 0] = FF(ii,jj-1,kk-1) ; /* voxels in a 3x3 neighborhood */
2508 nb[ 1] = FF(ii,jj ,kk-1) ;
2509 nb[ 2] = FF(ii,jj+1,kk-1) ;
2510 nb[ 3] = FF(ii,jj-1,kk ) ;
2511 nb[ 4] = FF(ii,jj ,kk ) ;
2512 nb[ 5] = FF(ii,jj+1,kk ) ;
2513 nb[ 6] = FF(ii,jj-1,kk+1) ;
2514 nb[ 7] = FF(ii,jj ,kk+1) ;
2515 nb[ 8] = FF(ii,jj+1,kk+1) ;
2516 for( ll=cc=0 ; ll < 9 ; ll++ ) if( nb[ll] != 0.0f ) cc++ ;
2517 if( cc < 2 ){ QQ(ii,jj,kk) = 0.0f ; nblast++ ; }
2518 }
2519 } } }
2520
2521 if( nblast > 0 ) memcpy( far , qar , sizeof(float)*nxyz) ;
2522 ntot += nblast ;
2523 }
2524
2525 free(qar) ;
2526 RETURN(ntot) ;
2527 }
2528
2529 #undef FF
2530 #undef QQ
2531