IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

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)

File:
1 edited

Legend:

Unmodified
Added
Removed

  • 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{
Note: See TracChangeset for help on using the changeset viewer.