/* ---------------------------------------------------------------------------- @COPYRIGHT : Copyright 1993,1994,1995 David MacDonald, McConnell Brain Imaging Centre, Montreal Neurological Institute, McGill University. Permission to use, copy, modify, and distribute this software and its documentation for any purpose and without fee is hereby granted, provided that the above copyright notice appear in all copies. The author and McGill University make no representations about the suitability of this software for any purpose. It is provided "as is" without express or implied warranty. ---------------------------------------------------------------------------- */ #include #ifndef lint static char rcsid[] = "$Header: /software/source/minc/cvsroot/minc/volume_io/Geometry/splines.c,v 1.6 2004/10/04 20:23:51 bert Exp $"; #endif /*--- Weighting functions which define the splines, all of which are interpolation splines, with the exception of the quadratic spline */ static Real constant_coefs[1][1] = { { 1.0 } }; static Real linear_coefs[2][2] = { { 1.0, 0.0 }, { -1.0, 1.0 } }; static Real quadratic_coefs[3][3] = { { 0.5, 0.5, 0.0 }, { -1.0, 1.0, 0.0 }, { 0.5, -1.0, 0.5 } }; static Real cubic_coefs[4][4] = { { 0.0, 1.0, 0.0, 0.0 }, { -0.5, 0.0, 0.5, 0.0 }, { 1.0, -2.5, 2.0, -0.5 }, { -0.5, 1.5, -1.5, 0.5 } }; /* ----------------------------- MNI Header ----------------------------------- @NAME : get_linear_spline_coefs @INPUT : @OUTPUT : coefs 2 by 2 array of coefficients @RETURNS : @DESCRIPTION: Passes back the basis matrix of the linear spline. @METHOD : @GLOBALS : @CALLS : @CREATED : Jan 21, 1995 David MacDonald @MODIFIED : ---------------------------------------------------------------------------- */ VIOAPI void get_linear_spline_coefs( Real **coefs ) { int i, j; for_less( i, 0, 2 ) for_less( j, 0, 2 ) coefs[i][j] = linear_coefs[i][j]; } /* ----------------------------- MNI Header ----------------------------------- @NAME : get_quadratic_spline_coefs @INPUT : @OUTPUT : coefs 3 by 3 array of coefficients @RETURNS : @DESCRIPTION: Passes back the basis matrix of the quadratic spline. @METHOD : @GLOBALS : @CALLS : @CREATED : Jan 21, 1995 David MacDonald @MODIFIED : ---------------------------------------------------------------------------- */ VIOAPI void get_quadratic_spline_coefs( Real **coefs ) { int i, j; for_less( i, 0, 3 ) for_less( j, 0, 3 ) coefs[i][j] = quadratic_coefs[i][j]; } /* ----------------------------- MNI Header ----------------------------------- @NAME : get_cubic_spline_coefs @INPUT : @OUTPUT : coefs 4 by 4 array of coefficients @RETURNS : @DESCRIPTION: Passes back the basis matrix of the cubic interpolating (Catmull-Romm) spline. @METHOD : @GLOBALS : @CALLS : @CREATED : Jan 21, 1995 David MacDonald @MODIFIED : ---------------------------------------------------------------------------- */ VIOAPI void get_cubic_spline_coefs( Real **coefs ) { int i, j; for_less( i, 0, 4 ) for_less( j, 0, 4 ) coefs[i][j] = cubic_coefs[i][j]; } /* ----------------------------- MNI Header ----------------------------------- @NAME : cubic_interpolate @INPUT : u - position to evaluate, between 0 and 1 v0 - four control vertices v1 v2 v3 @OUTPUT : @RETURNS : interpolated value @DESCRIPTION: Performs cubic interpolation, where a value of u = 0 returns v1, a value of u = 1 returns v2, and intermediate values smoothly interpolate. @METHOD : @GLOBALS : @CALLS : @CREATED : Jan 21, 1995 David MacDonald @MODIFIED : ---------------------------------------------------------------------------- */ VIOAPI Real cubic_interpolate( Real u, Real v0, Real v1, Real v2, Real v3 ) { Real coefs[4], value; coefs[0] = v0; coefs[1] = v1; coefs[2] = v2; coefs[3] = v3; evaluate_univariate_interpolating_spline( u, 4, coefs, 0, &value ); return( value ); } /* ----------------------------- MNI Header ----------------------------------- @NAME : evaluate_univariate_interpolating_spline @INPUT : u degree - 2,3,4 for linear, quadratic, or cubic coefs[degree] - control vertices n_derivs - number of derivatives to compute @OUTPUT : derivs - 1 + n_derivs values and derivatives @RETURNS : @DESCRIPTION: Passes back the interpolated value and n_derivs derivatives in the derivs array. @METHOD : @GLOBALS : @CALLS : @CREATED : Jan 21, 1995 David MacDonald @MODIFIED : ---------------------------------------------------------------------------- */ VIOAPI void evaluate_univariate_interpolating_spline( Real u, int degree, Real coefs[], int n_derivs, Real derivs[] ) { evaluate_interpolating_spline( 1, &u, degree, 1, coefs, n_derivs, derivs ); } /* ----------------------------- MNI Header ----------------------------------- @NAME : evaluate_bivariate_interpolating_spline @INPUT : u - position to evaluate v degree - 2,3,4 for linear, quadratic, or cubic coefs - control vertices, size degree * degree n_derivs - number of derivatives to compute @OUTPUT : derivs - (1 + n_derivs) * (1 + n_derivs) values and derivatives @RETURNS : @DESCRIPTION: Passes back the interpolated value and derivatives in the derivs array. derivs is a 1D array that is conceptually 2 dimensional, indexed by dx and dy, where dx and dy range from 0 to n_derivs, indicating which value or derivative. For example 0,0 refers to the interpolated value whereas 1,0 refers to the derivative of the function wrt u. @METHOD : @GLOBALS : @CALLS : @CREATED : Jan 21, 1995 David MacDonald @MODIFIED : ---------------------------------------------------------------------------- */ VIOAPI void evaluate_bivariate_interpolating_spline( Real u, Real v, int degree, Real coefs[], int n_derivs, Real derivs[] ) { Real positions[2]; positions[0] = u; positions[1] = v; evaluate_interpolating_spline( 2, positions, degree, 1, coefs, n_derivs, derivs ); } /* ----------------------------- MNI Header ----------------------------------- @NAME : evaluate_trivariate_interpolating_spline @INPUT : u - position to evaluate v w degree - 2,3,4 for linear, quadratic, or cubic coefs - control vertices, size degree * degree * degree n_derivs - number of derivatives to compute @OUTPUT : derivs - (1 + n_derivs) * (1 + n_derivs) * (1 + n_derivs) values and derivatives @RETURNS : @DESCRIPTION: Passes back the interpolated value and derivatives in the derivs array. derivs is a 1D array that is conceptually 3 dimensional, indexed by dx, dy, and dz, where dx, dy, and dz each range from 0 to n_derivs, indicating which value or derivative. For example 0,0,0 refers to the interpolated value whereas 1,0,1 refers to the derivative of the function wrt u and w. @METHOD : @GLOBALS : @CALLS : @CREATED : Jan 21, 1995 David MacDonald @MODIFIED : ---------------------------------------------------------------------------- */ VIOAPI void evaluate_trivariate_interpolating_spline( Real u, Real v, Real w, int degree, Real coefs[], int n_derivs, Real derivs[] ) { Real positions[3]; positions[0] = u; positions[1] = v; positions[2] = w; evaluate_interpolating_spline( 3, positions, degree, 1, coefs, n_derivs, derivs ); } #define MAX_DIMS 100 /* ----------------------------- MNI Header ----------------------------------- @NAME : evaluate_interpolating_spline @INPUT : n_dims - dimensionality of the spline, >= 1 parameters - u, v, w,... position of spline degree - 2,3,4 for linear, quadratic, or cubic n_values - number of values to interpolate at the point coefs - [n_values]*[degree]*[degree]... control vertices n_derivs - number of derivatives to compute @OUTPUT : derivs - (n_values) * (1 + n_derivs) * (1 + n_derivs) * ... values and derivatives @RETURNS : @DESCRIPTION: Passes back the interpolated value and derivatives in the derivs array. derivs is a 1D array that is conceptually multi-dimensional, indexed by v, dx, dy, dz, etc., where dx, dy, dz, etc. each range from 0 to n_derivs, and v ranges from 0 to n_values-1. For example, if n_dims is 3 and n_values is 4, then the 4D index of derivs[2,0,0,0] refers to the interpolated value of the 3rd component of the 4 valued function. derivs[1,0,1,1] refers to the derivative of the 2nd component of the 4 valued function with respect to v and w. @METHOD : @GLOBALS : @CALLS : @CREATED : Jan 21, 1995 David MacDonald @MODIFIED : ---------------------------------------------------------------------------- */ VIOAPI void evaluate_interpolating_spline( int n_dims, Real parameters[], int degree, int n_values, Real coefs[], int n_derivs, Real derivs[] ) { int d, degrees[MAX_DIMS], n_derivs_list[MAX_DIMS]; Real *bases[MAX_DIMS]; if( degree < 1 || degree > 4 ) { print_error( "evaluate_interpolating_spline: invalid degree: %d\n", degree ); return; } if( n_dims < 1 || n_dims > MAX_DIMS ) { print_error( "evaluate_interpolating_spline: invalid n dims: %d\n", n_dims ); return; } switch( degree ) { case 1: bases[0] = &constant_coefs[0][0]; break; case 2: bases[0] = &linear_coefs[0][0]; break; case 3: bases[0] = &quadratic_coefs[0][0]; break; case 4: bases[0] = &cubic_coefs[0][0]; break; } for_less( d, 1, n_dims ) bases[d] = bases[0]; for_less( d, 0, n_dims ) { degrees[d] = degree; n_derivs_list[d] = n_derivs; } spline_tensor_product( n_dims, parameters, degrees, bases, n_values, coefs, n_derivs_list, derivs ); }