1 /* ----------------------------------------------------------------------------
2 @COPYRIGHT :
3 Copyright 1993,1994,1995 David MacDonald,
4 McConnell Brain Imaging Centre,
5 Montreal Neurological Institute, McGill University.
6 Permission to use, copy, modify, and distribute this
7 software and its documentation for any purpose and without
8 fee is hereby granted, provided that the above copyright
9 notice appear in all copies. The author and McGill University
10 make no representations about the suitability of this
11 software for any purpose. It is provided "as is" without
12 express or implied warranty.
13 ---------------------------------------------------------------------------- */
14
15 #include <internal_volume_io.h>
16
17 /* ----------------------------- MNI Header -----------------------------------
18 @NAME : convert_voxel_to_value
19 @INPUT : volume
20 voxel
21 @OUTPUT :
22 @RETURNS : real value
23 @DESCRIPTION: Converts a voxel value to a real value.
24 @METHOD :
25 @GLOBALS :
26 @CALLS :
27 @CREATED : Sep. 1, 1995 David MacDonald
28 @MODIFIED :
29 ---------------------------------------------------------------------------- */
30
convert_voxel_to_value(Volume volume,Real voxel)31 VIOAPI Real convert_voxel_to_value(
32 Volume volume,
33 Real voxel )
34 {
35 if( volume->real_range_set )
36 return( volume->real_value_scale * voxel +
37 volume->real_value_translation );
38 else
39 return( voxel );
40 }
41
42 /* ----------------------------- MNI Header -----------------------------------
43 @NAME : convert_value_to_voxel
44 @INPUT : volume
45 value
46 @OUTPUT :
47 @RETURNS : voxel value
48 @DESCRIPTION: Converts a real value to a voxel value.
49 @METHOD :
50 @GLOBALS :
51 @CALLS :
52 @CREATED : Sep. 1, 1995 David MacDonald
53 @MODIFIED :
54 ---------------------------------------------------------------------------- */
55
convert_value_to_voxel(Volume volume,Real value)56 VIOAPI Real convert_value_to_voxel(
57 Volume volume,
58 Real value )
59 {
60 if( volume->real_range_set )
61 return( (value - volume->real_value_translation) /
62 volume->real_value_scale );
63 else
64 return( value );
65 }
66
67 /* ----------------------------- MNI Header -----------------------------------
68 @NAME : get_volume_voxel_value
69 @INPUT : volume
70 v0 voxel indices
71 v1
72 v2
73 v3
74 v4
75 @OUTPUT :
76 @RETURNS : Voxel value
77 @DESCRIPTION: Returns the voxel at the specified voxel index.
78 @METHOD :
79 @GLOBALS :
80 @CALLS :
81 @CREATED : Mar. 20, 1995 David MacDonald
82 @MODIFIED :
83 ---------------------------------------------------------------------------- */
84
get_volume_voxel_value(Volume volume,int v0,int v1,int v2,int v3,int v4)85 VIOAPI Real get_volume_voxel_value(
86 Volume volume,
87 int v0,
88 int v1,
89 int v2,
90 int v3,
91 int v4 )
92 {
93 Real voxel;
94
95 GET_VOXEL_TYPED( voxel, (Real), volume, v0, v1, v2, v3, v4 );
96
97 return( voxel );
98 }
99
100 /* ----------------------------- MNI Header -----------------------------------
101 @NAME : get_volume_real_value
102 @INPUT : volume
103 v0 voxel indices
104 v1
105 v2
106 v3
107 v4
108 @OUTPUT :
109 @RETURNS : Real value
110 @DESCRIPTION: Returns the volume real value at the specified voxel index.
111 @METHOD :
112 @GLOBALS :
113 @CALLS :
114 @CREATED : Mar. 20, 1995 David MacDonald
115 @MODIFIED :
116 ---------------------------------------------------------------------------- */
117
get_volume_real_value(Volume volume,int v0,int v1,int v2,int v3,int v4)118 VIOAPI Real get_volume_real_value(
119 Volume volume,
120 int v0,
121 int v1,
122 int v2,
123 int v3,
124 int v4 )
125 {
126 Real voxel, value;
127
128 voxel = get_volume_voxel_value( volume, v0, v1, v2, v3, v4 );
129
130 value = convert_voxel_to_value( volume, voxel );
131
132 return( value );
133 }
134
135 /* ----------------------------- MNI Header -----------------------------------
136 @NAME : set_volume_voxel_value
137 @INPUT : volume
138 v0 voxel indices
139 v1
140 v2
141 v3
142 v4
143 voxel
144 @OUTPUT :
145 @RETURNS :
146 @DESCRIPTION: Sets the voxel at the specified voxel index.
147 @METHOD :
148 @GLOBALS :
149 @CALLS :
150 @CREATED : Mar. 20, 1995 David MacDonald
151 @MODIFIED :
152 ---------------------------------------------------------------------------- */
153
set_volume_voxel_value(Volume volume,int v0,int v1,int v2,int v3,int v4,Real voxel)154 VIOAPI void set_volume_voxel_value(
155 Volume volume,
156 int v0,
157 int v1,
158 int v2,
159 int v3,
160 int v4,
161 Real voxel )
162 {
163 SET_VOXEL( volume, v0, v1, v2, v3, v4, voxel );
164 }
165
166 /* ----------------------------- MNI Header -----------------------------------
167 @NAME : set_volume_real_value
168 @INPUT : volume
169 v0 voxel indices
170 v1
171 v2
172 v3
173 v4
174 value
175 @OUTPUT :
176 @RETURNS :
177 @DESCRIPTION: Sets the volume real value at the specified voxel index.
178 @METHOD :
179 @GLOBALS :
180 @CALLS :
181 @CREATED : Mar. 20, 1995 David MacDonald
182 @MODIFIED :
183 ---------------------------------------------------------------------------- */
184
set_volume_real_value(Volume volume,int v0,int v1,int v2,int v3,int v4,Real value)185 VIOAPI void set_volume_real_value(
186 Volume volume,
187 int v0,
188 int v1,
189 int v2,
190 int v3,
191 int v4,
192 Real value )
193 {
194 Real voxel;
195 Data_types data_type;
196
197 voxel = convert_value_to_voxel( volume, value );
198
199 data_type = get_volume_data_type( volume );
200
201 if( data_type != FLOAT &&
202 data_type != DOUBLE )
203 {
204 voxel = (Real) ROUND( voxel );
205 }
206
207 set_volume_voxel_value( volume, v0, v1, v2, v3, v4, voxel );
208 }
209
210 /* ----------------------------- MNI Header -----------------------------------
211 @NAME : trilinear_interpolate
212 @INPUT : volume
213 voxel
214 outside_value
215 @OUTPUT : value
216 derivs
217 @RETURNS :
218 @DESCRIPTION: Computes the trilinear interpolation of the 8 coeficients, passing
219 the value back in the 'value' parameter. If the derivs parameter
220 is not null, then the 3 derivatives are also computed and
221 passed back.
222 @METHOD :
223 @GLOBALS :
224 @CALLS :
225 @CREATED : Mar. 20, 1995 David MacDonald
226 @MODIFIED : Jan. 8, 1996 D. MacDonald - now check for outside volume is
227 done inside this procedure
228 ---------------------------------------------------------------------------- */
229
trilinear_interpolate(Volume volume,Real voxel[],Real outside_value,Real * value,Real derivs[])230 static void trilinear_interpolate(
231 Volume volume,
232 Real voxel[],
233 Real outside_value,
234 Real *value,
235 Real derivs[] )
236 {
237 int c, i, j, k, *sizes, dx, dy, dz;
238 Real x, y, z, u, v, w;
239 Real coefs[8];
240 Real du00, du01, du10, du11, c00, c01, c10, c11, c0, c1, du0, du1;
241 Real dv0, dv1, dw, scale_factor;
242
243 sizes = &volume->array.sizes[0];
244
245 x = voxel[0];
246 y = voxel[1];
247 z = voxel[2];
248
249 if( x >= 0.0 && x < (Real) sizes[0]-1.0 &&
250 y >= 0.0 && y < (Real) sizes[1]-1.0 &&
251 z >= 0.0 && z < (Real) sizes[2]-1.0 )
252 {
253 i = (int) x;
254 j = (int) y;
255 k = (int) z;
256
257 GET_VOXEL_3D_TYPED( coefs[0], (Real), volume, i , j , k );
258 GET_VOXEL_3D_TYPED( coefs[1], (Real), volume, i , j , k+1 );
259 GET_VOXEL_3D_TYPED( coefs[2], (Real), volume, i , j+1, k );
260 GET_VOXEL_3D_TYPED( coefs[3], (Real), volume, i , j+1, k+1 );
261 GET_VOXEL_3D_TYPED( coefs[4], (Real), volume, i+1, j , k );
262 GET_VOXEL_3D_TYPED( coefs[5], (Real), volume, i+1, j , k+1 );
263 GET_VOXEL_3D_TYPED( coefs[6], (Real), volume, i+1, j+1, k );
264 GET_VOXEL_3D_TYPED( coefs[7], (Real), volume, i+1, j+1, k+1 );
265 }
266 else
267 {
268 outside_value = convert_value_to_voxel( volume, outside_value );
269
270 i = FLOOR( x );
271 j = FLOOR( y );
272 k = FLOOR( z );
273
274 c = 0;
275 for_less( dx, 0, 2 )
276 for_less( dy, 0, 2 )
277 for_less( dz, 0, 2 )
278 {
279 if( i + dx >= 0 && i + dx < sizes[0] &&
280 j + dy >= 0 && j + dy < sizes[1] &&
281 k + dz >= 0 && k + dz < sizes[2] )
282 {
283 GET_VOXEL_3D_TYPED( coefs[c], (Real), volume, i+dx, j+dy, k+dz);
284 }
285 else
286 coefs[c] = outside_value;
287 ++c;
288 }
289 }
290
291 u = x - (Real) i;
292 v = y - (Real) j;
293 w = z - (Real) k;
294
295 /*--- get the 4 differences in the u direction */
296
297 du00 = coefs[4] - coefs[0];
298 du01 = coefs[5] - coefs[1];
299 du10 = coefs[6] - coefs[2];
300 du11 = coefs[7] - coefs[3];
301
302 /*--- reduce to a 2D problem, by interpolating in the u direction */
303
304 c00 = coefs[0] + u * du00;
305 c01 = coefs[1] + u * du01;
306 c10 = coefs[2] + u * du10;
307 c11 = coefs[3] + u * du11;
308
309 /*--- get the 2 differences in the v direction for the 2D problem */
310
311 dv0 = c10 - c00;
312 dv1 = c11 - c01;
313
314 /*--- reduce 2D to a 1D problem, by interpolating in the v direction */
315
316 c0 = c00 + v * dv0;
317 c1 = c01 + v * dv1;
318
319 /*--- get the 1 difference in the w direction for the 1D problem */
320
321 dw = c1 - c0;
322
323 /*--- if the value is desired, interpolate in 1D to get the value */
324
325 if( value != NULL )
326 {
327 *value = convert_voxel_to_value( volume, c0 + w * dw );
328 }
329
330 /*--- if the derivatives are desired, compute them */
331
332 if( derivs != NULL )
333 {
334 if( volume->real_range_set )
335 scale_factor = volume->real_value_scale;
336 else
337 scale_factor = 1.0;
338
339 /*--- reduce the 2D u derivs to 1D */
340
341 du0 = INTERPOLATE( v, du00, du10 );
342 du1 = INTERPOLATE( v, du01, du11 );
343
344 /*--- interpolate the 1D problems in w, or for Z deriv, just use dw */
345
346 derivs[X] = scale_factor * INTERPOLATE( w, du0, du1 );
347 derivs[Y] = scale_factor * INTERPOLATE( w, dv0, dv1 );
348 derivs[Z] = scale_factor * dw;
349 }
350 }
351
352 /* ----------------------------- MNI Header -----------------------------------
353 @NAME : interpolate_volume
354 @INPUT : n_dims - number of dimensions to interpolate
355 parameters[] - 0 to 1 parameters for each dim
356 n_values - number of values to interpolate
357 degree - degree of interpolation, 4 == cubic
358 coefs - [degree*degree*degree... *n_values] coeficients
359 @OUTPUT : values - pass back values
360 first_deriv - pass first derivs [n_values][n_dims]
361 second_deriv - pass back values [n_values][n_dims][n_dims]
362 @RETURNS :
363 @DESCRIPTION: Computes the interpolation of the box specified by coefs and
364 its derivatives, if necessary.
365 @METHOD :
366 @GLOBALS :
367 @CALLS :
368 @CREATED : Mar. 20, 1995 David MacDonald
369 @MODIFIED :
370 ---------------------------------------------------------------------------- */
371
372 #define MAX_DERIV_SIZE 100
373
interpolate_volume(int n_dims,Real parameters[],int n_values,int degree,Real coefs[],Real values[],Real ** first_deriv,Real *** second_deriv)374 static void interpolate_volume(
375 int n_dims,
376 Real parameters[],
377 int n_values,
378 int degree,
379 Real coefs[],
380 Real values[],
381 Real **first_deriv,
382 Real ***second_deriv )
383 {
384 int v, d, d2, n_derivs, derivs_per_value, mult, mult2;
385 Real fixed_size_derivs[MAX_DERIV_SIZE];
386 Real *derivs;
387
388 /*--- determine how many derivatives should be computed */
389
390 if( second_deriv != NULL )
391 n_derivs = 2;
392 else if( first_deriv != NULL )
393 n_derivs = 1;
394 else
395 n_derivs = 0;
396
397 /*--- compute the total number of values, 1st and 2nd derivatives per val*/
398
399 derivs_per_value = 1;
400 for_less( d, 0, n_dims )
401 derivs_per_value *= 1 + n_derivs;
402
403 /*--- make storage for the spline routines to place the answers */
404
405 if( n_values * derivs_per_value <= MAX_DERIV_SIZE )
406 {
407 derivs = fixed_size_derivs;
408 }
409 else
410 {
411 ALLOC( derivs, n_values * derivs_per_value );
412 }
413
414 /*--- evaluate the interpolating spline */
415
416 evaluate_interpolating_spline( n_dims, parameters, degree, n_values, coefs,
417 n_derivs, derivs );
418
419 /*--- derivs is now a one dimensional array representing
420 derivs[n_values][1+n_derivs][1+n_derivs]...,
421 where derivs[0][0][0][0]... = the interpolated value for 1st comp,
422 derivs[1][0][0][0]... = the interpolated value for 2nd comp,
423 derivs[n_values-1][0][0][0]... = interpolated value for last,
424 derivs[0][1][0][0]... = derivative of 1st comp wrt x
425 derivs[0][0][1][0]... = derivative of 1st comp wrt y
426 derivs[1][1][0][0]... = derivative of 2nd comp wrt x, etc. */
427
428 if( values != NULL )
429 {
430 for_less( v, 0, n_values )
431 values[v] = derivs[v*derivs_per_value];
432 }
433
434 /*--- fill in the first derivatives, if required */
435
436 if( first_deriv != NULL )
437 {
438 mult = 1;
439 for_down( d, n_dims-1, 0 )
440 {
441 for_less( v, 0, n_values )
442 first_deriv[v][d] = derivs[mult + v*derivs_per_value];
443
444 mult *= 1 + n_derivs;
445 }
446 }
447
448 /*--- fill in the second derivatives, if required */
449
450 if( second_deriv != NULL )
451 {
452 mult = 1;
453 for_down( d, n_dims-1, 0 )
454 {
455 for_less( v, 0, n_values )
456 second_deriv[v][d][d] = derivs[2*mult + v*derivs_per_value];
457
458 mult *= 1 + n_derivs;
459 }
460
461 mult = 1;
462 for_down( d, n_dims-1, 0 )
463 {
464 mult2 = 1;
465
466 for_down( d2, n_dims-1, d+1 )
467 {
468 for_less( v, 0, n_values )
469 {
470 second_deriv[v][d][d2] =
471 derivs[mult+mult2+v*derivs_per_value];
472 second_deriv[v][d2][d] =
473 derivs[mult+mult2+v*derivs_per_value];
474 }
475
476 mult2 *= 1 + n_derivs;
477 }
478
479 mult *= 1 + n_derivs;
480 }
481 }
482
483 if( n_values * derivs_per_value > MAX_DERIV_SIZE )
484 {
485 FREE( derivs );
486 }
487 }
488
489 /* ----------------------------- MNI Header -----------------------------------
490 @NAME : extract_coefficients
491 @INPUT : volume
492 start
493 end
494 inc
495 @OUTPUT : coefs
496 @RETURNS :
497 @DESCRIPTION: Extracts the coefficients from a volume.
498 @METHOD :
499 @GLOBALS :
500 @CALLS :
501 @CREATED : Jun 21, 1995 David MacDonald
502 @MODIFIED :
503 ---------------------------------------------------------------------------- */
504
extract_coefficients(Volume volume,int start[],int end[],Real coefs[],int inc[])505 static void extract_coefficients(
506 Volume volume,
507 int start[],
508 int end[],
509 Real coefs[],
510 int inc[] )
511 {
512 int inc0, inc1, inc2, inc3, inc4;
513 int ind;
514 int start0, start1, start2, start3, start4;
515 int end0, end1, end2, end3, end4, n_dims;
516 int v0, v1, v2, v3, v4;
517
518 /*--- for speed, use non-array variables for the loops */
519
520 start0 = start[0];
521 start1 = start[1];
522 start2 = start[2];
523 start3 = start[3];
524 start4 = start[4];
525
526 end0 = end[0];
527 end1 = end[1];
528 end2 = end[2];
529 end3 = end[3];
530 end4 = end[4];
531
532 inc0 = inc[0];
533 inc1 = inc[1];
534 inc2 = inc[2];
535 inc3 = inc[3];
536 inc4 = inc[4];
537
538 /*--- adjust the inc stride for stepping through coefs to account
539 for the additions of the inner loops */
540
541 n_dims = get_volume_n_dimensions(volume);
542
543 if( n_dims >= 2 )
544 inc0 -= inc1 * (end1 - start1);
545 if( n_dims >= 3 )
546 inc1 -= inc2 * (end2 - start2);
547 if( n_dims >= 4 )
548 inc2 -= inc3 * (end3 - start3);
549 if( n_dims >= 5 )
550 inc3 -= inc4 * (end4 - start4);
551
552 /*--- get the coefs[] from the volume. For speed, do each dimension
553 separately */
554
555 ind = 0;
556
557 switch( n_dims )
558 {
559 case 1:
560 for_less( v0, start0, end0 )
561 {
562 GET_VALUE_1D_TYPED( coefs[ind], (Real), volume, v0 );
563 ind += inc0;
564 }
565 break;
566
567 case 2:
568 for_less( v0, start0, end0 )
569 {
570 for_less( v1, start1, end1 )
571 {
572 GET_VALUE_2D_TYPED( coefs[ind], (Real), volume, v0, v1 );
573 ind += inc1;
574 }
575 ind += inc0;
576 }
577 break;
578
579 case 3:
580 for_less( v0, start0, end0 )
581 {
582 for_less( v1, start1, end1 )
583 {
584 for_less( v2, start2, end2 )
585 {
586 GET_VALUE_3D_TYPED( coefs[ind], (Real), volume, v0, v1, v2);
587 ind += inc2;
588 }
589 ind += inc1;
590 }
591 ind += inc0;
592 }
593 break;
594
595 case 4:
596 for_less( v0, start0, end0 )
597 {
598 for_less( v1, start1, end1 )
599 {
600 for_less( v2, start2, end2 )
601 {
602 for_less( v3, start3, end3 )
603 {
604 GET_VALUE_4D_TYPED( coefs[ind], (Real),
605 volume, v0, v1, v2, v3 );
606 ind += inc3;
607 }
608 ind += inc2;
609 }
610 ind += inc1;
611 }
612 ind += inc0;
613 }
614 break;
615
616 case 5:
617 for_less( v0, start0, end0 )
618 {
619 for_less( v1, start1, end1 )
620 {
621 for_less( v2, start2, end2 )
622 {
623 for_less( v3, start3, end3 )
624 {
625 for_less( v4, start4, end4 )
626 {
627 GET_VALUE_5D_TYPED( coefs[ind], (Real), volume,
628 v0, v1, v2, v3, v4 );
629 ind += inc4;
630 }
631 ind += inc3;
632 }
633 ind += inc2;
634 }
635 ind += inc1;
636 }
637 ind += inc0;
638 }
639 break;
640 }
641 }
642
643 static Real interpolation_tolerance = 0.0;
644
645 /* ----------------------------- MNI Header -----------------------------------
646 @NAME : set_volume_interpolation_tolerance
647 @INPUT : tolerance
648 @OUTPUT :
649 @RETURNS :
650 @DESCRIPTION: Sets the tolerance which defines how close to a voxel centre we
651 must be in order to just pass back the voxel value, rather than
652 interpolating.
653 @METHOD :
654 @GLOBALS :
655 @CALLS :
656 @CREATED : Apr. 11, 1996 David MacDonald
657 @MODIFIED :
658 ---------------------------------------------------------------------------- */
659
set_volume_interpolation_tolerance(Real tolerance)660 VIOAPI void set_volume_interpolation_tolerance(
661 Real tolerance )
662 {
663 interpolation_tolerance = tolerance;
664 }
665
666 /* ----------------------------- MNI Header -----------------------------------
667 @NAME : evaluate_volume
668 @INPUT : volume
669 voxel
670 interpolating_dimensions - whether each dimension is interpolated
671 degrees_continuity
672 use_linear_at_edge
673 outside_value
674 @OUTPUT : values
675 first_deriv
676 second_deriv
677 @RETURNS :
678 @DESCRIPTION: Takes a voxel space position and evaluates the value within
679 the volume by nearest_neighbour, linear, quadratic, or
680 cubic interpolation. degrees_continuity == 2 corresponds to
681 cubic, 1 for quadratic, etc.
682 If first_deriv is not a null pointer, then the first derivatives
683 are passed back. Similarly for the second_deriv.
684 If use_linear_at_edge is TRUE, then near the boundaries, either
685 linear or nearest neighbour interpolation is used, even if cubic
686 is specified by the degrees_continuity.
687 If use_linear_at_edge is FALSE, then the 'outside_value' is used
688 to provide coefficients for outside the volume, and the degree
689 specified by degrees_continuity is used.
690
691 Each dimension may or may not be interpolated, specified by the
692 interpolating_dimensions parameter. For instance, a 4D volume
693 of x,y,z,RGB may be interpolated in 3D (x,y,z) for each of the
694 3 RGB components, with one call to evaluate_volume.
695 @CREATED : Mar 1993 David MacDonald
696 @MODIFIED :
697 ---------------------------------------------------------------------------- */
698
699 #define MAX_COEF_SPACE 1000
700
evaluate_volume(Volume volume,Real voxel[],BOOLEAN interpolating_dimensions[],int degrees_continuity,BOOLEAN use_linear_at_edge,Real outside_value,Real values[],Real ** first_deriv,Real *** second_deriv)701 VIOAPI int evaluate_volume(
702 Volume volume,
703 Real voxel[],
704 BOOLEAN interpolating_dimensions[],
705 int degrees_continuity,
706 BOOLEAN use_linear_at_edge,
707 Real outside_value,
708 Real values[],
709 Real **first_deriv,
710 Real ***second_deriv )
711 {
712 int inc[MAX_DIMENSIONS];
713 int start_index, spline_degree;
714 int next_d;
715 int n, v, d, dim, n_values, sizes[MAX_DIMENSIONS], n_dims;
716 int start[MAX_DIMENSIONS], n_interp_dims;
717 int end[MAX_DIMENSIONS];
718 int interp_dims[MAX_DIMENSIONS];
719 int n_coefs;
720 Real fraction[MAX_DIMENSIONS], bound, *coefs, pos;
721 Real fixed_size_coefs[MAX_COEF_SPACE];
722 BOOLEAN fully_inside, fully_outside, on_grid_point;
723
724 n_dims = get_volume_n_dimensions(volume);
725
726 /*--- check for the common case trilinear interpolation of 1 value */
727
728 if( n_dims == 3 && degrees_continuity == 0 && second_deriv == NULL &&
729 (interpolating_dimensions == NULL ||
730 interpolating_dimensions[0] &&
731 interpolating_dimensions[1] &&
732 interpolating_dimensions[2]) )
733 {
734 Real *deriv;
735
736 if( first_deriv == NULL )
737 deriv = NULL;
738 else
739 deriv = first_deriv[0];
740
741 trilinear_interpolate( volume, voxel, outside_value, &values[0],
742 deriv );
743
744 return( 1 );
745 }
746
747 /*--- check if the degrees continuity is between nearest neighbour
748 and cubic */
749
750 if( degrees_continuity < -1 || degrees_continuity > 2 )
751 {
752 print_error( "Warning: evaluate_volume(), degrees invalid: %d\n",
753 degrees_continuity );
754 degrees_continuity = 0;
755 }
756
757 get_volume_sizes( volume, sizes );
758
759 /*--- check if we are near a voxel centre, if so just use nearest neighbour
760 and avoid the expensive interpolation, unless we need derivatives */
761
762 if( interpolation_tolerance > 0.0 &&
763 first_deriv == NULL && second_deriv == NULL )
764 {
765 on_grid_point = TRUE;
766 for_less( d, 0, n_dims )
767 {
768 if( interpolating_dimensions == NULL || interpolating_dimensions[d])
769 {
770 pos = (Real) ROUND( voxel[d] );
771 if( voxel[d] < pos - interpolation_tolerance ||
772 voxel[d] > pos + interpolation_tolerance )
773 {
774 on_grid_point = FALSE;
775 break;
776 }
777 }
778 }
779 if( on_grid_point )
780 degrees_continuity = -1;
781 }
782
783 bound = (Real) degrees_continuity / 2.0;
784
785 /*--- if we must use linear interpolation near the boundaries, then
786 check if we are near the boundaries, and adjust the degrees_continuity
787 accordingly */
788
789 if( use_linear_at_edge && degrees_continuity >= 0 )
790 {
791 for_less( d, 0, n_dims )
792 {
793 if( interpolating_dimensions == NULL || interpolating_dimensions[d])
794 {
795 while( degrees_continuity >= 0 &&
796 (voxel[d] < bound ||
797 voxel[d] > (Real) sizes[d] - 1.0 - bound ||
798 bound == (Real) sizes[d] - 1.0 - bound ) )
799 {
800 --degrees_continuity;
801 if( degrees_continuity == 1 )
802 degrees_continuity = 0;
803 bound = (Real) degrees_continuity / 2.0;
804 }
805 }
806 }
807 }
808
809 /*--- now check which dimensions are being interpolated. Also, compute
810 how many values must be interpolated, which are all the values not
811 in the interpolated dimensions */
812
813 n_interp_dims = 0;
814 n_values = 1;
815 n_coefs = 1;
816 spline_degree = degrees_continuity + 2;
817
818 fully_inside = TRUE;
819 fully_outside = FALSE;
820
821 for_less( d, 0, n_dims )
822 {
823 if( interpolating_dimensions == NULL || interpolating_dimensions[d])
824 {
825 interp_dims[n_interp_dims] = d;
826 pos = voxel[d] - bound;
827 start[d] = FLOOR( pos );
828 fraction[n_interp_dims] = pos - (Real) start[d];
829
830 if( voxel[d] == (Real) sizes[d] - 1.0 - bound )
831 {
832 --start[d];
833 fraction[n_interp_dims] = 1.0;
834 }
835
836 end[d] = start[d] + spline_degree;
837 n_coefs *= spline_degree;
838
839 if( start[d] < 0 || end[d] > sizes[d] )
840 {
841 fully_inside = FALSE;
842
843 if( end[d] <= 0 || start[d] >= sizes[d] )
844 {
845 fully_outside = TRUE;
846 break;
847 }
848 }
849
850 ++n_interp_dims;
851 }
852 else
853 n_values *= sizes[d];
854 }
855
856 /*--- check if the first derivatives are uncomputable */
857
858 if( first_deriv != NULL && (fully_outside || degrees_continuity < 0) )
859 {
860 for_less( v, 0, n_values )
861 for_less( d, 0, n_interp_dims )
862 first_deriv[v][d] = 0.0;
863 }
864
865 /*--- check if the second derivatives are uncomputable */
866
867 if( second_deriv != NULL && (fully_outside || degrees_continuity < 1) )
868 {
869 for_less( v, 0, n_values )
870 for_less( d, 0, n_interp_dims )
871 for_less( dim, 0, n_interp_dims )
872 second_deriv[v][d][dim] = 0.0;
873 }
874
875 /*--- check if the values are uncomputable, i.e., outside volume */
876
877 if( fully_outside )
878 {
879 if( values != NULL )
880 {
881 for_less( v, 0, n_values )
882 values[v] = outside_value;
883 }
884
885 return( n_values );
886 }
887
888 /*--- add the non-interpolated dimensions to the list of dimensions, in
889 order, after the interpolated dimensions */
890
891 n = 0;
892 for_less( d, 0, n_dims )
893 {
894 if( interpolating_dimensions != NULL && !interpolating_dimensions[d] )
895 {
896 interp_dims[n_interp_dims+n] = d;
897 start[d] = 0;
898 end[d] = sizes[d];
899 ++n;
900 }
901 }
902
903 /*--- make room for the coeficients */
904
905 if( n_values * n_coefs > MAX_COEF_SPACE )
906 {
907 ALLOC( coefs, n_values * n_coefs );
908 }
909 else
910 coefs = fixed_size_coefs;
911
912 /*--- figure out the offset within coefs. If we are inside, the offset
913 is zero, since all coefs must be filled in. If we are partially
914 inside, set the offset to the first coef within the volume. */
915
916 if( !fully_inside )
917 {
918 /*--- compute the increments in the coefs[] array for each dimension,
919 in order to simulate a multidimensional array with a single dim
920 array, coefs */
921
922 inc[interp_dims[n_dims-1]] = 1;
923 for_down( d, n_dims-2, 0 )
924 {
925 next_d = interp_dims[d+1];
926 inc[interp_dims[d]] = inc[next_d] * (end[next_d] - start[next_d]);
927 }
928
929 start_index = 0;
930 for_less( d, 0, n_dims )
931 {
932 if( start[d] < 0 )
933 {
934 start_index += -start[d] * inc[d];
935 start[d] = 0;
936 }
937
938 if( end[d] > sizes[d] )
939 end[d] = sizes[d];
940 }
941
942 for_less( v, 0, n_values * n_coefs )
943 coefs[v] = outside_value;
944
945 /*--- get the necessary coeficients from the volume */
946
947 extract_coefficients( volume, start, end, &coefs[start_index], inc );
948 }
949 else
950 {
951 /*--- get the necessary coeficients from the volume */
952
953 for_less( d, n_dims, MAX_DIMENSIONS )
954 {
955 start[d] = 0;
956 end[d] = 0;
957 }
958
959 get_volume_value_hyperslab( volume,
960 start[0], start[1], start[2], start[3], start[4],
961 end[0] - start[0], end[1] - start[1],
962 end[2] - start[2], end[3] - start[3],
963 end[4] - start[4], coefs );
964 }
965
966 /*--- now that we have the coeficients, do the interpolation */
967
968 switch( degrees_continuity )
969 {
970 case -1: /*--- nearest neighbour interpolation */
971 for_less( v, 0, n_values )
972 values[v] = coefs[v];
973 break;
974
975 case 0:
976 case 1:
977 case 2:
978 interpolate_volume( n_interp_dims, fraction, n_values,
979 spline_degree, coefs,
980 values, first_deriv, second_deriv );
981 break;
982 }
983
984 if( n_values * n_coefs > MAX_COEF_SPACE )
985 {
986 FREE( coefs );
987 }
988
989 return( n_values );
990 }
991
992 /* ----------------------------- MNI Header -----------------------------------
993 @NAME : evaluate_volume_in_world
994 @INPUT : volume
995 x
996 y
997 z
998 degrees_continuity - 0 = linear, 2 = cubic
999 use_linear_at_edge
1000 outside_value
1001 @OUTPUT : values
1002 deriv_x
1003 deriv_y
1004 deriv_z
1005 deriv_xx
1006 deriv_xy
1007 deriv_xz
1008 deriv_yy
1009 deriv_yz
1010 deriv_zz
1011 @RETURNS :
1012 @DESCRIPTION: Takes a world space position and evaluates the value within
1013 the volume.
1014 If deriv_x is not a null pointer, then the 3 derivatives are
1015 passed back. If deriv_xx is not null, then the 6 second
1016 derivatives are passed back. If the volume is 3D, then only
1017 one value, and one derivative per deriv_x,etc. is passed back.
1018 If the volume has more than 3 dimensions, say 5 dimensions, with
1019 dimensions 3 and 4 being the non-spatial dimensions, then there
1020 will be sizes[3] * sizes[4] values passed back. The derivatives
1021 are converted to world space.
1022 @CREATED : Mar 1993 David MacDonald
1023 @MODIFIED :
1024 ---------------------------------------------------------------------------- */
1025
evaluate_volume_in_world(Volume volume,Real x,Real y,Real z,int degrees_continuity,BOOLEAN use_linear_at_edge,Real outside_value,Real values[],Real deriv_x[],Real deriv_y[],Real deriv_z[],Real deriv_xx[],Real deriv_xy[],Real deriv_xz[],Real deriv_yy[],Real deriv_yz[],Real deriv_zz[])1026 VIOAPI void evaluate_volume_in_world(
1027 Volume volume,
1028 Real x,
1029 Real y,
1030 Real z,
1031 int degrees_continuity,
1032 BOOLEAN use_linear_at_edge,
1033 Real outside_value,
1034 Real values[],
1035 Real deriv_x[],
1036 Real deriv_y[],
1037 Real deriv_z[],
1038 Real deriv_xx[],
1039 Real deriv_xy[],
1040 Real deriv_xz[],
1041 Real deriv_yy[],
1042 Real deriv_yz[],
1043 Real deriv_zz[] )
1044 {
1045 Real ignore;
1046 Real voxel[MAX_DIMENSIONS];
1047 Real **first_deriv, ***second_deriv;
1048 Real t[N_DIMENSIONS][MAX_DIMENSIONS];
1049 int c, d, dim, v, n_values, n_dims, axis;
1050 int sizes[MAX_DIMENSIONS], dims_interpolated[N_DIMENSIONS];
1051 BOOLEAN interpolating_dimensions[MAX_DIMENSIONS];
1052
1053 /*--- convert the world space to a voxel coordinate */
1054
1055 convert_world_to_voxel( volume, x, y, z, voxel );
1056 get_volume_sizes( volume, sizes );
1057
1058 /*--- initialize all dimensions to not being interpolated */
1059
1060 n_dims = get_volume_n_dimensions( volume );
1061 for_less( d, 0, n_dims )
1062 interpolating_dimensions[d] = FALSE;
1063
1064 /*--- set each spatial dimension to being interpolated */
1065
1066 for_less( d, 0, N_DIMENSIONS )
1067 {
1068 axis = volume->spatial_axes[d];
1069 if( axis < 0 )
1070 {
1071 // print_error(
1072 // "evaluate_volume_in_world(): must have 3 spatial axes.\n" );
1073 // return;
1074 } else {
1075 interpolating_dimensions[axis] = TRUE;
1076 }
1077
1078 }
1079
1080 /*--- compute the number of values, the product of the sizes of the
1081 non-interpolating dimensions */
1082
1083 n_values = 1;
1084 for_less( d, 0, n_dims )
1085 {
1086 if( !interpolating_dimensions[d] )
1087 n_values *= sizes[d];
1088 }
1089
1090 /*--- make room for the first derivative, if necessary */
1091
1092 if( deriv_x != NULL )
1093 {
1094 ALLOC2D( first_deriv, n_values, N_DIMENSIONS );
1095 }
1096 else
1097 first_deriv = NULL;
1098
1099 /*--- make room for the second derivative, if necessary */
1100
1101 if( deriv_xx != NULL )
1102 {
1103 ALLOC3D( second_deriv, n_values, N_DIMENSIONS, N_DIMENSIONS );
1104 }
1105 else
1106 second_deriv = NULL;
1107
1108 /*--- evaluate the volume and derivatives in voxel space */
1109
1110 n_values = evaluate_volume( volume, voxel, interpolating_dimensions,
1111 degrees_continuity, use_linear_at_edge, outside_value,
1112 values, first_deriv, second_deriv );
1113
1114 /*--- if the derivative is desired, convert the voxel derivative
1115 to world space */
1116
1117 if( deriv_x != NULL || deriv_xx != NULL )
1118 {
1119 /*--- figure out the dimensions interpolated, in order */
1120
1121 dim = 0;
1122 for_less( d, 0, n_dims )
1123 {
1124 if( interpolating_dimensions[d] )
1125 {
1126 dims_interpolated[dim] = d;
1127 ++dim;
1128 }
1129 }
1130 }
1131
1132 if( deriv_x != NULL )
1133 {
1134 for_less( v, 0, n_values ) /*--- convert the deriv of each value */
1135 {
1136 /*--- get the voxel coordinates of the first derivative */
1137
1138 for_less( c, 0, N_DIMENSIONS )
1139 voxel[dims_interpolated[c]] = first_deriv[v][c];
1140
1141 /*--- convert the voxel-space derivative to a world derivative */
1142
1143 convert_voxel_normal_vector_to_world( volume, voxel,
1144 &deriv_x[v], &deriv_y[v], &deriv_z[v] );
1145 }
1146
1147 FREE2D( first_deriv );
1148 }
1149
1150 /*--- if the derivative is desired, convert the voxel derivative
1151 to world space */
1152
1153 if( deriv_xx != (Real *) 0 )
1154 {
1155 for_less( v, 0, n_values ) /*--- convert the deriv of each value */
1156 {
1157 /*--- get the voxel coordinates of the first derivative */
1158
1159 for_less( dim, 0, N_DIMENSIONS )
1160 {
1161 for_less( c, 0, N_DIMENSIONS )
1162 voxel[dims_interpolated[c]] = second_deriv[v][dim][c];
1163
1164 /*--- convert the voxel-space derivative to a world derivative*/
1165
1166 convert_voxel_normal_vector_to_world( volume, voxel,
1167 &t[X][dims_interpolated[dim]],
1168 &t[Y][dims_interpolated[dim]],
1169 &t[Z][dims_interpolated[dim]] );
1170 }
1171
1172 /*--- now convert the results to world */
1173
1174 convert_voxel_normal_vector_to_world( volume, t[X],
1175 &deriv_xx[v], &ignore, &ignore );
1176
1177 convert_voxel_normal_vector_to_world( volume, t[Y],
1178 &deriv_xy[v], &deriv_yy[v], &ignore );
1179
1180 convert_voxel_normal_vector_to_world( volume, t[Z],
1181 &deriv_xz[v], &deriv_yz[v], &deriv_zz[v] );
1182 }
1183
1184 FREE3D( second_deriv );
1185 }
1186 }
1187