1 /* ----------------------------- MNI Header -----------------------------------
2 @NAME : mincresample
3 @INPUT : argc, argv - command line arguments
4 @OUTPUT : (none)
5 @RETURNS : error status
6 @DESCRIPTION: Program to resample a minc file along different spatial
7 coordinate axes.
8 @METHOD :
9 @GLOBALS :
10 @CALLS :
11 @CREATED : February 8, 1993 (Peter Neelin)
12 @MODIFIED :
13 * $Log: mincresample.c,v $
14 * Revision 6.23 2008-01-17 02:33:02 rotor
15 * * removed all rcsids
16 * * removed a bunch of ^L's that somehow crept in
17 * * removed old (and outdated) BUGS file
18 *
19 * Revision 6.22 2008/01/13 09:38:54 stever
20 * Avoid compiler warnings about functions and variables that are defined
21 * but not used. Remove some such functions and variables,
22 * conditionalize some, and move static declarations out of header files
23 * into C files.
24 *
25 * Revision 6.21 2008/01/12 19:08:15 stever
26 * Add __attribute__ ((unused)) to all rcsid variables.
27 *
28 * Revision 6.20 2006/07/28 18:19:46 baghdadi
29 * *** empty log message ***
30 *
31 * Revision 6.19 2006/07/28 17:54:57 baghdadi
32 * added vector dimension to list of excluded files (minc2.
33 *
34 * Revision 6.18 2005/08/26 21:07:17 bert
35 * Use #if rather than #ifdef with MINC2 symbol, and be sure to include config.h whereever MINC2 is used
36 *
37 * Revision 6.17 2005/07/13 21:34:24 bert
38 * Add sinc interpolant (ported from 1.X branch)
39 *
40 * Revision 6.16 2004/11/01 22:38:39 bert
41 * Eliminate all references to minc_def.h
42 *
43 * Revision 6.15 2004/04/30 19:53:40 bert
44 * Remove unused variable
45 *
46 * Revision 6.14 2004/04/30 18:52:49 bert
47 * Remove some unused variables
48 *
49 * Revision 6.13 2004/04/27 15:31:20 bert
50 * Added -2 option
51 *
52 * Revision 6.12 2003/09/18 15:01:33 bert
53 * Removed unnecessary brackets from initializer
54 *
55 * Revision 6.11 2002/11/06 13:32:23 jason
56 * Fixes to mincresample: setting the interpolation type is now done
57 * through an enum rather than function pointers.
58 *
59 * Revision 6.10 2002/10/30 13:53:02 jason
60 * Added a ARGV_LONG argument type to ParseArgv. Used that type for the nelements variable in mincresample
61 *
62 * Revision 6.9 2001/08/24 19:12:50 neelin
63 * Re-ordered variables so that image variable is last. This fixes a problem
64 * with an uninitialized processing variable that shows up when mincresample
65 * is linked with netcdf 3.5.0 and an older program linked with 2.3.2 dies
66 * when reading the mincresample output file.
67 * It also allows use of large volumes with 64-bit minc: all variables
68 * must have 32-bit offset, but one (the image) can extend beyond the
69 * 32-bit limit.
70 *
71 * Revision 6.8 2001/08/16 13:32:39 neelin
72 * Partial fix for valid_range of different type from image (problems
73 * arising from double to float conversion/rounding). NOT COMPLETE.
74 *
75 * Revision 6.7 2001/04/24 13:38:45 neelin
76 * Replaced NC_NAT with MI_ORIGINAL_TYPE.
77 *
78 * Revision 6.6 2001/04/17 18:40:22 neelin
79 * Modifications to work with NetCDF 3.x
80 * In particular, changed NC_LONG to NC_INT (and corresponding longs to ints).
81 * Changed NC_UNSPECIFIED to NC_NAT.
82 * A few fixes to the configure script.
83 *
84 * Revision 6.5 2001/04/10 12:44:34 neelin
85 * Re-introduced fix for fillvalues that was lost with last version
86 *
87 * Revision 6.4 2001/04/02 14:58:09 neelin
88 * Added -keep_real_range option to prevent rescaling of slices on output
89 *
90 * Revision 6.2 1999/10/19 14:45:27 neelin
91 * Fixed Log subsitutions for CVS
92 *
93 * Revision 6.1 1998/02/19 15:04:24 neelin
94 * Minor bug fixes.
95 *
96 * Revision 6.0 1997/09/12 13:23:21 neelin
97 * Release of minc version 0.6
98 *
99 * Revision 5.0 1997/08/21 13:24:22 neelin
100 * Release of minc version 0.5
101 *
102 * Revision 4.1 1997/08/13 15:41:12 neelin
103 * Fixed initialization problem that caused crashing under Linux.
104 *
105 * Revision 4.0 1997/05/07 19:59:42 neelin
106 * Release of minc version 0.4
107 *
108 * Revision 3.5 1996/01/31 15:22:02 neelin
109 * Fixed bug in transformation of input sampling.
110 *
111 * Revision 3.4 1996/01/30 14:10:24 neelin
112 * Added check for -transformation without -like, -tfm_input_sampling or
113 * -use_input_sampling because of change in behaviour.
114 *
115 * Revision 3.3 1995/12/12 19:15:35 neelin
116 * Added -spacetype, -talairach and -units options.
117 *
118 * Revision 3.2 1995/11/21 14:13:20 neelin
119 * Transform input sampling with transformation and use this as default.
120 * Added -tfm_input_sampling to specify above option.
121 * Added -use_input_sampling to get old behaviour (no longer the default).
122 * Added -origin option (to specify coordinate instead of start values).
123 * Added -standard_sampling option (to set standard values of start, step
124 * and direction cosines).
125 * Added -invert_transformation option.
126 *
127 * Revision 3.1 1995/11/07 15:04:02 neelin
128 * Modified argument parsing so that only one pass is done.
129 *
130 * Revision 3.0 1995/05/15 19:30:57 neelin
131 * Release of minc version 0.3
132 *
133 * Revision 2.3 1995/05/05 19:11:05 neelin
134 * Modified call to input_transform.
135 *
136 * Revision 2.2 1995/02/09 14:05:51 neelin
137 * Mods to make irix 5 lint happy.
138 *
139 * Revision 2.1 1995/02/08 19:31:47 neelin
140 * Moved ARGSUSED statements for irix 5 lint.
141 *
142 * Revision 2.0 1994/09/28 10:32:46 neelin
143 * Release of minc version 0.2
144 *
145 * Revision 1.16 94/09/28 10:32:33 neelin
146 * Pre-release
147 *
148 * Revision 1.15 94/03/17 14:12:09 neelin
149 * Exit with failure if no argument given for -transformation or -like.
150 * .,
151 *
152 * Revision 1.14 94/03/15 16:44:21 neelin
153 * Changed default from -clobber to -noclobber.
154 *
155 * Revision 1.13 93/11/04 15:13:13 neelin
156 * Added support for irregularly spaced dimensions.
157 *
158 * Revision 1.12 93/11/03 14:32:44 neelin
159 * Turn off fill for output file.
160 *
161 * Revision 1.11 93/11/03 12:32:17 neelin
162 * Change ncopen, nccreate and ncclose to miopen, micreate and miclose.
163 *
164 * Revision 1.10 93/11/02 11:23:06 neelin
165 * Handle imagemax/min potentially varying over slices (for vector data, etc.)
166 *
167 * Revision 1.9 93/10/12 12:47:50 neelin
168 * Use volume_io.h instead of def_mni.h
169 *
170 * Revision 1.8 93/09/16 09:56:36 neelin
171 * Added use of open_file_with_default_suffix in get_transformation to
172 * append appropriate suffix for xfm files.
173 *
174 * Revision 1.7 93/08/11 14:28:19 neelin
175 * Modified get_arginfo and check_imageminmax to modify type of volume (not
176 * file) so that output volume gets the input volume type by default when
177 * an icv is used on input.
178 *
179 * Revision 1.6 93/08/11 13:27:59 neelin
180 * Converted to use Dave MacDonald's General_transform code.
181 * Fixed bug in get_slice - for non-linear transformations coord was
182 * transformed, then used again as a starting coordinate.
183 * Handle files that have image-max/min that doesn't vary over slices.
184 * Handle files that have image-max/min varying over row/cols.
185 * Allow volume to extend to voxel edge for -nearest_neighbour interpolation.
186 * Handle out-of-range values (-fill values from a previous mincresample, for
187 * example).
188 * Save transformation file as a string attribute to processing variable.
189 *
190 @COPYRIGHT :
191 Copyright 1993 Peter Neelin, McConnell Brain Imaging Centre,
192 Montreal Neurological Institute, McGill University.
193 Permission to use, copy, modify, and distribute this
194 software and its documentation for any purpose and without
195 fee is hereby granted, provided that the above copyright
196 notice appear in all copies. The author and McGill University
197 make no representations about the suitability of this
198 software for any purpose. It is provided "as is" without
199 express or implied warranty.
200 ---------------------------------------------------------------------------- */
201
202 #if HAVE_CONFIG_H
203 #include "config.h"
204 #endif
205
206 #include <stdlib.h>
207 #include <stdio.h>
208 #include <float.h>
209 #include <limits.h>
210 #include <string.h>
211 #include <math.h>
212 #include <minc.h>
213 #include <ParseArgv.h>
214 #include <time_stamp.h>
215 #include <volume_io.h>
216 #include <convert_origin_to_start.h>
217 #include "mincresample.h"
218
219 /* Kludge to catch use of -transformation without -like and without
220 -tfm_input_sampling or -use_input_sampling */
221 #define TRANSFORM_CHANGE_KLUDGE
222 #ifdef TRANSFORM_CHANGE_KLUDGE
223 int Specified_like = FALSE;
224 int Specified_transform = FALSE;
225 #define SAMPLING_ACTION_NOT_SET (-1)
226 #endif
227
228 static void get_arginfo(int argc, char *argv[],
229 Program_Flags *program_flags,
230 VVolume *in_vol, VVolume *out_vol,
231 General_transform *transformation);
232 static void check_imageminmax(File_Info *fp, Volume_Data *volume);
233 static void get_file_info(char *filename, int initialized_volume_def,
234 Volume_Definition *volume_def,
235 File_Info *file_info);
236 static void get_args_volume_def(Volume_Definition *input_volume_def,
237 Volume_Definition *args_volume_def);
238 static void transform_volume_def(Transform_Info *transform_info,
239 Volume_Definition *input_volume_def,
240 Volume_Definition *transformed_volume_def);
241 static int is_zero_vector(double vector[]);
242 static void normalize_vector(double vector[]);
243 static void create_output_file(char *filename, int clobber,
244 Volume_Definition *volume_def,
245 File_Info *in_file,
246 File_Info *out_file,
247 char *tm_stamp,
248 Transform_Info *transform_info);
249 static void get_voxel_to_world_transf(Volume_Definition *volume_def,
250 General_transform *voxel_to_world);
251 static void irregular_transform_function(void *user_data,
252 Real x,
253 Real y,
254 Real z,
255 Real *x_trans,
256 Real *y_trans,
257 Real *z_trans);
258 static void irregular_inverse_transform_function(void *user_data,
259 Real x,
260 Real y,
261 Real z,
262 Real *x_trans,
263 Real *y_trans,
264 Real *z_trans);
265 static double get_default_range(char *what, nc_type datatype, int is_signed);
266 static void finish_up(VVolume *in_vol, VVolume *out_vol);
267 static int get_transformation(char *dst, char *key, char *nextArg);
268 static int get_model_file(char *dst, char *key, char *nextArg);
269 static int set_standard_sampling(char *dst, char *key, char *nextArg);
270 static int set_spacetype(char *dst, char *key, char *nextArg);
271 static int set_units(char *dst, char *key, char *nextArg);
272 static int get_axis_order(char *dst, char *key, char *nextArg);
273 static int get_fillvalue(char *dst, char *key, char *nextArg);
274
275 /* Main program */
276
main(int argc,char * argv[])277 int main(int argc, char *argv[])
278 {
279 VVolume in_vol_struct, out_vol_struct;
280 VVolume *in_vol = &in_vol_struct, *out_vol = &out_vol_struct;
281 General_transform transformation;
282 Program_Flags program_flags;
283
284 /* Get argument information */
285 get_arginfo(argc, argv, &program_flags, in_vol, out_vol, &transformation);
286
287 /* Do the resampling */
288 resample_volumes(&program_flags, in_vol, out_vol, &transformation);
289
290 /* Finish up */
291 finish_up(in_vol, out_vol);
292
293 exit(EXIT_SUCCESS);
294 }
295
296 /* ----------------------------- MNI Header -----------------------------------
297 @NAME : get_arginfo
298 @INPUT : argc - number of command-line arguments
299 argv - command-line arguments
300 @OUTPUT : program_flags - data for program execution
301 in_vol - description of input volume.
302 out_vol - description of output volume.
303 transformation - description of world transformation
304 @RETURNS : (nothing)
305 @DESCRIPTION: Routine to get information from arguments about input and
306 output files and transfomation. Sets up all structures
307 completely (including allocating space for data).
308 @METHOD :
309 @GLOBALS :
310 @CALLS :
311 @CREATED : February 8, 1993 (Peter Neelin)
312 @MODIFIED :
313 ---------------------------------------------------------------------------- */
get_arginfo(int argc,char * argv[],Program_Flags * program_flags,VVolume * in_vol,VVolume * out_vol,General_transform * transformation)314 static void get_arginfo(int argc, char *argv[],
315 Program_Flags *program_flags,
316 VVolume *in_vol, VVolume *out_vol,
317 General_transform *transformation)
318 {
319 /* Argument parsing information */
320 #ifdef TRANSFORM_CHANGE_KLUDGE
321 static int transform_input_sampling = SAMPLING_ACTION_NOT_SET;
322 #else
323 static int transform_input_sampling = TRUE;
324 #endif
325 static Arg_Data args={
326 FALSE, /* Clobber */
327 FALSE, /* Keep scale */
328 MI_ORIGINAL_TYPE, /* Flag that type not set */
329 INT_MIN, /* Flag that is_signed has not been set */
330 {-DBL_MAX, -DBL_MAX}, /* Flag that range not set */
331 FILL_DEFAULT, /* Flag indicating that fillvalue not set */
332 {NO_VALUE, NO_VALUE, NO_VALUE}, /* Flag indicating that origin not set */
333 {TRUE}, /* Verbose */
334 TRILINEAR, /* use trilinear interpolation by default */
335 {FALSE, NULL, NULL, 0, NULL}, /* Transformation info is empty at start.
336 Transformation must be set before invoking
337 argument parsing */
338 { /* Use flags to indicate that values not set */
339 {NO_AXIS, NO_AXIS, NO_AXIS}, /* Axis order */
340 {0, 0, 0}, /* nelements */
341 {0.0, 0.0, 0.0}, /* step */
342 {NO_VALUE, NO_VALUE, NO_VALUE}, /* start */
343 {{NO_VALUE, NO_VALUE, NO_VALUE}, /* Default direction cosines */
344 {NO_VALUE, NO_VALUE, NO_VALUE},
345 {NO_VALUE, NO_VALUE, NO_VALUE}},
346 {NULL, NULL, NULL}, /* Pointers to coordinate arrays */
347 {"", "", ""}, /* units */
348 {"", "", ""} /* spacetype */
349 },
350 FALSE /* MINC 2.0 format? */
351 };
352
353 static ArgvInfo argTable[] = {
354 #if MINC2
355 {"-2", ARGV_CONSTANT, (char *) TRUE,
356 (char *) &args.v2format,
357 "Produce a MINC 2.0 format output file."},
358 #endif /* MINC2 */
359 {"-clobber", ARGV_CONSTANT, (char *) TRUE,
360 (char *) &args.clobber,
361 "Overwrite existing file."},
362 {"-noclobber", ARGV_CONSTANT, (char *) FALSE,
363 (char *) &args.clobber,
364 "Do not overwrite existing file (default)."},
365 {"-verbose", ARGV_CONSTANT, (char *) TRUE,
366 (char *) &args.flags.verbose,
367 "Print out log messages as processing is being done (default).\n"},
368 {"-quiet", ARGV_CONSTANT, (char *) FALSE,
369 (char *) &args.flags.verbose,
370 "Do not print out any log messages.\n"},
371 {"-transformation", ARGV_FUNC, (char *) get_transformation,
372 (char *) &args.transform_info,
373 "File giving world transformation. (Default = identity)."},
374 {"-invert_transformation", ARGV_CONSTANT, (char *) TRUE,
375 (char *) &args.transform_info.invert_transform,
376 "Invert the transformation before using it.\n"},
377 {"-noinvert_transformation", ARGV_CONSTANT, (char *) FALSE,
378 (char *) &args.transform_info.invert_transform,
379 "Do not invert the transformation (default).\n"},
380 {"-tfm_input_sampling", ARGV_CONSTANT, (char *) TRUE,
381 (char *) &transform_input_sampling,
382 "Transform the input sampling with the transform (default).\n"},
383 {"-use_input_sampling", ARGV_CONSTANT, (char *) FALSE,
384 (char *) &transform_input_sampling,
385 "Use the input sampling without transforming (old behaviour).\n"},
386 {"-like", ARGV_FUNC, (char *) get_model_file,
387 (char *) &args.volume_def,
388 "Specifies a model file for the resampling."},
389 {"-standard_sampling", ARGV_FUNC, (char *) set_standard_sampling,
390 (char *) &args.volume_def,
391 "Set the sampling to standard values (step, start and dircos)."},
392 {"-spacetype", ARGV_FUNC, (char *) set_spacetype,
393 (char *) &args.volume_def,
394 "Set the spacetype attribute to a specified string."},
395 {"-talairach", ARGV_FUNC, (char *) set_spacetype,
396 (char *) &args.volume_def,
397 "Output is in Talairach space."},
398 {"-units", ARGV_FUNC, (char *) set_units,
399 (char *) &args.volume_def,
400 "Specify the units of the output sampling."},
401 {"-nelements", ARGV_LONG, (char *) 3,
402 (char *) args.volume_def.nelements,
403 "Number of elements along each dimension (X, Y, Z)"},
404 {"-xnelements", ARGV_LONG, (char *) 0,
405 (char *) &args.volume_def.nelements[X],
406 "Number of elements along the X dimension"},
407 {"-ynelements", ARGV_LONG, (char *) 0,
408 (char *) &args.volume_def.nelements[Y],
409 "Number of elements along the Y dimension"},
410 {"-znelements", ARGV_LONG, (char *) 0,
411 (char *) &args.volume_def.nelements[Z],
412 "Number of elements along the Z dimension"},
413 {"-size", ARGV_LONG, (char *) 3,
414 (char *) args.volume_def.nelements,
415 "synonym for -nelements)"},
416 {"-xsize", ARGV_LONG, (char *) 0,
417 (char *) &args.volume_def.nelements[X],
418 "synonym for -xnelements"},
419 {"-ysize", ARGV_LONG, (char *) 0,
420 (char *) &args.volume_def.nelements[Y],
421 "synonym for -ynelements"},
422 {"-zsize", ARGV_LONG, (char *) 0,
423 (char *) &args.volume_def.nelements[Z],
424 "synonym for -ynelements"},
425 {"-step", ARGV_FLOAT, (char *) 3,
426 (char *) args.volume_def.step,
427 "Step size along each dimension (X, Y, Z)"},
428 {"-xstep", ARGV_FLOAT, (char *) 0,
429 (char *) &args.volume_def.step[X],
430 "Step size along the X dimension"},
431 {"-ystep", ARGV_FLOAT, (char *) 0,
432 (char *) &args.volume_def.step[Y],
433 "Step size along the Y dimension"},
434 {"-zstep", ARGV_FLOAT, (char *) 0,
435 (char *) &args.volume_def.step[Z],
436 "Step size along the Z dimension"},
437 {"-start", ARGV_FLOAT, (char *) 3,
438 (char *) args.volume_def.start,
439 "Start point along each dimension (X, Y, Z)"},
440 {"-xstart", ARGV_FLOAT, (char *) 0,
441 (char *) &args.volume_def.start[X],
442 "Start point along the X dimension"},
443 {"-ystart", ARGV_FLOAT, (char *) 0,
444 (char *) &args.volume_def.start[Y],
445 "Start point along the Y dimension"},
446 {"-zstart", ARGV_FLOAT, (char *) 0,
447 (char *) &args.volume_def.start[Z],
448 "Start point along the Z dimension"},
449 {"-dircos", ARGV_FLOAT, (char *) 9,
450 (char *) args.volume_def.dircos,
451 "Direction cosines along each dimension (X, Y, Z)"},
452 {"-xdircos", ARGV_FLOAT, (char *) 3,
453 (char *) args.volume_def.dircos[X],
454 "Direction cosines along the X dimension"},
455 {"-ydircos", ARGV_FLOAT, (char *) 3,
456 (char *) args.volume_def.dircos[Y],
457 "Direction cosines along the Y dimension"},
458 {"-zdircos", ARGV_FLOAT, (char *) 3,
459 (char *) args.volume_def.dircos[Z],
460 "Direction cosines along the Z dimension"},
461 {"-origin", ARGV_FLOAT, (char *) 3, (char *) args.origin,
462 "Origin of first pixel in 3D space"},
463 {"-transverse", ARGV_FUNC, (char *) get_axis_order,
464 (char *) &args.volume_def,
465 "Write out transverse slices"},
466 {"-sagittal", ARGV_FUNC, (char *) get_axis_order,
467 (char *) &args.volume_def,
468 "Write out sagittal slices"},
469 {"-coronal", ARGV_FUNC, (char *) get_axis_order,
470 (char *) &args.volume_def,
471 "Write out coronal slices"},
472 {"-byte", ARGV_CONSTANT, (char *) NC_BYTE, (char *) &args.datatype,
473 "Write out byte data"},
474 {"-short", ARGV_CONSTANT, (char *) NC_SHORT, (char *) &args.datatype,
475 "Write out short integer data"},
476 {"-int", ARGV_CONSTANT, (char *) NC_INT, (char *) &args.datatype,
477 "Write out 32-bit integer data"},
478 {"-long", ARGV_CONSTANT, (char *) NC_INT, (char *) &args.datatype,
479 "Superseded by -int"},
480 {"-float", ARGV_CONSTANT, (char *) NC_FLOAT, (char *) &args.datatype,
481 "Write out single-precision floating-point data"},
482 {"-double", ARGV_CONSTANT, (char *) NC_DOUBLE, (char *) &args.datatype,
483 "Write out double-precision floating-point data"},
484 {"-signed", ARGV_CONSTANT, (char *) TRUE, (char *) &args.is_signed,
485 "Write signed integer data"},
486 {"-unsigned", ARGV_CONSTANT, (char *) FALSE, (char *) &args.is_signed,
487 "Write unsigned integer data"},
488 {"-range", ARGV_FLOAT, (char *) 2, (char *) args.vrange,
489 "Valid range for output data"},
490 {"-keep_real_range", ARGV_CONSTANT, (char *) TRUE,
491 (char *) &args.keep_real_range,
492 "Keep the real scale of the input volume"},
493 {"-nokeep_real_range", ARGV_CONSTANT, (char *) FALSE,
494 (char *) &args.keep_real_range,
495 "Do not keep the real scale of the data (default)"},
496 {"-nofill", ARGV_FUNC, (char *) get_fillvalue,
497 (char *) &args.fillvalue,
498 "Use value zero for points outside of input volume"},
499 {"-fill", ARGV_FUNC, (char *) get_fillvalue,
500 (char *) &args.fillvalue,
501 "Use a fill value for points outside of input volume"},
502 {"-fillvalue", ARGV_FLOAT, (char *) 0,
503 (char *) &args.fillvalue,
504 "Specify a fill value for points outside of input volume"},
505 {"-trilinear", ARGV_CONSTANT, (char *) TRILINEAR,
506 (char *) &args.interpolant_type,
507 "Do trilinear interpolation"},
508 {"-tricubic", ARGV_CONSTANT, (char *) TRICUBIC,
509 (char *) &args.interpolant_type,
510 "Do tricubic interpolation"},
511 {"-nearest_neighbour", ARGV_CONSTANT,
512 (char *) N_NEIGHBOUR,
513 (char *) &args.interpolant_type,
514 "Do nearest neighbour interpolation"},
515 {"-sinc", ARGV_CONSTANT,
516 (char *) WINDOWED_SINC,
517 (char *) &args.interpolant_type,
518 "Do windowed sinc interpolation"},
519 {"-width", ARGV_INT,
520 (char *) 1,
521 (char *) &sinc_half_width,
522 "Set half-width of sinc window (1-10)" },
523 {"-hanning", ARGV_CONSTANT,
524 (char *) SINC_WINDOW_HANNING,
525 (char *) &sinc_window_type,
526 "Set sinc window type to Hanning"},
527 {"-hamming", ARGV_CONSTANT,
528 (char *) SINC_WINDOW_HAMMING,
529 (char *) &sinc_window_type,
530 "Set sinc window type to Hamming"},
531 {NULL, ARGV_END, NULL, NULL, NULL}
532 };
533
534 /* Other variables */
535 int idim, index;
536 int out_vindex; /* Volume indices (0, 1 or 2) */
537 int out_findex; /* File indices (0 to ndims-1) */
538 long size, total_size;
539 char *infile, *outfile;
540 File_Info *fp;
541 char *tm_stamp, *pname;
542 Volume_Definition input_volume_def, transformed_volume_def;
543 General_transform input_transformation;
544 int cflags;
545
546 /* Initialize the transformation to identity */
547 create_linear_transform(&input_transformation, NULL);
548 args.transform_info.transformation = &input_transformation;
549
550 /* Get the time stamp */
551 tm_stamp = time_stamp(argc, argv);
552
553 /* Save the program name */
554 pname=argv[0];
555
556 /* Call ParseArgv */
557 if (ParseArgv(&argc, argv, argTable, 0) || (argc!=3)) {
558 (void) fprintf(stderr,
559 "\nUsage: %s [<options>] <infile> <outfile>\n", pname);
560 (void) fprintf(stderr,
561 " %s [-help]\n\n", pname);
562 exit(EXIT_FAILURE);
563 }
564 infile = argv[1];
565 outfile = argv[2];
566
567 #ifdef TRANSFORM_CHANGE_KLUDGE
568 if (Specified_transform &&
569 !Specified_like &&
570 (transform_input_sampling == SAMPLING_ACTION_NOT_SET)) {
571 (void) fprintf(stderr,
572 "Use -like, -tfm_input_sampling or "
573 "-use_input_sampling with -transformation\n");
574 exit(EXIT_FAILURE);
575 }
576 if (transform_input_sampling == SAMPLING_ACTION_NOT_SET) {
577 transform_input_sampling = TRUE;
578 }
579 #endif
580
581 /* Check for an inverted transform. This looks backwards because we
582 normally invert the transform. */
583 if (args.transform_info.invert_transform) {
584 copy_general_transform(args.transform_info.transformation,
585 transformation);
586 }
587 else {
588 create_inverse_general_transform(args.transform_info.transformation,
589 transformation);
590 }
591 args.transform_info.transformation = transformation;
592
593 /* Get rid of the input transformation */
594 delete_general_transform(&input_transformation);
595
596 /* Check input file for default argument information */
597 in_vol->file = malloc(sizeof(File_Info));
598 get_file_info(infile, FALSE, &input_volume_def, in_vol->file);
599 transform_volume_def((transform_input_sampling ?
600 &args.transform_info : NULL),
601 &input_volume_def,
602 &transformed_volume_def);
603 get_args_volume_def(&transformed_volume_def, &args.volume_def);
604
605 /* Check that direction cosines are normalized and look for origin
606 option */
607 for (idim=0; idim < WORLD_NDIMS; idim++) {
608 normalize_vector(args.volume_def.dircos[idim]);
609 if (is_zero_vector(args.volume_def.dircos[idim])) {
610 (void) fprintf(stderr, "Bad direction cosines.\n");
611 exit(EXIT_FAILURE);
612 }
613 }
614 if (args.origin[0] != NO_VALUE) {
615 if (convert_origin_to_start(args.origin,
616 args.volume_def.dircos[XCOORD],
617 args.volume_def.dircos[YCOORD],
618 args.volume_def.dircos[ZCOORD],
619 args.volume_def.start) != 0) {
620 (void) fprintf(stderr, "Error converting origin to start value: ");
621 (void) fprintf(stderr, "Bad direction cosines.\n");
622 exit(EXIT_FAILURE);
623 }
624 }
625
626 /* Save the voxel_to_world transformation information */
627 in_vol->voxel_to_world = malloc(sizeof(General_transform));
628 in_vol->world_to_voxel = malloc(sizeof(General_transform));
629 get_voxel_to_world_transf(&input_volume_def, in_vol->voxel_to_world);
630 create_inverse_general_transform(in_vol->voxel_to_world,
631 in_vol->world_to_voxel);
632
633 /* Get input volume data information */
634 in_vol->slice = NULL;
635 in_vol->volume = malloc(sizeof(Volume_Data));
636 in_vol->volume->datatype = in_vol->file->datatype;
637 in_vol->volume->is_signed = in_vol->file->is_signed;
638 in_vol->volume->vrange[0] = in_vol->file->vrange[0];
639 in_vol->volume->vrange[1] = in_vol->file->vrange[1];
640 if (args.fillvalue == FILL_DEFAULT) {
641 in_vol->volume->fillvalue = 0.0;
642 in_vol->volume->use_fill = TRUE;
643 }
644 else {
645 in_vol->volume->fillvalue = args.fillvalue;
646 in_vol->volume->use_fill = (args.fillvalue != -DBL_MAX);
647 }
648
649 /* set the function pointer defining the type of interpolation */
650 switch (args.interpolant_type ) {
651 case TRICUBIC:
652 in_vol->volume->interpolant = tricubic_interpolant;
653 break;
654 case TRILINEAR:
655 in_vol->volume->interpolant = trilinear_interpolant;
656 break;
657 case N_NEIGHBOUR:
658 in_vol->volume->interpolant = nearest_neighbour_interpolant;
659 break;
660 case WINDOWED_SINC:
661 in_vol->volume->interpolant = windowed_sinc_interpolant;
662 if (sinc_half_width < SINC_HALF_WIDTH_MIN ||
663 sinc_half_width > SINC_HALF_WIDTH_MAX) {
664 fprintf(stderr, "Invalid sinc half-window size %d\n",
665 sinc_half_width);
666 exit(EXIT_FAILURE);
667 }
668 break;
669 default:
670 (void) fprintf(stderr, "Error determining interpolation type\n");
671 exit(EXIT_FAILURE);
672 }
673
674 /* Check min/max variables */
675 fp = in_vol->file;
676 fp->using_icv=FALSE;
677 if ((fp->maxid != MI_ERROR) && (fp->minid != MI_ERROR) &&
678 (fp->datatype!=NC_FLOAT) && (fp->datatype!=NC_DOUBLE)) {
679 check_imageminmax(fp, in_vol->volume);
680 }
681
682 /* Get space for volume data */
683 total_size = 1;
684 for (idim=0; idim < WORLD_NDIMS; idim++) {
685 index = input_volume_def.axes[idim];
686 size = input_volume_def.nelements[idim];
687 total_size *= size;
688 in_vol->volume->size[index] = size;
689 }
690 in_vol->volume->data = malloc((size_t) total_size *
691 nctypelen(in_vol->volume->datatype));
692
693 /* Get space for slice scale and offset */
694 in_vol->volume->scale =
695 malloc(sizeof(double) * in_vol->volume->size[SLC_AXIS]);
696 in_vol->volume->offset =
697 malloc(sizeof(double) * in_vol->volume->size[SLC_AXIS]);
698
699 /* Save the program flags */
700 *program_flags = args.flags;
701
702 /* Set the default output file datatype */
703 if (args.datatype == MI_ORIGINAL_TYPE)
704 args.datatype = in_vol->file->datatype;
705
706 /* Explicitly force output files to have regular spacing */
707 for (idim=0; idim < WORLD_NDIMS; idim++) {
708 if (args.volume_def.coords[idim] != NULL) {
709 free(args.volume_def.coords[idim]);
710 args.volume_def.coords[idim] = NULL;
711 }
712 }
713
714 /* Check to see if sign and range have been explicitly set. If not set
715 them now */
716 if (args.is_signed == INT_MIN) {
717 if (args.datatype == in_vol->file->datatype)
718 args.is_signed = in_vol->file->is_signed;
719 else
720 args.is_signed = (args.datatype != NC_BYTE);
721 }
722 if (args.vrange[0] == -DBL_MAX) {
723 if ((args.datatype == in_vol->file->datatype) &&
724 (args.is_signed == in_vol->file->is_signed)) {
725 args.vrange[0] = in_vol->file->vrange[0];
726 args.vrange[1] = in_vol->file->vrange[1];
727 }
728 else {
729 args.vrange[0] = get_default_range(MIvalid_min, args.datatype,
730 args.is_signed);
731 args.vrange[1] = get_default_range(MIvalid_max, args.datatype,
732 args.is_signed);
733 }
734 }
735
736 /* Set up the file description for the output file */
737 out_vol->file = malloc(sizeof(File_Info));
738 out_vol->file->ndims = in_vol->file->ndims;
739 out_vol->file->datatype = args.datatype;
740 out_vol->file->is_signed = args.is_signed;
741 out_vol->file->vrange[0] = args.vrange[0];
742 out_vol->file->vrange[1] = args.vrange[1];
743 for (idim=0; idim < out_vol->file->ndims; idim++) {
744 out_vol->file->nelements[idim] = in_vol->file->nelements[idim];
745 out_vol->file->world_axes[idim] = in_vol->file->world_axes[idim];
746 }
747 out_vol->file->keep_real_range = args.keep_real_range;
748
749 /* Get space for output slice */
750 out_vol->volume = NULL;
751 out_vol->slice = malloc(sizeof(Slice_Data));
752
753 /* Loop through list of axes, getting size of volume and slice */
754 total_size = 1;
755 for (idim=0; idim < WORLD_NDIMS; idim++) {
756
757 /* Get the index for input and output volumes */
758 out_vindex = args.volume_def.axes[idim]; /* 0, 1 or 2 */
759 out_findex = in_vol->file->indices[out_vindex]; /* 0 to ndims-1 */
760 size = args.volume_def.nelements[idim];
761
762 /* Update output axes and indices and nelements */
763 out_vol->file->nelements[out_findex] = size;
764 out_vol->file->world_axes[out_findex] = idim;
765 out_vol->file->axes[idim] = out_vindex;
766 out_vol->file->indices[out_vindex] = out_findex;
767
768 /* Update slice size */
769 if (out_vindex != 0) {
770 out_vol->slice->size[out_vindex-1] = size;
771 total_size *= size;
772 }
773 }
774 out_vol->slice->data = malloc((size_t) total_size * sizeof(double));
775
776 /* Create the output file */
777 if (args.clobber) {
778 cflags = NC_CLOBBER;
779 }
780 else {
781 cflags = NC_NOCLOBBER;
782 }
783 #if MINC2
784 if (args.v2format) {
785 cflags |= MI2_CREATE_V2;
786 }
787 #endif /* MINC2 */
788 create_output_file(outfile, cflags, &args.volume_def,
789 in_vol->file, out_vol->file,
790 tm_stamp, &args.transform_info);
791
792 /* Save the voxel_to_world transformation information */
793 out_vol->voxel_to_world = malloc(sizeof(General_transform));
794 out_vol->world_to_voxel = malloc(sizeof(General_transform));
795 get_voxel_to_world_transf(&args.volume_def, out_vol->voxel_to_world);
796 create_inverse_general_transform(out_vol->voxel_to_world,
797 out_vol->world_to_voxel);
798
799 /* Free the time stamp */
800 free(tm_stamp);
801
802 }
803
804 /* ----------------------------- MNI Header -----------------------------------
805 @NAME : check_imageminmax
806 @INPUT : fp - pointer to file description
807 volume - pointer to volume description
808 @OUTPUT :
809 @RETURNS : (nothing)
810 @DESCRIPTION: Routine to check that MIimagemax and MIimagemin do not vary
811 over volume rows and columns. If they do, set up an icv to
812 handle it.
813 @METHOD :
814 @GLOBALS :
815 @CALLS :
816 @CREATED : August 5, 1993 (Peter Neelin)
817 @MODIFIED :
818 ---------------------------------------------------------------------------- */
check_imageminmax(File_Info * fp,Volume_Data * volume)819 static void check_imageminmax(File_Info *fp, Volume_Data *volume)
820 {
821 int ndims, idim, dim[MAX_VAR_DIMS], imgdim[MAX_VAR_DIMS];
822 int ivar, varid;
823
824 /* Get MIimage dimensions */
825 (void) ncvarinq(fp->mincid, fp->imgid, NULL, NULL, &ndims, imgdim, NULL);
826
827 /* Check MIimagemax/min dimensions */
828 for (ivar=0; ivar<2; ivar++) {
829 varid = (ivar==0 ? fp->maxid : fp->minid);
830 (void) ncvarinq(fp->mincid, varid, NULL, NULL, &ndims, dim, NULL);
831 for (idim=0; idim < ndims; idim++) {
832 if ((dim[idim] == imgdim[fp->indices[ROW_AXIS]]) ||
833 (dim[idim] == imgdim[fp->indices[COL_AXIS]])) {
834 fp->using_icv = TRUE;
835 }
836 } /* End loop over MIimagemax/min dimensions */
837 } /* End loop over variables MIimagemax/min */
838
839 /* Set up an icv if needed to handle values varying over slice dims. */
840 if (fp->using_icv) {
841
842 /* Change type to floating point so that there is no loss of
843 precision (except possibly for long values). */
844 if (volume->datatype != NC_DOUBLE)
845 volume->datatype = NC_FLOAT;
846 volume->is_signed = TRUE;
847
848 /* Create the icv */
849 fp->icvid = miicv_create();
850 (void) miicv_setint(fp->icvid, MI_ICV_TYPE, volume->datatype);
851 (void) miicv_setstr(fp->icvid, MI_ICV_SIGN,
852 (volume->is_signed ? MI_SIGNED : MI_UNSIGNED));
853 (void) miicv_setint(fp->icvid, MI_ICV_DO_NORM, TRUE);
854 (void) miicv_setint(fp->icvid, MI_ICV_DO_FILLVALUE, TRUE);
855 (void) miicv_attach(fp->icvid, fp->mincid, fp->imgid);
856
857 /* Get max and min for doing valid range checking */
858 (void) miicv_inqdbl(fp->icvid, MI_ICV_NORM_MIN, &volume->vrange[0]);
859 (void) miicv_inqdbl(fp->icvid, MI_ICV_NORM_MAX, &volume->vrange[1]);
860
861 }
862
863 }
864
865 /* ----------------------------- MNI Header -----------------------------------
866 @NAME : get_file_info
867 @INPUT : filename - name of file to read
868 initialized_volume_def - if TRUE, then volume_def is taken
869 as being properly initialized and arrays are freed if
870 non-NULL. Otherwise arrays are not freed.
871 @OUTPUT : volume_def - description of volume
872 file_info - description of file
873 @RETURNS : (nothing)
874 @DESCRIPTION: Routine to get information about the volume definition of
875 a minc file.
876 @METHOD :
877 @GLOBALS :
878 @CALLS :
879 @CREATED : February 9, 1993 (Peter Neelin)
880 @MODIFIED :
881 ---------------------------------------------------------------------------- */
get_file_info(char * filename,int initialized_volume_def,Volume_Definition * volume_def,File_Info * file_info)882 static void get_file_info(char *filename, int initialized_volume_def,
883 Volume_Definition *volume_def,
884 File_Info *file_info)
885 {
886 int dim[MAX_VAR_DIMS], dimid;
887 int axis_counter, idim, jdim, cur_axis;
888 int varndims, vardim[MAX_VAR_DIMS];
889 long varstart, varcount, dimlength;
890 char attstr[MI_MAX_ATTSTR_LEN];
891 char dimname[MAX_NC_NAME];
892 enum {UNKNOWN, REGULAR, IRREGULAR} coord_spacing;
893
894 /* Open the minc file */
895 file_info->mincid = miopen(filename, NC_NOWRITE);
896 file_info->name = filename;
897
898 /* Get variable identifiers */
899 file_info->imgid = ncvarid(file_info->mincid, MIimage);
900 ncopts = 0;
901 file_info->maxid = ncvarid(file_info->mincid, MIimagemax);
902 file_info->minid = ncvarid(file_info->mincid, MIimagemin);
903 ncopts = NC_VERBOSE | NC_FATAL;
904
905 /* Get information about datatype dimensions of variable */
906 (void) miget_datatype(file_info->mincid, file_info->imgid,
907 &file_info->datatype, &file_info->is_signed);
908
909 /* Get valid max and min */
910 (void) miget_valid_range(file_info->mincid, file_info->imgid,
911 file_info->vrange);
912
913 /* Get information about dimensions */
914 (void) ncvarinq(file_info->mincid, file_info->imgid, NULL, NULL,
915 &file_info->ndims, dim, NULL);
916
917 /* Set variables for keeping track of spatial dimensions */
918 axis_counter = 0; /* Keeps track of values for axes */
919
920 /* Initialize volume definition variables */
921 for (idim=0; idim < WORLD_NDIMS; idim++) {
922 volume_def->axes[idim] = NO_AXIS;
923 volume_def->step[idim] = 1.0;
924 volume_def->start[idim] = 0.0;
925 for (jdim=0; jdim < WORLD_NDIMS; jdim++) {
926 if (jdim==idim)
927 volume_def->dircos[idim][jdim] = 1.0;
928 else
929 volume_def->dircos[idim][jdim] = 0.0;
930 }
931 if (initialized_volume_def && (volume_def->coords[idim] != NULL)) {
932 free(volume_def->coords[idim]);
933 }
934 volume_def->coords[idim] = NULL;
935 (void) strcpy(volume_def->units[idim], "mm");
936 (void) strcpy(volume_def->spacetype[idim], MI_NATIVE);
937 }
938
939 /* Loop through dimensions, getting dimension information */
940
941 for (idim=0; idim < file_info->ndims; idim++) {
942
943 /* Get size of dimension */
944 (void) ncdiminq(file_info->mincid, dim[idim], dimname,
945 &file_info->nelements[idim]);
946
947 /* Check variable name */
948 cur_axis = NO_AXIS;
949 if (strcmp(dimname, MIxspace)==0)
950 cur_axis = XAXIS;
951 else if (strcmp(dimname, MIyspace)==0)
952 cur_axis = YAXIS;
953 else if (strcmp(dimname, MIzspace)==0)
954 cur_axis = ZAXIS;
955
956 /* Save world axis info */
957 file_info->world_axes[idim] = cur_axis;
958
959 /* Check for spatial dimension */
960 if (cur_axis == NO_AXIS) continue;
961
962 /* Set axis */
963 if (volume_def->axes[cur_axis] != NO_AXIS) {
964 (void) fprintf(stderr, "Repeated spatial dimension %s in file %s.\n",
965 dimname, filename);
966 exit(EXIT_FAILURE);
967 }
968 volume_def->axes[cur_axis] = axis_counter++;
969
970 /* Save spatial axis specific info */
971 file_info->axes[cur_axis] = volume_def->axes[cur_axis];
972 file_info->indices[volume_def->axes[cur_axis]] = idim;
973 volume_def->nelements[cur_axis] = file_info->nelements[idim];
974
975 /* Check for existence of variable */
976 ncopts = 0;
977 dimid = ncvarid(file_info->mincid, dimname);
978 ncopts = NC_VERBOSE | NC_FATAL;
979 if (dimid == MI_ERROR) continue;
980
981 /* Get attributes from variable */
982 ncopts = 0;
983 (void) miattget1(file_info->mincid, dimid, MIstep,
984 NC_DOUBLE, &volume_def->step[cur_axis]);
985 if (volume_def->step[cur_axis] == 0.0)
986 volume_def->step[cur_axis] = 1.0;
987 (void) miattget1(file_info->mincid, dimid, MIstart,
988 NC_DOUBLE, &volume_def->start[cur_axis]);
989 (void) miattget(file_info->mincid, dimid, MIdirection_cosines,
990 NC_DOUBLE, WORLD_NDIMS,
991 volume_def->dircos[cur_axis], NULL);
992 (void) miattgetstr(file_info->mincid, dimid, MIunits,
993 MI_MAX_ATTSTR_LEN, volume_def->units[cur_axis]);
994 (void) miattgetstr(file_info->mincid, dimid, MIspacetype,
995 MI_MAX_ATTSTR_LEN, volume_def->spacetype[cur_axis]);
996 ncopts = NC_VERBOSE | NC_FATAL;
997
998 /* Normalize the direction cosine */
999 normalize_vector(volume_def->dircos[cur_axis]);
1000
1001 /* Look for irregular coordinates for dimension variable */
1002 ncopts = 0;
1003 coord_spacing = UNKNOWN;
1004 dimlength = volume_def->nelements[cur_axis];
1005 if (miattgetstr(file_info->mincid, dimid, MIspacing, MI_MAX_ATTSTR_LEN,
1006 attstr) != NULL) {
1007 if (strcmp(attstr, MI_IRREGULAR) == 0)
1008 coord_spacing = IRREGULAR;
1009 else if (strcmp(attstr, MI_REGULAR) == 0)
1010 coord_spacing = REGULAR;
1011 }
1012 if (ncvarinq(file_info->mincid, dimid, NULL, NULL,
1013 &varndims, vardim, NULL) == MI_ERROR) {
1014 ncopts = NC_VERBOSE | NC_FATAL;
1015 continue;
1016 }
1017 if ((coord_spacing != REGULAR) &&
1018 (varndims == 1) && (vardim[0] == dim[idim])) {
1019 coord_spacing = IRREGULAR;
1020 }
1021 if ((coord_spacing == UNKNOWN) || (dimlength <= 1)) {
1022 coord_spacing = REGULAR;
1023 }
1024 if (coord_spacing == IRREGULAR) {
1025 volume_def->coords[cur_axis] = malloc(sizeof(double) * dimlength);
1026 varstart = 0;
1027 varcount = dimlength;
1028 if (mivarget(file_info->mincid, dimid, &varstart, &varcount,
1029 NC_DOUBLE, MI_SIGNED, volume_def->coords[cur_axis])
1030 == MI_ERROR) {
1031 ncopts = NC_VERBOSE | NC_FATAL;
1032 free(volume_def->coords[cur_axis]);
1033 volume_def->coords[cur_axis] = NULL;
1034 continue;
1035 }
1036 volume_def->start[cur_axis] = volume_def->coords[cur_axis][0];
1037 if (dimlength > 1) {
1038 volume_def->step[cur_axis] =
1039 (volume_def->coords[cur_axis][dimlength-1] -
1040 volume_def->coords[cur_axis][0]) /
1041 (dimlength - 1);
1042 if (volume_def->step[cur_axis] == 0.0)
1043 volume_def->step[cur_axis] = 1.0;
1044 }
1045 }
1046 ncopts = NC_VERBOSE | NC_FATAL;
1047
1048 } /* End of loop over dimensions */
1049
1050 /* Check that we have the correct number of spatial dimensions */
1051 if (axis_counter != WORLD_NDIMS) {
1052 (void) fprintf(stderr,
1053 "Incorrect number of spatial dimensions in file %s.\n",
1054 filename);
1055 exit(EXIT_FAILURE);
1056 }
1057
1058 return;
1059 }
1060
1061 /* ----------------------------- MNI Header -----------------------------------
1062 @NAME : get_args_volume_def
1063 @INPUT : input_volume_def - description of input volume
1064 @OUTPUT : args_volume_def - description of new output volume
1065 @RETURNS : (nothing)
1066 @DESCRIPTION: Routine to copy appropriate information from input volume
1067 definition to output volume definition, overriding values
1068 not set on the command line.
1069 @METHOD :
1070 @GLOBALS :
1071 @CALLS :
1072 @CREATED : November 7, 1995 (Peter Neelin)
1073 @MODIFIED :
1074 ---------------------------------------------------------------------------- */
get_args_volume_def(Volume_Definition * input_volume_def,Volume_Definition * args_volume_def)1075 static void get_args_volume_def(Volume_Definition *input_volume_def,
1076 Volume_Definition *args_volume_def)
1077 {
1078 int idim, jdim;
1079
1080 /* Loop over dimensions */
1081 for (idim=0; idim < WORLD_NDIMS; idim++) {
1082
1083 /* Copy values, as needed */
1084 if (args_volume_def->axes[idim] == NO_AXIS)
1085 args_volume_def->axes[idim] = input_volume_def->axes[idim];
1086 if (args_volume_def->nelements[idim] == 0)
1087 args_volume_def->nelements[idim] = input_volume_def->nelements[idim];
1088 if (args_volume_def->step[idim] == 0.0)
1089 args_volume_def->step[idim] = input_volume_def->step[idim];
1090 if (args_volume_def->start[idim] == NO_VALUE)
1091 args_volume_def->start[idim] = input_volume_def->start[idim];
1092 if (args_volume_def->dircos[idim][0] == NO_VALUE) {
1093 for (jdim=0; jdim < WORLD_NDIMS; jdim++)
1094 args_volume_def->dircos[idim][jdim] =
1095 input_volume_def->dircos[idim][jdim];
1096 }
1097 if (strlen(args_volume_def->units[idim]) == 0)
1098 (void) strcpy(args_volume_def->units[idim],
1099 input_volume_def->units[idim]);
1100 if (strlen(args_volume_def->spacetype[idim]) == 0)
1101 (void) strcpy(args_volume_def->spacetype[idim],
1102 input_volume_def->spacetype[idim]);
1103
1104 }
1105 }
1106
1107 /* ----------------------------- MNI Header -----------------------------------
1108 @NAME : transform_volume_def
1109 @INPUT : transform_info - description of output to input transform
1110 (if NULL, then don't transform).
1111 input_volume_def - description of input volume
1112 @OUTPUT : transformed_volume_def - description of new output volume
1113 @RETURNS : (nothing)
1114 @DESCRIPTION: Routine to copy appropriate information from input volume
1115 definition to a new volume definition, after transformation
1116 from with the output to input transform.
1117 @METHOD :
1118 @GLOBALS :
1119 @CALLS :
1120 @CREATED : November 7, 1995 (Peter Neelin)
1121 @MODIFIED :
1122 ---------------------------------------------------------------------------- */
transform_volume_def(Transform_Info * transform_info,Volume_Definition * input_volume_def,Volume_Definition * transformed_volume_def)1123 static void transform_volume_def(Transform_Info *transform_info,
1124 Volume_Definition *input_volume_def,
1125 Volume_Definition *transformed_volume_def)
1126 {
1127 int idim, jdim;
1128 Coord_Vector origin = {0, 0, 0};
1129 Coord_Vector transformed_origin, vector;
1130 double length, nelements;
1131
1132 /* Copy the volume definition */
1133 *transformed_volume_def = *input_volume_def;
1134
1135 /* Modify things if the transform has changed */
1136 if ((transform_info != NULL) && (transform_info->file_name != NULL)) {
1137
1138 /* Set up origin vector */
1139 for (idim=0; idim < WORLD_NDIMS; idim++)
1140 for (jdim=0; jdim < WORLD_NDIMS; jdim++)
1141 origin[idim] += input_volume_def->start[jdim] *
1142 input_volume_def->dircos[jdim][idim];
1143
1144 /* Transform origin vector */
1145 DO_INVERSE_TRANSFORM(transformed_origin,
1146 transform_info->transformation, origin);
1147
1148 /* Loop over dimensions */
1149 for (idim=0; idim < WORLD_NDIMS; idim++) {
1150
1151 /* Get number of elements */
1152 nelements = input_volume_def->nelements[idim] - 1;
1153 if (nelements < 1) nelements = 1.0;
1154
1155 /* Transform whole axis */
1156 VECTOR_SCALAR_MULT(vector, input_volume_def->dircos[idim],
1157 nelements * input_volume_def->step[idim]);
1158 VECTOR_ADD(vector, vector, origin);
1159 DO_INVERSE_TRANSFORM(vector, transform_info->transformation, vector);
1160 VECTOR_DIFF(vector, vector, transformed_origin);
1161
1162 /* Calculate length of transformed axis - new step value */
1163 length = sqrt(vector[XCOORD]*vector[XCOORD] +
1164 vector[YCOORD]*vector[YCOORD] +
1165 vector[ZCOORD]*vector[ZCOORD]);
1166 transformed_volume_def->step[idim] = length / nelements;
1167
1168 /* Calculate new direction cosines */
1169 for (jdim=0; jdim < WORLD_NDIMS; jdim++) {
1170 transformed_volume_def->dircos[idim][jdim] =
1171 (length > 0.0 ? vector[jdim] / length :
1172 (jdim == idim ? 1.0 : 0.0));
1173 }
1174
1175 /* Make sure that direction cosines point the right way */
1176 if (transformed_volume_def->dircos[idim][idim] < 0.0) {
1177 VECTOR_SCALAR_MULT(transformed_volume_def->dircos[idim],
1178 transformed_volume_def->dircos[idim],
1179 -1.0);
1180 transformed_volume_def->step[idim] *= -1.0;
1181 }
1182 }
1183
1184 /* Calculate the start along each axis - check for bad dircos. */
1185 if (convert_origin_to_start(transformed_origin,
1186 transformed_volume_def->dircos[XCOORD],
1187 transformed_volume_def->dircos[YCOORD],
1188 transformed_volume_def->dircos[ZCOORD],
1189 transformed_volume_def->start) != 0) {
1190
1191 /* If dircos are bad, set them to default */
1192 for (idim=0; idim < WORLD_NDIMS; idim++) {
1193 for (jdim=0; jdim < WORLD_NDIMS; jdim++) {
1194 transformed_volume_def->dircos[idim][jdim] =
1195 ((idim==jdim) ? 1.0 : 0.0);
1196 }
1197 }
1198
1199 /* Try to convert point again */
1200 if (convert_origin_to_start(transformed_origin,
1201 transformed_volume_def->dircos[XCOORD],
1202 transformed_volume_def->dircos[YCOORD],
1203 transformed_volume_def->dircos[ZCOORD],
1204 transformed_volume_def->start) != 0) {
1205 (void) fprintf(stderr,
1206 "Serious problem converting origin to start!\n");
1207 exit(EXIT_FAILURE);
1208 }
1209 }
1210
1211 }
1212
1213 }
1214
1215 /* ----------------------------- MNI Header -----------------------------------
1216 @NAME : is_zero_vector
1217 @INPUT : vector - 3D vector
1218 @OUTPUT : (none)
1219 @RETURNS : 1 if vector has zero length, 0 otherwise
1220 @DESCRIPTION: Routine to check for a zero length vector.
1221 @METHOD :
1222 @GLOBALS :
1223 @CALLS :
1224 @CREATED : November 9, 1995 (Peter Neelin)
1225 @MODIFIED :
1226 ---------------------------------------------------------------------------- */
is_zero_vector(double vector[])1227 static int is_zero_vector(double vector[])
1228 {
1229 return ((vector[0] == 0.0) &&
1230 (vector[1] == 0.0) &&
1231 (vector[2] == 0.0));
1232 }
1233
1234 /* ----------------------------- MNI Header -----------------------------------
1235 @NAME : normalize_vector
1236 @INPUT : vector - 3D vector
1237 @OUTPUT : (none)
1238 @RETURNS : (nothing)
1239 @DESCRIPTION: Routine to normalize a vector
1240 @METHOD :
1241 @GLOBALS :
1242 @CALLS :
1243 @CREATED : November 9, 1995 (Peter Neelin)
1244 @MODIFIED :
1245 ---------------------------------------------------------------------------- */
normalize_vector(double vector[])1246 static void normalize_vector(double vector[])
1247 {
1248 int idim;
1249 double magnitude;
1250
1251 /* Normalize the direction cosine */
1252 magnitude = 0.0;
1253 for (idim=0; idim < WORLD_NDIMS; idim++) {
1254 magnitude += (vector[idim] * vector[idim]);
1255 }
1256 magnitude = sqrt(magnitude);
1257 if (magnitude > 0.0) {
1258 for (idim=0; idim < WORLD_NDIMS; idim++) {
1259 vector[idim] /= magnitude;
1260 }
1261 }
1262
1263 }
1264
1265 /* ----------------------------- MNI Header -----------------------------------
1266 @NAME : create_output_file
1267 @INPUT : filename - name of file to create
1268 cflags - flag creation flags
1269 volume_def - description of volume
1270 in_file - description of input file
1271 out_file - description of output file
1272 tm_stamp - string describing program invocation
1273 transformation - transformation to be applied to data
1274 @OUTPUT : (nothing)
1275 @RETURNS : (nothing)
1276 @DESCRIPTION: Routine to create an minc output file and set up its parameters
1277 properly.
1278 @METHOD :
1279 @GLOBALS :
1280 @CALLS :
1281 @CREATED : February 9, 1993 (Peter Neelin)
1282 @MODIFIED :
1283 ---------------------------------------------------------------------------- */
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)1284 static void create_output_file(char *filename, int cflags,
1285 Volume_Definition *volume_def,
1286 File_Info *in_file,
1287 File_Info *out_file,
1288 char *tm_stamp,
1289 Transform_Info *transform_info)
1290 {
1291 int ndims, in_dims[MAX_VAR_DIMS], out_dims[MAX_VAR_DIMS];
1292 int out_maxmin_dims[MAX_VAR_DIMS];
1293 char dimname[MAX_NC_NAME];
1294 int nmaxmin_dims, nimage_dims, axis, idim, dimid, varid, itrans;
1295 int att_length;
1296 nc_type datatype;
1297 int nexcluded, excluded_vars[10];
1298 char *string;
1299 int dim_exists, is_volume_dimension;
1300 int in_index, out_index;
1301
1302 /* Save the file name */
1303 out_file->name = filename;
1304
1305 /* Create the list of excluded variables */
1306 nexcluded = 0;
1307 ncopts = 0;
1308 if ((varid=ncvarid(in_file->mincid, MIxspace)) != MI_ERROR)
1309 excluded_vars[nexcluded++] = varid;
1310 if ((varid=ncvarid(in_file->mincid, MIyspace)) != MI_ERROR)
1311 excluded_vars[nexcluded++] = varid;
1312 if ((varid=ncvarid(in_file->mincid, MIzspace)) != MI_ERROR)
1313 excluded_vars[nexcluded++] = varid;
1314 if ((varid=ncvarid(in_file->mincid, MIimage)) != MI_ERROR)
1315 excluded_vars[nexcluded++] = varid;
1316 if ((varid=ncvarid(in_file->mincid, MIimagemax)) != MI_ERROR)
1317 excluded_vars[nexcluded++] = varid;
1318 if ((varid=ncvarid(in_file->mincid, MIimagemin)) != MI_ERROR)
1319 excluded_vars[nexcluded++] = varid;
1320 #if MINC2
1321 if ((varid=ncvarid(in_file->mincid, MIvector_dimension)) != MI_ERROR)
1322 excluded_vars[nexcluded++] = varid;
1323 #endif /* MINC2 */
1324 ncopts = NC_VERBOSE | NC_FATAL;
1325
1326 /* Create the file */
1327 out_file->mincid = micreate(filename, cflags);
1328
1329 /* Copy all other variable definitions */
1330 (void) micopy_all_var_defs(in_file->mincid, out_file->mincid,
1331 nexcluded, excluded_vars);
1332
1333 /* Add the time stamp */
1334 ncopts=0;
1335 if ((ncattinq(out_file->mincid, NC_GLOBAL, MIhistory, &datatype,
1336 &att_length) == MI_ERROR) ||
1337 (datatype != NC_CHAR))
1338 att_length = 0;
1339 att_length += strlen(tm_stamp) + 1;
1340 string = malloc(att_length);
1341 string[0] = '\0';
1342 (void) miattgetstr(out_file->mincid, NC_GLOBAL, MIhistory, att_length,
1343 string);
1344 ncopts=NC_VERBOSE | NC_FATAL;
1345 (void) strcat(string, tm_stamp);
1346 (void) miattputstr(out_file->mincid, NC_GLOBAL, MIhistory, string);
1347 free(string);
1348
1349 /* Get the dimension ids from the input file */
1350 (void) ncvarinq(in_file->mincid, in_file->imgid, NULL, NULL,
1351 &ndims, in_dims, NULL);
1352
1353 /* Check for screw-up on number of dimensions */
1354 if (ndims != out_file->ndims) {
1355 (void) fprintf(stderr,
1356 "Error in number of dimensions for output file.\n");
1357 exit(EXIT_FAILURE);
1358 }
1359
1360 /* Create the dimensions for the image variable */
1361 for (out_index=0; out_index<ndims; out_index++) {
1362
1363 /* Check to see if this is a volume dimension */
1364 is_volume_dimension = (out_file->world_axes[out_index] != NO_AXIS);
1365
1366 /* Get the input index */
1367 if (!is_volume_dimension)
1368 in_index = out_index;
1369 else {
1370 axis = out_file->world_axes[out_index];
1371 if ((axis<0) || (axis>=WORLD_NDIMS)) {
1372 (void) fprintf(stderr,
1373 "Error creating dimensions for output file.\n");
1374 exit(EXIT_FAILURE);
1375 }
1376 in_index = in_file->indices[in_file->axes[axis]];
1377 }
1378
1379 /* Get the dimension name from the input file */
1380 (void) ncdiminq(in_file->mincid, in_dims[in_index], dimname, NULL);
1381
1382 /* Check to see if the dimension already exists */
1383 ncopts = 0;
1384 out_dims[out_index] = ncdimid(out_file->mincid, dimname);
1385 ncopts = NC_VERBOSE | NC_FATAL;
1386 dim_exists = (out_dims[out_index] != MI_ERROR);
1387
1388 /* If we have a volume dimension and it exists already with the wrong
1389 size, then we must rename it */
1390 if (is_volume_dimension && dim_exists &&
1391 (out_file->nelements[out_index] != in_file->nelements[in_index])) {
1392 string = malloc(MAX_NC_NAME);
1393 ncopts = 0;
1394 idim = 0;
1395 do {
1396 (void) sprintf(string, "%s%d", dimname, idim);
1397 idim++;
1398 } while (ncdimid(out_file->mincid, string) != MI_ERROR);
1399 ncopts = NC_VERBOSE | NC_FATAL;
1400 (void) ncdimrename(out_file->mincid, out_dims[out_index], string);
1401 free(string);
1402 out_dims[out_index] = ncdimdef(out_file->mincid, dimname,
1403 out_file->nelements[out_index]);
1404 }
1405 else if (!dim_exists)
1406 out_dims[out_index] = ncdimdef(out_file->mincid, dimname,
1407 out_file->nelements[out_index]);
1408
1409 /* If this is a volume dimension, create a variable */
1410 if (is_volume_dimension) {
1411
1412 /* Create the variable */
1413 dimid = micreate_group_variable(out_file->mincid, dimname);
1414 (void) miattputdbl(out_file->mincid, dimid, MIstep,
1415 volume_def->step[axis]);
1416 (void) miattputdbl(out_file->mincid, dimid, MIstart,
1417 volume_def->start[axis]);
1418 (void) ncattput(out_file->mincid, dimid, MIdirection_cosines,
1419 NC_DOUBLE, WORLD_NDIMS, volume_def->dircos[axis]);
1420 (void) miattputstr(out_file->mincid, dimid, MIunits,
1421 volume_def->units[axis]);
1422 (void) miattputstr(out_file->mincid, dimid, MIspacetype,
1423 volume_def->spacetype[axis]);
1424
1425 } /* If volume dimension */
1426 } /* Loop over dimensions */
1427
1428 /* Create the image max and min variables. These do not vary over
1429 the volume rows and columns even if they are not image dimensions,
1430 so we have to copy down the elements of the array (excluding image
1431 dimensions). Compute number of slices (row,columns) in an image
1432 (minc file def'n) and number of images in the file. */
1433 nimage_dims = 2;
1434 ncdiminq(out_file->mincid, out_dims[ndims-1], dimname, NULL);
1435 if (strcmp(dimname, MIvector_dimension)==0)
1436 nimage_dims++;
1437 nmaxmin_dims = 0;
1438 out_file->slices_per_image = 1;
1439 out_file->images_per_file = 1;
1440 for (idim=0; idim<ndims; idim++) {
1441 if ((idim != out_file->indices[COL_AXIS]) &&
1442 (idim != out_file->indices[ROW_AXIS])) {
1443 if (idim < ndims-nimage_dims) {
1444 out_maxmin_dims[nmaxmin_dims] = out_dims[idim];
1445 nmaxmin_dims++;
1446 out_file->images_per_file *= out_file->nelements[idim];
1447 }
1448 else {
1449 out_file->slices_per_image *= out_file->nelements[idim];
1450 }
1451 }
1452 }
1453 out_file->do_slice_renormalization =
1454 ((out_file->datatype != NC_FLOAT) &&
1455 (out_file->datatype != NC_DOUBLE) &&
1456 (out_file->slices_per_image > 1));
1457
1458 /* Create the variables */
1459 out_file->maxid = micreate_std_variable(out_file->mincid, MIimagemax,
1460 NC_DOUBLE, nmaxmin_dims,
1461 out_maxmin_dims);
1462 if (in_file->maxid != MI_ERROR)
1463 (void) micopy_all_atts(in_file->mincid, in_file->maxid,
1464 out_file->mincid, out_file->maxid);
1465 out_file->minid = micreate_std_variable(out_file->mincid, MIimagemin,
1466 NC_DOUBLE, nmaxmin_dims,
1467 out_maxmin_dims);
1468 if (in_file->minid != MI_ERROR)
1469 (void) micopy_all_atts(in_file->mincid, in_file->minid,
1470 out_file->mincid, out_file->minid);
1471
1472 /* Add transformation information to image processing variable if
1473 a transformation is given on the command line */
1474
1475 if (transform_info->file_name != NULL) {
1476
1477 ncopts = 0;
1478
1479 /* Get id of processing variable (create it if needed) */
1480 varid = ncvarid(out_file->mincid, PROCESSING_VAR);
1481 if (varid == MI_ERROR) {
1482 varid = ncvardef(out_file->mincid, PROCESSING_VAR, NC_INT, 0, NULL);
1483 (void) miadd_child(out_file->mincid,
1484 ncvarid(out_file->mincid, MIrootvariable), varid);
1485 }
1486
1487 /* Look for an unused transformation attribute */
1488 string = malloc(MI_MAX_ATTSTR_LEN);
1489 itrans = 0;
1490 do {
1491 (void) sprintf(string, "transformation%d-filename", itrans);
1492 itrans++;
1493 } while (ncattinq(out_file->mincid, varid, string,
1494 NULL, NULL) != MI_ERROR);
1495 itrans--;
1496
1497 /* Reset error handling */
1498 ncopts = NC_VERBOSE | NC_FATAL;
1499
1500 /* Add the attributes describing the transformation */
1501 (void) miattputstr(out_file->mincid, varid, string,
1502 transform_info->file_name);
1503 (void) sprintf(string, "transformation%d-filedata", itrans);
1504 (void) miattputstr(out_file->mincid, varid, string,
1505 transform_info->file_contents);
1506 if (transform_info->invert_transform) {
1507 (void) sprintf(string, "transformation%d-inverted", itrans);
1508 (void) miattputstr(out_file->mincid, varid, string, MI_TRUE);
1509 }
1510
1511 free(string);
1512
1513 } /* If transform specified on command line */
1514
1515 /* Create the image variable */
1516 out_file->imgid = micreate_std_variable(out_file->mincid, MIimage,
1517 out_file->datatype,
1518 ndims, out_dims);
1519 (void) micopy_all_atts(in_file->mincid, in_file->imgid,
1520 out_file->mincid, out_file->imgid);
1521 (void) miattputstr(out_file->mincid, out_file->imgid, MIcomplete,
1522 MI_FALSE);
1523 (void) miset_valid_range(out_file->mincid, out_file->imgid,
1524 out_file->vrange);
1525 if (out_file->is_signed)
1526 (void) miattputstr(out_file->mincid, out_file->imgid,
1527 MIsigntype, MI_SIGNED);
1528 else
1529 (void) miattputstr(out_file->mincid, out_file->imgid,
1530 MIsigntype, MI_UNSIGNED);
1531
1532 /* Get into data mode */
1533 (void) ncsetfill(out_file->mincid, NC_NOFILL);
1534 (void) ncendef(out_file->mincid);
1535
1536 /* Copy all the other data */
1537 (void) micopy_all_var_values(in_file->mincid, out_file->mincid,
1538 nexcluded, excluded_vars);
1539
1540 /* Create and attach an icv */
1541 out_file->using_icv = TRUE;
1542 out_file->icvid = miicv_create();
1543 (void) miicv_setint(out_file->icvid, MI_ICV_TYPE, NC_DOUBLE);
1544 (void) miicv_setint(out_file->icvid, MI_ICV_DO_NORM, TRUE);
1545 (void) miicv_setint(out_file->icvid, MI_ICV_USER_NORM, TRUE);
1546 (void) miicv_attach(out_file->icvid, out_file->mincid, out_file->imgid);
1547
1548 return;
1549 }
1550
1551 /* ----------------------------- MNI Header -----------------------------------
1552 @NAME : get_voxel_to_world_transf
1553 @INPUT : volume_def - description of volume
1554 @OUTPUT : voxel_to_world - transformation
1555 @RETURNS : (nothing)
1556 @DESCRIPTION: Routine to convert a Volume definition specification of sampling
1557 to a voxel-to-world transformation
1558 @METHOD :
1559 @GLOBALS :
1560 @CALLS :
1561 @CREATED : February 9, 1993 (Peter Neelin)
1562 @MODIFIED :
1563 ---------------------------------------------------------------------------- */
get_voxel_to_world_transf(Volume_Definition * volume_def,General_transform * voxel_to_world)1564 static void get_voxel_to_world_transf(Volume_Definition *volume_def,
1565 General_transform *voxel_to_world)
1566 {
1567 int idim, jdim, cur_dim, vol_axis;
1568 int irregular;
1569 long ielement, dimlength;
1570 Transform matrix;
1571 Irregular_Transform_Data *irreg_transf_data;
1572 General_transform linear_transf, irreg_transf;
1573
1574 /* Make an identity matrix */
1575 make_identity_transform(&matrix);
1576
1577 /* Loop over rows of matrix */
1578 for (idim=0; idim<WORLD_NDIMS; idim++) {
1579
1580 /* Loop over columns of matrix */
1581 for (jdim=0; jdim<WORLD_NDIMS; jdim++) {
1582 cur_dim = volume_def->axes[jdim];
1583
1584 /* Get rotation/scale components of matrix */
1585 Transform_elem(matrix, idim, cur_dim) =
1586 volume_def->step[jdim] * volume_def->dircos[jdim][idim];
1587
1588 /* Get translation components */
1589 Transform_elem(matrix, idim, VOL_NDIMS) +=
1590 volume_def->start[jdim] * volume_def->dircos[jdim][idim];
1591 }
1592 }
1593
1594 /* Save the general transform */
1595 create_linear_transform(voxel_to_world, &matrix);
1596
1597 /* Check for an irregularly spaced dimension */
1598 irregular = FALSE;
1599 for (idim=0; idim < WORLD_NDIMS; idim++) {
1600 if (volume_def->coords[idim] != NULL)
1601 irregular = TRUE;
1602 }
1603
1604 /* If we have an irregularly spaced dimension, then create the appropriate
1605 transform */
1606 if (irregular) {
1607 irreg_transf_data = malloc(sizeof(*irreg_transf_data));
1608
1609 /* Loop through the axes */
1610 for (idim=0; idim < WORLD_NDIMS; idim++) {
1611
1612 vol_axis = volume_def->axes[idim];
1613 irreg_transf_data->last_index[vol_axis] = 0;
1614
1615 /* Check whether the axis is irregularly spaced or not */
1616 if (volume_def->coords[idim] == NULL) {
1617 irreg_transf_data->nelements[vol_axis] = 0;
1618 irreg_transf_data->coords[vol_axis] = NULL;
1619 }
1620 else {
1621
1622 /* If irregular then get the number of elements and
1623 allocate space */
1624 dimlength = volume_def->nelements[idim];
1625 irreg_transf_data->nelements[vol_axis] = dimlength;
1626 irreg_transf_data->coords[vol_axis] =
1627 malloc(sizeof(double) * dimlength);
1628
1629 /* Normalize the coordinates to give first coord 0 and
1630 last coord n-1 (so that we can concat with the linear
1631 transform already created */
1632 for (ielement=0; ielement < dimlength; ielement++) {
1633 irreg_transf_data->coords[vol_axis][ielement] =
1634 (volume_def->coords[idim][ielement] -
1635 volume_def->start[idim]) /
1636 volume_def->step[idim];
1637 }
1638 }
1639 }
1640
1641 /* Create an irregular transform and free the data (we only free
1642 the Irregular_Transform_Data structure, not the coord pointers,
1643 since these arrays are not copied) */
1644 create_user_transform(&irreg_transf, (void *) irreg_transf_data,
1645 sizeof(*irreg_transf_data),
1646 irregular_transform_function,
1647 irregular_inverse_transform_function);
1648 free(irreg_transf_data);
1649
1650 /* Concatenate the linear transform with the irregular transform */
1651 copy_general_transform(voxel_to_world, &linear_transf);
1652 delete_general_transform(voxel_to_world);
1653 concat_general_transforms(&irreg_transf, &linear_transf,
1654 voxel_to_world);
1655
1656 }
1657
1658 return;
1659 }
1660
1661 /* ----------------------------- MNI Header -----------------------------------
1662 @NAME : irregular_transform_function
1663 @INPUT : user_data - pointer to user data
1664 x, y, z - coordinate to transform
1665 @OUTPUT : x_trans, y_trans, z_trans - resulting coordinate
1666 @RETURNS : (nothin)
1667 @DESCRIPTION: Routine to transform irregularly spaced coordinate to a
1668 regular spacing.
1669 @METHOD :
1670 @GLOBALS :
1671 @CALLS :
1672 @CREATED : November 4, 1993 (Peter Neelin)
1673 @MODIFIED :
1674 ---------------------------------------------------------------------------- */
irregular_transform_function(void * user_data,Real x,Real y,Real z,Real * x_trans,Real * y_trans,Real * z_trans)1675 static void irregular_transform_function(void *user_data,
1676 Real x,
1677 Real y,
1678 Real z,
1679 Real *x_trans,
1680 Real *y_trans,
1681 Real *z_trans)
1682 {
1683 Irregular_Transform_Data *irreg_transf_data;
1684 int idim;
1685 long dimlength, index;
1686 double coord[VOL_NDIMS], coord_transf[VOL_NDIMS];
1687 double frac;
1688
1689 /* Get the pointer to the data */
1690 irreg_transf_data = (Irregular_Transform_Data *) user_data;
1691
1692 /* Get the coordinate */
1693 coord[0] = x; coord[1] = y; coord[2] = z;
1694
1695 /* Loop through axes, renormalizing coordinate */
1696 for (idim=0; idim < VOL_NDIMS; idim++) {
1697 dimlength = irreg_transf_data->nelements[idim];
1698 if (dimlength <= 1) {
1699 coord_transf[idim] = coord[idim];
1700 }
1701 else {
1702
1703 /* For irregular dimension, do linear interpolation between coords */
1704 index = (long) coord[idim];
1705 if (index < 0) index = 0;
1706 if (index > dimlength-2) index = dimlength-2;
1707 frac = coord[idim] - index;
1708 coord_transf[idim] =
1709 (1.0 - frac) * irreg_transf_data->coords[idim][index] +
1710 frac * irreg_transf_data->coords[idim][index + 1];
1711
1712 }
1713 }
1714
1715 /* Save the coordinates */
1716 *x_trans = coord_transf[0];
1717 *y_trans = coord_transf[1];
1718 *z_trans = coord_transf[2];
1719
1720 return;
1721 }
1722
1723 /* ----------------------------- MNI Header -----------------------------------
1724 @NAME : irregular_inverse_transform_function
1725 @INPUT : user_data - pointer to user data
1726 x, y, z - coordinate to transform
1727 @OUTPUT : x_trans, y_trans, z_trans - resulting coordinate
1728 @RETURNS : (nothin)
1729 @DESCRIPTION: Routine to transform irregularly spaced coordinate to a
1730 regular spacing.
1731 @METHOD :
1732 @GLOBALS :
1733 @CALLS :
1734 @CREATED : November 4, 1993 (Peter Neelin)
1735 @MODIFIED :
1736 ---------------------------------------------------------------------------- */
irregular_inverse_transform_function(void * user_data,Real x,Real y,Real z,Real * x_trans,Real * y_trans,Real * z_trans)1737 static void irregular_inverse_transform_function(void *user_data,
1738 Real x,
1739 Real y,
1740 Real z,
1741 Real *x_trans,
1742 Real *y_trans,
1743 Real *z_trans)
1744 {
1745 Irregular_Transform_Data *irreg_transf_data;
1746 int idim, not_found;
1747 long dimlength, index;
1748 double coord[VOL_NDIMS], coord_transf[VOL_NDIMS];
1749 double first, next, step, frac;
1750
1751 /* Get the pointer to the data */
1752 irreg_transf_data = (Irregular_Transform_Data *) user_data;
1753
1754 /* Get the coordinate */
1755 coord[0] = x; coord[1] = y; coord[2] = z;
1756
1757 /* Loop through axes, renormalizing coordinate */
1758 for (idim=0; idim < VOL_NDIMS; idim++) {
1759 dimlength = irreg_transf_data->nelements[idim];
1760 if (dimlength <= 1) {
1761 coord_transf[idim] = coord[idim];
1762 }
1763 else {
1764
1765 /* Search for the closest index (checking for out-of-range values) */
1766 index = irreg_transf_data->last_index[idim];
1767 if (index < 0) index = 0;
1768 if (index > dimlength-2) index = dimlength-2;
1769 not_found = TRUE;
1770 while (not_found) {
1771 if (coord[idim] < irreg_transf_data->coords[idim][index]) {
1772 if (index > 0) index--;
1773 else {
1774 index = 0;
1775 not_found = FALSE;
1776 }
1777 }
1778 else if (coord[idim] > irreg_transf_data->coords[idim][index+1]) {
1779 if (index < dimlength-2) index++;
1780 else {
1781 index = dimlength-2;
1782 not_found = FALSE;
1783 }
1784 }
1785 else {
1786 not_found = FALSE;
1787 }
1788
1789 }
1790 irreg_transf_data->last_index[idim] = index;
1791
1792 /* Do linear interpolation */
1793 first = irreg_transf_data->coords[idim][index];
1794 next = irreg_transf_data->coords[idim][index + 1];
1795 step = next - first;
1796 if (step == 0.0)
1797 frac = 0.0;
1798 else
1799 frac = (coord[idim] - first) / step;
1800 coord_transf[idim] =
1801 (1.0 - frac) * index + frac * (index + 1);
1802
1803 }
1804 }
1805
1806 /* Save the coordinates */
1807 *x_trans = coord_transf[0];
1808 *y_trans = coord_transf[1];
1809 *z_trans = coord_transf[2];
1810
1811 return;
1812 }
1813
1814 /* ----------------------------- MNI Header -----------------------------------
1815 @NAME : get_default_range
1816 @INPUT : what - MIvalid_min means get default min, MIvalid_min means
1817 get default min
1818 datatype - type of variable
1819 is_signed - TRUE if variable is signed
1820 @OUTPUT : (none)
1821 @RETURNS : default maximum or minimum for the datatype
1822 @DESCRIPTION: Return the defaults maximum or minimum for a given datatype
1823 and sign.
1824 @METHOD :
1825 @GLOBALS :
1826 @CALLS :
1827 @CREATED : August 10, 1992 (Peter Neelin)
1828 @MODIFIED : February 10, 1993 (Peter Neelin)
1829 - ripped off from MINC code
1830 ---------------------------------------------------------------------------- */
get_default_range(char * what,nc_type datatype,int is_signed)1831 static double get_default_range(char *what, nc_type datatype, int is_signed)
1832 {
1833 double limit;
1834
1835 if (strcmp(what, MIvalid_max)==0) {
1836 switch (datatype) {
1837 case NC_INT: limit = (is_signed) ? INT_MAX : UINT_MAX; break;
1838 case NC_SHORT: limit = (is_signed) ? SHRT_MAX : USHRT_MAX; break;
1839 case NC_BYTE: limit = (is_signed) ? SCHAR_MAX : UCHAR_MAX; break;
1840 default: limit = DEFAULT_MAX; break;
1841 }
1842 }
1843 else if (strcmp(what, MIvalid_min)==0) {
1844 switch (datatype) {
1845 case NC_INT: limit = (is_signed) ? INT_MIN : 0; break;
1846 case NC_SHORT: limit = (is_signed) ? SHRT_MIN : 0; break;
1847 case NC_BYTE: limit = (is_signed) ? SCHAR_MIN : 0; break;
1848 default: limit = DEFAULT_MIN; break;
1849 }
1850 }
1851 else {
1852 limit = 0.0;
1853 }
1854
1855 return limit;
1856 }
1857
1858
1859 /* ----------------------------- MNI Header -----------------------------------
1860 @NAME : finish_up
1861 @INPUT : in_vol - input volume
1862 out_vol - output volume
1863 @OUTPUT : (nothing)
1864 @RETURNS : (nothing)
1865 @DESCRIPTION: Routine to finish up at end of program, closing files, etc.
1866 @METHOD :
1867 @GLOBALS :
1868 @CALLS :
1869 @CREATED : February 15, 1993 (Peter Neelin)
1870 @MODIFIED :
1871 ---------------------------------------------------------------------------- */
finish_up(VVolume * in_vol,VVolume * out_vol)1872 static void finish_up(VVolume *in_vol, VVolume *out_vol)
1873 {
1874 File_Info *in_file, *out_file;
1875
1876 /* Get file info pointers */
1877 in_file = in_vol->file;
1878 out_file = out_vol->file;
1879
1880 /* Close the output file */
1881 (void) miattputstr(out_file->mincid, out_file->imgid, MIcomplete, MI_TRUE);
1882 if (out_file->using_icv) {
1883 (void) miicv_free(out_file->icvid);
1884 }
1885 (void) miclose(out_file->mincid);
1886
1887 /* Close the input file */
1888 if (in_file->using_icv) {
1889 (void) miicv_free(in_file->icvid);
1890 }
1891 (void) miclose(in_file->mincid);
1892
1893 return;
1894 }
1895
1896 /* ----------------------------- MNI Header -----------------------------------
1897 @NAME : get_transformation
1898 @INPUT : dst - Pointer to client data from argument table
1899 key - argument key
1900 nextArg - argument following key
1901 @OUTPUT : (nothing)
1902 @RETURNS : TRUE so that ParseArgv will discard nextArg, unless there
1903 is no following argument.
1904 @DESCRIPTION: Routine called by ParseArgv to read in a transformation file
1905 @METHOD :
1906 @GLOBALS :
1907 @CALLS :
1908 @CREATED : February 15, 1993 (Peter Neelin)
1909 @MODIFIED :
1910 ---------------------------------------------------------------------------- */
get_transformation(char * dst,char * key,char * nextArg)1911 static int get_transformation(char *dst, char *key, char *nextArg)
1912 /* ARGSUSED */
1913 {
1914 Transform_Info *transform_info;
1915 General_transform *transformation;
1916 FILE *fp;
1917 int ch, index;
1918
1919 /* Check for following argument */
1920 if (nextArg == NULL) {
1921 (void) fprintf(stderr,
1922 "\"%s\" option requires an additional argument\n",
1923 key);
1924 exit(EXIT_FAILURE);
1925 }
1926
1927 /* Get pointer to transform info structure */
1928 transform_info = (Transform_Info *) dst;
1929
1930 /* Save file name */
1931 transform_info->file_name = nextArg;
1932 transformation = transform_info->transformation;
1933
1934 /* Open the file */
1935 if (strcmp(nextArg, "-") == 0) {
1936 /* Create a temporary for standard input */
1937 fp=tmpfile();
1938 if (fp==NULL) {
1939 (void) fprintf(stderr, "Error opening temporary file.\n");
1940 exit(EXIT_FAILURE);
1941 }
1942 while ((ch=getc(stdin))!=EOF) (void) putc(ch, fp);
1943 rewind(fp);
1944 }
1945 else {
1946 if (open_file_with_default_suffix(nextArg,
1947 get_default_transform_file_suffix(),
1948 READ_FILE, ASCII_FORMAT, &fp) != OK) {
1949 (void) fprintf(stderr, "Error opening transformation file %s.\n",
1950 nextArg);
1951 exit(EXIT_FAILURE);
1952 }
1953 }
1954
1955 /* Read in the file for later use */
1956 if (transform_info->file_contents == NULL) {
1957 transform_info->file_contents = malloc(TRANSFORM_BUFFER_INCREMENT);
1958 transform_info->buffer_length = TRANSFORM_BUFFER_INCREMENT;
1959 }
1960 for (index = 0; (ch=getc(fp)) != EOF; index++) {
1961 if (index >= transform_info->buffer_length-1) {
1962 transform_info->buffer_length += TRANSFORM_BUFFER_INCREMENT;
1963 transform_info->file_contents =
1964 realloc(transform_info->file_contents,
1965 transform_info->buffer_length);
1966 }
1967 transform_info->file_contents[index] = (char) ch;
1968 }
1969 transform_info->file_contents[index] = '\0';
1970 rewind(fp);
1971
1972 /* Get rid of the old transformation */
1973 delete_general_transform(transformation);
1974
1975 /* Read the file */
1976 if (input_transform(fp, transform_info->file_name,
1977 transformation)!=OK) {
1978 (void) fprintf(stderr, "Error reading transformation file.\n");
1979 exit(EXIT_FAILURE);
1980 }
1981 (void) close_file(fp);
1982
1983 #ifdef TRANSFORM_CHANGE_KLUDGE
1984 Specified_transform = TRUE;
1985 #endif
1986
1987 return TRUE;
1988 }
1989
1990 /* ----------------------------- MNI Header -----------------------------------
1991 @NAME : get_model_file
1992 @INPUT : dst - Pointer to client data from argument table
1993 key - argument key
1994 nextArg - argument following key
1995 @OUTPUT : (nothing)
1996 @RETURNS : TRUE so that ParseArgv will discard nextArg unless there
1997 is no following argument.
1998 @DESCRIPTION: Routine called by ParseArgv to read in a model file (-like)
1999 @METHOD :
2000 @GLOBALS :
2001 @CALLS :
2002 @CREATED : February 15, 1993 (Peter Neelin)
2003 @MODIFIED :
2004 ---------------------------------------------------------------------------- */
get_model_file(char * dst,char * key,char * nextArg)2005 static int get_model_file(char *dst, char *key, char *nextArg)
2006 /* ARGSUSED */
2007 {
2008 Volume_Definition *volume_def;
2009 File_Info file;
2010
2011 /* Check for following argument */
2012 if (nextArg == NULL) {
2013 (void) fprintf(stderr,
2014 "\"%s\" option requires an additional argument\n",
2015 key);
2016 exit(EXIT_FAILURE);
2017 }
2018
2019 /* Get pointer to volume definition structure */
2020 volume_def = (Volume_Definition *) dst;
2021
2022 /* Get file information */
2023 get_file_info(nextArg, TRUE, volume_def, &file);
2024
2025 /* Close the file */
2026 (void) miclose(file.mincid);
2027
2028 #ifdef TRANSFORM_CHANGE_KLUDGE
2029 Specified_like = TRUE;
2030 #endif
2031
2032 return TRUE;
2033 }
2034
2035 /* ----------------------------- MNI Header -----------------------------------
2036 @NAME : set_standard_sampling
2037 @INPUT : dst - Pointer to client data from argument table
2038 key - argument key
2039 nextArg - argument following key
2040 @OUTPUT : (nothing)
2041 @RETURNS : FALSE so that ParseArgv will not discard nextArg.
2042 @DESCRIPTION: Routine called by ParseArgv to set the sampling to standard
2043 values (sets only step, start and direction cosines).
2044 @METHOD :
2045 @GLOBALS :
2046 @CALLS :
2047 @CREATED : November 14, 1995 (Peter Neelin)
2048 @MODIFIED :
2049 ---------------------------------------------------------------------------- */
set_standard_sampling(char * dst,char * key,char * nextArg)2050 static int set_standard_sampling(char *dst, char *key, char *nextArg)
2051 /* ARGSUSED */
2052 {
2053 Volume_Definition *volume_def;
2054 int idim, jdim;
2055
2056 /* Get pointer to volume definition structure */
2057 volume_def = (Volume_Definition *) dst;
2058
2059 /* Set sensible values */
2060 for (idim=0; idim < WORLD_NDIMS; idim++) {
2061 volume_def->step[idim] = 1.0;
2062 volume_def->start[idim] = 0.0;
2063 for (jdim=0; jdim < WORLD_NDIMS; jdim++) {
2064 volume_def->dircos[idim][jdim] = (idim == jdim ? 1.0 : 0.0);
2065 }
2066 }
2067
2068 return FALSE;
2069 }
2070
2071 /* ----------------------------- MNI Header -----------------------------------
2072 @NAME : set_spacetype
2073 @INPUT : dst - Pointer to client data from argument table
2074 key - argument key
2075 nextArg - argument following key
2076 @OUTPUT : (nothing)
2077 @RETURNS : TRUE if nextArg should be discarded, FALSE otherwise
2078 @DESCRIPTION: Routine called by ParseArgv to set the space type of the
2079 output sampling.
2080 @METHOD :
2081 @GLOBALS :
2082 @CALLS :
2083 @CREATED : December 12, 1995 (Peter Neelin)
2084 @MODIFIED :
2085 ---------------------------------------------------------------------------- */
set_spacetype(char * dst,char * key,char * nextArg)2086 static int set_spacetype(char *dst, char *key, char *nextArg)
2087 /* ARGSUSED */
2088 {
2089 Volume_Definition *volume_def;
2090 char *spacetype;
2091 int idim;
2092 int return_value;
2093
2094 /* Get pointer to client data */
2095 volume_def = (Volume_Definition *) dst;
2096
2097 /* Check key for spacetype */
2098 return_value = FALSE;
2099 if (strcmp(key, "-talairach") == 0) {
2100 spacetype = MI_TALAIRACH;
2101 }
2102 else {
2103
2104 /* Check for following argument */
2105 if (nextArg == NULL) {
2106 (void) fprintf(stderr,
2107 "\"%s\" option requires an additional argument\n",
2108 key);
2109 exit(EXIT_FAILURE);
2110 }
2111
2112 spacetype = nextArg;
2113
2114 return_value = TRUE;
2115 }
2116
2117 /* Copy the strings */
2118 for (idim=0; idim < WORLD_NDIMS; idim++) {
2119 (void) strncpy(volume_def->spacetype[idim], spacetype,
2120 MI_MAX_ATTSTR_LEN);
2121 volume_def->spacetype[idim][MI_MAX_ATTSTR_LEN-1] = '\0';
2122 }
2123
2124 return return_value;
2125
2126 }
2127
2128 /* ----------------------------- MNI Header -----------------------------------
2129 @NAME : set_units
2130 @INPUT : dst - Pointer to client data from argument table
2131 key - argument key
2132 nextArg - argument following key
2133 @OUTPUT : (nothing)
2134 @RETURNS : TRUE if nextArg should be discarded, FALSE otherwise
2135 @DESCRIPTION: Routine called by ParseArgv to set the units of the
2136 output sampling.
2137 @METHOD :
2138 @GLOBALS :
2139 @CALLS :
2140 @CREATED : December 12, 1995 (Peter Neelin)
2141 @MODIFIED :
2142 ---------------------------------------------------------------------------- */
set_units(char * dst,char * key,char * nextArg)2143 static int set_units(char *dst, char *key, char *nextArg)
2144 /* ARGSUSED */
2145 {
2146 Volume_Definition *volume_def;
2147 char *units;
2148 int idim;
2149
2150 /* Get pointer to client data */
2151 volume_def = (Volume_Definition *) dst;
2152
2153 /* Check for following argument */
2154 if (nextArg == NULL) {
2155 (void) fprintf(stderr,
2156 "\"%s\" option requires an additional argument\n",
2157 key);
2158 exit(EXIT_FAILURE);
2159 }
2160 units = nextArg;
2161
2162 /* Copy the strings */
2163 for (idim=0; idim < WORLD_NDIMS; idim++) {
2164 (void) strncpy(volume_def->units[idim], units, MI_MAX_ATTSTR_LEN);
2165 volume_def->units[idim][MI_MAX_ATTSTR_LEN-1] = '\0';
2166 }
2167
2168 return TRUE;
2169
2170 }
2171
2172 /* ----------------------------- MNI Header -----------------------------------
2173 @NAME : get_axis_order
2174 @INPUT : dst - Pointer to client data from argument table
2175 key - argument key
2176 nextArg - argument following key
2177 @OUTPUT : (nothing)
2178 @RETURNS : FALSE so that ParseArgv will not discard nextArg
2179 @DESCRIPTION: Routine called by ParseArgv to get the axis order
2180 @METHOD :
2181 @GLOBALS :
2182 @CALLS :
2183 @CREATED : February 15, 1993 (Peter Neelin)
2184 @MODIFIED :
2185 ---------------------------------------------------------------------------- */
get_axis_order(char * dst,char * key,char * nextArg)2186 static int get_axis_order(char *dst, char *key, char *nextArg)
2187 /* ARGSUSED */
2188 {
2189 Volume_Definition *volume_def;
2190
2191 /* Get pointer to client data */
2192 volume_def = (Volume_Definition *) dst;
2193
2194 /* Check key for order */
2195 if (strcmp(key, "-transverse") == 0) {
2196 volume_def->axes[Z] = 0;
2197 volume_def->axes[Y] = 1;
2198 volume_def->axes[X] = 2;
2199 }
2200 else if (strcmp(key, "-sagittal") == 0) {
2201 volume_def->axes[X] = 0;
2202 volume_def->axes[Z] = 1;
2203 volume_def->axes[Y] = 2;
2204 }
2205 else if (strcmp(key, "-coronal") == 0) {
2206 volume_def->axes[Y] = 0;
2207 volume_def->axes[Z] = 1;
2208 volume_def->axes[X] = 2;
2209 }
2210
2211 return FALSE;
2212 }
2213
2214 /* ----------------------------- MNI Header -----------------------------------
2215 @NAME : get_fillvalue
2216 @INPUT : dst - Pointer to client data from argument table
2217 key - argument key
2218 nextArg - argument following key
2219 @OUTPUT : (nothing)
2220 @RETURNS : FALSE so that ParseArgv will not discard nextArg
2221 @DESCRIPTION: Routine called by ParseArgv to set the fill value
2222 @METHOD :
2223 @GLOBALS :
2224 @CALLS :
2225 @CREATED : February 15, 1993 (Peter Neelin)
2226 @MODIFIED :
2227 ---------------------------------------------------------------------------- */
get_fillvalue(char * dst,char * key,char * nextArg)2228 static int get_fillvalue(char *dst, char *key, char *nextArg)
2229 /* ARGSUSED */
2230 {
2231 double *dptr;
2232
2233 /* Get pointer to client data */
2234 dptr = (double *) dst;
2235
2236 /* Check key for fill value to set */
2237 if (strcmp(key, "-fill") == 0) {
2238 *dptr = -DBL_MAX;
2239 }
2240 else if (strcmp(key, "-nofill") == 0) {
2241 *dptr = FILL_DEFAULT;
2242 }
2243
2244 return FALSE;
2245 }
2246