Mercurial > hg > minc-tools
changeset 772:32437b7bfce4
*** empty log message ***
| author | david <david> |
|---|---|
| date | Wed, 22 Feb 1995 14:42:48 +0000 |
| parents | 61cf793fe122 |
| children | e1e0f33f27d3 |
| files | volume_io/Volumes/evaluate.c |
| diffstat | 1 files changed, 378 insertions(+), 811 deletions(-) [+] |
line wrap: on
line diff
--- a/volume_io/Volumes/evaluate.c +++ b/volume_io/Volumes/evaluate.c @@ -1,22 +1,282 @@ #include <internal_volume_io.h> -#include <vols.h> -#include <numerical.h> + +private void interpolate_volume( + int n_dims, + Real parameters[], + int n_values, + int degree, + Real coefs[], + Real values[], + Real **first_deriv, + Real ***second_deriv ) +{ + int v, d, d2, n_derivs, derivs_per_value, mult, mult2; + Real *derivs; + + if( second_deriv != NULL ) + n_derivs = 2; + else if( first_deriv != NULL ) + n_derivs = 1; + else + n_derivs = 0; + + derivs_per_value = 1; + for_less( d, 0, n_dims ) + derivs_per_value *= 1 + n_derivs; + + ALLOC( derivs, n_values * derivs_per_value ); + + evaluate_interpolating_spline( n_dims, parameters, degree, n_values, coefs, + n_derivs, derivs ); + + /* --- derivs is now a one dimensional array representing + derivs[n_values][1+n_derivs][1+n_derivs]... */ + + if( values != NULL ) + { + for_less( v, 0, n_values ) + values[v] = derivs[v*derivs_per_value]; + } -private void trilinear_interpolate_volume( - Volume volume, - Real x, - Real y, - Real z, - BOOLEAN x_use_higher, - BOOLEAN y_use_higher, - BOOLEAN z_use_higher, - Real *value, - Real *deriv_x, - Real *deriv_y, - Real *deriv_z ); + if( first_deriv != NULL ) + { + mult = 1; + for_down( d, n_dims-1, 0 ) + { + for_less( v, 0, n_values ) + first_deriv[v][d] = derivs[mult + v*derivs_per_value]; + + mult *= 1 + n_derivs; + } + } + + if( second_deriv != NULL ) + { + mult = 1; + for_down( d, n_dims-1, 0 ) + { + for_less( v, 0, n_values ) + second_deriv[v][d][d] = derivs[2*mult + v*derivs_per_value]; + + mult *= 1 + n_derivs; + } + + mult = 1; + for_down( d, n_dims-1, 0 ) + { + mult2 = 1; + + for_down( d2, n_dims-1, d+1 ) + { + for_less( v, 0, n_values ) + { + second_deriv[v][d][d2] = + derivs[mult+mult2+v*derivs_per_value]; + second_deriv[v][d2][d] = + derivs[mult+mult2+v*derivs_per_value]; + } + + mult2 *= 1 + n_derivs; + } + + mult *= 1 + n_derivs; + } + } + + FREE( derivs ); +} /* ----------------------------- MNI Header ----------------------------------- -@NAME : evaluate_3D_volume_in_world +@NAME : evaluate_volume +@INPUT : volume + x + y + z +@OUTPUT : value + deriv_x + deriv_y + deriv_z +@RETURNS : +@DESCRIPTION: Takes a voxel space position and evaluates the value within + the volume by nearest_neighbour, linear, quadratic, or + cubic interpolation. + If first_deriv is not a null pointer, then the first derivatives + are passed back. Similarly for the second_deriv. +@CREATED : Mar 1993 David MacDonald +@MODIFIED : +---------------------------------------------------------------------------- */ + +public int evaluate_volume( + Volume volume, + Real voxel[], + BOOLEAN interpolating_dimensions[], + int degrees_continuity, + Real values[], + Real **first_deriv, + Real ***second_deriv, + BOOLEAN *inside ) +{ + int n, v, d, dim, n_values, sizes[MAX_DIMENSIONS], n_dims, ind; + int start[MAX_DIMENSIONS], n_interp_dims; + int end[MAX_DIMENSIONS]; + int interp_dims[MAX_DIMENSIONS]; + int n_coefs; + int vi[MAX_DIMENSIONS]; + Real value, fraction[MAX_DIMENSIONS], bound, *coefs; + + if( degrees_continuity > 2 ) + { + print( "Warning: evaluate_volume(), degrees too large: %d\n", + degrees_continuity ); + degrees_continuity = 2; + } + + n_dims = get_volume_n_dimensions(volume); + get_volume_sizes( volume, sizes ); + n_values = 1; + n_interp_dims = 0; + + for_less( d, 0, n_dims ) + { + if( interpolating_dimensions == NULL || interpolating_dimensions[d] ) + { + if( degrees_continuity > -2 && + (voxel[d] < (Real) degrees_continuity * 0.5 || + voxel[d] > (Real)(sizes[d]-1) - (Real)degrees_continuity*0.5) ) + { + --degrees_continuity; + for( ; degrees_continuity > -2; --degrees_continuity ) + { + if( voxel[d] >= (Real) degrees_continuity * 0.5 && + voxel[d] <= (Real)(sizes[d]-1) - + (Real)degrees_continuity * 0.5 ) + { + break; + } + } + } + ++n_interp_dims; + } + else + { + n_values *= sizes[d]; + } + } + + if( inside != NULL ) + *inside = degrees_continuity >= -1; + + if( degrees_continuity < -1 ) + { + value = get_volume_real_min( volume ); + + for_less( v, 0, n_values ) + values[v] = value; + } + + if( degrees_continuity < 0 && first_deriv != NULL ) + { + for_less( v, 0, n_values ) + for_less( d, 0, n_interp_dims ) + first_deriv[v][d] = 0.0; + } + + if( degrees_continuity < 1 && second_deriv != NULL ) + { + for_less( v, 0, n_values ) + for_less( d, 0, n_interp_dims ) + for_less( dim, 0, n_interp_dims ) + second_deriv[v][d][dim] = 0.0; + } + + if( degrees_continuity > -2 ) + { + bound = degrees_continuity / 2.0; + n_interp_dims = 0; + n_coefs = 1; + + for_less( d, 0, n_dims ) + { + if( interpolating_dimensions == NULL || interpolating_dimensions[d]) + { + interp_dims[n_interp_dims] = d; + n_coefs *= 2 + degrees_continuity; + if( voxel[d] >= (Real) sizes[d] - 1.0 - bound ) + { + start[n_interp_dims] = sizes[d] - degrees_continuity - 2; + fraction[n_interp_dims] = 1.0; + } + else + { + if( voxel[d] < bound ) + voxel[d] = bound; + + start[n_interp_dims] = (int) ( voxel[d] - bound ); + fraction[n_interp_dims] = FRACTION( voxel[d] - bound ); + } + + end[n_interp_dims] = start[d] + degrees_continuity + 2; + ++n_interp_dims; + } + } + + n = 0; + for_less( d, 0, n_dims ) + { + if( interpolating_dimensions != NULL && + !interpolating_dimensions[d] ) + { + interp_dims[n_interp_dims+n] = d; + start[n_interp_dims+n] = 0; + end[n_interp_dims+n] = sizes[d]; + ++n; + } + } + + for_less( d, n_dims, MAX_DIMENSIONS ) + { + start[d] = 0; + end[d] = 1; + interp_dims[d] = d; + } + + ALLOC( coefs, n_values * n_coefs ); + + ind = 0; + for_less( vi[interp_dims[0]], start[0], end[0] ) + for_less( vi[interp_dims[1]], start[1], end[1] ) + for_less( vi[interp_dims[2]], start[2], end[2] ) + for_less( vi[interp_dims[3]], start[3], end[3] ) + for_less( vi[interp_dims[4]], start[4], end[4] ) + { + GET_VALUE( coefs[ind], volume, vi[0], vi[1], vi[2], vi[3], vi[4] ); + ++ind; + } + + switch( degrees_continuity ) + { + case -1: + for_less( v, 0, n_values ) + values[v] = coefs[v]; + break; + + case 0: + case 1: + case 2: + interpolate_volume( n_interp_dims, fraction, n_values, + degrees_continuity + 2, coefs, + values, first_deriv, second_deriv ); + break; + } + + FREE( coefs ); + } + + return( n_values ); +} + +/* ----------------------------- MNI Header ----------------------------------- +@NAME : evaluate_volume_in_world @INPUT : volume x y @@ -35,833 +295,140 @@ @MODIFIED : ---------------------------------------------------------------------------- */ -public void evaluate_3D_volume_in_world( +public BOOLEAN evaluate_volume_in_world( Volume volume, Real x, Real y, Real z, int degrees_continuity, - Real *value, - 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 ) + 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[] ) { - Real ignore, dxx, dxy, dxz, dyy, dyz, dzz; + Real ignore; Real voxel[MAX_DIMENSIONS]; - Real txx, txy, txz; - Real tyx, tyy, tyz; - Real tzx, tzy, tzz; - - if( get_volume_n_dimensions(volume) != 3 ) - { - handle_internal_error( - "evaluate_3D_volume_in_world: volume must be 3D.\n" ); - } - - convert_world_to_voxel( volume, x, y, z, voxel ); - - evaluate_3D_volume( volume, voxel[X], voxel[Y], voxel[Z], - degrees_continuity, value, deriv_x, deriv_y, deriv_z, - deriv_xx, deriv_xy, deriv_xz, - deriv_yy, deriv_yz, deriv_zz ); - - if( deriv_x != (Real *) 0 ) - { - convert_voxel_normal_vector_to_world( volume, - *deriv_x, *deriv_y, *deriv_z, - deriv_x, deriv_y, deriv_z ); - } - - if( deriv_xx != (Real *) 0 ) - { - dxx = *deriv_xx; - dxy = *deriv_xy; - dxz = *deriv_xz; - dyy = *deriv_yy; - dyz = *deriv_yz; - dzz = *deriv_zz; - convert_voxel_normal_vector_to_world( volume, - dxx, dxy, dxz, - &txx, &txy, &txz ); - convert_voxel_normal_vector_to_world( volume, - dxy, dyy, dyz, - &tyx, &tyy, &tyz ); - convert_voxel_normal_vector_to_world( volume, - dxz, dyz, dzz, - &tzx, &tzy, &tzz ); - - convert_voxel_normal_vector_to_world( volume, - txx, tyx, tzx, - deriv_xx, &ignore, &ignore ); - convert_voxel_normal_vector_to_world( volume, - txy, tyy, tzy, - deriv_xy, deriv_yy, &ignore ); - convert_voxel_normal_vector_to_world( volume, - txz, tyz, tzz, - deriv_xz, deriv_yz, deriv_zz ); - } -} - -public void special_evaluate_3D_volume_in_world( - Volume volume, - Real x, - Real y, - Real z, - BOOLEAN x_use_higher, - BOOLEAN y_use_higher, - BOOLEAN z_use_higher, - Real *value, - Real *deriv_x, - Real *deriv_y, - Real *deriv_z ) -{ - Real voxel[MAX_DIMENSIONS]; + Real **first_deriv, ***second_deriv; + Real t[MAX_DIMENSIONS][MAX_DIMENSIONS]; + int c, d, dim, v, n_values, n_interp_dims, n_dims, axis; int sizes[MAX_DIMENSIONS]; + BOOLEAN interpolating_dimensions[MAX_DIMENSIONS], inside; convert_world_to_voxel( volume, x, y, z, voxel ); + get_volume_sizes( volume, sizes ); - get_volume_sizes( volume, sizes ); - if( x < -0.5 || x > (Real) sizes[X] - 0.5 || - y < -0.5 || y > (Real) sizes[Y] - 0.5 || - z < -0.5 || z > (Real) sizes[Z] - 0.5 ) + n_dims = get_volume_n_dimensions( volume ); + for_less( d, 0, n_dims ) + interpolating_dimensions[d] = FALSE; + + for_less( d, 0, N_DIMENSIONS ) { - *value = get_volume_voxel_min( volume ); - if( deriv_x != (Real *) NULL ) + axis = volume->spatial_axes[d]; + if( axis < 0 ) { - *deriv_x = 0.0; - *deriv_y = 0.0; - *deriv_z = 0.0; + print("evaluate_volume_in_world(): must have 3 spatial axes.\n"); + return( FALSE ); } - return; + + interpolating_dimensions[axis] = TRUE; + } + + + n_values = 1; + for_less( d, 0, n_dims ) + { + if( !interpolating_dimensions[d] ) + n_values *= sizes[d]; } - trilinear_interpolate_volume( volume, voxel[X], voxel[Y], voxel[Z], - x_use_higher, y_use_higher, z_use_higher, - value, deriv_x, deriv_y, deriv_z ); - - if( deriv_x != (Real *) 0 ) + if( deriv_x != NULL ) { - convert_voxel_normal_vector_to_world( volume, - *deriv_x, *deriv_y, *deriv_z, - deriv_x, deriv_y, deriv_z ); + ALLOC2D( first_deriv, n_values, N_DIMENSIONS ); } -} + else + first_deriv = NULL; -/* ----------------------------- MNI Header ----------------------------------- -@NAME : evaluate_3D_slice_in_world -@INPUT : volume - x - y - z -@OUTPUT : value - deriv_x - deriv_y - deriv_xx - deriv_xy - deriv_yy -@RETURNS : -@DESCRIPTION: Takes a world space position and evaluates the value within - the volume by bilinear interpolation within the slice. - If deriv_x is not a null pointer, then the derivatives are passed - back. -@CREATED : Mar 1993 David MacDonald -@MODIFIED : ----------------------------------------------------------------------------- */ + if( deriv_xx != NULL ) + { + ALLOC3D( second_deriv, n_values, N_DIMENSIONS, N_DIMENSIONS ); + } + else + second_deriv = NULL; + + n_values = evaluate_volume( volume, voxel, interpolating_dimensions, + degrees_continuity, values, first_deriv, second_deriv, + &inside ); -public void evaluate_3D_slice_in_world( - Volume volume, - Real x, - Real y, - Real z, - Real *value, - Real *deriv_x, - Real *deriv_y, - Real *deriv_xx, - Real *deriv_xy, - Real *deriv_yy ) -{ - Real ignore, dxx, dxy, dyy; - Real voxel[MAX_DIMENSIONS]; - Real txx, txy, txz; - Real tyx, tyy, tyz; - Real tzx, tzy, tzz; + if( deriv_x != NULL ) + { + for_less( v, 0, n_values ) + { + n_interp_dims = 0; + for_less( c, 0, n_dims ) + { + if( interpolating_dimensions[c] ) + { + voxel[c] = first_deriv[v][n_interp_dims]; + ++n_interp_dims; + } + else + voxel[c] = 0.0; + } - if( get_volume_n_dimensions(volume) != 3 ) - { - handle_internal_error( - "evaluate_3D_slice_in_world: volume must be 3D.\n" ); - } - - convert_world_to_voxel( volume, x, y, z, voxel ); - - evaluate_3D_slice( volume, voxel[X], voxel[Y], voxel[Z], - value, deriv_x, deriv_y, - deriv_xx, deriv_xy, deriv_yy ); - - if( deriv_x != (Real *) 0 ) - { - convert_voxel_normal_vector_to_world( volume, - *deriv_x, *deriv_y, 0.0, - deriv_x, deriv_y, &ignore ); + convert_voxel_normal_vector_to_world( volume, voxel, + &deriv_x[v], &deriv_y[v], &deriv_z[v] ); + } } if( deriv_xx != (Real *) 0 ) { - dxx = *deriv_xx; - dxy = *deriv_xy; - dyy = *deriv_yy; - convert_voxel_normal_vector_to_world( volume, - dxx, dxy, 0.0, - &txx, &txy, &txz ); - convert_voxel_normal_vector_to_world( volume, - dxy, dyy, 0.0, - &tyx, &tyy, &tyz ); - convert_voxel_normal_vector_to_world( volume, - 0.0, 0.0, 0.0, - &tzx, &tzy, &tzz ); - - convert_voxel_normal_vector_to_world( volume, - txx, tyx, tzx, - deriv_xx, &ignore, &ignore ); - convert_voxel_normal_vector_to_world( volume, - txy, tyy, tzy, - deriv_xy, deriv_yy, &ignore ); - } -} - -/* ----------------------------- MNI Header ----------------------------------- -@NAME : triconstant_interpolate_volume -@INPUT : volume - x - y - z -@OUTPUT : value - deriv_x - deriv_y - deriv_z -@RETURNS : -@DESCRIPTION: Returns the value within the volume, assuming constant voxels, - (step function). Derivative is approximated by neighbours. -@METHOD : -@GLOBALS : -@CALLS : -@CREATED : 1993 David MacDonald -@MODIFIED : ----------------------------------------------------------------------------- */ - -private void triconstant_interpolate_volume( - Volume volume, - Real x, - Real y, - Real z, - Real *value, - Real *deriv_x, - Real *deriv_y, - Real *deriv_z ) -{ - int i, j, k, sizes[MAX_DIMENSIONS]; - Real prev, next, dx, dy, dz, voxel_value; - - get_volume_sizes( volume, sizes ); - - i = ROUND( x ); - if( i == sizes[X] ) - i = sizes[X]-1; - j = ROUND( y ); - if( j == sizes[Y] ) - j = sizes[Y]-1; - k = ROUND( z ); - if( k == sizes[Z] ) - k = sizes[Z]-1; - - GET_VALUE_3D( voxel_value, volume, i, j, k ); - - if( value != (Real *) 0 ) - *value = voxel_value; - - if( deriv_x != (Real *) NULL ) - { - /* --- get derivative wrt x */ - - dx = 0; - if( i == 0 ) - prev = voxel_value; - else + for_less( v, 0, n_values ) { - GET_VALUE_3D( prev, volume, i-1, j, k ); - ++dx; - } - - if( i == sizes[X]-1 ) - next = voxel_value; - else - { - GET_VALUE_3D( next, volume, i+1, j, k ); - ++dx; - } - - if( dx == 0 ) - *deriv_x = 0.0; - else - *deriv_x = (next - prev) / (Real) dx; - - /* --- get derivative wrt y */ - - dy = 0; - if( j == 0 ) - prev = voxel_value; - else - { - GET_VALUE_3D( prev, volume, i, j-1, k ); - ++dy; - } - - if( j == sizes[Y]-1 ) - next = voxel_value; - else - { - GET_VALUE_3D( next, volume, i, j+1, k ); - ++dy; - } + for_less( dim, 0, N_DIMENSIONS ) + { + n_interp_dims = 0; + for_less( c, 0, n_dims ) + { + if( interpolating_dimensions[c] ) + { + voxel[c] = second_deriv[v][dim][n_interp_dims]; + ++n_interp_dims; + } + else + voxel[c] = 0.0; + } - if( dy == 0 ) - *deriv_y = 0.0; - else - *deriv_y = (next - prev) / (Real) dy; - - /* --- get derivative wrt z */ - - dz = 0; - if( k == 0 ) - prev = voxel_value; - else - { - GET_VALUE_3D( prev, volume, i, j, k-1 ); - ++dz; - } - - if( k == sizes[Z]-1 ) - next = voxel_value; - else - { - GET_VALUE_3D( next, volume, i, j, k+1 ); - ++dz; - } - - if( dz == 0 ) - *deriv_z = 0.0; - else - *deriv_z = (next - prev) / (Real) dz; - } -} - -/* ----------------------------- MNI Header ----------------------------------- -@NAME : trilinear_interpolate_volume -@INPUT : volume - x - y - z -@OUTPUT : value - deriv_x - deriv_y - deriv_z -@RETURNS : -@DESCRIPTION: Returns the value within the volume, assuming trilinear - interpolation. -@METHOD : -@GLOBALS : -@CALLS : -@CREATED : 1993 David MacDonald -@MODIFIED : ----------------------------------------------------------------------------- */ + convert_voxel_normal_vector_to_world( volume, voxel, + &t[dim][volume->spatial_axes[X]], + &t[dim][volume->spatial_axes[Y]], + &t[dim][volume->spatial_axes[Z]] ); + } + + for_less( c, 0, n_dims ) + voxel[c] = t[c][volume->spatial_axes[X]]; -private void trilinear_interpolate_volume( - Volume volume, - Real x, - Real y, - Real z, - BOOLEAN x_use_higher, - BOOLEAN y_use_higher, - BOOLEAN z_use_higher, - Real *value, - Real *deriv_x, - Real *deriv_y, - Real *deriv_z ) -{ - int i, j, k, sizes[MAX_DIMENSIONS]; - Real u, v, w; - Real c000, c001, c010, c011, c100, c101, c110, c111; - Real c00, c01, c10, c11; - Real c0, c1; - Real du00, du01, du10, du11, du0, du1; - Real dv0, dv1; - - get_volume_sizes( volume, sizes ); + convert_voxel_normal_vector_to_world( volume, voxel, + &deriv_xx[v], &ignore, &ignore ); + + for_less( c, 0, n_dims ) + voxel[c] = t[c][volume->spatial_axes[Y]]; - if( x == (Real) (sizes[X]-1) ) - { - i = sizes[X]-2; - u = 1.0; - } - else - { - i = (int) x; - u = FRACTION( x ); - if( !x_use_higher && u == 0.0 && i > 0 ) - { - --i; - u = 1.0; - } - } + convert_voxel_normal_vector_to_world( volume, voxel, + &deriv_xy[v], &deriv_yy[v], &ignore ); + + for_less( c, 0, n_dims ) + voxel[c] = t[c][volume->spatial_axes[Z]]; - if( y == (Real) (sizes[Y]-1) ) - { - j = sizes[Y]-2; - v = 1.0; - } - else - { - j = (int) y; - v = FRACTION( y ); - if( !y_use_higher && v == 0.0 && j > 0 ) - { - --j; - v = 1.0; - } - } - - if( z == (Real) (sizes[Z]-1) ) - { - k = sizes[Z]-2; - w = 1.0; - } - else - { - k = (int) z; - w = FRACTION( z ); - if( !z_use_higher && w == 0.0 && k > 0 ) - { - --k; - w = 1.0; + convert_voxel_normal_vector_to_world( volume, voxel, + &deriv_xz[v], &deriv_yz[v], &deriv_zz[v] ); } } - GET_VALUE_3D( c000, volume, i, j, k ); - GET_VALUE_3D( c001, volume, i, j, k+1 ); - GET_VALUE_3D( c010, volume, i, j+1, k ); - GET_VALUE_3D( c011, volume, i, j+1, k+1 ); - GET_VALUE_3D( c100, volume, i+1, j, k ); - GET_VALUE_3D( c101, volume, i+1, j, k+1 ); - GET_VALUE_3D( c110, volume, i+1, j+1, k ); - GET_VALUE_3D( c111, volume, i+1, j+1, k+1 ); - - du00 = c100 - c000; - du01 = c101 - c001; - du10 = c110 - c010; - du11 = c111 - c011; - - c00 = c000 + u * du00; - c01 = c001 + u * du01; - c10 = c010 + u * du10; - c11 = c011 + u * du11; - - dv0 = c10 - c00; - dv1 = c11 - c01; - - c0 = c00 + v * dv0; - c1 = c01 + v * dv1; - - if( value != (Real *) 0 ) - *value = INTERPOLATE( w, c0, c1 ); - - if( deriv_x != (Real *) 0 ) - { - du0 = INTERPOLATE( v, du00, du10 ); - du1 = INTERPOLATE( v, du01, du11 ); - - *deriv_x = INTERPOLATE( w, du0, du1 ); - *deriv_y = INTERPOLATE( w, dv0, dv1 ); - *deriv_z = (c1 - c0); - } -} - -private void trivar_interpolate_volume( - Volume volume, - Real x, - Real y, - Real z, - int degree, - Real *value, - 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 ) -{ - int i, j, k, tu, tv, tw, ind, n_derivs; - Real u, v, w, bound; - Real coefs[4*4*4]; - Real derivs[3*3*3]; - int sizes[MAX_DIMENSIONS]; - - get_volume_sizes( volume, sizes ); - - bound = (degree - 2) / 2.0; - - if( x >= (Real) sizes[X] - 1.0 - bound ) - { - i = sizes[X] - degree; - u = 1.0; - } - else - { - if( x < bound ) - x = bound; - - i = (int) ( x - bound ); - u = FRACTION( x - bound ); - } - - if( y >= (Real) sizes[Y] - 1.0 - bound ) - { - j = sizes[Y] - degree; - v = 1.0; - } - else - { - if( y < bound ) - y = bound; - - j = (int) ( y - bound ); - v = FRACTION( y - bound ); - } - - if( z >= (Real) sizes[Z] - 1.0 - bound ) - { - k = sizes[Z] - degree; - w = 1.0; - } - else - { - if( z < bound ) - z = bound; - - k = (int) ( z - bound ); - w = FRACTION( z - bound ); - } - - ind = 0; - for_less( tu, 0, degree ) - for_less( tv, 0, degree ) - for_less( tw, 0, degree ) - { - GET_VALUE_3D( coefs[ind], volume, i+tu, j+tv, k+tw ); - ++ind; - } - - if( deriv_xx != NULL ) - n_derivs = 2; - else if( deriv_x != NULL ) - n_derivs = 1; - else - n_derivs = 0; - - evaluate_trivariate_catmull_spline( u, v, w, degree, coefs, n_derivs, - derivs ); - - if( deriv_xx != NULL ) - { - *deriv_xx = derivs[IJK(2,0,0,n_derivs+1,n_derivs+1)]; - *deriv_xy = derivs[IJK(1,1,0,n_derivs+1,n_derivs+1)]; - *deriv_xz = derivs[IJK(1,0,1,n_derivs+1,n_derivs+1)]; - *deriv_yy = derivs[IJK(0,2,0,n_derivs+1,n_derivs+1)]; - *deriv_yz = derivs[IJK(0,1,1,n_derivs+1,n_derivs+1)]; - *deriv_zz = derivs[IJK(0,0,2,n_derivs+1,n_derivs+1)]; - } - - if( deriv_x != NULL ) - { - *deriv_x = derivs[IJK(1,0,0,n_derivs+1,n_derivs+1)]; - *deriv_y = derivs[IJK(0,1,0,n_derivs+1,n_derivs+1)]; - *deriv_z = derivs[IJK(0,0,1,n_derivs+1,n_derivs+1)]; - } - - if( value != NULL ) - *value = derivs[IJK(0,0,0,n_derivs+1,n_derivs+1)]; + return( inside ); } - -private void bicubic_interpolate_volume( - Volume volume, - Real x, - Real y, - Real z, - Real *value, - Real *deriv_x, - Real *deriv_y, - Real *deriv_xx, - Real *deriv_xy, - Real *deriv_yy ) -{ - int i, j, k, tu, tv, ind, n_derivs; - Real u, v; - Real coefs[4*4]; - Real derivs[3*3]; - int sizes[MAX_DIMENSIONS]; - - get_volume_sizes( volume, sizes ); - - if( x >= (Real) sizes[X] - 1.5 ) - { - i = sizes[X]-2; - u = 1.0; - } - else - { - i = (int) x; - u = FRACTION( x ); - } - - if( y >= (Real) sizes[Y] - 1.5 ) - { - j = sizes[Y]-2; - v = 1.0; - } - else - { - j = (int) y; - v = FRACTION( y ); - } - - if( z == (Real) sizes[Z] - 1.5 ) - { - k = sizes[Z]-2; - } - else - { - k = (int) z; - } - - ind = 0; - for_less( tu, 0, 4 ) - for_less( tv, 0, 4 ) - { - GET_VALUE_3D( coefs[ind], volume, i-1+tu, j-1+tv, k ); - ++ind; - } - - if( deriv_xx != NULL ) - n_derivs = 2; - else if( deriv_x != NULL ) - n_derivs = 1; - else - n_derivs = 0; - - evaluate_bivariate_catmull_spline( u, v, 4, coefs, n_derivs, derivs ); - - if( deriv_xx != NULL ) - { - *deriv_xx = derivs[IJ(2,0,n_derivs+1)]; - *deriv_xy = derivs[IJ(1,1,n_derivs+1)]; - *deriv_yy = derivs[IJ(0,2,n_derivs+1)]; - } - - if( deriv_x != NULL ) - { - *deriv_x = derivs[IJ(1,0,n_derivs+1)]; - *deriv_y = derivs[IJ(0,1,n_derivs+1)]; - } - - if( value != NULL ) - *value = derivs[IJ(0,0,n_derivs+1)]; -} - -/* ----------------------------- MNI Header ----------------------------------- -@NAME : evaluate_3D_volume -@INPUT : volume - x - y - z -@OUTPUT : value - deriv_x - deriv_y - deriv_z -@RETURNS : -@DESCRIPTION: Takes a voxel space position and evaluates the value within - the volume by trilinear interpolation. - If deriv_x is not a null pointer, then the 3 derivatives are passed - back. -@CREATED : Mar 1993 David MacDonald -@MODIFIED : ----------------------------------------------------------------------------- */ - -public void evaluate_3D_volume( - Volume volume, - Real x, - Real y, - Real z, - int degrees_continuity, - Real *value, - 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 ) -{ - int sizes[MAX_DIMENSIONS]; - - if( get_volume_n_dimensions(volume) != 3 ) - { - handle_internal_error( "evaluate_3D_volume: volume must be 3D.\n" ); - } - - get_volume_sizes( volume, sizes ); - - if( x < -0.5 || x > (Real) sizes[X] - 0.5 || - y < -0.5 || y > (Real) sizes[Y] - 0.5 || - z < -0.5 || z > (Real) sizes[Z] - 0.5 ) - { - *value = get_volume_voxel_min( volume ); - if( deriv_x != (Real *) NULL ) - { - *deriv_x = 0.0; - *deriv_y = 0.0; - *deriv_z = 0.0; - } - if( deriv_xx != (Real *) NULL ) - { - *deriv_xx = 0.0; - *deriv_xy = 0.0; - *deriv_xz = 0.0; - *deriv_yy = 0.0; - *deriv_yz = 0.0; - *deriv_zz = 0.0; - } - return; - } - - if( x < (Real) degrees_continuity * 0.5 || - x > (Real) (sizes[X]-1) - (Real) degrees_continuity * 0.5 || - y < (Real) degrees_continuity * 0.5 || - y > (Real) (sizes[Y]-1) - (Real) degrees_continuity * 0.5 || - z < (Real) degrees_continuity * 0.5 || - z > (Real) (sizes[Z]-1) - (Real) degrees_continuity * 0.5 ) - { - degrees_continuity = -1; - } - - if( degrees_continuity < 1 && deriv_xx != (Real *) NULL ) - { - *deriv_xx = 0.0; - *deriv_xy = 0.0; - *deriv_xz = 0.0; - *deriv_yy = 0.0; - *deriv_yz = 0.0; - *deriv_zz = 0.0; - } - - switch( degrees_continuity ) - { - case -1: - triconstant_interpolate_volume( volume, x, y, z, value, - deriv_x, deriv_y, deriv_z ); - break; - - case 0: - trilinear_interpolate_volume( volume, x, y, z, - TRUE, TRUE, TRUE, value, - deriv_x, deriv_y, deriv_z ); - break; - - case 1: - case 2: - trivar_interpolate_volume( volume, x, y, z, degrees_continuity+2, - value, - deriv_x, deriv_y, deriv_z, - deriv_xx, deriv_xy, deriv_xz, - deriv_yy, deriv_yz, deriv_zz ); - break; - - default: - handle_internal_error( "evaluate_3D_volume: invalid continuity" ); - } -} - -/* ----------------------------- MNI Header ----------------------------------- -@NAME : evaluate_3D_slice -@INPUT : volume - x - y - z -@OUTPUT : value - deriv_x - deriv_y - deriv_xx - deriv_xy - deriv_yy -@RETURNS : -@DESCRIPTION: Takes a voxel position and evaluates the value within - the volume by bilinear interpolation within the slice. - If deriv_x is not a null pointer, then the - derivatives are passed back. -@CREATED : Mar 1993 David MacDonald -@MODIFIED : ----------------------------------------------------------------------------- */ - -public void evaluate_3D_slice( - Volume volume, - Real x, - Real y, - Real z, - Real *value, - Real *deriv_x, - Real *deriv_y, - Real *deriv_xx, - Real *deriv_xy, - Real *deriv_yy ) -{ - int sizes[MAX_DIMENSIONS]; - - if( get_volume_n_dimensions(volume) != 3 ) - { - handle_internal_error( "evaluate_3D_slice: volume must be 3D.\n" ); - } - - get_volume_sizes( volume, sizes ); - - if( x < 1.0 || x > (Real) sizes[X] - 2.0 || - y < 1.0 || y > (Real) sizes[Y] - 2.0 || - z < 1.0 || z > (Real) sizes[Z] - 2.0 ) - { - *value = get_volume_voxel_min( volume ); - if( deriv_x != (Real *) NULL ) - { - *deriv_x = 0.0; - *deriv_y = 0.0; - } - if( deriv_xx != (Real *) NULL ) - { - *deriv_xx = 0.0; - *deriv_xy = 0.0; - *deriv_yy = 0.0; - } - return; - } - - bicubic_interpolate_volume( volume, x, y, z, value, - deriv_x, deriv_y, - deriv_xx, deriv_xy, deriv_yy ); -}