IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Oct 14, 2021, 9:51:58 AM (5 years ago)
Author:
eugene
Message:

redo 1D chebyshevs to used macros for each order (up to 9)

File:
1 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/ipp-dev-20210817/psLib/src/math/psPolynomial.c

    r41831 r41839  
    257257}
    258258
     259# define CHEB_EVAL_0(OUT,IN) {OUT = 1.0;}
     260# define CHEB_EVAL_1(OUT,IN) {                       OUT = IN; }
     261# define CHEB_EVAL_2(OUT,IN) {psF64 X2 = PS_SQR(IN); OUT = 2.0*X2 - 1.0; }
     262# define CHEB_EVAL_3(OUT,IN) {psF64 X2 = PS_SQR(IN); OUT = IN*(4.0*X2 - 3.0); }
     263# define CHEB_EVAL_4(OUT,IN) {psF64 X2 = PS_SQR(IN); OUT = X2*(8.0*X2 - 8.0) + 1.0; }
     264# define CHEB_EVAL_5(OUT,IN) {psF64 X2 = PS_SQR(IN); OUT = IN *(X2*(16.0*X2 - 20.0) + 5.0); }
     265# define CHEB_EVAL_6(OUT,IN) {psF64 X2 = PS_SQR(IN); OUT = X2*(X2*(32.0*X2 - 48.0) + 18.0) - 1.0; }
     266# define CHEB_EVAL_7(OUT,IN) {psF64 X2 = PS_SQR(IN); OUT = IN *(X2*(X2*(64.0*X2 - 112.0) + 56.0) - 7.0); }
     267# define CHEB_EVAL_8(OUT,IN) {psF64 X2 = PS_SQR(IN); OUT = X2*(X2*(X2*(128.0*X2 - 256.0) + 160.0) - 32.0) + 1.0; }
     268# define CHEB_EVAL_9(OUT,IN) {psF64 X2 = PS_SQR(IN); OUT = IN *(X2*(X2*(X2*(256.0*X2 - 576.0) + 432.0) - 129.0) + 9.0); }
     269
    259270/** This function generates a vector containing the values of a Chebyshev polynomial of
    260271    the given order evaluated at the coordinates given by the input vector, i.e., this
    … …  
    274285    switch (order) {
    275286      case 0:
    276         for (int i = 0; i < vec->n; i++) {                                                   out->data.F64[i] = 1.0; } break;
     287        for (int i = 0; i < vec->n; i++) { CHEB_EVAL_0(out->data.F64[i], vec->data.F64[i]); } break;
    277288      case 1:
    278         for (int i = 0; i < vec->n; i++) { psF64 x = vec->data.F64[i];                       out->data.F64[i] = x; } break;
     289        for (int i = 0; i < vec->n; i++) { CHEB_EVAL_1(out->data.F64[i], vec->data.F64[i]); } break;
    279290      case 2:
    280         for (int i = 0; i < vec->n; i++) { psF64 x = vec->data.F64[i]; psF64 x2 = PS_SQR(x); out->data.F64[i] = 2.0*x2 - 1.0; } break;
     291        for (int i = 0; i < vec->n; i++) { CHEB_EVAL_2(out->data.F64[i], vec->data.F64[i]); } break;
    281292      case 3:
    282         for (int i = 0; i < vec->n; i++) { psF64 x = vec->data.F64[i]; psF64 x2 = PS_SQR(x); out->data.F64[i] = x*(4.0*x2 - 3.0); } break;
     293        for (int i = 0; i < vec->n; i++) { CHEB_EVAL_3(out->data.F64[i], vec->data.F64[i]); } break;
    283294      case 4:
    284         for (int i = 0; i < vec->n; i++) { psF64 x = vec->data.F64[i]; psF64 x2 = PS_SQR(x); out->data.F64[i] = x2*(8.0*x2 - 8.0) + 1.0; } break;
     295        for (int i = 0; i < vec->n; i++) { CHEB_EVAL_4(out->data.F64[i], vec->data.F64[i]); } break;
    285296      case 5:
    286         for (int i = 0; i < vec->n; i++) { psF64 x = vec->data.F64[i]; psF64 x2 = PS_SQR(x); out->data.F64[i] = x *(x2*(16.0*x2 - 20.0) + 5.0); } break;
     297        for (int i = 0; i < vec->n; i++) { CHEB_EVAL_5(out->data.F64[i], vec->data.F64[i]); } break;
    287298      case 6:
    288         for (int i = 0; i < vec->n; i++) { psF64 x = vec->data.F64[i]; psF64 x2 = PS_SQR(x); out->data.F64[i] = x2*(x2*(32.0*x2 - 48.0) + 18.0) - 1.0; } break;
     299        for (int i = 0; i < vec->n; i++) { CHEB_EVAL_6(out->data.F64[i], vec->data.F64[i]); } break;
    289300      case 7:
    290         for (int i = 0; i < vec->n; i++) { psF64 x = vec->data.F64[i]; psF64 x2 = PS_SQR(x); out->data.F64[i] = x *(x2*(x2*(64.0*x2 - 112.0) + 56.0) - 7.0); } break;
     301        for (int i = 0; i < vec->n; i++) { CHEB_EVAL_7(out->data.F64[i], vec->data.F64[i]); } break;
    291302      case 8:
    292         for (int i = 0; i < vec->n; i++) { psF64 x = vec->data.F64[i]; psF64 x2 = PS_SQR(x); out->data.F64[i] = x2*(x2*(x2*(128.0*x2 - 256.0) + 160.0) - 32.0) + 1.0; } break;
     303        for (int i = 0; i < vec->n; i++) { CHEB_EVAL_8(out->data.F64[i], vec->data.F64[i]); } break;
    293304      case 9:
    294         for (int i = 0; i < vec->n; i++) { psF64 x = vec->data.F64[i]; psF64 x2 = PS_SQR(x); out->data.F64[i] = x *(x2*(x2*(x2*(256.0*x2 - 576.0) + 432.0) - 129.0) + 9.0); } break;
     305        for (int i = 0; i < vec->n; i++) { CHEB_EVAL_9(out->data.F64[i], vec->data.F64[i]); } break;
    295306      default:
    296307        psWarning ("Chebyshev orders higher than 9 are not yet coded\n");
    … …  
    332343}
    333344
    334 // XXX: You can do this without having to psAlloc() vector d.
    335 // XXX: How does the mask vector affect Clenshaw's formula?
    336 // NOTE: We assume that x is scaled between -1.0 and 1.0;
    337 // XXX: Create a faster version for low-order Chebyshevs.
    338 static psF64 chebPolynomial1DEval(
    339     psF64 x,
    340     const psPolynomial1D* poly)
    341 {
    342     PS_ASSERT_DOUBLE_WITHIN_RANGE(x, -1.0, 1.0, NAN);
     345static psF64 chebPolynomial1DEval(psF64 x, const psPolynomial1D* poly) {
     346
    343347    PS_ASSERT_INT_LARGER_THAN_OR_EQUAL(poly->nX, 0, NAN);
     348
     349    psF64 xNorm = x*poly->scale[0] + poly->zero[0];
     350
     351    psF64 polySum = 0.0;
     352
     353    for (int ix = 0; ix <= poly->nX; ix++) {
     354        if (poly->coeffMask[ix] & PS_POLY_MASK_SET) continue;
     355        psF64 xCheb = NAN;
     356        switch (ix) {
     357          case 0: CHEB_EVAL_0 (xCheb, xNorm); break;
     358          case 1: CHEB_EVAL_1 (xCheb, xNorm); break;
     359          case 2: CHEB_EVAL_2 (xCheb, xNorm); break;
     360          case 3: CHEB_EVAL_3 (xCheb, xNorm); break;
     361          case 4: CHEB_EVAL_4 (xCheb, xNorm); break;
     362          case 5: CHEB_EVAL_5 (xCheb, xNorm); break;
     363          case 6: CHEB_EVAL_6 (xCheb, xNorm); break;
     364          case 7: CHEB_EVAL_7 (xCheb, xNorm); break;
     365          case 8: CHEB_EVAL_8 (xCheb, xNorm); break;
     366          case 9: CHEB_EVAL_9 (xCheb, xNorm); break;
     367          default:
     368            break;
     369        }
     370        polySum += poly->coeff[ix] * xCheb;
     371    }
     372    return polySum;
     373}
     374
     375/*** version 1 is a general case and could be used for Norder > 9.  ***/
     376# ifdef CHEB_VERSION_1
     377void oldcode_1(void) {
    344378    psVector *d;
     379    psF64 tmp = 0.0;
    345380
    346381    unsigned int nTerms = 1 + poly->nX;
    347382    unsigned int i;
    348     psF64 tmp = 0.0;
    349383
    350384    // Special case where the Chebyshev poly is constant.
    … …  
    367401    }
    368402
    369     if (1) {
    370         // General case where the Chebyshev poly has 2 or more terms.
    371         d = psVectorAlloc(nTerms, PS_TYPE_F64);
    372         if (!(poly->coeffMask[nTerms-1] & PS_POLY_MASK_SET)) {
    373             d->data.F64[nTerms-1] = poly->coeff[nTerms-1];
    374         } else {
    375             d->data.F64[nTerms-1] = 0.0;
    376         }
    377 
    378         d->data.F64[nTerms-2] = (2.0 * x * d->data.F64[nTerms-1]);
    379         if (!(poly->coeffMask[nTerms-2] & PS_POLY_MASK_SET)) {
    380             d->data.F64[nTerms-2] += poly->coeff[nTerms-2];
    381         }
    382 
    383         for (i=nTerms-3;i>=1;i--) {
    384             d->data.F64[i] = (2.0 * x * d->data.F64[i+1]) - (d->data.F64[i+2]);
    385             if (!(poly->coeffMask[i] & PS_POLY_MASK_SET)) {
    386                 d->data.F64[i] += poly->coeff[i];
    387             }
    388         }
    389 
    390         tmp = (x * d->data.F64[1]) - (d->data.F64[2]);
    391         if (!(poly->coeffMask[0] & PS_POLY_MASK_SET)) {
    392             tmp += (0.5 * poly->coeff[0]);
    393         }
    394         psFree(d);
     403    // General case where the Chebyshev poly has 2 or more terms.
     404    d = psVectorAlloc(nTerms, PS_TYPE_F64);
     405    if (!(poly->coeffMask[nTerms-1] & PS_POLY_MASK_SET)) {
     406        d->data.F64[nTerms-1] = poly->coeff[nTerms-1];
    395407    } else {
    396         // XXX: This is old code that does not use Clenshaw's formula.  Get rid of it.
    397         psPolynomial1D **chebPolys = p_psCreateChebyshevPolys(1 + poly->nX);
    398 
    399         tmp = 0.0;
    400         for (psS32 i=0;i<(1 + poly->nX);i++) {
    401             tmp+= (poly->coeff[i] * psPolynomial1DEval(chebPolys[i], x));
    402         }
    403         tmp-= (poly->coeff[0]/2.0);
    404 
    405         for (psS32 i=0;i<(1 + poly->nX);i++) {
    406             psFree(chebPolys[i]);
    407         }
    408         psFree(chebPolys);
    409     }
     408        d->data.F64[nTerms-1] = 0.0;
     409    }
     410
     411    d->data.F64[nTerms-2] = (2.0 * x * d->data.F64[nTerms-1]);
     412    if (!(poly->coeffMask[nTerms-2] & PS_POLY_MASK_SET)) {
     413        d->data.F64[nTerms-2] += poly->coeff[nTerms-2];
     414    }
     415
     416    for (i=nTerms-3;i>=1;i--) {
     417        d->data.F64[i] = (2.0 * x * d->data.F64[i+1]) - (d->data.F64[i+2]);
     418        if (!(poly->coeffMask[i] & PS_POLY_MASK_SET)) {
     419            d->data.F64[i] += poly->coeff[i];
     420        }
     421    }
     422
     423    tmp = (x * d->data.F64[1]) - (d->data.F64[2]);
     424    if (!(poly->coeffMask[0] & PS_POLY_MASK_SET)) {
     425        tmp += (0.5 * poly->coeff[0]);
     426    }
     427    psFree(d);
     428}
     429# endif
     430
     431/*** version 0 should be removed when version 2 is ready ***/
     432# ifdef CHEB_VERSION_0
     433void oldcode_0(void) {
     434    // XXX: This is old code that does not use Clenshaw's formula.  Get rid of it.
     435    psPolynomial1D **chebPolys = p_psCreateChebyshevPolys(1 + poly->nX);
     436
     437    tmp = 0.0;
     438    for (psS32 i=0;i<(1 + poly->nX);i++) {
     439        tmp+= (poly->coeff[i] * psPolynomial1DEval(chebPolys[i], x));
     440    }
     441    tmp-= (poly->coeff[0]/2.0);
     442
     443    for (psS32 i=0;i<(1 + poly->nX);i++) {
     444        psFree(chebPolys[i]);
     445    }
     446    psFree(chebPolys);
    410447
    411448    return(tmp);
    412449}
     450# endif
    413451
    414452static psF64 ordPolynomial2DEval(psF64 x,
    … …  
    442480                                  const psPolynomial2D* poly)
    443481{
    444   // XXX transform x,y to chebyshev range
    445   // PS_ASSERT_DOUBLE_WITHIN_RANGE(x, -1.0, 1.0, 0.0);
    446   // PS_ASSERT_DOUBLE_WITHIN_RANGE(y, -1.0, 1.0, 0.0);
    447482    PS_ASSERT_POLY_NON_NULL(poly, NAN);
    448 
    449     unsigned int loop_x = 0;
    450     unsigned int loop_y = 0;
    451     unsigned int i = 0;
    452     psF64 polySum = 0.0;
    453     psPolynomial1D* *chebPolys = NULL;
    454     unsigned int maxChebyPoly = 0;
    455483
    456484    psF64 xNorm = x*poly->scale[0] + poly->zero[0];
    457485    psF64 yNorm = y*poly->scale[1] + poly->zero[1];
    458486
    459     // Determine how many Chebyshev polynomials
    460     // are needed, then create them.
    461     maxChebyPoly = poly->nX;
    462     if (poly->nY > maxChebyPoly) {
    463         maxChebyPoly = poly->nY;
    464     }
    465     chebPolys = p_psCreateChebyshevPolys(maxChebyPoly + 1);
    466 
    467     for (loop_x = 0; loop_x < (1 + poly->nX); loop_x++) {
    468         for (loop_y = 0; loop_y < (1 + poly->nY); loop_y++) {
    469             if (!(poly->coeffMask[loop_x][loop_y] & PS_POLY_MASK_SET)) {
    470                 polySum += poly->coeff[loop_x][loop_y] *
    471                            psPolynomial1DEval(chebPolys[loop_x], xNorm) *
    472                            psPolynomial1DEval(chebPolys[loop_y], yNorm);
    473             }
    474         }
    475     }
    476     for (i=0;i<maxChebyPoly+1;i++) {
    477         psFree(chebPolys[i]);
    478     }
    479     psFree(chebPolys);
     487    psF64 polySum = 0.0;
     488
     489    // XXX this could be quicker if we saved the N xvalues are re-used the resuls
     490    for (int ix = 0; ix <= poly->nX; ix++) {
     491        psF64 xCheb = NAN;
     492        switch (ix) {
     493          case 0: CHEB_EVAL_0 (xCheb, xNorm); break;
     494          case 1: CHEB_EVAL_1 (xCheb, xNorm); break;
     495          case 2: CHEB_EVAL_2 (xCheb, xNorm); break;
     496          case 3: CHEB_EVAL_3 (xCheb, xNorm); break;
     497          case 4: CHEB_EVAL_4 (xCheb, xNorm); break;
     498          case 5: CHEB_EVAL_5 (xCheb, xNorm); break;
     499          case 6: CHEB_EVAL_6 (xCheb, xNorm); break;
     500          case 7: CHEB_EVAL_7 (xCheb, xNorm); break;
     501          case 8: CHEB_EVAL_8 (xCheb, xNorm); break;
     502          case 9: CHEB_EVAL_9 (xCheb, xNorm); break;
     503          default:
     504            break;
     505        }
     506        for (int iy = 0; iy <= poly->nY; iy++) {
     507            if (poly->coeffMask[ix][iy] & PS_POLY_MASK_SET) continue;
     508            psF64 yCheb = NAN;
     509            switch (iy) {
     510              case 0: CHEB_EVAL_0 (yCheb, yNorm); break;
     511              case 1: CHEB_EVAL_1 (yCheb, yNorm); break;
     512              case 2: CHEB_EVAL_2 (yCheb, yNorm); break;
     513              case 3: CHEB_EVAL_3 (yCheb, yNorm); break;
     514              case 4: CHEB_EVAL_4 (yCheb, yNorm); break;
     515              case 5: CHEB_EVAL_5 (yCheb, yNorm); break;
     516              case 6: CHEB_EVAL_6 (yCheb, yNorm); break;
     517              case 7: CHEB_EVAL_7 (yCheb, yNorm); break;
     518              case 8: CHEB_EVAL_8 (yCheb, yNorm); break;
     519              case 9: CHEB_EVAL_9 (yCheb, yNorm); break;
     520              default:
     521                break;
     522            }
     523            polySum += poly->coeff[ix][iy] * xCheb * yCheb;
     524        }
     525    }
    480526    return(polySum);
    481527}
    … …  
    9561002}
    9571003
     1004/* note these functions accept unscaled values and apply the scaling saved on poly */
    9581005psF64 psPolynomial1DEval(const psPolynomial1D* poly,
    9591006                         psF64 x)
    … …  
    9631010    if (poly->type == PS_POLYNOMIAL_ORD) {
    9641011        return(ordPolynomial1DEval(x, poly));
    965     } else if (poly->type == PS_POLYNOMIAL_CHEB) {
     1012    }
     1013    if (poly->type == PS_POLYNOMIAL_CHEB) {
    9661014        return(chebPolynomial1DEval(x, poly));
    967     } else {
    968         psError(PS_ERR_BAD_PARAMETER_TYPE, true,
    969                 _("Unknown polynomial type 0x%x found.  Evaluation failed."),
    970                 poly->type);
    971     }
     1015    }
     1016    psError(PS_ERR_BAD_PARAMETER_TYPE, true,
     1017            _("Unknown polynomial type 0x%x found.  Evaluation failed."),
     1018            poly->type);
     1019
    9721020    return(NAN);
    9731021}
    … …  
    10131061    if (poly->type == PS_POLYNOMIAL_ORD) {
    10141062        return(ordPolynomial2DEval(x, y, poly));
    1015     } else if (poly->type == PS_POLYNOMIAL_CHEB) {
     1063    }
     1064    if (poly->type == PS_POLYNOMIAL_CHEB) {
    10161065        return(chebPolynomial2DEval(x, y, poly));
    1017     } else {
    1018         psError(PS_ERR_BAD_PARAMETER_TYPE, true,
    1019                 _("Unknown polynomial type 0x%x found.  Evaluation failed."),
    1020                 poly->type);
    1021     }
     1066    }
     1067    psError(PS_ERR_BAD_PARAMETER_TYPE, true,
     1068            _("Unknown polynomial type 0x%x found.  Evaluation failed."),
     1069            poly->type);
    10221070    return(NAN);
    10231071}
    … …  
    10771125            for (int jy = 0; jy < nYterm; jy++) {
    10781126                psVector *jyCheb = yPolySet->data[jy];
     1127                if (poly->coeffMask[jx][jy] & PS_POLY_MASK_SET) continue;
    10791128                sum += poly->coeff[jx][jy] * jxCheb->data.F64[i] * jyCheb->data.F64[i];
    10801129            }
Note: See TracChangeset for help on using the changeset viewer.