Changeset 42336

Timestamp:

Jan 30, 2023, 9:51:20 AM (3 years ago)

Author:

eugene

Message:

merge from branches/eam_branches/psLib.20230123: rework of psSpline (cleaner names, do not convert to generic polynomials, allow defined 1st deriv or 0.0 2nd deriv for boundary conditions)

Location:

trunk/psLib

Files:

• . (modified) (1 prop)
• src/math/psPolynomial.h (modified) (2 diffs)
• src/math/psSpline.c (modified) (9 diffs)
• src/math/psSpline.h (modified) (5 diffs)
• test/math/mana.spline.pro (copied) (copied from branches/eam_branches/psLib.20230123/test/math/mana.spline.pro )
• test/math/tap_psSpline1D.c (modified) (12 diffs)

trunk/psLib/src/math/psPolynomial.h

@@ -353 +353 @@
 } \
 
+// XXX warning: this is fragile if NAME contains an external 'i'
 #define PS_POLY_PRINT_1D(NAME) \
 printf("Poly %s: (nX) is (%d)\n", #NAME, NAME->nX);\
@@ -359 +360 @@
 }\
 
+// XXX warning: this is fragile if NAME contains an external 'i' or 'j'
 #define PS_POLY_PRINT_2D(NAME) \
 printf("Poly %s: (nX, nY) is (%d, %d)\n", #NAME, NAME->nX, NAME->nY);\

trunk/psLib/src/math/psSpline.c

@@ -1 +1 @@
 /** @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 @@
 #endif
 
-/*****************************************************************************/
-/*  INCLUDE FILES                                                            */
-/*****************************************************************************/
 #include <stdio.h>
 #include <stdbool.h>
@@ -36 +24 @@
 #include "psMathUtils.h"
 
-/*****************************************************************************/
-/* DEFINE STATEMENTS                                                         */
-/*****************************************************************************/
-
-/*****************************************************************************/
-/* TYPE DEFINITIONS                                                          */
-/*****************************************************************************/
-
-/*****************************************************************************/
-/* GLOBAL VARIABLES                                                          */
-/*****************************************************************************/
-
-/*****************************************************************************/
-/* FILE STATIC VARIABLES                                                     */
-/*****************************************************************************/
-
-/*****************************************************************************/
-/* FUNCTION IMPLEMENTATION - LOCAL                                           */
-/*****************************************************************************/
-
-static void spline1DFree(psSpline1D *tmpSpline)
-{
-    if (tmpSpline == NULL) {
-        return;
-    }
-
-    if (tmpSpline->spline != NULL) {
-        for (psS32 i=0;i<tmpSpline->n;i++) {
-            psFree((tmpSpline->spline)[i]);
-        }
-        psFree(tmpSpline->spline);
-    }
-
-    if (tmpSpline->p_psDeriv2 != NULL) {
-        psFree(tmpSpline->p_psDeriv2);
-    }
-
-    psFree(tmpSpline->knots);
+static void psSpline1DFree(psSpline1D *tmpSpline)
+{
+    if (tmpSpline == NULL) return;
+
+    psFree(tmpSpline->xKnots);
+    psFree(tmpSpline->yKnots);
+    psFree(tmpSpline->d2yKnots);
 
     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);
+    if (!mySpline->xKnots) return;
+    if (!mySpline->yKnots) return;
+    if (!mySpline->d2yKnots) return;
     for (psS32 i = 0 ; i < mySpline->n ; i++) {
-        PS_POLY1D_PRINT(mySpline->spline[i]);
-    }
-    PS_VECTOR_PRINT_F32(mySpline->knots);
-}
+        printf("(x, y, d2y) : %f %f %f\n", mySpline->xKnots[i], mySpline->yKnots[i], mySpline->d2yKnots[i]);
+    }
+}
+*/
 
 /*****************************************************************************
-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
+
+EAM 2023.01.19 : the comment above is wrong: the code below implements the 
+splines with *only* the 1st derivatives at the end points set to 0.0.  the 
+2nd derivatives are constrained by the equations for the splines.  it is not 
+possible to specify both the endpoint 1st and 2nd derivaties: there would be 
+too many constraints for the number of free parameters.
+
+It is not clear that choosing to set the end point 1st derivatives to 0.0 is the best
+option.  setting the 2nd derivatives to zero allow for linear extrapolation.
+
  *****************************************************************************/
-#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
+
+void calculateSecondDerivs(
+    const psSpline1D *mySpline,
+    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__);
-    if (psTraceGetLevel("psLib.math") >= 6) {
-        p_psVectorPrint(1, (psVector *) x, "x");
-        p_psVectorPrint(1, (psVector *) y, "y");
-    }
-    psS32 n = y->n;
-    psF32 *u = (psF32 *) psAlloc(n * sizeof(psF32));
-    psF32 *derivs2 = (psF32 *) psAlloc(n * sizeof(psF32));
-    psF32 *X = (psF32 *) & (x->data.F32[0]);
-    psF32 *Y = (psF32 *) & (y->data.F32[0]);
-    //
-    // The second derivatives at the endpoints, undefined in the SDR,
-    // are set in psAssert.h: PS_LEFT_SPLINE_DERIV, PS_RIGHT_SPLINE_DERIV.
-    //
-    derivs2[0] = -0.5;
-    u[0]= (3.0/(X[1]-X[0])) * ((Y[1]-Y[0])/(X[1]-X[0]) - PS_LEFT_SPLINE_DERIV);
-
-    for (psS32 i=1;i<=(n-2);i++) {
-        psF32 sig = (X[i] - X[i-1]) / (X[i+1] - X[i-1]);
-        psF32 p = sig * derivs2[i-1] + 2.0;
-        derivs2[i] = (sig - 1.0) / p;
-        u[i] = ((Y[i+1] - Y[i])/(X[i+1]-X[i])) - ((Y[i]-Y[i-1])/(X[i]-X[i-1]));
-        u[i] = ((6.0 * u[i] / (X[i+1] - X[i-1])) - (sig * u[i-1])) / p;
+
+    psS32 n = mySpline->n; // n is the number of knots
+    psF32 *u   = (psF32 *) psAlloc(n * sizeof(psF32));
+
+    psF32 *X   = mySpline->xKnots;
+    psF32 *Y   = mySpline->yKnots;
+    psF32 *d2y = mySpline->d2yKnots;
+
+    if (isfinite(dyLower)) {
+      d2y[0] = -0.5;
+      u[0]   = (3.0/(X[1]-X[0])) * ((Y[1]-Y[0])/(X[1]-X[0]) - dyLower);
+    } else {
+      d2y[0] = 0.0;
+      u[0]   = 0.0;
+    }
+
+    for (psS32 i = 1; i < n - 1; i++) {
+        psF32 dX = (X[i] - X[i-1]) / (X[i+1] - X[i-1]);
+        psF32 dY = dX * d2y[i-1] + 2.0;
+        d2y[i] = (dX - 1.0) / dY;
+        u[i]   = ((Y[i+1] - Y[i])/(X[i+1] - X[i])) - ((Y[i] - Y[i-1])/(X[i] - X[i-1]));
+        u[i]   = ((6.0 * u[i] / (X[i+1] - X[i-1])) - (dX * u[i-1])) / dY;
 
         psTrace("psLib.math", 6, "X[%d] is %f\n", i, X[i]);
@@ -153 +107 @@
     }
 
-    psF32 qn = 0.5;
-    u[n-1] = (3.0/(X[n-1]-X[n-2])) * (PS_RIGHT_SPLINE_DERIV - (Y[n-1]-Y[n-2])/(X[n-1]-X[n-2]));
-    derivs2[n-1] = (u[n-1] - (qn * u[n-2])) / ((qn * derivs2[n-2]) + 1.0);
-
-    for (psS32 k=(n-2);k>=0;k--) {
-        derivs2[k] = derivs2[k] * derivs2[k+1] + u[k];
-        psTrace("psLib.math", 6, "derivs2[%d] is %f\n", k, derivs2[k]);
+    if (isfinite(dyUpper)) {
+      psF32 qn = 0.5;
+      u[n-1] = (3.0/(X[n-1]-X[n-2])) * (dyUpper - (Y[n-1]-Y[n-2])/(X[n-1]-X[n-2]));
+      d2y[n-1] = (u[n-1] - (qn * u[n-2])) / ((qn * d2y[n-2]) + 1.0);
+    } else {
+      d2y[n-1] = 0;
+    }
+
+    for (psS32 k = n-2; k >= 0; k--) {
+        d2y[k] = d2y[k] * d2y[k+1] + u[k];
+        psTrace("psLib.math", 6, "derivs2[%d] is %f\n", k, d2y[k]);
     }
     psFree(u);
     psTrace("psLib.math", 4, "---- %s() end ----\n", __func__);
-    return(derivs2);
+    return;
 }
 
@@ -174 +132 @@
 {
     PS_ASSERT_PTR(ptr, false);
-    return ( psMemGetDeallocator(ptr) == (psFreeFunc)spline1DFree );
+    return ( psMemGetDeallocator(ptr) == (psFreeFunc) psSpline1DFree );
 }
 
@@ -180 +138 @@
 {
     psSpline1D *tmpSpline = (psSpline1D *) psAlloc(sizeof(psSpline1D));
+
     tmpSpline->n = 0;
-    tmpSpline->spline = NULL;
-    tmpSpline->knots = NULL;
-    tmpSpline->p_psDeriv2 = NULL;
-    psMemSetDeallocator(tmpSpline, (psFreeFunc) spline1DFree);
+    tmpSpline->xMin = NAN;
+    tmpSpline->xMax = NAN;
+    tmpSpline->xDel = NAN;
+    tmpSpline->equalSpacing = false;
+
+    tmpSpline->xKnots   = NULL;
+    tmpSpline->yKnots   = NULL;
+    tmpSpline->d2yKnots = NULL;
+    psMemSetDeallocator(tmpSpline, (psFreeFunc) psSpline1DFree);
 
     return(tmpSpline);
 }
 
+/** Create an empty 1D spline **/
+psSpline1D *psSpline1DCreate(
+    int nKnots)                 ///< number of knots
+{
+    psSpline1D *spline = psSpline1DAlloc();
+    spline->n = nKnots; // number of knots
+
+    spline->xKnots   = (psF32 *) psAlloc( spline->n * sizeof(psF32));
+    spline->yKnots   = (psF32 *) psAlloc( spline->n * sizeof(psF32));
+    spline->d2yKnots = (psF32 *) psAlloc( spline->n * sizeof(psF32));
+
+    // x & y can both be F32 or F64. should knots be F64?
+    for (psS32 i = 0 ; i < spline->n ; i++) {
+        spline->xKnots[i]   = NAN;
+        spline->yKnots[i]   = NAN;
+        spline->d2yKnots[i] = NAN;
+    }
+    return(spline);
+}
+
 /*****************************************************************************
-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 +193 @@
 XXX: What types must be supported?
  *****************************************************************************/
-psSpline1D *psVectorFitSpline1D(
+psSpline1D *psSpline1DFitVector(
     const psVector* x,                  ///< Ordinates.
-    const psVector* y)                  ///< Coordinates.
+    const psVector* y,                  ///< Coordinates.
+    psF32 dyLower,                      ///< 1st derivative at lower bound
+    psF32 dyUpper)                      ///< 1st derivative at upper bound
 {
     psTrace("psLib.math", 3, "---- %s() begin ----\n", __func__);
+    PS_ASSERT_VECTOR_NON_NULL(x, NULL);
     PS_ASSERT_VECTOR_NON_NULL(y, NULL);
+    PS_ASSERT_VECTORS_SIZE_EQUAL(x, y, NULL);
+    PS_ASSERT_VECTOR_TYPE_F32_OR_F64(x, NULL);
     PS_ASSERT_VECTOR_TYPE_F32_OR_F64(y, NULL);
     PS_ASSERT_LONG_LARGER_THAN_OR_EQUAL(y->n, (long)2, NULL);
-    psS32 numSplines = (y->n)-1;
-    psTrace("psLib.math", 5, "numSplines is %d\n", numSplines);
-
-    //
-    // Create the psSpline1D struct.
-    //
-    psSpline1D *spline = psSpline1DAlloc();
-    spline->n = numSplines;
-    spline->spline = (psPolynomial1D **) psAlloc(numSplines * sizeof(psPolynomial1D *));
-    for (psS32 i=0;i<numSplines;i++) {
-        spline->spline[i] = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 3);
-    }
-
-    //
-    // The following code ensures that xPtr and yPtr points to a psF32 psVector.
-    //
-    // XXX: Use the vector copy and create routines here:
-    //
-
-    spline->knots = psVectorAlloc(y->n, PS_TYPE_F32);
-    if (x != NULL) {
-        PS_ASSERT_VECTOR_NON_NULL(x, NULL);
-        PS_ASSERT_VECTORS_SIZE_EQUAL(x, y, NULL);
-        PS_ASSERT_VECTOR_TYPE_F32_OR_F64(x, NULL);
-        if (x->type.type == PS_TYPE_F32) {
-            for (psS32 i = 0 ; i < x->n ; i++) {
-                spline->knots->data.F32[i] = x->data.F32[i];
-            }
-        } else if (x->type.type == PS_TYPE_F64) {
-            for (psS32 i = 0 ; i < x->n ; i++) {
-                spline->knots->data.F32[i] = (psF32) x->data.F64[i];
-            }
-        }
-    } else {
-        for (psS32 i = 0 ; i < y->n ; i++) {
-            spline->knots->data.F32[i] = (psF32) i;
-        }
-    }
-    psVector *xPtr = spline->knots;
-
-    psVector *yPtr = NULL;
-    // Convert y to F32 if necessary.
-    if (PS_TYPE_F64 == y->type.type) {
-        yPtr = psVectorCopy(NULL, y, PS_TYPE_F32);
-    } else {
-        yPtr = (psVector *) y;
-    }
-
-    //
+
+    psSpline1D *spline = psSpline1DCreate(y->n);
+
+    // x & y can both be F32 or F64. should knots be F64?
+    for (psS32 i = 0 ; i < spline->n ; i++) {
+        spline->xKnots[i] = (x->type.type == PS_TYPE_F32) ? x->data.F32[i] : x->data.F64[i];
+        spline->yKnots[i] = (y->type.type == PS_TYPE_F32) ? y->data.F32[i] : y->data.F64[i];
+    }
+
     // Generate the second derivatives at each data point.
-    //
-    spline->p_psDeriv2 = calculateSecondDerivs(xPtr, yPtr);
-
-    //
-    // We generate the coefficients of the spline polynomials.  I can't
-    // concisely explain how this code works.  See above function comments
-    // and Numerical Recipes in C.
-    //
-    for (psS32 i=0 ; i < numSplines ; i++) {
-        psF32 H = xPtr->data.F32[i+1] - xPtr->data.F32[i];
-        if (fabs(H) <= FLT_EPSILON) {
-            psError(PS_ERR_UNKNOWN, false, "x data points are not distinct (%d %d) (%f %f).\n",
-                    i, i+1, xPtr->data.F32[i], xPtr->data.F32[i+1]);
-        }
-        psTrace("psLib.math", 6, "x data (%f - %f) (%f)\n", xPtr->data.F32[i], xPtr->data.F32[i+1], H);
-        //
-        // ******** Calculate 0-order term ********
-        //
-        // From (1)
-        spline->spline[i]->coeff[0] = yPtr->data.F32[i] * xPtr->data.F32[i+1]/H;
-        // From (2)
-        spline->spline[i]->coeff[0]-= (yPtr->data.F32[i+1] * xPtr->data.F32[i])/H;
-        // From (3)
-        psF32 tmp = (xPtr->data.F32[i+1] * xPtr->data.F32[i+1] * xPtr->data.F32[i+1]) / (H * H * H);
-        tmp-= xPtr->data.F32[i+1] / H;
-        tmp*= spline->p_psDeriv2[i] * H * H / 6.0;
-        spline->spline[i]->coeff[0]+= tmp;
-        // From (4)
-        tmp = -(xPtr->data.F32[i] * xPtr->data.F32[i] * xPtr->data.F32[i]) / (H * H * H);
-        tmp+= xPtr->data.F32[i] / H;
-        tmp*= spline->p_psDeriv2[i+1] * H * H / 6.0;
-        spline->spline[i]->coeff[0]+= tmp;
-
-        //
-        // ******** Calculate 1-order term ********
-        //
-        // From (1)
-        spline->spline[i]->coeff[1] = -(yPtr->data.F32[i]) / H;
-        // From (2)
-        spline->spline[i]->coeff[1]+= yPtr->data.F32[i+1] / H;
-        // From (3)
-        tmp = -3.0 * xPtr->data.F32[i+1] * xPtr->data.F32[i+1] / (H * H * H);
-        tmp+= (1.0 / H);
-        tmp*= spline->p_psDeriv2[i] * H * H / 6.0;
-        spline->spline[i]->coeff[1]+= tmp;
-        // From (4)
-        tmp = 3.0 * xPtr->data.F32[i] * xPtr->data.F32[i] / (H * H * H);
-        tmp-= 1.0 / H;
-        tmp*= spline->p_psDeriv2[i+1] * H * H / 6.0;
-        spline->spline[i]->coeff[1]+= tmp;
-
-        //
-        // ******** Calculate 2-order term ********
-        //
-        // From (3)
-        spline->spline[i]->coeff[2] = spline->p_psDeriv2[i] * 3.0 * xPtr->data.F32[i+1] / (6.0 * H);
-        // From (4)
-        spline->spline[i]->coeff[2]-= spline->p_psDeriv2[i+1] * 3.0 * xPtr->data.F32[i] / (6.0 * H);
-
-        //
-        // ******** Calculate 3-order term ********
-        //
-        // From (3)
-        spline->spline[i]->coeff[3] = -spline->p_psDeriv2[i] / (6.0 * H);
-        // From (4)
-        spline->spline[i]->coeff[3]+=  spline->p_psDeriv2[i+1] / (6.0 * H);
-
-        psTrace("psLib.math", 6, "(spline->spline[%u])->coeff[0] is %f\n", i, spline->spline[i]->coeff[0]);
-        psTrace("psLib.math", 6, "(spline->spline[%u])->coeff[1] is %f\n", i, spline->spline[i]->coeff[1]);
-        psTrace("psLib.math", 6, "(spline->spline[%u])->coeff[2] is %f\n", i, spline->spline[i]->coeff[2]);
-        psTrace("psLib.math", 6, "(spline->spline[%u])->coeff[3] is %f\n", i, spline->spline[i]->coeff[3]);
-    }
-
-    if (PS_TYPE_F64 == y->type.type) {
-        psFree(yPtr);
-    }
+    calculateSecondDerivs(spline, dyLower, dyUpper);
+
     psTrace("psLib.math", 3, "---- %s() end ----\n", __func__);
     return(spline);
 }
 
-
-/*****************************************************************************
-psSpline1DEval(): this routine takes an existing spline of arbitrary order
-and an independent x value.  It determines which spline that x corresponds
-to by doing a bracket disection on the knots of the spline data structure
-(vectorBinDisectF32()).  Then it evaluates the spline at that x location
-by a call to the 1D polynomial functions.
-
-XXX: The spline eval functions require input and output to be F32.  however
-     the spline fit functions require F32 and F64.
-
-XXX: This only works if spline->knots if psF32.  Must we add support for psU32 and
-psF64?
- *****************************************************************************/
+psPolynomial1D *psSpline1DToPoly (psSpline1D *spline, int n) {
+    PS_ASSERT_INT_LESS_THAN(n, spline->n - 1, NULL);
+
+    // convert the cubic spline coeffs to a polynomial. See above function comments and
+    // Numerical Recipes in C.
+
+    psF32 *xKnots   = spline->xKnots;
+    psF32 *yKnots   = spline->yKnots;
+    psF32 *d2yKnots = spline->d2yKnots;
+
+    psF32 H = xKnots[n+1] - xKnots[n];
+    if (fabs(H) <= FLT_EPSILON) {
+        psError(PS_ERR_UNKNOWN, false, "x data points are not distinct (%d %d) (%f %f).\n",
+                n, n+1, xKnots[n], xKnots[n+1]);
+    }
+
+    psPolynomial1D *myPoly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 3);
+
+    // ******** Calculate 0-order term ********
+    // From (1)
+    myPoly->coeff[0] = yKnots[n] * xKnots[n+1]/H;
+    // From (2)
+    myPoly->coeff[0]-= (yKnots[n+1] * xKnots[n])/H;
+    // From (3)
+    psF32 tmp = (xKnots[n+1] * xKnots[n+1] * xKnots[n+1]) / (H * H * H);
+    tmp-= xKnots[n+1] / H;
+    tmp*= d2yKnots[n] * H * H / 6.0;
+    myPoly->coeff[0]+= tmp;
+    // From (4)
+    tmp = -(xKnots[n] * xKnots[n] * xKnots[n]) / (H * H * H);
+    tmp+= xKnots[n] / H;
+    tmp*= d2yKnots[n+1] * H * H / 6.0;
+    myPoly->coeff[0]+= tmp;
+
+    //
+    // ******** Calculate 1-order term ********
+    //
+    // From (1)
+    myPoly->coeff[1] = -(yKnots[n]) / H;
+    // From (2)
+    myPoly->coeff[1]+= yKnots[n+1] / H;
+    // From (3)
+    tmp = -3.0 * xKnots[n+1] * xKnots[n+1] / (H * H * H);
+    tmp+= (1.0 / H);
+    tmp*= d2yKnots[n] * H * H / 6.0;
+    myPoly->coeff[1]+= tmp;
+    // From (4)
+    tmp = 3.0 * xKnots[n] * xKnots[n] / (H * H * H);
+    tmp-= 1.0 / H;
+    tmp*= d2yKnots[n+1] * H * H / 6.0;
+    myPoly->coeff[1]+= tmp;
+
+    //
+    // ******** Calculate 2-order term ********
+    //
+    // From (3)
+    myPoly->coeff[2] = d2yKnots[n] * 3.0 * xKnots[n+1] / (6.0 * H);
+    // From (4)
+    myPoly->coeff[2]-= d2yKnots[n+1] * 3.0 * xKnots[n] / (6.0 * H);
+
+    //
+    // ******** Calculate 3-order term ********
+    //
+    // From (3)
+    myPoly->coeff[3] = -d2yKnots[n] / (6.0 * H);
+    // From (4)
+    myPoly->coeff[3]+=  d2yKnots[n+1] / (6.0 * H);
+
+    psTrace("psLib.math", 6, "(spline->spline[%u])->coeff[0] is %f\n", n, myPoly->coeff[0]);
+    psTrace("psLib.math", 6, "(spline->spline[%u])->coeff[1] is %f\n", n, myPoly->coeff[1]);
+    psTrace("psLib.math", 6, "(spline->spline[%u])->coeff[2] is %f\n", n, myPoly->coeff[2]);
+    psTrace("psLib.math", 6, "(spline->spline[%u])->coeff[3] is %f\n", n, myPoly->coeff[3]);
+
+    return myPoly;
+}
+
+// given an already-constructed spline, check/assert that the
+// knot spacing is equal.  this allows some optimization
+bool psSpline1DisEqualSpacing (psSpline1D *spline) {
+
+    PS_ASSERT_PTR_NON_NULL(spline, false);
+    PS_ASSERT_PTR_NON_NULL(spline->xKnots, false);
+    PS_ASSERT_PTR_NON_NULL(spline->yKnots, false);
+    PS_ASSERT_PTR_NON_NULL(spline->d2yKnots, false);
+    
+    // if the spline has equally-spaced xKnots, the values of the
+    // xKnots can be predicted from the first, last, and delta values
+
+    int n = spline->n;
+    spline->xMax = spline->xKnots[n-1];
+    spline->xMin = spline->xKnots[0];
+    spline->xDel = (spline->xMax - spline->xMin) / (n - 1);
+    
+    // check that the xKnots actually follow this spacing:
+
+    for (int i = 1; i < n - 1; i++) {
+        float xValue = spline->xMin + i*spline->xDel;
+        fprintf (stderr, "%d %f - %f = %f\n", i, spline->xKnots[i], xValue, spline->xKnots[i] - xValue);
+    }
+
+    spline->equalSpacing = true;
+    return true;
+}
+
+// XXX EAM : changing implementation to use yKnot, d2yKnot instead of polynomials
 float psSpline1DEval(
     const psSpline1D *spline,
@@ -368 +334 @@
     PS_ASSERT_PTR_NON_NULL(spline, NAN);
     PS_ASSERT_INT_NONNEGATIVE(spline->n, NAN);
-    PS_ASSERT_VECTOR_TYPE(spline->knots, PS_TYPE_F32, NAN);
+    PS_ASSERT_PTR_NON_NULL(spline->xKnots, 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));
-}
-
+
+    // XXX this should be linear extrapolation at the high or low ends
+    if (x < spline->xKnots[0])   return psSpline1DEval_Segment(spline,   0, x);
+    if (x > spline->xKnots[n-1]) return psSpline1DEval_Segment(spline, n-2, x);
+
+    if (spline->equalSpacing) {
+        int bin = (x - spline->xMin) / spline->xDel;
+        bin = PS_MIN(PS_MAX(bin, 0), n - 2);
+        return psSpline1DEval_Segment(spline, bin, x);
+    }
+
+    /* find correct element in array (x must be sorted) */
+    int lo = 0;
+    int hi = n-1;
+    while (hi - lo > 1) {
+        int i = 0.5*(hi+lo);
+        if (spline->xKnots[i] > x) {
+            hi = i;
+        } else {
+            lo = i;
+        }
+    }
+
+    return psSpline1DEval_Segment(spline, lo, x);
+}
+
+// evaluate the spline at the given coordinate using the specified segment
+// (YMMV if you use the wrong segment!)
+float psSpline1DEval_Segment(
+    const psSpline1D *spline,
+    int n,
+    float x)
+{
+    psTrace("psLib.math", 3, "---- %s() begin ----\n", __func__);
+    PS_ASSERT_PTR_NON_NULL(spline, NAN);
+    PS_ASSERT_PTR_NON_NULL(spline->xKnots, NAN);
+    PS_ASSERT_PTR_NON_NULL(spline->yKnots, NAN);
+    PS_ASSERT_PTR_NON_NULL(spline->d2yKnots, NAN);
+    PS_ASSERT_INT_NONNEGATIVE(spline->n, NAN);
+    PS_ASSERT_INT_LESS_THAN(n, spline->n - 1, NAN);
+
+    psF32 dX = spline->xKnots[n+1] - spline->xKnots[n];
+    psF32 A  = (spline->xKnots[n+1] - x) / dX;
+    psF32 B  = (x - spline->xKnots[n]) / dX;
+
+    psF32 value = A*spline->yKnots[n] + B*spline->yKnots[n+1] + ((A*A*A - A)*spline->d2yKnots[n] + (B*B*B - B)*spline->d2yKnots[n+1])*(dX*dX) / 6.0;
+    return value;
+}
 
 /*****************************************************************************
+ returns a vector of the same type as the input (x)
  *****************************************************************************/
 psVector *psSpline1DEvalVector(
@@ -405 +394 @@
 {
     psTrace("psLib.math", 3, "---- %s() begin ----\n", __func__);
-    PS_ASSERT_PTR_NON_NULL(spline, NULL);
-    PS_ASSERT_VECTOR_TYPE(spline->knots, PS_TYPE_F32, NULL);
+    PS_ASSERT_PTR_NON_NULL(spline,           NULL);
+    PS_ASSERT_PTR_NON_NULL(spline->xKnots,   NULL);
+    PS_ASSERT_PTR_NON_NULL(spline->yKnots,   NULL);
+    PS_ASSERT_PTR_NON_NULL(spline->d2yKnots, NULL);
+    PS_ASSERT_INT_NONNEGATIVE(spline->n,     NULL);
+
     PS_ASSERT_VECTOR_NON_NULL(x, NULL);
     PS_ASSERT_VECTOR_TYPE_F32_OR_F64(x, NULL);
-    if (psTraceGetLevel("psLib.math") >= 6) {
-        PS_VECTOR_PRINT_F32(x);
-        PS_PRINT_SPLINE2((psSpline1D *) spline);
-    }
-
-    psVector *tmpVector = psVectorAlloc(x->n, PS_TYPE_F32);
+
+    psVector *tmpVector = psVectorAlloc(x->n, x->type.type);
     if (x->type.type == PS_TYPE_F32) {
         for (psS32 i=0;i<x->n;i++) {
             tmpVector->data.F32[i] = psSpline1DEval(spline, x->data.F32[i]);
         }
-    } else if (x->type.type == PS_TYPE_F64) {
+    } else {
         for (psS32 i=0;i<x->n;i++) {
-            tmpVector->data.F32[i] = psSpline1DEval(spline, (psF32) x->data.F64[i]);
+            tmpVector->data.F64[i] = psSpline1DEval(spline, (psF32) x->data.F64[i]);
         }
     }
-
+    
     psTrace("psLib.math", 3, "---- %s() end ----\n", __func__);
     return(tmpVector);

trunk/psLib/src/math/psSpline.h

@@ -5 +5 @@
  * and evaluate splines.
  *
- * @author GLG, MHPCC
+ * @author GLG (MHPCC)
+ * reworked by EAM (IfA)
  *
  * @version $Revision: 1.62 $ $Name: not supported by cvs2svn $
@@ -27 +28 @@
 #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.
+    unsigned int n;    ///< The number of knots
+    psF32 *xKnots; ///< x-coordinate of the knots
+    psF32 *yKnots; ///< y-coordinate of the knots
+    psF32 *d2yKnots; ///< 2nd derivative of y at the knots
+    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
  */
 psSpline1D *psSpline1DAlloc(void) PS_ATTR_MALLOC;
+
+/** Create an empty 1D spline **/
+psSpline1D *psSpline1DCreate(
+    int nKnots);                        ///< number of knots
 
 /** Evaluates 1-D spline polynomials at a specific coordinate.
@@ -53 +61 @@
  */
 float psSpline1DEval(
-    const psSpline1D *spline,          ///< Coefficients for spline polynomials
+    const psSpline1D *spline,          ///< spline pointer
+    float x                            ///< location at which to evaluate
+);
+
+/** Evaluates 1-D spline polynomials at a specific coordinate.
+ *
+ *  @return float    result of spline polynomials evaluated at given location
+ */
+float psSpline1DEval_Segment(
+    const psSpline1D *spline,          ///< spline pointer
+    int n,                             // choice of segment
     float x                            ///< location at which to evaluate
 );
@@ -71 +89 @@
  *  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
+    const psVector* y,                  ///< Coordinates.
+    psF32 dyLower,                      ///< 1st derivative at lower bound
+    psF32 dyUpper                       ///< 1st derivative at upper bound
 );
 
@@ -87 +109 @@
     psPtr ptr                          ///< the pointer whose type to check
 );
+
+// check for equal spacing and set internal boolean if true
+bool psSpline1DisEqualSpacing (psSpline1D *spline);
+
+// convert the cubic spline elements to a simply ordinary polynomial for segment n
+psPolynomial1D *psSpline1DToPoly (psSpline1D *spline, int n);
 
 /*****************************************************************************

trunk/psLib/test/math/tap_psSpline1D.c

@@ -1 @@
-/* @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 +35 @@
     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 +85 @@
 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 +107 @@
         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 +146 @@
                     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 +161 @@
                         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 +183 @@
 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 +196 @@
         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 +286 @@
     }
 
-    // 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 +299 @@
         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 +309 @@
         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 +328 @@
         ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
     }
-
 
     ok(genericF32Test(1, myFunc00, false), "Generic, simple mapping, F32 test. 1 spline");