IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 42319


Ignore:
Timestamp:
Jan 23, 2023, 7:10:15 AM (3 years ago)
Author:
eugene
Message:

workin on spline implementation

Location:
branches/eam_branches/psLib.20230123
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/psLib.20230123/src/math/psSpline.c

    r23486 r42319  
    11/** @file psSpline.c
    2 *
    3 *  @brief Contains basic spline allocation, deallocation, fitting,
    4 *         and evaluation routines.
    5 *
    6 *  This file contains the routines that allocate, free, and evaluate splines.
    7 *
    8 *  @version $Revision: 1.158 $ $Name: not supported by cvs2svn $
    9 *  @date $Date: 2007-03-14 02:36:28 $
    10 *
    11 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
     2    re-written by EAM 2023.01.23
    123*/
    134
    … …  
    167#endif
    178
    18 /*****************************************************************************/
    19 /*  INCLUDE FILES                                                            */
    20 /*****************************************************************************/
    219#include <stdio.h>
    2210#include <stdbool.h>
    … …  
    3624#include "psMathUtils.h"
    3725
    38 /*****************************************************************************/
    39 /* DEFINE STATEMENTS                                                         */
    40 /*****************************************************************************/
    41 
    42 /*****************************************************************************/
    43 /* TYPE DEFINITIONS                                                          */
    44 /*****************************************************************************/
    45 
    46 /*****************************************************************************/
    47 /* GLOBAL VARIABLES                                                          */
    48 /*****************************************************************************/
    49 
    50 /*****************************************************************************/
    51 /* FILE STATIC VARIABLES                                                     */
    52 /*****************************************************************************/
    53 
    54 /*****************************************************************************/
    55 /* FUNCTION IMPLEMENTATION - LOCAL                                           */
    56 /*****************************************************************************/
    57 
    58 static void spline1DFree(psSpline1D *tmpSpline)
    59 {
    60     if (tmpSpline == NULL) {
    61         return;
    62     }
     26static void psSpline1DFree(psSpline1D *tmpSpline)
     27{
     28    if (tmpSpline == NULL) return;
    6329
    6430    if (tmpSpline->spline != NULL) {
    … …  
    7844}
    7945
    80 
    81 static void PS_POLY1D_PRINT(
    82     psPolynomial1D *poly)
    83 {
    84     printf("-------------- PS_POLY1D_PRINT() --------------\n");
    85     printf("poly->nX is %d\n", poly->nX);
    86     for (psS32 i = 0 ; i < (1 + poly->nX) ; i++) {
    87         printf("poly->coeff[%d] is %f\n", i, poly->coeff[i]);
    88     }
    89 }
    90 
    9146static void PS_PRINT_SPLINE2(psSpline1D *mySpline)
    9247{
    … …  
    9449    printf("mySpline->n is %d\n", mySpline->n);
    9550    for (psS32 i = 0 ; i < mySpline->n ; i++) {
    96         PS_POLY1D_PRINT(mySpline->spline[i]);
    97     }
    98     PS_VECTOR_PRINT_F32(mySpline->knots);
     51        // print the knots
     52    }
    9953}
    10054
    10155/*****************************************************************************
    102 CalculateSecondDerivs(): Given a set of x/y vectors corresponding to a
    103 tabulated function at n points, this routine calculates the second
    104 derivatives of the interpolating cubic splines at those n points.
    105 
    106 The first and second derivatives at the endpoints were previously undefined in
    107 the SDR.  From bugzilla #???, they are required to be 0.0, implementing natural
    108 splines.
    109 
    110 Endpoints are defined by
    111     PS_LEFT_SPLINE_DERIV
    112     PS_RIGHT_SPLINE_DERIV
    113 
    114 This routine assumes that vectors x and y are of the appropriate types/sizes
    115 (F32).
    116 
    117 XXX: use recycled vectors for internal data.
    118 XXX: do an F64 version?
     56CalculateSecondDerivs(): Given a set of x,y vectors corresponding to the spline knots, this
     57routine calculates the second derivatives of the interpolating cubic splines at those n points.
     58
     59The boundary conditions may be:
     60* first derivative specified OR
     61* second derivatives = 0
     62
     63One or the other condition may be true for each of the upper and lower bounds
     64
     65NOTE: vectors must be F32
     66
    11967 *****************************************************************************/
    120 #define PS_LEFT_SPLINE_DERIV 0.0
    121 #define PS_RIGHT_SPLINE_DERIV 0.0
     68
    12269static psF32 *calculateSecondDerivs(
    123     const psVector* x,                  ///< Ordinates
    124     const psVector* y)                  ///< Coordinates
     70    const psVector* x, ///< Ordinates
     71    const psVector* y, ///< Coordinates
     72    psF32 dyLower; // if not NAN, lower-bound 1st derivative is defined
     73    psF32 dyUpper; // if not NAN, lower-bound 1st derivative is defined
     74    )                 
    12575{
    12676    psTrace("psLib.math", 4, "---- %s() begin ----\n", __func__);
    … …  
    190140
    191141/*****************************************************************************
    192 psVectorFitSpline1D(): given a set of x/y vectors, this routine generates the
     142psSpline1DFitVector(): given a set of x/y vectors, this routine generates the
    193143linear or cublic splines which satisfy those data points.
    194144
    … …  
    209159XXX: What types must be supported?
    210160 *****************************************************************************/
    211 psSpline1D *psVectorFitSpline1D(
     161psSpline1D *psSpline1DFitVector(
    212162    const psVector* x,                  ///< Ordinates.
    213163    const psVector* y)                  ///< Coordinates.
    … …  
    361311psF64?
    362312 *****************************************************************************/
    363 float psSpline1DEval(
    364     const psSpline1D *spline,
     313float psSpline1DEval_Old(
     314    const psSpline1D_Old *spline,
    365315    float x)
    366316{
    … …  
    397347}
    398348
     349// XXX EAM : changing implementation to use yKnot, d2yKnot instead of polynomials
     350float psSpline1DEval_New(
     351    const psSpline1D *spline,
     352    float x)
     353{
     354    psTrace("psLib.math", 3, "---- %s() begin ----\n", __func__);
     355    PS_ASSERT_PTR_NON_NULL(spline, NAN);
     356    PS_ASSERT_INT_NONNEGATIVE(spline->n, NAN);
     357    PS_ASSERT_VECTOR_TYPE(spline->knots, PS_TYPE_F32, NAN);
     358
     359    psS32 n = spline->n;
     360    if ((x < spline->knots->data.F32[0]) || (x > spline->knots->data.F32[spline->knots->n-1])) {
     361        // If x is outside the range of spline->knots, generate a warning
     362        // message, then return the left, or right, endpoint.
     363        psLogMsg(__func__, PS_LOG_WARN,
     364                 "psSpline1DEval(): x ordinate (%f) is outside the spline range (%f - %f) (%d).",
     365                 x, spline->knots->data.F32[0], spline->knots->data.F32[n-1], n);
     366
     367        psS32 binNum = (x < spline->knots->data.F32[0]) ? 0 : n-1;
     368        psTrace("psLib.math", 3, "---- %s() end ----\n", __func__);
     369        return(psPolynomial1DEval(spline->spline[binNum], x));
     370    }
     371
     372    psScalar tmpScalar;
     373    tmpScalar.type.type = PS_TYPE_F32;
     374    tmpScalar.data.F32 = x;
     375    psVectorBinaryDisectResult result;
     376    psS32 binNum = psVectorBinaryDisect(&result, spline->knots, &tmpScalar);
     377    if (result != PS_BINARY_DISECT_PASS) {
     378        psError(PS_ERR_UNKNOWN, false, "Could not perform bin dissection on spline->knots.\n");
     379        return(NAN);
     380    }
     381
     382    psTrace("psLib.math", 3, "---- %s() end ----\n", __func__);
     383    return(psPolynomial1DEval(spline->spline[binNum], x));
     384}
     385
    399386
    400387/*****************************************************************************

  • branches/eam_branches/psLib.20230123/src/math/psSpline.h

    r23486 r42319  
    55 * and evaluate splines.
    66 *
    7  * @author GLG, MHPCC
     7 * @author GLG (MHPCC)
     8 * reworked by EAM (IfA)
    89 *
    910 * @version $Revision: 1.62 $ $Name: not supported by cvs2svn $
    … …  
    2728#include "psPolynomial.h"
    2829
    29 #define PS_LEFT_SPLINE_DERIV 0.0
    30 #define PS_RIGHT_SPLINE_DERIV 0.0
    31 
    32 /** One-Dimensional Spline */
     30/** One-Dimensional Spline
     31    This structure represents a 1D cubic spline.  Note the option for
     32    equally-spaced knots allows a quick selection of the correct spline
     33    segment.  The values (xMin, xMax, xDel) are stored for ease of access.
     34 */
    3335typedef struct
    3436{
    3537    unsigned int n;                    ///< The number of spline pieces
    36     psPolynomial1D **spline;           ///< An array of n pointers to the spline polynomials
    37     psVector *knots;                   ///< The boundaries between each spline piece.  Size is n+1.
    38     psF32 *p_psDeriv2;                 ///< For cubic splines, the second derivative at each domain point.  Size is n+1.
     38    psF32 *xKnots; ///< x-coordinate of the knots (n+1)
     39    psF32 *yKnots; ///< y-coordinate of the knots (n+1)
     40    psF32 *d2yKnots; ///< 2nd derivative of y at the knots (n+1)
     41    bool equalSpacing; // if knots are equally spaced, the seqment choice can be optimized
     42    psF32 xMin; // for equally-spaced knots, the value at the lower bound (xKnots[0])
     43    psF32 xMax; // for equally-spaced knots, the value at the upper bound (xKnots[n])
     44    psF32 xDel; // for equally-spaced knots, the spacing (xKnots[n] - xKnots[0])/n
    3945}
    4046psSpline1D;
    4147
    42 /** Allocates a psSpline1D structure
    43  *
    44  *  Allocator for psSpline1D.
     48/** Allocator for psSpline1D.
    4549 *
    4650 *  @return psSpline1D*    new 1-D spline struct
    … …  
    5357 */
    5458float psSpline1DEval(
    55     const psSpline1D *spline,          ///< Coefficients for spline polynomials
     59    const psSpline1D *spline,          ///< spline pointer
    5660    float x                            ///< location at which to evaluate
    5761);
    … …  
    7175 *  generates the linear splines which satisfy those data points.
    7276 *
     77 *  XXX EAM: add option to generate / select a subset of knots from the data
     78 *
    7379 *  @return psSpline1D*:  the calculated one-dimensional splines
    7480 */
    75 psSpline1D *psVectorFitSpline1D(
     81psSpline1D *psSpline1DFitVector(
    7682    const psVector* x,                 ///< Ordinates (or NULL to just use the indices)
    7783    const psVector* y                  ///< Coordinates

  • branches/eam_branches/psLib.20230123/test/math/tap_psSpline1D.c

    r13124 r42319  
    1 /* @file  tst_psImageManip.c
    2 *
    3 *  @brief This file will contain tests for all of the public psLib functions
     1/* @file  tap_psImageManip.c
     2*
     3*  @brief This file contains tests for all of the public psLib functions
    44*         that implement 1-D spline functionality:
    55* psSpline1DAlloc()
    66* psSpline1DEval()
    77* psSpline1DEvalVector()
    8 * psVectorFitSpline1D()
     8* psSpline1DFitVector()
    99*
    1010*         This file is composed of the tests formerly in tst_psFunc02.c,
    … …  
    3535    psF32 expect)
    3636{
     37    // NOTE: this returns NAN if 'expect' == 0.0
     38    // NOTE 2: this is not testing a percent, but fractional error
    3739    if ((fabs(actual - expect) / fabs(expect)) > ERROR_TOLERANCE_PERCENT) {
    3840        return(true);
    … …  
    8385typedef psF64 (*mappingFuncF64)(psF64 x);
    8486
     87/* EAM 2023.01.22 : these tests are not well considered.  They generate a set of N+1 points
     88   to use as knots at integer values from 0 to N, then evaluate the spline half-way between
     89   the knots.   the function above is a cubic, so a cubic spline should fit it perfectly, which
     90   is fine.  Perhaps this test suite should use a pre-defined collection of data points?
     91 */
     92
    8593bool genericF32Test(psS32 NumSplines, mappingFuncF32 func, bool xNull)
    8694{
    87     // We test the psVectorFitSpline1D, psSpline1DEval() functions.  F32 version.
     95    // We test the psSpline1DFitVector, psSpline1DEval() functions.  F32 version.
    8896    bool testStatus = true;
    8997    {
     98        // Generate the vector data
    9099        psMemId id = psMemGetId();
    91100        psVector *xF32 = psVectorAlloc(NumSplines+1, PS_TYPE_F32);
    … …  
    98107        psSpline1D *tmpSpline = NULL;
    99108        if (!xNull) {
    100             tmpSpline = psVectorFitSpline1D(xF32, yF32);
     109            tmpSpline = psSpline1DFitVector(xF32, yF32);
    101110        } else {
    102             tmpSpline = psVectorFitSpline1D(NULL, yF32);
    103         }
     111            tmpSpline = psSpline1DFitVector(NULL, yF32);
     112        }
     113
    104114        if(tmpSpline == NULL) {
    105             diag("psVectorFitSpline1D() returned NULL");
     115            diag("psSpline1DFitVector() returned NULL");
    106116            testStatus = false;
    107117        } else {
    108118            if (tmpSpline->n != NumSplines) {
    109                 diag("psVectorFitSpline1D() did not properly set the psSpline1D->n member");
     119                diag("psSpline1DFitVector() did not properly set the psSpline1D->n member");
    110120                testStatus = false;
    111121            }
    112122            if(tmpSpline->spline == NULL) {
    113                 diag("psVectorFitSpline1D() returned a NULL psSpline1D->spline member.");
     123                diag("psSpline1DFitVector() returned a NULL psSpline1D->spline member.");
    114124                testStatus = false;
    115125            }
    116126            for (psS32 i = 0 ; i < NumSplines ; i++) {
    117127                if (tmpSpline->spline[i] == NULL) {
    118                     diag("psVectorFitSpline1D() returned a NULL psSpline1D->spline[%d] member.", i);
     128                    diag("psSpline1DFitVector() returned a NULL psSpline1D->spline[%d] member.", i);
    119129                    testStatus = false;
    120130                }
    121131            }
    122132            if (tmpSpline->knots == NULL) {
    123                 diag("psVectorFitSpline1D() returned a NULL psSpline1D->knots member");
     133                diag("psSpline1DFitVector() returned a NULL psSpline1D->knots member");
    124134                testStatus = false;
    125135            }
    126136            if (tmpSpline->p_psDeriv2 == NULL) {
    127                 diag("psVectorFitSpline1D()returned a NULL psSpline1D->p_psDeriv2 member");
     137                diag("psSpline1DFitVector()returned a NULL psSpline1D->p_psDeriv2 member");
    128138                testStatus = false;
    129139            }
    … …  
    136146                    testStatus = false;
    137147                    diag("TEST ERROR: f(%f) is %f, should be %f", x, y, myFunc00(x));
     148                    // XXX EAM : the truth value above should be 'func' not myFunc00
    138149                }
    139150            }
    … …  
    150161                        diag("TEST ERROR: f(%f) is %f, should be %f", xF32->data.F32[i],
    151162                             yF32Test->data.F32[i], myFunc00(xF32->data.F32[i]));
     163                    // XXX EAM : the truth value above should be 'func' not myFunc00
    152164                    }
    153165                }
    … …  
    171183bool genericF64Test(psS32 NumSplines, mappingFuncF64 func, bool xNull)
    172184{
    173     // We test the psVectorFitSpline1D, psSpline1DEval() functions.  F64 version.
     185    // We test the psSpline1DFitVector, psSpline1DEval() functions.  F64 version.
    174186    bool testStatus = true;
    175187    {
    … …  
    184196        psSpline1D *tmpSpline = NULL;
    185197        if (!xNull) {
    186             tmpSpline = psVectorFitSpline1D(xF64, yF64);
     198            tmpSpline = psSpline1DFitVector(xF64, yF64);
    187199        } else {
    188             tmpSpline = psVectorFitSpline1D(NULL, yF64);
     200            tmpSpline = psSpline1DFitVector(NULL, yF64);
    189201        }
    190202        if(tmpSpline == NULL) {
    191             diag("psVectorFitSpline1D() returned NULL");
     203            diag("psSpline1DFitVector() returned NULL");
    192204            testStatus = false;
    193205        } else {
    194206            if (tmpSpline->n != NumSplines) {
    195                 diag("psVectorFitSpline1D() did not properly set the psSpline1D->n member");
     207                diag("psSpline1DFitVector() did not properly set the psSpline1D->n member");
    196208                testStatus = false;
    197209            }
    198210            if(tmpSpline->spline == NULL) {
    199                 diag("psVectorFitSpline1D() returned a NULL psSpline1D->spline member.");
     211                diag("psSpline1DFitVector() returned a NULL psSpline1D->spline member.");
    200212                testStatus = false;
    201213            }
    202214            for (psS32 i = 0 ; i < NumSplines ; i++) {
    203215                if (tmpSpline->spline[i] == NULL) {
    204                     diag("psVectorFitSpline1D() returned a NULL psSpline1D->spline[%d] member.", i);
     216                    diag("psSpline1DFitVector() returned a NULL psSpline1D->spline[%d] member.", i);
    205217                    testStatus = false;
    206218                }
    207219            }
    208220            if (tmpSpline->knots == NULL) {
    209                 diag("psVectorFitSpline1D() returned a NULL psSpline1D->knots member");
     221                diag("psSpline1DFitVector() returned a NULL psSpline1D->knots member");
    210222                testStatus = false;
    211223            }
    212224            if (tmpSpline->p_psDeriv2 == NULL) {
    213                 diag("psVectorFitSpline1D()returned a NULL psSpline1D->p_psDeriv2 member");
     225                diag("psSpline1DFitVector()returned a NULL psSpline1D->p_psDeriv2 member");
    214226                testStatus = false;
    215227            }
    … …  
    274286    }
    275287
    276     // psSplineEvalTest_sub(): Call psVectorFitSpline1D with NULL arguments.
    277     {
    278         psMemId id = psMemGetId();
    279         psSpline1D *tmpSpline = psVectorFitSpline1D(NULL, NULL);
    280         ok(tmpSpline == NULL, "psVectorFitSpline1D() returns NULL with NULL arguments");
     288    // psSplineEvalTest_sub(): Call psSpline1DFitVector with NULL arguments.
     289    {
     290        psMemId id = psMemGetId();
     291        psSpline1D *tmpSpline = psSpline1DFitVector(NULL, NULL);
     292        ok(tmpSpline == NULL, "psSpline1DFitVector() returns NULL with NULL arguments");
    281293        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
    282294    }
    … …  
    287299        psMemId id = psMemGetId();
    288300        float y = psSpline1DEval(NULL, 0.0);
    289         ok(isnan(y), "psSpline1DEval() returned NAN will NULL input spline");
     301        ok(isnan(y), "psSpline1DEval() returned NAN with NULL input spline");
    290302        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
    291303    }
    … …  
    297309        psVector *x = psVectorAlloc(10, PS_TYPE_F32);
    298310        psVector *y = psSpline1DEvalVector(NULL, x);
    299         ok(y == NULL, "psSpline1DEvalVector() returned NAN will NULL input spline");
     311        ok(y == NULL, "psSpline1DEvalVector() returned NULL with NULL input spline");
    300312        psFree(x);
    301313        psFree(y);
    … …  
    316328        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
    317329    }
    318 
    319330
    320331    ok(genericF32Test(1, myFunc00, false), "Generic, simple mapping, F32 test. 1 spline");
Note: See TracChangeset for help on using the changeset viewer.