Changeset 41896

Timestamp:

Nov 4, 2021, 6:10:51 PM (5 years ago)

Author:

eugene

Message:

merge changes from eam_branches/ipp-dev-20210817 (fix chebyshevs 1D, 2D, set DB field to NULL for inf, psFitsTableNew)

Location:

trunk/psLib

Files:

• src/astro/psCoord.c (modified) (9 diffs)
• src/astro/psCoord.h (modified) (1 diff)
• src/db/psDB.c (modified) (1 diff)
• src/fits/psFitsTableNew.c (modified) (2 diffs)
• src/fits/psFitsTableNew.h (modified) (2 diffs)
• src/math/psMinimizePolyFit.c (modified) (4 diffs)
• src/math/psPolynomial.c (modified) (12 diffs)
• src/math/psPolynomial.h (modified) (5 diffs)
• test/math/tap_psPolyFit2DCheb.c (copied) (copied from branches/eam_branches/ipp-dev-20210817/psLib/test/math/tap_psPolyFit2DCheb.c )

trunk/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);

trunk/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;
 

trunk/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;
     }

trunk/psLib/src/fits/psFitsTableNew.c

@@ -137 +137 @@
 
     return table;
+}
+
+psFitsTable *psFitsTableCreate (psArray *tableColumns, int numRows) {
+
+    psFitsTable *table = psFitsTableAlloc (tableColumns->n, numRows);
+
+    // now define the table colums
+    for (int i = 0; i < tableColumns->n; i++) {
+        psFitsTableColumn *column = &table->columns[i];
+
+        psFitsTableColumn *myColumn = tableColumns->data[i];
+
+        column->name = psStringCopy (myColumn->name);
+        column->type = myColumn->type;
+        column->vectorType = myColumn->vectorType;
+        column->elementSize = myColumn->elementSize;
+        
+        psMetadataAddS32(table->index, PS_LIST_TAIL, column->name, 0, NULL, i);
+
+        switch (column->type) {
+          case PS_DATA_S8:  column->data.S8  = psAlloc(sizeof(psS8) *table->numRows); break;
+          case PS_DATA_S16: column->data.S16 = psAlloc(sizeof(psS16)*table->numRows); break;
+          case PS_DATA_S32: column->data.S32 = psAlloc(sizeof(psS32)*table->numRows); break;
+          case PS_DATA_S64: column->data.S64 = psAlloc(sizeof(psS64)*table->numRows); break;
+          case PS_DATA_U8:  column->data.U8  = psAlloc(sizeof(psU8) *table->numRows); break;
+          case PS_DATA_U16: column->data.U16 = psAlloc(sizeof(psU16)*table->numRows); break;
+          case PS_DATA_U32: column->data.U32 = psAlloc(sizeof(psU32)*table->numRows); break;
+          case PS_DATA_U64: column->data.U64 = psAlloc(sizeof(psU64)*table->numRows); break;
+          case PS_DATA_F32: column->data.F32 = psAlloc(sizeof(psF32)*table->numRows); break;
+          case PS_DATA_F64: column->data.F64 = psAlloc(sizeof(psF64)*table->numRows); break;
+          default: break;
+        }
+    }
+    return table;
+}
+
+static void freeColumn(psFitsTableColumn *column)
+{
+    // add any element frees here
+    return;
+}
+
+psFitsTableColumn *psFitsTableColumnAlloc (char *name, psDataType type) {
+
+    psFitsTableColumn *column = psAlloc(sizeof(psFitsTableColumn));
+    psMemSetDeallocator(column, (psFreeFunc) freeColumn);
+
+    column->name = psStringCopy (name);
+    column->type = type;
+    column->vectorType  = 0;
+    column->elementSize = 1;
+    column->data.S8 = NULL; // no data yet allocated, place-holder only
+    return column;
+}
+
+bool psFitsTableColumnAdd (psArray *tableColumns, char *name, psDataType type) {
+    psFitsTableColumn *column = psFitsTableColumnAlloc (name, type);
+    psArrayAdd (tableColumns, 10, column);
+    return true;
 }
 
@@ -397 +456 @@
 psFitsTableGetNumTYPE(Bool, BOOL);
 
+#define psFitsTableSetNumTYPE(TYPE, NAME) \
+    bool psFitsTableSet##NAME(psFitsTable *table, int row, const char *name, ps##TYPE value) { \
+                                                                        \
+    if (row >= table->numRows || row < 0) { return false; }             \
+                                                                        \
+    bool mdok;                                                          \
+    int col = psMetadataLookupS32(&mdok, table->index, name);           \
+    if (!mdok) { return false; }                                        \
+    if (col < 0) { return false; }                                      \
+                                                                        \
+    psFitsTableColumn *column = &table->columns[col];                   \
+    switch (column->type) {                                             \
+      case PS_DATA_S8:                                                  \
+        column->data.S8[row] = (psS8) value;                            \
+        break;                                                          \
+      case PS_DATA_S16:                                                 \
+        column->data.S16[row] = (psS16) value;                          \
+        break;                                                          \
+      case PS_DATA_S32:                                                 \
+        column->data.S32[row] = (psS32) value;                          \
+        break;                                                          \
+      case PS_DATA_S64:                                                 \
+        column->data.S64[row] = (psS64) value;                          \
+        break;                                                          \
+      case PS_DATA_U8:                                                  \
+        column->data.U8[row] = (psU8) value;                            \
+        break;                                                          \
+      case PS_DATA_U16:                                                 \
+        column->data.U16[row] = (psU16) value;                          \
+        break;                                                          \
+      case PS_DATA_U32:                                                 \
+        column->data.U32[row] = (psU32) value;                          \
+        break;                                                          \
+      case PS_DATA_U64:                                                 \
+        column->data.U64[row] = (psU64) value;                          \
+        break;                                                          \
+      case PS_DATA_F32:                                                 \
+        column->data.F32[row] = (psF32) value;                          \
+        break;                                                          \
+      case PS_DATA_F64:                                                 \
+        column->data.F64[row] = (psF64) value;                          \
+        break;                                                          \
+      case PS_DATA_BOOL:                                                \
+        if (value) { column->data.F64[row] = 1; }                       \
+        else { column->data.F64[row] = 0; }                             \
+        break;                                                          \
+      default:                                                          \
+        /* if you get to this point, the value is not a number. */      \
+        return false;                                                   \
+        break;                                                          \
+    }                                                                   \
+    return true;                                                        \
+}
+
+psFitsTableSetNumTYPE(S8, S8);
+psFitsTableSetNumTYPE(U8, U8);
+psFitsTableSetNumTYPE(S16, S16);
+psFitsTableSetNumTYPE(U16, U16);
+psFitsTableSetNumTYPE(S32, S32);
+psFitsTableSetNumTYPE(U32, U32);
+psFitsTableSetNumTYPE(S64, S64);
+psFitsTableSetNumTYPE(U64, U64);
+psFitsTableSetNumTYPE(F32, F32);
+psFitsTableSetNumTYPE(F64, F64);
+psFitsTableSetNumTYPE(Bool, BOOL);
+
 // remove rows from table that are marked in the mask array as censored.
 bool psFitsTableCensor(psFitsTable *table, bool *censorMask)

trunk/psLib/src/fits/psFitsTableNew.h

@@ -55 +55 @@
 // Get value for given row and column name
 psBool psFitsTableGetBool(bool *status, psFitsTable *table, int row, const char* name);
-psS8 psFitsTableGetS8(bool *status, psFitsTable *table, int row, const char* name);
-psU8 psFitsTableGetU8(bool *status, psFitsTable *table, int row, const char* name);
+psS8  psFitsTableGetS8(bool *status, psFitsTable *table, int row, const char* name);
+psU8  psFitsTableGetU8(bool *status, psFitsTable *table, int row, const char* name);
 psS16 psFitsTableGetS16(bool *status, psFitsTable *table, int row, const char* name);
 psU16 psFitsTableGetU16(bool *status, psFitsTable *table, int row, const char* name);
@@ -63 +63 @@
 psS64 psFitsTableGetS64(bool *status, psFitsTable *table, int row, const char* name);
 psU64 psFitsTableGetU64(bool *status, psFitsTable *table, int row, const char* name);
-
 psF32 psFitsTableGetF32(bool *status, psFitsTable *table, int row, const char* name);
 psF64 psFitsTableGetF64(bool *status, psFitsTable *table, int row, const char* name);
 
+psFitsTable *psFitsTableCreate (psArray *tableColumns, int numRows);
+bool psFitsTableColumnAdd (psArray *tableColumns, char *name, psDataType type);
+psFitsTableColumn *psFitsTableColumnAlloc (char *name, psDataType type);
+
+bool psFitsTableSetS8 (psFitsTable *table, int row, const char* name, psS8  value);
+bool psFitsTableSetU8 (psFitsTable *table, int row, const char* name, psU8  value);
+bool psFitsTableSetS16(psFitsTable *table, int row, const char* name, psS16 value);
+bool psFitsTableSetU16(psFitsTable *table, int row, const char* name, psU16 value);
+bool psFitsTableSetS32(psFitsTable *table, int row, const char* name, psS32 value);
+bool psFitsTableSetU32(psFitsTable *table, int row, const char* name, psU32 value);
+bool psFitsTableSetS64(psFitsTable *table, int row, const char* name, psS64 value);
+bool psFitsTableSetU64(psFitsTable *table, int row, const char* name, psU64 value);
+bool psFitsTableSetF32(psFitsTable *table, int row, const char* name, psF32 value);
+bool psFitsTableSetF64(psFitsTable *table, int row, const char* name, psF64 value);
+
 #endif

trunk/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");

trunk/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;
+}
+
+# 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
+    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++) { CHEB_EVAL_0(out->data.F64[i], vec->data.F64[i]); } break;
+      case 1:
+        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++) { CHEB_EVAL_2(out->data.F64[i], vec->data.F64[i]); } break;
+      case 3:
+        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++) { CHEB_EVAL_4(out->data.F64[i], vec->data.F64[i]); } break;
+      case 5:
+        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++) { CHEB_EVAL_6(out->data.F64[i], vec->data.F64[i]); } break;
+      case 7:
+        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++) { CHEB_EVAL_8(out->data.F64[i], vec->data.F64[i]); } break;
+      case 9:
+        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");
+        psFree (out);
+        return NULL;
+    }
+
+    return out;
+}
+
 /*****************************************************************************
     Polynomial coefficients will be accessed in [w][x][y][z] fashion.
@@ -233 +343 @@
 }
 
-// XXX: You can do this without having to psAlloc() vector d.
-// XXX: How does the mask vector effect Crenshaw'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.
@@ -268 +401 @@
     }
 
-    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,
@@ -343 +480 @@
                                   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);
     PS_ASSERT_POLY_NON_NULL(poly, NAN);
 
-    unsigned int loop_x = 0;
-    unsigned int loop_y = 0;
-    unsigned int i = 0;
+    psF64 xNorm = x*poly->scale[0] + poly->zero[0];
+    psF64 yNorm = y*poly->scale[1] + poly->zero[1];
+
     psF64 polySum = 0.0;
-    psPolynomial1D* *chebPolys = NULL;
-    unsigned int maxChebyPoly = 0;
-
-    // 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], x) *
-                           psPolynomial1DEval(chebPolys[loop_y], y);
-            }
-        }
-    }
-    for (i=0;i<maxChebyPoly+1;i++) {
-        psFree(chebPolys[i]);
-    }
-    psFree(chebPolys);
+
+    // XXX this could be quicker if we saved the N xvalues are re-used the resuls
+    for (int ix = 0; ix <= poly->nX; 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);
 }
@@ -603 +752 @@
     }
 
+    // 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 +793 @@
             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 +919 @@
             }
         }
+    }
+
+    // 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 +991 @@
     }
 
+    // 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);
 }
 
+/* note these functions accept unscaled values and apply the scaling saved on poly */
 psF64 psPolynomial1DEval(const psPolynomial1D* poly,
                          psF64 x)
@@ -830 +1010 @@
     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);
 }
 
 // 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)
@@ -878 +1061 @@
     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);
+}
+
+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];
+                if (poly->coeffMask[jx][jy] & PS_POLY_MASK_SET) continue;
+                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;
 }
 
@@ -901 +1153 @@
     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;
 }
 

trunk/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);
 
 /*****************************************************************************