IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Sep 29, 2005, 9:40:59 AM (21 years ago)
Author:
gusciora
Message:

....

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r5155 r5179  
    77*  splines.
    88*
    9 *  @version $Revision: 1.129 $ $Name: not supported by cvs2svn $
    10 *  @date $Date: 2005-09-27 23:16:59 $
     9*  @version $Revision: 1.130 $ $Name: not supported by cvs2svn $
     10*  @date $Date: 2005-09-29 19:40:59 $
    1111*
    1212*  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    491491    PS_ASSERT_VECTOR_NON_NULL(y, NULL);
    492492    PS_ASSERT_VECTOR_TYPE_F32_OR_F64(y, NULL);
    493     PS_ASSERT_VECTOR_NON_NULL(x, NULL);
    494     PS_ASSERT_VECTOR_TYPE_F32_OR_F64(x, NULL);
    495     PS_ASSERT_VECTORS_SIZE_EQUAL(x, y, NULL);
    496493    psS32 numSplines = (y->n)-1;
    497 
    498494    psTrace(__func__, 5, "numSplines is %d\n", numSplines);
    499495    //
    500496    // The following code ensures that xPtr and yPtr points to a psF32 psVector.
    501497    //
    502     psVector *xF32 = NULL;
    503498    psVector *xPtr = NULL;
    504     if (PS_TYPE_F64 == x->type.type) {
    505         xF32 = psVectorAlloc(y->n, PS_TYPE_F32);
    506         for (psS32 i = 0 ; i < x->n ; i++) {
    507             xF32->data.F32[i] = (psF32) x->data.F64[i];
     499    if (x != NULL) {
     500        // Convert x to F32 if necessary.
     501        PS_ASSERT_VECTOR_NON_NULL(x, NULL);
     502        PS_ASSERT_VECTOR_TYPE_F32_OR_F64(x, NULL);
     503        PS_ASSERT_VECTORS_SIZE_EQUAL(x, y, NULL);
     504        if (PS_TYPE_F64 == x->type.type) {
     505            xPtr = psVectorAlloc(y->n, PS_TYPE_F32);
     506            for (psS32 i = 0 ; i < x->n ; i++) {
     507                xPtr->data.F32[i] = (psF32) x->data.F64[i];
     508            }
     509        } else if (PS_TYPE_F32 == x->type.type) {
     510            xPtr = (psVector *) x;
    508511        }
    509         xPtr = xF32;
    510     } else if (PS_TYPE_F32 == x->type.type) {
    511         xPtr = (psVector *) x;
    512     }
    513 
    514     psVector *yF32 = NULL;
     512    } else {
     513        // Allocate an index vector for x
     514        xPtr = psVectorCreate(NULL, 0.0, (psF64) y->n, 1.0, PS_TYPE_F32);
     515        printf("x, y is (%d, %d)\n", (psS32) xPtr->n, (psS32) y->n);
     516        xPtr = xPtr;
     517    }
     518
    515519    psVector *yPtr = NULL;
     520    // Convert y to F32 if necessary.
    516521    if (PS_TYPE_F64 == y->type.type) {
    517         yF32 = psVectorAlloc(y->n, PS_TYPE_F32);
     522        yPtr = psVectorAlloc(y->n, PS_TYPE_F32);
    518523        for (psS32 i = 0 ; i < y->n ; i++) {
    519             yF32->data.F32[i] = (psF32) y->data.F64[i];
     524            yPtr->data.F32[i] = (psF32) y->data.F64[i];
    520525        }
    521         yPtr = yF32;
    522526    } else {
    523527        yPtr = (psVector *) y;
    … …  
    531535    spline->spline = (psPolynomial1D **) psAlloc(numSplines * sizeof(psPolynomial1D *));
    532536    for (psS32 i=0;i<numSplines;i++) {
    533         (spline->spline)[i] = psPolynomial1DAlloc(3, PS_POLYNOMIAL_ORD);
    534     }
    535     //
    536     // XXX: Ensure that the knots are distinct, and monotonic.
    537     // XXX: Use a vector dup macro, or function here.
    538     //
    539     spline->knots = psVectorAlloc(x->n, PS_TYPE_F32);
    540     for (psS32 i = 0 ; i < x->n ; i++) {
    541         spline->knots->data.F32[i] = xPtr->data.F32[i];
    542     }
     537        spline->spline[i] = psPolynomial1DAlloc(3, PS_POLYNOMIAL_ORD);
     538    }
     539    spline->knots = psVectorCopy(NULL, xPtr, PS_TYPE_F32);
    543540    //
    544541    // Generate the second derivatives at each data point.
    … …  
    554551        psF32 H = xPtr->data.F32[i+1] - xPtr->data.F32[i];
    555552        if (fabs(H) <= FLT_EPSILON) {
    556             printf("XXX: Generate error: x data points are not distinct (%d %d).\n", i, i+1);
     553            psError(PS_ERR_UNKNOWN, false, "x data points are not distinct (%d %d).\n", i, i+1);
    557554        }
    558555        psTrace(__func__, 6, "x data (%f - %f) (%f)\n", xPtr->data.F32[i], xPtr->data.F32[i+1], H);
    559 
    560556        //
    561557        // ******** Calculate 0-order term ********
    562558        //
    563559        // From (1)
    564         (spline->spline[i])->coeff[0] = (yPtr->data.F32[i] * xPtr->data.F32[i+1]/H);
     560        spline->spline[i]->coeff[0] = yPtr->data.F32[i] * xPtr->data.F32[i+1]/H;
    565561        // From (2)
    566         ((spline->spline[i])->coeff[0])-= ((yPtr->data.F32[i+1] * xPtr->data.F32[i])/H);
     562        spline->spline[i]->coeff[0]-= (yPtr->data.F32[i+1] * xPtr->data.F32[i])/H;
    567563        // From (3)
    568564        psF32 tmp = (xPtr->data.F32[i+1] * xPtr->data.F32[i+1] * xPtr->data.F32[i+1]) / (H * H * H);
    569         tmp-= (xPtr->data.F32[i+1] / H);
    570         tmp*= (spline->p_psDeriv2)[i] * H * H / 6.0;
    571         ((spline->spline[i])->coeff[0])+= tmp;
     565        tmp-= xPtr->data.F32[i+1] / H;
     566        tmp*= spline->p_psDeriv2[i] * H * H / 6.0;
     567        spline->spline[i]->coeff[0]+= tmp;
    572568        // From (4)
    573569        tmp = -(xPtr->data.F32[i] * xPtr->data.F32[i] * xPtr->data.F32[i]) / (H * H * H);
    574         tmp+= (xPtr->data.F32[i] / H);
    575         tmp*= (spline->p_psDeriv2)[i+1] * H * H / 6.0;
    576         ((spline->spline[i])->coeff[0])+= tmp;
     570        tmp+= xPtr->data.F32[i] / H;
     571        tmp*= spline->p_psDeriv2[i+1] * H * H / 6.0;
     572        (spline->spline[i]->coeff[0])+= tmp;
    577573
    578574        //
    … …  
    580576        //
    581577        // From (1)
    582         (spline->spline[i])->coeff[1] = -(yPtr->data.F32[i]) / H;
     578        spline->spline[i]->coeff[1] = -(yPtr->data.F32[i]) / H;
    583579        // From (2)
    584         ((spline->spline[i])->coeff[1])+= (yPtr->data.F32[i+1] / H);
     580        spline->spline[i]->coeff[1]+= yPtr->data.F32[i+1] / H;
    585581        // From (3)
    586         tmp = -3.0 * (xPtr->data.F32[i+1] * xPtr->data.F32[i+1]) / (H * H * H);
     582        tmp = -3.0 * xPtr->data.F32[i+1] * xPtr->data.F32[i+1] / (H * H * H);
    587583        tmp+= (1.0 / H);
    588         tmp*= ((spline->p_psDeriv2)[i]) * H * H / 6.0;
    589         ((spline->spline[i])->coeff[1])+= tmp;
     584        tmp*= spline->p_psDeriv2[i] * H * H / 6.0;
     585        spline->spline[i]->coeff[1]+= tmp;
    590586        // From (4)
    591         tmp = 3.0 * (xPtr->data.F32[i] * xPtr->data.F32[i]) / (H * H * H);
    592         tmp-= (1.0 / H);
    593         tmp*= ((spline->p_psDeriv2)[i+1]) * H * H / 6.0;
    594         ((spline->spline[i])->coeff[1])+= tmp;
     587        tmp = 3.0 * xPtr->data.F32[i] * xPtr->data.F32[i] / (H * H * H);
     588        tmp-= 1.0 / H;
     589        tmp*= spline->p_psDeriv2[i+1] * H * H / 6.0;
     590        spline->spline[i]->coeff[1]+= tmp;
    595591
    596592        //
    … …  
    598594        //
    599595        // From (3)
    600         (spline->spline[i])->coeff[2] = ((spline->p_psDeriv2)[i]) * 3.0 * xPtr->data.F32[i+1] / (6.0 * H);
     596        spline->spline[i]->coeff[2] = ((spline->p_psDeriv2)[i]) * 3.0 * xPtr->data.F32[i+1] / (6.0 * H);
    601597        // From (4)
    602         ((spline->spline[i])->coeff[2])-= (((spline->p_psDeriv2)[i+1]) * 3.0 * xPtr->data.F32[i] / (6.0 * H));
     598        (spline->spline[i]->coeff[2])-= (((spline->p_psDeriv2)[i+1]) * 3.0 * xPtr->data.F32[i] / (6.0 * H));
    603599
    604600        //
    … …  
    606602        //
    607603        // From (3)
    608         (spline->spline[i])->coeff[3] = -((spline->p_psDeriv2)[i]) / (6.0 * H);
     604        spline->spline[i]->coeff[3] = -((spline->p_psDeriv2)[i]) / (6.0 * H);
    609605        // From (4)
    610         ((spline->spline[i])->coeff[3])+=  ((spline->p_psDeriv2)[i+1]) / (6.0 * H);
    611 
    612         psTrace(__func__, 6, "(spline->spline[%u])->coeff[0] is %f\n", i, (spline->spline[i])->coeff[0]);
    613         psTrace(__func__, 6, "(spline->spline[%u])->coeff[1] is %f\n", i, (spline->spline[i])->coeff[1]);
    614         psTrace(__func__, 6, "(spline->spline[%u])->coeff[2] is %f\n", i, (spline->spline[i])->coeff[2]);
    615         psTrace(__func__, 6, "(spline->spline[%u])->coeff[3] is %f\n", i, (spline->spline[i])->coeff[3]);
    616     }
    617 
    618     if (PS_TYPE_F64 == x->type.type) {
    619         psFree(xF32);
     606        (spline->spline[i]->coeff[3])+=  ((spline->p_psDeriv2)[i+1]) / (6.0 * H);
     607
     608        psTrace(__func__, 6, "(spline->spline[%u])->coeff[0] is %f\n", i, spline->spline[i]->coeff[0]);
     609        psTrace(__func__, 6, "(spline->spline[%u])->coeff[1] is %f\n", i, spline->spline[i]->coeff[1]);
     610        psTrace(__func__, 6, "(spline->spline[%u])->coeff[2] is %f\n", i, spline->spline[i]->coeff[2]);
     611        psTrace(__func__, 6, "(spline->spline[%u])->coeff[3] is %f\n", i, spline->spline[i]->coeff[3]);
     612    }
     613
     614    if ((x == NULL) || (PS_TYPE_F64 == x->type.type)) {
     615        psFree(xPtr);
    620616    }
    621617    if (PS_TYPE_F64 == y->type.type) {
    622         psFree(yF32);
     618        psFree(yPtr);
    623619    }
    624620    psTrace(__func__, 3, "---- %s() end ----\n", __func__);
    … …  
    637633     the spline fit functions require F32 and F64.
    638634 
    639 XXX: This only works if spline0>knots if psF32.  Must we add support for psU32 and
     635XXX: This only works if spline->knots if psF32.  Must we add support for psU32 and
    640636psF64?
    641637 *****************************************************************************/
Note: See TracChangeset for help on using the changeset viewer.