Changeset 41831 for branches/eam_branches/ipp-dev-20210817

Timestamp:

Oct 11, 2021, 11:33:33 AM (5 years ago)

Author:

eugene

Message:

add chebyshev polynomials to 2D fitting; add chebyshev support functions; use chebyshev polynomials for inverse transformation

Location:

branches/eam_branches/ipp-dev-20210817/psLib

Files:

• src/astro/psCoord.c (modified) (9 diffs)
• src/astro/psCoord.h (modified) (1 diff)
• src/db/psDB.c (modified) (1 diff)
• src/math/psMinimizePolyFit.c (modified) (4 diffs)
• src/math/psPolynomial.c (modified) (13 diffs)
• src/math/psPolynomial.h (modified) (5 diffs)
• test/math/tap_psPolyFit2DCheb.c (added)

branches/eam_branches/ipp-dev-20210817/psLib/src/astro/psCoord.c

@@ -40 +40 @@
 #include "psMinimizePolyFit.h"
 
+# define TEST_SAVE_INVERSE_TRANSFORM 0
+
+# if (TEST_SAVE_INVERSE_TRANSFORM)
+# include "psBinaryOp.h"
+# endif
+
 # define ELIXIR_CODE 1
 
@@ -86 +92 @@
     }
 
-    psPlaneTransform *out = psPlaneTransformAlloc(1, 1);
+    // since the output polynomial is 1st order, a Chebyshev is not really useful
+    psPlaneTransform *out = psPlaneTransformAlloc(1, 1, PS_POLYNOMIAL_ORD); 
 
     psF64 r12 = 0.0;
@@ -191 +198 @@
 }
 
-psPlaneTransform* psPlaneTransformAlloc(int order1, int order2)
+psPlaneTransform* psPlaneTransformAlloc(int order1, int order2, psPolynomialType type)
 {
     PS_ASSERT_INT_NONNEGATIVE(order1, NULL);
@@ -198 +205 @@
     psPlaneTransform *pt = psAlloc(sizeof(psPlaneTransform));
 
-    pt->x = psPolynomial2DAlloc(PS_POLYNOMIAL_ORD, order1, order2);
-    pt->y = psPolynomial2DAlloc(PS_POLYNOMIAL_ORD, order1, order2);
+    pt->x = psPolynomial2DAlloc(type, order1, order2);
+    pt->y = psPolynomial2DAlloc(type, order1, order2);
 
     psMemSetDeallocator(pt, (psFreeFunc) planeTransformFree);
@@ -797 +804 @@
     }
 
+    // both polynomials in both input transforms must match type -- and for now be ORD
+    psAssert (trans1->x->type == PS_POLYNOMIAL_ORD, "fix for CHEB");
+    psAssert (trans1->y->type == PS_POLYNOMIAL_ORD, "fix for CHEB");
+    psAssert (trans2->x->type == PS_POLYNOMIAL_ORD, "fix for CHEB");
+    psAssert (trans2->y->type == PS_POLYNOMIAL_ORD, "fix for CHEB");
+
     //
     // Determine the size of the new psPlaneTransform.
@@ -814 +827 @@
     psPlaneTransform *myPT = NULL;
     if (out == NULL) {
-        myPT = psPlaneTransformAlloc(orderX, orderY);
+        myPT = psPlaneTransformAlloc(orderX, orderY, PS_POLYNOMIAL_ORD);
     } else {
         if ((out->x->nX == orderX) &&
@@ -834 +847 @@
         } else {
             psFree(out);
-            myPT = psPlaneTransformAlloc(orderX, orderY);
+            myPT = psPlaneTransformAlloc(orderX, orderY, PS_POLYNOMIAL_ORD);
         }
     }
@@ -1037 +1050 @@
     psPlaneTransform *myPT = NULL;
     if (out == NULL) {
-        myPT = psPlaneTransformAlloc(order, order);
+      myPT = psPlaneTransformAlloc(order, order, PS_POLYNOMIAL_CHEB);
     } else {
       // the user has supplied a model with a specific order : fit that order
@@ -1071 +1084 @@
     result &= psVectorFitPolynomial2D(myPT->x, NULL, 0, xOut, NULL, xIn, yIn);
     result &= psVectorFitPolynomial2D(myPT->y, NULL, 0, yOut, NULL, xIn, yIn);
+
+# if (TEST_SAVE_INVERSE_TRANSFORM)
+
+    psVector *xFit = psPolynomial2DEvalVector (myPT->x, xIn, yIn);
+    psVector *xRes = (psVector *) psBinaryOp (NULL, xOut, "-", xFit);
+
+    psVector *yFit = psPolynomial2DEvalVector (myPT->y, xIn, yIn);
+    psVector *yRes = (psVector *) psBinaryOp (NULL, yOut, "-", yFit);
+
+    static int nOut = 0;
+    char filename[1024];
+    snprintf (filename, 1024, "test.fit.%03d.ply", nOut); 
+    FILE *fp = fopen (filename, "w");
+    
+    fprintf (fp, " ---- xFit ---- \n");
+
+    for (int ix = 0; ix < myPT->x->nX + 1; ix++) {
+      for (int iy = 0; iy < myPT->x->nY + 1; iy++) {
+        fprintf (fp, "%18.12e ", myPT->x->coeff[ix][iy]);
+      }
+      fprintf (fp, "\n");
+    }
+    fprintf (fp, " ---- yFit ---- \n");
+
+    for (int ix = 0; ix < myPT->y->nX + 1; ix++) {
+      for (int iy = 0; iy < myPT->y->nY + 1; iy++) {
+        fprintf (fp, "%18.12e ", myPT->y->coeff[ix][iy]);
+      }
+      fprintf (fp, "\n");
+    }
+    fclose (fp);
+    
+    snprintf (filename, 1024, "test.fit.%03d.dat", nOut); nOut ++;
+    FILE *f1 = fopen (filename, "w");
+    for (int i = 0; i < xFit->n; i++) {
+      fprintf (f1, "%d : %f %f : %f %f : %f %f : %f %f\n", i,
+               xIn->data.F64[i], yIn->data.F64[i],
+               xOut->data.F64[i], yOut->data.F64[i],
+               xFit->data.F64[i], yFit->data.F64[i],
+               xRes->data.F64[i], yRes->data.F64[i]);
+    }
+    fclose (f1);
+
+    psStats *myStats = psStatsAlloc (PS_STAT_SAMPLE_STDEV);
+    psVectorStats (myStats, xRes, NULL, NULL, 0); float dX = myStats->sampleStdev;
+    psVectorStats (myStats, yRes, NULL, NULL, 0); float dY = myStats->sampleStdev;
+    fprintf (stderr, "xRes Sigma: %f  --  yRes Sigma %f\n", dX, dY);
+
+    psFree (myStats);
+    psFree (xFit);
+    psFree (yFit);
+    psFree (xRes);
+    psFree (yRes);
+
+# endif
 
     psFree(inCoord);

branches/eam_branches/ipp-dev-20210817/psLib/src/astro/psCoord.h

@@ -178 +178 @@
 
 psPlaneTransform* psPlaneTransformAlloc(
-    int order1,                        ///< The order of the x term in the transform.
-    int order2                         ///< The order of the y term in the transform.
+    int order1,                       ///< The order of the x term in the transform.
+    int order2,                       ///< The order of the y term in the transform.
+    psPolynomialType type             ///< The polynomial type (ORD or CHEB) for this transform
 ) PS_ATTR_MALLOC;
 

branches/eam_branches/ipp-dev-20210817/psLib/src/db/psDB.c

@@ -2875 +2875 @@
 
     switch (pType) {
-    case PS_DATA_S8:
+      case PS_DATA_S8:
         isNaN = PS_IS_NAN(psS8, data, PS_MAX_S8);
         break;
-    case PS_DATA_S16:
+      case PS_DATA_S16:
         isNaN = PS_IS_NAN(psS16, data, PS_MAX_S16);
         break;
-    case PS_DATA_S32:
+      case PS_DATA_S32:
         isNaN = PS_IS_NAN(psS32, data, PS_MAX_S32);
         break;
-    case PS_DATA_S64:
+      case PS_DATA_S64:
         isNaN = PS_IS_NAN(psS64, data, PS_MAX_S64);
         break;
-    case PS_DATA_U8:
+      case PS_DATA_U8:
         isNaN = PS_IS_NAN(psU8, data, PS_MAX_U8);
         break;
-    case PS_DATA_U16:
+      case PS_DATA_U16:
         isNaN = PS_IS_NAN(psU16, data, PS_MAX_U16);
         break;
-    case PS_DATA_U32:
+      case PS_DATA_U32:
         isNaN = PS_IS_NAN(psU32, data, PS_MAX_U32);
         break;
-    case PS_DATA_U64:
+      case PS_DATA_U64:
         isNaN = PS_IS_NAN(psU64, data, PS_MAX_U64);
         break;
-    case PS_DATA_F32:
-      isNaN = isnan(*((psF32 *) data));
-      break;
-    case PS_DATA_F64:
-      isNaN = isnan(*((psF64 *) data));
-      break;
-    case PS_DATA_BOOL:
-        // XXX: what is NaN for a bool?
-        isNaN = PS_IS_NAN(psU8, data, PS_MAX_U8);
+      case PS_DATA_F32:
+        isNaN = !isfinite(*((psF32 *) data)); // trap nan, +inf, -inf
+        break;
+      case PS_DATA_F64:
+        isNaN = !isfinite(*((psF64 *) data)); // trap nan, +inf, -inf
+        break;
+      case PS_DATA_BOOL:
+        isNaN = PS_IS_NAN(psU8, data, PS_MAX_U8); // probably meaningless
         break;
     }

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

@@ -35 +35 @@
 
 #include "psMinimizePolyFit.h"
+#include "psAbort.h"
 #include "psAssert.h"
 #include "psMinimizeLMM.h"  // For Gauss-Jordan routines
@@ -1225 +1226 @@
 
 /******************************************************************************
+VectorFitPolynomial2DCheb(myPoly, *mask, maskValue, *f, *fErr, *x, *y): This is
+a private routine which will fit a 2-D polynomial to a set of (x, y)-(f)
+pairs.  All non-NULL vectors must be of type PS_TYPE_F64.
+ 
+ *****************************************************************************/
+static bool VectorFitPolynomial2DCheb(
+    psPolynomial2D* myPoly,
+    const psVector *f,
+    const psVector *x,
+    const psVector *y)
+{
+    psTrace("psLib.math", 4, "---- %s() begin ----\n", __func__);
+    PS_ASSERT_POLY_NON_NULL(myPoly, false);
+    PS_ASSERT_INT_NONNEGATIVE(myPoly->nX, false);
+    PS_ASSERT_INT_NONNEGATIVE(myPoly->nY, false);
+    PS_ASSERT_VECTOR_NON_NULL(f, false);
+    PS_ASSERT_VECTOR_TYPE(f, PS_TYPE_F64, false);
+    PS_ASSERT_VECTOR_NON_NULL(x, false);
+    PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_F64, false);
+    PS_ASSERT_VECTORS_SIZE_EQUAL(f, x, false);
+    PS_ASSERT_VECTOR_NON_NULL(y, false);
+    PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F64, false);
+    PS_ASSERT_VECTORS_SIZE_EQUAL(f, y, false);
+
+    // Number of polynomial terms
+    int nXterm = 1 + myPoly->nX;      // Number of terms in x
+    int nYterm = 1 + myPoly->nY;      // Number of terms in y
+    int nTerm = nXterm * nYterm;      // Total number of terms
+    if (nXterm > 9) {
+        psError(PS_ERR_UNKNOWN, false, "failed 2D chebyshev fit: orders higher than 9 are not yet coded\n");
+        return false;
+    }
+    if (nYterm > 9) {
+        psError(PS_ERR_UNKNOWN, false, "failed 2D chebyshev fit: orders higher than 9 are not yet coded\n");
+        return false;
+    }
+
+    // determine scale factors
+    if (!psChebyshevSetScale (myPoly, x, 0)) { psError(PS_ERR_UNKNOWN, false, "failed 2D chebyshev fit.\n"); return false; }
+    if (!psChebyshevSetScale (myPoly, y, 1)) { psError(PS_ERR_UNKNOWN, false, "failed 2D chebyshev fit.\n"); return false; }
+
+    // generate normalized vectors
+    psVector *xNorm = psChebyshevNormVector (myPoly, x, 0);
+    psVector *yNorm = psChebyshevNormVector (myPoly, y, 1);
+
+    // generate the N cheb polynomials based on xNorm, yNorm
+    psArray *xPolySet = psArrayAlloc (nXterm);
+    for (int i = 0; i < nXterm; i++) {
+        xPolySet->data[i] = psChebyshevPolyVector (xNorm, i);
+    }
+    psArray *yPolySet = psArrayAlloc (nYterm);
+    for (int i = 0; i < nYterm; i++) {
+        yPolySet->data[i] = psChebyshevPolyVector (yNorm, i);
+    }
+
+    psF64 *fData = f->data.F64;         // Dereference f
+
+    psImage *A = psImageAlloc(nTerm, nTerm, PS_TYPE_F64); // Least-squares matrix
+    psVector *B = psVectorAlloc(nTerm, PS_TYPE_F64); // Least-squares vector
+
+    // Initialize data structures (should not be able to fail)
+    psAssert (psImageInit(A, 0.0), "Could initialize data structures A");
+    psAssert (psVectorInit(B, 0.0), "Could initialize data structures B");
+
+    // Dereference stuff, to make the loop go faster
+    psF64 **matrix = A->data.F64;       // Dereference the least-squares matrix
+    psF64 *vector = B->data.F64;        // Dereference the least-squares vector
+
+    // loop over all elements of the data vector
+    for (int k = 0; k < x->n; k++) {
+
+        if (!finite(fData[k])) continue;
+    
+        // XXX can we only calculate the upper diagonal?
+        int nelem = 0;
+        for (int jx = 0; jx < nXterm; jx++) {
+            psVector *jxCheb = xPolySet->data[jx];
+            for (int jy = 0; jy < nYterm; jy++) {
+                psVector *jyCheb = yPolySet->data[jy];
+                psF64 chebValue = jxCheb->data.F64[k] * jyCheb->data.F64[k];
+                
+                vector[nelem] += fData[k] * chebValue;
+
+                int melem = 0;
+                for (int kx = 0; kx < nXterm; kx++) {
+                    psVector *kxCheb = xPolySet->data[kx];
+                    for (int ky = 0; ky < nYterm; ky++) {
+                        psVector *kyCheb = yPolySet->data[ky];
+                        matrix[nelem][melem] += chebValue * kxCheb->data.F64[k]*kyCheb->data.F64[k];
+                        melem++;
+                    }
+                }
+                nelem++;
+            }
+        }
+    }
+
+    if (psTraceGetLevel("psLib.math") >= 4) {
+        printf("Least-squares vector:\n");
+        for (int i = 0; i < nTerm; i++) {
+            printf("%f ", B->data.F64[i]);
+        }
+        printf("\n");
+        printf("Least-squares matrix:\n");
+        for (int i = 0; i < nTerm; i++) {
+            for (int j = 0; j < nTerm; j++) {
+                printf("%f ", A->data.F64[i][j]);
+            }
+            printf("\n");
+        }
+    }
+
+    bool status = false;
+    if (USE_GAUSS_JORDAN) {
+        status = psMatrixGJSolve(A, B);
+    } else {
+        status = psMatrixLUSolve(A, B);
+    }
+    if (!status) {
+        psError(PS_ERR_UNKNOWN, false, "Could not solve linear equations.\n");
+        goto escape;
+    } 
+
+    // unroll the result:
+    int nelem = 0;
+    for (int jx = 0; jx < nXterm; jx++) {
+        for (int jy = 0; jy < nYterm; jy++) {
+            myPoly->coeff[jx][jy]    = B->data.F64[nelem];
+            myPoly->coeffErr[jx][jy] = sqrt(A->data.F64[nelem][nelem]);
+            nelem ++;
+        }
+    }
+    psFree(A);
+    psFree(B);
+
+    psFree (xNorm);
+    psFree (yNorm);
+    psFree (xPolySet);
+    psFree (yPolySet);
+
+    return true;
+
+escape:
+    psFree (A);
+    psFree (B);
+    return false;
+}
+
+/******************************************************************************
 psVectorFitPolynomial2D():  This routine fits a 2D polynomial of arbitrary
 degree (specified in poly) to the data points (x, y)-(f) and returns that
@@ -1240 +1390 @@
 {
     PS_ASSERT_POLY_NON_NULL(poly, false);
-    PS_ASSERT_POLY_TYPE(poly, PS_POLYNOMIAL_ORD, false);
+    // PS_ASSERT_POLY_TYPE(poly, PS_POLYNOMIAL_ORD, false);
 
     PS_ASSERT_VECTOR_NON_NULL(f, false);
@@ -1277 +1427 @@
         break;
     case PS_POLYNOMIAL_CHEB:
-        if (mask != NULL) {
-            psLogMsg(__func__, PS_LOG_WARN, "WARNING: ignoring mask and maskValue with Chebyshev polynomials.\n");
-        }
-        psError(PS_ERR_UNKNOWN, true, "2-D Chebyshev polynomial vector fitting has not been implemented.  Returning NULL.\n");
-        result = false;
-        break;
+      if (mask != NULL) {
+          psLogMsg(__func__, PS_LOG_WARN, "WARNING: ignoring mask and maskValue with Chebyshev polynomials.\n");
+      }
+      if (fErr != NULL) {
+          psLogMsg(__func__, PS_LOG_WARN, "WARNING: ignoring error values for Chebyshev polynomials.\n");
+      }
+      result = VectorFitPolynomial2DCheb(poly, f64, x64, y64);
+      if (!result) {
+          psError(PS_ERR_UNKNOWN, true, "Could not fit polynomial.  Returning NULL.\n");
+      }
+      break;
     default:
         psError(PS_ERR_UNKNOWN, true, "Incorrect polynomial type.  Returning NULL.\n");

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

@@ -36 +36 @@
 #include "psLogMsg.h"
 #include "psPolynomial.h"
+#include "psAbort.h"
 #include "psAssert.h"
 
@@ -203 +204 @@
 
 
+/** This function calculates the appropriate scaling factors needed to normalize the
+ * input vector to the range -1 : +1.  These are stored on the polynomial in the given
+ * direction.
+ */
+bool psChebyshevSetScale (psPolynomial2D* myPoly, const psVector *vec, int dir) {
+
+    psAssert ((dir == 0) || (dir == 1), "invalid direction %d\n", dir);
+
+    // find the min and max of the vector
+    psF64 minValue = NAN;
+    psF64 maxValue = NAN;
+
+    for (int i = 0; i < vec->n; i++) {
+        if (isnan(vec->data.F64[i])) continue;
+        if (isnan(minValue)) { minValue = vec->data.F64[i]; }
+        if (isnan(maxValue)) { maxValue = vec->data.F64[i]; }
+        minValue = PS_MIN(minValue, vec->data.F64[i]);
+        maxValue = PS_MAX(maxValue, vec->data.F64[i]);
+    }
+    if (minValue == maxValue) {
+        psWarning ("insufficient data range to determine scale factors\n");
+        return false;
+    }
+
+    myPoly->scale[dir] = 2.0 / (maxValue - minValue);
+    myPoly->zero[dir]  = 1 - myPoly->scale[dir] * maxValue;
+    return true;
+}
+
+/** This function generates a normalized vector in the range -1 : +1 based on the input
+    vector using the scale factors stored in myPoly in the given direction.
+ */
+psVector *psChebyshevNormVector (const psPolynomial2D* myPoly, const psVector *vec, int dir) {
+
+    psVector *norm = psVectorAlloc (vec->n, PS_TYPE_F64);
+
+    if (vec->type.type == PS_TYPE_F64) {
+        for (int i = 0; i < vec->n; i++) {
+            norm->data.F64[i] = vec->data.F64[i]*myPoly->scale[dir] + myPoly->zero[dir];
+        }
+        return norm;
+    }
+    if (vec->type.type == PS_TYPE_F32) {
+        for (int i = 0; i < vec->n; i++) {
+            norm->data.F64[i] = vec->data.F32[i]*myPoly->scale[dir] + myPoly->zero[dir];
+        }
+        return norm;
+    }
+
+    psError(PS_ERR_UNKNOWN, true, "invalid type for chebyshev polynomial");
+    return NULL;
+}
+
+/** 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
+    function returns the vector T^n (x_i) where x_i is the input vector of values and n is
+    the polynomial order.
+ */
+psVector *psChebyshevPolyVector (const psVector *vec, int order) {
+
+    if (order > 9) {
+        psWarning ("Chebyshev orders higher than 9 are not yet coded\n");
+        return NULL;
+    }
+
+    psVector *out = psVectorAlloc (vec->n, PS_TYPE_F64);
+
+    // easy but non-general implementation
+    switch (order) {
+      case 0:
+        for (int i = 0; i < vec->n; i++) {                                                   out->data.F64[i] = 1.0; } break;
+      case 1:
+        for (int i = 0; i < vec->n; i++) { psF64 x = vec->data.F64[i];                       out->data.F64[i] = x; } 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;
+      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;
+      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;
+      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;
+      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;
+      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;
+      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;
+      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;
+      default:
+        psWarning ("Chebyshev orders higher than 9 are not yet coded\n");
+        psFree (out);
+        return NULL;
+    }
+
+    return out;
+}
+
 /*****************************************************************************
     Polynomial coefficients will be accessed in [w][x][y][z] fashion.
@@ -234 +333 @@
 
 // XXX: You can do this without having to psAlloc() vector d.
-// XXX: How does the mask vector effect Crenshaw's formula?
+// 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.
@@ -343 +442 @@
                                   const psPolynomial2D* poly)
 {
-    PS_ASSERT_DOUBLE_WITHIN_RANGE(x, -1.0, 1.0, 0.0);
-    PS_ASSERT_DOUBLE_WITHIN_RANGE(y, -1.0, 1.0, 0.0);
+  // 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);
 
@@ -354 +454 @@
     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.
@@ -366 +469 @@
             if (!(poly->coeffMask[loop_x][loop_y] & PS_POLY_MASK_SET)) {
                 polySum += poly->coeff[loop_x][loop_y] *
-                           psPolynomial1DEval(chebPolys[loop_x], x) *
-                           psPolynomial1DEval(chebPolys[loop_y], y);
+                           psPolynomial1DEval(chebPolys[loop_x], xNorm) *
+                           psPolynomial1DEval(chebPolys[loop_y], yNorm);
             }
         }
@@ -603 +706 @@
     }
 
+    // scale & zero are used for Chebyshev polynomials to define the relationship between
+    // the independent variables and the normalized version with range -1 : +1.  These
+    // must be determined for a specific data set.
+    newPoly->scale[0] = NAN;
+    newPoly->zero[0]  = NAN;
+
     return(newPoly);
 }
@@ -638 +747 @@
             newPoly->coeffMask[x][y] = PS_POLY_MASK_NONE;
         }
+    }
+
+    // scale & zero are used for Chebyshev polynomials to define the relationship between
+    // the independent variables and the normalized version with range -1 : +1.  These
+    // must be determined for a specific data set.
+    for (int i = 0; i < 2; i++) {
+      newPoly->scale[i] = NAN;
+      newPoly->zero[i]  = NAN;
     }
 
@@ -756 +873 @@
             }
         }
+    }
+
+    // scale & zero are used for Chebyshev polynomials to define the relationship between
+    // the independent variables and the normalized version with range -1 : +1.  These
+    // must be determined for a specific data set.
+    for (int i = 0; i < 3; i++) {
+      newPoly->scale[i] = NAN;
+      newPoly->zero[i]  = NAN;
     }
 
@@ -820 +945 @@
     }
 
+    // scale & zero are used for Chebyshev polynomials to define the relationship between
+    // the independent variables and the normalized version with range -1 : +1.  These
+    // must be determined for a specific data set.
+    for (int i = 0; i < 4; i++) {
+      newPoly->scale[i] = NAN;
+      newPoly->zero[i]  = NAN;
+    }
+
     return(newPoly);
 }
@@ -841 +974 @@
 
 // this function must accept F32 and F64 input x vectors
+// EAM XXX these functions seem inefficiently implemented with many nested function calls.
+// they might benefit from unrolling.
 psVector *psPolynomial1DEvalVector(const psPolynomial1D *poly,
                                    const psVector *x)
@@ -888 +1023 @@
 }
 
+psVector *psPolynomial2DEvalChebVector(const psPolynomial2D *poly,
+                                       const psVector *x,
+                                       const psVector *y)
+{
+
+    if (!isfinite(poly->scale[0]) || !isfinite(poly->zero[0]) || !isfinite(poly->scale[1]) || !isfinite(poly->zero[1])) {
+        // re-calculate if not already determined?  
+        psError(PS_ERR_UNKNOWN, true, "normalization scales are not set for chebyshev polynomial");
+        return (NULL);
+    }
+
+    // Number of polynomial terms
+    int nXterm = 1 + poly->nX;      // Number of terms in x
+    int nYterm = 1 + poly->nY;      // Number of terms in y
+    if (nXterm > 9) {
+        psError(PS_ERR_UNKNOWN, false, "failed 2D chebyshev fit: orders higher than 9 are not yet coded\n");
+        return NULL;
+    }
+    if (nYterm > 9) {
+        psError(PS_ERR_UNKNOWN, false, "failed 2D chebyshev fit: orders higher than 9 are not yet coded\n");
+        return NULL;
+    }
+
+    // Generate normalized vectors for the range -1 : +1.  These functions cast to psF64
+    psVector *xNorm = psChebyshevNormVector (poly, x, 0);
+    psVector *yNorm = psChebyshevNormVector (poly, y, 1);
+    
+    // Generate the N cheb polynomials based on xNorm, yNorm
+    psArray *xPolySet = psArrayAlloc (nXterm);
+    for (int i = 0; i < nXterm; i++) {
+        xPolySet->data[i] = psChebyshevPolyVector (xNorm, i);
+    }
+    psArray *yPolySet = psArrayAlloc (nYterm);
+    for (int i = 0; i < nYterm; i++) {
+        yPolySet->data[i] = psChebyshevPolyVector (yNorm, i);
+    }
+
+    psVector *out = psVectorAlloc (x->n, PS_TYPE_F64);
+
+    psF64 *xData = xNorm->data.F64;
+    psF64 *yData = yNorm->data.F64;
+    psF64 *fData = out->data.F64;
+
+    // loop over all elements of the data vector
+    for (int i = 0; i < x->n; i++) {
+
+        if (!finite(xData[i])) {fData[i] = NAN; continue; }
+        if (!finite(yData[i])) {fData[i] = NAN; continue; }
+
+        psF64 sum = 0.0;
+        for (int jx = 0; jx < nXterm; jx++) {
+            psVector *jxCheb = xPolySet->data[jx];
+            for (int jy = 0; jy < nYterm; jy++) {
+                psVector *jyCheb = yPolySet->data[jy];
+                sum += poly->coeff[jx][jy] * jxCheb->data.F64[i] * jyCheb->data.F64[i];
+            }
+        }
+        fData[i] = sum;
+    }
+
+    psFree (xPolySet);
+    psFree (yPolySet);
+    psFree (xNorm);
+    psFree (yNorm);
+
+    return out;
+}
+
 // this function must support input data types of F32 and F64
 // all input vectors data types must match (all F32 or all F64)
@@ -901 +1104 @@
     PS_ASSERT_VECTOR_TYPE_F32_OR_F64(y, NULL);
 
-    psVector *tmp;
-    unsigned int vecLen=x->n;
+    unsigned int vecLen = x->n;
+
+    // input vector types must match
+    if (y->type.type != x->type.type) {
+        psError(PS_ERR_UNKNOWN, true, "type mismatch in data vectors");
+        return (NULL);
+    }
 
     // Determine the length of the output vector to by the minimum of the x,y vectors
-    if (y->n < vecLen) {
-        vecLen = y->n;
+    // XXX shouldn't we require x & y to have the same length?  seems meaningless otherwise
+    if (y->n != vecLen) {
+        psError(PS_ERR_UNKNOWN, true, "length mismatch in data vectors");
+        return (NULL);
+    }
+
+    if (poly->type == PS_POLYNOMIAL_CHEB) {
+        psVector *out = psPolynomial2DEvalChebVector (poly, x, y);
+        return out;
     }
 
     switch (x->type.type) {
-    case PS_TYPE_F32:
-        if (y->type.type != x->type.type) {
-            psError(PS_ERR_UNKNOWN, true, "type mismatch in data vectors");
-            return (NULL);
-        }
-
-        // Create output vector to return
-        tmp = psVectorAlloc(vecLen, PS_TYPE_F32);
-
-        // Evaluate the polynomial at the specified points
-        for (unsigned int i=0; i<vecLen; i++) {
-            tmp->data.F32[i] = psPolynomial2DEval(poly,x->data.F32[i],y->data.F32[i]);
-        }
-        break;
-    case PS_TYPE_F64:
-        if (y->type.type != x->type.type) {
-            psError(PS_ERR_UNKNOWN, true, "type mismatch in data vectors");
-            return (NULL);
-        }
-
-        // Create output vector to return
-        tmp = psVectorAlloc(vecLen, PS_TYPE_F64);
-
-        // Evaluate the polynomial at the specified points
-        for (unsigned int i=0; i<vecLen; i++) {
-            tmp->data.F64[i] = psPolynomial2DEval(poly,x->data.F64[i],y->data.F64[i]);
-        }
-        break;
-    default:
+      case PS_TYPE_F32: {
+          // Create output vector to return
+          psVector *out = psVectorAlloc(vecLen, PS_TYPE_F32);
+
+          // Evaluate the polynomial at the specified points
+          for (unsigned int i = 0; i < vecLen; i++) {
+              out->data.F32[i] = psPolynomial2DEval(poly,x->data.F32[i],y->data.F32[i]);
+          }
+          return out;
+      }
+      case PS_TYPE_F64: {
+          // Create output vector to return
+          psVector *out = psVectorAlloc(vecLen, PS_TYPE_F64);
+
+          // Evaluate the polynomial at the specified points
+          for (unsigned int i = 0; i < vecLen; i++) {
+              out->data.F64[i] = psPolynomial2DEval(poly,x->data.F64[i],y->data.F64[i]);
+          }
+          return out;
+      }
+      default:
         psError(PS_ERR_UNKNOWN, false, "invalid input data type.\n");
         return (NULL);
     }
-    // Return output vector
-    return(tmp);
+    psAbort ("impossible");
+    return NULL;
 }
 

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

@@ -68 +68 @@
     psF64 *coeffErr;                    ///< Error in coefficients
     psMaskType *coeffMask;              ///< Coefficient mask
+    double scale[1];                    ///< Chebyshev scale factor
+    double zero[1];                     ///< Chebyshev zero point
 }
 psPolynomial1D;
@@ -80 +82 @@
     psF64 **coeffErr;                   ///< Error in coefficients
     psMaskType **coeffMask;                  ///< Coefficients mask
+    double scale[2];                    ///< Chebyshev scale factor
+    double zero[2];                     ///< Chebyshev zero point
 }
 psPolynomial2D;
@@ -93 +97 @@
     psF64 ***coeffErr;                  ///< Error in coefficients
     psMaskType ***coeffMask;                 ///< Coefficients mask
+    double scale[3];                    ///< Chebyshev scale factor
+    double zero[3];                     ///< Chebyshev zero point
 }
 psPolynomial3D;
@@ -107 +113 @@
     psF64 ****coeffErr;                 ///< Error in coefficients
     psMaskType ****coeffMask;           ///< Coefficients mask
+    double scale[4];                    ///< Chebyshev scale factor
+    double zero[4];                     ///< Chebyshev zero point
 }
 psPolynomial4D;
@@ -305 +313 @@
 p_chebyPolys;
 
+
+// chebyshev support functions:
+bool psChebyshevSetScale (psPolynomial2D* myPoly, const psVector *vec, int dir);
+psVector *psChebyshevNormVector (const psPolynomial2D* myPoly, const psVector *vec, int dir);
+psVector *psChebyshevPolyVector (const psVector *vec, int order);
 
 /*****************************************************************************