Index: /branches/eam_branches/psLib.20230123/src/math/psSpline.c =================================================================== --- /branches/eam_branches/psLib.20230123/src/math/psSpline.c (revision 42318) +++ /branches/eam_branches/psLib.20230123/src/math/psSpline.c (revision 42319) @@ -1,13 +1,4 @@ /** @file psSpline.c -* -* @brief Contains basic spline allocation, deallocation, fitting, -* and evaluation routines. -* -* This file contains the routines that allocate, free, and evaluate splines. -* -* @version $Revision: 1.158 $ $Name: not supported by cvs2svn $ -* @date $Date: 2007-03-14 02:36:28 $ -* -* Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii + re-written by EAM 2023.01.23 */ @@ -16,7 +7,4 @@ #endif -/*****************************************************************************/ -/* INCLUDE FILES */ -/*****************************************************************************/ #include #include @@ -36,29 +24,7 @@ #include "psMathUtils.h" -/*****************************************************************************/ -/* DEFINE STATEMENTS */ -/*****************************************************************************/ - -/*****************************************************************************/ -/* TYPE DEFINITIONS */ -/*****************************************************************************/ - -/*****************************************************************************/ -/* GLOBAL VARIABLES */ -/*****************************************************************************/ - -/*****************************************************************************/ -/* FILE STATIC VARIABLES */ -/*****************************************************************************/ - -/*****************************************************************************/ -/* FUNCTION IMPLEMENTATION - LOCAL */ -/*****************************************************************************/ - -static void spline1DFree(psSpline1D *tmpSpline) -{ - if (tmpSpline == NULL) { - return; - } +static void psSpline1DFree(psSpline1D *tmpSpline) +{ + if (tmpSpline == NULL) return; if (tmpSpline->spline != NULL) { @@ -78,15 +44,4 @@ } - -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) { @@ -94,33 +49,28 @@ 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); + // print the knots + } } /***************************************************************************** -CalculateSecondDerivs(): Given a set of x/y vectors corresponding to a -tabulated function at n points, this routine calculates the second -derivatives of the interpolating cubic splines at those n points. - -The first and second derivatives at the endpoints were previously undefined in -the SDR. From bugzilla #???, they are required to be 0.0, implementing natural -splines. - -Endpoints are defined by - PS_LEFT_SPLINE_DERIV - PS_RIGHT_SPLINE_DERIV - -This routine assumes that vectors x and y are of the appropriate types/sizes -(F32). - -XXX: use recycled vectors for internal data. -XXX: do an F64 version? +CalculateSecondDerivs(): Given a set of x,y vectors corresponding to the spline knots, this +routine calculates the second derivatives of the interpolating cubic splines at those n points. + +The boundary conditions may be: +* first derivative specified OR +* second derivatives = 0 + +One or the other condition may be true for each of the upper and lower bounds + +NOTE: vectors must be F32 + *****************************************************************************/ -#define PS_LEFT_SPLINE_DERIV 0.0 -#define PS_RIGHT_SPLINE_DERIV 0.0 + static psF32 *calculateSecondDerivs( - const psVector* x, ///< Ordinates - const psVector* y) ///< Coordinates + const psVector* x, ///< Ordinates + const psVector* y, ///< Coordinates + psF32 dyLower; // if not NAN, lower-bound 1st derivative is defined + psF32 dyUpper; // if not NAN, lower-bound 1st derivative is defined + ) { psTrace("psLib.math", 4, "---- %s() begin ----\n", __func__); @@ -190,5 +140,5 @@ /***************************************************************************** -psVectorFitSpline1D(): given a set of x/y vectors, this routine generates the +psSpline1DFitVector(): given a set of x/y vectors, this routine generates the linear or cublic splines which satisfy those data points. @@ -209,5 +159,5 @@ XXX: What types must be supported? *****************************************************************************/ -psSpline1D *psVectorFitSpline1D( +psSpline1D *psSpline1DFitVector( const psVector* x, ///< Ordinates. const psVector* y) ///< Coordinates. @@ -361,6 +311,6 @@ psF64? *****************************************************************************/ -float psSpline1DEval( - const psSpline1D *spline, +float psSpline1DEval_Old( + const psSpline1D_Old *spline, float x) { @@ -397,4 +347,41 @@ } +// XXX EAM : changing implementation to use yKnot, d2yKnot instead of polynomials +float psSpline1DEval_New( + const psSpline1D *spline, + float x) +{ + psTrace("psLib.math", 3, "---- %s() begin ----\n", __func__); + PS_ASSERT_PTR_NON_NULL(spline, NAN); + PS_ASSERT_INT_NONNEGATIVE(spline->n, NAN); + PS_ASSERT_VECTOR_TYPE(spline->knots, PS_TYPE_F32, NAN); + + psS32 n = spline->n; + if ((x < spline->knots->data.F32[0]) || (x > spline->knots->data.F32[spline->knots->n-1])) { + // If x is outside the range of spline->knots, generate a warning + // message, then return the left, or right, endpoint. + psLogMsg(__func__, PS_LOG_WARN, + "psSpline1DEval(): x ordinate (%f) is outside the spline range (%f - %f) (%d).", + x, spline->knots->data.F32[0], spline->knots->data.F32[n-1], n); + + psS32 binNum = (x < spline->knots->data.F32[0]) ? 0 : n-1; + psTrace("psLib.math", 3, "---- %s() end ----\n", __func__); + return(psPolynomial1DEval(spline->spline[binNum], x)); + } + + psScalar tmpScalar; + tmpScalar.type.type = PS_TYPE_F32; + tmpScalar.data.F32 = x; + psVectorBinaryDisectResult result; + psS32 binNum = psVectorBinaryDisect(&result, spline->knots, &tmpScalar); + if (result != PS_BINARY_DISECT_PASS) { + psError(PS_ERR_UNKNOWN, false, "Could not perform bin dissection on spline->knots.\n"); + return(NAN); + } + + psTrace("psLib.math", 3, "---- %s() end ----\n", __func__); + return(psPolynomial1DEval(spline->spline[binNum], x)); +} + /***************************************************************************** Index: /branches/eam_branches/psLib.20230123/src/math/psSpline.h =================================================================== --- /branches/eam_branches/psLib.20230123/src/math/psSpline.h (revision 42318) +++ /branches/eam_branches/psLib.20230123/src/math/psSpline.h (revision 42319) @@ -5,5 +5,6 @@ * and evaluate splines. * - * @author GLG, MHPCC + * @author GLG (MHPCC) + * reworked by EAM (IfA) * * @version $Revision: 1.62 $ $Name: not supported by cvs2svn $ @@ -27,20 +28,23 @@ #include "psPolynomial.h" -#define PS_LEFT_SPLINE_DERIV 0.0 -#define PS_RIGHT_SPLINE_DERIV 0.0 - -/** One-Dimensional Spline */ +/** One-Dimensional Spline + This structure represents a 1D cubic spline. Note the option for + equally-spaced knots allows a quick selection of the correct spline + segment. The values (xMin, xMax, xDel) are stored for ease of access. + */ typedef struct { unsigned int n; ///< The number of spline pieces - psPolynomial1D **spline; ///< An array of n pointers to the spline polynomials - psVector *knots; ///< The boundaries between each spline piece. Size is n+1. - psF32 *p_psDeriv2; ///< For cubic splines, the second derivative at each domain point. Size is n+1. + psF32 *xKnots; ///< x-coordinate of the knots (n+1) + psF32 *yKnots; ///< y-coordinate of the knots (n+1) + psF32 *d2yKnots; ///< 2nd derivative of y at the knots (n+1) + bool equalSpacing; // if knots are equally spaced, the seqment choice can be optimized + psF32 xMin; // for equally-spaced knots, the value at the lower bound (xKnots[0]) + psF32 xMax; // for equally-spaced knots, the value at the upper bound (xKnots[n]) + psF32 xDel; // for equally-spaced knots, the spacing (xKnots[n] - xKnots[0])/n } psSpline1D; -/** Allocates a psSpline1D structure - * - * Allocator for psSpline1D. +/** Allocator for psSpline1D. * * @return psSpline1D* new 1-D spline struct @@ -53,5 +57,5 @@ */ float psSpline1DEval( - const psSpline1D *spline, ///< Coefficients for spline polynomials + const psSpline1D *spline, ///< spline pointer float x ///< location at which to evaluate ); @@ -71,7 +75,9 @@ * generates the linear splines which satisfy those data points. * + * XXX EAM: add option to generate / select a subset of knots from the data + * * @return psSpline1D*: the calculated one-dimensional splines */ -psSpline1D *psVectorFitSpline1D( +psSpline1D *psSpline1DFitVector( const psVector* x, ///< Ordinates (or NULL to just use the indices) const psVector* y ///< Coordinates Index: /branches/eam_branches/psLib.20230123/test/math/tap_psSpline1D.c =================================================================== --- /branches/eam_branches/psLib.20230123/test/math/tap_psSpline1D.c (revision 42318) +++ /branches/eam_branches/psLib.20230123/test/math/tap_psSpline1D.c (revision 42319) @@ -1,10 +1,10 @@ -/* @file tst_psImageManip.c -* -* @brief This file will contain tests for all of the public psLib functions +/* @file tap_psImageManip.c +* +* @brief This file contains tests for all of the public psLib functions * that implement 1-D spline functionality: * psSpline1DAlloc() * psSpline1DEval() * psSpline1DEvalVector() -* psVectorFitSpline1D() +* psSpline1DFitVector() * * This file is composed of the tests formerly in tst_psFunc02.c, @@ -35,4 +35,6 @@ psF32 expect) { + // NOTE: this returns NAN if 'expect' == 0.0 + // NOTE 2: this is not testing a percent, but fractional error if ((fabs(actual - expect) / fabs(expect)) > ERROR_TOLERANCE_PERCENT) { return(true); @@ -83,9 +85,16 @@ typedef psF64 (*mappingFuncF64)(psF64 x); +/* EAM 2023.01.22 : these tests are not well considered. They generate a set of N+1 points + to use as knots at integer values from 0 to N, then evaluate the spline half-way between + the knots. the function above is a cubic, so a cubic spline should fit it perfectly, which + is fine. Perhaps this test suite should use a pre-defined collection of data points? + */ + bool genericF32Test(psS32 NumSplines, mappingFuncF32 func, bool xNull) { - // We test the psVectorFitSpline1D, psSpline1DEval() functions. F32 version. + // We test the psSpline1DFitVector, psSpline1DEval() functions. F32 version. bool testStatus = true; { + // Generate the vector data psMemId id = psMemGetId(); psVector *xF32 = psVectorAlloc(NumSplines+1, PS_TYPE_F32); @@ -98,32 +107,33 @@ psSpline1D *tmpSpline = NULL; if (!xNull) { - tmpSpline = psVectorFitSpline1D(xF32, yF32); + tmpSpline = psSpline1DFitVector(xF32, yF32); } else { - tmpSpline = psVectorFitSpline1D(NULL, yF32); - } + tmpSpline = psSpline1DFitVector(NULL, yF32); + } + if(tmpSpline == NULL) { - diag("psVectorFitSpline1D() returned NULL"); + diag("psSpline1DFitVector() returned NULL"); testStatus = false; } else { if (tmpSpline->n != NumSplines) { - diag("psVectorFitSpline1D() did not properly set the psSpline1D->n member"); + diag("psSpline1DFitVector() did not properly set the psSpline1D->n member"); testStatus = false; } if(tmpSpline->spline == NULL) { - diag("psVectorFitSpline1D() returned a NULL psSpline1D->spline member."); + diag("psSpline1DFitVector() returned a NULL psSpline1D->spline member."); testStatus = false; } for (psS32 i = 0 ; i < NumSplines ; i++) { if (tmpSpline->spline[i] == NULL) { - diag("psVectorFitSpline1D() returned a NULL psSpline1D->spline[%d] member.", i); + diag("psSpline1DFitVector() returned a NULL psSpline1D->spline[%d] member.", i); testStatus = false; } } if (tmpSpline->knots == NULL) { - diag("psVectorFitSpline1D() returned a NULL psSpline1D->knots member"); + diag("psSpline1DFitVector() returned a NULL psSpline1D->knots member"); testStatus = false; } if (tmpSpline->p_psDeriv2 == NULL) { - diag("psVectorFitSpline1D()returned a NULL psSpline1D->p_psDeriv2 member"); + diag("psSpline1DFitVector()returned a NULL psSpline1D->p_psDeriv2 member"); testStatus = false; } @@ -136,4 +146,5 @@ testStatus = false; diag("TEST ERROR: f(%f) is %f, should be %f", x, y, myFunc00(x)); + // XXX EAM : the truth value above should be 'func' not myFunc00 } } @@ -150,4 +161,5 @@ diag("TEST ERROR: f(%f) is %f, should be %f", xF32->data.F32[i], yF32Test->data.F32[i], myFunc00(xF32->data.F32[i])); + // XXX EAM : the truth value above should be 'func' not myFunc00 } } @@ -171,5 +183,5 @@ bool genericF64Test(psS32 NumSplines, mappingFuncF64 func, bool xNull) { - // We test the psVectorFitSpline1D, psSpline1DEval() functions. F64 version. + // We test the psSpline1DFitVector, psSpline1DEval() functions. F64 version. bool testStatus = true; { @@ -184,32 +196,32 @@ psSpline1D *tmpSpline = NULL; if (!xNull) { - tmpSpline = psVectorFitSpline1D(xF64, yF64); + tmpSpline = psSpline1DFitVector(xF64, yF64); } else { - tmpSpline = psVectorFitSpline1D(NULL, yF64); + tmpSpline = psSpline1DFitVector(NULL, yF64); } if(tmpSpline == NULL) { - diag("psVectorFitSpline1D() returned NULL"); + diag("psSpline1DFitVector() returned NULL"); testStatus = false; } else { if (tmpSpline->n != NumSplines) { - diag("psVectorFitSpline1D() did not properly set the psSpline1D->n member"); + diag("psSpline1DFitVector() did not properly set the psSpline1D->n member"); testStatus = false; } if(tmpSpline->spline == NULL) { - diag("psVectorFitSpline1D() returned a NULL psSpline1D->spline member."); + diag("psSpline1DFitVector() returned a NULL psSpline1D->spline member."); testStatus = false; } for (psS32 i = 0 ; i < NumSplines ; i++) { if (tmpSpline->spline[i] == NULL) { - diag("psVectorFitSpline1D() returned a NULL psSpline1D->spline[%d] member.", i); + diag("psSpline1DFitVector() returned a NULL psSpline1D->spline[%d] member.", i); testStatus = false; } } if (tmpSpline->knots == NULL) { - diag("psVectorFitSpline1D() returned a NULL psSpline1D->knots member"); + diag("psSpline1DFitVector() returned a NULL psSpline1D->knots member"); testStatus = false; } if (tmpSpline->p_psDeriv2 == NULL) { - diag("psVectorFitSpline1D()returned a NULL psSpline1D->p_psDeriv2 member"); + diag("psSpline1DFitVector()returned a NULL psSpline1D->p_psDeriv2 member"); testStatus = false; } @@ -274,9 +286,9 @@ } - // psSplineEvalTest_sub(): Call psVectorFitSpline1D with NULL arguments. - { - psMemId id = psMemGetId(); - psSpline1D *tmpSpline = psVectorFitSpline1D(NULL, NULL); - ok(tmpSpline == NULL, "psVectorFitSpline1D() returns NULL with NULL arguments"); + // psSplineEvalTest_sub(): Call psSpline1DFitVector with NULL arguments. + { + psMemId id = psMemGetId(); + psSpline1D *tmpSpline = psSpline1DFitVector(NULL, NULL); + ok(tmpSpline == NULL, "psSpline1DFitVector() returns NULL with NULL arguments"); ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks"); } @@ -287,5 +299,5 @@ psMemId id = psMemGetId(); float y = psSpline1DEval(NULL, 0.0); - ok(isnan(y), "psSpline1DEval() returned NAN will NULL input spline"); + ok(isnan(y), "psSpline1DEval() returned NAN with NULL input spline"); ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks"); } @@ -297,5 +309,5 @@ psVector *x = psVectorAlloc(10, PS_TYPE_F32); psVector *y = psSpline1DEvalVector(NULL, x); - ok(y == NULL, "psSpline1DEvalVector() returned NAN will NULL input spline"); + ok(y == NULL, "psSpline1DEvalVector() returned NULL with NULL input spline"); psFree(x); psFree(y); @@ -316,5 +328,4 @@ ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks"); } - ok(genericF32Test(1, myFunc00, false), "Generic, simple mapping, F32 test. 1 spline");