Mercurial > hg > minc-tools
changeset 769:df6fbc2fc74e
Initial revision
author | david <david> |
---|---|
date | Tue, 21 Feb 1995 12:58:33 +0000 |
parents | 775ad72771e6 |
children | 5bf93c6b3d1d |
files | volume_io/Geometry/splines.c volume_io/Geometry/tensors.c |
diffstat | 2 files changed, 294 insertions(+), 0 deletions(-) [+] |
line wrap: on
line diff
new file mode 100644 --- /dev/null +++ b/volume_io/Geometry/splines.c @@ -0,0 +1,159 @@ +#include <internal_volume_io.h> + +private Real linear_coefs[2][2] = { + { 1.0, 0.0 }, + { -1.0, 1.0 } + }; +private Real quadratic_coefs[3][3] = { + { 0.5, 0.5, 0.0 }, + { -1.0, 1.0, 0.0 }, + { 0.5, -1.0, 0.5 } + }; + +private 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 } + }; + +public 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]; +} + +public 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]; +} + +public 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]; +} + +public 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 ); +} + +public 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 ); +} + +public 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 ); +} + +public 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 + +public 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 < 2 || degree > 4 ) + { + print( "evaluate_interpolating_spline: invalid degree: %d\n", degree ); + return; + } + + if( n_dims < 1 || n_dims > MAX_DIMS ) + { + print( "evaluate_interpolating_spline: invalid n dims: %d\n", n_dims ); + return; + } + + switch( degree ) + { + 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 ); +}
new file mode 100755 --- /dev/null +++ b/volume_io/Geometry/tensors.c @@ -0,0 +1,135 @@ +#include <internal_volume_io.h> + +private void multiply_matrices( + int x1, + int y1, + Real m1[], + int sa1, + int sb1, + int y2, + Real m2[], + int sa2, + int sb2, + Real prod[], + int sap, + int sbp ) +{ + int i, j, k; + Real sum, *m1_ptr, *m2_ptr, *prod_ptr; + + for_less( i, 0, x1 ) + { + prod_ptr = &prod[i*sap]; + for_less( j, 0, y2 ) + { + sum = 0.0; + m1_ptr = &m1[i*sa1]; + m2_ptr = &m2[j*sb2]; + for_less( k, 0, y1 ) + { + sum += (*m1_ptr) * (*m2_ptr); + m1_ptr += sb1; + m2_ptr += sa2; + } + *prod_ptr = sum; + prod_ptr += sbp; + } + } +} + +#define MAX_DEGREE 4 +#define MAX_TOTAL_VALUES 4000 + +public void spline_tensor_product( + int n_dims, + Real positions[], + int degrees[], /* [n_dims] */ + Real *bases[], /* [n_dims][degress[dim]*degrees[dim]] */ + int n_values, + Real coefs[], /* [n_values*degrees[0]*degrees[1]*...] */ + int n_derivs[], /* [n_dims] */ + Real results[] ) /* [n_values*n_derivs[0]*n_derivs[1]*...] */ +{ + int i, deriv, d, k, total_values, src; + int ind, prev_ind; + Real us[MAX_DEGREE*MAX_DEGREE], weights[MAX_DEGREE*MAX_DEGREE]; + Real *tmp_results[2], *r; + Real static_tmp_results[2*MAX_TOTAL_VALUES]; + + /*--- check arguments */ + + total_values = n_values; + for_less( d, 0, n_dims ) + { + if( degrees[d] < 2 || degrees[d] > MAX_DEGREE ) + { + print( + "spline_tensor_product: Degree %d is not be between 2 and %d\n", + degrees[d], MAX_DEGREE ); + return; + } + total_values *= degrees[d]; + } + + if( total_values > MAX_TOTAL_VALUES ) + { + print( "Spline size too large for static memory.\n" ); + return; + } + + tmp_results[0] = &static_tmp_results[0]; + tmp_results[1] = &static_tmp_results[total_values]; + + for_less( i, 0, total_values ) + tmp_results[0][i] = coefs[i]; + + src = 0; + + /*--- do each dimension */ + + for_less( d, 0, n_dims ) + { + us[0] = 1.0; + for_less( k, 1, degrees[d] ) + us[k] = us[k-1] * positions[d]; + + ind = degrees[d]; + for_inclusive( deriv, 1, n_derivs[d] ) + { + for_less( k, 0, deriv ) + { + us[ind] = 0.0; + ++ind; + } + + prev_ind = IJ( deriv-1, deriv-1, degrees[d] ); + for_less( k, deriv, degrees[d] ) + { + us[ind] = us[prev_ind] * (Real) k; + ++ind; + ++prev_ind; + } + } + + multiply_matrices( 1 + n_derivs[d], degrees[d], us, degrees[d], 1, + degrees[d], bases[d], degrees[d], 1, + weights, degrees[d], 1 ); + + total_values = n_values; + for_less( i, 0, d ) + total_values *= 1 + n_derivs[i]; + for_less( i, d+1, n_dims ) + total_values *= degrees[i]; + + if( d == n_dims-1 ) + r = results; + else + r = tmp_results[1-src]; + + multiply_matrices( 1 + n_derivs[d], degrees[d], weights, degrees[d], 1, + total_values, tmp_results[src], total_values, 1, + r, 1, 1 + n_derivs[d] ); + + src = 1 - src; + } +}