1 /* ----------------------------- MNI Header -----------------------------------
2 @NAME : resample_volumes.c
3 @DESCRIPTION: File containing routines to resample minc volumes.
4 @METHOD :
5 @GLOBALS :
6 @CREATED : February 8, 1993 (Peter Neelin)
7 @MODIFIED :
8 * $Log: resample_volumes.c,v $
9 * Revision 6.12 2008-01-17 02:33:02 rotor
10 * * removed all rcsids
11 * * removed a bunch of ^L's that somehow crept in
12 * * removed old (and outdated) BUGS file
13 *
14 * Revision 6.11 2008/01/13 09:49:00 stever
15 * Add do_Ncubic_interpolation declaration, moved here from mincresample.h.
16 *
17 * Revision 6.10 2008/01/13 09:38:54 stever
18 * Avoid compiler warnings about functions and variables that are defined
19 * but not used. Remove some such functions and variables,
20 * conditionalize some, and move static declarations out of header files
21 * into C files.
22 *
23 * Revision 6.9 2008/01/12 19:08:15 stever
24 * Add __attribute__ ((unused)) to all rcsid variables.
25 *
26 * Revision 6.8 2008/01/11 07:17:08 stever
27 * Remove unused variables.
28 *
29 * Revision 6.7 2007/12/12 20:55:26 rotor
30 * * added a bunch of bug fixes from Claude.
31 *
32 * Revision 6.6 2006/05/19 00:13:39 bert
33 * Add config.h
34 *
35 * Revision 6.5 2005/07/13 21:34:25 bert
36 * Add sinc interpolant (ported from 1.X branch)
37 *
38 * Revision 6.4 2004/11/01 22:38:39 bert
39 * Eliminate all references to minc_def.h
40 *
41 * Revision 6.3 2001/08/16 13:32:39 neelin
42 * Partial fix for valid_range of different type from image (problems
43 * arising from double to float conversion/rounding). NOT COMPLETE.
44 *
45 * Revision 6.2 2001/04/02 14:58:10 neelin
46 * Added -keep_real_range option to prevent rescaling of slices on output
47 *
48 * Revision 6.1 1999/10/19 14:45:28 neelin
49 * Fixed Log subsitutions for CVS
50 *
51 * Revision 6.0 1997/09/12 13:23:21 neelin
52 * Release of minc version 0.6
53 *
54 * Revision 5.0 1997/08/21 13:24:22 neelin
55 * Release of minc version 0.5
56 *
57 * Revision 4.0 1997/05/07 19:59:42 neelin
58 * Release of minc version 0.4
59 *
60 * Revision 3.0 1995/05/15 19:30:57 neelin
61 * Release of minc version 0.3
62 *
63 * Revision 2.0 1994/09/28 10:32:50 neelin
64 * Release of minc version 0.2
65 *
66 * Revision 1.15 94/09/28 10:32:38 neelin
67 * Pre-release
68 *
69 * Revision 1.14 94/03/15 10:58:47 neelin
70 * Fixed tricubic interpolation (wasn't initializing variable for fillvalue
71 * detection).
72 *
73 * Revision 1.13 94/03/01 14:33:48 neelin
74 * Fixed error in calculating valid minimum for float volumes.
75 *
76 * Revision 1.12 93/11/02 11:23:52 neelin
77 * Handle imagemax/min potentially varying over slices (for vector data, etc.)
78 *
79 * Revision 1.11 93/10/20 14:06:36 neelin
80 * Modified tri-linear interpolation to allow volumes to extend epsilon
81 * beyond the first voxel.
82 * Added code to handle volume dimensions of size one (in tri-linear
83 * interpolation).
84 *
85 * Revision 1.10 93/10/15 13:48:22 neelin
86 * Removed include of recipes.h
87 *
88 * Revision 1.9 93/10/12 12:48:08 neelin
89 * Use volume_io.h instead of def_mni.h
90 *
91 * Revision 1.8 93/08/11 14:29:53 neelin
92 * Use volume->datatype in load_volume instead of file->datatype.
93 * Loop for slice max/min is from 0, not 1 in load_volume.
94 *
95 * Revision 1.7 93/08/11 13:34:15 neelin
96 * Converted to use Dave MacDonald's General_transform code.
97 * Fixed bug in get_slice - for non-linear transformations coord was
98 * transformed, then used again as a starting coordinate.
99 * Handle files that have image-max/min that doesn't vary over slices.
100 * Handle files that have image-max/min varying over row/cols.
101 * Allow volume to extend to voxel edge for -nearest_neighbour interpolation.
102 * Handle out-of-range values (-fill values from a previous mincresample, for
103 * example).
104 * Save transformation file as a string attribute to processing variable.
105 *
106 @COPYRIGHT :
107 Copyright 1993 Peter Neelin, McConnell Brain Imaging Centre,
108 Montreal Neurological Institute, McGill University.
109 Permission to use, copy, modify, and distribute this
110 software and its documentation for any purpose and without
111 fee is hereby granted, provided that the above copyright
112 notice appear in all copies. The author and McGill University
113 make no representations about the suitability of this
114 software for any purpose. It is provided "as is" without
115 express or implied warranty.
116 ---------------------------------------------------------------------------- */
117
118 #if HAVE_CONFIG_H
119 #include "config.h"
120 #endif
121
122 #include <stdlib.h>
123 #include <stdio.h>
124 #include <float.h>
125 #include <string.h>
126 #include <math.h>
127 #include <minc.h>
128 #include <volume_io.h>
129 #include "mincresample.h"
130
131 static void load_volume(File_Info *file, long start[], long count[],
132 Volume_Data *volume);
133 static void get_slice(long slice_num, VVolume *in_vol, VVolume *out_vol,
134 General_transform *transformation,
135 double *minimum, double *maximum);
136 static void renormalize_slices(Program_Flags *program_flags, VVolume *out_vol,
137 double slice_min[], double slice_max[]);
138 static int do_Ncubic_interpolation(Volume_Data *volume,
139 long index[], int cur_dim,
140 double frac[], double *result);
141
142
143
144 /* ----------------------------- MNI Header -----------------------------------
145 @NAME : resample_volumes
146 @INPUT : program_flags - data for program execution
147 in_vol - description of input volume
148 out_vol - description of output volume
149 transformation - description of world transformation
150 @OUTPUT : (none)
151 @RETURNS : (none)
152 @DESCRIPTION: Resamples in_vol into file specified by out_vol using given
153 world transformation.
154 @METHOD :
155 @GLOBALS :
156 @CALLS :
157 @CREATED : February 8, 1993 (Peter Neelin)
158 @MODIFIED :
159 ---------------------------------------------------------------------------- */
resample_volumes(Program_Flags * program_flags,VVolume * in_vol,VVolume * out_vol,General_transform * transformation)160 void resample_volumes(Program_Flags *program_flags,
161 VVolume *in_vol, VVolume *out_vol,
162 General_transform *transformation)
163 {
164 long in_start[MAX_VAR_DIMS], in_count[MAX_VAR_DIMS], in_end[MAX_VAR_DIMS];
165 long out_start[MAX_VAR_DIMS], out_count[MAX_VAR_DIMS];
166 long mm_start[MAX_VAR_DIMS]; /* Vector for min/max variables */
167 long nslice, islice, slice_count;
168 int idim, index, slice_index;
169 double maximum, minimum, valid_range[2];
170 double *slice_max, *slice_min;
171 File_Info *ifp,*ofp;
172
173 /* Set pointers to file information */
174 ifp = in_vol->file;
175 ofp = out_vol->file;
176
177 /* Allocate slice min/max arrays if needed */
178 if (ofp->do_slice_renormalization) {
179 slice_min = malloc(ofp->images_per_file * ofp->slices_per_image *
180 sizeof(double));
181 slice_max = malloc(ofp->images_per_file * ofp->slices_per_image *
182 sizeof(double));
183 }
184
185 /* Set input file start, count and end vectors for reading a volume
186 at a time */
187 (void) miset_coords(ifp->ndims, (long) 0, in_start);
188 (void) miset_coords(ifp->ndims, (long) 1, in_count);
189 for (idim=0; idim < ifp->ndims; idim++) {
190 in_end[idim] = ifp->nelements[idim];
191 }
192 for (idim=0; idim < VOL_NDIMS; idim++) {
193 index = ifp->indices[idim];
194 in_count[index] = ifp->nelements[index];
195 }
196
197 /* Set output file count for writing a slice and get the number of
198 output slices */
199 (void) miset_coords(ifp->ndims, (long) 1, out_count);
200 for (idim=0; idim < VOL_NDIMS; idim++) {
201 index = ofp->indices[idim];
202 if (idim==0) {
203 slice_index = index;
204 nslice = ofp->nelements[index];
205 }
206 else {
207 out_count[index] = ofp->nelements[index];
208 }
209 }
210
211 /* Initialize global max and min */
212 valid_range[0] = DBL_MAX;
213 valid_range[1] = -DBL_MAX;
214
215 /* Initialize file max/min slice count */
216 slice_count = 0;
217
218 /* Print log message */
219 if (program_flags->verbose) {
220 (void) fprintf(stderr, "Transforming slices:");
221 (void) fflush(stderr);
222 }
223
224 /* Loop over input volumes */
225
226 while (in_start[0] < in_end[0]) {
227
228 /* Copy the start vector */
229 for (idim=0; idim < ifp->ndims; idim++)
230 out_start[idim] = in_start[idim];
231
232 /* Read in the volume */
233 load_volume(ifp, in_start, in_count, in_vol->volume);
234
235 /* Loop over slices */
236 for (islice=0; islice < nslice; islice++) {
237
238 /* Print log message */
239 if (program_flags->verbose) {
240 (void) fprintf(stderr, ".");
241 (void) fflush(stderr);
242 }
243
244 /* Set slice number in out_start */
245 out_start[slice_index] = islice;
246
247 /* Get the slice */
248 get_slice(islice, in_vol, out_vol, transformation,
249 &minimum, &maximum);
250
251 /* Check whether we are keep the input range */
252 if (ofp->keep_real_range) {
253 minimum = in_vol->volume->real_range[0];
254 maximum = in_vol->volume->real_range[1];
255 }
256
257 /* Update global max and min */
258 if (maximum > valid_range[1]) valid_range[1] = maximum;
259 if (minimum < valid_range[0]) valid_range[0] = minimum;
260
261 /* Write the max, min and slice */
262 (void) mivarput1(ofp->mincid, ofp->maxid,
263 mitranslate_coords(ofp->mincid,
264 ofp->imgid, out_start,
265 ofp->maxid, mm_start),
266 NC_DOUBLE, NULL, &maximum);
267 (void) mivarput1(ofp->mincid, ofp->minid,
268 mitranslate_coords(ofp->mincid,
269 ofp->imgid, out_start,
270 ofp->minid, mm_start),
271 NC_DOUBLE, NULL, &minimum);
272 (void) miicv_put(ofp->icvid, out_start, out_count,
273 out_vol->slice->data);
274
275 /* Save the max, min if needed */
276 if (ofp->do_slice_renormalization) {
277 slice_max[slice_count] = maximum;
278 slice_min[slice_count] = minimum;
279 }
280
281 /* Increment slice count */
282 slice_count++;
283
284 } /* End loop over slices */
285
286 /* Increment in_start counter */
287 idim = ofp->ndims-1;
288 in_start[idim] += in_count[idim];
289 while ( (idim>0) && (in_start[idim] >= in_end[idim])) {
290 in_start[idim] = 0;
291 idim--;
292 in_start[idim] += in_count[idim];
293 }
294
295 } /* End loop over volumes */
296
297 /* Print end of log message */
298 if (program_flags->verbose) {
299 (void) fprintf(stderr, "Done\n");
300 (void) fflush(stderr);
301 }
302
303 /* If output volume is floating point, write out global max and min */
304 if ((ofp->datatype == NC_FLOAT) || (ofp->datatype == NC_DOUBLE)) {
305 (void) miset_valid_range(ofp->mincid, ofp->imgid, valid_range);
306 }
307
308 /* Recompute slices and free vectors, if needed */
309 if (ofp->do_slice_renormalization) {
310 renormalize_slices(program_flags, out_vol, slice_min, slice_max);
311 free(slice_min);
312 free(slice_max);
313 }
314
315 }
316
317 /* ----------------------------- MNI Header -----------------------------------
318 @NAME : load_volume
319 @INPUT : file - description of input file
320 start - index of start of volume in minc file
321 count - vector size of volume in minc file
322 volume - description of volume data
323 @OUTPUT : volume - contains new volume data and scales and offsets
324 @RETURNS : (none)
325 @DESCRIPTION: Loads a volume from a minc file.
326 @METHOD :
327 @GLOBALS :
328 @CALLS :
329 @CREATED : February 10, 1993 (Peter Neelin)
330 @MODIFIED :
331 ---------------------------------------------------------------------------- */
load_volume(File_Info * file,long start[],long count[],Volume_Data * volume)332 void load_volume(File_Info *file, long start[], long count[],
333 Volume_Data *volume)
334 {
335 long nread, islice, mm_start[MAX_VAR_DIMS], mm_count[MAX_VAR_DIMS];
336 int varid, ivar, idim, ndims;
337 double *values, maximum, minimum, denom;
338
339 /* Load the file */
340 if (file->using_icv) {
341 (void) miicv_get(file->icvid, start, count, volume->data);
342 }
343 else {
344 (void) ncvarget(file->mincid, file->imgid,
345 start, count, volume->data);
346 }
347
348 /* Read the max and min from the file into the scale and offset variables
349 (maxima into scale and minima into offset) if datatype is not
350 floating point */
351
352 /* Don't do it ourselves if the icv is doing it for us */
353 if (file->using_icv) {
354 for (islice=0; islice < volume->size[SLC_AXIS]; islice++) {
355 (void) miicv_inqdbl(file->icvid, MI_ICV_NORM_MAX,
356 &volume->scale[islice]);
357 (void) miicv_inqdbl(file->icvid, MI_ICV_NORM_MIN,
358 &volume->offset[islice]);
359 }
360 }
361
362 /* If either max/min variable doesn't exist, or if type is floating
363 point, then set max to 1 and min to 0 */
364 else if ((file->maxid == MI_ERROR) || (file->minid == MI_ERROR) ||
365 (volume->datatype == NC_FLOAT) ||
366 (volume->datatype == NC_DOUBLE)) {
367 for (islice=0; islice < volume->size[SLC_AXIS]; islice++) {
368 volume->scale[islice] = DEFAULT_MAX;
369 volume->offset[islice] = DEFAULT_MIN;
370 }
371 }
372
373 /* Otherwise the imagemax/min variables exist - use them for integer
374 types */
375 else {
376
377 for (ivar=0; ivar<2; ivar++) {
378
379 /* Set up variables */
380 varid = (ivar == 0 ? file->maxid : file->minid);
381 values = (ivar == 0 ? volume->scale : volume->offset);
382
383 /* Read max or min */
384 (void) mivarget(file->mincid, varid,
385 mitranslate_coords(file->mincid, file->imgid, start,
386 varid,
387 miset_coords(MAX_VAR_DIMS, 0L, mm_start)),
388 mitranslate_coords(file->mincid, file->imgid, count,
389 varid,
390 miset_coords(MAX_VAR_DIMS, 1L, mm_count)),
391 NC_DOUBLE, NULL, values);
392
393 /* Check for number of values read */
394 (void) ncvarinq(file->mincid, varid, NULL, NULL, &ndims, NULL, NULL);
395 for (nread=1, idim=0; idim<ndims; idim++) nread *= mm_count[idim];
396
397 /* If only one value read, then copy it for each slice */
398 if (nread==1) {
399 for (islice=1; islice < volume->size[SLC_AXIS]; islice++) {
400 values[islice] = values[0];
401 }
402 }
403 else if (nread != volume->size[SLC_AXIS]) {
404 (void) fprintf(stderr,
405 "Program bug while reading image max/min\n");
406 exit(EXIT_FAILURE);
407 }
408
409 } /* Loop over max/min variables */
410
411 } /* If max/min variables both exist */
412
413 /* Get the real range of this volume (offset is min, scale is max) */
414 volume->real_range[0] = volume->offset[0];
415 volume->real_range[1] = volume->scale[0];
416 for (islice=1; islice < volume->size[SLC_AXIS]; islice++) {
417 if (volume->offset[islice] < volume->real_range[0])
418 volume->real_range[0] = volume->offset[islice];
419 if (volume->scale[islice] > volume->real_range[1])
420 volume->real_range[1] = volume->scale[islice];
421 }
422
423 /* Calculate the scale and offset */
424 for (islice=0; islice < volume->size[SLC_AXIS]; islice++) {
425
426 /* If the variables type is floating point, then don't scale */
427 if ((volume->datatype==NC_FLOAT) || (volume->datatype==NC_DOUBLE)) {
428 volume->scale[islice] = 1.0;
429 volume->offset[islice] = 0.0;
430 }
431
432 /* Otherwise, calculate scale and offset */
433 else {
434 maximum = volume->scale[islice];
435 minimum = volume->offset[islice];
436 denom = file->vrange[1] - file->vrange[0];
437 if (denom != 0.0) {
438 volume->scale[islice] = (maximum - minimum) / denom;
439 }
440 else {
441 volume->scale[islice] = 0.0;
442 }
443 volume->offset[islice] = minimum -
444 file->vrange[0] * volume->scale[islice];
445 }
446
447 } /* End of loop through slices */
448 }
449
450 /* ----------------------------- MNI Header -----------------------------------
451 @NAME : get_slice
452 @INPUT : in_vol - description of input volume
453 out_vol - description of output volume
454 transformation - description of world transformation
455 @OUTPUT : out_vol - slice field contains new slice
456 minimum - slice minimum (excluding data from outside volume)
457 maximum - slice maximum (excluding data from outside volume)
458 @RETURNS : (none)
459 @DESCRIPTION: Resamples current volume of in_vol into slice in out_vol
460 using given world transformation.
461 @METHOD :
462 @GLOBALS :
463 @CALLS :
464 @CREATED : February 8, 1993 (Peter Neelin)
465 @MODIFIED :
466 ---------------------------------------------------------------------------- */
get_slice(long slice_num,VVolume * in_vol,VVolume * out_vol,General_transform * transformation,double * minimum,double * maximum)467 void get_slice(long slice_num, VVolume *in_vol, VVolume *out_vol,
468 General_transform *transformation,
469 double *minimum, double *maximum)
470 {
471 Slice_Data *slice;
472 Volume_Data *volume;
473 double *dptr;
474 long irow, icol;
475 int all_linear;
476 int idim;
477
478 /* Coordinate vectors for stepping through slice */
479 Coord_Vector zero = {0, 0, 0};
480 Coord_Vector column = {0, 0, 1};
481 Coord_Vector row = {0, 1, 0};
482 Coord_Vector start = {0, 0, 0}; /* start[SLICE] set later to slice_num */
483 Coord_Vector coord, transf_coord;
484
485 /* Transformation stuff */
486 General_transform total_transf, temp_transf;
487
488 /* Get slice and volume pointers */
489 volume = in_vol->volume;
490 slice = out_vol->slice;
491
492 /* Concatenate transforms */
493 concat_general_transforms(out_vol->voxel_to_world,
494 transformation, &temp_transf);
495 concat_general_transforms(&temp_transf, in_vol->world_to_voxel,
496 &total_transf);
497 delete_general_transform(&temp_transf);
498
499 /* Check for complete linear transformation */
500 all_linear = (get_transform_type(&total_transf) == LINEAR);
501
502 /* Transform vectors for linear transformation */
503 start[SLICE] = slice_num;
504 if (all_linear) {
505 DO_TRANSFORM(zero, &total_transf, zero);
506 DO_TRANSFORM(column, &total_transf, column);
507 DO_TRANSFORM(row, &total_transf, row);
508 DO_TRANSFORM(start, &total_transf, start);
509 }
510
511 /* Make sure that row and column are vectors and not points */
512 VECTOR_DIFF(row, row, zero);
513 VECTOR_DIFF(column, column, zero);
514
515 /* Initialize maximum and minimum */
516 *maximum = -DBL_MAX;
517 *minimum = DBL_MAX;
518
519 /* Loop over rows of slice */
520
521 for (irow=0; irow < slice->size[SLICE_ROW]; irow++) {
522
523 /* Set starting coordinate of row */
524 VECTOR_SCALAR_MULT(coord, row, irow);
525 VECTOR_ADD(coord, coord, start);
526
527 /* Loop over columns */
528
529 dptr = slice->data + irow*slice->size[SLICE_COL];
530 for (icol=0; icol < slice->size[SLICE_COL]; icol++) {
531
532 /* If transformation is not completely linear, then transform
533 voxel to world, world to world and world to voxel, as needed */
534 for (idim=0; idim<WORLD_NDIMS; idim++)
535 transf_coord[idim]=coord[idim];
536 if (!all_linear) {
537 DO_TRANSFORM(transf_coord, &total_transf, transf_coord);
538 }
539
540 /* Do interpolation */
541 if (INTERPOLATE(volume, transf_coord, dptr) || volume->use_fill) {
542 if (*dptr > *maximum) *maximum = *dptr;
543 if (*dptr < *minimum) *minimum = *dptr;
544 }
545
546 /* Increment coordinate */
547 VECTOR_ADD(coord, coord, column);
548
549 /* Increment slice pointer */
550 dptr++;
551
552 } /* Loop over columns */
553 } /* Loop over rows */
554
555 if ((*maximum == -DBL_MAX) && (*minimum == DBL_MAX)) {
556 *minimum = 0.0;
557 *maximum = SMALL_VALUE;
558 }
559 else if (*maximum <= *minimum) {
560 if (*minimum == 0.0)
561 *maximum = SMALL_VALUE;
562 else if (*minimum < 0.0)
563 *maximum = 0.0;
564 else
565 *maximum = 2.0 * (*minimum);
566 }
567
568 /* Delete the transformation */
569 delete_general_transform(&total_transf);
570
571 }
572
573 /* ----------------------------- MNI Header -----------------------------------
574 @NAME : trilinear_interpolant
575 @INPUT : volume - pointer to volume data
576 coord - point at which volume should be interpolated in voxel
577 units (with 0 being first point of the volume).
578 @OUTPUT : result - interpolated value.
579 @RETURNS : TRUE if coord is within the volume, FALSE otherwise.
580 @DESCRIPTION: Routine to interpolate a volume at a point with tri-linear
581 interpolation.
582 @METHOD :
583 @GLOBALS :
584 @CALLS :
585 @CREATED : February 10, 1993 (Peter Neelin)
586 @MODIFIED :
587 ---------------------------------------------------------------------------- */
trilinear_interpolant(Volume_Data * volume,Coord_Vector coord,double * result)588 int trilinear_interpolant(Volume_Data *volume,
589 Coord_Vector coord, double *result)
590 {
591 long slcind, rowind, colind, slcmax, rowmax, colmax;
592 long slcnext, rownext, colnext;
593 static double f0, f1, f2, r0, r1, r2, r1r2, r1f2, f1r2, f1f2;
594 static double v000, v001, v010, v011, v100, v101, v110, v111;
595
596 /* Check that the coordinate is inside the volume */
597 slcmax = volume->size[SLC_AXIS] - 1;
598 rowmax = volume->size[ROW_AXIS] - 1;
599 colmax = volume->size[COL_AXIS] - 1;
600 if ((coord[SLICE] < -VOXEL_COORD_EPS) ||
601 (coord[SLICE] > slcmax+VOXEL_COORD_EPS) ||
602 (coord[ROW] < -VOXEL_COORD_EPS) ||
603 (coord[ROW] > rowmax+VOXEL_COORD_EPS) ||
604 (coord[COLUMN] < -VOXEL_COORD_EPS) ||
605 (coord[COLUMN] > colmax+VOXEL_COORD_EPS)) {
606 *result = volume->fillvalue;
607 return FALSE;
608 }
609
610 /* Get the whole part of the coordinate */
611 slcind = (long) coord[SLICE];
612 rowind = (long) coord[ROW];
613 colind = (long) coord[COLUMN];
614 if (slcind >= slcmax-1) slcind = slcmax-1;
615 if (rowind >= rowmax-1) rowind = rowmax-1;
616 if (colind >= colmax-1) colind = colmax-1;
617
618 /* Get the next voxel up */
619 slcnext = slcind+1;
620 rownext = rowind+1;
621 colnext = colind+1;
622
623 /* Check for case of dimension of length one */
624 if (slcmax == 0) {
625 slcind = 0;
626 slcnext = 0;
627 }
628 if (rowmax == 0) {
629 rowind = 0;
630 rownext = 0;
631 }
632 if (colmax == 0) {
633 colind = 0;
634 colnext = 0;
635 }
636
637 /* Get the relevant voxels */
638 VOLUME_VALUE(volume, slcind , rowind , colind , v000);
639 VOLUME_VALUE(volume, slcind , rowind , colnext, v001);
640 VOLUME_VALUE(volume, slcind , rownext, colind , v010);
641 VOLUME_VALUE(volume, slcind , rownext, colnext, v011);
642 VOLUME_VALUE(volume, slcnext, rowind , colind , v100);
643 VOLUME_VALUE(volume, slcnext, rowind , colnext, v101);
644 VOLUME_VALUE(volume, slcnext, rownext, colind , v110);
645 VOLUME_VALUE(volume, slcnext, rownext, colnext, v111);
646
647 /* Check that the values are not fill values */
648 if ((v000 < volume->vrange[0]) || (v000 > volume->vrange[1]) ||
649 (v001 < volume->vrange[0]) || (v001 > volume->vrange[1]) ||
650 (v010 < volume->vrange[0]) || (v010 > volume->vrange[1]) ||
651 (v011 < volume->vrange[0]) || (v011 > volume->vrange[1]) ||
652 (v100 < volume->vrange[0]) || (v100 > volume->vrange[1]) ||
653 (v101 < volume->vrange[0]) || (v101 > volume->vrange[1]) ||
654 (v110 < volume->vrange[0]) || (v110 > volume->vrange[1]) ||
655 (v111 < volume->vrange[0]) || (v111 > volume->vrange[1])) {
656 *result = volume->fillvalue;
657 return FALSE;
658 }
659
660 /* Get the fraction parts */
661 f0 = coord[SLICE] - slcind;
662 f1 = coord[ROW] - rowind;
663 f2 = coord[COLUMN] - colind;
664 r0 = 1.0 - f0;
665 r1 = 1.0 - f1;
666 r2 = 1.0 - f2;
667
668 /* Do the interpolation */
669 r1r2 = r1 * r2;
670 r1f2 = r1 * f2;
671 f1r2 = f1 * r2;
672 f1f2 = f1 * f2;
673 *result =
674 r0 * (volume->scale[slcind] *
675 (r1r2 * v000 +
676 r1f2 * v001 +
677 f1r2 * v010 +
678 f1f2 * v011) + volume->offset[slcind]);
679 *result +=
680 f0 * (volume->scale[slcind+1] *
681 (r1r2 * v100 +
682 r1f2 * v101 +
683 f1r2 * v110 +
684 f1f2 * v111) + volume->offset[slcind+1]);
685
686 return TRUE;
687
688 }
689
690 /* ----------------------------- MNI Header -----------------------------------
691 @NAME : tricubic_interpolant
692 @INPUT : volume - pointer to volume data
693 coord - point at which volume should be interpolated in voxel
694 units (with 0 being first point of the volume).
695 @OUTPUT : result - interpolated value.
696 @RETURNS : TRUE if coord is within the volume, FALSE otherwise.
697 @DESCRIPTION: Routine to interpolate a volume at a point with tri-cubic
698 interpolation.
699 @METHOD :
700 @GLOBALS :
701 @CALLS :
702 @CREATED : February 12, 1993 (Peter Neelin)
703 @MODIFIED :
704 ---------------------------------------------------------------------------- */
tricubic_interpolant(Volume_Data * volume,Coord_Vector coord,double * result)705 int tricubic_interpolant(Volume_Data *volume,
706 Coord_Vector coord, double *result)
707 {
708 long slcind, rowind, colind, slcmax, rowmax, colmax, index[VOL_NDIMS];
709 double frac[VOL_NDIMS];
710
711 /* Check that the coordinate is inside the volume */
712 slcmax = volume->size[SLC_AXIS] - 1;
713 rowmax = volume->size[ROW_AXIS] - 1;
714 colmax = volume->size[COL_AXIS] - 1;
715
716 /* Identify a slice or a volume. Slice must be in x-y. */
717 if( slcmax == 0 ) {
718 /* 2-d slice */
719 if( (coord[ROW] < 0) || (coord[ROW] > rowmax) ||
720 (coord[COLUMN] < 0) || (coord[COLUMN] > colmax) ) {
721 *result = volume->fillvalue;
722 return FALSE;
723 }
724 } else {
725 /* 3-d volume */
726 if( (coord[SLICE] < 0) || (coord[SLICE] > slcmax) ||
727 (coord[ROW] < 0) || (coord[ROW] > rowmax) ||
728 (coord[COLUMN] < 0) || (coord[COLUMN] > colmax) ) {
729 *result = volume->fillvalue;
730 return FALSE;
731 }
732 }
733
734 /* Get the whole and fractional part of the coordinate */
735 slcind = (long)coord[SLICE];
736 rowind = (long)coord[ROW];
737 colind = (long)coord[COLUMN];
738 frac[0] = coord[SLICE] - slcind;
739 frac[1] = coord[ROW] - rowind;
740 frac[2] = coord[COLUMN] - colind;
741 slcind--;
742 rowind--;
743 colind--;
744
745 /* Check for dimension of image. */
746
747 /* Check for edges in the 2-d plane - do linear interpolation at edges */
748 /* Note: Spline stencil is right-sided, so ok to start at 0. */
749 if( slcmax == 0 && rowmax > 0 && colmax > 0 ) {
750 if( (rowind > rowmax-3) || (rowind < 0) ||
751 (colind > colmax-3) || (colind < 0)) {
752 return trilinear_interpolant(volume, coord, result);
753 } else {
754 slcind = 0; /* there is only slice 0 */
755 index[0]=slcind; index[1]=rowind; index[2]=colind;
756 if( do_Ncubic_interpolation(volume, index, 1, frac, result) ) {
757 /* scaling not done for a slice in do_Ncubic_interpolation, */
758 /* only done for a volume */
759 *result = (*result) * volume->scale[slcind] + volume->offset[slcind];
760 return TRUE;
761 } else {
762 return FALSE;
763 }
764 }
765 }
766
767 /* Check for edges in the 3-d volume - do linear interpolation at edges */
768 /* Note: Spline stencil is right-sided, so ok to start at 0. */
769 if ((slcind > slcmax-3) || (slcind < 0) ||
770 (rowind > rowmax-3) || (rowind < 0) ||
771 (colind > colmax-3) || (colind < 0)) {
772 return trilinear_interpolant(volume, coord, result);
773 } else {
774 index[0]=slcind; index[1]=rowind; index[2]=colind;
775 /* Do the interpolation and return its value */
776 return do_Ncubic_interpolation(volume, index, 0, frac, result);
777 }
778
779 }
780
781 /* ----------------------------- MNI Header -----------------------------------
782 @NAME : do_Ncubic_interpolation
783 @INPUT : volume - pointer to volume data
784 index - indices to start point in volume
785 (bottom of 4x4x4 cube for interpolation)
786 cur_dim - dimension to be interpolated (0 = volume, 1 = slice,
787 2 = line)
788 @OUTPUT : result - interpolated value.
789 @RETURNS : TRUE if coord is within the volume, FALSE otherwise.
790 @DESCRIPTION: Routine to interpolate a volume, slice or line (specified by
791 cur_dim).
792 @METHOD :
793 @GLOBALS :
794 @CALLS :
795 @CREATED : February 12, 1993 (Peter Neelin)
796 @MODIFIED :
797 ---------------------------------------------------------------------------- */
do_Ncubic_interpolation(Volume_Data * volume,long index[],int cur_dim,double frac[],double * result)798 int do_Ncubic_interpolation(Volume_Data *volume,
799 long index[], int cur_dim,
800 double frac[], double *result)
801 {
802 long base_index;
803 double v0, v1, v2, v3, u;
804 int found_fillvalue;
805
806 /* Save index that we will change */
807 base_index = index[cur_dim];
808
809 /* If last dimension, then just get the values */
810 found_fillvalue = FALSE;
811 if (cur_dim == VOL_NDIMS-1) {
812 VOLUME_VALUE(volume, index[0] ,index[1], index[2], v0);
813 index[cur_dim]++;
814 VOLUME_VALUE(volume, index[0] ,index[1], index[2], v1);
815 index[cur_dim]++;
816 VOLUME_VALUE(volume, index[0] ,index[1], index[2], v2);
817 index[cur_dim]++;
818 VOLUME_VALUE(volume, index[0] ,index[1], index[2], v3);
819
820 /* Check for fillvalues */
821 if ((v0 < volume->vrange[0]) || (v0 > volume->vrange[1]) ||
822 (v1 < volume->vrange[0]) || (v1 > volume->vrange[1]) ||
823 (v2 < volume->vrange[0]) || (v2 > volume->vrange[1]) ||
824 (v3 < volume->vrange[0]) || (v3 > volume->vrange[1])) {
825 found_fillvalue = TRUE;
826 }
827
828 }
829
830 /* Otherwise, recurse */
831 else {
832 if (!do_Ncubic_interpolation(volume, index, cur_dim+1, frac, &v0)) {
833 found_fillvalue = TRUE;
834 }
835 index[cur_dim]++;
836 if (!do_Ncubic_interpolation(volume, index, cur_dim+1, frac, &v1)) {
837 found_fillvalue = TRUE;
838 }
839 index[cur_dim]++;
840 if (!do_Ncubic_interpolation(volume, index, cur_dim+1, frac, &v2)) {
841 found_fillvalue = TRUE;
842 }
843 index[cur_dim]++;
844 if (!do_Ncubic_interpolation(volume, index, cur_dim+1, frac, &v3)) {
845 found_fillvalue = TRUE;
846 }
847 }
848
849 /* Restore index */
850 index[cur_dim] = base_index;
851
852 /* Check for fill value found */
853 if (found_fillvalue) {
854 *result = volume->fillvalue;
855 return FALSE;
856 }
857
858 /* Scale values for slices */
859 if (cur_dim==0) {
860 v0 = v0 * volume->scale[base_index ] + volume->offset[base_index ];
861 v1 = v1 * volume->scale[base_index+1] + volume->offset[base_index+1];
862 v2 = v2 * volume->scale[base_index+2] + volume->offset[base_index+2];
863 v3 = v3 * volume->scale[base_index+3] + volume->offset[base_index+3];
864 }
865
866 /* Get fraction */
867 u = frac[cur_dim];
868
869 /* Do tricubic interpolation (code from Dave MacDonald).
870 Gives v1 and v2 at u = 0 and 1 and gives continuity of intensity
871 and first derivative. */
872 *result =
873 ( (v1) + (u) * (
874 0.5 * ((v2)-(v0)) + (u) * (
875 (v0) - 2.5 * (v1) + 2.0 * (v2) - 0.5 * (v3) + (u) * (
876 -0.5 * (v0) + 1.5 * (v1) - 1.5 * (v2) + 0.5 * (v3) )
877 )
878 )
879 );
880
881 return TRUE;
882 }
883
884 /************************************************************************
885 * Windowed Sinc Interpolant
886 *
887 * The technique used is borrowed from Neil Thacker et al., "Improved
888 * Quality of Re-sliced MR Images Using Re-normalized Sinc
889 * Interpolation", Journal of Magnetic Resonance Imaging 10:582-588
890 * (1999).
891 *
892 * Any bugs are of course my own fault!
893 *
894 * -bert
895 */
896
897 int sinc_half_width = SINC_HALF_WIDTH_MAX / 2;
898
899 enum sinc_interpolant_window_t sinc_window_type = SINC_WINDOW_HANNING;
900
901 /* basic windowed sinc function */
902
903 static double
windowed_sinc(double delta)904 windowed_sinc(double delta)
905 {
906 double phase;
907 double sinc;
908 double window;
909
910 /* Calculate the sinc function.
911 */
912 phase = delta * M_PI;
913
914 if (phase == 0.0) {
915 sinc = 1.0;
916 }
917 else {
918 sinc = sin(phase) / phase;
919 }
920
921 switch (sinc_window_type) {
922 case SINC_WINDOW_HANNING:
923 /* Calculate the Hanning window.
924 */
925 window = 0.50 + 0.50 * cos(phase / (1.0 + sinc_half_width));
926 break;
927
928 case SINC_WINDOW_HAMMING:
929 /* Calculate the Hamming window.
930 */
931 window = 0.54 + 0.46 * cos(phase / (1.0 + sinc_half_width));
932 break;
933
934 default:
935 window = 1.0; /* No window */
936 break;
937 }
938 return (sinc * window);
939 }
940
941 /* Floating-point (unscaled) multiply/accumulate operations */
942 #define SINC_FRND result += *pix_ptr++ * *win_ptr++
943
944 /* Unroll those loops!! */
945 #define SINC_FMAC \
946 result = 0.0; \
947 switch (sinc_half_width) { \
948 case 10: SINC_FRND; SINC_FRND; \
949 case 9: SINC_FRND; SINC_FRND; \
950 case 8: SINC_FRND; SINC_FRND; \
951 case 7: SINC_FRND; SINC_FRND; \
952 case 6: SINC_FRND; SINC_FRND; \
953 case 5: SINC_FRND; SINC_FRND; \
954 case 4: SINC_FRND; SINC_FRND; \
955 case 3: SINC_FRND; SINC_FRND; \
956 case 2: SINC_FRND; SINC_FRND; \
957 case 1: SINC_FRND; SINC_FRND; \
958 SINC_FRND; /* Do the leftover */ \
959 }
960
961 /* Integer (scaled) multiply/accumulate operations */
962 #define SINC_IRND result += ((slope * *pix_ptr++) + intercept) * *win_ptr++
963
964 /* Unroll the loops!! */
965 #define SINC_IMAC \
966 result = 0.0; \
967 switch (sinc_half_width) { \
968 case 10: SINC_IRND; SINC_IRND; \
969 case 9: SINC_IRND; SINC_IRND; \
970 case 8: SINC_IRND; SINC_IRND; \
971 case 7: SINC_IRND; SINC_IRND; \
972 case 6: SINC_IRND; SINC_IRND; \
973 case 5: SINC_IRND; SINC_IRND; \
974 case 4: SINC_IRND; SINC_IRND; \
975 case 3: SINC_IRND; SINC_IRND; \
976 case 2: SINC_IRND; SINC_IRND; \
977 case 1: SINC_IRND; SINC_IRND; \
978 SINC_IRND; /* Do the leftover */ \
979 }
980
981 static double
sinc_mac_d(Volume_Data * volume,int z,int y,int x,double * win_ptr)982 sinc_mac_d(Volume_Data *volume, int z, int y, int x, double *win_ptr)
983 {
984 double result;
985 long offset;
986 double *pix_ptr;
987
988 offset = (z * volume->size[ROW_AXIS] + y) * volume->size[COL_AXIS] + x;
989
990 pix_ptr = (double *) volume->data + offset;
991
992 SINC_FMAC;
993
994 return (result);
995 }
996
997 static double
sinc_mac_f(Volume_Data * volume,int z,int y,int x,double * win_ptr)998 sinc_mac_f(Volume_Data *volume, int z, int y, int x, double *win_ptr)
999 {
1000 double result;
1001 long offset;
1002 float *pix_ptr;
1003
1004 offset = (z * volume->size[ROW_AXIS] + y) * volume->size[COL_AXIS] + x;
1005
1006 pix_ptr = (float *) volume->data + offset;
1007
1008 SINC_FMAC;
1009
1010 return (result);
1011 }
1012
1013 static double
sinc_mac_uc(Volume_Data * volume,int z,int y,int x,double * win_ptr)1014 sinc_mac_uc(Volume_Data *volume, int z, int y, int x, double *win_ptr)
1015 {
1016 double result;
1017 double slope = volume->scale[z];
1018 double intercept = volume->offset[z];
1019 long offset;
1020 unsigned char *pix_ptr;
1021
1022 offset = (z * volume->size[ROW_AXIS] + y) * volume->size[COL_AXIS] + x;
1023
1024 pix_ptr = (unsigned char *) volume->data + offset;
1025
1026 SINC_IMAC;
1027
1028 return (result);
1029 }
1030
1031 static double
sinc_mac_sc(Volume_Data * volume,int z,int y,int x,double * win_ptr)1032 sinc_mac_sc(Volume_Data *volume, int z, int y, int x, double *win_ptr)
1033 {
1034 double result;
1035 double slope = volume->scale[z];
1036 double intercept = volume->offset[z];
1037 long offset;
1038 char *pix_ptr;
1039
1040 offset = (z * volume->size[ROW_AXIS] + y) * volume->size[COL_AXIS] + x;
1041
1042 pix_ptr = (char *) volume->data + offset;
1043
1044 SINC_IMAC;
1045
1046 return (result);
1047 }
1048
1049 static double
sinc_mac_us(Volume_Data * volume,int z,int y,int x,double * win_ptr)1050 sinc_mac_us(Volume_Data *volume, int z, int y, int x, double *win_ptr)
1051 {
1052 double result;
1053 double slope = volume->scale[z];
1054 double intercept = volume->offset[z];
1055 long offset;
1056 unsigned short *pix_ptr;
1057
1058 offset = (z * volume->size[ROW_AXIS] + y) * volume->size[COL_AXIS] + x;
1059
1060 pix_ptr = (unsigned short *) volume->data + offset;
1061
1062 SINC_IMAC;
1063
1064 return (result);
1065 }
1066
1067 static double
sinc_mac_ss(Volume_Data * volume,int z,int y,int x,double * win_ptr)1068 sinc_mac_ss(Volume_Data *volume, int z, int y, int x, double *win_ptr)
1069 {
1070 double result;
1071 double slope = volume->scale[z];
1072 double intercept = volume->offset[z];
1073 long offset;
1074 short *pix_ptr;
1075
1076 offset = (z * volume->size[ROW_AXIS] + y) * volume->size[COL_AXIS] + x;
1077
1078 pix_ptr = (short *) volume->data + offset;
1079
1080 SINC_IMAC;
1081
1082 return (result);
1083 }
1084
1085 static double
sinc_mac_ui(Volume_Data * volume,int z,int y,int x,double * win_ptr)1086 sinc_mac_ui(Volume_Data *volume, int z, int y, int x, double *win_ptr)
1087 {
1088 double result;
1089 double slope = volume->scale[z];
1090 double intercept = volume->offset[z];
1091 long offset;
1092 unsigned int *pix_ptr;
1093
1094 offset = (z * volume->size[ROW_AXIS] + y) * volume->size[COL_AXIS] + x;
1095
1096 pix_ptr = (unsigned int *) volume->data + offset;
1097
1098 SINC_IMAC;
1099
1100 return (result);
1101 }
1102
1103 static double
sinc_mac_si(Volume_Data * volume,int z,int y,int x,double * win_ptr)1104 sinc_mac_si(Volume_Data *volume, int z, int y, int x, double *win_ptr)
1105 {
1106 double result;
1107 double slope = volume->scale[z];
1108 double intercept = volume->offset[z];
1109 long offset;
1110 int *pix_ptr;
1111
1112 offset = (z * volume->size[ROW_AXIS] + y) * volume->size[COL_AXIS] + x;
1113
1114 pix_ptr = (int *) volume->data + offset;
1115
1116 SINC_IMAC;
1117
1118 return (result);
1119 }
1120
1121 /* ----------------------------- MNI Header -----------------------------------
1122 @NAME : windowed_sinc_interpolant
1123 @INPUT : volume - pointer to volume data
1124 coord - point at which volume should be interpolated in voxel
1125 units (with 0 being first point of the volume).
1126 @OUTPUT : result - interpolated value.
1127 @RETURNS : TRUE if coord is within the volume, FALSE otherwise.
1128 @DESCRIPTION: Routine to interpolate a volume at a point with windowed
1129 sinc interpolation.
1130 @METHOD :
1131 @GLOBALS :
1132 @CALLS :
1133 @CREATED : July 11 2005 (Robert Vincent)
1134 @MODIFIED :
1135 ---------------------------------------------------------------------------- */
1136 int
windowed_sinc_interpolant(Volume_Data * volume,Coord_Vector coord,double * result)1137 windowed_sinc_interpolant(Volume_Data *volume, Coord_Vector coord,
1138 double *result)
1139 {
1140 double zt, yt, xt;
1141 double zf, yf, xf;
1142 int zi, yi, xi;
1143 int i, j;
1144 double new;
1145 double tmp;
1146 long slcmax, rowmax, colmax;
1147 double zw[SINC_HALF_WIDTH_MAX * 2 + 1];
1148 double yw[SINC_HALF_WIDTH_MAX * 2 + 1];
1149 double xw[SINC_HALF_WIDTH_MAX * 2 + 1];
1150
1151 slcmax = volume->size[SLC_AXIS] - 1;
1152 rowmax = volume->size[ROW_AXIS] - 1;
1153 colmax = volume->size[COL_AXIS] - 1;
1154
1155 if ((coord[SLICE] < 0) || (coord[SLICE] > slcmax) ||
1156 (coord[ROW] < 0) || (coord[ROW] > rowmax) ||
1157 (coord[COLUMN] < 0) || (coord[COLUMN] > colmax)) {
1158 *result = volume->fillvalue;
1159 return FALSE;
1160 }
1161
1162 zi = (int) coord[SLICE];
1163 yi = (int) coord[ROW];
1164 xi = (int) coord[COLUMN];
1165
1166 /* Check for edges - do linear interpolation at edges */
1167 if ((zi > slcmax-sinc_half_width) || (zi < sinc_half_width) ||
1168 (yi > rowmax-sinc_half_width) || (yi < sinc_half_width) ||
1169 (xi > colmax-sinc_half_width) || (xi < sinc_half_width)) {
1170 return trilinear_interpolant(volume, coord, result);
1171 }
1172
1173 /* Calculate fractional part of the coordinate.
1174 */
1175 zf = coord[SLICE] - zi;
1176 yf = coord[ROW] - yi;
1177 xf = coord[COLUMN] - xi;
1178
1179 /* Initialize the totals.
1180 */
1181 zt = 0.0;
1182 yt = 0.0;
1183 xt = 0.0;
1184
1185 /* Generate the three windowed sinc functions.
1186 */
1187 for (i = -sinc_half_width; i <= sinc_half_width; i++) {
1188 tmp = windowed_sinc(zf - i);
1189 zw[i + sinc_half_width] = tmp;
1190 zt += tmp;
1191
1192 tmp = windowed_sinc(yf - i);
1193 yw[i + sinc_half_width] = tmp;
1194 yt += tmp;
1195
1196 tmp = windowed_sinc(xf - i);
1197 xw[i + sinc_half_width] = tmp;
1198 xt += tmp;
1199 }
1200
1201 /* Now calculate the new value.
1202 */
1203 new = 0.0;
1204 for (i = -sinc_half_width; i <= sinc_half_width; i++) {
1205 for (j = -sinc_half_width; j <= sinc_half_width; j++) {
1206 switch (volume->datatype) {
1207 case NC_BYTE:
1208 if (volume->is_signed) {
1209 tmp = sinc_mac_sc(volume,
1210 zi + i,
1211 yi + j,
1212 xi - sinc_half_width,
1213 xw);
1214 }
1215 else {
1216 tmp = sinc_mac_uc(volume,
1217 zi + i,
1218 yi + j,
1219 xi - sinc_half_width,
1220 xw);
1221 }
1222 break;
1223 case NC_SHORT:
1224 if (volume->is_signed) {
1225 tmp = sinc_mac_ss(volume,
1226 zi + i,
1227 yi + j,
1228 xi - sinc_half_width,
1229 xw);
1230 }
1231 else {
1232 tmp = sinc_mac_us(volume,
1233 zi + i,
1234 yi + j,
1235 xi - sinc_half_width,
1236 xw);
1237 }
1238 break;
1239 case NC_INT:
1240 if (volume->is_signed) {
1241 tmp = sinc_mac_si(volume,
1242 zi + i,
1243 yi + j,
1244 xi - sinc_half_width,
1245 xw);
1246 }
1247 else {
1248 tmp = sinc_mac_ui(volume,
1249 zi + i,
1250 yi + j,
1251 xi - sinc_half_width,
1252 xw);
1253 }
1254 break;
1255 case NC_FLOAT:
1256 tmp = sinc_mac_f(volume,
1257 zi + i,
1258 yi + j,
1259 xi - sinc_half_width,
1260 xw);
1261 break;
1262 case NC_DOUBLE:
1263 tmp = sinc_mac_d(volume,
1264 zi + i,
1265 yi + j,
1266 xi - sinc_half_width,
1267 xw);
1268 break;
1269
1270 default:
1271 fprintf(stderr, "UNHANDLED TYPE!!!\n");
1272 break;
1273 }
1274 new += zw[i + sinc_half_width] * yw[j + sinc_half_width] * tmp;
1275 }
1276 }
1277 *result = (new / (zt * yt * xt));
1278 return TRUE;
1279 }
1280
1281
1282 /* ----------------------------- MNI Header -----------------------------------
1283 @NAME : nearest_neighbour_interpolant
1284 @INPUT : volume - pointer to volume data
1285 coord - point at which volume should be interpolated in voxel
1286 units (with 0 being first point of the volume).
1287 @OUTPUT : result - interpolated value.
1288 @RETURNS : TRUE if coord is within the volume, FALSE otherwise.
1289 @DESCRIPTION: Routine to interpolate a volume at a point with nearest
1290 neighbour interpolation. Allows the coord to be outside
1291 the volume by up to 1/2 a pixel.
1292 @METHOD :
1293 @GLOBALS :
1294 @CALLS :
1295 @CREATED : February 12, 1993 (Peter Neelin)
1296 @MODIFIED :
1297 ---------------------------------------------------------------------------- */
nearest_neighbour_interpolant(Volume_Data * volume,Coord_Vector coord,double * result)1298 int nearest_neighbour_interpolant(Volume_Data *volume,
1299 Coord_Vector coord, double *result)
1300 {
1301 long slcind, rowind, colind, slcmax, rowmax, colmax;
1302
1303 /* Check that the coordinate is inside the volume */
1304 slcmax = volume->size[SLC_AXIS] - 1;
1305 rowmax = volume->size[ROW_AXIS] - 1;
1306 colmax = volume->size[COL_AXIS] - 1;
1307 slcind = ROUND(coord[SLICE]);
1308 rowind = ROUND(coord[ROW]);
1309 colind = ROUND(coord[COLUMN]);
1310 if ((slcind < 0) || (slcind > slcmax) ||
1311 (rowind < 0) || (rowind > rowmax) ||
1312 (colind < 0) || (colind > colmax)) {
1313 *result = volume->fillvalue;
1314 return FALSE;
1315 }
1316
1317 /* Get the value */
1318 VOLUME_VALUE(volume, slcind , rowind , colind , *result);
1319
1320 /* Check for fillvalue on input */
1321 if ((*result < volume->vrange[0]) || (*result > volume->vrange[1])) {
1322 *result = volume->fillvalue;
1323 return FALSE;
1324 }
1325
1326 *result = volume->scale[slcind] * (*result) + volume->offset[slcind];
1327
1328 return TRUE;
1329
1330 }
1331
1332 /* ----------------------------- MNI Header -----------------------------------
1333 @NAME : renormalize_slices
1334 @INPUT : ofp - output file pointer
1335 slice_min - array of slice minima
1336 slice_max - array of slice maxima
1337 @OUTPUT : (none)
1338 @RETURNS : (nothing)
1339 @DESCRIPTION: Routine to loop through the output file and renormalize the
1340 slices.
1341 @METHOD :
1342 @GLOBALS :
1343 @CALLS :
1344 @CREATED : October 29, 1993 (Peter Neelin)
1345 @MODIFIED :
1346 ---------------------------------------------------------------------------- */
renormalize_slices(Program_Flags * program_flags,VVolume * out_vol,double slice_min[],double slice_max[])1347 static void renormalize_slices(Program_Flags *program_flags, VVolume *out_vol,
1348 double slice_min[], double slice_max[])
1349 {
1350 File_Info *ofp;
1351 long start[MAX_VAR_DIMS], count[MAX_VAR_DIMS], end[MAX_VAR_DIMS];
1352 long mm_start[MAX_VAR_DIMS];
1353 long nslice, image, islice, ivolume, image_slice, slice_count;
1354 int idim, slice_index, index;
1355 double *image_maximum, *image_minimum;
1356
1357 /* Set pointer to file information */
1358 ofp = out_vol->file;
1359
1360 /* Set output file start, count and end vectors for stepping
1361 a slice at a time */
1362 (void) miset_coords(ofp->ndims, (long) 0, start);
1363 (void) miset_coords(ofp->ndims, (long) 1, count);
1364 for (idim=0; idim < ofp->ndims; idim++) {
1365 end[idim] = ofp->nelements[idim];
1366 }
1367 for (idim=0; idim < VOL_NDIMS; idim++) {
1368 index = ofp->indices[idim];
1369 if (idim==0) {
1370 slice_index = index;
1371 nslice = ofp->nelements[index];
1372 count[slice_index] = 1;
1373 }
1374 else {
1375 count[index] = ofp->nelements[index];
1376 }
1377 }
1378
1379 /* Find the max/min for each image */
1380 image_maximum = malloc( sizeof(double) * ofp->images_per_file);
1381 image_minimum = malloc( sizeof(double) * ofp->images_per_file);
1382 for (image=0; image < ofp->images_per_file; image++) {
1383 image_maximum[image] = -DBL_MAX;
1384 image_minimum[image] = DBL_MAX;
1385 }
1386 slice_count = 0;
1387 for (ivolume=0; ivolume < ofp->images_per_file / nslice; ivolume++) {
1388 for (image_slice=0; image_slice<ofp->slices_per_image; image_slice++) {
1389 for (islice=0; islice < nslice; islice++) {
1390 image = ivolume * nslice + islice;
1391 image_maximum[image] =
1392 MAX(image_maximum[image], slice_max[slice_count]);
1393 image_minimum[image] =
1394 MIN(image_minimum[image], slice_min[slice_count]);
1395 slice_count++;
1396 }
1397 }
1398 }
1399
1400 /* Initialize file max/min slice count */
1401 slice_count = 0;
1402
1403 /* Print log message */
1404 if (program_flags->verbose) {
1405 (void) fprintf(stderr, "Renormalizing slices:");
1406 (void) fflush(stderr);
1407 }
1408
1409 /* Loop over output volumes */
1410
1411 while (start[0] < end[0]) {
1412
1413 /* Loop over slices */
1414 for (islice=0; islice < nslice; islice++) {
1415
1416 /* Print log message */
1417 if (program_flags->verbose) {
1418 (void) fprintf(stderr, ".");
1419 (void) fflush(stderr);
1420 }
1421
1422 /* Set slice number in out_start */
1423 start[slice_index] = islice;
1424
1425 /* Figure out which image, slice, volume, etc we are dealing with */
1426 islice = slice_count % nslice;
1427 ivolume = slice_count / nslice / ofp->slices_per_image;
1428 image = ivolume * nslice + islice;
1429
1430 /* Get the slice min/max start coordinate */
1431 (void) mitranslate_coords(ofp->mincid,
1432 ofp->imgid, start,
1433 ofp->maxid, mm_start),
1434
1435 /* Write out the slice max and min */
1436 (void) mivarput1(ofp->mincid, ofp->maxid, mm_start,
1437 NC_DOUBLE, NULL, &slice_max[slice_count]);
1438 (void) mivarput1(ofp->mincid, ofp->minid, mm_start,
1439 NC_DOUBLE, NULL, &slice_min[slice_count]);
1440
1441 /* Read in the slice */
1442 (void) miicv_get(ofp->icvid, start, count, out_vol->slice->data);
1443
1444 /* Write the image max, min and slice */
1445 (void) mivarput1(ofp->mincid, ofp->maxid, mm_start,
1446 NC_DOUBLE, NULL, &image_maximum[image]);
1447 (void) mivarput1(ofp->mincid, ofp->minid, mm_start,
1448 NC_DOUBLE, NULL, &image_minimum[image]);
1449 (void) miicv_put(ofp->icvid, start, count, out_vol->slice->data);
1450
1451 /* Increment slice count */
1452 slice_count++;
1453
1454 } /* End loop over slices */
1455
1456 /* Set slice counter to ensure that we move to the next volume */
1457 start[slice_index] = end[slice_index] - 1;
1458
1459 /* Increment in_start counter */
1460 idim = ofp->ndims-1;
1461 start[idim] += count[idim];
1462 while ( (idim>0) && (start[idim] >= end[idim])) {
1463 start[idim] = 0;
1464 idim--;
1465 start[idim] += count[idim];
1466 }
1467
1468 } /* End loop over volumes */
1469
1470 /* Print end of log message */
1471 if (program_flags->verbose) {
1472 (void) fprintf(stderr, "Done\n");
1473 (void) fflush(stderr);
1474 }
1475
1476 /* Free the image max/min arrays */
1477 free(image_maximum);
1478 free(image_minimum);
1479
1480 return;
1481 }
1482
1483