/* ----------------------------- MNI Header ----------------------------------- @NAME : mincmath @INPUT : argc, argv - command line arguments @OUTPUT : (none) @RETURNS : status @DESCRIPTION: Program to do simple math on minc files @METHOD : @GLOBALS : @CALLS : @CREATED : April 28, 1995 (Peter Neelin) @MODIFIED : * $Log: mincmath.c,v $ * Revision 6.9 2004/12/14 23:40:07 bert * Get rid of c99 compilation problems * * Revision 6.8 2004/12/03 21:56:51 bert * Include config.h * * Revision 6.7 2004/11/01 22:38:39 bert * Eliminate all references to minc_def.h * * Revision 6.6 2004/04/27 15:31:45 bert * Added -2 option * * Revision 6.5 2001/04/24 13:38:44 neelin * Replaced NC_NAT with MI_ORIGINAL_TYPE. * * Revision 6.4 2001/04/17 18:40:22 neelin * Modifications to work with NetCDF 3.x * In particular, changed NC_LONG to NC_INT (and corresponding longs to ints). * Changed NC_UNSPECIFIED to NC_NAT. * A few fixes to the configure script. * * Revision 6.3 2001/01/15 14:58:43 neelin * Modified description of -segment option. * * Revision 6.2 2000/07/07 13:19:54 neelin * Added option -filelist to read file names from a file. This gets around * command-line length limits. * * Revision 6.1 1999/10/19 14:45:26 neelin * Fixed Log subsitutions for CVS * * Revision 6.0 1997/09/12 13:24:17 neelin * Release of minc version 0.6 * * Revision 5.0 1997/08/21 13:25:16 neelin * Release of minc version 0.5 * * Revision 4.0 1997/05/07 20:02:04 neelin * Release of minc version 0.4 * * Revision 3.5 1997/04/24 13:48:51 neelin * Fixed handling of invalid or uninitialized data for cumulative operations. * Added options -illegal_value and -count_valid. * * Revision 3.4 1997/04/23 19:34:56 neelin * Added options -maximum, -minimum, -abs. * * Revision 3.3 1996/01/17 21:24:06 neelin * Added -exp and -log options. * * Revision 3.2 1996/01/16 13:29:31 neelin * Added -invert option. * * Revision 3.1 1995/12/13 16:22:24 neelin * Added -check_dimensions and -nocheck_dimensions options. * * Revision 3.0 1995/05/15 19:32:42 neelin * Release of minc version 0.3 * * Revision 1.2 1995/05/03 16:13:46 neelin * Changed default for -copy/-nocopy to depend on number of input files. * * Revision 1.1 1995/05/03 13:19:56 neelin * Initial revision * ---------------------------------------------------------------------------- */ #ifndef lint static char rcsid[]="$Header: /software/source/minc/progs/mincmath/mincmath.c,v 6.9 2004/12/14 23:40:07 bert Exp $"; #endif #define _GNU_SOURCE 1 #include "config.h" #include #include #include #include #include #include #include #include #include #include /* Constants */ #ifndef TRUE # define TRUE 1 # define FALSE 0 #endif /* Data values for invalid data and for uninitialized data */ #define INVALID_DATA -DBL_MAX #define UNINITIALIZED_DATA DBL_MAX /* Values for representing default case for command-line options */ #define DEFAULT_DBL DBL_MAX #define DEFAULT_BOOL -1 /* Typedefs */ typedef enum { UNSPECIFIED_OP = 0, ADD_OP, SUB_OP, MULT_OP, DIV_OP, SQRT_OP, SQUARE_OP, SCALE_OP, CLAMP_OP, SEGMENT_OP, NSEGMENT_OP, PERCENTDIFF_OP, EQ_OP, NE_OP, GT_OP, GE_OP, LT_OP, LE_OP, AND_OP, OR_OP, NOT_OP, ISNAN_OP, NISNAN_OP, INVERT_OP, EXP_OP, LOG_OP, MAX_OP, MIN_OP, ABS_OP, COUNT_OP } Operation; typedef enum { ILLEGAL_NUMOP, UNARY_NUMOP, BINARY_NUMOP, NARY_NUMOP } Num_Operands; /* Table that matches [Operation][Number of constants (0,1,2)] to a number of volume operands */ Num_Operands OperandTable[][3] = { ILLEGAL_NUMOP, ILLEGAL_NUMOP, ILLEGAL_NUMOP, /* UNSPECIFIED_OP */ NARY_NUMOP, UNARY_NUMOP, ILLEGAL_NUMOP, /* ADD_OP */ BINARY_NUMOP, UNARY_NUMOP, ILLEGAL_NUMOP, /* SUB_OP */ NARY_NUMOP, UNARY_NUMOP, ILLEGAL_NUMOP, /* MULT_OP */ BINARY_NUMOP, UNARY_NUMOP, ILLEGAL_NUMOP, /* DIV_OP */ UNARY_NUMOP, ILLEGAL_NUMOP, ILLEGAL_NUMOP, /* SQRT_OP */ UNARY_NUMOP, ILLEGAL_NUMOP, ILLEGAL_NUMOP, /* SQUARE_OP */ ILLEGAL_NUMOP, UNARY_NUMOP, UNARY_NUMOP, /* SCALE_OP */ ILLEGAL_NUMOP, ILLEGAL_NUMOP, UNARY_NUMOP, /* CLAMP_OP */ ILLEGAL_NUMOP, ILLEGAL_NUMOP, UNARY_NUMOP, /* SEGMENT_OP */ ILLEGAL_NUMOP, ILLEGAL_NUMOP, UNARY_NUMOP, /* NSEGMENT_OP */ BINARY_NUMOP, BINARY_NUMOP, ILLEGAL_NUMOP, /* PERCENTDIFF_OP */ BINARY_NUMOP, UNARY_NUMOP, ILLEGAL_NUMOP, /* EQ_OP */ BINARY_NUMOP, UNARY_NUMOP, ILLEGAL_NUMOP, /* NE_OP */ BINARY_NUMOP, UNARY_NUMOP, ILLEGAL_NUMOP, /* GT_OP */ BINARY_NUMOP, UNARY_NUMOP, ILLEGAL_NUMOP, /* GE_OP */ BINARY_NUMOP, UNARY_NUMOP, ILLEGAL_NUMOP, /* LT_OP */ BINARY_NUMOP, UNARY_NUMOP, ILLEGAL_NUMOP, /* LE_OP */ NARY_NUMOP, ILLEGAL_NUMOP, ILLEGAL_NUMOP, /* AND_OP */ NARY_NUMOP, ILLEGAL_NUMOP, ILLEGAL_NUMOP, /* OR_OP */ UNARY_NUMOP, ILLEGAL_NUMOP, ILLEGAL_NUMOP, /* NOT_OP */ UNARY_NUMOP, ILLEGAL_NUMOP, ILLEGAL_NUMOP, /* ISNAN_OP */ UNARY_NUMOP, ILLEGAL_NUMOP, ILLEGAL_NUMOP, /* NISNAN_OP */ UNARY_NUMOP, UNARY_NUMOP, ILLEGAL_NUMOP, /* INVERT_OP */ UNARY_NUMOP, UNARY_NUMOP, UNARY_NUMOP, /* EXP_OP */ UNARY_NUMOP, UNARY_NUMOP, UNARY_NUMOP, /* LOG_OP */ NARY_NUMOP, ILLEGAL_NUMOP, ILLEGAL_NUMOP, /* MAX_OP */ NARY_NUMOP, ILLEGAL_NUMOP, ILLEGAL_NUMOP, /* MIN_OP */ UNARY_NUMOP, ILLEGAL_NUMOP, ILLEGAL_NUMOP, /* ABS_OP */ NARY_NUMOP, ILLEGAL_NUMOP, ILLEGAL_NUMOP, /* COUNT_OP */ ILLEGAL_NUMOP, ILLEGAL_NUMOP, ILLEGAL_NUMOP /* nothing */ }; /* Structure for window information */ typedef struct { Operation operation; Num_Operands num_operands; int num_constants; double constants[2]; int propagate_nan; double illegal_value; } Math_Data; /* Function prototypes */ static void 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); static void 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); static void start_math(void *caller_data, long num_voxels, int output_num_buffers, int output_vector_length, double *output_data[], Loop_Info *loop_info); static void end_math(void *caller_data, long num_voxels, int output_num_buffers, int output_vector_length, double *output_data[], Loop_Info *loop_info); static char **read_file_names(char *filelist, int *num_files); /* Argument variables */ int clobber = FALSE; int verbose = TRUE; int debug = FALSE; nc_type datatype = MI_ORIGINAL_TYPE; int is_signed = FALSE; double valid_range[2] = {0.0, 0.0}; int copy_all_header = DEFAULT_BOOL; char *loop_dimension = NULL; int max_buffer_size_in_kb = 4 * 1024; double constant = DEFAULT_DBL; double constant2[2] = {DEFAULT_DBL, DEFAULT_DBL}; Operation operation = UNSPECIFIED_OP; int propagate_nan = TRUE; int use_nan_for_illegal_values = TRUE; double value_for_illegal_operations = DEFAULT_DBL; int check_dim_info = TRUE; char *filelist = NULL; #ifdef MINC2 int minc2_format = FALSE; #endif /* MINC2 defined */ /* Argument table */ ArgvInfo argTable[] = { {NULL, ARGV_HELP, (char *) NULL, (char *) NULL, "General options:"}, #ifdef MINC2 {"-2", ARGV_CONSTANT, (char *) TRUE, (char *) &minc2_format, "Produce a MINC 2.0 format output file"}, #endif /* MINC2 defined */ {"-clobber", ARGV_CONSTANT, (char *) TRUE, (char *) &clobber, "Overwrite existing file."}, {"-noclobber", ARGV_CONSTANT, (char *) FALSE, (char *) &clobber, "Don't overwrite existing file (default)."}, {"-no_clobber", ARGV_CONSTANT, (char *) FALSE, (char *) &clobber, "Synonym for -noclobber."}, {"-verbose", ARGV_CONSTANT, (char *) TRUE, (char *) &verbose, "Print out log messages (default)."}, {"-quiet", ARGV_CONSTANT, (char *) FALSE, (char *) &verbose, "Do not print out log messages."}, {"-debug", ARGV_CONSTANT, (char *) TRUE, (char *) &debug, "Print out debugging messages."}, {"-filelist", ARGV_STRING, (char *) 1, (char *) &filelist, "Specify the name of a file containing input file names (- for stdin)."}, {"-copy_header", ARGV_CONSTANT, (char *) TRUE, (char *) ©_all_header, "Copy all of the header from the first file."}, {"-nocopy_header", ARGV_CONSTANT, (char *) FALSE, (char *) ©_all_header, "Do not copy all of the header from the first file."}, {"-filetype", ARGV_CONSTANT, (char *) MI_ORIGINAL_TYPE, (char *) &datatype, "Use data type of first file (default)."}, {"-byte", ARGV_CONSTANT, (char *) NC_BYTE, (char *) &datatype, "Write out byte data."}, {"-short", ARGV_CONSTANT, (char *) NC_SHORT, (char *) &datatype, "Write out short integer data."}, {"-int", ARGV_CONSTANT, (char *) NC_INT, (char *) &datatype, "Write out 32-bit integer data."}, {"-long", ARGV_CONSTANT, (char *) NC_INT, (char *) &datatype, "Superseded by -int."}, {"-float", ARGV_CONSTANT, (char *) NC_FLOAT, (char *) &datatype, "Write out single-precision floating-point data."}, {"-double", ARGV_CONSTANT, (char *) NC_DOUBLE, (char *) &datatype, "Write out double-precision floating-point data."}, {"-signed", ARGV_CONSTANT, (char *) TRUE, (char *) &is_signed, "Write signed integer data."}, {"-unsigned", ARGV_CONSTANT, (char *) FALSE, (char *) &is_signed, "Write unsigned integer data (default if type specified)."}, {"-range", ARGV_FLOAT, (char *) 2, (char *) valid_range, "Valid range for output data."}, {"-max_buffer_size_in_kb", ARGV_INT, (char *) 1, (char *) &max_buffer_size_in_kb, "Specify the maximum size of the internal buffers (in kbytes)."}, {"-dimension", ARGV_STRING, (char *) 1, (char *) &loop_dimension, "Specify a dimension along which we wish to perform a calculation."}, {"-check_dimensions", ARGV_CONSTANT, (char *) TRUE, (char *) &check_dim_info, "Check that files have matching dimensions (default)."}, {"-nocheck_dimensions", ARGV_CONSTANT, (char *) FALSE, (char *) &check_dim_info, "Do not check that files have matching dimensions."}, {"-ignore_nan", ARGV_CONSTANT, (char *) FALSE, (char *) &propagate_nan, "Ignore invalid data (NaN) for accumulations."}, {"-propagate_nan", ARGV_CONSTANT, (char *) TRUE, (char *) &propagate_nan, "Invalid data in any file at a voxel produces a NaN (default)."}, {"-nan", ARGV_CONSTANT, (char *) TRUE, (char *) &use_nan_for_illegal_values, "Output NaN when an illegal operation is done (default)."}, {"-zero", ARGV_CONSTANT, (char *) FALSE, (char *) &use_nan_for_illegal_values, "Output zero when an illegal operation is done."}, {"-illegal_value", ARGV_FLOAT, (char *) 1, (char *) &value_for_illegal_operations, "Value to write out when an illegal operation is done."}, {NULL, ARGV_HELP, (char *) NULL, (char *) NULL, "Options for specifying constants:"}, {"-constant", ARGV_FLOAT, (char *) 1, (char *) &constant, "Specify a constant argument."}, {"-const", ARGV_FLOAT, (char *) 1, (char *) &constant, "Synonym for -constant."}, {"-const2", ARGV_FLOAT, (char *) 2, (char *) constant2, "Specify two constant arguments."}, {NULL, ARGV_HELP, (char *) NULL, (char *) NULL, "Operations:"}, {"-add", ARGV_CONSTANT, (char *) ADD_OP, (char *) &operation, "Add N volumes or volume + constant."}, {"-sub", ARGV_CONSTANT, (char *) SUB_OP, (char *) &operation, "Subtract 2 volumes or volume - constant."}, {"-mult", ARGV_CONSTANT, (char *) MULT_OP, (char *) &operation, "Multiply N volumes or volume * constant."}, {"-div", ARGV_CONSTANT, (char *) DIV_OP, (char *) &operation, "Divide 2 volumes or volume / constant."}, {"-invert", ARGV_CONSTANT, (char *) INVERT_OP, (char *) &operation, "Calculate 1/x at each voxel (use -constant for c/x)."}, {"-sqrt", ARGV_CONSTANT, (char *) SQRT_OP, (char *) &operation, "Take square root of a volume."}, {"-square", ARGV_CONSTANT, (char *) SQUARE_OP, (char *) &operation, "Take square of a volume."}, {"-abs", ARGV_CONSTANT, (char *) ABS_OP, (char *) &operation, "Take absolute value of a volume."}, {"-max", ARGV_CONSTANT, (char *) MAX_OP, (char *) &operation, "Synonym for -maximum."}, {"-maximum", ARGV_CONSTANT, (char *) MAX_OP, (char *) &operation, "Find maximum of N volumes."}, {"-minimum", ARGV_CONSTANT, (char *) MIN_OP, (char *) &operation, "Find minimum of N volumes."}, {"-exp", ARGV_CONSTANT, (char *) EXP_OP, (char *) &operation, "Calculate c2*exp(c1*x). The constants c1 and c2 default to 1."}, {"-log", ARGV_CONSTANT, (char *) LOG_OP, (char *) &operation, "Calculate log(x/c2)/c1. The constants c1 and c2 default to 1."}, {"-scale", ARGV_CONSTANT, (char *) SCALE_OP, (char *) &operation, "Scale a volume: volume * c1 + c2."}, {"-clamp", ARGV_CONSTANT, (char *) CLAMP_OP, (char *) &operation, "Clamp a volume to lie between two values."}, {"-segment", ARGV_CONSTANT, (char *) SEGMENT_OP, (char *) &operation, "Segment a volume using range of -const2: within range = 1, outside range = 0."}, {"-nsegment", ARGV_CONSTANT, (char *) NSEGMENT_OP, (char *) &operation, "Opposite of -segment: within range = 0, outside range = 1."}, {"-percentdiff", ARGV_CONSTANT, (char *) PERCENTDIFF_OP, (char *) &operation, "Percent difference between 2 volumes, thresholded (const def=0.0)."}, {"-pd", ARGV_CONSTANT, (char *) PERCENTDIFF_OP, (char *) &operation, "Synonym for -percentdiff."}, {"-eq", ARGV_CONSTANT, (char *) EQ_OP, (char *) &operation, "Test for integer vol1 == vol2 or vol1 == const."}, {"-ne", ARGV_CONSTANT, (char *) NE_OP, (char *) &operation, "Test for integer vol1 != vol2 or vol1 != const."}, {"-gt", ARGV_CONSTANT, (char *) GT_OP, (char *) &operation, "Test for vol1 > vol2 or vol1 > const."}, {"-ge", ARGV_CONSTANT, (char *) GE_OP, (char *) &operation, "Test for vol1 >= vol2 or vol1 >= const."}, {"-lt", ARGV_CONSTANT, (char *) LT_OP, (char *) &operation, "Test for vol1 < vol2 or vol1 < const."}, {"-le", ARGV_CONSTANT, (char *) LE_OP, (char *) &operation, "Test for vol1 <= vol2 or vol1 <= const."}, {"-and", ARGV_CONSTANT, (char *) AND_OP, (char *) &operation, "Calculate vol1 && vol2 (&& ...)."}, {"-or", ARGV_CONSTANT, (char *) OR_OP, (char *) &operation, "Calculate vol1 || vol2 (|| ...)."}, {"-not", ARGV_CONSTANT, (char *) NOT_OP, (char *) &operation, "Calculate !vol1."}, {"-isnan", ARGV_CONSTANT, (char *) ISNAN_OP, (char *) &operation, "Test for NaN values in vol1."}, {"-nisnan", ARGV_CONSTANT, (char *) NISNAN_OP, (char *) &operation, "Negation of -isnan."}, {"-count_valid", ARGV_CONSTANT, (char *) COUNT_OP, (char *) &operation, "Count the number of valid values in N volumes."}, {NULL, ARGV_END, NULL, NULL, NULL} }; /* Main program */ int main(int argc, char *argv[]) { char **infiles, **outfiles; int nfiles, nout; char *arg_string; Math_Data math_data; Loop_Options *loop_options; char *pname; int num_constants; Num_Operands num_operands; VoxelFunction math_function; /* Save time stamp and args */ arg_string = time_stamp(argc, argv); /* Get arguments */ pname = argv[0]; if (ParseArgv(&argc, argv, argTable, 0) || (argc < 2)) { (void) fprintf(stderr, "\nUsage: %s [options] [ ...] \n", pname); (void) fprintf(stderr, " %s -help\n\n", pname); exit(EXIT_FAILURE); } nout = 1; outfiles = &argv[argc-1]; /* Get the list of input files either from the command line or from a file, or report an error if both are specified */ nfiles = argc - 2; if (filelist == NULL) { infiles = &argv[1]; } else if (nfiles <= 0) { infiles = read_file_names(filelist, &nfiles); if (infiles == NULL) { (void) fprintf(stderr, "Error reading in file names from file \"%s\"\n", filelist); exit(EXIT_FAILURE); } } else { (void) fprintf(stderr, "Do not specify both -filelist and input file names\n"); exit(EXIT_FAILURE); } /* Make sure that we have something to process */ if (nfiles == 0) { (void) fprintf(stderr, "No input files specified\n"); exit(EXIT_FAILURE); } /* Handle special case of COUNT_OP - it always assume -ignore_nan and -zero */ if (operation == COUNT_OP) { propagate_nan = FALSE; value_for_illegal_operations = 0.0; } /* Check that the arguments make sense */ if ((constant != DEFAULT_DBL) && (constant2[0] != DEFAULT_DBL)) { (void) fprintf(stderr, "%s: Specify only one of -constant or -const2\n", pname); exit(EXIT_FAILURE); } if (constant != DEFAULT_DBL) num_constants = 1; else if (constant2[0] != DEFAULT_DBL) num_constants = 2; else num_constants = 0; num_operands = OperandTable[operation][num_constants]; if (num_operands == ILLEGAL_NUMOP) { (void) fprintf(stderr, "%s: Operation and constants do not match.\n", pname); exit(EXIT_FAILURE); } if ((num_operands != NARY_NUMOP) && (loop_dimension != NULL)) { (void) fprintf(stderr, "%s: Use -dimension only for cumulative ops.\n", pname); exit(EXIT_FAILURE); } if ((num_operands == UNARY_NUMOP) && (nfiles != 1)) { (void) fprintf(stderr, "%s: Expected only one input file.\n", pname); exit(EXIT_FAILURE); } if ((num_operands == BINARY_NUMOP) && (nfiles != 2)) { (void) fprintf(stderr, "%s: Expected two input files.\n", pname); exit(EXIT_FAILURE); } if ((num_operands == NARY_NUMOP) && (nfiles < 1) && (loop_dimension == NULL)) { (void) fprintf(stderr, "%s: Expected at least one input files.\n", pname); exit(EXIT_FAILURE); } /* Set default copy_all_header according to number of input files */ if (copy_all_header == DEFAULT_BOOL) copy_all_header = (nfiles == 1); /* Set up math data structure */ math_data.operation = operation; math_data.num_operands = num_operands; math_data.propagate_nan = propagate_nan; math_data.num_constants = num_constants; switch (num_constants) { case 1: math_data.constants[0] = constant; break; case 2: math_data.constants[0] = constant2[0]; math_data.constants[1] = constant2[1]; break; } if (value_for_illegal_operations != DEFAULT_DBL) math_data.illegal_value = value_for_illegal_operations; else if (use_nan_for_illegal_values) math_data.illegal_value = INVALID_DATA; else math_data.illegal_value = 0.0; /* Do math */ loop_options = create_loop_options(); set_loop_verbose(loop_options, verbose); set_loop_clobber(loop_options, clobber); #ifdef MINC2 set_loop_v2format(loop_options, minc2_format); #endif /* MINC2 defined */ set_loop_datatype(loop_options, datatype, is_signed, valid_range[0], valid_range[1]); if (num_operands == NARY_NUMOP) { math_function = accum_math; set_loop_accumulate(loop_options, TRUE, 0, start_math, end_math); } else { math_function = do_math; } set_loop_copy_all_header(loop_options, copy_all_header); set_loop_dimension(loop_options, loop_dimension); set_loop_buffer_size(loop_options, (long) 1024 * max_buffer_size_in_kb); set_loop_check_dim_info(loop_options, check_dim_info); voxel_loop(nfiles, infiles, nout, outfiles, arg_string, loop_options, math_function, (void *) &math_data); free_loop_options(loop_options); exit(EXIT_SUCCESS); } /* ----------------------------- MNI Header ----------------------------------- @NAME : do_math @INPUT : Standard for voxel loop @OUTPUT : Standard for voxel loop @RETURNS : (nothing) @DESCRIPTION: Routine doing math operations. @METHOD : @GLOBALS : @CALLS : @CREATED : April 25, 1995 (Peter Neelin) @MODIFIED : ---------------------------------------------------------------------------- */ static void 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) /* ARGSUSED */ { Math_Data *math_data; long ivox; double value1, value2; double illegal_value; Operation operation; int num_constants, iconst; double constants[2]; /* Get pointer to window info */ math_data = (Math_Data *) caller_data; /* Check arguments */ if ((input_num_buffers > 2) || (output_num_buffers != 1) || (output_vector_length != input_vector_length)) { (void) fprintf(stderr, "Bad arguments to do_math!\n"); exit(EXIT_FAILURE); } /* Get info */ operation = math_data->operation; num_constants = math_data->num_constants; for (iconst=0; iconst < sizeof(constants)/sizeof(constants[0]); iconst++) { if (iconst < num_constants) constants[iconst] = math_data->constants[iconst]; else if ((operation == INVERT_OP) || (operation == EXP_OP) || (operation == LOG_OP)) constants[iconst] = 1.0; else constants[iconst] = 0.0; } illegal_value = math_data->illegal_value; /* Set default second value */ value2 = constants[0]; /* Loop through the voxels */ for (ivox=0; ivox < num_voxels*input_vector_length; ivox++) { value1 = input_data[0][ivox]; if (input_num_buffers == 2) value2 = input_data[1][ivox]; if ((value1 == INVALID_DATA) || (value2 == INVALID_DATA)) { switch(operation) { case ISNAN_OP: output_data[0][ivox] = 1.0; break; case NISNAN_OP: output_data[0][ivox] = 0.0; break; default: output_data[0][ivox] = INVALID_DATA; break; } } else { switch (operation) { case ADD_OP: output_data[0][ivox] = value1 + value2; break; case SUB_OP: output_data[0][ivox] = value1 - value2; break; case MULT_OP: output_data[0][ivox] = value1 * value2; break; case DIV_OP: if (value2 != 0.0) output_data[0][ivox] = value1 / value2; else output_data[0][ivox] = illegal_value; break; case INVERT_OP: if (value1 == 0.0) output_data[0][ivox] = illegal_value; else output_data[0][ivox] = value2 / value1; break; case SQRT_OP: if (value1 < 0.0) output_data[0][ivox] = illegal_value; else output_data[0][ivox] = sqrt(value1); break; case SQUARE_OP: output_data[0][ivox] = value1 * value1; break; case ABS_OP: if (value1 < 0.0) output_data[0][ivox] = -value1; else output_data[0][ivox] = value1; break; case EXP_OP: output_data[0][ivox] = constants[1] * exp(value1 * constants[0]); break; case LOG_OP: if ((value1 <= 0.0) || (constants[1] <= 0.0) || (constants[0] == 0.0)) output_data[0][ivox] = illegal_value; else output_data[0][ivox] = log(value1/constants[1])/constants[0]; break; case SCALE_OP: output_data[0][ivox] = value1 * constants[0] + constants[1]; break; case CLAMP_OP: if (value1 < constants[0]) value1 = constants[0]; else if (value1 > constants[1]) value1 = constants[1]; output_data[0][ivox] = value1; break; case SEGMENT_OP: if ((value1 < constants[0]) || (value1 > constants[1])) output_data[0][ivox] = 0.0; else output_data[0][ivox] = 1.0; break; case NSEGMENT_OP: if ((value1 < constants[0]) || (value1 > constants[1])) output_data[0][ivox] = 1.0; else output_data[0][ivox] = 0.0; break; case PERCENTDIFF_OP: if ((value1 < constants[0]) || (value1 == 0.0)) output_data[0][ivox] = illegal_value; else { output_data[0][ivox] = 100.0 * (value1 - value2) / value1; } break; case EQ_OP: output_data[0][ivox] = (((rint(value1)-rint(value2)) == 0.0) ? 1.0 : 0.0); break; case NE_OP: output_data[0][ivox] = (((rint(value1)-rint(value2)) != 0.0) ? 1.0 : 0.0); break; case GT_OP: output_data[0][ivox] = value1 > value2; break; case GE_OP: output_data[0][ivox] = value1 >= value2; break; case LT_OP: output_data[0][ivox] = value1 < value2; break; case LE_OP: output_data[0][ivox] = value1 <= value2; break; case AND_OP: output_data[0][ivox] = (((rint(value1) != 0.0) && (rint(value2) != 0.0)) ? 1.0 : 0.0); break; case OR_OP: output_data[0][ivox] = (((rint(value1) != 0.0) || (rint(value2) != 0.0)) ? 1.0 : 0.0); break; case NOT_OP: output_data[0][ivox] = ((rint(value1) == 0.0) ? 1.0 : 0.0); break; case ISNAN_OP: output_data[0][ivox] = 0.0; /* To get here, value is not nan */ break; case NISNAN_OP: output_data[0][ivox] = 1.0; break; default: (void) fprintf(stderr, "Bad op in do_math!\n"); exit(EXIT_FAILURE); } } } return; } /* ----------------------------- MNI Header ----------------------------------- @NAME : accum_math @INPUT : Standard for voxel loop @OUTPUT : Standard for voxel loop @RETURNS : (nothing) @DESCRIPTION: Routine for doing accumulation math operations. @METHOD : @GLOBALS : @CALLS : @CREATED : April 25, 1995 (Peter Neelin) @MODIFIED : ---------------------------------------------------------------------------- */ static void 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) /* ARGSUSED */ { Math_Data *math_data; long ivox; double value, oldvalue; Operation operation; int propagate_nan; /* Get pointer to window info */ math_data = (Math_Data *) caller_data; /* Check arguments */ if ((input_num_buffers != 1) || (output_num_buffers != 1) || (output_vector_length != input_vector_length)) { (void) fprintf(stderr, "Bad arguments to accum_math!\n"); exit(EXIT_FAILURE); } /* Get info */ operation = math_data->operation; propagate_nan = math_data->propagate_nan; /* Loop through the voxels */ for (ivox=0; ivox < num_voxels*input_vector_length; ivox++) { /* Get previous value and the next value */ oldvalue = output_data[0][ivox]; value = input_data[0][ivox]; /* If the new data is invalid, then either mark the output as invalid or ignore it */ if (value == INVALID_DATA) { if (propagate_nan) { output_data[0][ivox] = INVALID_DATA; } } /* If we haven't set anything yet, then just copy the new value */ else if (oldvalue == UNINITIALIZED_DATA) { output_data[0][ivox] = input_data[0][ivox]; } /* Do the operation if the old data and the new data are valid */ else if (oldvalue != INVALID_DATA) { switch (operation) { case ADD_OP: output_data[0][ivox] = oldvalue + value; break; case MULT_OP: output_data[0][ivox] = oldvalue * value; break; case AND_OP: output_data[0][ivox] = (((oldvalue != 0.0) && (rint(value) != 0.0)) ? 1.0 : 0.0); break; case OR_OP: output_data[0][ivox] = (((oldvalue != 0.0) || (rint(value) != 0.0)) ? 1.0 : 0.0); break; case MAX_OP: if (value > oldvalue) output_data[0][ivox] = value; break; case MIN_OP: if (value < oldvalue) output_data[0][ivox] = value; break; case COUNT_OP: output_data[0][ivox]++; break; default: (void) fprintf(stderr, "Bad op in accum_math!\n"); exit(EXIT_FAILURE); } } } /* Loop over voxels */ return; } /* ----------------------------- MNI Header ----------------------------------- @NAME : start_math @INPUT : Standard for voxel loop @OUTPUT : Standard for voxel loop @RETURNS : (nothing) @DESCRIPTION: Start routine for math accumulation. @METHOD : @GLOBALS : @CALLS : @CREATED : April 25, 1995 (Peter Neelin) @MODIFIED : ---------------------------------------------------------------------------- */ static void start_math(void *caller_data, long num_voxels, int output_num_buffers, int output_vector_length, double *output_data[], Loop_Info *loop_info) /* ARGSUSED */ { Math_Data *math_data; long ivox; /* Get pointer to window info */ math_data = (Math_Data *) caller_data; /* Check arguments */ if (output_num_buffers != 1) { (void) fprintf(stderr, "Bad arguments to start_math!\n"); exit(EXIT_FAILURE); } /* Get info */ operation = math_data->operation; /* Loop through the voxels, marking them all as uninitialized. We treat COUNT_OP as a special case since it always has a value. This is especially important to prevent it from going through the code in accum_math for handling the first valid voxel which just assigns the first value. */ for (ivox=0; ivox < num_voxels*output_vector_length; ivox++) { switch (operation) { case COUNT_OP: output_data[0][ivox] = 0.0; break; default: output_data[0][ivox] = UNINITIALIZED_DATA; break; } } return; } /* ----------------------------- MNI Header ----------------------------------- @NAME : end_math @INPUT : Standard for voxel loop @OUTPUT : Standard for voxel loop @RETURNS : (nothing) @DESCRIPTION: Start routine for math accumulation. @METHOD : @GLOBALS : @CALLS : @CREATED : April 25, 1995 (Peter Neelin) @MODIFIED : ---------------------------------------------------------------------------- */ static void end_math(void *caller_data, long num_voxels, int output_num_buffers, int output_vector_length, double *output_data[], Loop_Info *loop_info) /* ARGSUSED */ { Math_Data *math_data; long ivox; double value; double illegal_value; /* Get pointer to window info */ math_data = (Math_Data *) caller_data; /* Check arguments */ if (output_num_buffers != 1) { (void) fprintf(stderr, "Bad arguments to end_math!\n"); exit(EXIT_FAILURE); } /* Get info */ operation = math_data->operation; illegal_value = math_data->illegal_value; /* Loop through the voxels, checking for uninitialized values */ for (ivox=0; ivox < num_voxels*output_vector_length; ivox++) { value = output_data[0][ivox]; if ((value == UNINITIALIZED_DATA) || (value == INVALID_DATA)) { output_data[0][ivox] = illegal_value; } } return; } /* ----------------------------- MNI Header ----------------------------------- @NAME : read_file_names @INPUT : filelist - name of file from which to read names @OUTPUT : num_files - number of files read in @RETURNS : Pointer to a NULL-terminated array of file names @DESCRIPTION: Reads in a list of file names from file filelist or stdin if "-" is specified. Returns NULL if an error occurs. If no error occurs, then a pointer to an empty array is returned and num_files is zero. @METHOD : @GLOBALS : @CALLS : @CREATED : March 8, 1995 (Peter Neelin) @MODIFIED : ---------------------------------------------------------------------------- */ static char **read_file_names(char *filelist, int *num_files) { #define FILE_NAME_ALLOC_SIZE 10 char **files; int array_size; int nfiles; FILE *fp; char line[PATH_MAX+1]; int length; /* Open the file */ if (strcmp(filelist, "-") == 0) { fp = stdin; } else { fp = fopen(filelist, "r"); if (fp == NULL) { (void) fprintf(stderr, "Error opening file \"%s\"\n", filelist); return NULL; } } /* Allocate an initial array and NULL-terminate it */ array_size = FILE_NAME_ALLOC_SIZE; files = malloc(sizeof(*files) * array_size); if (files == NULL) { (void) fprintf(stderr, "Error allocating memory\n"); return NULL; } nfiles = 0; files[nfiles] = NULL; /* Read in file names */ while (fgets(line, sizeof(line)/sizeof(line[0]), fp) != NULL) { /* Remove a trailing newline and check that there is a name */ length = strlen(line); if ((length > 0) && (line[length-1] == '\n')) { line[length-1] = '\0'; length--; } if (length == 0) continue; /* Make room for names if needed */ while (nfiles >= array_size-1) { array_size += FILE_NAME_ALLOC_SIZE; files = realloc(files, sizeof(*files) * array_size); if (files == NULL) { (void) fprintf(stderr, "Error allocating memory\n"); return NULL; } } /* Save the name, making sure that the list is NULL-terminated */ files[nfiles] = strdup(line); if (files[nfiles] == NULL) { (void) fprintf(stderr, "Error allocating memory\n"); return NULL; } nfiles++; files[nfiles] = NULL; } /* Close the file */ (void) fclose(fp); /* Return the number of files */ *num_files = nfiles; return files; }