1 /* ----------------------------- MNI Header -----------------------------------
2 @NAME : mincmath
3 @INPUT : argc, argv - command line arguments
4 @OUTPUT : (none)
5 @RETURNS : status
6 @DESCRIPTION: Program to do simple math on minc files
7 @METHOD :
8 @GLOBALS :
9 @CALLS :
10 @CREATED : April 28, 1995 (Peter Neelin)
11 @MODIFIED :
12 * $Log: mincmath.c,v $
13 * Revision 6.15 2009-01-20 11:58:13 rotor
14 * * CMakeLists.txt: updated version
15 * * Updated Changelog to include releases
16 * * Warning cleanups below
17 * * conversion/dcm2mnc/minc_file.c: fixed printf type
18 * * conversion/dcm2mnc/siemens_to_dicom.c: fixed printf type
19 * * conversion/ecattominc/machine_indep.c: added string.h and fixed
20 * 2 fprintf missing format args
21 * * conversion/micropet/upet2mnc.c: fixed two fprintf format args
22 * * conversion/minctoecat/ecat_write.c: added string.h
23 * * conversion/minctoecat/minctoecat.c: added missing argument to fprintf
24 * * conversion/nifti1/mnc2nii.c: fixed incorrect printf type
25 * * progs/mincview/invert_raw_image.c: added fwrite checking
26 *
27 * Revision 6.14 2008/01/17 02:33:02 rotor
28 * * removed all rcsids
29 * * removed a bunch of ^L's that somehow crept in
30 * * removed old (and outdated) BUGS file
31 *
32 * Revision 6.13 2008/01/13 04:30:28 stever
33 * Add braces around static initializers.
34 *
35 * Revision 6.12 2008/01/12 19:08:15 stever
36 * Add __attribute__ ((unused)) to all rcsid variables.
37 *
38 * Revision 6.11 2007/12/11 12:43:01 rotor
39 * * added static to all global variables in main programs to avoid linking
40 * problems with libraries (compress in mincconvert and libz for example)
41 *
42 * Revision 6.10 2005/08/26 21:07:17 bert
43 * Use #if rather than #ifdef with MINC2 symbol, and be sure to include config.h whereever MINC2 is used
44 *
45 * Revision 6.9 2004/12/14 23:40:07 bert
46 * Get rid of c99 compilation problems
47 *
48 * Revision 6.8 2004/12/03 21:56:51 bert
49 * Include config.h
50 *
51 * Revision 6.7 2004/11/01 22:38:39 bert
52 * Eliminate all references to minc_def.h
53 *
54 * Revision 6.6 2004/04/27 15:31:45 bert
55 * Added -2 option
56 *
57 * Revision 6.5 2001/04/24 13:38:44 neelin
58 * Replaced NC_NAT with MI_ORIGINAL_TYPE.
59 *
60 * Revision 6.4 2001/04/17 18:40:22 neelin
61 * Modifications to work with NetCDF 3.x
62 * In particular, changed NC_LONG to NC_INT (and corresponding longs to ints).
63 * Changed NC_UNSPECIFIED to NC_NAT.
64 * A few fixes to the configure script.
65 *
66 * Revision 6.3 2001/01/15 14:58:43 neelin
67 * Modified description of -segment option.
68 *
69 * Revision 6.2 2000/07/07 13:19:54 neelin
70 * Added option -filelist to read file names from a file. This gets around
71 * command-line length limits.
72 *
73 * Revision 6.1 1999/10/19 14:45:26 neelin
74 * Fixed Log subsitutions for CVS
75 *
76 * Revision 6.0 1997/09/12 13:24:17 neelin
77 * Release of minc version 0.6
78 *
79 * Revision 5.0 1997/08/21 13:25:16 neelin
80 * Release of minc version 0.5
81 *
82 * Revision 4.0 1997/05/07 20:02:04 neelin
83 * Release of minc version 0.4
84 *
85 * Revision 3.5 1997/04/24 13:48:51 neelin
86 * Fixed handling of invalid or uninitialized data for cumulative operations.
87 * Added options -illegal_value and -count_valid.
88 *
89 * Revision 3.4 1997/04/23 19:34:56 neelin
90 * Added options -maximum, -minimum, -abs.
91 *
92 * Revision 3.3 1996/01/17 21:24:06 neelin
93 * Added -exp and -log options.
94 *
95 * Revision 3.2 1996/01/16 13:29:31 neelin
96 * Added -invert option.
97 *
98 * Revision 3.1 1995/12/13 16:22:24 neelin
99 * Added -check_dimensions and -nocheck_dimensions options.
100 *
101 * Revision 3.0 1995/05/15 19:32:42 neelin
102 * Release of minc version 0.3
103 *
104 * Revision 1.2 1995/05/03 16:13:46 neelin
105 * Changed default for -copy/-nocopy to depend on number of input files.
106 *
107 * Revision 1.1 1995/05/03 13:19:56 neelin
108 * Initial revision
109 *
110 ---------------------------------------------------------------------------- */
111
112 #define _GNU_SOURCE 1
113 #if HAVE_CONFIG_H
114 #include "config.h"
115 #endif
116 #include <stdlib.h>
117 #include <stdio.h>
118 #include <string.h>
119 #include <float.h>
120 #include <limits.h>
121 #include <math.h>
122 #include <minc.h>
123 #include <ParseArgv.h>
124 #include <time_stamp.h>
125 #include <voxel_loop.h>
126 #include "read_file_names.h"
127
128 /* Constants */
129
130 #ifndef TRUE
131 # define TRUE 1
132 # define FALSE 0
133 #endif
134
135 /* Data values for invalid data and for uninitialized data */
136 #define INVALID_DATA -DBL_MAX
137 #define UNINITIALIZED_DATA DBL_MAX
138
139 /* Values for representing default case for command-line options */
140 #define DEFAULT_DBL DBL_MAX
141 #define DEFAULT_BOOL -1
142
143 /* Typedefs */
144 typedef enum {
145 UNSPECIFIED_OP = 0, ADD_OP, SUB_OP, MULT_OP, DIV_OP, SQRT_OP, SQUARE_OP,
146 SCALE_OP, CLAMP_OP, SEGMENT_OP, NSEGMENT_OP, PERCENTDIFF_OP,
147 EQ_OP, NE_OP, GT_OP, GE_OP, LT_OP, LE_OP, AND_OP, OR_OP, NOT_OP,
148 ISNAN_OP, NISNAN_OP, INVERT_OP, EXP_OP, LOG_OP, MAX_OP, MIN_OP, ABS_OP,
149 COUNT_OP
150 } Operation;
151
152 typedef enum {
153 ILLEGAL_NUMOP, UNARY_NUMOP, BINARY_NUMOP, NARY_NUMOP
154 } Num_Operands;
155
156 /* Table that matches [Operation][Number of constants (0,1,2)] to a number
157 of volume operands */
158 Num_Operands OperandTable[][3] = {
159 { ILLEGAL_NUMOP, ILLEGAL_NUMOP, ILLEGAL_NUMOP }, /* UNSPECIFIED_OP */
160 { NARY_NUMOP, UNARY_NUMOP, ILLEGAL_NUMOP }, /* ADD_OP */
161 { BINARY_NUMOP, UNARY_NUMOP, ILLEGAL_NUMOP }, /* SUB_OP */
162 { NARY_NUMOP, UNARY_NUMOP, ILLEGAL_NUMOP }, /* MULT_OP */
163 { BINARY_NUMOP, UNARY_NUMOP, ILLEGAL_NUMOP }, /* DIV_OP */
164 { UNARY_NUMOP, ILLEGAL_NUMOP, ILLEGAL_NUMOP }, /* SQRT_OP */
165 { UNARY_NUMOP, ILLEGAL_NUMOP, ILLEGAL_NUMOP }, /* SQUARE_OP */
166 { ILLEGAL_NUMOP, UNARY_NUMOP, UNARY_NUMOP }, /* SCALE_OP */
167 { ILLEGAL_NUMOP, ILLEGAL_NUMOP, UNARY_NUMOP }, /* CLAMP_OP */
168 { ILLEGAL_NUMOP, ILLEGAL_NUMOP, UNARY_NUMOP }, /* SEGMENT_OP */
169 { ILLEGAL_NUMOP, ILLEGAL_NUMOP, UNARY_NUMOP }, /* NSEGMENT_OP */
170 { BINARY_NUMOP, BINARY_NUMOP, ILLEGAL_NUMOP }, /* PERCENTDIFF_OP */
171 { BINARY_NUMOP, UNARY_NUMOP, ILLEGAL_NUMOP }, /* EQ_OP */
172 { BINARY_NUMOP, UNARY_NUMOP, ILLEGAL_NUMOP }, /* NE_OP */
173 { BINARY_NUMOP, UNARY_NUMOP, ILLEGAL_NUMOP }, /* GT_OP */
174 { BINARY_NUMOP, UNARY_NUMOP, ILLEGAL_NUMOP }, /* GE_OP */
175 { BINARY_NUMOP, UNARY_NUMOP, ILLEGAL_NUMOP }, /* LT_OP */
176 { BINARY_NUMOP, UNARY_NUMOP, ILLEGAL_NUMOP }, /* LE_OP */
177 { NARY_NUMOP, ILLEGAL_NUMOP, ILLEGAL_NUMOP }, /* AND_OP */
178 { NARY_NUMOP, ILLEGAL_NUMOP, ILLEGAL_NUMOP }, /* OR_OP */
179 { UNARY_NUMOP, ILLEGAL_NUMOP, ILLEGAL_NUMOP }, /* NOT_OP */
180 { UNARY_NUMOP, ILLEGAL_NUMOP, ILLEGAL_NUMOP }, /* ISNAN_OP */
181 { UNARY_NUMOP, ILLEGAL_NUMOP, ILLEGAL_NUMOP }, /* NISNAN_OP */
182 { UNARY_NUMOP, UNARY_NUMOP, ILLEGAL_NUMOP }, /* INVERT_OP */
183 { UNARY_NUMOP, UNARY_NUMOP, UNARY_NUMOP }, /* EXP_OP */
184 { UNARY_NUMOP, UNARY_NUMOP, UNARY_NUMOP }, /* LOG_OP */
185 { NARY_NUMOP, ILLEGAL_NUMOP, ILLEGAL_NUMOP }, /* MAX_OP */
186 { NARY_NUMOP, ILLEGAL_NUMOP, ILLEGAL_NUMOP }, /* MIN_OP */
187 { UNARY_NUMOP, ILLEGAL_NUMOP, ILLEGAL_NUMOP }, /* ABS_OP */
188 { NARY_NUMOP, ILLEGAL_NUMOP, ILLEGAL_NUMOP }, /* COUNT_OP */
189 { ILLEGAL_NUMOP, ILLEGAL_NUMOP, ILLEGAL_NUMOP } /* nothing */
190 };
191
192 /* Structure for window information */
193 typedef struct {
194 Operation operation;
195 Num_Operands num_operands;
196 int num_constants;
197 double constants[2];
198 int propagate_nan;
199 double illegal_value;
200 } Math_Data;
201
202 /* Function prototypes */
203 static void do_math(void *caller_data, long num_voxels,
204 int input_num_buffers, int input_vector_length,
205 double *input_data[],
206 int output_num_buffers, int output_vector_length,
207 double *output_data[],
208 Loop_Info *loop_info);
209 static void accum_math(void *caller_data, long num_voxels,
210 int input_num_buffers, int input_vector_length,
211 double *input_data[],
212 int output_num_buffers, int output_vector_length,
213 double *output_data[],
214 Loop_Info *loop_info);
215 static void start_math(void *caller_data, long num_voxels,
216 int output_num_buffers, int output_vector_length,
217 double *output_data[],
218 Loop_Info *loop_info);
219 static void end_math(void *caller_data, long num_voxels,
220 int output_num_buffers, int output_vector_length,
221 double *output_data[],
222 Loop_Info *loop_info);
223
224 /* Argument variables */
225 static int clobber = FALSE;
226 static int verbose = TRUE;
227 static int debug = FALSE;
228 static nc_type datatype = MI_ORIGINAL_TYPE;
229 static int is_signed = FALSE;
230 static double valid_range[2] = {0.0, 0.0};
231 static int copy_all_header = DEFAULT_BOOL;
232 static char *loop_dimension = NULL;
233 static int max_buffer_size_in_kb = 4 * 1024;
234 static double constant = DEFAULT_DBL;
235 static double constant2[2] = {DEFAULT_DBL, DEFAULT_DBL};
236 static Operation operation = UNSPECIFIED_OP;
237 static int propagate_nan = TRUE;
238 static int use_nan_for_illegal_values = TRUE;
239 static double value_for_illegal_operations = DEFAULT_DBL;
240 static int check_dim_info = TRUE;
241 static char *filelist = NULL;
242 #if MINC2
243 static int minc2_format = FALSE;
244 #endif /* MINC2 */
245
246 /* Argument table */
247 static ArgvInfo argTable[] = {
248 {NULL, ARGV_HELP, (char *) NULL, (char *) NULL,
249 "General options:"},
250 #if MINC2
251 {"-2", ARGV_CONSTANT, (char *) TRUE, (char *) &minc2_format,
252 "Produce a MINC 2.0 format output file"},
253 #endif /* MINC2 */
254 {"-clobber", ARGV_CONSTANT, (char *) TRUE, (char *) &clobber,
255 "Overwrite existing file."},
256 {"-noclobber", ARGV_CONSTANT, (char *) FALSE, (char *) &clobber,
257 "Don't overwrite existing file (default)."},
258 {"-no_clobber", ARGV_CONSTANT, (char *) FALSE, (char *) &clobber,
259 "Synonym for -noclobber."},
260 {"-verbose", ARGV_CONSTANT, (char *) TRUE, (char *) &verbose,
261 "Print out log messages (default)."},
262 {"-quiet", ARGV_CONSTANT, (char *) FALSE, (char *) &verbose,
263 "Do not print out log messages."},
264 {"-debug", ARGV_CONSTANT, (char *) TRUE, (char *) &debug,
265 "Print out debugging messages."},
266 {"-filelist", ARGV_STRING, (char *) 1, (char *) &filelist,
267 "Specify the name of a file containing input file names (- for stdin)."},
268 {"-copy_header", ARGV_CONSTANT, (char *) TRUE, (char *) ©_all_header,
269 "Copy all of the header from the first file."},
270 {"-nocopy_header", ARGV_CONSTANT, (char *) FALSE, (char *) ©_all_header,
271 "Do not copy all of the header from the first file."},
272 {"-filetype", ARGV_CONSTANT, (char *) MI_ORIGINAL_TYPE, (char *) &datatype,
273 "Use data type of first file (default)."},
274 {"-byte", ARGV_CONSTANT, (char *) NC_BYTE, (char *) &datatype,
275 "Write out byte data."},
276 {"-short", ARGV_CONSTANT, (char *) NC_SHORT, (char *) &datatype,
277 "Write out short integer data."},
278 {"-int", ARGV_CONSTANT, (char *) NC_INT, (char *) &datatype,
279 "Write out 32-bit integer data."},
280 {"-long", ARGV_CONSTANT, (char *) NC_INT, (char *) &datatype,
281 "Superseded by -int."},
282 {"-float", ARGV_CONSTANT, (char *) NC_FLOAT, (char *) &datatype,
283 "Write out single-precision floating-point data."},
284 {"-double", ARGV_CONSTANT, (char *) NC_DOUBLE, (char *) &datatype,
285 "Write out double-precision floating-point data."},
286 {"-signed", ARGV_CONSTANT, (char *) TRUE, (char *) &is_signed,
287 "Write signed integer data."},
288 {"-unsigned", ARGV_CONSTANT, (char *) FALSE, (char *) &is_signed,
289 "Write unsigned integer data (default if type specified)."},
290 {"-range", ARGV_FLOAT, (char *) 2, (char *) valid_range,
291 "Valid range for output data."},
292 {"-max_buffer_size_in_kb", ARGV_INT, (char *) 1,
293 (char *) &max_buffer_size_in_kb,
294 "Specify the maximum size of the internal buffers (in kbytes)."},
295 {"-dimension", ARGV_STRING, (char *) 1, (char *) &loop_dimension,
296 "Specify a dimension along which we wish to perform a calculation."},
297 {"-check_dimensions", ARGV_CONSTANT, (char *) TRUE,
298 (char *) &check_dim_info,
299 "Check that files have matching dimensions (default)."},
300 {"-nocheck_dimensions", ARGV_CONSTANT, (char *) FALSE,
301 (char *) &check_dim_info,
302 "Do not check that files have matching dimensions."},
303 {"-ignore_nan", ARGV_CONSTANT, (char *) FALSE, (char *) &propagate_nan,
304 "Ignore invalid data (NaN) for accumulations."},
305 {"-propagate_nan", ARGV_CONSTANT, (char *) TRUE, (char *) &propagate_nan,
306 "Invalid data in any file at a voxel produces a NaN (default)."},
307 {"-nan", ARGV_CONSTANT, (char *) TRUE,
308 (char *) &use_nan_for_illegal_values,
309 "Output NaN when an illegal operation is done (default)."},
310 {"-zero", ARGV_CONSTANT, (char *) FALSE,
311 (char *) &use_nan_for_illegal_values,
312 "Output zero when an illegal operation is done."},
313 {"-illegal_value", ARGV_FLOAT, (char *) 1,
314 (char *) &value_for_illegal_operations,
315 "Value to write out when an illegal operation is done."},
316 {NULL, ARGV_HELP, (char *) NULL, (char *) NULL,
317 "Options for specifying constants:"},
318 {"-constant", ARGV_FLOAT, (char *) 1, (char *) &constant,
319 "Specify a constant argument."},
320 {"-const", ARGV_FLOAT, (char *) 1, (char *) &constant,
321 "Synonym for -constant."},
322 {"-const2", ARGV_FLOAT, (char *) 2, (char *) constant2,
323 "Specify two constant arguments."},
324 {NULL, ARGV_HELP, (char *) NULL, (char *) NULL,
325 "Operations:"},
326 {"-add", ARGV_CONSTANT, (char *) ADD_OP, (char *) &operation,
327 "Add N volumes or volume + constant."},
328 {"-sub", ARGV_CONSTANT, (char *) SUB_OP, (char *) &operation,
329 "Subtract 2 volumes or volume - constant."},
330 {"-mult", ARGV_CONSTANT, (char *) MULT_OP, (char *) &operation,
331 "Multiply N volumes or volume * constant."},
332 {"-div", ARGV_CONSTANT, (char *) DIV_OP, (char *) &operation,
333 "Divide 2 volumes or volume / constant."},
334 {"-invert", ARGV_CONSTANT, (char *) INVERT_OP, (char *) &operation,
335 "Calculate 1/x at each voxel (use -constant for c/x)."},
336 {"-sqrt", ARGV_CONSTANT, (char *) SQRT_OP, (char *) &operation,
337 "Take square root of a volume."},
338 {"-square", ARGV_CONSTANT, (char *) SQUARE_OP, (char *) &operation,
339 "Take square of a volume."},
340 {"-abs", ARGV_CONSTANT, (char *) ABS_OP, (char *) &operation,
341 "Take absolute value of a volume."},
342 {"-max", ARGV_CONSTANT, (char *) MAX_OP, (char *) &operation,
343 "Synonym for -maximum."},
344 {"-maximum", ARGV_CONSTANT, (char *) MAX_OP, (char *) &operation,
345 "Find maximum of N volumes."},
346 {"-minimum", ARGV_CONSTANT, (char *) MIN_OP, (char *) &operation,
347 "Find minimum of N volumes."},
348 {"-exp", ARGV_CONSTANT, (char *) EXP_OP, (char *) &operation,
349 "Calculate c2*exp(c1*x). The constants c1 and c2 default to 1."},
350 {"-log", ARGV_CONSTANT, (char *) LOG_OP, (char *) &operation,
351 "Calculate log(x/c2)/c1. The constants c1 and c2 default to 1."},
352 {"-scale", ARGV_CONSTANT, (char *) SCALE_OP, (char *) &operation,
353 "Scale a volume: volume * c1 + c2."},
354 {"-clamp", ARGV_CONSTANT, (char *) CLAMP_OP, (char *) &operation,
355 "Clamp a volume to lie between two values."},
356 {"-segment", ARGV_CONSTANT, (char *) SEGMENT_OP, (char *) &operation,
357 "Segment a volume using range of -const2: within range = 1, outside range = 0."},
358 {"-nsegment", ARGV_CONSTANT, (char *) NSEGMENT_OP, (char *) &operation,
359 "Opposite of -segment: within range = 0, outside range = 1."},
360 {"-percentdiff", ARGV_CONSTANT, (char *) PERCENTDIFF_OP,
361 (char *) &operation,
362 "Percent difference between 2 volumes, thresholded (const def=0.0)."},
363 {"-pd", ARGV_CONSTANT, (char *) PERCENTDIFF_OP, (char *) &operation,
364 "Synonym for -percentdiff."},
365 {"-eq", ARGV_CONSTANT, (char *) EQ_OP, (char *) &operation,
366 "Test for integer vol1 == vol2 or vol1 == const."},
367 {"-ne", ARGV_CONSTANT, (char *) NE_OP, (char *) &operation,
368 "Test for integer vol1 != vol2 or vol1 != const."},
369 {"-gt", ARGV_CONSTANT, (char *) GT_OP, (char *) &operation,
370 "Test for vol1 > vol2 or vol1 > const."},
371 {"-ge", ARGV_CONSTANT, (char *) GE_OP, (char *) &operation,
372 "Test for vol1 >= vol2 or vol1 >= const."},
373 {"-lt", ARGV_CONSTANT, (char *) LT_OP, (char *) &operation,
374 "Test for vol1 < vol2 or vol1 < const."},
375 {"-le", ARGV_CONSTANT, (char *) LE_OP, (char *) &operation,
376 "Test for vol1 <= vol2 or vol1 <= const."},
377 {"-and", ARGV_CONSTANT, (char *) AND_OP, (char *) &operation,
378 "Calculate vol1 && vol2 (&& ...)."},
379 {"-or", ARGV_CONSTANT, (char *) OR_OP, (char *) &operation,
380 "Calculate vol1 || vol2 (|| ...)."},
381 {"-not", ARGV_CONSTANT, (char *) NOT_OP, (char *) &operation,
382 "Calculate !vol1."},
383 {"-isnan", ARGV_CONSTANT, (char *) ISNAN_OP, (char *) &operation,
384 "Test for NaN values in vol1."},
385 {"-nisnan", ARGV_CONSTANT, (char *) NISNAN_OP, (char *) &operation,
386 "Negation of -isnan."},
387 {"-count_valid", ARGV_CONSTANT, (char *) COUNT_OP, (char *) &operation,
388 "Count the number of valid values in N volumes."},
389 {NULL, ARGV_END, NULL, NULL, NULL}
390 };
391
392 /* Main program */
393
main(int argc,char * argv[])394 int main(int argc, char *argv[])
395 {
396 char **infiles, **outfiles;
397 int nfiles, nout;
398 char *arg_string;
399 Math_Data math_data;
400 Loop_Options *loop_options;
401 char *pname;
402 int num_constants;
403 Num_Operands num_operands;
404 VoxelFunction math_function;
405
406 /* Save time stamp and args */
407 arg_string = time_stamp(argc, argv);
408
409 /* Get arguments */
410 pname = argv[0];
411 if (ParseArgv(&argc, argv, argTable, 0) || (argc < 2)) {
412 (void) fprintf(stderr,
413 "\nUsage: %s [options] [<in1.mnc> ...] <out.mnc>\n",
414 pname);
415 (void) fprintf(stderr,
416 " %s -help\n\n", pname);
417 exit(EXIT_FAILURE);
418 }
419 nout = 1;
420 outfiles = &argv[argc-1];
421
422 /* Get the list of input files either from the command line or
423 from a file, or report an error if both are specified */
424 nfiles = argc - 2;
425 if (filelist == NULL) {
426 infiles = &argv[1];
427 }
428 else if (nfiles <= 0) {
429 infiles = read_file_names(filelist, &nfiles);
430 if (infiles == NULL) {
431 (void) fprintf(stderr,
432 "Error reading in file names from file \"%s\"\n",
433 filelist);
434 exit(EXIT_FAILURE);
435 }
436 }
437 else {
438 (void) fprintf(stderr,
439 "Do not specify both -filelist and input file names\n");
440 exit(EXIT_FAILURE);
441 }
442
443 /* Make sure that we have something to process */
444 if (nfiles == 0) {
445 (void) fprintf(stderr, "No input files specified\n");
446 exit(EXIT_FAILURE);
447 }
448
449 /* Handle special case of COUNT_OP - it always assume -ignore_nan and
450 -zero */
451 if (operation == COUNT_OP) {
452 propagate_nan = FALSE;
453 value_for_illegal_operations = 0.0;
454 }
455
456 /* Check that the arguments make sense */
457 if ((constant != DEFAULT_DBL) && (constant2[0] != DEFAULT_DBL)) {
458 (void) fprintf(stderr, "%s: Specify only one of -constant or -const2\n",
459 pname);
460 exit(EXIT_FAILURE);
461 }
462 if (constant != DEFAULT_DBL)
463 num_constants = 1;
464 else if (constant2[0] != DEFAULT_DBL)
465 num_constants = 2;
466 else
467 num_constants = 0;
468 num_operands = OperandTable[operation][num_constants];
469 if (num_operands == ILLEGAL_NUMOP) {
470 (void) fprintf(stderr, "%s: Operation and constants do not match.\n",
471 pname);
472 exit(EXIT_FAILURE);
473 }
474 if ((num_operands != NARY_NUMOP) && (loop_dimension != NULL)) {
475 (void) fprintf(stderr, "%s: Use -dimension only for cumulative ops.\n",
476 pname);
477 exit(EXIT_FAILURE);
478 }
479 if ((num_operands == UNARY_NUMOP) && (nfiles != 1)) {
480 (void) fprintf(stderr, "%s: Expected only one input file.\n", pname);
481 exit(EXIT_FAILURE);
482 }
483 if ((num_operands == BINARY_NUMOP) && (nfiles != 2)) {
484 (void) fprintf(stderr, "%s: Expected two input files.\n", pname);
485 exit(EXIT_FAILURE);
486 }
487 if ((num_operands == NARY_NUMOP) && (nfiles < 1) &&
488 (loop_dimension == NULL)) {
489 (void) fprintf(stderr, "%s: Expected at least one input files.\n",
490 pname);
491 exit(EXIT_FAILURE);
492 }
493
494 /* Set default copy_all_header according to number of input files */
495 if (copy_all_header == DEFAULT_BOOL)
496 copy_all_header = (nfiles == 1);
497
498 /* Set up math data structure */
499 math_data.operation = operation;
500 math_data.num_operands = num_operands;
501 math_data.propagate_nan = propagate_nan;
502 math_data.num_constants = num_constants;
503 switch (num_constants) {
504 case 1:
505 math_data.constants[0] = constant;
506 break;
507 case 2:
508 math_data.constants[0] = constant2[0];
509 math_data.constants[1] = constant2[1];
510 break;
511 }
512 if (value_for_illegal_operations != DEFAULT_DBL)
513 math_data.illegal_value = value_for_illegal_operations;
514 else if (use_nan_for_illegal_values)
515 math_data.illegal_value = INVALID_DATA;
516 else
517 math_data.illegal_value = 0.0;
518
519 /* Do math */
520 loop_options = create_loop_options();
521 set_loop_verbose(loop_options, verbose);
522 set_loop_clobber(loop_options, clobber);
523 #if MINC2
524 set_loop_v2format(loop_options, minc2_format);
525 #endif /* MINC2 */
526 set_loop_datatype(loop_options, datatype, is_signed,
527 valid_range[0], valid_range[1]);
528 if (num_operands == NARY_NUMOP) {
529 math_function = accum_math;
530 set_loop_accumulate(loop_options, TRUE, 0, start_math, end_math);
531 }
532 else {
533 math_function = do_math;
534 }
535 set_loop_copy_all_header(loop_options, copy_all_header);
536 set_loop_dimension(loop_options, loop_dimension);
537 set_loop_buffer_size(loop_options, (long) 1024 * max_buffer_size_in_kb);
538 set_loop_check_dim_info(loop_options, check_dim_info);
539 voxel_loop(nfiles, infiles, nout, outfiles, arg_string, loop_options,
540 math_function, (void *) &math_data);
541 free_loop_options(loop_options);
542
543 exit(EXIT_SUCCESS);
544 }
545
546 /* ----------------------------- MNI Header -----------------------------------
547 @NAME : do_math
548 @INPUT : Standard for voxel loop
549 @OUTPUT : Standard for voxel loop
550 @RETURNS : (nothing)
551 @DESCRIPTION: Routine doing math operations.
552 @METHOD :
553 @GLOBALS :
554 @CALLS :
555 @CREATED : April 25, 1995 (Peter Neelin)
556 @MODIFIED :
557 ---------------------------------------------------------------------------- */
do_math(void * caller_data,long num_voxels,int input_num_buffers,int input_vector_length,double * input_data[],int output_num_buffers,int output_vector_length,double * output_data[],Loop_Info * loop_info)558 static void do_math(void *caller_data, long num_voxels,
559 int input_num_buffers, int input_vector_length,
560 double *input_data[],
561 int output_num_buffers, int output_vector_length,
562 double *output_data[],
563 Loop_Info *loop_info)
564 /* ARGSUSED */
565 {
566 Math_Data *math_data;
567 long ivox;
568 double value1, value2;
569 double illegal_value;
570 Operation operation;
571 int num_constants, iconst;
572 double constants[2];
573
574 /* Get pointer to window info */
575 math_data = (Math_Data *) caller_data;
576
577 /* Check arguments */
578 if ((input_num_buffers > 2) || (output_num_buffers != 1) ||
579 (output_vector_length != input_vector_length)) {
580 (void) fprintf(stderr, "Bad arguments to do_math!\n");
581 exit(EXIT_FAILURE);
582 }
583
584 /* Get info */
585 operation = math_data->operation;
586 num_constants = math_data->num_constants;
587 for (iconst=0; iconst < sizeof(constants)/sizeof(constants[0]); iconst++) {
588 if (iconst < num_constants)
589 constants[iconst] = math_data->constants[iconst];
590 else if ((operation == INVERT_OP) ||
591 (operation == EXP_OP) ||
592 (operation == LOG_OP))
593 constants[iconst] = 1.0;
594 else
595 constants[iconst] = 0.0;
596 }
597 illegal_value = math_data->illegal_value;
598
599 /* Set default second value */
600 value2 = constants[0];
601
602 /* Loop through the voxels */
603 for (ivox=0; ivox < num_voxels*input_vector_length; ivox++) {
604 value1 = input_data[0][ivox];
605 if (input_num_buffers == 2)
606 value2 = input_data[1][ivox];
607 if ((value1 == INVALID_DATA) || (value2 == INVALID_DATA)) {
608 switch(operation) {
609 case ISNAN_OP:
610 output_data[0][ivox] = 1.0;
611 break;
612 case NISNAN_OP:
613 output_data[0][ivox] = 0.0;
614 break;
615 default:
616 output_data[0][ivox] = INVALID_DATA;
617 break;
618 }
619 }
620 else {
621 switch (operation) {
622 case ADD_OP:
623 output_data[0][ivox] = value1 + value2; break;
624 case SUB_OP:
625 output_data[0][ivox] = value1 - value2; break;
626 case MULT_OP:
627 output_data[0][ivox] = value1 * value2; break;
628 case DIV_OP:
629 if (value2 != 0.0)
630 output_data[0][ivox] = value1 / value2;
631 else
632 output_data[0][ivox] = illegal_value;
633 break;
634 case INVERT_OP:
635 if (value1 == 0.0)
636 output_data[0][ivox] = illegal_value;
637 else
638 output_data[0][ivox] = value2 / value1;
639 break;
640 case SQRT_OP:
641 if (value1 < 0.0)
642 output_data[0][ivox] = illegal_value;
643 else
644 output_data[0][ivox] = sqrt(value1);
645 break;
646 case SQUARE_OP:
647 output_data[0][ivox] = value1 * value1; break;
648 case ABS_OP:
649 if (value1 < 0.0)
650 output_data[0][ivox] = -value1;
651 else
652 output_data[0][ivox] = value1;
653 break;
654 case EXP_OP:
655 output_data[0][ivox] = constants[1] * exp(value1 * constants[0]);
656 break;
657 case LOG_OP:
658 if ((value1 <= 0.0) || (constants[1] <= 0.0) ||
659 (constants[0] == 0.0))
660 output_data[0][ivox] = illegal_value;
661 else
662 output_data[0][ivox] = log(value1/constants[1])/constants[0];
663 break;
664 case SCALE_OP:
665 output_data[0][ivox] = value1 * constants[0] + constants[1]; break;
666 case CLAMP_OP:
667 if (value1 < constants[0])
668 value1 = constants[0];
669 else if (value1 > constants[1])
670 value1 = constants[1];
671 output_data[0][ivox] = value1;
672 break;
673 case SEGMENT_OP:
674 if ((value1 < constants[0]) || (value1 > constants[1]))
675 output_data[0][ivox] = 0.0;
676 else
677 output_data[0][ivox] = 1.0;
678 break;
679 case NSEGMENT_OP:
680 if ((value1 < constants[0]) || (value1 > constants[1]))
681 output_data[0][ivox] = 1.0;
682 else
683 output_data[0][ivox] = 0.0;
684 break;
685 case PERCENTDIFF_OP:
686 if ((value1 < constants[0]) || (value1 == 0.0))
687 output_data[0][ivox] = illegal_value;
688 else {
689 output_data[0][ivox] = 100.0 * (value1 - value2) / value1;
690 }
691 break;
692 case EQ_OP:
693 output_data[0][ivox] = (((rint(value1)-rint(value2)) == 0.0)
694 ? 1.0 : 0.0);
695 break;
696 case NE_OP:
697 output_data[0][ivox] = (((rint(value1)-rint(value2)) != 0.0)
698 ? 1.0 : 0.0);
699 break;
700 case GT_OP:
701 output_data[0][ivox] = value1 > value2; break;
702 case GE_OP:
703 output_data[0][ivox] = value1 >= value2; break;
704 case LT_OP:
705 output_data[0][ivox] = value1 < value2; break;
706 case LE_OP:
707 output_data[0][ivox] = value1 <= value2; break;
708 case AND_OP:
709 output_data[0][ivox] =
710 (((rint(value1) != 0.0) && (rint(value2) != 0.0)) ? 1.0 : 0.0);
711 break;
712 case OR_OP:
713 output_data[0][ivox] =
714 (((rint(value1) != 0.0) || (rint(value2) != 0.0)) ? 1.0 : 0.0);
715 break;
716 case NOT_OP:
717 output_data[0][ivox] = ((rint(value1) == 0.0) ? 1.0 : 0.0);
718 break;
719 case ISNAN_OP:
720 output_data[0][ivox] = 0.0; /* To get here, value is not nan */
721 break;
722 case NISNAN_OP:
723 output_data[0][ivox] = 1.0;
724 break;
725 default:
726 (void) fprintf(stderr, "Bad op in do_math!\n");
727 exit(EXIT_FAILURE);
728 }
729 }
730 }
731
732 return;
733 }
734
735 /* ----------------------------- MNI Header -----------------------------------
736 @NAME : accum_math
737 @INPUT : Standard for voxel loop
738 @OUTPUT : Standard for voxel loop
739 @RETURNS : (nothing)
740 @DESCRIPTION: Routine for doing accumulation math operations.
741 @METHOD :
742 @GLOBALS :
743 @CALLS :
744 @CREATED : April 25, 1995 (Peter Neelin)
745 @MODIFIED :
746 ---------------------------------------------------------------------------- */
accum_math(void * caller_data,long num_voxels,int input_num_buffers,int input_vector_length,double * input_data[],int output_num_buffers,int output_vector_length,double * output_data[],Loop_Info * loop_info)747 static void accum_math(void *caller_data, long num_voxels,
748 int input_num_buffers, int input_vector_length,
749 double *input_data[],
750 int output_num_buffers, int output_vector_length,
751 double *output_data[],
752 Loop_Info *loop_info)
753 /* ARGSUSED */
754 {
755 Math_Data *math_data;
756 long ivox;
757 double value, oldvalue;
758 Operation operation;
759 int propagate_nan;
760
761 /* Get pointer to window info */
762 math_data = (Math_Data *) caller_data;
763
764 /* Check arguments */
765 if ((input_num_buffers != 1) || (output_num_buffers != 1) ||
766 (output_vector_length != input_vector_length)) {
767 (void) fprintf(stderr, "Bad arguments to accum_math!\n");
768 exit(EXIT_FAILURE);
769 }
770
771 /* Get info */
772 operation = math_data->operation;
773 propagate_nan = math_data->propagate_nan;
774
775 /* Loop through the voxels */
776 for (ivox=0; ivox < num_voxels*input_vector_length; ivox++) {
777
778 /* Get previous value and the next value */
779 oldvalue = output_data[0][ivox];
780 value = input_data[0][ivox];
781
782 /* If the new data is invalid, then either mark the output as invalid
783 or ignore it */
784 if (value == INVALID_DATA) {
785 if (propagate_nan) {
786 output_data[0][ivox] = INVALID_DATA;
787 }
788 }
789
790 /* If we haven't set anything yet, then just copy the new value */
791 else if (oldvalue == UNINITIALIZED_DATA) {
792 output_data[0][ivox] = input_data[0][ivox];
793 }
794
795 /* Do the operation if the old data and the new data are valid */
796 else if (oldvalue != INVALID_DATA) {
797 switch (operation) {
798 case ADD_OP:
799 output_data[0][ivox] = oldvalue + value;
800 break;
801 case MULT_OP:
802 output_data[0][ivox] = oldvalue * value;
803 break;
804 case AND_OP:
805 output_data[0][ivox] =
806 (((oldvalue != 0.0) && (rint(value) != 0.0)) ? 1.0 : 0.0);
807 break;
808 case OR_OP:
809 output_data[0][ivox] =
810 (((oldvalue != 0.0) || (rint(value) != 0.0)) ? 1.0 : 0.0);
811 break;
812 case MAX_OP:
813 if (value > oldvalue)
814 output_data[0][ivox] = value;
815 break;
816 case MIN_OP:
817 if (value < oldvalue)
818 output_data[0][ivox] = value;
819 break;
820 case COUNT_OP:
821 output_data[0][ivox]++;
822 break;
823 default:
824 (void) fprintf(stderr, "Bad op in accum_math!\n");
825 exit(EXIT_FAILURE);
826 }
827 }
828
829 } /* Loop over voxels */
830
831 return;
832 }
833
834 /* ----------------------------- MNI Header -----------------------------------
835 @NAME : start_math
836 @INPUT : Standard for voxel loop
837 @OUTPUT : Standard for voxel loop
838 @RETURNS : (nothing)
839 @DESCRIPTION: Start routine for math accumulation.
840 @METHOD :
841 @GLOBALS :
842 @CALLS :
843 @CREATED : April 25, 1995 (Peter Neelin)
844 @MODIFIED :
845 ---------------------------------------------------------------------------- */
start_math(void * caller_data,long num_voxels,int output_num_buffers,int output_vector_length,double * output_data[],Loop_Info * loop_info)846 static void start_math(void *caller_data, long num_voxels,
847 int output_num_buffers, int output_vector_length,
848 double *output_data[],
849 Loop_Info *loop_info)
850 /* ARGSUSED */
851 {
852 Math_Data *math_data;
853 long ivox;
854
855 /* Get pointer to window info */
856 math_data = (Math_Data *) caller_data;
857
858 /* Check arguments */
859 if (output_num_buffers != 1) {
860 (void) fprintf(stderr, "Bad arguments to start_math!\n");
861 exit(EXIT_FAILURE);
862 }
863
864 /* Get info */
865 operation = math_data->operation;
866
867 /* Loop through the voxels, marking them all as uninitialized. We treat
868 COUNT_OP as a special case since it always has a value. This is
869 especially important to prevent it from going through
870 the code in accum_math for handling the first valid voxel which
871 just assigns the first value. */
872 for (ivox=0; ivox < num_voxels*output_vector_length; ivox++) {
873 switch (operation) {
874 case COUNT_OP:
875 output_data[0][ivox] = 0.0;
876 break;
877 default:
878 output_data[0][ivox] = UNINITIALIZED_DATA;
879 break;
880 }
881 }
882
883 return;
884 }
885
886 /* ----------------------------- MNI Header -----------------------------------
887 @NAME : end_math
888 @INPUT : Standard for voxel loop
889 @OUTPUT : Standard for voxel loop
890 @RETURNS : (nothing)
891 @DESCRIPTION: Start routine for math accumulation.
892 @METHOD :
893 @GLOBALS :
894 @CALLS :
895 @CREATED : April 25, 1995 (Peter Neelin)
896 @MODIFIED :
897 ---------------------------------------------------------------------------- */
end_math(void * caller_data,long num_voxels,int output_num_buffers,int output_vector_length,double * output_data[],Loop_Info * loop_info)898 static void end_math(void *caller_data, long num_voxels,
899 int output_num_buffers, int output_vector_length,
900 double *output_data[],
901 Loop_Info *loop_info)
902 /* ARGSUSED */
903 {
904 Math_Data *math_data;
905 long ivox;
906 double value;
907 double illegal_value;
908
909 /* Get pointer to window info */
910 math_data = (Math_Data *) caller_data;
911
912 /* Check arguments */
913 if (output_num_buffers != 1) {
914 (void) fprintf(stderr, "Bad arguments to end_math!\n");
915 exit(EXIT_FAILURE);
916 }
917
918 /* Get info */
919 operation = math_data->operation;
920 illegal_value = math_data->illegal_value;
921
922 /* Loop through the voxels, checking for uninitialized values */
923 for (ivox=0; ivox < num_voxels*output_vector_length; ivox++) {
924 value = output_data[0][ivox];
925 if ((value == UNINITIALIZED_DATA) || (value == INVALID_DATA)) {
926 output_data[0][ivox] = illegal_value;
927 }
928 }
929
930 return;
931 }
932