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 #ifdef HAVE_CONFIG_H
15 #include "config.h"
16 #endif /*HAVE_CONFIG_H*/
17
18 #include <internal_volume_io.h>
19
20 #define DEGREES_CONTINUITY 2 /* -1 = Nearest; 0 = Linear; 1 = Quadratic; 2 = Cubic interpolation */
21 #define SPLINE_DEGREE ((DEGREES_CONTINUITY) + 2)
22
23 #define N_COMPONENTS VIO_N_DIMENSIONS /* displacement vector has 3 components */
24
25 #define FOUR_DIMS 4
26
27 #ifdef USE_NEWTONS_METHOD
28 #define INVERSE_FUNCTION_TOLERANCE 0.01
29 #define INVERSE_DELTA_TOLERANCE 1.0e-5
30 #define MAX_INVERSE_ITERATIONS 20
31 #endif
32
33 static void evaluate_grid_volume(
34 VIO_Volume volume,
35 VIO_Real x,
36 VIO_Real y,
37 VIO_Real z,
38 int degrees_continuity,
39 VIO_Real values[],
40 VIO_Real deriv_x[],
41 VIO_Real deriv_y[],
42 VIO_Real deriv_z[] );
43
44 /* ----------------------------- MNI Header -----------------------------------
45 @NAME : grid_transform_point
46 @INPUT : transform
47 x
48 y
49 z
50 @OUTPUT : x_transformed
51 y_transformed
52 z_transformed
53 @RETURNS :
54 @DESCRIPTION: Applies a grid transform to the point
55 @METHOD :
56 @GLOBALS :
57 @CALLS :
58 @CREATED : Feb. 21, 1995 David MacDonald
59 @MODIFIED :
60 ---------------------------------------------------------------------------- */
61
grid_transform_point(VIO_General_transform * transform,VIO_Real x,VIO_Real y,VIO_Real z,VIO_Real * x_transformed,VIO_Real * y_transformed,VIO_Real * z_transformed)62 VIOAPI VIO_Status grid_transform_point(
63 VIO_General_transform *transform,
64 VIO_Real x,
65 VIO_Real y,
66 VIO_Real z,
67 VIO_Real *x_transformed,
68 VIO_Real *y_transformed,
69 VIO_Real *z_transformed )
70 {
71 VIO_Real displacements[N_COMPONENTS];
72 VIO_Volume volume;
73
74 /* --- the volume that defines the transform is an offset vector,
75 so evaluate the volume at the given position and add the
76 resulting offset to the given position */
77
78 if(!transform->displacement_volume)
79 return VIO_ERROR;
80
81 volume = (VIO_Volume) transform->displacement_volume;
82
83 evaluate_grid_volume( volume, x, y, z, DEGREES_CONTINUITY, displacements,
84 NULL, NULL, NULL );
85
86 *x_transformed = x + displacements[VIO_X];
87 *y_transformed = y + displacements[VIO_Y];
88 *z_transformed = z + displacements[VIO_Z];
89
90 return VIO_OK;
91 }
92
93 #ifdef USE_NEWTONS_METHOD
94 /* ----------------------------- MNI Header -----------------------------------
95 @NAME : forward_function
96 @INPUT : function_data - contains transform info
97 parameters - x,y,z position
98 @OUTPUT : values - where x,y,z, maps to
99 derivatives - the 3 by 3 derivatives of the mapping
100 @RETURNS :
101 @DESCRIPTION: This function does the same thing as grid_transform_point(),
102 but also gets derivatives. This function is passed to the
103 newton function solution routine to perform the inverse mapping
104 of the grid transformation.
105 @METHOD :
106 @GLOBALS :
107 @CALLS :
108 @CREATED : Feb. , 1995 David MacDonald
109 @MODIFIED :
110 ---------------------------------------------------------------------------- */
111
forward_function(void * function_data,VIO_Real parameters[],VIO_Real values[],VIO_Real ** derivatives)112 static void forward_function(
113 void *function_data,
114 VIO_Real parameters[],
115 VIO_Real values[],
116 VIO_Real **derivatives )
117 {
118 int c;
119 VIO_General_transform *transform;
120 VIO_Real deriv_x[N_COMPONENTS], deriv_y[N_COMPONENTS];
121 VIO_Real deriv_z[N_COMPONENTS];
122 VIO_Volume volume;
123
124 transform = (VIO_General_transform *) function_data;
125
126 /* --- store the offset vector in values[0-2] */
127
128 volume = (VIO_Volume) transform->displacement_volume;
129
130 evaluate_grid_volume( volume, parameters[X], parameters[Y], parameters[Z],
131 DEGREES_CONTINUITY, values,
132 deriv_x, deriv_y, deriv_z );
133
134 for_less( c, 0, N_COMPONENTS )
135 {
136 values[c] += parameters[c]; /* to get x',y',z', add offset to x,y,z */
137
138 /*--- given the derivatives of the offset, compute the
139 derivatives of (x,y,z) + offset, with respect to x,y,z */
140
141 derivatives[c][X] = deriv_x[c];
142 derivatives[c][Y] = deriv_y[c];
143 derivatives[c][Z] = deriv_z[c];
144
145 derivatives[c][c] += 1.0; /* deriv of (x,y,z) w.r.t. x or y or z */
146 }
147 }
148
149 /* ----------------------------- MNI Header -----------------------------------
150 @NAME : grid_inverse_transform_point
151 @INPUT : transform
152 x
153 y
154 z
155 @OUTPUT : x_transformed
156 y_transformed
157 z_transformed
158 @RETURNS : VIO_Status
159 @DESCRIPTION: Applies the inverse grid transform to the point. This is done
160 by using newton-rhapson steps to find the point which maps to
161 the parameters (x,y,z).
162 @METHOD :
163 @GLOBALS :
164 @CALLS :
165 @CREATED : Feb. 21, 1995 David MacDonald
166 @MODIFIED :
167 ---------------------------------------------------------------------------- */
168
169 /* ---------------------------------------------------------------------------
170
171 There are two different versions of the grid inverse function. I would
172 have hoped that my version worked best, since it uses first derivatives and
173 Newton's method. However, Louis' version seems to work better, perhaps since
174 it matches the code he uses in minctracc to generate the grid transforms.
175
176 - David MacDonald
177
178 ---------------------------------------------------------------------------- */
179
grid_inverse_transform_point(VIO_General_transform * transform,VIO_Real x,VIO_Real y,VIO_Real z,VIO_Real * x_transformed,VIO_Real * y_transformed,VIO_Real * z_transformed)180 VIOAPI VIO_Status grid_inverse_transform_point(
181 VIO_General_transform *transform,
182 VIO_Real x,
183 VIO_Real y,
184 VIO_Real z,
185 VIO_Real *x_transformed,
186 VIO_Real *y_transformed,
187 VIO_Real *z_transformed )
188 {
189 VIO_Real solution[VIO_N_DIMENSIONS];
190 VIO_Real initial_guess[VIO_N_DIMENSIONS];
191 VIO_Real desired_values[VIO_N_DIMENSIONS];
192
193 /* --- fill in the initial guess */
194
195 initial_guess[X] = x;
196 initial_guess[Y] = y;
197 initial_guess[Z] = z;
198
199 /* --- define what the desired function values are */
200
201 desired_values[X] = x;
202 desired_values[Y] = y;
203 desired_values[Z] = z;
204
205 /* --- find the x,y,z that are mapped to the desired values */
206
207 if( newton_root_find( VIO_N_DIMENSIONS, forward_function,
208 (void *) transform,
209 initial_guess, desired_values,
210 solution, INVERSE_FUNCTION_TOLERANCE,
211 INVERSE_DELTA_TOLERANCE, MAX_INVERSE_ITERATIONS ))
212 {
213 *x_transformed = solution[X];
214 *y_transformed = solution[Y];
215 *z_transformed = solution[Z];
216 return VIO_OK;
217 }
218 else /* --- if no solution found, not sure what is reasonable to return */
219 {
220 *x_transformed = x;
221 *y_transformed = y;
222 *z_transformed = z;
223 return VIO_ERROR;
224 }
225 return VIO_ERROR;
226 }
227 #endif
228
229 /* ----------------------------- MNI Header -----------------------------------
230 @NAME : grid_inverse_transform_point_with_input_steps
231 @INPUT : transform
232 x
233 y
234 z
235 input_volume_steps
236 @OUTPUT : x_transformed
237 y_transformed
238 z_transformed
239 @RETURNS :
240 @DESCRIPTION: Transforms the point by the inverse of the grid transform.
241 Approximates the solution using a simple iterative step
242 method.
243 @METHOD :
244 @GLOBALS :
245 @CALLS :
246 @CREATED : 1993? Louis Collins
247 @MODIFIED : 1994 David MacDonald
248 @MODIFIED : 2013 June 10, Matthijs van Eede, added the possibility to pass
249 along the step sizes of the input file which is
250 being resampled to determine the appropriate error
251 margin (ftol)
252 ---------------------------------------------------------------------------- */
253
grid_inverse_transform_point_with_input_steps(VIO_General_transform * transform,VIO_Real x,VIO_Real y,VIO_Real z,VIO_Real * input_volume_steps,VIO_Real * x_transformed,VIO_Real * y_transformed,VIO_Real * z_transformed)254 VIOAPI VIO_Status grid_inverse_transform_point_with_input_steps(
255 VIO_General_transform *transform,
256 VIO_Real x,
257 VIO_Real y,
258 VIO_Real z,
259 VIO_Real *input_volume_steps,
260 VIO_Real *x_transformed,
261 VIO_Real *y_transformed,
262 VIO_Real *z_transformed )
263 {
264 #define NUMBER_TRIES 10
265 int tries;
266 VIO_Real best_x, best_y, best_z;
267 VIO_Real tx, ty, tz;
268 VIO_Real gx, gy, gz;
269 VIO_Real error_x, error_y, error_z, error, smallest_e;
270 VIO_Real ftol;
271 VIO_Status status=VIO_ERROR;
272 int sizes[VIO_MAX_DIMENSIONS];
273 VIO_Real steps[VIO_MAX_DIMENSIONS];
274 VIO_Volume volume;
275 short d, vector_dim = -1;
276 int i;
277
278 if((status=grid_transform_point( transform, x, y, z, &tx, &ty, &tz ))!=VIO_OK)
279 return status;
280 tx = x - (tx - x);
281 ty = y - (ty - y);
282 tz = z - (tz - z);
283
284 if((status=grid_transform_point( transform, tx, ty, tz, &gx, &gy, &gz ))!=VIO_OK)
285 return status;
286
287 error_x = x - gx;
288 error_y = y - gy;
289 error_z = z - gz;
290
291 tries = 0;
292
293 smallest_e = VIO_FABS(error_x) + VIO_FABS(error_y) + VIO_FABS(error_z);
294 best_x = tx;
295 best_y = ty;
296 best_z = tz;
297
298 // Adapt ftol to grid step sizes. For 1mm stx volume with grid 4mm, we
299 // are using ftol=0.05 (=4mm/80). For histology data at grid 0.125mm,
300 // then use ftol=0.125/80=0.0015625, which is fine on 0.01mm volume.
301 // ftol = 0.05; // good for MNI space at 1mm (grid 2mm or 4mm)
302 // ftol = 0.001; // acceptable for big brain slice (grid at 0.125, voxel at 0.01mm)
303
304 // This is ok too:
305 // Make the error a fraction of the initial residual.
306 // ftol = 0.05 * smallest_e + 0.0001;
307
308
309 volume = (VIO_Volume) transform->displacement_volume;
310 get_volume_sizes( volume, sizes );
311 get_volume_separations( volume, steps );
312
313 /*--- find which of 4 dimensions is the vector dimension */
314
315 for_less( vector_dim, 0, FOUR_DIMS ) {
316 for_less( d, 0, VIO_N_DIMENSIONS ) {
317 if( volume->spatial_axes[d] == vector_dim ) break;
318 }
319 if( d == VIO_N_DIMENSIONS ) break;
320 }
321
322 // If we have information about the step sizes of the input volume that is
323 // being resampled, base the tolerance on those instead of the step sizes
324 // of the deformation grid; there can be a significant difference.
325
326 ftol = -1.0;
327 if( input_volume_steps != NULL){
328 for_less( i, 0, 3 ) {
329 if( ftol < 0 ) ftol = input_volume_steps[i];
330 if( input_volume_steps[i] < ftol ) ftol = input_volume_steps[i];
331 }
332 }
333 else {
334 for_less( d, 0, FOUR_DIMS ) {
335 if( d == vector_dim ) continue;
336 if( sizes[d] == 1 ) continue;
337 if( ftol < 0 ) ftol = steps[d];
338 if( steps[d] < ftol ) ftol = steps[d];
339 }
340 }
341
342 ftol = ftol / 80.0;
343 if( ftol > 0.05 ) ftol = 0.05; // just to be sure for large grids
344
345 while( ++tries < NUMBER_TRIES && smallest_e > ftol ) {
346 tx += 0.95 * error_x;
347 ty += 0.95 * error_y;
348 tz += 0.95 * error_z;
349
350 if((status=grid_transform_point( transform, tx, ty, tz, &gx, &gy, &gz ))!=VIO_OK)
351 return status;
352
353 error_x = x - gx;
354 error_y = y - gy;
355 error_z = z - gz;
356
357 error = VIO_FABS(error_x) + VIO_FABS(error_y) + VIO_FABS(error_z);
358
359 if( error < smallest_e ) {
360 smallest_e = error;
361 best_x = tx;
362 best_y = ty;
363 best_z = tz;
364 }
365 }
366
367 *x_transformed = best_x;
368 *y_transformed = best_y;
369 *z_transformed = best_z;
370 return VIO_OK;
371 }
372
373 /* ----------------------------- MNI Header -----------------------------------
374 @NAME : grid_inverse_transform_point
375 @INPUT : transform
376 x
377 y
378 z
379 @OUTPUT : x_transformed
380 y_transformed
381 z_transformed
382 @RETURNS :
383 @DESCRIPTION: Transforms the point by the inverse of the grid transform.
384 Approximates the solution using a simple iterative step
385 method.
386 @METHOD :
387 @GLOBALS :
388 @CALLS :
389 @CREATED : 1993? Louis Collins
390 @MODIFIED : 1994 David MacDonald
391 @MODIFIED :
392 ---------------------------------------------------------------------------- */
393
grid_inverse_transform_point(VIO_General_transform * transform,VIO_Real x,VIO_Real y,VIO_Real z,VIO_Real * x_transformed,VIO_Real * y_transformed,VIO_Real * z_transformed)394 VIOAPI VIO_Status grid_inverse_transform_point(
395 VIO_General_transform *transform,
396 VIO_Real x,
397 VIO_Real y,
398 VIO_Real z,
399 VIO_Real *x_transformed,
400 VIO_Real *y_transformed,
401 VIO_Real *z_transformed )
402 {
403 return grid_inverse_transform_point_with_input_steps(transform, x, y, z,
404 NULL,
405 x_transformed, y_transformed, z_transformed );
406 }
407
408 /* ----------------------------- MNI Header -----------------------------------
409 @NAME : evaluate_grid_volume
410 @INPUT : volume
411 voxel
412 degrees_continuity
413 @OUTPUT : values
414 derivs (if non-NULL)
415 @RETURNS :
416 @DESCRIPTION: Takes a voxel space position and evaluates the value within
417 the volume by nearest_neighbour, linear, quadratic, or
418 cubic interpolation. Rather than use the generic evaluate_volume
419 function, this special purpose function is a bit faster.
420 @CREATED : Mar. 16, 1995 David MacDonald
421 @MODIFIED :
422 ---------------------------------------------------------------------------- */
423
evaluate_grid_volume(VIO_Volume volume,VIO_Real x,VIO_Real y,VIO_Real z,int degrees_continuity,VIO_Real values[],VIO_Real deriv_x[],VIO_Real deriv_y[],VIO_Real deriv_z[])424 static void evaluate_grid_volume(
425 VIO_Volume volume,
426 VIO_Real x,
427 VIO_Real y,
428 VIO_Real z,
429 int degrees_continuity,
430 VIO_Real values[],
431 VIO_Real deriv_x[],
432 VIO_Real deriv_y[],
433 VIO_Real deriv_z[] )
434 {
435 VIO_Real voxel[VIO_MAX_DIMENSIONS], voxel_vector[VIO_MAX_DIMENSIONS];
436 int inc0, inc1, inc2, inc3, inc[VIO_MAX_DIMENSIONS], derivs_per_value;
437 int ind0, vector_dim;
438 int start0, start1, start2, start3, inc_so_far;
439 int end0, end1, end2, end3;
440 int v0, v1, v2, v3;
441 int v, d, id, sizes[VIO_MAX_DIMENSIONS];
442 int start[VIO_MAX_DIMENSIONS];
443 int end[VIO_MAX_DIMENSIONS];
444 VIO_Real fraction[VIO_MAX_DIMENSIONS], bound, pos;
445 VIO_Real coefs[SPLINE_DEGREE*SPLINE_DEGREE*SPLINE_DEGREE*N_COMPONENTS];
446 VIO_Real values_derivs[N_COMPONENTS + N_COMPONENTS * VIO_N_DIMENSIONS];
447 int is_2dslice = -1;
448
449
450 convert_world_to_voxel( volume, x, y, z, voxel );
451
452 if( get_volume_n_dimensions(volume) != FOUR_DIMS )
453 handle_internal_error( "evaluate_grid_volume" );
454
455 /*--- find which of 4 dimensions is the vector dimension */
456
457 for_less( vector_dim, 0, FOUR_DIMS ) {
458 for_less( d, 0, VIO_N_DIMENSIONS ) {
459 if( volume->spatial_axes[d] == vector_dim )
460 break;
461 }
462 if( d == VIO_N_DIMENSIONS )
463 break;
464 }
465
466 get_volume_sizes( volume, sizes );
467
468 /*--- if a 2-d slice, do best interpolation in the plane */
469
470 for_less( d, 0, FOUR_DIMS ) {
471 if( d == vector_dim ) continue;
472 if( sizes[d] == 1 ) {
473 is_2dslice = d;
474 }
475 }
476
477 bound = (VIO_Real) degrees_continuity / 2.0;
478
479 /*--- if near the edges, reduce the degrees of continuity.
480 This is very important. Doing cubic (with a shifted
481 stencil) near a boundary will cause trouble because
482 the stencil needs to be centered. This is why quadratic
483 is also disabled (not symmetric stencil). CL. */
484
485 for_less( d, 0, FOUR_DIMS ) {
486 if( d == is_2dslice ) continue;
487 if( d == vector_dim ) continue;
488 while( degrees_continuity >= -1 &&
489 (voxel[d] < bound ||
490 voxel[d] > (VIO_Real) sizes[d] - 1.0 - bound ||
491 bound == (VIO_Real) sizes[d] - 1.0 - bound ) ) {
492 --degrees_continuity;
493 if( degrees_continuity == 1 )
494 degrees_continuity = 0;
495 bound = (VIO_Real) degrees_continuity / 2.0;
496 }
497 }
498
499 /*--- check to fill in the first derivative */
500
501 if( degrees_continuity < 0 && deriv_x != NULL ) {
502 for_less( v, 0, N_COMPONENTS ) {
503 deriv_x[v] = 0.0;
504 deriv_y[v] = 0.0;
505 deriv_z[v] = 0.0;
506 }
507 }
508
509 /*--- check if outside */
510
511 for( d = 0; d < FOUR_DIMS; d++ ) {
512 if( d == vector_dim ) continue;
513 if( voxel[d] < -0.5 || voxel[d] > sizes[d]-0.5 ) {
514 for_less( v, 0, N_COMPONENTS ) {
515 values[v] = 0.0;
516 }
517 return;
518 }
519 }
520
521 /*--- determine the starting positions in the volume to grab control
522 vertices */
523
524 id = 0;
525 for_less( d, 0, FOUR_DIMS ) {
526 if( d == vector_dim ) continue;
527 if( d == is_2dslice ) {
528 start[d] = 0;
529 end[d] = 1;
530 } else {
531 pos = voxel[d] - bound;
532 start[d] = VIO_FLOOR( pos );
533 if( start[d] < 0 ) {
534 start[d] = 0;
535 } else if( start[d]+degrees_continuity+1 >= sizes[d] ) {
536 start[d] = sizes[d] - degrees_continuity - 2;
537 }
538 end[d] = start[d] + degrees_continuity + 2;
539 fraction[id] = pos - (double) start[d];
540 ++id;
541 }
542 }
543
544 /*--- create the strides */
545
546 start[vector_dim] = 0;
547 end[vector_dim] = N_COMPONENTS;
548
549 inc_so_far = N_COMPONENTS;
550 for_down( d, FOUR_DIMS-1, 0 ) {
551 if( d != vector_dim ) {
552 inc[d] = inc_so_far;
553 inc_so_far *= ( end[d] - start[d] );
554 }
555 }
556
557 /*--- copy stride arrays to variables for speed */
558
559 inc[vector_dim] = 1;
560
561 start0 = start[0];
562 start1 = start[1];
563 start2 = start[2];
564 start3 = start[3];
565
566 end0 = end[0];
567 end1 = end[1];
568 end2 = end[2];
569 end3 = end[3];
570
571 inc0 = inc[0] - inc[1] * (end1 - start1);
572 inc1 = inc[1] - inc[2] * (end2 - start2);
573 inc2 = inc[2] - inc[3] * (end3 - start3);
574 inc3 = inc[3];
575
576 /*--- extract values from volume */
577
578 ind0 = 0;
579
580 for_less( v0, start0, end0 ) {
581 for_less( v1, start1, end1 ) {
582 for_less( v2, start2, end2 ) {
583 for_less( v3, start3, end3 ) {
584 GET_VALUE_4D_TYPED( coefs[ind0], (VIO_Real), volume,
585 v0, v1, v2, v3 );
586 ind0 += inc3;
587 }
588 ind0 += inc2;
589 }
590 ind0 += inc1;
591 }
592 ind0 += inc0;
593 }
594
595 /*--- interpolate values */
596
597 if( degrees_continuity == -1 ) {
598 for_less( v, 0, N_COMPONENTS )
599 values[v] = coefs[v];
600 } else {
601 if( is_2dslice == -1 ) {
602 evaluate_interpolating_spline( VIO_N_DIMENSIONS, fraction,
603 degrees_continuity + 2,
604 N_COMPONENTS, coefs, 0, values_derivs );
605 } else {
606 evaluate_interpolating_spline( VIO_N_DIMENSIONS-1, fraction,
607 degrees_continuity + 2,
608 N_COMPONENTS, coefs, 0, values_derivs );
609 }
610
611 /*--- extract values and derivatives from values_derivs */
612
613 if( deriv_x != NULL )
614 derivs_per_value = 8;
615 else
616 derivs_per_value = 1;
617
618 for_less( v, 0, N_COMPONENTS ) {
619 values[v] = values_derivs[v*derivs_per_value];
620 }
621
622 if( deriv_x != NULL )
623 {
624 for_less( v, 0, N_COMPONENTS )
625 {
626 id = 0;
627 for_less( d, 0, FOUR_DIMS )
628 {
629 if( d != vector_dim )
630 {
631 voxel_vector[d] = values_derivs[v*8 + (4>>id)];
632 ++id;
633 }
634 else
635 voxel_vector[d] = 0.0;
636 }
637
638 convert_voxel_normal_vector_to_world( volume, voxel_vector,
639 &deriv_x[v], &deriv_y[v], &deriv_z[v] );
640 }
641 }
642 }
643 }
644