IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 5066


Ignore:
Timestamp:
Sep 19, 2005, 9:53:13 AM (21 years ago)
Author:
drobbin
Message:

Made changes to poly struct (psMaskType) and spline fxns (unsigned int's)

Location:
trunk/psLib/src/math
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • trunk/psLib/src/math/psPolynomial.c

    r4991 r5066  
    77*  polynomials.  It also contains a Gaussian functions.
    88*
    9 *  @version $Revision: 1.120 $ $Name: not supported by cvs2svn $
    10 *  @date $Date: 2005-09-11 22:18:40 $
     9*  @version $Revision: 1.121 $ $Name: not supported by cvs2svn $
     10*  @date $Date: 2005-09-19 19:53:13 $
    1111*
    1212*  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    162162outer coefficients of the Chebyshev polynomials.
    163163 *****************************************************************************/
    164 static psPolynomial1D **createChebyshevPolys(psS32 maxChebyPoly)
     164static psPolynomial1D **createChebyshevPolys(unsigned int maxChebyPoly)
    165165{
    166166    PS_ASSERT_INT_NONNEGATIVE(maxChebyPoly, NULL);
    … …  
    169169
    170170    chebPolys = (psPolynomial1D **) psAlloc(maxChebyPoly * sizeof(psPolynomial1D *));
    171     for (psS32 i = 0; i < maxChebyPoly; i++) {
     171    for (unsigned int i = 0; i < maxChebyPoly; i++) {
    172172        chebPolys[i] = psPolynomial1DAlloc(i + 1, PS_POLYNOMIAL_ORD);
    173173    }
    … …  
    191191    } else {
    192192        // XXX: Code this.
    193         printf("WARNING: %d-order chebyshev polynomials not correctly implemented.\n", maxChebyPoly);
     193        printf("WARNING: %u-order chebyshev polynomials not correctly implemented.\n", maxChebyPoly);
    194194    }
    195195
    … …  
    200200    Polynomial coefficients will be accessed in [w][x][y][z] fashion.
    201201 *****************************************************************************/
    202 static psF64 ordPolynomial1DEval(psF64 x, const psPolynomial1D* poly)
    203 {
    204     psS32 loop_x = 0;
     202static psF64 ordPolynomial1DEval(psF64 x,
     203                                 const psPolynomial1D* poly)
     204{
     205    unsigned int loop_x = 0;
    205206    psF64 polySum = 0.0;
    206207    psF64 xSum = 1.0;
    … …  
    209210            "---- Calling ordPolynomial1DEval(%lf)\n", x);
    210211    psTrace(".psLib.dataManip.psPolynomial.ordPolynomial1DEval", 4,
    211             "Polynomial order is %d\n", poly->n);
     212            "Polynomial order is %u\n", poly->n);
    212213    for (loop_x = 0; loop_x < poly->n; loop_x++) {
    213214        psTrace(".psLib.dataManip.psPolynomial.ordPolynomial1DEval", 4,
    214                 "Polynomial coeff[%d] is %lf\n", loop_x, poly->coeff[loop_x]);
     215                "Polynomial coeff[%u] is %lf\n", loop_x, poly->coeff[loop_x]);
    215216    }
    216217
    … …  
    230231// XXX: How does the mask vector effect Crenshaw's formula?
    231232// XXX: We assume that x is scaled between -1.0 and 1.0;
    232 static psF64 chebPolynomial1DEval(psF64 x, const psPolynomial1D* poly)
     233static psF64 chebPolynomial1DEval(psF64 x,
     234                                  const psPolynomial1D* poly)
    233235{
    234236    PS_ASSERT_DOUBLE_WITHIN_RANGE(x, -1.0, 1.0, 0.0);
    235237    // XXX: Create a macro for this in psConstants.h
    236238    if (poly->n < 1) {
    237         psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Error: Chebyshev polynomial as order %d.", poly->n);
     239        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Error: Chebyshev polynomial as order %u.", poly->n);
    238240        return(NAN);
    239241    }
    240242    psVector *d;
    241     psS32 n = poly->n;
    242     psS32 i;
     243    unsigned int n = poly->n;
     244    unsigned int i;
    243245    psF64 tmp = 0.0;
    244246
    … …  
    318320    PS_ASSERT_POLY_NON_NULL(poly, NAN);
    319321
    320     psS32 loop_x = 0;
    321     psS32 loop_y = 0;
     322    unsigned int loop_x = 0;
     323    unsigned int loop_y = 0;
    322324    psF64 polySum = 0.0;
    323325    psF64 xSum = 1.0;
    … …  
    338340}
    339341
    340 static psF64 chebPolynomial2DEval(psF64 x, psF64 y, const psPolynomial2D* poly)
     342static psF64 chebPolynomial2DEval(psF64 x,
     343                                  psF64 y,
     344                                  const psPolynomial2D* poly)
    341345{
    342346    PS_ASSERT_DOUBLE_WITHIN_RANGE(x, -1.0, 1.0, 0.0);
    … …  
    344348    PS_ASSERT_POLY_NON_NULL(poly, NAN);
    345349
    346     psS32 loop_x = 0;
    347     psS32 loop_y = 0;
    348     psS32 i = 0;
     350    unsigned int loop_x = 0;
     351    unsigned int loop_y = 0;
     352    unsigned int i = 0;
    349353    psF64 polySum = 0.0;
    350354    psPolynomial1D* *chebPolys = NULL;
    351     psS32 maxChebyPoly = 0;
     355    unsigned int maxChebyPoly = 0;
    352356
    353357    // Determine how many Chebyshev polynomials
    … …  
    375379}
    376380
    377 static psF64 ordPolynomial3DEval(psF64 x, psF64 y, psF64 z, const psPolynomial3D* poly)
    378 {
    379     psS32 loop_x = 0;
    380     psS32 loop_y = 0;
    381     psS32 loop_z = 0;
     381static psF64 ordPolynomial3DEval(psF64 x,
     382                                 psF64 y,
     383                                 psF64 z,
     384                                 const psPolynomial3D* poly)
     385{
     386    unsigned int loop_x = 0;
     387    unsigned int loop_y = 0;
     388    unsigned int loop_z = 0;
    382389    psF64 polySum = 0.0;
    383390    psF64 xSum = 1.0;
    … …  
    403410}
    404411
    405 static psF64 chebPolynomial3DEval(psF64 x, psF64 y, psF64 z, const psPolynomial3D* poly)
     412static psF64 chebPolynomial3DEval(psF64 x,
     413                                  psF64 y,
     414                                  psF64 z,
     415                                  const psPolynomial3D* poly)
    406416{
    407417    PS_ASSERT_DOUBLE_WITHIN_RANGE(x, -1.0, 1.0, 0.0);
    408418    PS_ASSERT_DOUBLE_WITHIN_RANGE(y, -1.0, 1.0, 0.0);
    409419    PS_ASSERT_DOUBLE_WITHIN_RANGE(z, -1.0, 1.0, 0.0);
    410     psS32 loop_x = 0;
    411     psS32 loop_y = 0;
    412     psS32 loop_z = 0;
    413     psS32 i = 0;
     420    unsigned int loop_x = 0;
     421    unsigned int loop_y = 0;
     422    unsigned int loop_z = 0;
     423    unsigned int i = 0;
    414424    psF64 polySum = 0.0;
    415425    psPolynomial1D* *chebPolys = NULL;
    416     psS32 maxChebyPoly = 0;
     426    unsigned int maxChebyPoly = 0;
    417427
    418428    // Determine how many Chebyshev polynomials
    … …  
    447457}
    448458
    449 static psF64 ordPolynomial4DEval(psF64 x, psF64 y, psF64 z, psF64 t, const psPolynomial4D* poly)
    450 {
    451     psS32 loop_x = 0;
    452     psS32 loop_y = 0;
    453     psS32 loop_z = 0;
    454     psS32 loop_t = 0;
     459static psF64 ordPolynomial4DEval(psF64 x,
     460                                 psF64 y,
     461                                 psF64 z,
     462                                 psF64 t,
     463                                 const psPolynomial4D* poly)
     464{
     465    unsigned int loop_x = 0;
     466    unsigned int loop_y = 0;
     467    unsigned int loop_z = 0;
     468    unsigned int loop_t = 0;
    455469    psF64 polySum = 0.0;
    456470    psF64 xSum = 1.0;
    … …  
    481495}
    482496
    483 static psF64 chebPolynomial4DEval(psF64 x, psF64 y, psF64 z, psF64 t, const psPolynomial4D* poly)
     497static psF64 chebPolynomial4DEval(psF64 x,
     498                                  psF64 y,
     499                                  psF64 z,
     500                                  psF64 t,
     501                                  const psPolynomial4D* poly)
    484502{
    485503    PS_ASSERT_DOUBLE_WITHIN_RANGE(x, -1.0, 1.0, 0.0);
    … …  
    487505    PS_ASSERT_DOUBLE_WITHIN_RANGE(z, -1.0, 1.0, 0.0);
    488506    PS_ASSERT_DOUBLE_WITHIN_RANGE(t, -1.0, 1.0, 0.0);
    489     psS32 loop_x = 0;
    490     psS32 loop_y = 0;
    491     psS32 loop_z = 0;
    492     psS32 loop_t = 0;
    493     psS32 i = 0;
     507    unsigned int loop_x = 0;
     508    unsigned int loop_y = 0;
     509    unsigned int loop_z = 0;
     510    unsigned int loop_t = 0;
     511    unsigned int i = 0;
    494512    psF64 polySum = 0.0;
    495513    psPolynomial1D* *chebPolys = NULL;
    496     psS32 maxChebyPoly = 0;
     514    unsigned int maxChebyPoly = 0;
    497515
    498516    // Determine how many Chebyshev polynomials
    … …  
    542560    evaluated Gaussian is: \f[ exp(-\frac{(x-mean)^2}{2\sigma^2}) \f]
    543561 *****************************************************************************/
    544 float psGaussian(float x, float mean, float sigma, bool normal)
     562float psGaussian(float x,
     563                 float mean,
     564                 float sigma,
     565                 bool normal)
    545566{
    546567    psF32 tmp = 1.0;
    … …  
    568589 *****************************************************************************/
    569590#define PS_XXX_GAUSSIAN_SEED 1995
    570 psVector* p_psGaussianDev(psF32 mean, psF32 sigma, psS32 Npts)
     591psVector* p_psGaussianDev(psF32 mean,
     592                          psF32 sigma,
     593                          unsigned int Npts)
    571594{
    572595    PS_ASSERT_INT_NONNEGATIVE(Npts, NULL);
    … …  
    575598    psRandom *r = psRandomAlloc(PS_RANDOM_TAUS, PS_XXX_GAUSSIAN_SEED);
    576599    psVector* gauss = psVectorAlloc(Npts, PS_TYPE_F32);
    577     for (psS32 i = 0; i < Npts; i++) {
     600    for (unsigned int i = 0; i < Npts; i++) {
    578601        gauss->data.F32[i] = mean + p_psRandomGaussian(r, sigma);
    579602    }
    … …  
    586609    This routine must allocate memory for the polynomial structures.
    587610 *****************************************************************************/
    588 psPolynomial1D* psPolynomial1DAlloc(int n,
     611psPolynomial1D* psPolynomial1DAlloc(unsigned int n,
    589612                                    psPolynomialType type)
    590613{
    591614    PS_ASSERT_INT_POSITIVE(n, NULL);
    592615
    593     int i = 0;
     616    unsigned int i = 0;
    594617    psPolynomial1D* newPoly = NULL;
    595618
    … …  
    601624    newPoly->coeff = psAlloc(n * sizeof(psF64));
    602625    newPoly->coeffErr = psAlloc(n * sizeof(psF64));
    603     newPoly->mask = (char *)psAlloc(n * sizeof(char));
     626    newPoly->mask = (psMaskType *)psAlloc(n * sizeof(psMaskType));
    604627    for (i = 0; i < n; i++) {
    605628        newPoly->coeff[i] = 0.0;
    … …  
    611634}
    612635
    613 psPolynomial2D* psPolynomial2DAlloc( int nX,  int nY,
     636psPolynomial2D* psPolynomial2DAlloc( unsigned int nX,
     637                                     unsigned int nY,
    614638                                     psPolynomialType type)
    615639{
    … …  
    617641    PS_ASSERT_INT_POSITIVE(nY, NULL);
    618642
    619     int x = 0;
    620     int y = 0;
     643    unsigned int x = 0;
     644    unsigned int y = 0;
    621645    psPolynomial2D* newPoly = NULL;
    622646
    … …  
    630654    newPoly->coeff = psAlloc(nX * sizeof(psF64 *));
    631655    newPoly->coeffErr = psAlloc(nX * sizeof(psF64 *));
    632     newPoly->mask = (char **)psAlloc(nX * sizeof(char *));
     656    newPoly->mask = (psMaskType **)psAlloc(nX * sizeof(psMaskType *));
    633657    for (x = 0; x < nX; x++) {
    634658        newPoly->coeff[x] = psAlloc(nY * sizeof(psF64));
    635659        newPoly->coeffErr[x] = psAlloc(nY * sizeof(psF64));
    636         newPoly->mask[x] = (char *)psAlloc(nY * sizeof(char));
     660        newPoly->mask[x] = (psMaskType *)psAlloc(nY * sizeof(psMaskType));
    637661    }
    638662    for (x = 0; x < nX; x++) {
    … …  
    647671}
    648672
    649 psPolynomial3D* psPolynomial3DAlloc( int nX,  int nY,  int nZ,
     673psPolynomial3D* psPolynomial3DAlloc( unsigned int nX,
     674                                     unsigned int nY,
     675                                     unsigned int nZ,
    650676                                     psPolynomialType type)
    651677{
    … …  
    654680    PS_ASSERT_INT_POSITIVE(nZ, NULL);
    655681
    656     psS32 x = 0;
    657     psS32 y = 0;
    658     psS32 z = 0;
     682    unsigned int x = 0;
     683    unsigned int y = 0;
     684    unsigned int z = 0;
    659685    psPolynomial3D* newPoly = NULL;
    660686
    … …  
    669695    newPoly->coeff = psAlloc(nX * sizeof(psF64 **));
    670696    newPoly->coeffErr = psAlloc(nX * sizeof(psF64 **));
    671     newPoly->mask = (char ***)psAlloc(nX * sizeof(char **));
     697    newPoly->mask = (psMaskType ***)psAlloc(nX * sizeof(psMaskType **));
    672698    for (x = 0; x < nX; x++) {
    673699        newPoly->coeff[x] = psAlloc(nY * sizeof(psF64 *));
    674700        newPoly->coeffErr[x] = psAlloc(nY * sizeof(psF64 *));
    675         newPoly->mask[x] = (char **)psAlloc(nY * sizeof(char *));
     701        newPoly->mask[x] = (psMaskType **)psAlloc(nY * sizeof(psMaskType *));
    676702        for (y = 0; y < nY; y++) {
    677703            newPoly->coeff[x][y] = psAlloc(nZ * sizeof(psF64));
    678704            newPoly->coeffErr[x][y] = psAlloc(nZ * sizeof(psF64));
    679             newPoly->mask[x][y] = (char *)psAlloc(nZ * sizeof(char));
     705            newPoly->mask[x][y] = (psMaskType *)psAlloc(nZ * sizeof(psMaskType));
    680706        }
    681707    }
    … …  
    693719}
    694720
    695 psPolynomial4D* psPolynomial4DAlloc( int nX,  int nY,  int nZ,  int nT,
     721psPolynomial4D* psPolynomial4DAlloc( unsigned int nX,
     722                                     unsigned int nY,
     723                                     unsigned int nZ,
     724                                     unsigned int nT,
    696725                                     psPolynomialType type)
    697726{
    … …  
    701730    PS_ASSERT_INT_POSITIVE(nT, NULL);
    702731
    703     psS32 x = 0;
    704     psS32 y = 0;
    705     psS32 z = 0;
    706     psS32 t = 0;
     732    unsigned int x = 0;
     733    unsigned int y = 0;
     734    unsigned int z = 0;
     735    unsigned int t = 0;
    707736    psPolynomial4D* newPoly = NULL;
    708737
    … …  
    718747    newPoly->coeff = psAlloc(nX * sizeof(psF64 ***));
    719748    newPoly->coeffErr = psAlloc(nX * sizeof(psF64 ***));
    720     newPoly->mask = (char ****)psAlloc(nX * sizeof(char ***));
     749    newPoly->mask = (psMaskType ****)psAlloc(nX * sizeof(psMaskType ***));
    721750    for (x = 0; x < nX; x++) {
    722751        newPoly->coeff[x] = psAlloc(nY * sizeof(psF64 **));
    723752        newPoly->coeffErr[x] = psAlloc(nY * sizeof(psF64 **));
    724         newPoly->mask[x] = (char ***)psAlloc(nY * sizeof(char **));
     753        newPoly->mask[x] = (psMaskType ***)psAlloc(nY * sizeof(psMaskType **));
    725754        for (y = 0; y < nY; y++) {
    726755            newPoly->coeff[x][y] = psAlloc(nZ * sizeof(psF64 *));
    727756            newPoly->coeffErr[x][y] = psAlloc(nZ * sizeof(psF64 *));
    728             newPoly->mask[x][y] = (char **)psAlloc(nZ * sizeof(char *));
     757            newPoly->mask[x][y] = (psMaskType **)psAlloc(nZ * sizeof(psMaskType *));
    729758            for (z = 0; z < nZ; z++) {
    730759                newPoly->coeff[x][y][z] = psAlloc(nT * sizeof(psF64));
    731760                newPoly->coeffErr[x][y][z] = psAlloc(nT * sizeof(psF64));
    732                 newPoly->mask[x][y][z] = (char *)psAlloc(nT * sizeof(char));
     761                newPoly->mask[x][y][z] = (psMaskType *)psAlloc(nT * sizeof(psMaskType));
    733762            }
    734763        }
    … …  
    749778}
    750779
    751 psF64 psPolynomial1DEval(const psPolynomial1D* poly, psF64 x)
     780psF64 psPolynomial1DEval(const psPolynomial1D* poly,
     781                         psF64 x)
    752782{
    753783    PS_ASSERT_POLY_NON_NULL(poly, NAN);
    … …  
    775805
    776806    tmp = psVectorAlloc(x->n, PS_TYPE_F64);
    777     for (psS32 i=0;i<x->n;i++) {
     807    for (unsigned int i=0;i<x->n;i++) {
    778808        tmp->data.F64[i] = psPolynomial1DEval(poly, x->data.F64[i]);
    779809    }
    … …  
    782812}
    783813
    784 psF64 psPolynomial2DEval(const psPolynomial2D* poly, psF64 x, psF64 y)
     814psF64 psPolynomial2DEval(const psPolynomial2D* poly,
     815                         psF64 x,
     816                         psF64 y)
    785817{
    786818    PS_ASSERT_POLY_NON_NULL(poly, NAN);
    … …  
    810842
    811843    psVector *tmp;
    812     psS32 vecLen=x->n;
     844    unsigned int vecLen=x->n;
    813845
    814846    // Determine the length of the output vector to by the minimum of the x,y vectors
    … …  
    821853
    822854    // Evaluate the polynomial at the specified points
    823     for (psS32 i=0; i<vecLen; i++) {
     855    for (unsigned int i=0; i<vecLen; i++) {
    824856        tmp->data.F64[i] = psPolynomial2DEval(poly,x->data.F64[i],y->data.F64[i]);
    825857    }
    … …  
    829861}
    830862
    831 psF64 psPolynomial3DEval(const psPolynomial3D* poly, psF64 x, psF64 y, psF64 z)
     863psF64 psPolynomial3DEval(const psPolynomial3D* poly,
     864                         psF64 x,
     865                         psF64 y,
     866                         psF64 z)
    832867{
    833868    PS_ASSERT_POLY_NON_NULL(poly, NAN);
    … …  
    860895
    861896    psVector *tmp;
    862     psS32 vecLen=x->n;
     897    unsigned int vecLen=x->n;
    863898
    864899    // Determine the length of output vector from min of the input vectors
    … …  
    874909
    875910    // Evaluate polynomial
    876     for (psS32 i = 0; i < vecLen; i++) {
     911    for (unsigned int i = 0; i < vecLen; i++) {
    877912        tmp->data.F64[i] = psPolynomial3DEval(poly,
    878913                                              x->data.F64[i],
    … …  
    885920}
    886921
    887 psF64 psPolynomial4DEval(const psPolynomial4D* poly, psF64 x, psF64 y, psF64 z, psF64 t)
     922psF64 psPolynomial4DEval(const psPolynomial4D* poly,
     923                         psF64 x,
     924                         psF64 y,
     925                         psF64 z,
     926                         psF64 t)
    888927{
    889928    PS_ASSERT_POLY_NON_NULL(poly, NAN);
    … …  
    918957
    919958    psVector *tmp;
    920     psS32 vecLen=x->n;
     959    unsigned int vecLen=x->n;
    921960
    922961    // Determine output vector size from min of input vectors
    … …  
    935974
    936975    // Evaluate polynomial
    937     for (psS32 i = 0; i < vecLen; i++) {
     976    for (unsigned int i = 0; i < vecLen; i++) {
    938977        tmp->data.F64[i] = psPolynomial4DEval(poly,
    939978                                              x->data.F64[i],

  • trunk/psLib/src/math/psPolynomial.h

    r4969 r5066  
    1111 *  @author GLG, MHPCC
    1212 *
    13  *  @version $Revision: 1.55 $ $Name: not supported by cvs2svn $
    14  *  @date $Date: 2005-09-08 00:02:48 $
     13 *  @version $Revision: 1.56 $ $Name: not supported by cvs2svn $
     14 *  @date $Date: 2005-09-19 19:53:13 $
    1515 *
    1616 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    5656    psF32 mean,                        ///< The mean of the Gaussian
    5757    psF32 sigma,                       ///< The sigma of the Gaussian
    58     psS32 Npts                         ///< The size of the vector
     58    unsigned int Npts                  ///< The size of the vector
    5959);
    6060
    … …  
    7373{
    7474    psPolynomialType type;             ///< Polynomial type
    75     int n;                             ///< Number of terms
     75    unsigned int n;                    ///< Number of terms
    7676    psF64 *coeff;                      ///< Coefficients
    7777    psF64 *coeffErr;                   ///< Error in coefficients
    78     char *mask;                        ///< Coefficient mask
     78    psMaskType *mask;                  ///< Coefficient mask
    7979}
    8080psPolynomial1D;
    … …  
    8484{
    8585    psPolynomialType type;             ///< Polynomial type
    86     int nX;                            ///< Number of terms in x
    87     int nY;                            ///< Number of terms in y
     86    unsigned int nX;                   ///< Number of terms in x
     87    unsigned int nY;                   ///< Number of terms in y
    8888    psF64 **coeff;                     ///< Coefficients
    8989    psF64 **coeffErr;                  ///< Error in coefficients
    90     char **mask;                       ///< Coefficients mask
     90    psMaskType **mask;                 ///< Coefficients mask
    9191}
    9292psPolynomial2D;
    … …  
    9696{
    9797    psPolynomialType type;             ///< Polynomial type
    98     int nX;                           ///< Number of terms in x
    99     int nY;                            ///< Number of terms in y
    100     int nZ;                           ///< Number of terms in z
     98    unsigned int nX;                   ///< Number of terms in x
     99    unsigned int nY;                   ///< Number of terms in y
     100    unsigned int nZ;                   ///< Number of terms in z
    101101    psF64 ***coeff;                    ///< Coefficients
    102102    psF64 ***coeffErr;                 ///< Error in coefficients
    103     char ***mask;                      ///< Coefficients mask
     103    psMaskType ***mask;                ///< Coefficients mask
    104104}
    105105psPolynomial3D;
    … …  
    109109{
    110110    psPolynomialType type;             ///< Polynomial type
    111     int nX;                            ///< Number of terms in x
    112     int nY;                            ///< Number of terms in y
    113     int nZ;                            ///< Number of terms in z
    114     int nT;                            ///< Number of terms in t
     111    unsigned int nX;                   ///< Number of terms in x
     112    unsigned int nY;                   ///< Number of terms in y
     113    unsigned int nZ;                   ///< Number of terms in z
     114    unsigned int nT;                   ///< Number of terms in t
    115115    psF64 ****coeff;                   ///< Coefficients
    116116    psF64 ****coeffErr;                ///< Error in coefficients
    117     char ****mask;                     ///< Coefficients mask
     117    psMaskType ****mask;               ///< Coefficients mask
    118118}
    119119psPolynomial4D;
    … …  
    125125 */
    126126psPolynomial1D* psPolynomial1DAlloc(
    127     int n,                             ///< Number of terms
     127    unsigned int n,                    ///< Number of terms
    128128    psPolynomialType type              ///< Polynomial Type
    129129);
    … …  
    134134 */
    135135psPolynomial2D* psPolynomial2DAlloc(
    136     int nX,                   ///< Number of terms in x
    137     int nY,                   ///< Number of terms in y
     136    unsigned int nX,                   ///< Number of terms in x
     137    unsigned int nY,                   ///< Number of terms in y
    138138    psPolynomialType type              ///< Polynomial Type
    139139);
    … …  
    144144 */
    145145psPolynomial3D* psPolynomial3DAlloc(
    146     int nX,                            ///< Number of terms in x
    147     int nY,                            ///< Number of terms in y
    148     int nZ,                            ///< Number of terms in z
     146    unsigned int nX,                   ///< Number of terms in x
     147    unsigned int nY,                   ///< Number of terms in y
     148    unsigned int nZ,                   ///< Number of terms in z
    149149    psPolynomialType type              ///< Polynomial Type
    150150);
    … …  
    155155 */
    156156psPolynomial4D* psPolynomial4DAlloc(
    157     int nX,                            ///< Number of terms in x
    158     int nY,                            ///< Number of terms in y
    159     int nZ,                            ///< Number of terms in z
    160     int nT,                            ///< Number of terms in t
     157    unsigned int nX,                   ///< Number of terms in x
     158    unsigned int nY,                   ///< Number of terms in y
     159    unsigned int nZ,                   ///< Number of terms in z
     160    unsigned int nT,                   ///< Number of terms in t
    161161    psPolynomialType type              ///< Polynomial Type
    162162);

  • trunk/psLib/src/math/psSpline.c

    r4991 r5066  
    77*  splines.
    88*
    9 *  @version $Revision: 1.123 $ $Name: not supported by cvs2svn $
    10 *  @date $Date: 2005-09-11 22:18:40 $
     9*  @version $Revision: 1.124 $ $Name: not supported by cvs2svn $
     10*  @date $Date: 2005-09-19 19:53:13 $
    1111*
    1212*
    … …  
    4343/*****************************************************************************/
    4444static void spline1DFree(psSpline1D *tmpSpline);
    45 static psS32 vectorBinDisectF32(psF32 *bins,psS32 numBins,psF32 x);
    46 static psS32 vectorBinDisectS32(psS32 *bins,psS32 numBins,psS32 x);
     45static unsigned int vectorBinDisectF32(psF32 *bins,psS32 numBins,psF32 x);
     46static unsigned int vectorBinDisectS32(psS32 *bins,psS32 numBins,psS32 x);
    4747
    4848/*****************************************************************************/
    … …  
    6969static void spline1DFree(psSpline1D *tmpSpline)
    7070{
    71     psS32 i;
     71    unsigned int i;
    7272
    7373    if (tmpSpline == NULL) {
    … …  
    101101static psF32 fullInterpolate1D##TYPE(ps##TYPE *domain, \
    102102                                     ps##TYPE *range, \
    103                                      psS32 n, \
     103                                     unsigned int n, \
    104104                                     ps##TYPE x) \
    105105{ \
    106106    \
    107     psS32 i; \
    108     psS32 m; \
     107    unsigned int i; \
     108    unsigned int m; \
    109109    static psVector *p = NULL; \
    110110    p = psVectorRecycle(p, n, PS_TYPE_##TYPE); \
    … …  
    169169static psF32 interpolate1DF32(psF32 *domain,
    170170                              psF32 *range,
    171                               psS32 n,
    172                               psS32 order,
     171                              unsigned int n,
     172                              unsigned int order,
    173173                              psF32 x)
    174174{
    … …  
    178178
    179179    psS32 binNum;
    180     psS32 numIntPoints = order+1;
     180    unsigned int numIntPoints = order+1;
    181181    psS32 origin;
    182182
    … …  
    229229            "---- calculateSecondDerivs() begin ----\n");
    230230
    231     psS32 i;
    232     psS32 k;
     231    unsigned int i;
     232    unsigned int k;
    233233    psF32 sig;
    234234    psF32 p;
    235     psS32 n = y->n;
     235    unsigned int n = y->n;
    236236    psF32 *u = (psF32 *) psAlloc(n * sizeof(psF32));
    237237    psF32 *derivs2 = (psF32 *) psAlloc(n * sizeof(psF32));
    … …  
    253253
    254254        psTrace(".psLib.dataManip.calculateSecondDerivs", 6,
    255                 "X[%d] is %f\n", i, X[i]);
     255                "X[%u] is %f\n", i, X[i]);
    256256        psTrace(".psLib.dataManip.calculateSecondDerivs", 6,
    257                 "Y[%d] is %f\n", i, Y[i]);
     257                "Y[%u] is %f\n", i, Y[i]);
    258258        psTrace(".psLib.dataManip.calculateSecondDerivs", 6,
    259                 "u[%d] is %f\n", i, u[i]);
     259                "u[%u] is %f\n", i, u[i]);
    260260    }
    261261
    … …  
    268268
    269269        psTrace(".psLib.dataManip.calculateSecondDerivs", 6,
    270                 "derivs2[%d] is %f\n", k, derivs2[k]);
     270                "derivs2[%u] is %f\n", k, derivs2[k]);
    271271    }
    272272
    … …  
    314314XXX: Assumes mySpline->knots is psF32.  Must add psU32 and psF64.
    315315 *****************************************************************************/
    316 psSpline1D *psVectorFitSpline1D(psSpline1D *mySpline,     ///< The spline which will be generated.
     316psSpline1D *psVectorFitSpline1D(psSpline1D *spline,     ///< The spline which will be generated.
    317317                                const psVector* x,        ///< Ordinates (or NULL to just use the indices)
    318318                                const psVector* y,        ///< Coordinates
    … …  
    321321    PS_ASSERT_VECTOR_NON_NULL(y, NULL);
    322322    PS_ASSERT_VECTOR_TYPE_F32_OR_F64(y, NULL);
    323     if (mySpline != NULL) {
    324         PS_ASSERT_VECTOR_TYPE(mySpline->knots, PS_TYPE_F32, NULL);
     323    if (spline != NULL) {
     324        PS_ASSERT_VECTOR_TYPE(spline->knots, PS_TYPE_F32, NULL);
    325325    }
    326326    psTrace(".psLib.dataManip.psVectorFitSpline1D", 4,
    327327            "---- psVectorFitSpline1D() begin ----\n");
    328     psS32 numSplines = (y->n)-1;
     328    unsigned int numSplines = (y->n)-1;
    329329    psF32 tmp;
    330330    psF32 H;
    331     psS32 i;
     331    unsigned int i;
    332332    psF32 slope;
    333333    psVector *x32 = NULL;
    … …  
    364364        This can not be implemented until SDR states what order spline should be
    365365        created.
    366         Should we error if mySpline is not NULL?
     366        Should we error if spline is not NULL?
    367367        Should we error if mySPline is not NULL?
    368368    */
    369     if (mySpline == NULL) {
    370         mySpline = psSpline1DAllocGeneric(x32, 3);
    371     }
    372     PS_ASSERT_PTR_NON_NULL(mySpline, NULL);
    373     PS_ASSERT_INT_NONNEGATIVE(mySpline->n, NULL);
    374 
    375     if (y32->n != (1 + mySpline->n)) {
     369    if (spline == NULL) {
     370        spline = psSpline1DAllocGeneric(x32, 3);
     371    }
     372    PS_ASSERT_PTR_NON_NULL(spline, NULL);
     373    PS_ASSERT_INT_NONNEGATIVE(spline->n, NULL);
     374
     375    if (y32->n != (1 + spline->n)) {
    376376        psError(PS_ERR_BAD_PARAMETER_SIZE, true,
    377                 "data size / spline size mismatch (%d %d)\n",
    378                 y32->n, mySpline->n);
     377                "data size / spline size mismatch (%u %u)\n",
     378                y32->n, spline->n);
    379379        return(NULL);
    380380    }
    … …  
    382382    // If these are linear splines, which means their polynomials will have
    383383    // two coefficients, then we do the simple calculation.
    384     if (2 == (mySpline->spline[0])->n) {
     384    if (2 == (spline->spline[0])->n) {
     385        for (i=0;i<spline->n;i++) {
     386            slope = (y32->data.F32[i+1] - y32->data.F32[i]) /
     387                    (spline->knots->data.F32[i+1] - spline->knots->data.F32[i]);
     388            (spline->spline[i])->coeff[0] = y32->data.F32[i] -
     389                                            (slope * spline->knots->data.F32[i]);
     390
     391            (spline->spline[i])->coeff[1] = slope;
     392            psTrace(".psLib.dataManip.psMinimize.psVectorFitSpline1D", 4,
     393                    "---- spline %u coeffs are (%f, %f)\n", i,
     394                    (spline->spline[i])->coeff[0],
     395                    (spline->spline[i])->coeff[1]);
     396        }
     397        psTrace(".psLib.dataManip.psMinimize.psVectorFitSpline1D", 4,
     398                "---- Exiting psVectorFitSpline1D()()\n");
     399        return((psSpline1D *) spline);
     400    }
     401
     402    // Check if these are cubic splines (n==4).  If not, psError.
     403    if (4 != (spline->spline[0])->n) {
     404        psError(PS_ERR_BAD_PARAMETER_SIZE, true,
     405                "Don't know how to generate %u-order splines.",
     406                (spline->spline[0])->n-1);
     407        return(NULL);
     408    }
     409
     410    // If we get here, then we know these are cubic splines.  We first
     411    // generate the second derivatives at each data point.
     412    spline->p_psDeriv2 = calculateSecondDerivs(x32, y32);
     413    for (i=0;i<y32->n;i++)
     414        psTrace(".psLib.dataManip.psVectorFitSpline1D", 6,
     415                "Second deriv[%u] is %f\n", i, spline->p_psDeriv2[i]);
     416
     417    // We generate the coefficients of the spline polynomials.  I can't
     418    // concisely explain how this code works.  See above function comments
     419    // and Numerical Recipes in C.
     420    for (i=0;i<numSplines;i++) {
     421        H = x32->data.F32[i+1] - x32->data.F32[i];
     422        psTrace(".psLib.dataManip.psVectorFitSpline1D", 4,
     423                "x data (%f - %f) (%f)\n",
     424                x32->data.F32[i],
     425                x32->data.F32[i+1], H);
     426        //
     427        // ******** Calculate 0-order term ********
     428        //
     429        // From (1)
     430        (spline->spline[i])->coeff[0] = (y32->data.F32[i] * x32->data.F32[i+1]/H);
     431        // From (2)
     432        ((spline->spline[i])->coeff[0])-= ((y32->data.F32[i+1] * x32->data.F32[i])/H);
     433        // From (3)
     434        tmp = (x32->data.F32[i+1] * x32->data.F32[i+1] * x32->data.F32[i+1]) / (H * H * H);
     435        tmp-= (x32->data.F32[i+1] / H);
     436        tmp*= (spline->p_psDeriv2)[i] * H * H / 6.0;
     437        ((spline->spline[i])->coeff[0])+= tmp;
     438        // From (4)
     439        tmp = -(x32->data.F32[i] * x32->data.F32[i] * x32->data.F32[i]) / (H * H * H);
     440        tmp+= (x32->data.F32[i] / H);
     441        tmp*= (spline->p_psDeriv2)[i+1] * H * H / 6.0;
     442        ((spline->spline[i])->coeff[0])+= tmp;
     443
     444        //
     445        // ******** Calculate 1-order term ********
     446        //
     447        // From (1)
     448        (spline->spline[i])->coeff[1] = -(y32->data.F32[i]) / H;
     449        // From (2)
     450        ((spline->spline[i])->coeff[1])+= (y32->data.F32[i+1] / H);
     451        // From (3)
     452        tmp = -3.0 * (x32->data.F32[i+1] * x32->data.F32[i+1]) / (H * H * H);
     453        tmp+= (1.0 / H);
     454        tmp*= ((spline->p_psDeriv2)[i]) * H * H / 6.0;
     455        ((spline->spline[i])->coeff[1])+= tmp;
     456        // From (4)
     457        tmp = 3.0 * (x32->data.F32[i] * x32->data.F32[i]) / (H * H * H);
     458        tmp-= (1.0 / H);
     459        tmp*= ((spline->p_psDeriv2)[i+1]) * H * H / 6.0;
     460        ((spline->spline[i])->coeff[1])+= tmp;
     461
     462        //
     463        // ******** Calculate 2-order term ********
     464        //
     465        // From (3)
     466        (spline->spline[i])->coeff[2] = ((spline->p_psDeriv2)[i]) * 3.0 * x32->data.F32[i+1] / (6.0 * H);
     467        // From (4)
     468        ((spline->spline[i])->coeff[2])-= (((spline->p_psDeriv2)[i+1]) * 3.0 * x32->data.F32[i] / (6.0 * H));
     469
     470        //
     471        // ******** Calculate 3-order term ********
     472        //
     473        // From (3)
     474        (spline->spline[i])->coeff[3] = -((spline->p_psDeriv2)[i]) / (6.0 * H);
     475        // From (4)
     476        ((spline->spline[i])->coeff[3])+=  ((spline->p_psDeriv2)[i+1]) / (6.0 * H);
     477
     478        psTrace(".psLib.dataManip.psVectorFitSpline1D", 6,
     479                "(spline->spline[%u])->coeff[0] is %f\n", i, (spline->spline[i])->coeff[0]);
     480        psTrace(".psLib.dataManip.psVectorFitSpline1D", 6,
     481                "(spline->spline[%u])->coeff[1] is %f\n", i, (spline->spline[i])->coeff[1]);
     482        psTrace(".psLib.dataManip.psVectorFitSpline1D", 6,
     483                "(spline->spline[%u])->coeff[2] is %f\n", i, (spline->spline[i])->coeff[2]);
     484        psTrace(".psLib.dataManip.psVectorFitSpline1D", 6,
     485                "(spline->spline[%u])->coeff[3] is %f\n", i, (spline->spline[i])->coeff[3]);
     486
     487    }
     488
     489    psTrace(".psLib.dataManip.psVectorFitSpline1D", 4,
     490            "---- psVectorFitSpline1D() end ----\n");
     491    return(spline);
     492}
     493
     494
     495
     496
     497
     498
     499
     500
     501
     502/*****************************************************************************
     503psVectorFitSpline1DNEW(): given a psSpline1D data structure and a set of x/y
     504 
     505xF32 and yF32 are internal psVectors which are used to hold the psF32 versions
     506of the input data, if necessary.  xPtr and yPtr are pointers to either xF32 or
     507the x argument.  All computation is done on xPtr and yPtr.  xF32 and yF32 will
     508simply be psFree() at the end.
     509 
     510XXX: nKnots makes no sense.  This number is always equal to the size of the x
     511an y vectors.
     512 
     513XXX: How do we specify the spline order?  For now, order=3.
     514 *****************************************************************************/
     515#define PS_XXX_SPLINE_ORDER 3
     516psSpline1D *psVectorFitSpline1DNEW(const psVector* x,        ///< Ordinates (or NULL to just use the indices)
     517                                   const psVector* y,        ///< Coordinates
     518                                   unsigned int nKnots)
     519{
     520    psTrace(".psLib.dataManip.psVectorFitSpline1D", 4, "---- psVectorFitSpline1DNEW() begin ----\n");
     521    PS_ASSERT_VECTOR_NON_NULL(y, NULL);
     522    PS_ASSERT_VECTOR_TYPE_F32_OR_F64(y, NULL);
     523    //    PS_ASSERT_INT_EQUAL(y->n, nKnots);
     524
     525    //
     526    // The following code ensures that xPtr points to a psF32 version of the
     527    // ordinate data.
     528    //
     529    psVector *xF32 = NULL;
     530    psVector *xPtr = NULL;
     531    if (x != NULL) {
     532        PS_ASSERT_VECTORS_SIZE_EQUAL(x, y, NULL);
     533        if (PS_TYPE_F64 == x->type.type) {
     534            xF32 = psVectorAlloc(y->n, PS_TYPE_F32);
     535            for (unsigned int i = 0 ; i < x->n ; i++) {
     536                xF32->data.F32[i] = (psF32) x->data.F64[i];
     537            }
     538            xPtr = xF32;
     539        } else if (PS_TYPE_F32 == x->type.type) {
     540            xPtr = (psVector *) x;
     541        } else {
     542            psError(PS_ERR_UNKNOWN, true, "psVector x is wrong type.\n");
     543            return(NULL);
     544        }
     545    } else {
     546        // If x==NULL, create an x32 vector with x values set to (0:n).
     547        xF32 = psVectorAlloc(y->n, PS_TYPE_F32);
     548        for (unsigned int i = 0 ; i < x->n ; i++) {
     549            xF32->data.F32[i] = (psF32) i;
     550        }
     551        xPtr = xF32;
     552    }
     553
     554    //
     555    // If y is of type psF64, then create a new vector yF32 and convert the
     556    // y elements.  Regardless of y's type, we create a yPtr which will be
     557    // used in the remainder of this function.
     558    //
     559    psVector *yF32 = NULL;
     560    psVector *yPtr = NULL;
     561    if (PS_TYPE_F64 == y->type.type) {
     562        yF32 = psVectorAlloc(y->n, PS_TYPE_F32);
     563        for (unsigned int i = 0 ; i < y->n ; i++) {
     564            yF32->data.F32[i] = (psF32) y->data.F64[i];
     565        }
     566        yPtr = yF32;
     567    } else {
     568        yPtr = (psVector *) y;
     569    }
     570
     571    psSpline1D *mySpline = psSpline1DAllocGeneric(xPtr, PS_XXX_SPLINE_ORDER);
     572
     573    unsigned int numSplines = nKnots - 1;
     574    psF32 tmp;
     575    psF32 H;
     576    unsigned int i;
     577    psF32 slope;
     578    // XXX: get rid of x32 and y32 (this is from old code)
     579    psVector *x32 = xPtr;
     580    psVector *y32 = yPtr;
     581
     582    // If these are linear splines, which means their polynomials will have
     583    // two coefficients, then we do the simple calculation.
     584    if (1 == PS_XXX_SPLINE_ORDER) {
    385585        for (i=0;i<mySpline->n;i++) {
    386586            slope = (y32->data.F32[i+1] - y32->data.F32[i]) /
    … …  
    391591            (mySpline->spline[i])->coeff[1] = slope;
    392592            psTrace(".psLib.dataManip.psMinimize.psVectorFitSpline1D", 4,
    393                     "---- mySpline %d coeffs are (%f, %f)\n", i,
     593                    "---- mySpline %u coeffs are (%f, %f)\n", i,
    394594                    (mySpline->spline[i])->coeff[0],
    395595                    (mySpline->spline[i])->coeff[1]);
    … …  
    400600    }
    401601
     602    //
    402603    // Check if these are cubic splines (n==4).  If not, psError.
    403     if (4 != (mySpline->spline[0])->n) {
     604    //
     605    if (3 != PS_XXX_SPLINE_ORDER) {
    404606        psError(PS_ERR_BAD_PARAMETER_SIZE, true,
    405                 "Don't know how to generate %d-order splines.",
    406                 (mySpline->spline[0])->n-1);
     607                "Don't know how to generate %d-order splines.", PS_XXX_SPLINE_ORDER);
    407608        return(NULL);
    408609    }
    409610
     611    //
    410612    // If we get here, then we know these are cubic splines.  We first
    411613    // generate the second derivatives at each data point.
     614    //
    412615    mySpline->p_psDeriv2 = calculateSecondDerivs(x32, y32);
    413616    for (i=0;i<y32->n;i++)
    414617        psTrace(".psLib.dataManip.psVectorFitSpline1D", 6,
    415                 "Second deriv[%d] is %f\n", i, mySpline->p_psDeriv2[i]);
    416 
     618                "Second deriv[%u] is %f\n", i, mySpline->p_psDeriv2[i]);
     619
     620    //
    417621    // We generate the coefficients of the spline polynomials.  I can't
    418622    // concisely explain how this code works.  See above function comments
    419623    // and Numerical Recipes in C.
     624    //
    420625    for (i=0;i<numSplines;i++) {
    421626        H = x32->data.F32[i+1] - x32->data.F32[i];
    … …  
    477682
    478683        psTrace(".psLib.dataManip.psVectorFitSpline1D", 6,
    479                 "(mySpline->spline[%d])->coeff[0] is %f\n", i, (mySpline->spline[i])->coeff[0]);
     684                "(mySpline->spline[%u])->coeff[0] is %f\n", i, (mySpline->spline[i])->coeff[0]);
    480685        psTrace(".psLib.dataManip.psVectorFitSpline1D", 6,
    481                 "(mySpline->spline[%d])->coeff[1] is %f\n", i, (mySpline->spline[i])->coeff[1]);
     686                "(mySpline->spline[%u])->coeff[1] is %f\n", i, (mySpline->spline[i])->coeff[1]);
    482687        psTrace(".psLib.dataManip.psVectorFitSpline1D", 6,
    483                 "(mySpline->spline[%d])->coeff[2] is %f\n", i, (mySpline->spline[i])->coeff[2]);
     688                "(mySpline->spline[%u])->coeff[2] is %f\n", i, (mySpline->spline[i])->coeff[2]);
    484689        psTrace(".psLib.dataManip.psVectorFitSpline1D", 6,
    485                 "(mySpline->spline[%d])->coeff[3] is %f\n", i, (mySpline->spline[i])->coeff[3]);
    486 
    487     }
    488 
    489     psTrace(".psLib.dataManip.psVectorFitSpline1D", 4,
    490             "---- psVectorFitSpline1D() end ----\n");
    491     return(mySpline);
    492 }
    493 
    494 
    495 
    496 
    497 
    498 
    499 
    500 
    501 
    502 /*****************************************************************************
    503 psVectorFitSpline1DNEW(): given a psSpline1D data structure and a set of x/y
    504  
    505 xF32 and yF32 are internal psVectors which are used to hold the psF32 versions
    506 of the input data, if necessary.  xPtr and yPtr are pointers to either xF32 or
    507 the x argument.  All computation is done on xPtr and yPtr.  xF32 and yF32 will
    508 simply be psFree() at the end.
    509  
    510 XXX: nKnots makes no sense.  This number is always equal to the size of the x
    511 an y vectors.
    512  
    513 XXX: How do we specify the spline order?  For now, order=3.
    514  *****************************************************************************/
    515 #define PS_XXX_SPLINE_ORDER 3
    516 psSpline1D *psVectorFitSpline1DNEW(const psVector* x,        ///< Ordinates (or NULL to just use the indices)
    517                                    const psVector* y,        ///< Coordinates
    518                                    int nKnots)
    519 {
    520     psTrace(".psLib.dataManip.psVectorFitSpline1D", 4, "---- psVectorFitSpline1DNEW() begin ----\n");
    521     PS_ASSERT_VECTOR_NON_NULL(y, NULL);
    522     PS_ASSERT_VECTOR_TYPE_F32_OR_F64(y, NULL);
    523     //    PS_ASSERT_INT_EQUAL(y->n, nKnots);
    524 
    525     //
    526     // The following code ensures that xPtr points to a psF32 version of the
    527     // ordinate data.
    528     //
    529     psVector *xF32 = NULL;
    530     psVector *xPtr = NULL;
    531     if (x != NULL) {
    532         PS_ASSERT_VECTORS_SIZE_EQUAL(x, y, NULL);
    533         if (PS_TYPE_F64 == x->type.type) {
    534             xF32 = psVectorAlloc(y->n, PS_TYPE_F32);
    535             for (psS32 i = 0 ; i < x->n ; i++) {
    536                 xF32->data.F32[i] = (psF32) x->data.F64[i];
    537             }
    538             xPtr = xF32;
    539         } else if (PS_TYPE_F32 == x->type.type) {
    540             xPtr = (psVector *) x;
    541         } else {
    542             psError(PS_ERR_UNKNOWN, true, "psVector x is wrong type.\n");
    543             return(NULL);
    544         }
    545     } else {
    546         // If x==NULL, create an x32 vector with x values set to (0:n).
    547         xF32 = psVectorAlloc(y->n, PS_TYPE_F32);
    548         for (psS32 i = 0 ; i < x->n ; i++) {
    549             xF32->data.F32[i] = (psF32) i;
    550         }
    551         xPtr = xF32;
    552     }
    553 
    554     //
    555     // If y is of type psF64, then create a new vector yF32 and convert the
    556     // y elements.  Regardless of y's type, we create a yPtr which will be
    557     // used in the remainder of this function.
    558     //
    559     psVector *yF32 = NULL;
    560     psVector *yPtr = NULL;
    561     if (PS_TYPE_F64 == y->type.type) {
    562         yF32 = psVectorAlloc(y->n, PS_TYPE_F32);
    563         for (psS32 i = 0 ; i < y->n ; i++) {
    564             yF32->data.F32[i] = (psF32) y->data.F64[i];
    565         }
    566         yPtr = yF32;
    567     } else {
    568         yPtr = (psVector *) y;
    569     }
    570 
    571     psSpline1D *mySpline = psSpline1DAllocGeneric(xPtr, PS_XXX_SPLINE_ORDER);
    572 
    573     psS32 numSplines = nKnots - 1;
    574     psF32 tmp;
    575     psF32 H;
    576     psS32 i;
    577     psF32 slope;
    578     // XXX: get rid of x32 and y32 (this is from old code)
    579     psVector *x32 = xPtr;
    580     psVector *y32 = yPtr;
    581 
    582     // If these are linear splines, which means their polynomials will have
    583     // two coefficients, then we do the simple calculation.
    584     if (1 == PS_XXX_SPLINE_ORDER) {
    585         for (i=0;i<mySpline->n;i++) {
    586             slope = (y32->data.F32[i+1] - y32->data.F32[i]) /
    587                     (mySpline->knots->data.F32[i+1] - mySpline->knots->data.F32[i]);
    588             (mySpline->spline[i])->coeff[0] = y32->data.F32[i] -
    589                                               (slope * mySpline->knots->data.F32[i]);
    590 
    591             (mySpline->spline[i])->coeff[1] = slope;
    592             psTrace(".psLib.dataManip.psMinimize.psVectorFitSpline1D", 4,
    593                     "---- mySpline %d coeffs are (%f, %f)\n", i,
    594                     (mySpline->spline[i])->coeff[0],
    595                     (mySpline->spline[i])->coeff[1]);
    596         }
    597         psTrace(".psLib.dataManip.psMinimize.psVectorFitSpline1D", 4,
    598                 "---- Exiting psVectorFitSpline1D()()\n");
    599         return((psSpline1D *) mySpline);
    600     }
    601 
    602     //
    603     // Check if these are cubic splines (n==4).  If not, psError.
    604     //
    605     if (3 != PS_XXX_SPLINE_ORDER) {
    606         psError(PS_ERR_BAD_PARAMETER_SIZE, true,
    607                 "Don't know how to generate %d-order splines.", PS_XXX_SPLINE_ORDER);
    608         return(NULL);
    609     }
    610 
    611     //
    612     // If we get here, then we know these are cubic splines.  We first
    613     // generate the second derivatives at each data point.
    614     //
    615     mySpline->p_psDeriv2 = calculateSecondDerivs(x32, y32);
    616     for (i=0;i<y32->n;i++)
    617         psTrace(".psLib.dataManip.psVectorFitSpline1D", 6,
    618                 "Second deriv[%d] is %f\n", i, mySpline->p_psDeriv2[i]);
    619 
    620     //
    621     // We generate the coefficients of the spline polynomials.  I can't
    622     // concisely explain how this code works.  See above function comments
    623     // and Numerical Recipes in C.
    624     //
    625     for (i=0;i<numSplines;i++) {
    626         H = x32->data.F32[i+1] - x32->data.F32[i];
    627         psTrace(".psLib.dataManip.psVectorFitSpline1D", 4,
    628                 "x data (%f - %f) (%f)\n",
    629                 x32->data.F32[i],
    630                 x32->data.F32[i+1], H);
    631         //
    632         // ******** Calculate 0-order term ********
    633         //
    634         // From (1)
    635         (mySpline->spline[i])->coeff[0] = (y32->data.F32[i] * x32->data.F32[i+1]/H);
    636         // From (2)
    637         ((mySpline->spline[i])->coeff[0])-= ((y32->data.F32[i+1] * x32->data.F32[i])/H);
    638         // From (3)
    639         tmp = (x32->data.F32[i+1] * x32->data.F32[i+1] * x32->data.F32[i+1]) / (H * H * H);
    640         tmp-= (x32->data.F32[i+1] / H);
    641         tmp*= (mySpline->p_psDeriv2)[i] * H * H / 6.0;
    642         ((mySpline->spline[i])->coeff[0])+= tmp;
    643         // From (4)
    644         tmp = -(x32->data.F32[i] * x32->data.F32[i] * x32->data.F32[i]) / (H * H * H);
    645         tmp+= (x32->data.F32[i] / H);
    646         tmp*= (mySpline->p_psDeriv2)[i+1] * H * H / 6.0;
    647         ((mySpline->spline[i])->coeff[0])+= tmp;
    648 
    649         //
    650         // ******** Calculate 1-order term ********
    651         //
    652         // From (1)
    653         (mySpline->spline[i])->coeff[1] = -(y32->data.F32[i]) / H;
    654         // From (2)
    655         ((mySpline->spline[i])->coeff[1])+= (y32->data.F32[i+1] / H);
    656         // From (3)
    657         tmp = -3.0 * (x32->data.F32[i+1] * x32->data.F32[i+1]) / (H * H * H);
    658         tmp+= (1.0 / H);
    659         tmp*= ((mySpline->p_psDeriv2)[i]) * H * H / 6.0;
    660         ((mySpline->spline[i])->coeff[1])+= tmp;
    661         // From (4)
    662         tmp = 3.0 * (x32->data.F32[i] * x32->data.F32[i]) / (H * H * H);
    663         tmp-= (1.0 / H);
    664         tmp*= ((mySpline->p_psDeriv2)[i+1]) * H * H / 6.0;
    665         ((mySpline->spline[i])->coeff[1])+= tmp;
    666 
    667         //
    668         // ******** Calculate 2-order term ********
    669         //
    670         // From (3)
    671         (mySpline->spline[i])->coeff[2] = ((mySpline->p_psDeriv2)[i]) * 3.0 * x32->data.F32[i+1] / (6.0 * H);
    672         // From (4)
    673         ((mySpline->spline[i])->coeff[2])-= (((mySpline->p_psDeriv2)[i+1]) * 3.0 * x32->data.F32[i] / (6.0 * H));
    674 
    675         //
    676         // ******** Calculate 3-order term ********
    677         //
    678         // From (3)
    679         (mySpline->spline[i])->coeff[3] = -((mySpline->p_psDeriv2)[i]) / (6.0 * H);
    680         // From (4)
    681         ((mySpline->spline[i])->coeff[3])+=  ((mySpline->p_psDeriv2)[i+1]) / (6.0 * H);
    682 
    683         psTrace(".psLib.dataManip.psVectorFitSpline1D", 6,
    684                 "(mySpline->spline[%d])->coeff[0] is %f\n", i, (mySpline->spline[i])->coeff[0]);
    685         psTrace(".psLib.dataManip.psVectorFitSpline1D", 6,
    686                 "(mySpline->spline[%d])->coeff[1] is %f\n", i, (mySpline->spline[i])->coeff[1]);
    687         psTrace(".psLib.dataManip.psVectorFitSpline1D", 6,
    688                 "(mySpline->spline[%d])->coeff[2] is %f\n", i, (mySpline->spline[i])->coeff[2]);
    689         psTrace(".psLib.dataManip.psVectorFitSpline1D", 6,
    690                 "(mySpline->spline[%d])->coeff[3] is %f\n", i, (mySpline->spline[i])->coeff[3]);
     690                "(mySpline->spline[%u])->coeff[3] is %f\n", i, (mySpline->spline[i])->coeff[3]);
    691691
    692692    }
    … …  
    726726XXX: What should be the default type for knots be?  psF32 is assumed.
    727727 *****************************************************************************/
    728 psSpline1D *psSpline1DAlloc(int numSplines,
    729                             int order,
     728psSpline1D *psSpline1DAlloc(unsigned int numSplines,
     729                            unsigned int order,
    730730                            float min,
    731731                            float max)
    … …  
    743743    //
    744744    tmpSpline->spline = (psPolynomial1D **) psAlloc(numSplines * sizeof(psPolynomial1D *));
    745     for (psS32 i=0;i<numSplines;i++) {
     745    for (unsigned int i=0;i<numSplines;i++) {
    746746        (tmpSpline->spline)[i] = psPolynomial1DAlloc(order+1, PS_POLYNOMIAL_ORD);
    747747    }
    … …  
    756756    psF32 width = (max - min) / ((psF32) numSplines);
    757757    tmpSpline->knots->data.F32[0] = min;
    758     for (psS32 i=1;i<numSplines;i++) {
     758    for (unsigned int i=1;i<numSplines;i++) {
    759759        tmpSpline->knots->data.F32[i] = min + (width * (psF32) i);
    760760    }
    … …  
    773773
    774774    if (in->type.type == PS_TYPE_F32) {
    775         for (psS32 i = 0 ; i < in->n ; i++) {
     775        for (unsigned int i = 0 ; i < in->n ; i++) {
    776776            out->data.F32[i] = in->data.F32[i];
    777777        }
    778778    } else if (in->type.type == PS_TYPE_F64) {
    779         for (psS32 i = 0 ; i < in->n ; i++) {
     779        for (unsigned int i = 0 ; i < in->n ; i++) {
    780780            out->data.F64[i] = in->data.F64[i];
    781781        }
    … …  
    791791 *****************************************************************************/
    792792psSpline1D *psSpline1DAllocGeneric(const psVector *bounds,
    793                                    int order)
     793                                   unsigned int order)
    794794{
    795795    PS_ASSERT_VECTOR_NON_NULL(bounds, NULL);
    … …  
    799799
    800800    psSpline1D *tmpSpline = (psSpline1D *) psAlloc(sizeof(psSpline1D));
    801     psS32 numSplines = bounds->n - 1;
     801    unsigned int numSplines = bounds->n - 1;
    802802    tmpSpline->n = numSplines;
    803803
    … …  
    807807    //
    808808    tmpSpline->spline = (psPolynomial1D **) psAlloc(numSplines * sizeof(psPolynomial1D *));
    809     for (psS32 i=0;i<numSplines;i++) {
     809    for (unsigned int i=0;i<numSplines;i++) {
    810810        (tmpSpline->spline)[i] = psPolynomial1DAlloc(order+1, PS_POLYNOMIAL_ORD);
    811811    }
    … …  
    818818    // XXX:Ensure that the knots are monotonic.
    819819    //
    820     for (psS32 i=0;i<bounds->n-1;i++) {
     820    for (unsigned int i=0;i<bounds->n-1;i++) {
    821821        if (FLT_EPSILON >= fabs(bounds->data.F32[i+1]-bounds->data.F32[i])) {
    822822            psError(PS_ERR_UNKNOWN, true, "data points must be distinct ([%d] %f %f)\n", i, bounds->data.F32[i], bounds->data.F32[i+1]);
    … …  
    839839 *****************************************************************************/
    840840#define FUNC_MACRO_VECTOR_BIN_DISECT(TYPE) \
    841 static psS32 vectorBinDisect##TYPE(ps##TYPE *bins, \
    842                                    psS32 numBins, \
    843                                    ps##TYPE x) \
     841static unsigned int vectorBinDisect##TYPE(ps##TYPE *bins, \
     842        psS32 numBins, \
     843        ps##TYPE x) \
    844844{ \
    845845    psS32 min; \
    … …  
    906906XXX: Assert that the psVector and psScalar have the same type.
    907907 *****************************************************************************/
    908 psS32 p_psVectorBinDisect(psVector *bins,
    909                           psScalar *x)
     908unsigned int p_psVectorBinDisect(psVector *bins,
     909                                 psScalar *x)
    910910{
    911911    PS_ASSERT_VECTOR_NON_NULL(bins, -4);
    … …  
    972972psScalar *p_psVectorInterpolate(psVector *domain,
    973973                                psVector *range,
    974                                 int order,
     974                                unsigned int order,
    975975                                psScalar *x)
    976976{

  • trunk/psLib/src/math/psSpline.h

    r4991 r5066  
    1010 *  @author GLG, MHPCC
    1111 *
    12  *  @version $Revision: 1.55 $ $Name: not supported by cvs2svn $
    13  *  @date $Date: 2005-09-11 22:18:40 $
     12 *  @version $Revision: 1.56 $ $Name: not supported by cvs2svn $
     13 *  @date $Date: 2005-09-19 19:53:13 $
    1414 *
    1515 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    3535typedef struct
    3636{
    37     int n;                             ///< The number of spline pieces
     37    unsigned int n;                    ///< The number of spline pieces
    3838    psPolynomial1D **spline;           ///< An array of n pointers to the spline polynomials
    3939    psVector *knots;                   ///< The boundaries between each spline piece.  Size is n+1.
    … …  
    5151 */
    5252psSpline1D *psSpline1DAlloc(
    53     int n,                             ///< Number of spline polynomials
    54     int order,                         ///< Order of spline polynomials
     53    unsigned int n,                    ///< Number of spline polynomials
     54    unsigned int order,                ///< Order of spline polynomials
    5555    float min,                         ///< Lower boundary value of spline polynomials
    5656    float max                          ///< Upper boundary value of spline polynomials
    … …  
    6565psSpline1D *psSpline1DAllocGeneric(
    6666    const psVector *bounds,            ///< Bounds for spline polynomials
    67     int order                          ///< Order of spline polynomials
     67    unsigned int order                 ///< Order of spline polynomials
    6868);
    6969
    … …  
    9191 *  @return psS32    corresponding index number of specified value
    9292 */
    93 psS32 p_psVectorBinDisect(
     93unsigned int p_psVectorBinDisect(
    9494    psVector *bins,                    ///< Array of non-decreasing values
    9595    psScalar *x                        ///< Target value to find
    … …  
    104104    psVector *domain,                  ///< Domain (x coords) for interpolation
    105105    psVector *range,                   ///< Range (y coords) for interpolation
    106     int order,                         ///< Order of interpolation function
     106    unsigned int order,                ///< Order of interpolation function
    107107    psScalar *x                        ///< Location at which to evaluate
    108108);
    … …  
    126126 */
    127127psSpline1D *psVectorFitSpline1D(
    128     psSpline1D *mySpline,              ///< The spline which will be generated.
     128    psSpline1D *spline,                ///< The spline which will be generated.
    129129    const psVector* x,                 ///< Ordinates (or NULL to just use the indices)
    130130    const psVector* y,                 ///< Coordinates
    … …  
    135135    const psVector* x,                 ///< Ordinates (or NULL to just use the indices)
    136136    const psVector* y,                 ///< Coordinates
    137     int nKnots
     137    unsigned int nKnots                ///< Number of Knots
    138138);
    139139
Note: See TracChangeset for help on using the changeset viewer.