1 /* ----------------------------- MNI Header ----------------------------------- 2 @NAME : mincresample.h 3 @DESCRIPTION: Header file for mincresample.c 4 @METHOD : 5 @GLOBALS : 6 @CALLS : 7 @CREATED : February 8, 1993 (Peter Neelin) 8 @MODIFIED : 9 * $Log: mincresample.h,v $ 10 * Revision 6.8 2008-01-13 09:38:54 stever 11 * Avoid compiler warnings about functions and variables that are defined 12 * but not used. Remove some such functions and variables, 13 * conditionalize some, and move static declarations out of header files 14 * into C files. 15 * 16 * Revision 6.7 2005/07/13 21:34:25 bert 17 * Add sinc interpolant (ported from 1.X branch) 18 * 19 * Revision 6.6 2004/11/01 22:38:39 bert 20 * Eliminate all references to minc_def.h 21 * 22 * Revision 6.5 2004/04/27 15:31:20 bert 23 * Added -2 option 24 * 25 * Revision 6.4 2002/11/06 13:32:23 jason 26 * Fixes to mincresample: setting the interpolation type is now done 27 * through an enum rather than function pointers. 28 * 29 * Revision 6.3 2001/04/17 18:40:23 neelin 30 * Modifications to work with NetCDF 3.x 31 * In particular, changed NC_LONG to NC_INT (and corresponding longs to ints). 32 * Changed NC_UNSPECIFIED to NC_NAT. 33 * A few fixes to the configure script. 34 * 35 * Revision 6.2 2001/04/02 14:58:09 neelin 36 * Added -keep_real_range option to prevent rescaling of slices on output 37 * 38 * Revision 6.1 1999/10/19 14:45:27 neelin 39 * Fixed Log subsitutions for CVS 40 * 41 * Revision 6.0 1997/09/12 13:23:21 neelin 42 * Release of minc version 0.6 43 * 44 * Revision 5.0 1997/08/21 13:24:22 neelin 45 * Release of minc version 0.5 46 * 47 * Revision 4.1 1997/08/13 15:41:12 neelin 48 * Fixed initialization problem that caused crashing under Linux. 49 * 50 * Revision 4.0 1997/05/07 19:59:42 neelin 51 * Release of minc version 0.4 52 * 53 * Revision 3.4 1996/01/31 15:22:02 neelin 54 * Fixed bug in transformation of input sampling. 55 * 56 * Revision 3.3 1995/12/12 19:15:35 neelin 57 * Added -spacetype, -talairach and -units options. 58 * 59 * Revision 3.2 1995/11/21 14:13:20 neelin 60 * Transform input sampling with transformation and use this as default. 61 * Added -tfm_input_sampling to specify above option. 62 * Added -use_input_sampling to get old behaviour (no longer the default). 63 * Added -origin option (to specify coordinate instead of start values). 64 * Added -standard_sampling option (to set standard values of start, step 65 * and direction cosines). 66 * Added -invert_transformation option. 67 * 68 * Revision 3.1 1995/11/07 15:04:02 neelin 69 * Modified argument parsing so that only one pass is done. 70 * 71 * Revision 3.0 1995/05/15 19:30:57 neelin 72 * Release of minc version 0.3 73 * 74 * Revision 2.0 1994/09/28 10:32:48 neelin 75 * Release of minc version 0.2 76 * 77 * Revision 1.11 94/09/28 10:32:40 neelin 78 * Pre-release 79 * 80 * Revision 1.10 93/11/04 15:13:40 neelin 81 * Added support for irregularly spaced dimensions. 82 * 83 * Revision 1.9 93/11/02 11:23:56 neelin 84 * Handle imagemax/min potentially varying over slices (for vector data, etc.) 85 * 86 * Revision 1.8 93/10/20 14:05:42 neelin 87 * Added VOXEL_COORD_EPS - an epsilon for doing voxel coordinate comparisons. 88 * 89 * Revision 1.7 93/08/11 14:31:50 neelin 90 * Changed prototype for check_imageminmax. 91 * 92 * Revision 1.6 93/08/11 13:34:20 neelin 93 * Converted to use Dave MacDonald's General_transform code. 94 * Fixed bug in get_slice - for non-linear transformations coord was 95 * transformed, then used again as a starting coordinate. 96 * Handle files that have image-max/min that doesn't vary over slices. 97 * Handle files that have image-max/min varying over row/cols. 98 * Allow volume to extend to voxel edge for -nearest_neighbour interpolation. 99 * Handle out-of-range values (-fill values from a previous mincresample, for 100 * example). 101 * Save transformation file as a string attribute to processing variable. 102 * 103 @COPYRIGHT : 104 Copyright 1993 Peter Neelin, McConnell Brain Imaging Centre, 105 Montreal Neurological Institute, McGill University. 106 Permission to use, copy, modify, and distribute this 107 software and its documentation for any purpose and without 108 fee is hereby granted, provided that the above copyright 109 notice appear in all copies. The author and McGill University 110 make no representations about the suitability of this 111 software for any purpose. It is provided "as is" without 112 express or implied warranty. 113 ---------------------------------------------------------------------------- */ 114 115 /* Constants used in program */ 116 117 /* Number of dimensions for various things */ 118 #define VOL_NDIMS 3 /* Number of volume dimensions */ 119 #define WORLD_NDIMS 3 /* Number of world spatial dimensions */ 120 #define SLICE_NDIMS 2 /* Number of slice dimensions */ 121 #define MAT_NDIMS WORLD_NDIMS+1 /* Number of dims for homogenous matrix */ 122 /* For referring to world axes in arrays subscripted by WORLD_NDIMS */ 123 #define NO_AXIS -1 124 #define XAXIS 0 125 #define YAXIS 1 126 #define ZAXIS 2 127 /* For referring to volume axes in arrays subscripted by VOL_NDIMS */ 128 #define SLC_AXIS 0 129 #define ROW_AXIS 1 130 #define COL_AXIS 2 131 /* For referring to slice axes in arrays subscripted by SLICE_NDIMS */ 132 #define SLICE_ROW 0 133 #define SLICE_COL 1 134 /* For referring to world coordinates in Coord_Vector */ 135 #define XCOORD 0 136 #define YCOORD 1 137 #define ZCOORD 2 138 /* For referring to volume coordinates in Coord_Vector */ 139 #define SLICE 0 140 #define ROW 1 141 #define COLUMN 2 142 /* Various constants */ 143 #define NO_VALUE DBL_MAX /* Constant to flag fact that value not set */ 144 #define DEFAULT_MAX 1.0 145 #define DEFAULT_MIN 0.0 146 #define FILL_DEFAULT DBL_MAX /* Fillvalue indicating -nofill */ 147 #define SMALL_VALUE (100.0*FLT_MIN) /* A small floating-point value */ 148 #define VOXEL_COORD_EPS (100.0*FLT_EPSILON) /* Epsilon for voxel coords */ 149 #define TRANSFORM_BUFFER_INCREMENT 256 150 #define PROCESSING_VAR "processing" 151 #define TEMP_IMAGE_VAR "mincresample-temporary-image" 152 #ifndef TRUE 153 # define TRUE 1 154 # define FALSE 0 155 #endif 156 157 /* Types used in program */ 158 enum Interpolant_type { TRILINEAR, TRICUBIC, N_NEIGHBOUR, WINDOWED_SINC }; 159 160 typedef double Coord_Vector[WORLD_NDIMS]; 161 162 typedef struct { 163 char *name; 164 int mincid; 165 int imgid; 166 int maxid; 167 int minid; 168 int ndims; 169 nc_type datatype; 170 int is_signed; 171 double vrange[2]; /* [0]=min, [1]=max */ 172 long nelements[MAX_VAR_DIMS]; /* Size of each dimension */ 173 int world_axes[MAX_VAR_DIMS]; /* Relates variable index to X, Y, Z 174 or NO_AXIS */ 175 int indices[VOL_NDIMS]; /* Indices of volume dimenions (subscripted 176 from slowest to fastest) */ 177 int axes[WORLD_NDIMS]; /* Relates world X,Y,Z (index) to dimension 178 order (value=0,1,2; 0=slowest varying) */ 179 int using_icv; /* True if we are using an icv to read data */ 180 int icvid; /* Id of icv (if used) */ 181 long slices_per_image; /* Number of volume slices (row, column) per 182 minc file image */ 183 long images_per_file; /* Number of minc file images in the file */ 184 int do_slice_renormalization; /* Flag indicating that we need to 185 loop through the data a second time, 186 recomputing the slices to normalize 187 images properly */ 188 int keep_real_range; /* Flag indicating whether we should keep 189 the real range of the input data or not */ 190 } File_Info; 191 192 typedef struct { 193 long size[SLICE_NDIMS]; /* Size of each dimension */ 194 double *data; /* Pointer to slice data */ 195 } Slice_Data; 196 197 typedef struct Volume_Data_Struct Volume_Data; 198 typedef int (*Interpolating_Function) 199 (Volume_Data *volume, Coord_Vector coord, double *result); 200 struct Volume_Data_Struct { 201 nc_type datatype; /* Type of data in volume */ 202 int is_signed; /* Sign of data (TRUE if signed) */ 203 int use_fill; /* TRUE if fill values should be used in 204 calculation of output image max/min */ 205 double fillvalue; /* Value to return when out of bounds */ 206 double vrange[2]; /* [0]=min, [1]=max */ 207 double real_range[2]; /* Real min and max for current volume */ 208 int size[VOL_NDIMS]; /* Size of each dimension */ 209 void *data; /* Pointer to volume data */ 210 double *scale; /* Pointer to array of scales for slices */ 211 double *offset; /* Pointer to array of offsets for slices */ 212 Interpolating_Function interpolant; /* Function Pointer */ 213 }; 214 215 typedef struct { 216 File_Info *file; /* Information about associated file */ 217 Volume_Data *volume; /* Volume data for (input volume) */ 218 Slice_Data *slice; /* Slice data for (output volume) */ 219 General_transform *voxel_to_world; 220 General_transform *world_to_voxel; 221 } VVolume; 222 223 typedef struct { 224 int axes[WORLD_NDIMS]; /* Relates world X,Y,Z (index) to dimension 225 order (value=0,1,2; 0=slowest varying) */ 226 long nelements[WORLD_NDIMS]; /* These are subscripted by X, Y and Z */ 227 double step[WORLD_NDIMS]; 228 double start[WORLD_NDIMS]; 229 double dircos[WORLD_NDIMS][WORLD_NDIMS]; 230 double *coords[WORLD_NDIMS]; 231 char units[WORLD_NDIMS][MI_MAX_ATTSTR_LEN]; 232 char spacetype[WORLD_NDIMS][MI_MAX_ATTSTR_LEN]; 233 } Volume_Definition; 234 235 typedef struct { 236 int verbose; 237 } Program_Flags; 238 239 typedef struct { 240 int invert_transform; 241 char *file_name; 242 char *file_contents; 243 long buffer_length; 244 General_transform *transformation; 245 } Transform_Info; 246 247 typedef struct { 248 int clobber; 249 int keep_real_range; 250 nc_type datatype; 251 int is_signed; 252 double vrange[2]; 253 double fillvalue; 254 double origin[3]; 255 Program_Flags flags; 256 enum Interpolant_type interpolant_type; /* Type of interpolation */ 257 Transform_Info transform_info; 258 Volume_Definition volume_def; 259 int v2format; /* If non-zero, create a MINC 2.0 output */ 260 } Arg_Data; 261 262 typedef struct { 263 long last_index[VOL_NDIMS]; 264 long nelements[VOL_NDIMS]; 265 double *coords[VOL_NDIMS]; 266 } Irregular_Transform_Data; 267 268 /* Macros used in program */ 269 270 #define DO_TRANSFORM(result, transformation, coord) \ 271 general_transform_point(transformation, \ 272 coord[XCOORD], coord[YCOORD], coord[ZCOORD], \ 273 &result[XCOORD], &result[YCOORD], &result[ZCOORD]) 274 275 #define DO_INVERSE_TRANSFORM(result, transformation, coord) \ 276 general_inverse_transform_point(transformation, \ 277 coord[XCOORD], coord[YCOORD], coord[ZCOORD], \ 278 &result[XCOORD], &result[YCOORD], &result[ZCOORD]) 279 280 #define IS_LINEAR(transformation) \ 281 (get_transform_type(transformation)==LINEAR) 282 283 #define VECTOR_COPY(result, first) { \ 284 result[XCOORD] = first[XCOORD]; \ 285 result[YCOORD] = first[YCOORD]; \ 286 result[ZCOORD] = first[ZCOORD]; \ 287 } 288 289 #define VECTOR_DIFF(result, first, second) { \ 290 result[XCOORD] = first[XCOORD] - second[XCOORD]; \ 291 result[YCOORD] = first[YCOORD] - second[YCOORD]; \ 292 result[ZCOORD] = first[ZCOORD] - second[ZCOORD]; \ 293 } 294 295 #define VECTOR_ADD(result, first, second) { \ 296 result[XCOORD] = first[XCOORD] + second[XCOORD]; \ 297 result[YCOORD] = first[YCOORD] + second[YCOORD]; \ 298 result[ZCOORD] = first[ZCOORD] + second[ZCOORD]; \ 299 } 300 301 #define VECTOR_SCALAR_MULT(result, vector, scalar) { \ 302 result[XCOORD] = vector[XCOORD] * (scalar); \ 303 result[YCOORD] = vector[YCOORD] * (scalar); \ 304 result[ZCOORD] = vector[ZCOORD] * (scalar); \ 305 } 306 307 #ifdef INTERPOLATE 308 # undef INTERPOLATE 309 #endif 310 311 #define INTERPOLATE(volume, coord, result) \ 312 (*volume->interpolant) (volume, coord, result) 313 314 #define VOLUME_VALUE(volume, slcind, rowind, colind, value) \ 315 { \ 316 long offset; \ 317 \ 318 offset = ((slcind)*volume->size[ROW_AXIS] + \ 319 (rowind))*volume->size[COL_AXIS] + (colind); \ 320 switch (volume->datatype) { \ 321 case NC_BYTE: \ 322 if (volume->is_signed) \ 323 value = *((signed char *) volume->data + offset); \ 324 else \ 325 value = *((unsigned char *) volume->data + offset); \ 326 break; \ 327 case NC_SHORT: \ 328 if (volume->is_signed) \ 329 value = *((signed short *) volume->data + offset); \ 330 else \ 331 value = *((unsigned short *) volume->data + offset); \ 332 break; \ 333 case NC_INT: \ 334 if (volume->is_signed) \ 335 value = *((signed int *) volume->data + offset); \ 336 else \ 337 value = *((unsigned int *) volume->data + offset); \ 338 break; \ 339 case NC_FLOAT: \ 340 value = *((float *) volume->data + offset); \ 341 break; \ 342 case NC_DOUBLE: \ 343 value = *((double *) volume->data + offset); \ 344 break; \ 345 } \ 346 } 347 348 /* Function prototypes */ 349 350 extern void resample_volumes(Program_Flags *program_flags, 351 VVolume *in_vol, VVolume *out_vol, 352 General_transform *transformation); 353 extern int trilinear_interpolant(Volume_Data *volume, 354 Coord_Vector coord, double *result); 355 extern int tricubic_interpolant(Volume_Data *volume, 356 Coord_Vector coord, double *result); 357 extern int nearest_neighbour_interpolant(Volume_Data *volume, 358 Coord_Vector coord, double *result); 359 extern int windowed_sinc_interpolant(Volume_Data *volume, 360 Coord_Vector coord, double *result); 361 362 #define SINC_HALF_WIDTH_MAX 10 363 #define SINC_HALF_WIDTH_MIN 1 364 365 enum sinc_interpolant_window_t { 366 SINC_WINDOW_NONE, 367 SINC_WINDOW_HANNING, 368 SINC_WINDOW_HAMMING 369 }; 370 371 extern enum sinc_interpolant_window_t sinc_window_type; 372 extern int sinc_half_width; 373