Index: /branches/eam_branches/ipp-dev-20210817/psLib/src/math/psPolynomial.c =================================================================== --- /branches/eam_branches/ipp-dev-20210817/psLib/src/math/psPolynomial.c (revision 41838) +++ /branches/eam_branches/ipp-dev-20210817/psLib/src/math/psPolynomial.c (revision 41839) @@ -257,4 +257,15 @@ } +# define CHEB_EVAL_0(OUT,IN) {OUT = 1.0;} +# define CHEB_EVAL_1(OUT,IN) { OUT = IN; } +# define CHEB_EVAL_2(OUT,IN) {psF64 X2 = PS_SQR(IN); OUT = 2.0*X2 - 1.0; } +# define CHEB_EVAL_3(OUT,IN) {psF64 X2 = PS_SQR(IN); OUT = IN*(4.0*X2 - 3.0); } +# define CHEB_EVAL_4(OUT,IN) {psF64 X2 = PS_SQR(IN); OUT = X2*(8.0*X2 - 8.0) + 1.0; } +# define CHEB_EVAL_5(OUT,IN) {psF64 X2 = PS_SQR(IN); OUT = IN *(X2*(16.0*X2 - 20.0) + 5.0); } +# define CHEB_EVAL_6(OUT,IN) {psF64 X2 = PS_SQR(IN); OUT = X2*(X2*(32.0*X2 - 48.0) + 18.0) - 1.0; } +# define CHEB_EVAL_7(OUT,IN) {psF64 X2 = PS_SQR(IN); OUT = IN *(X2*(X2*(64.0*X2 - 112.0) + 56.0) - 7.0); } +# 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; } +# 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); } + /** This function generates a vector containing the values of a Chebyshev polynomial of the given order evaluated at the coordinates given by the input vector, i.e., this @@ -274,23 +285,23 @@ switch (order) { case 0: - for (int i = 0; i < vec->n; i++) { out->data.F64[i] = 1.0; } break; + for (int i = 0; i < vec->n; i++) { CHEB_EVAL_0(out->data.F64[i], vec->data.F64[i]); } break; case 1: - for (int i = 0; i < vec->n; i++) { psF64 x = vec->data.F64[i]; out->data.F64[i] = x; } break; + for (int i = 0; i < vec->n; i++) { CHEB_EVAL_1(out->data.F64[i], vec->data.F64[i]); } break; case 2: - 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; + for (int i = 0; i < vec->n; i++) { CHEB_EVAL_2(out->data.F64[i], vec->data.F64[i]); } break; case 3: - 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; + for (int i = 0; i < vec->n; i++) { CHEB_EVAL_3(out->data.F64[i], vec->data.F64[i]); } break; case 4: - 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; + for (int i = 0; i < vec->n; i++) { CHEB_EVAL_4(out->data.F64[i], vec->data.F64[i]); } break; case 5: - 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; + for (int i = 0; i < vec->n; i++) { CHEB_EVAL_5(out->data.F64[i], vec->data.F64[i]); } break; case 6: - 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; + for (int i = 0; i < vec->n; i++) { CHEB_EVAL_6(out->data.F64[i], vec->data.F64[i]); } break; case 7: - 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; + for (int i = 0; i < vec->n; i++) { CHEB_EVAL_7(out->data.F64[i], vec->data.F64[i]); } break; case 8: - 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; + for (int i = 0; i < vec->n; i++) { CHEB_EVAL_8(out->data.F64[i], vec->data.F64[i]); } break; case 9: - 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; + for (int i = 0; i < vec->n; i++) { CHEB_EVAL_9(out->data.F64[i], vec->data.F64[i]); } break; default: psWarning ("Chebyshev orders higher than 9 are not yet coded\n"); @@ -332,19 +343,42 @@ } -// XXX: You can do this without having to psAlloc() vector d. -// XXX: How does the mask vector affect Clenshaw's formula? -// NOTE: We assume that x is scaled between -1.0 and 1.0; -// XXX: Create a faster version for low-order Chebyshevs. -static psF64 chebPolynomial1DEval( - psF64 x, - const psPolynomial1D* poly) -{ - PS_ASSERT_DOUBLE_WITHIN_RANGE(x, -1.0, 1.0, NAN); +static psF64 chebPolynomial1DEval(psF64 x, const psPolynomial1D* poly) { + PS_ASSERT_INT_LARGER_THAN_OR_EQUAL(poly->nX, 0, NAN); + + psF64 xNorm = x*poly->scale[0] + poly->zero[0]; + + psF64 polySum = 0.0; + + for (int ix = 0; ix <= poly->nX; ix++) { + if (poly->coeffMask[ix] & PS_POLY_MASK_SET) continue; + psF64 xCheb = NAN; + switch (ix) { + case 0: CHEB_EVAL_0 (xCheb, xNorm); break; + case 1: CHEB_EVAL_1 (xCheb, xNorm); break; + case 2: CHEB_EVAL_2 (xCheb, xNorm); break; + case 3: CHEB_EVAL_3 (xCheb, xNorm); break; + case 4: CHEB_EVAL_4 (xCheb, xNorm); break; + case 5: CHEB_EVAL_5 (xCheb, xNorm); break; + case 6: CHEB_EVAL_6 (xCheb, xNorm); break; + case 7: CHEB_EVAL_7 (xCheb, xNorm); break; + case 8: CHEB_EVAL_8 (xCheb, xNorm); break; + case 9: CHEB_EVAL_9 (xCheb, xNorm); break; + default: + break; + } + polySum += poly->coeff[ix] * xCheb; + } + return polySum; +} + +/*** version 1 is a general case and could be used for Norder > 9. ***/ +# ifdef CHEB_VERSION_1 +void oldcode_1(void) { psVector *d; + psF64 tmp = 0.0; unsigned int nTerms = 1 + poly->nX; unsigned int i; - psF64 tmp = 0.0; // Special case where the Chebyshev poly is constant. @@ -367,48 +401,52 @@ } - if (1) { - // General case where the Chebyshev poly has 2 or more terms. - d = psVectorAlloc(nTerms, PS_TYPE_F64); - if (!(poly->coeffMask[nTerms-1] & PS_POLY_MASK_SET)) { - d->data.F64[nTerms-1] = poly->coeff[nTerms-1]; - } else { - d->data.F64[nTerms-1] = 0.0; - } - - d->data.F64[nTerms-2] = (2.0 * x * d->data.F64[nTerms-1]); - if (!(poly->coeffMask[nTerms-2] & PS_POLY_MASK_SET)) { - d->data.F64[nTerms-2] += poly->coeff[nTerms-2]; - } - - for (i=nTerms-3;i>=1;i--) { - d->data.F64[i] = (2.0 * x * d->data.F64[i+1]) - (d->data.F64[i+2]); - if (!(poly->coeffMask[i] & PS_POLY_MASK_SET)) { - d->data.F64[i] += poly->coeff[i]; - } - } - - tmp = (x * d->data.F64[1]) - (d->data.F64[2]); - if (!(poly->coeffMask[0] & PS_POLY_MASK_SET)) { - tmp += (0.5 * poly->coeff[0]); - } - psFree(d); + // General case where the Chebyshev poly has 2 or more terms. + d = psVectorAlloc(nTerms, PS_TYPE_F64); + if (!(poly->coeffMask[nTerms-1] & PS_POLY_MASK_SET)) { + d->data.F64[nTerms-1] = poly->coeff[nTerms-1]; } else { - // XXX: This is old code that does not use Clenshaw's formula. Get rid of it. - psPolynomial1D **chebPolys = p_psCreateChebyshevPolys(1 + poly->nX); - - tmp = 0.0; - for (psS32 i=0;i<(1 + poly->nX);i++) { - tmp+= (poly->coeff[i] * psPolynomial1DEval(chebPolys[i], x)); - } - tmp-= (poly->coeff[0]/2.0); - - for (psS32 i=0;i<(1 + poly->nX);i++) { - psFree(chebPolys[i]); - } - psFree(chebPolys); - } + d->data.F64[nTerms-1] = 0.0; + } + + d->data.F64[nTerms-2] = (2.0 * x * d->data.F64[nTerms-1]); + if (!(poly->coeffMask[nTerms-2] & PS_POLY_MASK_SET)) { + d->data.F64[nTerms-2] += poly->coeff[nTerms-2]; + } + + for (i=nTerms-3;i>=1;i--) { + d->data.F64[i] = (2.0 * x * d->data.F64[i+1]) - (d->data.F64[i+2]); + if (!(poly->coeffMask[i] & PS_POLY_MASK_SET)) { + d->data.F64[i] += poly->coeff[i]; + } + } + + tmp = (x * d->data.F64[1]) - (d->data.F64[2]); + if (!(poly->coeffMask[0] & PS_POLY_MASK_SET)) { + tmp += (0.5 * poly->coeff[0]); + } + psFree(d); +} +# endif + +/*** version 0 should be removed when version 2 is ready ***/ +# ifdef CHEB_VERSION_0 +void oldcode_0(void) { + // XXX: This is old code that does not use Clenshaw's formula. Get rid of it. + psPolynomial1D **chebPolys = p_psCreateChebyshevPolys(1 + poly->nX); + + tmp = 0.0; + for (psS32 i=0;i<(1 + poly->nX);i++) { + tmp+= (poly->coeff[i] * psPolynomial1DEval(chebPolys[i], x)); + } + tmp-= (poly->coeff[0]/2.0); + + for (psS32 i=0;i<(1 + poly->nX);i++) { + psFree(chebPolys[i]); + } + psFree(chebPolys); return(tmp); } +# endif static psF64 ordPolynomial2DEval(psF64 x, @@ -442,40 +480,48 @@ const psPolynomial2D* poly) { - // XXX transform x,y to chebyshev range - // PS_ASSERT_DOUBLE_WITHIN_RANGE(x, -1.0, 1.0, 0.0); - // PS_ASSERT_DOUBLE_WITHIN_RANGE(y, -1.0, 1.0, 0.0); PS_ASSERT_POLY_NON_NULL(poly, NAN); - - unsigned int loop_x = 0; - unsigned int loop_y = 0; - unsigned int i = 0; - psF64 polySum = 0.0; - psPolynomial1D* *chebPolys = NULL; - unsigned int maxChebyPoly = 0; psF64 xNorm = x*poly->scale[0] + poly->zero[0]; psF64 yNorm = y*poly->scale[1] + poly->zero[1]; - // Determine how many Chebyshev polynomials - // are needed, then create them. - maxChebyPoly = poly->nX; - if (poly->nY > maxChebyPoly) { - maxChebyPoly = poly->nY; - } - chebPolys = p_psCreateChebyshevPolys(maxChebyPoly + 1); - - for (loop_x = 0; loop_x < (1 + poly->nX); loop_x++) { - for (loop_y = 0; loop_y < (1 + poly->nY); loop_y++) { - if (!(poly->coeffMask[loop_x][loop_y] & PS_POLY_MASK_SET)) { - polySum += poly->coeff[loop_x][loop_y] * - psPolynomial1DEval(chebPolys[loop_x], xNorm) * - psPolynomial1DEval(chebPolys[loop_y], yNorm); - } - } - } - for (i=0;inX; ix++) { + psF64 xCheb = NAN; + switch (ix) { + case 0: CHEB_EVAL_0 (xCheb, xNorm); break; + case 1: CHEB_EVAL_1 (xCheb, xNorm); break; + case 2: CHEB_EVAL_2 (xCheb, xNorm); break; + case 3: CHEB_EVAL_3 (xCheb, xNorm); break; + case 4: CHEB_EVAL_4 (xCheb, xNorm); break; + case 5: CHEB_EVAL_5 (xCheb, xNorm); break; + case 6: CHEB_EVAL_6 (xCheb, xNorm); break; + case 7: CHEB_EVAL_7 (xCheb, xNorm); break; + case 8: CHEB_EVAL_8 (xCheb, xNorm); break; + case 9: CHEB_EVAL_9 (xCheb, xNorm); break; + default: + break; + } + for (int iy = 0; iy <= poly->nY; iy++) { + if (poly->coeffMask[ix][iy] & PS_POLY_MASK_SET) continue; + psF64 yCheb = NAN; + switch (iy) { + case 0: CHEB_EVAL_0 (yCheb, yNorm); break; + case 1: CHEB_EVAL_1 (yCheb, yNorm); break; + case 2: CHEB_EVAL_2 (yCheb, yNorm); break; + case 3: CHEB_EVAL_3 (yCheb, yNorm); break; + case 4: CHEB_EVAL_4 (yCheb, yNorm); break; + case 5: CHEB_EVAL_5 (yCheb, yNorm); break; + case 6: CHEB_EVAL_6 (yCheb, yNorm); break; + case 7: CHEB_EVAL_7 (yCheb, yNorm); break; + case 8: CHEB_EVAL_8 (yCheb, yNorm); break; + case 9: CHEB_EVAL_9 (yCheb, yNorm); break; + default: + break; + } + polySum += poly->coeff[ix][iy] * xCheb * yCheb; + } + } return(polySum); } @@ -956,4 +1002,5 @@ } +/* note these functions accept unscaled values and apply the scaling saved on poly */ psF64 psPolynomial1DEval(const psPolynomial1D* poly, psF64 x) @@ -963,11 +1010,12 @@ if (poly->type == PS_POLYNOMIAL_ORD) { return(ordPolynomial1DEval(x, poly)); - } else if (poly->type == PS_POLYNOMIAL_CHEB) { + } + if (poly->type == PS_POLYNOMIAL_CHEB) { return(chebPolynomial1DEval(x, poly)); - } else { - psError(PS_ERR_BAD_PARAMETER_TYPE, true, - _("Unknown polynomial type 0x%x found. Evaluation failed."), - poly->type); - } + } + psError(PS_ERR_BAD_PARAMETER_TYPE, true, + _("Unknown polynomial type 0x%x found. Evaluation failed."), + poly->type); + return(NAN); } @@ -1013,11 +1061,11 @@ if (poly->type == PS_POLYNOMIAL_ORD) { return(ordPolynomial2DEval(x, y, poly)); - } else if (poly->type == PS_POLYNOMIAL_CHEB) { + } + if (poly->type == PS_POLYNOMIAL_CHEB) { return(chebPolynomial2DEval(x, y, poly)); - } else { - psError(PS_ERR_BAD_PARAMETER_TYPE, true, - _("Unknown polynomial type 0x%x found. Evaluation failed."), - poly->type); - } + } + psError(PS_ERR_BAD_PARAMETER_TYPE, true, + _("Unknown polynomial type 0x%x found. Evaluation failed."), + poly->type); return(NAN); } @@ -1077,4 +1125,5 @@ for (int jy = 0; jy < nYterm; jy++) { psVector *jyCheb = yPolySet->data[jy]; + if (poly->coeffMask[jx][jy] & PS_POLY_MASK_SET) continue; sum += poly->coeff[jx][jy] * jxCheb->data.F64[i] * jyCheb->data.F64[i]; }