Index: trunk/psLib/src/math/psSpline.c =================================================================== --- trunk/psLib/src/math/psSpline.c (revision 6204) +++ trunk/psLib/src/math/psSpline.c (revision 6305) @@ -4,9 +4,8 @@ * and evaluation routines. * -* This file will hold the functions for allocated, freeing, and evaluating -* splines. +* This file contains the routines that allocate, free, and evaluate splines. * -* @version $Revision: 1.134 $ $Name: not supported by cvs2svn $ -* @date $Date: 2006-01-26 21:10:22 $ +* @version $Revision: 1.135 $ $Name: not supported by cvs2svn $ +* @date $Date: 2006-02-02 21:09:08 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -31,4 +30,5 @@ #include "psConstants.h" #include "psErrorText.h" +#include "psMathUtils.h" /*****************************************************************************/ @@ -72,4 +72,25 @@ return; +} + + +static void PS_POLY1D_PRINT( + psPolynomial1D *poly) +{ + printf("-------------- PS_POLY1D_PRINT() --------------\n"); + printf("poly->nX is %d\n", poly->nX); + for (psS32 i = 0 ; i < (1 + poly->nX) ; i++) { + printf("poly->coeff[%d] is %f\n", i, poly->coeff[i]); + } +} + +static void PS_PRINT_SPLINE2(psSpline1D *mySpline) +{ + printf("-------------- PS_PRINT_SPLINE2() --------------\n"); + printf("mySpline->n is %d\n", mySpline->n); + for (psS32 i = 0 ; i < mySpline->n ; i++) { + PS_POLY1D_PRINT(mySpline->spline[i]); + } + PS_VECTOR_PRINT_F32(mySpline->knots); } @@ -139,310 +160,4 @@ psTrace(__func__, 4, "---- %s() end ----\n", __func__); return(derivs2); -} - - -/***************************************************************************** -vectorBinDisectTYPE(): This is a macro for a private function which takes as -input a vector an array of data as well as a single value for that data. The -input vector values are assumed to be non-decreasing (v[i-1] <= v[i] for all -i). This routine does a binary disection of the vector and returns "i" such -that (v[i] <= x <= v[i+1). If x lies outside the range of v[], then this -routine prints a warning message and returns (-2 or -1). - *****************************************************************************/ -#define FUNC_MACRO_VECTOR_BIN_DISECT(TYPE) \ -static psS32 vectorBinDisect##TYPE(ps##TYPE *bins, \ - psS32 numBins, \ - ps##TYPE x) \ -{ \ - psS32 min; \ - psS32 max; \ - psS32 mid; \ - psTrace(__func__, 4, "---- %s() begin ----\n", __func__); \ - if (x < bins[0]) { \ - psLogMsg(__func__, PS_LOG_WARN, \ - "vectorBinDisect%s(): ordinate %f is outside vector range (%f - %f).", \ - #TYPE, x, bins[0], bins[numBins-1]); \ - return(-2); \ - } \ - if (x > bins[numBins-1]) { \ - psLogMsg(__func__, PS_LOG_WARN, \ - "vectorBinDisect%s(): ordinate %f is outside vector range (%f - %f).", \ - #TYPE, x, bins[0], bins[numBins-1]); \ - return(-1); \ - } \ - \ - min = 0; \ - max = numBins-2; \ - mid = ((max+1)-min)/2; \ - while (min != max) { \ - psTrace(__func__, 6, "(min, mid, max) is (%u, %u, %u): (x, bins) is (%f, %f)\n", min, mid, max, x, bins[mid]); \ - \ - if (x == bins[mid]) { \ - psTrace(__func__, 4, "---- %s(%d) end ----\n", __func__, mid); \ - return(mid); \ - } else if (x < bins[mid]) { \ - max = mid-1; \ - } else { \ - min = mid; \ - } \ - mid = ((max+1)+min)/2; \ - } \ - psTrace(__func__, 4, "---- %s(%d) end ----\n", __func__, min); \ - return(min); \ -} \ - -FUNC_MACRO_VECTOR_BIN_DISECT(S8) -FUNC_MACRO_VECTOR_BIN_DISECT(S16) -FUNC_MACRO_VECTOR_BIN_DISECT(S32) -FUNC_MACRO_VECTOR_BIN_DISECT(S64) -FUNC_MACRO_VECTOR_BIN_DISECT(U8) -FUNC_MACRO_VECTOR_BIN_DISECT(U16) -FUNC_MACRO_VECTOR_BIN_DISECT(U32) -FUNC_MACRO_VECTOR_BIN_DISECT(U64) -FUNC_MACRO_VECTOR_BIN_DISECT(F32) -FUNC_MACRO_VECTOR_BIN_DISECT(F64) - -/***************************************************************************** -p_psVectorBinDisect(): A wrapper to the above p_psVectorBinDisect(). - *****************************************************************************/ -psS32 p_psVectorBinDisect( - psVector *bins, - psScalar *x) -{ - PS_ASSERT_VECTOR_NON_NULL(bins, -4); - PS_ASSERT_VECTOR_NON_EMPTY(bins, -4); - PS_ASSERT_PTR_NON_NULL(x, -6); - PS_ASSERT_PTR_TYPE_EQUAL(x, bins, -3); - char* strType; - - switch (x->type.type) { - case PS_TYPE_U8: - return(vectorBinDisectU8(bins->data.U8, bins->n, x->data.U8)); - case PS_TYPE_U16: - return(vectorBinDisectU16(bins->data.U16, bins->n, x->data.U16)); - case PS_TYPE_U32: - return(vectorBinDisectU32(bins->data.U32, bins->n, x->data.U32)); - case PS_TYPE_U64: - return(vectorBinDisectU64(bins->data.U64, bins->n, x->data.U64)); - case PS_TYPE_S8: - return(vectorBinDisectS8(bins->data.S8, bins->n, x->data.S8)); - case PS_TYPE_S16: - return(vectorBinDisectS16(bins->data.S16, bins->n, x->data.S16)); - case PS_TYPE_S32: - return(vectorBinDisectS32(bins->data.S32, bins->n, x->data.S32)); - case PS_TYPE_S64: - return(vectorBinDisectS64(bins->data.S64, bins->n, x->data.S64)); - case PS_TYPE_F32: - return(vectorBinDisectF32(bins->data.F32, bins->n, x->data.F32)); - case PS_TYPE_F64: - return(vectorBinDisectF64(bins->data.F64, bins->n, x->data.F64)); - case PS_TYPE_C32: - PS_TYPE_NAME(strType,x->type.type); - psError(PS_ERR_BAD_PARAMETER_TYPE, - PS_ERRORTEXT_psSpline_TYPE_NOT_SUPPORTED, - strType); - return 0; - case PS_TYPE_C64: - PS_TYPE_NAME(strType,x->type.type); - psError(PS_ERR_BAD_PARAMETER_TYPE, - PS_ERRORTEXT_psSpline_TYPE_NOT_SUPPORTED, - strType); - return 0; - case PS_TYPE_BOOL: - PS_TYPE_NAME(strType,x->type.type); - psError(PS_ERR_BAD_PARAMETER_TYPE, - PS_ERRORTEXT_psSpline_TYPE_NOT_SUPPORTED, - strType); - return 0; - } - return(-3); -} - -/***************************************************************************** -fullInterpolate1DF32(): This routine will take as input n-element floating -point arrays domain and range, and the x value, assumed to lie with the -domain vector. It produces as output the (n-1)-order LaGrange interpolated -value of x. - -XXX: do we error check for non-distinct domain values? - -XXX: This is a somewhat general function. There's no reason to put it in this -file. - -XXX: Sort all these interpolation functions out. Why are there so many? - -XXX: The fullInterpolate1D macros are not used anywhere. - *****************************************************************************/ -#define FUNC_MACRO_FULL_INTERPOLATE_1D(TYPE) \ -static psF32 fullInterpolate1D##TYPE(ps##TYPE *domain, \ - ps##TYPE *range, \ - unsigned int n, \ - ps##TYPE x) \ -{ \ - \ - unsigned int i; \ - unsigned int m; \ - static psVector *p = NULL; \ - p = psVectorRecycle(p, n, PS_TYPE_##TYPE); \ - p_psMemSetPersistent(p, true); \ - p_psMemSetPersistent(p->data.TYPE, true); \ - \ - psTrace(__func__, 4, "---- %s() begin %u-order at x=%f) (%d data points) ----\n", __func__, n-1, x, n); \ - for (i=0;idata.TYPE[i] = range[i]; \ - psTrace(__func__, 6, "p->data.TYPE[%u] is %f\n", i, p->data.TYPE[i]); \ - \ - } \ - \ - /* From NR, during each iteration of the m loop, we are computing the \ - p_{i ... i+m} terms. \ - */ \ - for (m=1;mdata.TYPE[i] = (((x-domain[i+m]) * p->data.TYPE[i]) + \ - ((domain[i]-x) * p->data.TYPE[i+1])) / \ - (domain[i] - domain[i+m]); \ - psTrace(__func__, 6, "p->data.TYPE[%u] is %f\n", i, p->data.TYPE[i]); \ - } \ - } \ - psTrace(__func__, 4, "---- %s(....) end ----\n", __func__); \ - return(p->data.TYPE[0]); \ -} \ - -/* XXX: Do this correctly. -FUNC_MACRO_FULL_INTERPOLATE_1D(U8) -FUNC_MACRO_FULL_INTERPOLATE_1D(U16) -FUNC_MACRO_FULL_INTERPOLATE_1D(U32) -FUNC_MACRO_FULL_INTERPOLATE_1D(U64) -FUNC_MACRO_FULL_INTERPOLATE_1D(S8) -FUNC_MACRO_FULL_INTERPOLATE_1D(S16) -FUNC_MACRO_FULL_INTERPOLATE_1D(S32) -FUNC_MACRO_FULL_INTERPOLATE_1D(S64) -FUNC_MACRO_FULL_INTERPOLATE_1D(F64) -*/ -FUNC_MACRO_FULL_INTERPOLATE_1D(F32) - - -/***************************************************************************** -interpolate1DF32(): this is the base 1-D flat memory routine to perform -LaGrange interpolation. - *****************************************************************************/ -static psF32 interpolate1DF32( - psF32 *domain, - psF32 *range, - unsigned int n, - unsigned int order, - psF32 x) -{ - psTrace(__func__, 4, "---- %s() begin ----\n", __func__); - PS_ASSERT_PTR_NON_NULL(domain, NAN) - PS_ASSERT_PTR_NON_NULL(range, NAN) - // XXX: Check valid values for n, order, and x? - - psS32 binNum; - psS32 numIntPoints = order+1; - psS32 origin; - - binNum = vectorBinDisectF32(domain, n, x); - - if (0 == numIntPoints%2) { - origin = binNum - ((numIntPoints/2) - 1); - } else { - origin = binNum - (numIntPoints/2); - if ((x-domain[binNum]) > (domain[binNum+1]-x)) { - // x is closer to binNum+1. - origin = 1 + (binNum - (numIntPoints/2)); - } - } - if (origin < 0) { - origin = 0; - } - if ((origin + numIntPoints) > n) { - origin = n - numIntPoints; - } - - psTrace(__func__, 4, "---- %s(....) end ----\n", __func__); - return(fullInterpolate1DF32(&domain[origin], &range[origin], order+1, x)); -} - - -/***************************************************************************** -p_psVectorInterpolate(): This routine will take as input psVectors domain and -range, and the x value, assumed to lie with the domain vector. It produces -as output the LaGrange interpolated value of a polynomial of the specified -order around the point x. - -XXX: This stuff does not currently work with a mask. - -XXX: add another psScalar argument for the result (so you don't have to - allocate it each time). - -XXX: The VectorCopy routines seg fault when I declare range32 as static. - *****************************************************************************/ -psScalar *p_psVectorInterpolate( - psVector *domain, - psVector *range, - psS32 order, - psScalar *x) -{ - psTrace(__func__, 4, "---- %s() begin ----\n", __func__); - PS_ASSERT_VECTOR_NON_NULL(domain, NULL); - PS_ASSERT_VECTOR_NON_NULL(range, NULL); - PS_ASSERT_PTR_NON_NULL(x, NULL); - PS_ASSERT_INT_NONNEGATIVE(order, NULL); - PS_ASSERT_VECTORS_SIZE_EQUAL(domain, range, NULL); - PS_ASSERT_PTR_TYPE_EQUAL(domain, range, NULL); - PS_ASSERT_PTR_TYPE_EQUAL(domain, x, NULL); - - psVector *range32 = NULL; - psVector *domain32 = NULL; - if (order > (domain->n - 1)) { - psError(PS_ERR_BAD_PARAMETER_SIZE, true, - PS_ERRORTEXT_psSpline_NOT_ENOUGH_DATAPOINTS, - order); - return(NULL); - } - - if (x->type.type == PS_TYPE_F32) { - psTrace(__func__, 4, "---- %s(NULL) end ----\n", __func__); - return(psScalarAlloc(interpolate1DF32(domain->data.F32, - range->data.F32, - domain->n, - order, - x->data.F32), PS_TYPE_F32)); - } else if (x->type.type == PS_TYPE_F64) { - // XXX: use recycled vectors here. - range32 = psVectorCopy(range32, range, PS_TYPE_F32); - domain32 = psVectorCopy(domain32, domain, PS_TYPE_F32); - - psScalar *tmpScalar = psScalarAlloc((psF64) - interpolate1DF32(domain32->data.F32, - range32->data.F32, - domain32->n, - order, - (psF32) x->data.F64), PS_TYPE_F64); - psFree(range32); - psFree(domain32); - - psTrace(__func__, 4, "---- %s() end ----\n", __func__); - // XXX: Convert data type to F64? - return(tmpScalar); - - } else { - char* strType; - PS_TYPE_NAME(strType,x->type.type); - psError(PS_ERR_BAD_PARAMETER_TYPE, - PS_ERRORTEXT_psSpline_TYPE_NOT_SUPPORTED, - strType); - } - - psTrace(__func__, 4, "---- %s(NULL) end ----\n", __func__); - return(NULL); } @@ -627,59 +342,4 @@ } -void PS_POLY1D_PRINT( - psPolynomial1D *poly) -{ - printf("-------------- PS_POLY1D_PRINT() --------------\n"); - printf("poly->nX is %d\n", poly->nX); - for (psS32 i = 0 ; i < (1 + poly->nX) ; i++) { - printf("poly->coeff[%d] is %f\n", i, poly->coeff[i]); - } -} - -void PS_PRINT_SPLINE2(psSpline1D *mySpline) -{ - printf("-------------- PS_PRINT_SPLINE2() --------------\n"); - printf("mySpline->n is %d\n", mySpline->n); - for (psS32 i = 0 ; i < mySpline->n ; i++) { - PS_POLY1D_PRINT(mySpline->spline[i]); - } - PS_VECTOR_PRINT_F32(mySpline->knots); -} - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - /***************************************************************************** @@ -706,5 +366,9 @@ psS32 n = spline->n; - psS32 binNum = vectorBinDisectF32(spline->knots->data.F32, (spline->n)+1, x); + psScalar tmpScalar; + tmpScalar.type.type = PS_TYPE_F32; + tmpScalar.data.F32 = x; + psS32 binNum = p_psVectorBinDisect(spline->knots, &tmpScalar); + if (binNum < 0) { // @@ -746,5 +410,5 @@ if (psTraceGetLevel(__func__) >= 6) { PS_VECTOR_PRINT_F32(x); - // XXX: print spline + PS_PRINT_SPLINE2((psSpline1D *) spline); }