/* ----------------------------- MNI Header -----------------------------------
@NAME       : mincresample
@INPUT      : argc, argv - command line arguments
@OUTPUT     : (none)
@RETURNS    : error status
@DESCRIPTION: Program to resample a minc file along different spatial
              coordinate axes.
@METHOD     : 
@GLOBALS    : 
@CALLS      : 
@CREATED    : February 8, 1993 (Peter Neelin)
@MODIFIED   : 
 * $Log: mincresample.c,v $
 * Revision 6.23  2008-01-17 02:33:02  rotor
 *  * removed all rcsids
 *  * removed a bunch of ^L's that somehow crept in
 *  * removed old (and outdated) BUGS file
 *
 * Revision 6.22  2008/01/13 09:38:54  stever
 * Avoid compiler warnings about functions and variables that are defined
 * but not used.  Remove some such functions and variables,
 * conditionalize some, and move static declarations out of header files
 * into C files.
 *
 * Revision 6.21  2008/01/12 19:08:15  stever
 * Add __attribute__ ((unused)) to all rcsid variables.
 *
 * Revision 6.20  2006/07/28 18:19:46  baghdadi
 * *** empty log message ***
 *
 * Revision 6.19  2006/07/28 17:54:57  baghdadi
 * added vector dimension to list of excluded files (minc2.
 *
 * Revision 6.18  2005/08/26 21:07:17  bert
 * Use #if rather than #ifdef with MINC2 symbol, and be sure to include config.h whereever MINC2 is used
 *
 * Revision 6.17  2005/07/13 21:34:24  bert
 * Add sinc interpolant (ported from 1.X branch)
 *
 * Revision 6.16  2004/11/01 22:38:39  bert
 * Eliminate all references to minc_def.h
 *
 * Revision 6.15  2004/04/30 19:53:40  bert
 * Remove unused variable
 *
 * Revision 6.14  2004/04/30 18:52:49  bert
 * Remove some unused variables
 *
 * Revision 6.13  2004/04/27 15:31:20  bert
 * Added -2 option
 *
 * Revision 6.12  2003/09/18 15:01:33  bert
 * Removed unnecessary brackets from initializer
 *
 * Revision 6.11  2002/11/06 13:32:23  jason
 * Fixes to mincresample: setting the interpolation type is now done
 * through an enum rather than function pointers.
 *
 * Revision 6.10  2002/10/30 13:53:02  jason
 * Added a ARGV_LONG argument type to ParseArgv. Used that type for the nelements variable in mincresample
 *
 * Revision 6.9  2001/08/24 19:12:50  neelin
 * Re-ordered variables so that image variable is last. This fixes a problem
 * with an uninitialized processing variable that shows up when mincresample
 * is linked with netcdf 3.5.0 and an older program linked with 2.3.2 dies
 * when reading the mincresample output file.
 * It also allows use of large volumes with 64-bit minc: all variables
 * must have 32-bit offset, but one (the image) can extend beyond the
 * 32-bit limit.
 *
 * Revision 6.8  2001/08/16 13:32:39  neelin
 * Partial fix for valid_range of different type from image (problems
 * arising from double to float conversion/rounding). NOT COMPLETE.
 *
 * Revision 6.7  2001/04/24 13:38:45  neelin
 * Replaced NC_NAT with MI_ORIGINAL_TYPE.
 *
 * Revision 6.6  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.5  2001/04/10 12:44:34  neelin
 * Re-introduced fix for fillvalues that was lost with last version
 *
 * Revision 6.4  2001/04/02 14:58:09  neelin
 * Added -keep_real_range option to prevent rescaling of slices on output
 *
 * Revision 6.2  1999/10/19 14:45:27  neelin
 * Fixed Log subsitutions for CVS
 *
 * Revision 6.1  1998/02/19 15:04:24  neelin
 * Minor bug fixes.
 *
 * Revision 6.0  1997/09/12  13:23:21  neelin
 * Release of minc version 0.6
 *
 * Revision 5.0  1997/08/21  13:24:22  neelin
 * Release of minc version 0.5
 *
 * Revision 4.1  1997/08/13  15:41:12  neelin
 * Fixed initialization problem that caused crashing under Linux.
 *
 * Revision 4.0  1997/05/07  19:59:42  neelin
 * Release of minc version 0.4
 *
 * Revision 3.5  1996/01/31  15:22:02  neelin
 * Fixed bug in transformation of input sampling.
 *
 * Revision 3.4  1996/01/30  14:10:24  neelin
 * Added check for -transformation without -like, -tfm_input_sampling or
 * -use_input_sampling because of change in behaviour.
 *
 * Revision 3.3  1995/12/12  19:15:35  neelin
 * Added -spacetype, -talairach and -units options.
 *
 * Revision 3.2  1995/11/21  14:13:20  neelin
 * Transform input sampling with transformation and use this as default.
 * Added -tfm_input_sampling to specify above option.
 * Added -use_input_sampling to get old behaviour (no longer the default).
 * Added -origin option (to specify coordinate instead of start values).
 * Added -standard_sampling option (to set standard values of start, step
 * and direction cosines).
 * Added -invert_transformation option.
 *
 * Revision 3.1  1995/11/07  15:04:02  neelin
 * Modified argument parsing so that only one pass is done.
 *
 * Revision 3.0  1995/05/15  19:30:57  neelin
 * Release of minc version 0.3
 *
 * Revision 2.3  1995/05/05  19:11:05  neelin
 * Modified call to input_transform.
 *
 * Revision 2.2  1995/02/09  14:05:51  neelin
 * Mods to make irix 5 lint happy.
 *
 * Revision 2.1  1995/02/08  19:31:47  neelin
 * Moved ARGSUSED statements for irix 5 lint.
 *
 * Revision 2.0  1994/09/28  10:32:46  neelin
 * Release of minc version 0.2
 *
 * Revision 1.16  94/09/28  10:32:33  neelin
 * Pre-release
 * 
 * Revision 1.15  94/03/17  14:12:09  neelin
 * Exit with failure if no argument given for -transformation or -like.
 * .,
 * 
 * Revision 1.14  94/03/15  16:44:21  neelin
 * Changed default from -clobber to -noclobber.
 * 
 * Revision 1.13  93/11/04  15:13:13  neelin
 * Added support for irregularly spaced dimensions.
 * 
 * Revision 1.12  93/11/03  14:32:44  neelin
 * Turn off fill for output file.
 * 
 * Revision 1.11  93/11/03  12:32:17  neelin
 * Change ncopen, nccreate and ncclose to miopen, micreate and miclose.
 * 
 * Revision 1.10  93/11/02  11:23:06  neelin
 * Handle imagemax/min potentially varying over slices (for vector data, etc.)
 * 
 * Revision 1.9  93/10/12  12:47:50  neelin
 * Use volume_io.h instead of def_mni.h
 * 
 * Revision 1.8  93/09/16  09:56:36  neelin
 * Added use of open_file_with_default_suffix in get_transformation to
 * append appropriate suffix for xfm files.
 * 
 * Revision 1.7  93/08/11  14:28:19  neelin
 * Modified get_arginfo and check_imageminmax to modify type of volume (not
 * file) so that output volume gets the input volume type by default when
 * an icv is used on input.
 * 
 * Revision 1.6  93/08/11  13:27:59  neelin
 * Converted to use Dave MacDonald's General_transform code.
 * Fixed bug in get_slice - for non-linear transformations coord was
 * transformed, then used again as a starting coordinate.
 * Handle files that have image-max/min that doesn't vary over slices.
 * Handle files that have image-max/min varying over row/cols.
 * Allow volume to extend to voxel edge for -nearest_neighbour interpolation.
 * Handle out-of-range values (-fill values from a previous mincresample, for
 * example).
 * Save transformation file as a string attribute to processing variable.
 * 
@COPYRIGHT  :
              Copyright 1993 Peter Neelin, McConnell Brain Imaging Centre, 
              Montreal Neurological Institute, McGill University.
              Permission to use, copy, modify, and distribute this
              software and its documentation for any purpose and without
              fee is hereby granted, provided that the above copyright
              notice appear in all copies.  The author and McGill University
              make no representations about the suitability of this
              software for any purpose.  It is provided "as is" without
              express or implied warranty.
---------------------------------------------------------------------------- */

#if HAVE_CONFIG_H
#include "config.h"
#endif

#include <stdlib.h>
#include <stdio.h>
#include <float.h>
#include <limits.h>
#include <string.h>
#include <math.h>
#include <minc.h>
#include <ParseArgv.h>
#include <time_stamp.h>
#include <volume_io.h>
#include <convert_origin_to_start.h>
#include "mincresample.h"

/* Kludge to catch use of -transformation without -like and without 
   -tfm_input_sampling or -use_input_sampling */
#define TRANSFORM_CHANGE_KLUDGE
#ifdef TRANSFORM_CHANGE_KLUDGE
int Specified_like = FALSE;
int Specified_transform = FALSE;
#define SAMPLING_ACTION_NOT_SET (-1)
#endif

static void get_arginfo(int argc, char *argv[],
                        Program_Flags *program_flags,
                        VVolume *in_vol, VVolume *out_vol,
                        General_transform *transformation);
static void check_imageminmax(File_Info *fp, Volume_Data *volume);
static void get_file_info(char *filename, int initialized_volume_def,
                          Volume_Definition *volume_def,
                          File_Info *file_info);
static void get_args_volume_def(Volume_Definition *input_volume_def,
                                Volume_Definition *args_volume_def);
static void transform_volume_def(Transform_Info *transform_info,
                                 Volume_Definition *input_volume_def,
                                 Volume_Definition *transformed_volume_def);
static int is_zero_vector(double vector[]);
static void normalize_vector(double vector[]);
static void create_output_file(char *filename, int clobber,
                               Volume_Definition *volume_def,
                               File_Info *in_file,
                               File_Info *out_file,
                               char *tm_stamp,
                               Transform_Info *transform_info);
static void get_voxel_to_world_transf(Volume_Definition *volume_def,
                                      General_transform *voxel_to_world);
static void irregular_transform_function(void *user_data,
                                         Real x,
                                         Real y,
                                         Real z,
                                         Real *x_trans,
                                         Real *y_trans,
                                         Real *z_trans);
static void irregular_inverse_transform_function(void *user_data,
                                                 Real x,
                                                 Real y,
                                                 Real z,
                                                 Real *x_trans,
                                                 Real *y_trans,
                                                 Real *z_trans);
static double get_default_range(char *what, nc_type datatype, int is_signed);
static void finish_up(VVolume *in_vol, VVolume *out_vol);
static int get_transformation(char *dst, char *key, char *nextArg);
static int get_model_file(char *dst, char *key, char *nextArg);
static int set_standard_sampling(char *dst, char *key, char *nextArg);
static int set_spacetype(char *dst, char *key, char *nextArg);
static int set_units(char *dst, char *key, char *nextArg);
static int get_axis_order(char *dst, char *key, char *nextArg);
static int get_fillvalue(char *dst, char *key, char *nextArg);

/* Main program */

int main(int argc, char *argv[])
{
   VVolume in_vol_struct, out_vol_struct;
   VVolume *in_vol = &in_vol_struct, *out_vol = &out_vol_struct;
   General_transform transformation;
   Program_Flags program_flags;

   /* Get argument information */
   get_arginfo(argc, argv, &program_flags, in_vol, out_vol, &transformation);

   /* Do the resampling */
   resample_volumes(&program_flags, in_vol, out_vol, &transformation);

   /* Finish up */
   finish_up(in_vol, out_vol);

   exit(EXIT_SUCCESS);
}

/* ----------------------------- MNI Header -----------------------------------
@NAME       : get_arginfo
@INPUT      : argc - number of command-line arguments
              argv - command-line arguments
@OUTPUT     : program_flags - data for program execution
              in_vol - description of input volume.
              out_vol - description of output volume.
              transformation - description of world transformation
@RETURNS    : (nothing)
@DESCRIPTION: Routine to get information from arguments about input and 
              output files and transfomation. Sets up all structures
              completely (including allocating space for data).
@METHOD     : 
@GLOBALS    : 
@CALLS      : 
@CREATED    : February 8, 1993 (Peter Neelin)
@MODIFIED   : 
---------------------------------------------------------------------------- */
static void get_arginfo(int argc, char *argv[],
                        Program_Flags *program_flags,
                        VVolume *in_vol, VVolume *out_vol, 
                        General_transform *transformation)
{
   /* Argument parsing information */
#ifdef TRANSFORM_CHANGE_KLUDGE
   static int transform_input_sampling = SAMPLING_ACTION_NOT_SET;
#else
   static int transform_input_sampling = TRUE;
#endif
   static Arg_Data args={
      FALSE,                  /* Clobber */
      FALSE,                  /* Keep scale */
      MI_ORIGINAL_TYPE,         /* Flag that type not set */
      INT_MIN,                /* Flag that is_signed has not been set */
      {-DBL_MAX, -DBL_MAX},   /* Flag that range not set */
      FILL_DEFAULT,           /* Flag indicating that fillvalue not set */
      {NO_VALUE, NO_VALUE, NO_VALUE}, /* Flag indicating that origin not set */
      {TRUE},                 /* Verbose */
      TRILINEAR,              /* use trilinear interpolation by default */
      {FALSE, NULL, NULL, 0, NULL}, /* Transformation info is empty at start.
                                 Transformation must be set before invoking
                                 argument parsing */
      {                       /* Use flags to indicate that values not set */
          {NO_AXIS, NO_AXIS, NO_AXIS},          /* Axis order */
          {0, 0, 0},          /* nelements */
          {0.0, 0.0, 0.0},    /* step */
          {NO_VALUE, NO_VALUE, NO_VALUE},    /* start */
          {{NO_VALUE, NO_VALUE, NO_VALUE},   /* Default direction cosines */
           {NO_VALUE, NO_VALUE, NO_VALUE},
           {NO_VALUE, NO_VALUE, NO_VALUE}},
          {NULL, NULL, NULL}, /* Pointers to coordinate arrays */
          {"", "", ""},       /* units */
          {"", "", ""}        /* spacetype */
      },
      FALSE			/* MINC 2.0 format? */
   };

   static ArgvInfo argTable[] = {
#if MINC2
       {"-2", ARGV_CONSTANT, (char *) TRUE, 
	(char *) &args.v2format,
       "Produce a MINC 2.0 format output file."},
#endif /* MINC2 */
      {"-clobber", ARGV_CONSTANT, (char *) TRUE, 
          (char *) &args.clobber,
          "Overwrite existing file."},
      {"-noclobber", ARGV_CONSTANT, (char *) FALSE, 
          (char *) &args.clobber,
          "Do not overwrite existing file (default)."},
      {"-verbose", ARGV_CONSTANT, (char *) TRUE,
          (char *) &args.flags.verbose,
          "Print out log messages as processing is being done (default).\n"},
      {"-quiet", ARGV_CONSTANT, (char *) FALSE,
          (char *) &args.flags.verbose,
          "Do not print out any log messages.\n"},
      {"-transformation", ARGV_FUNC, (char *) get_transformation, 
          (char *) &args.transform_info,
          "File giving world transformation. (Default = identity)."},
      {"-invert_transformation", ARGV_CONSTANT, (char *) TRUE,
          (char *) &args.transform_info.invert_transform,
          "Invert the transformation before using it.\n"},
      {"-noinvert_transformation", ARGV_CONSTANT, (char *) FALSE,
          (char *) &args.transform_info.invert_transform,
          "Do not invert the transformation (default).\n"},
      {"-tfm_input_sampling", ARGV_CONSTANT, (char *) TRUE,
          (char *) &transform_input_sampling,
          "Transform the input sampling with the transform (default).\n"},
      {"-use_input_sampling", ARGV_CONSTANT, (char *) FALSE,
          (char *) &transform_input_sampling,
          "Use the input sampling without transforming (old behaviour).\n"},
      {"-like", ARGV_FUNC, (char *) get_model_file, 
          (char *) &args.volume_def,
          "Specifies a model file for the resampling."},
      {"-standard_sampling", ARGV_FUNC, (char *) set_standard_sampling, 
          (char *) &args.volume_def,
          "Set the sampling to standard values (step, start and dircos)."},
      {"-spacetype", ARGV_FUNC, (char *) set_spacetype, 
          (char *) &args.volume_def,
          "Set the spacetype attribute to a specified string."},
      {"-talairach", ARGV_FUNC, (char *) set_spacetype, 
          (char *) &args.volume_def,
          "Output is in Talairach space."},
      {"-units", ARGV_FUNC, (char *) set_units, 
          (char *) &args.volume_def,
          "Specify the units of the output sampling."},
      {"-nelements", ARGV_LONG, (char *) 3, 
          (char *) args.volume_def.nelements,
          "Number of elements along each dimension (X, Y, Z)"},
      {"-xnelements", ARGV_LONG, (char *) 0, 
          (char *) &args.volume_def.nelements[X],
          "Number of elements along the X dimension"},
      {"-ynelements", ARGV_LONG, (char *) 0, 
          (char *) &args.volume_def.nelements[Y],
          "Number of elements along the Y dimension"},
      {"-znelements", ARGV_LONG, (char *) 0, 
          (char *) &args.volume_def.nelements[Z],
          "Number of elements along the Z dimension"},
      {"-step", ARGV_FLOAT, (char *) 3, 
          (char *) args.volume_def.step,
          "Step size along each dimension (X, Y, Z)"},
      {"-xstep", ARGV_FLOAT, (char *) 0, 
          (char *) &args.volume_def.step[X],
          "Step size along the X dimension"},
      {"-ystep", ARGV_FLOAT, (char *) 0, 
          (char *) &args.volume_def.step[Y],
          "Step size along the Y dimension"},
      {"-zstep", ARGV_FLOAT, (char *) 0, 
          (char *) &args.volume_def.step[Z],
          "Step size along the Z dimension"},
      {"-start", ARGV_FLOAT, (char *) 3, 
          (char *) args.volume_def.start,
          "Start point along each dimension (X, Y, Z)"},
      {"-xstart", ARGV_FLOAT, (char *) 0, 
          (char *) &args.volume_def.start[X],
          "Start point along the X dimension"},
      {"-ystart", ARGV_FLOAT, (char *) 0, 
          (char *) &args.volume_def.start[Y],
          "Start point along the Y dimension"},
      {"-zstart", ARGV_FLOAT, (char *) 0, 
          (char *) &args.volume_def.start[Z],
          "Start point along the Z dimension"},
      {"-dircos", ARGV_FLOAT, (char *) 9, 
          (char *) args.volume_def.dircos,
          "Direction cosines along each dimension (X, Y, Z)"},
      {"-xdircos", ARGV_FLOAT, (char *) 3, 
          (char *) args.volume_def.dircos[X],
          "Direction cosines along the X dimension"},
      {"-ydircos", ARGV_FLOAT, (char *) 3, 
          (char *) args.volume_def.dircos[Y],
          "Direction cosines along the Y dimension"},
      {"-zdircos", ARGV_FLOAT, (char *) 3, 
          (char *) args.volume_def.dircos[Z],
          "Direction cosines along the Z dimension"},
      {"-origin", ARGV_FLOAT, (char *) 3, (char *) args.origin,
          "Origin of first pixel in 3D space"},
      {"-transverse", ARGV_FUNC, (char *) get_axis_order, 
          (char *) &args.volume_def,
          "Write out transverse slices"},
      {"-sagittal", ARGV_FUNC, (char *) get_axis_order, 
          (char *) &args.volume_def,
          "Write out sagittal slices"},
      {"-coronal", ARGV_FUNC, (char *) get_axis_order, 
          (char *) &args.volume_def,
          "Write out coronal slices"},
      {"-byte", ARGV_CONSTANT, (char *) NC_BYTE, (char *) &args.datatype,
          "Write out byte data"},
      {"-short", ARGV_CONSTANT, (char *) NC_SHORT, (char *) &args.datatype,
          "Write out short integer data"},
      {"-int", ARGV_CONSTANT, (char *) NC_INT, (char *) &args.datatype,
          "Write out 32-bit integer data"},
      {"-long", ARGV_CONSTANT, (char *) NC_INT, (char *) &args.datatype,
          "Superseded by -int"},
      {"-float", ARGV_CONSTANT, (char *) NC_FLOAT, (char *) &args.datatype,
          "Write out single-precision floating-point data"},
      {"-double", ARGV_CONSTANT, (char *) NC_DOUBLE, (char *) &args.datatype,
          "Write out double-precision floating-point data"},
      {"-signed", ARGV_CONSTANT, (char *) TRUE, (char *) &args.is_signed,
          "Write signed integer data"},
      {"-unsigned", ARGV_CONSTANT, (char *) FALSE, (char *) &args.is_signed,
          "Write unsigned integer data"},
      {"-range", ARGV_FLOAT, (char *) 2, (char *) args.vrange,
          "Valid range for output data"},
      {"-keep_real_range", ARGV_CONSTANT, (char *) TRUE, 
          (char *) &args.keep_real_range,
          "Keep the real scale of the input volume"},
      {"-nokeep_real_range", ARGV_CONSTANT, (char *) FALSE, 
          (char *) &args.keep_real_range,
          "Do not keep the real scale of the data (default)"},
      {"-nofill", ARGV_FUNC, (char *) get_fillvalue, 
          (char *) &args.fillvalue,
          "Use value zero for points outside of input volume"},
      {"-fill", ARGV_FUNC, (char *) get_fillvalue, 
          (char *) &args.fillvalue,
          "Use a fill value for points outside of input volume"},
      {"-fillvalue", ARGV_FLOAT, (char *) 0, 
          (char *) &args.fillvalue,
          "Specify a fill value for points outside of input volume"},
      {"-trilinear", ARGV_CONSTANT, (char *) TRILINEAR, 
          (char *) &args.interpolant_type,
          "Do trilinear interpolation"},
      {"-tricubic", ARGV_CONSTANT, (char *) TRICUBIC, 
          (char *) &args.interpolant_type,
          "Do tricubic interpolation"},
      {"-nearest_neighbour", ARGV_CONSTANT, 
          (char *) N_NEIGHBOUR, 
          (char *) &args.interpolant_type,
          "Do nearest neighbour interpolation"},
      {"-sinc", ARGV_CONSTANT,
       (char *) WINDOWED_SINC,
       (char *) &args.interpolant_type,
       "Do windowed sinc interpolation"},
      {"-width", ARGV_INT,
       (char *) 1,
       (char *) &sinc_half_width,
      "Set half-width of sinc window (1-10)" },
      {"-hanning", ARGV_CONSTANT,
       (char *) SINC_WINDOW_HANNING,
       (char *) &sinc_window_type,
       "Set sinc window type to Hanning"},
      {"-hamming", ARGV_CONSTANT,
       (char *) SINC_WINDOW_HAMMING,
       (char *) &sinc_window_type,
       "Set sinc window type to Hamming"},
      {NULL, ARGV_END, NULL, NULL, NULL}
   };

   /* Other variables */
   int idim, index;
   int out_vindex;              /* Volume indices (0, 1 or 2) */
   int out_findex;              /* File indices (0 to ndims-1) */
   long size, total_size;
   char *infile, *outfile;
   File_Info *fp;
   char *tm_stamp, *pname;
   Volume_Definition input_volume_def, transformed_volume_def;
   General_transform input_transformation;
   int cflags;

   /* Initialize the transformation to identity */
   create_linear_transform(&input_transformation, NULL);
   args.transform_info.transformation = &input_transformation;

   /* Get the time stamp */
   tm_stamp = time_stamp(argc, argv);

   /* Save the program name */
   pname=argv[0];

   /* Call ParseArgv */
   if (ParseArgv(&argc, argv, argTable, 0) || (argc!=3)) {
      (void) fprintf(stderr, 
                     "\nUsage: %s [<options>] <infile> <outfile>\n", pname);
      (void) fprintf(stderr,   
                     "       %s [-help]\n\n", pname);
      exit(EXIT_FAILURE);
   }
   infile = argv[1];
   outfile = argv[2];

#ifdef TRANSFORM_CHANGE_KLUDGE
   if (Specified_transform && 
       !Specified_like &&
       (transform_input_sampling == SAMPLING_ACTION_NOT_SET)) {
      (void) fprintf(stderr, 
                     "Use -like, -tfm_input_sampling or "
                     "-use_input_sampling with -transformation\n");
      exit(EXIT_FAILURE);
   }
   if (transform_input_sampling == SAMPLING_ACTION_NOT_SET) {
      transform_input_sampling = TRUE;
   }
#endif

   /* Check for an inverted transform. This looks backwards because we 
      normally invert the transform. */
   if (args.transform_info.invert_transform) {
      copy_general_transform(args.transform_info.transformation,
                             transformation);
   }
   else {
      create_inverse_general_transform(args.transform_info.transformation,
                                       transformation);
   }
   args.transform_info.transformation = transformation;

   /* Get rid of the input transformation */
   delete_general_transform(&input_transformation);

   /* Check input file for default argument information */
   in_vol->file = malloc(sizeof(File_Info));
   get_file_info(infile, FALSE, &input_volume_def, in_vol->file);
   transform_volume_def((transform_input_sampling ? 
                         &args.transform_info : NULL), 
                        &input_volume_def, 
                        &transformed_volume_def);
   get_args_volume_def(&transformed_volume_def, &args.volume_def);

   /* Check that direction cosines are normalized and look for origin 
      option */
   for (idim=0; idim < WORLD_NDIMS; idim++) {
      normalize_vector(args.volume_def.dircos[idim]);
      if (is_zero_vector(args.volume_def.dircos[idim])) {
         (void) fprintf(stderr, "Bad direction cosines.\n");
         exit(EXIT_FAILURE);
      }
   }
   if (args.origin[0] != NO_VALUE) {
      if (convert_origin_to_start(args.origin,
                                  args.volume_def.dircos[XCOORD],
                                  args.volume_def.dircos[YCOORD],
                                  args.volume_def.dircos[ZCOORD],
                                  args.volume_def.start) != 0) {
         (void) fprintf(stderr, "Error converting origin to start value: ");
         (void) fprintf(stderr, "Bad direction cosines.\n");
         exit(EXIT_FAILURE);
      }
   }

   /* Save the voxel_to_world transformation information */
   in_vol->voxel_to_world = malloc(sizeof(General_transform));
   in_vol->world_to_voxel = malloc(sizeof(General_transform));
   get_voxel_to_world_transf(&input_volume_def, in_vol->voxel_to_world);
   create_inverse_general_transform(in_vol->voxel_to_world,
                                    in_vol->world_to_voxel);

   /* Get input volume data information */
   in_vol->slice = NULL;
   in_vol->volume = malloc(sizeof(Volume_Data));
   in_vol->volume->datatype = in_vol->file->datatype;
   in_vol->volume->is_signed = in_vol->file->is_signed;
   in_vol->volume->vrange[0] = in_vol->file->vrange[0];
   in_vol->volume->vrange[1] = in_vol->file->vrange[1];
   if (args.fillvalue == FILL_DEFAULT) {
      in_vol->volume->fillvalue = 0.0;
      in_vol->volume->use_fill = TRUE;
   }
   else {
      in_vol->volume->fillvalue = args.fillvalue;
      in_vol->volume->use_fill = (args.fillvalue != -DBL_MAX);
   }

   /* set the function pointer defining the type of interpolation */
   switch (args.interpolant_type ) {
   case TRICUBIC:
     in_vol->volume->interpolant = tricubic_interpolant;
     break;
   case TRILINEAR:
     in_vol->volume->interpolant = trilinear_interpolant;
     break;
   case N_NEIGHBOUR:
     in_vol->volume->interpolant = nearest_neighbour_interpolant;
     break;
   case WINDOWED_SINC:
     in_vol->volume->interpolant = windowed_sinc_interpolant;
     if (sinc_half_width < SINC_HALF_WIDTH_MIN ||
         sinc_half_width > SINC_HALF_WIDTH_MAX) {
         fprintf(stderr, "Invalid sinc half-window size %d\n", 
                 sinc_half_width);
         exit(EXIT_FAILURE);
     }
     break;
   default:
     (void) fprintf(stderr, "Error determining interpolation type\n");
     exit(EXIT_FAILURE);
   }

   /* Check min/max variables */
   fp = in_vol->file;
   fp->using_icv=FALSE;
   if ((fp->maxid != MI_ERROR) && (fp->minid != MI_ERROR) &&
       (fp->datatype!=NC_FLOAT) && (fp->datatype!=NC_DOUBLE)) {
      check_imageminmax(fp, in_vol->volume);
   }

   /* Get space for volume data */
   total_size = 1;
   for (idim=0; idim < WORLD_NDIMS; idim++) {
      index = input_volume_def.axes[idim];
      size = input_volume_def.nelements[idim];
      total_size *= size;
      in_vol->volume->size[index] = size;
   }
   in_vol->volume->data = malloc((size_t) total_size * 
                                 nctypelen(in_vol->volume->datatype));

   /* Get space for slice scale and offset */
   in_vol->volume->scale = 
      malloc(sizeof(double) * in_vol->volume->size[SLC_AXIS]);
   in_vol->volume->offset = 
      malloc(sizeof(double) * in_vol->volume->size[SLC_AXIS]);

   /* Save the program flags */
   *program_flags = args.flags;

   /* Set the default output file datatype */
   if (args.datatype == MI_ORIGINAL_TYPE)
      args.datatype = in_vol->file->datatype;

   /* Explicitly force output files to have regular spacing */
   for (idim=0; idim < WORLD_NDIMS; idim++) {
      if (args.volume_def.coords[idim] != NULL) {
         free(args.volume_def.coords[idim]);
         args.volume_def.coords[idim] = NULL;
      }
   }

   /* Check to see if sign and range have been explicitly set. If not set
      them now */
   if (args.is_signed == INT_MIN) {
      if (args.datatype == in_vol->file->datatype)
         args.is_signed = in_vol->file->is_signed;
      else
         args.is_signed = (args.datatype != NC_BYTE);
   }
   if (args.vrange[0] == -DBL_MAX) {
      if ((args.datatype == in_vol->file->datatype) &&
          (args.is_signed == in_vol->file->is_signed)) {
         args.vrange[0] = in_vol->file->vrange[0];
         args.vrange[1] = in_vol->file->vrange[1];
      }
      else {
         args.vrange[0] = get_default_range(MIvalid_min, args.datatype, 
                                            args.is_signed);
         args.vrange[1] = get_default_range(MIvalid_max, args.datatype, 
                                            args.is_signed);
      }
   }

   /* Set up the file description for the output file */
   out_vol->file = malloc(sizeof(File_Info));
   out_vol->file->ndims = in_vol->file->ndims;
   out_vol->file->datatype = args.datatype;
   out_vol->file->is_signed = args.is_signed;
   out_vol->file->vrange[0] = args.vrange[0];
   out_vol->file->vrange[1] = args.vrange[1];
   for (idim=0; idim < out_vol->file->ndims; idim++) {
      out_vol->file->nelements[idim] = in_vol->file->nelements[idim];
      out_vol->file->world_axes[idim] = in_vol->file->world_axes[idim];
   }
   out_vol->file->keep_real_range = args.keep_real_range;

   /* Get space for output slice */
   out_vol->volume = NULL;
   out_vol->slice = malloc(sizeof(Slice_Data));

   /* Loop through list of axes, getting size of volume and slice */
   total_size = 1;
   for (idim=0; idim < WORLD_NDIMS; idim++) {
      
      /* Get the index for input and output volumes */
      out_vindex = args.volume_def.axes[idim];    /* 0, 1 or 2 */
      out_findex = in_vol->file->indices[out_vindex];   /* 0 to ndims-1 */
      size = args.volume_def.nelements[idim];

      /* Update output axes and indices and nelements */
      out_vol->file->nelements[out_findex] = size;
      out_vol->file->world_axes[out_findex] = idim;
      out_vol->file->axes[idim] = out_vindex;
      out_vol->file->indices[out_vindex] = out_findex;

      /* Update slice size */
      if (out_vindex != 0) {
         out_vol->slice->size[out_vindex-1] = size;
         total_size *= size;
      }
   }
   out_vol->slice->data = malloc((size_t) total_size * sizeof(double));

   /* Create the output file */
   if (args.clobber) {
       cflags = NC_CLOBBER;
   }
   else {
       cflags = NC_NOCLOBBER;
   }
#if MINC2
   if (args.v2format) {
       cflags |= MI2_CREATE_V2;
   }
#endif /* MINC2 */
   create_output_file(outfile, cflags, &args.volume_def, 
                      in_vol->file, out_vol->file,
                      tm_stamp, &args.transform_info);
   
   /* Save the voxel_to_world transformation information */
   out_vol->voxel_to_world = malloc(sizeof(General_transform));
   out_vol->world_to_voxel = malloc(sizeof(General_transform));
   get_voxel_to_world_transf(&args.volume_def, out_vol->voxel_to_world);
   create_inverse_general_transform(out_vol->voxel_to_world,
                                    out_vol->world_to_voxel);

   /* Free the time stamp */
   free(tm_stamp);

}

/* ----------------------------- MNI Header -----------------------------------
@NAME       : check_imageminmax
@INPUT      : fp - pointer to file description
              volume - pointer to volume description
@OUTPUT     : 
@RETURNS    : (nothing)
@DESCRIPTION: Routine to check that MIimagemax and MIimagemin do not vary
              over volume rows and columns. If they do, set up an icv to
              handle it.
@METHOD     : 
@GLOBALS    : 
@CALLS      : 
@CREATED    : August 5, 1993 (Peter Neelin)
@MODIFIED   : 
---------------------------------------------------------------------------- */
static void check_imageminmax(File_Info *fp, Volume_Data *volume)
{
   int ndims, idim, dim[MAX_VAR_DIMS], imgdim[MAX_VAR_DIMS];
   int ivar, varid;

   /* Get MIimage dimensions */
   (void) ncvarinq(fp->mincid, fp->imgid, NULL, NULL, &ndims, imgdim, NULL);

   /* Check MIimagemax/min dimensions */
   for (ivar=0; ivar<2; ivar++) {
      varid = (ivar==0 ? fp->maxid : fp->minid);
      (void) ncvarinq(fp->mincid, varid, NULL, NULL, &ndims, dim, NULL);
      for (idim=0; idim < ndims; idim++) {
         if ((dim[idim] == imgdim[fp->indices[ROW_AXIS]]) ||
             (dim[idim] == imgdim[fp->indices[COL_AXIS]])) {
            fp->using_icv = TRUE;
         }
      }        /* End loop over MIimagemax/min dimensions */
   }        /* End loop over variables MIimagemax/min */

   /* Set up an icv if needed to handle values varying over slice dims. */
   if (fp->using_icv) {

      /* Change type to floating point so that there is no loss of 
         precision (except possibly for long values). */
      if (volume->datatype != NC_DOUBLE)
         volume->datatype = NC_FLOAT;
      volume->is_signed = TRUE;

      /* Create the icv */
      fp->icvid = miicv_create();
      (void) miicv_setint(fp->icvid, MI_ICV_TYPE, volume->datatype);
      (void) miicv_setstr(fp->icvid, MI_ICV_SIGN, 
                          (volume->is_signed ? MI_SIGNED : MI_UNSIGNED));
      (void) miicv_setint(fp->icvid, MI_ICV_DO_NORM, TRUE);
      (void) miicv_setint(fp->icvid, MI_ICV_DO_FILLVALUE, TRUE);
      (void) miicv_attach(fp->icvid, fp->mincid, fp->imgid);

      /* Get max and min for doing valid range checking */
      (void) miicv_inqdbl(fp->icvid, MI_ICV_NORM_MIN, &volume->vrange[0]);
      (void) miicv_inqdbl(fp->icvid, MI_ICV_NORM_MAX, &volume->vrange[1]);
         
   }

}

/* ----------------------------- MNI Header -----------------------------------
@NAME       : get_file_info
@INPUT      : filename - name of file to read
              initialized_volume_def - if TRUE, then volume_def is taken 
                 as being properly initialized and arrays are freed if
                 non-NULL. Otherwise arrays are not freed.
@OUTPUT     : volume_def - description of volume
              file_info - description of file
@RETURNS    : (nothing)
@DESCRIPTION: Routine to get information about the volume definition of
              a minc file. 
@METHOD     : 
@GLOBALS    : 
@CALLS      : 
@CREATED    : February 9, 1993 (Peter Neelin)
@MODIFIED   : 
---------------------------------------------------------------------------- */
static void get_file_info(char *filename, int initialized_volume_def, 
                          Volume_Definition *volume_def,
                          File_Info *file_info)
{
   int dim[MAX_VAR_DIMS], dimid;
   int axis_counter, idim, jdim, cur_axis;
   int varndims, vardim[MAX_VAR_DIMS];
   long varstart, varcount, dimlength;
   char attstr[MI_MAX_ATTSTR_LEN];
   char dimname[MAX_NC_NAME];
   enum {UNKNOWN, REGULAR, IRREGULAR} coord_spacing;

   /* Open the minc file */
   file_info->mincid = miopen(filename, NC_NOWRITE);
   file_info->name = filename;

   /* Get variable identifiers */
   file_info->imgid = ncvarid(file_info->mincid, MIimage);
   ncopts = 0;
   file_info->maxid = ncvarid(file_info->mincid, MIimagemax);
   file_info->minid = ncvarid(file_info->mincid, MIimagemin);
   ncopts = NC_VERBOSE | NC_FATAL;

   /* Get information about datatype dimensions of variable */
   (void) miget_datatype(file_info->mincid, file_info->imgid, 
                         &file_info->datatype, &file_info->is_signed);

   /* Get valid max and min */
   (void) miget_valid_range(file_info->mincid, file_info->imgid, 
                            file_info->vrange);

   /* Get information about dimensions */
   (void) ncvarinq(file_info->mincid, file_info->imgid, NULL, NULL,
                   &file_info->ndims, dim, NULL);

   /* Set variables for keeping track of spatial dimensions */
   axis_counter = 0;                   /* Keeps track of values for axes */

   /* Initialize volume definition variables */
   for (idim=0; idim < WORLD_NDIMS; idim++) {
      volume_def->axes[idim] = NO_AXIS;
      volume_def->step[idim] = 1.0;
      volume_def->start[idim] = 0.0;
      for (jdim=0; jdim < WORLD_NDIMS; jdim++) {
         if (jdim==idim)
            volume_def->dircos[idim][jdim] = 1.0;
         else
            volume_def->dircos[idim][jdim] = 0.0;
      }
      if (initialized_volume_def && (volume_def->coords[idim] != NULL)) {
         free(volume_def->coords[idim]);
      }
      volume_def->coords[idim] = NULL;
      (void) strcpy(volume_def->units[idim], "mm");
      (void) strcpy(volume_def->spacetype[idim], MI_NATIVE);
   }

   /* Loop through dimensions, getting dimension information */

   for (idim=0; idim < file_info->ndims; idim++) {

      /* Get size of dimension */
      (void) ncdiminq(file_info->mincid, dim[idim], dimname, 
                      &file_info->nelements[idim]);

      /* Check variable name */
      cur_axis = NO_AXIS;
      if (strcmp(dimname, MIxspace)==0)
         cur_axis = XAXIS;
      else if (strcmp(dimname, MIyspace)==0)
         cur_axis = YAXIS;
      else if (strcmp(dimname, MIzspace)==0)
         cur_axis = ZAXIS;

      /* Save world axis info */
      file_info->world_axes[idim] = cur_axis;

      /* Check for spatial dimension */
      if (cur_axis == NO_AXIS) continue;

      /* Set axis */
      if (volume_def->axes[cur_axis] != NO_AXIS) {
         (void) fprintf(stderr, "Repeated spatial dimension %s in file %s.\n",
                 dimname, filename);
         exit(EXIT_FAILURE);
      }
      volume_def->axes[cur_axis] = axis_counter++;

      /* Save spatial axis specific info */
      file_info->axes[cur_axis] = volume_def->axes[cur_axis];
      file_info->indices[volume_def->axes[cur_axis]] = idim;
      volume_def->nelements[cur_axis] = file_info->nelements[idim];

      /* Check for existence of variable */
      ncopts = 0;
      dimid = ncvarid(file_info->mincid, dimname);
      ncopts = NC_VERBOSE | NC_FATAL;
      if (dimid == MI_ERROR) continue;
             
      /* Get attributes from variable */
      ncopts = 0;
      (void) miattget1(file_info->mincid, dimid, MIstep, 
                       NC_DOUBLE, &volume_def->step[cur_axis]);
      if (volume_def->step[cur_axis] == 0.0)
         volume_def->step[cur_axis] = 1.0;
      (void) miattget1(file_info->mincid, dimid, MIstart, 
                       NC_DOUBLE, &volume_def->start[cur_axis]);
      (void) miattget(file_info->mincid, dimid, MIdirection_cosines, 
                      NC_DOUBLE, WORLD_NDIMS, 
                      volume_def->dircos[cur_axis], NULL);
      (void) miattgetstr(file_info->mincid, dimid, MIunits, 
                         MI_MAX_ATTSTR_LEN, volume_def->units[cur_axis]);
      (void) miattgetstr(file_info->mincid, dimid, MIspacetype, 
                         MI_MAX_ATTSTR_LEN, volume_def->spacetype[cur_axis]);
      ncopts = NC_VERBOSE | NC_FATAL;

      /* Normalize the direction cosine */
      normalize_vector(volume_def->dircos[cur_axis]);

      /* Look for irregular coordinates for dimension variable */
      ncopts = 0;
      coord_spacing = UNKNOWN;
      dimlength = volume_def->nelements[cur_axis];
      if (miattgetstr(file_info->mincid, dimid, MIspacing, MI_MAX_ATTSTR_LEN,
                       attstr) != NULL) {
         if (strcmp(attstr, MI_IRREGULAR) == 0)
            coord_spacing = IRREGULAR;
         else if (strcmp(attstr, MI_REGULAR) == 0)
            coord_spacing = REGULAR;
      }
      if (ncvarinq(file_info->mincid, dimid, NULL, NULL, 
                   &varndims, vardim, NULL) == MI_ERROR) {
         ncopts = NC_VERBOSE | NC_FATAL;
         continue;
      }
      if ((coord_spacing != REGULAR) && 
          (varndims == 1) && (vardim[0] == dim[idim])) {
         coord_spacing = IRREGULAR;
      }
      if ((coord_spacing == UNKNOWN) || (dimlength <= 1)) {
         coord_spacing = REGULAR;
      }
      if (coord_spacing == IRREGULAR) {
         volume_def->coords[cur_axis] = malloc(sizeof(double) * dimlength);
         varstart = 0;
         varcount = dimlength;
         if (mivarget(file_info->mincid, dimid, &varstart, &varcount,
                      NC_DOUBLE, MI_SIGNED, volume_def->coords[cur_axis])
                   == MI_ERROR) {
            ncopts = NC_VERBOSE | NC_FATAL;
            free(volume_def->coords[cur_axis]);
            volume_def->coords[cur_axis] = NULL;
            continue;
         }
         volume_def->start[cur_axis] = volume_def->coords[cur_axis][0];
         if (dimlength > 1) {
            volume_def->step[cur_axis] = 
               (volume_def->coords[cur_axis][dimlength-1] - 
                            volume_def->coords[cur_axis][0]) /
                               (dimlength - 1);
            if (volume_def->step[cur_axis] == 0.0)
               volume_def->step[cur_axis] = 1.0;
         }
      }
      ncopts = NC_VERBOSE | NC_FATAL;

   }   /* End of loop over dimensions */

   /* Check that we have the correct number of spatial dimensions */
   if (axis_counter != WORLD_NDIMS) {
      (void) fprintf(stderr, 
                     "Incorrect number of spatial dimensions in file %s.\n",
                     filename);
         exit(EXIT_FAILURE);
   }

   return;
}

/* ----------------------------- MNI Header -----------------------------------
@NAME       : get_args_volume_def
@INPUT      : input_volume_def - description of input volume
@OUTPUT     : args_volume_def - description of new output volume
@RETURNS    : (nothing)
@DESCRIPTION: Routine to copy appropriate information from input volume 
              definition to output volume definition, overriding values 
              not set on the command line.
@METHOD     : 
@GLOBALS    : 
@CALLS      : 
@CREATED    : November 7, 1995 (Peter Neelin)
@MODIFIED   : 
---------------------------------------------------------------------------- */
static void get_args_volume_def(Volume_Definition *input_volume_def,
                                Volume_Definition *args_volume_def)
{
   int idim, jdim;

   /* Loop over dimensions */
   for (idim=0; idim < WORLD_NDIMS; idim++) {

      /* Copy values, as needed */
      if (args_volume_def->axes[idim] == NO_AXIS)
         args_volume_def->axes[idim] = input_volume_def->axes[idim];
      if (args_volume_def->nelements[idim] == 0)
         args_volume_def->nelements[idim] = input_volume_def->nelements[idim];
      if (args_volume_def->step[idim] == 0.0)
         args_volume_def->step[idim] = input_volume_def->step[idim];
      if (args_volume_def->start[idim] == NO_VALUE)
         args_volume_def->start[idim] = input_volume_def->start[idim];
      if (args_volume_def->dircos[idim][0] == NO_VALUE) {
         for (jdim=0; jdim < WORLD_NDIMS; jdim++)
            args_volume_def->dircos[idim][jdim] = 
               input_volume_def->dircos[idim][jdim];
      }
      if (strlen(args_volume_def->units[idim]) == 0)
         (void) strcpy(args_volume_def->units[idim],
                        input_volume_def->units[idim]);
      if (strlen(args_volume_def->spacetype[idim]) == 0)
         (void) strcpy(args_volume_def->spacetype[idim],
                        input_volume_def->spacetype[idim]);

   }
}

/* ----------------------------- MNI Header -----------------------------------
@NAME       : transform_volume_def
@INPUT      : transform_info - description of output to input transform
                 (if NULL, then don't transform).
              input_volume_def - description of input volume
@OUTPUT     : transformed_volume_def - description of new output volume
@RETURNS    : (nothing)
@DESCRIPTION: Routine to copy appropriate information from input volume 
              definition to a new volume definition, after transformation 
              from with the output to input transform.
@METHOD     : 
@GLOBALS    : 
@CALLS      : 
@CREATED    : November 7, 1995 (Peter Neelin)
@MODIFIED   : 
---------------------------------------------------------------------------- */
static void transform_volume_def(Transform_Info *transform_info,
                                 Volume_Definition *input_volume_def,
                                 Volume_Definition *transformed_volume_def)
{
   int idim, jdim;
   Coord_Vector origin = {0, 0, 0};
   Coord_Vector transformed_origin, vector;
   double length, nelements;

   /* Copy the volume definition */
   *transformed_volume_def = *input_volume_def;

   /* Modify things if the transform has changed */
   if ((transform_info != NULL) && (transform_info->file_name != NULL)) {

      /* Set up origin vector */
      for (idim=0; idim < WORLD_NDIMS; idim++)
         for (jdim=0; jdim < WORLD_NDIMS; jdim++)
            origin[idim] += input_volume_def->start[jdim] *
               input_volume_def->dircos[jdim][idim];

      /* Transform origin vector */
      DO_INVERSE_TRANSFORM(transformed_origin, 
                           transform_info->transformation, origin);

      /* Loop over dimensions */
      for (idim=0; idim < WORLD_NDIMS; idim++) {

         /* Get number of elements */
         nelements = input_volume_def->nelements[idim] - 1;
         if (nelements < 1) nelements = 1.0;

         /* Transform whole axis */
         VECTOR_SCALAR_MULT(vector, input_volume_def->dircos[idim],
                            nelements * input_volume_def->step[idim]);
         VECTOR_ADD(vector, vector, origin);
         DO_INVERSE_TRANSFORM(vector, transform_info->transformation, vector);
         VECTOR_DIFF(vector, vector, transformed_origin);

         /* Calculate length of transformed axis - new step value */
         length = sqrt(vector[XCOORD]*vector[XCOORD] +
                       vector[YCOORD]*vector[YCOORD] +
                       vector[ZCOORD]*vector[ZCOORD]);
         transformed_volume_def->step[idim] = length / nelements;

         /* Calculate new direction cosines */
         for (jdim=0; jdim < WORLD_NDIMS; jdim++) {
            transformed_volume_def->dircos[idim][jdim] = 
               (length > 0.0 ? vector[jdim] / length : 
                (jdim == idim ? 1.0 : 0.0));
         }

         /* Make sure that direction cosines point the right way */
         if (transformed_volume_def->dircos[idim][idim] < 0.0) {
            VECTOR_SCALAR_MULT(transformed_volume_def->dircos[idim],
                               transformed_volume_def->dircos[idim],
                               -1.0);
            transformed_volume_def->step[idim] *= -1.0;
         }
      }

      /* Calculate the start along each axis - check for bad dircos. */
      if (convert_origin_to_start(transformed_origin,
                                  transformed_volume_def->dircos[XCOORD],
                                  transformed_volume_def->dircos[YCOORD],
                                  transformed_volume_def->dircos[ZCOORD],
                                  transformed_volume_def->start) != 0) {

         /* If dircos are bad, set them to default */
         for (idim=0; idim < WORLD_NDIMS; idim++) {
            for (jdim=0; jdim < WORLD_NDIMS; jdim++) {
               transformed_volume_def->dircos[idim][jdim] = 
                  ((idim==jdim) ? 1.0 : 0.0);
            }
         }

         /* Try to convert point again */
         if (convert_origin_to_start(transformed_origin,
                                     transformed_volume_def->dircos[XCOORD],
                                     transformed_volume_def->dircos[YCOORD],
                                     transformed_volume_def->dircos[ZCOORD],
                                     transformed_volume_def->start) != 0) {
            (void) fprintf(stderr, 
                           "Serious problem converting origin to start!\n");
            exit(EXIT_FAILURE);
         }
      }

   }

}

/* ----------------------------- MNI Header -----------------------------------
@NAME       : is_zero_vector
@INPUT      : vector - 3D vector
@OUTPUT     : (none)
@RETURNS    : 1 if vector has zero length, 0 otherwise
@DESCRIPTION: Routine to check for a zero length vector.
@METHOD     : 
@GLOBALS    : 
@CALLS      : 
@CREATED    : November 9, 1995 (Peter Neelin)
@MODIFIED   : 
---------------------------------------------------------------------------- */
static int is_zero_vector(double vector[])
{
   return ((vector[0] == 0.0) &&
           (vector[1] == 0.0) &&
           (vector[2] == 0.0));
}

/* ----------------------------- MNI Header -----------------------------------
@NAME       : normalize_vector
@INPUT      : vector - 3D vector
@OUTPUT     : (none)
@RETURNS    : (nothing)
@DESCRIPTION: Routine to normalize a vector
@METHOD     : 
@GLOBALS    : 
@CALLS      : 
@CREATED    : November 9, 1995 (Peter Neelin)
@MODIFIED   : 
---------------------------------------------------------------------------- */
static void normalize_vector(double vector[])
{
   int idim;
   double magnitude;

   /* Normalize the direction cosine */
   magnitude = 0.0;
   for (idim=0; idim < WORLD_NDIMS; idim++) {
      magnitude += (vector[idim] * vector[idim]);
   }
   magnitude = sqrt(magnitude);
   if (magnitude > 0.0) {
      for (idim=0; idim < WORLD_NDIMS; idim++) {
         vector[idim] /= magnitude;
      }
   }

}

/* ----------------------------- MNI Header -----------------------------------
@NAME       : create_output_file
@INPUT      : filename - name of file to create
              cflags - flag creation flags
              volume_def - description of volume
              in_file - description of input file
              out_file - description of output file
              tm_stamp - string describing program invocation
              transformation - transformation to be applied to data
@OUTPUT     : (nothing) 
@RETURNS    : (nothing)
@DESCRIPTION: Routine to create an minc output file and set up its parameters
              properly.
@METHOD     : 
@GLOBALS    : 
@CALLS      : 
@CREATED    : February 9, 1993 (Peter Neelin)
@MODIFIED   : 
---------------------------------------------------------------------------- */
static void create_output_file(char *filename, int cflags, 
                               Volume_Definition *volume_def,
                               File_Info *in_file,
                               File_Info *out_file,
                               char *tm_stamp, 
                               Transform_Info *transform_info)
{
   int ndims, in_dims[MAX_VAR_DIMS], out_dims[MAX_VAR_DIMS];
   int out_maxmin_dims[MAX_VAR_DIMS];
   char dimname[MAX_NC_NAME];
   int nmaxmin_dims, nimage_dims, axis, idim, dimid, varid, itrans;
   int att_length;
   nc_type datatype;
   int nexcluded, excluded_vars[10];
   char *string;
   int dim_exists, is_volume_dimension;
   int in_index, out_index;

   /* Save the file name */
   out_file->name = filename;

   /* Create the list of excluded variables */
   nexcluded = 0;
   ncopts = 0;
   if ((varid=ncvarid(in_file->mincid, MIxspace)) != MI_ERROR)
      excluded_vars[nexcluded++] = varid;
   if ((varid=ncvarid(in_file->mincid, MIyspace)) != MI_ERROR)
      excluded_vars[nexcluded++] = varid;
   if ((varid=ncvarid(in_file->mincid, MIzspace)) != MI_ERROR)
      excluded_vars[nexcluded++] = varid;
   if ((varid=ncvarid(in_file->mincid, MIimage)) != MI_ERROR)
      excluded_vars[nexcluded++] = varid;
   if ((varid=ncvarid(in_file->mincid, MIimagemax)) != MI_ERROR)
      excluded_vars[nexcluded++] = varid;
   if ((varid=ncvarid(in_file->mincid, MIimagemin)) != MI_ERROR)
      excluded_vars[nexcluded++] = varid;
#if MINC2
   if ((varid=ncvarid(in_file->mincid, MIvector_dimension)) != MI_ERROR)
      excluded_vars[nexcluded++] = varid;
#endif /* MINC2 */
   ncopts = NC_VERBOSE | NC_FATAL;

   /* Create the file */
   out_file->mincid = micreate(filename, cflags);
 
   /* Copy all other variable definitions */
   (void) micopy_all_var_defs(in_file->mincid, out_file->mincid, 
                              nexcluded, excluded_vars);

   /* Add the time stamp */
   ncopts=0;
   if ((ncattinq(out_file->mincid, NC_GLOBAL, MIhistory, &datatype,
                 &att_length) == MI_ERROR) ||
       (datatype != NC_CHAR))
      att_length = 0;
   att_length += strlen(tm_stamp) + 1;
   string = malloc(att_length);
   string[0] = '\0';
   (void) miattgetstr(out_file->mincid, NC_GLOBAL, MIhistory, att_length, 
                      string);
   ncopts=NC_VERBOSE | NC_FATAL;
   (void) strcat(string, tm_stamp);
   (void) miattputstr(out_file->mincid, NC_GLOBAL, MIhistory, string);
   free(string);

   /* Get the dimension ids from the input file */
   (void) ncvarinq(in_file->mincid, in_file->imgid, NULL, NULL, 
                   &ndims, in_dims, NULL);

   /* Check for screw-up on number of dimensions */
   if (ndims != out_file->ndims) {
      (void) fprintf(stderr, 
                     "Error in number of dimensions for output file.\n");
      exit(EXIT_FAILURE);
   }

   /* Create the dimensions for the image variable */
   for (out_index=0; out_index<ndims; out_index++) {

      /* Check to see if this is a volume dimension */
      is_volume_dimension = (out_file->world_axes[out_index] != NO_AXIS);

      /* Get the input index */
      if (!is_volume_dimension) 
         in_index = out_index;
      else {
         axis = out_file->world_axes[out_index];
         if ((axis<0) || (axis>=WORLD_NDIMS)) {
            (void) fprintf(stderr,
                           "Error creating dimensions for output file.\n");
            exit(EXIT_FAILURE);
         }
         in_index = in_file->indices[in_file->axes[axis]];
      }

      /* Get the dimension name from the input file */
      (void) ncdiminq(in_file->mincid, in_dims[in_index], dimname, NULL);

      /* Check to see if the dimension already exists */
      ncopts = 0;
      out_dims[out_index] = ncdimid(out_file->mincid, dimname);
      ncopts = NC_VERBOSE | NC_FATAL;
      dim_exists = (out_dims[out_index] != MI_ERROR);

      /* If we have a volume dimension and it exists already with the wrong
         size, then we must rename it */
      if (is_volume_dimension && dim_exists && 
          (out_file->nelements[out_index] != in_file->nelements[in_index])) {
         string = malloc(MAX_NC_NAME);
         ncopts = 0;
         idim = 0;
         do {
            (void) sprintf(string, "%s%d", dimname, idim);
            idim++;
         } while (ncdimid(out_file->mincid, string) != MI_ERROR);
         ncopts = NC_VERBOSE | NC_FATAL;
         (void) ncdimrename(out_file->mincid, out_dims[out_index], string);
         free(string);
         out_dims[out_index] = ncdimdef(out_file->mincid, dimname, 
                                        out_file->nelements[out_index]);
      }
      else if (!dim_exists)
         out_dims[out_index] = ncdimdef(out_file->mincid, dimname, 
                                        out_file->nelements[out_index]);

      /* If this is a volume dimension, create a variable */
      if (is_volume_dimension) {

         /* Create the variable */
         dimid = micreate_group_variable(out_file->mincid, dimname);
         (void) miattputdbl(out_file->mincid, dimid, MIstep, 
                            volume_def->step[axis]);
         (void) miattputdbl(out_file->mincid, dimid, MIstart, 
                            volume_def->start[axis]);
         (void) ncattput(out_file->mincid, dimid, MIdirection_cosines, 
                         NC_DOUBLE, WORLD_NDIMS, volume_def->dircos[axis]);
         (void) miattputstr(out_file->mincid, dimid, MIunits, 
                            volume_def->units[axis]);
         (void) miattputstr(out_file->mincid, dimid, MIspacetype, 
                            volume_def->spacetype[axis]);

      }       /* If volume dimension */
   }       /* Loop over dimensions */

   /* Create the image max and min variables. These do not vary over
      the volume rows and columns even if they are not image dimensions,
      so we have to copy down the elements of the array (excluding image
      dimensions). Compute number of slices (row,columns) in an image
      (minc file def'n) and number of images in the file. */
   nimage_dims = 2;
   ncdiminq(out_file->mincid, out_dims[ndims-1], dimname, NULL);
   if (strcmp(dimname, MIvector_dimension)==0)
      nimage_dims++;
   nmaxmin_dims = 0;
   out_file->slices_per_image = 1;
   out_file->images_per_file = 1;
   for (idim=0; idim<ndims; idim++) {
      if ((idim != out_file->indices[COL_AXIS]) &&
          (idim != out_file->indices[ROW_AXIS])) {
         if (idim < ndims-nimage_dims) {
            out_maxmin_dims[nmaxmin_dims] = out_dims[idim];
            nmaxmin_dims++;
            out_file->images_per_file *= out_file->nelements[idim];
         }
         else {
            out_file->slices_per_image *= out_file->nelements[idim];
         }
      }
   }
   out_file->do_slice_renormalization =
      ((out_file->datatype != NC_FLOAT) && 
       (out_file->datatype != NC_DOUBLE) &&
       (out_file->slices_per_image > 1));

   /* Create the variables */
   out_file->maxid = micreate_std_variable(out_file->mincid, MIimagemax,
                                           NC_DOUBLE, nmaxmin_dims, 
                                           out_maxmin_dims);
   if (in_file->maxid != MI_ERROR)
      (void) micopy_all_atts(in_file->mincid, in_file->maxid,
                             out_file->mincid, out_file->maxid);
   out_file->minid = micreate_std_variable(out_file->mincid, MIimagemin,
                                           NC_DOUBLE, nmaxmin_dims, 
                                           out_maxmin_dims);
   if (in_file->minid != MI_ERROR)
      (void) micopy_all_atts(in_file->mincid, in_file->minid,
                             out_file->mincid, out_file->minid);

   /* Add transformation information to image processing variable if 
      a transformation is given on the command line */

   if (transform_info->file_name != NULL) {

      ncopts = 0;

      /* Get id of processing variable (create it if needed) */
      varid = ncvarid(out_file->mincid, PROCESSING_VAR);
      if (varid == MI_ERROR) {
         varid = ncvardef(out_file->mincid, PROCESSING_VAR, NC_INT, 0, NULL);
         (void) miadd_child(out_file->mincid, 
                            ncvarid(out_file->mincid, MIrootvariable), varid);
      }

      /* Look for an unused transformation attribute */
      string = malloc(MI_MAX_ATTSTR_LEN);
      itrans = 0;
      do {
         (void) sprintf(string, "transformation%d-filename", itrans);
         itrans++;
      } while (ncattinq(out_file->mincid, varid, string,
                        NULL, NULL) != MI_ERROR);
      itrans--;

      /* Reset error handling */
      ncopts = NC_VERBOSE | NC_FATAL;

      /* Add the attributes describing the transformation */
      (void) miattputstr(out_file->mincid, varid, string, 
                         transform_info->file_name);
      (void) sprintf(string, "transformation%d-filedata", itrans);
      (void) miattputstr(out_file->mincid, varid, string,
                         transform_info->file_contents);
      if (transform_info->invert_transform) {
         (void) sprintf(string, "transformation%d-inverted", itrans);
         (void) miattputstr(out_file->mincid, varid, string, MI_TRUE);
      }

      free(string);

   }         /* If transform specified on command line */

   /* Create the image variable */
   out_file->imgid = micreate_std_variable(out_file->mincid, MIimage, 
                                           out_file->datatype,
                                           ndims, out_dims);
   (void) micopy_all_atts(in_file->mincid, in_file->imgid,
                          out_file->mincid, out_file->imgid);
   (void) miattputstr(out_file->mincid, out_file->imgid, MIcomplete,
                      MI_FALSE);
   (void) miset_valid_range(out_file->mincid, out_file->imgid, 
                            out_file->vrange);
   if (out_file->is_signed)
      (void) miattputstr(out_file->mincid, out_file->imgid,
                         MIsigntype, MI_SIGNED);
   else
      (void) miattputstr(out_file->mincid, out_file->imgid,
                         MIsigntype, MI_UNSIGNED);

   /* Get into data mode */
   (void) ncsetfill(out_file->mincid, NC_NOFILL);
   (void) ncendef(out_file->mincid);

   /* Copy all the other data */
   (void) micopy_all_var_values(in_file->mincid, out_file->mincid,
                                nexcluded, excluded_vars);

   /* Create and attach an icv */
   out_file->using_icv = TRUE;
   out_file->icvid = miicv_create();
   (void) miicv_setint(out_file->icvid, MI_ICV_TYPE, NC_DOUBLE);
   (void) miicv_setint(out_file->icvid, MI_ICV_DO_NORM, TRUE);
   (void) miicv_setint(out_file->icvid, MI_ICV_USER_NORM, TRUE);
   (void) miicv_attach(out_file->icvid, out_file->mincid, out_file->imgid);

   return;
}

/* ----------------------------- MNI Header -----------------------------------
@NAME       : get_voxel_to_world_transf
@INPUT      : volume_def - description of volume
@OUTPUT     : voxel_to_world - transformation
@RETURNS    : (nothing)
@DESCRIPTION: Routine to convert a Volume definition specification of sampling
              to a voxel-to-world transformation
@METHOD     : 
@GLOBALS    : 
@CALLS      : 
@CREATED    : February 9, 1993 (Peter Neelin)
@MODIFIED   : 
---------------------------------------------------------------------------- */
static void get_voxel_to_world_transf(Volume_Definition *volume_def, 
                                      General_transform *voxel_to_world)
{
   int idim, jdim, cur_dim, vol_axis;
   int irregular;
   long ielement, dimlength;
   Transform matrix;
   Irregular_Transform_Data *irreg_transf_data;
   General_transform linear_transf, irreg_transf;

   /* Make an identity matrix */
   make_identity_transform(&matrix);

   /* Loop over rows of matrix */
   for (idim=0; idim<WORLD_NDIMS; idim++) {

      /* Loop over columns of matrix */
      for (jdim=0; jdim<WORLD_NDIMS; jdim++) {
         cur_dim = volume_def->axes[jdim];

         /* Get rotation/scale components of matrix */
         Transform_elem(matrix, idim, cur_dim) = 
            volume_def->step[jdim] * volume_def->dircos[jdim][idim];

         /* Get translation components */
         Transform_elem(matrix, idim, VOL_NDIMS) +=
            volume_def->start[jdim] * volume_def->dircos[jdim][idim];
      }
   }

   /* Save the general transform */
   create_linear_transform(voxel_to_world, &matrix);

   /* Check for an irregularly spaced dimension */
   irregular = FALSE;
   for (idim=0; idim < WORLD_NDIMS; idim++) {
      if (volume_def->coords[idim] != NULL)
         irregular = TRUE;
   }

   /* If we have an irregularly spaced dimension, then create the appropriate
      transform */
   if (irregular) {
      irreg_transf_data = malloc(sizeof(*irreg_transf_data));

      /* Loop through the axes */
      for (idim=0; idim < WORLD_NDIMS; idim++) {

         vol_axis = volume_def->axes[idim];
         irreg_transf_data->last_index[vol_axis] = 0;

         /* Check whether the axis is irregularly spaced or not */
         if (volume_def->coords[idim] == NULL) {
            irreg_transf_data->nelements[vol_axis] = 0;
            irreg_transf_data->coords[vol_axis] = NULL;
         }
         else {

            /* If irregular then get the number of elements and
               allocate space */
            dimlength = volume_def->nelements[idim];
            irreg_transf_data->nelements[vol_axis] = dimlength;
            irreg_transf_data->coords[vol_axis] =
               malloc(sizeof(double) * dimlength);

            /* Normalize the coordinates to give first coord 0 and
               last coord n-1 (so that we can concat with the linear
               transform already created */
            for (ielement=0; ielement < dimlength; ielement++) {
               irreg_transf_data->coords[vol_axis][ielement] = 
                  (volume_def->coords[idim][ielement] - 
                              volume_def->start[idim]) /
                                 volume_def->step[idim];
            }
         }
      }

      /* Create an irregular transform and free the data (we only free
         the Irregular_Transform_Data structure, not the coord pointers,
         since these arrays are not copied) */
      create_user_transform(&irreg_transf, (void *) irreg_transf_data,
                            sizeof(*irreg_transf_data),
                            irregular_transform_function,
                            irregular_inverse_transform_function);
      free(irreg_transf_data);

      /* Concatenate the linear transform with the irregular transform */
      copy_general_transform(voxel_to_world, &linear_transf);
      delete_general_transform(voxel_to_world);
      concat_general_transforms(&irreg_transf, &linear_transf,
                                voxel_to_world);

   }

   return;
}

/* ----------------------------- MNI Header -----------------------------------
@NAME       : irregular_transform_function
@INPUT      : user_data - pointer to user data
              x, y, z   - coordinate to transform
@OUTPUT     : x_trans, y_trans, z_trans - resulting coordinate
@RETURNS    : (nothin)
@DESCRIPTION: Routine to transform irregularly spaced coordinate to a 
              regular spacing.
@METHOD     : 
@GLOBALS    : 
@CALLS      : 
@CREATED    : November 4, 1993 (Peter Neelin)
@MODIFIED   : 
---------------------------------------------------------------------------- */
static void irregular_transform_function(void *user_data,
                                         Real x,
                                         Real y,
                                         Real z,
                                         Real *x_trans,
                                         Real *y_trans,
                                         Real *z_trans)
{
   Irregular_Transform_Data *irreg_transf_data;
   int idim;
   long dimlength, index;
   double coord[VOL_NDIMS], coord_transf[VOL_NDIMS];
   double frac;

   /* Get the pointer to the data */
   irreg_transf_data = (Irregular_Transform_Data *) user_data;

   /* Get the coordinate */
   coord[0] = x; coord[1] = y; coord[2] = z;

   /* Loop through axes, renormalizing coordinate */
   for (idim=0; idim < VOL_NDIMS; idim++) {
      dimlength = irreg_transf_data->nelements[idim];
      if (dimlength <= 1) {
         coord_transf[idim] = coord[idim];
      }
      else {

         /* For irregular dimension, do linear interpolation between coords */
         index = (long) coord[idim];
         if (index < 0) index = 0;
         if (index > dimlength-2) index = dimlength-2;
         frac = coord[idim] - index;
         coord_transf[idim] = 
            (1.0 - frac) * irreg_transf_data->coords[idim][index] + 
            frac * irreg_transf_data->coords[idim][index + 1];

      }
   }

   /* Save the coordinates */
   *x_trans = coord_transf[0];
   *y_trans = coord_transf[1];
   *z_trans = coord_transf[2];

   return;
}

/* ----------------------------- MNI Header -----------------------------------
@NAME       : irregular_inverse_transform_function
@INPUT      : user_data - pointer to user data
              x, y, z   - coordinate to transform
@OUTPUT     : x_trans, y_trans, z_trans - resulting coordinate
@RETURNS    : (nothin)
@DESCRIPTION: Routine to transform irregularly spaced coordinate to a 
              regular spacing.
@METHOD     : 
@GLOBALS    : 
@CALLS      : 
@CREATED    : November 4, 1993 (Peter Neelin)
@MODIFIED   : 
---------------------------------------------------------------------------- */
static void irregular_inverse_transform_function(void *user_data,
                                                 Real x,
                                                 Real y,
                                                 Real z,
                                                 Real *x_trans,
                                                 Real *y_trans,
                                                 Real *z_trans)
{
   Irregular_Transform_Data *irreg_transf_data;
   int idim, not_found;
   long dimlength, index;
   double coord[VOL_NDIMS], coord_transf[VOL_NDIMS];
   double first, next, step, frac;

   /* Get the pointer to the data */
   irreg_transf_data = (Irregular_Transform_Data *) user_data;

   /* Get the coordinate */
   coord[0] = x; coord[1] = y; coord[2] = z;

   /* Loop through axes, renormalizing coordinate */
   for (idim=0; idim < VOL_NDIMS; idim++) {
      dimlength = irreg_transf_data->nelements[idim];
      if (dimlength <= 1) {
         coord_transf[idim] = coord[idim];
      }
      else {

         /* Search for the closest index (checking for out-of-range values) */
         index = irreg_transf_data->last_index[idim];
         if (index < 0) index = 0;
         if (index > dimlength-2) index = dimlength-2;
         not_found = TRUE;
         while (not_found) {
            if (coord[idim] < irreg_transf_data->coords[idim][index]) {
               if (index > 0) index--;
               else {
                  index = 0;
                  not_found = FALSE;
               }
            }
            else if (coord[idim] > irreg_transf_data->coords[idim][index+1]) {
               if (index < dimlength-2) index++;
               else {
                  index = dimlength-2;
                  not_found = FALSE;
               }
            }
            else {
               not_found = FALSE;
            }
            
         }
         irreg_transf_data->last_index[idim] = index;

         /* Do linear interpolation */
         first = irreg_transf_data->coords[idim][index];
         next = irreg_transf_data->coords[idim][index + 1];
         step = next - first;
         if (step == 0.0)
            frac = 0.0;
         else
            frac = (coord[idim] - first) / step;
         coord_transf[idim] = 
            (1.0 - frac) * index + frac * (index + 1);

      }
   }

   /* Save the coordinates */
   *x_trans = coord_transf[0];
   *y_trans = coord_transf[1];
   *z_trans = coord_transf[2];

   return;
}

/* ----------------------------- MNI Header -----------------------------------
@NAME       : get_default_range
@INPUT      : what     - MIvalid_min means get default min, MIvalid_min means 
                 get default min
              datatype - type of variable
              is_signed   - TRUE if variable is signed
@OUTPUT     : (none)
@RETURNS    : default maximum or minimum for the datatype
@DESCRIPTION: Return the defaults maximum or minimum for a given datatype
              and sign.
@METHOD     : 
@GLOBALS    : 
@CALLS      : 
@CREATED    : August 10, 1992 (Peter Neelin)
@MODIFIED   : February 10, 1993 (Peter Neelin)
                 - ripped off from MINC code
---------------------------------------------------------------------------- */
static double get_default_range(char *what, nc_type datatype, int is_signed)
{
   double limit;

   if (strcmp(what, MIvalid_max)==0) {
      switch (datatype) {
      case NC_INT:   limit = (is_signed) ? INT_MAX : UINT_MAX; break;
      case NC_SHORT: limit = (is_signed) ? SHRT_MAX : USHRT_MAX; break;
      case NC_BYTE:  limit = (is_signed) ? SCHAR_MAX : UCHAR_MAX; break;
      default: limit = DEFAULT_MAX; break;
      }
   }
   else if (strcmp(what, MIvalid_min)==0) {
      switch (datatype) {
      case NC_INT:  limit = (is_signed) ? INT_MIN : 0; break;
      case NC_SHORT: limit = (is_signed) ? SHRT_MIN : 0; break;
      case NC_BYTE:  limit = (is_signed) ? SCHAR_MIN : 0; break;
      default: limit = DEFAULT_MIN; break;
      }
   }
   else {
      limit = 0.0;
   }

   return limit;
}


/* ----------------------------- MNI Header -----------------------------------
@NAME       : finish_up
@INPUT      : in_vol - input volume
              out_vol - output volume
@OUTPUT     : (nothing) 
@RETURNS    : (nothing)
@DESCRIPTION: Routine to finish up at end of program, closing files, etc.
@METHOD     : 
@GLOBALS    : 
@CALLS      : 
@CREATED    : February 15, 1993 (Peter Neelin)
@MODIFIED   : 
---------------------------------------------------------------------------- */
static void finish_up(VVolume *in_vol, VVolume *out_vol)
{
   File_Info *in_file, *out_file;

   /* Get file info pointers */
   in_file = in_vol->file;
   out_file = out_vol->file;

   /* Close the output file */
   (void) miattputstr(out_file->mincid, out_file->imgid, MIcomplete, MI_TRUE);
   if (out_file->using_icv) {
      (void) miicv_free(out_file->icvid);
   }
   (void) miclose(out_file->mincid);

   /* Close the input file */
   if (in_file->using_icv) {
      (void) miicv_free(in_file->icvid);
   }
   (void) miclose(in_file->mincid);

   return;
}

/* ----------------------------- MNI Header -----------------------------------
@NAME       : get_transformation
@INPUT      : dst - Pointer to client data from argument table
              key - argument key
              nextArg - argument following key
@OUTPUT     : (nothing) 
@RETURNS    : TRUE so that ParseArgv will discard nextArg, unless there
              is no following argument.
@DESCRIPTION: Routine called by ParseArgv to read in a transformation file
@METHOD     : 
@GLOBALS    : 
@CALLS      : 
@CREATED    : February 15, 1993 (Peter Neelin)
@MODIFIED   : 
---------------------------------------------------------------------------- */
static int get_transformation(char *dst, char *key, char *nextArg)
     /* ARGSUSED */
{
   Transform_Info *transform_info;
   General_transform *transformation;
   FILE *fp;
   int ch, index;

   /* Check for following argument */
   if (nextArg == NULL) {
      (void) fprintf(stderr, 
                     "\"%s\" option requires an additional argument\n",
                     key);
      exit(EXIT_FAILURE);
   }

   /* Get pointer to transform info structure */
   transform_info = (Transform_Info *) dst;

   /* Save file name */
   transform_info->file_name = nextArg;
   transformation = transform_info->transformation;

   /* Open the file */
   if (strcmp(nextArg, "-") == 0) {
      /* Create a temporary for standard input */
      fp=tmpfile();
      if (fp==NULL) {
         (void) fprintf(stderr, "Error opening temporary file.\n");
         exit(EXIT_FAILURE);
      }
      while ((ch=getc(stdin))!=EOF) (void) putc(ch, fp);
      rewind(fp);
   }
   else {
      if (open_file_with_default_suffix(nextArg,
                    get_default_transform_file_suffix(),
                    READ_FILE, ASCII_FORMAT, &fp) != OK) {
         (void) fprintf(stderr, "Error opening transformation file %s.\n",
                        nextArg);
         exit(EXIT_FAILURE);
      }
   }

   /* Read in the file for later use */
   if (transform_info->file_contents == NULL) {
      transform_info->file_contents = malloc(TRANSFORM_BUFFER_INCREMENT);
      transform_info->buffer_length = TRANSFORM_BUFFER_INCREMENT;
   }
   for (index = 0; (ch=getc(fp)) != EOF; index++) {
      if (index >= transform_info->buffer_length-1) {
         transform_info->buffer_length += TRANSFORM_BUFFER_INCREMENT;
         transform_info->file_contents = 
            realloc(transform_info->file_contents, 
                    transform_info->buffer_length);
      }
      transform_info->file_contents[index] = (char) ch;
   }
   transform_info->file_contents[index] = '\0';
   rewind(fp);

   /* Get rid of the old transformation */
   delete_general_transform(transformation);

   /* Read the file */
   if (input_transform(fp, transform_info->file_name, 
                       transformation)!=OK) {
      (void) fprintf(stderr, "Error reading transformation file.\n");
      exit(EXIT_FAILURE);
   }
   (void) close_file(fp);

#ifdef TRANSFORM_CHANGE_KLUDGE
   Specified_transform = TRUE;
#endif

   return TRUE;
}

/* ----------------------------- MNI Header -----------------------------------
@NAME       : get_model_file
@INPUT      : dst - Pointer to client data from argument table
              key - argument key
              nextArg - argument following key
@OUTPUT     : (nothing) 
@RETURNS    : TRUE so that ParseArgv will discard nextArg unless there
              is no following argument.
@DESCRIPTION: Routine called by ParseArgv to read in a model file (-like)
@METHOD     : 
@GLOBALS    : 
@CALLS      : 
@CREATED    : February 15, 1993 (Peter Neelin)
@MODIFIED   : 
---------------------------------------------------------------------------- */
static int get_model_file(char *dst, char *key, char *nextArg)
     /* ARGSUSED */
{
   Volume_Definition *volume_def;
   File_Info file;

   /* Check for following argument */
   if (nextArg == NULL) {
      (void) fprintf(stderr, 
                     "\"%s\" option requires an additional argument\n",
                     key);
      exit(EXIT_FAILURE);
   }

   /* Get pointer to volume definition structure */
   volume_def = (Volume_Definition *) dst;

   /* Get file information */
   get_file_info(nextArg, TRUE, volume_def, &file);

   /* Close the file */
   (void) miclose(file.mincid);

#ifdef TRANSFORM_CHANGE_KLUDGE
   Specified_like = TRUE;
#endif

   return TRUE;
}

/* ----------------------------- MNI Header -----------------------------------
@NAME       : set_standard_sampling
@INPUT      : dst - Pointer to client data from argument table
              key - argument key
              nextArg - argument following key
@OUTPUT     : (nothing) 
@RETURNS    : FALSE so that ParseArgv will not discard nextArg.
@DESCRIPTION: Routine called by ParseArgv to set the sampling to standard
              values (sets only step, start and direction cosines).
@METHOD     : 
@GLOBALS    : 
@CALLS      : 
@CREATED    : November 14, 1995 (Peter Neelin)
@MODIFIED   : 
---------------------------------------------------------------------------- */
static int set_standard_sampling(char *dst, char *key, char *nextArg)
     /* ARGSUSED */
{
   Volume_Definition *volume_def;
   int idim, jdim;

   /* Get pointer to volume definition structure */
   volume_def = (Volume_Definition *) dst;

   /* Set sensible values */
   for (idim=0; idim < WORLD_NDIMS; idim++) {
      volume_def->step[idim] = 1.0;
      volume_def->start[idim] = 0.0;
      for (jdim=0; jdim < WORLD_NDIMS; jdim++) {
         volume_def->dircos[idim][jdim] = (idim == jdim ? 1.0 : 0.0);
      }
   }

   return FALSE;
}

/* ----------------------------- MNI Header -----------------------------------
@NAME       : set_spacetype
@INPUT      : dst - Pointer to client data from argument table
              key - argument key
              nextArg - argument following key
@OUTPUT     : (nothing) 
@RETURNS    : TRUE if nextArg should be discarded, FALSE otherwise
@DESCRIPTION: Routine called by ParseArgv to set the space type of the
              output sampling.
@METHOD     : 
@GLOBALS    : 
@CALLS      : 
@CREATED    : December 12, 1995 (Peter Neelin)
@MODIFIED   : 
---------------------------------------------------------------------------- */
static int set_spacetype(char *dst, char *key, char *nextArg)
     /* ARGSUSED */
{
   Volume_Definition *volume_def;
   char *spacetype;
   int idim;
   int return_value;

   /* Get pointer to client data */
   volume_def = (Volume_Definition *) dst;

   /* Check key for spacetype */
   return_value = FALSE;
   if (strcmp(key, "-talairach") == 0) {
      spacetype = MI_TALAIRACH;
   }
   else {

      /* Check for following argument */
      if (nextArg == NULL) {
         (void) fprintf(stderr, 
                        "\"%s\" option requires an additional argument\n",
                        key);
         exit(EXIT_FAILURE);
      }

      spacetype = nextArg;

      return_value = TRUE;
   }

   /* Copy the strings */
   for (idim=0; idim < WORLD_NDIMS; idim++) {
      (void) strncpy(volume_def->spacetype[idim], spacetype, 
                     MI_MAX_ATTSTR_LEN);
      volume_def->spacetype[idim][MI_MAX_ATTSTR_LEN-1] = '\0';
   }

   return return_value;

}

/* ----------------------------- MNI Header -----------------------------------
@NAME       : set_units
@INPUT      : dst - Pointer to client data from argument table
              key - argument key
              nextArg - argument following key
@OUTPUT     : (nothing) 
@RETURNS    : TRUE if nextArg should be discarded, FALSE otherwise
@DESCRIPTION: Routine called by ParseArgv to set the units of the
              output sampling.
@METHOD     : 
@GLOBALS    : 
@CALLS      : 
@CREATED    : December 12, 1995 (Peter Neelin)
@MODIFIED   : 
---------------------------------------------------------------------------- */
static int set_units(char *dst, char *key, char *nextArg)
     /* ARGSUSED */
{
   Volume_Definition *volume_def;
   char *units;
   int idim;

   /* Get pointer to client data */
   volume_def = (Volume_Definition *) dst;

   /* Check for following argument */
   if (nextArg == NULL) {
      (void) fprintf(stderr, 
                     "\"%s\" option requires an additional argument\n",
                     key);
      exit(EXIT_FAILURE);
   }
   units = nextArg;

   /* Copy the strings */
   for (idim=0; idim < WORLD_NDIMS; idim++) {
      (void) strncpy(volume_def->units[idim], units, MI_MAX_ATTSTR_LEN);
      volume_def->units[idim][MI_MAX_ATTSTR_LEN-1] = '\0';
   }

   return TRUE;

}

/* ----------------------------- MNI Header -----------------------------------
@NAME       : get_axis_order
@INPUT      : dst - Pointer to client data from argument table
              key - argument key
              nextArg - argument following key
@OUTPUT     : (nothing) 
@RETURNS    : FALSE so that ParseArgv will not discard nextArg
@DESCRIPTION: Routine called by ParseArgv to get the axis order
@METHOD     : 
@GLOBALS    : 
@CALLS      : 
@CREATED    : February 15, 1993 (Peter Neelin)
@MODIFIED   : 
---------------------------------------------------------------------------- */
static int get_axis_order(char *dst, char *key, char *nextArg)
     /* ARGSUSED */
{
   Volume_Definition *volume_def;

   /* Get pointer to client data */
   volume_def = (Volume_Definition *) dst;

   /* Check key for order */
   if (strcmp(key, "-transverse") == 0) {
      volume_def->axes[Z] = 0;
      volume_def->axes[Y] = 1;
      volume_def->axes[X] = 2;
   }
   else if (strcmp(key, "-sagittal") == 0) {
      volume_def->axes[X] = 0;
      volume_def->axes[Z] = 1;
      volume_def->axes[Y] = 2;
   }
   else if (strcmp(key, "-coronal") == 0) {
      volume_def->axes[Y] = 0;
      volume_def->axes[Z] = 1;
      volume_def->axes[X] = 2;
   }

   return FALSE;
}

/* ----------------------------- MNI Header -----------------------------------
@NAME       : get_fillvalue
@INPUT      : dst - Pointer to client data from argument table
              key - argument key
              nextArg - argument following key
@OUTPUT     : (nothing) 
@RETURNS    : FALSE so that ParseArgv will not discard nextArg
@DESCRIPTION: Routine called by ParseArgv to set the fill value
@METHOD     : 
@GLOBALS    : 
@CALLS      : 
@CREATED    : February 15, 1993 (Peter Neelin)
@MODIFIED   : 
---------------------------------------------------------------------------- */
static int get_fillvalue(char *dst, char *key, char *nextArg)
     /* ARGSUSED */
{
   double *dptr;

   /* Get pointer to client data */
   dptr = (double *) dst;

   /* Check key for fill value to set */
   if (strcmp(key, "-fill") == 0) {
      *dptr = -DBL_MAX;
   }
   else if (strcmp(key, "-nofill") == 0) {
      *dptr = FILL_DEFAULT;
   }

   return FALSE;
}