Index: /trunk/psLib/src/astro/psCoord.c =================================================================== --- /trunk/psLib/src/astro/psCoord.c (revision 41895) +++ /trunk/psLib/src/astro/psCoord.c (revision 41896) @@ -40,4 +40,10 @@ #include "psMinimizePolyFit.h" +# define TEST_SAVE_INVERSE_TRANSFORM 0 + +# if (TEST_SAVE_INVERSE_TRANSFORM) +# include "psBinaryOp.h" +# endif + # define ELIXIR_CODE 1 @@ -86,5 +92,6 @@ } - 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,5 +198,5 @@ } -psPlaneTransform* psPlaneTransformAlloc(int order1, int order2) +psPlaneTransform* psPlaneTransformAlloc(int order1, int order2, psPolynomialType type) { PS_ASSERT_INT_NONNEGATIVE(order1, NULL); @@ -198,6 +205,6 @@ 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,4 +804,10 @@ } + // 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,5 +827,5 @@ psPlaneTransform *myPT = NULL; if (out == NULL) { - myPT = psPlaneTransformAlloc(orderX, orderY); + myPT = psPlaneTransformAlloc(orderX, orderY, PS_POLYNOMIAL_ORD); } else { if ((out->x->nX == orderX) && @@ -834,5 +847,5 @@ } else { psFree(out); - myPT = psPlaneTransformAlloc(orderX, orderY); + myPT = psPlaneTransformAlloc(orderX, orderY, PS_POLYNOMIAL_ORD); } } @@ -1037,5 +1050,5 @@ 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,4 +1084,59 @@ 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); Index: /trunk/psLib/src/astro/psCoord.h =================================================================== --- /trunk/psLib/src/astro/psCoord.h (revision 41895) +++ /trunk/psLib/src/astro/psCoord.h (revision 41896) @@ -178,6 +178,7 @@ 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; Index: /trunk/psLib/src/db/psDB.c =================================================================== --- /trunk/psLib/src/db/psDB.c (revision 41895) +++ /trunk/psLib/src/db/psDB.c (revision 41896) @@ -2875,37 +2875,36 @@ 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; } Index: /trunk/psLib/src/fits/psFitsTableNew.c =================================================================== --- /trunk/psLib/src/fits/psFitsTableNew.c (revision 41895) +++ /trunk/psLib/src/fits/psFitsTableNew.c (revision 41896) @@ -137,4 +137,63 @@ 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,4 +456,70 @@ 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) Index: /trunk/psLib/src/fits/psFitsTableNew.h =================================================================== --- /trunk/psLib/src/fits/psFitsTableNew.h (revision 41895) +++ /trunk/psLib/src/fits/psFitsTableNew.h (revision 41896) @@ -55,6 +55,6 @@ // 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,7 +63,21 @@ 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 Index: /trunk/psLib/src/math/psMinimizePolyFit.c =================================================================== --- /trunk/psLib/src/math/psMinimizePolyFit.c (revision 41895) +++ /trunk/psLib/src/math/psMinimizePolyFit.c (revision 41896) @@ -35,4 +35,5 @@ #include "psMinimizePolyFit.h" +#include "psAbort.h" #include "psAssert.h" #include "psMinimizeLMM.h" // For Gauss-Jordan routines @@ -1225,4 +1226,153 @@ /****************************************************************************** +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,5 +1390,5 @@ { 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,10 +1427,15 @@ 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"); Index: /trunk/psLib/src/math/psPolynomial.c =================================================================== --- /trunk/psLib/src/math/psPolynomial.c (revision 41895) +++ /trunk/psLib/src/math/psPolynomial.c (revision 41896) @@ -36,4 +36,5 @@ #include "psLogMsg.h" #include "psPolynomial.h" +#include "psAbort.h" #include "psAssert.h" @@ -203,4 +204,113 @@ +/** 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,19 +343,42 @@ } -// 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,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, @@ -343,36 +480,48 @@ 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;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); } @@ -603,4 +752,10 @@ } + // 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,4 +793,12 @@ 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,4 +919,12 @@ } } + } + + // 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,7 +991,16 @@ } + // 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,15 +1010,18 @@ 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,12 +1061,81 @@ 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,47 +1153,51 @@ 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; idata.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; idata.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; } Index: /trunk/psLib/src/math/psPolynomial.h =================================================================== --- /trunk/psLib/src/math/psPolynomial.h (revision 41895) +++ /trunk/psLib/src/math/psPolynomial.h (revision 41896) @@ -68,4 +68,6 @@ psF64 *coeffErr; ///< Error in coefficients psMaskType *coeffMask; ///< Coefficient mask + double scale[1]; ///< Chebyshev scale factor + double zero[1]; ///< Chebyshev zero point } psPolynomial1D; @@ -80,4 +82,6 @@ psF64 **coeffErr; ///< Error in coefficients psMaskType **coeffMask; ///< Coefficients mask + double scale[2]; ///< Chebyshev scale factor + double zero[2]; ///< Chebyshev zero point } psPolynomial2D; @@ -93,4 +97,6 @@ psF64 ***coeffErr; ///< Error in coefficients psMaskType ***coeffMask; ///< Coefficients mask + double scale[3]; ///< Chebyshev scale factor + double zero[3]; ///< Chebyshev zero point } psPolynomial3D; @@ -107,4 +113,6 @@ psF64 ****coeffErr; ///< Error in coefficients psMaskType ****coeffMask; ///< Coefficients mask + double scale[4]; ///< Chebyshev scale factor + double zero[4]; ///< Chebyshev zero point } psPolynomial4D; @@ -305,4 +313,9 @@ 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); /***************************************************************************** Index: /trunk/psLib/test/math/tap_psPolyFit2DCheb.c =================================================================== --- /trunk/psLib/test/math/tap_psPolyFit2DCheb.c (revision 41896) +++ /trunk/psLib/test/math/tap_psPolyFit2DCheb.c (revision 41896) @@ -0,0 +1,84 @@ +#include +#include +#include +//#include "tap.h" +//#include "pstap.h" + +int main(int argc, char **argv) +{ + psLogSetFormat("HLNM"); + psLogSetLevel(PS_LOG_INFO); + // plan_tests(88); + + if (argc != 4) { + fprintf (stderr, "USAGE: tap_psPolyFit2DCheb (input) (Nx) (Ny)\n"); + exit (2); + } + + FILE *f = fopen (argv[1], "r"); + int Nx = atoi(argv[2]); + int Ny = atoi(argv[3]); + + // data file should contain x y p + // determine the best fit 2D polynomial from x y to p + psVector *x = psVectorAllocEmpty(1000, PS_TYPE_F64); + psVector *y = psVectorAllocEmpty(1000, PS_TYPE_F64); + psVector *p = psVectorAllocEmpty(1000, PS_TYPE_F64); + + double xv, yv, pv; + + while (fscanf (f, "%lf %lf %lf", &xv, &yv, &pv) == 3) { + psVectorAppend (x, xv); + psVectorAppend (y, yv); + psVectorAppend (p, pv); + } + + // Fit an ordinary polynomial + { + psPolynomial2D *polyOrd = psPolynomial2DAlloc(PS_POLYNOMIAL_ORD, Nx, Ny); + int result = psVectorFitPolynomial2D(polyOrd, NULL, 0, p, NULL, x, y); + + fprintf (stdout, " ---- Ord ---- \n"); + + for (int ix = 0; ix < polyOrd->nX + 1; ix++) { + for (int iy = 0; iy < polyOrd->nY + 1; iy++) { + fprintf (stdout, "%18.12e ", polyOrd->coeff[ix][iy]); + } + fprintf (stdout, "\n"); + } + fprintf (stdout, "\n"); + + psVector *P = psPolynomial2DEvalVector (polyOrd, x, y); + + FILE *g1 = fopen ("ord.dat", "w"); + for (int i = 0; i < x->n; i++) { + fprintf (g1, "%lf %lf %lf %lf\n", x->data.F64[i], y->data.F64[i], p->data.F64[i], P->data.F64[i]); + } + fclose (g1); + } + + // Fit a chebychev polynomial + { + psPolynomial2D *polyCheb = psPolynomial2DAlloc(PS_POLYNOMIAL_CHEB, Nx, Ny); + + int result = psVectorFitPolynomial2D(polyCheb, NULL, 0, p, NULL, x, y); + + fprintf (stdout, " ---- Cheb ---- \n"); + + for (int ix = 0; ix < polyCheb->nX + 1; ix++) { + for (int iy = 0; iy < polyCheb->nY + 1; iy++) { + fprintf (stdout, "%18.12e ", polyCheb->coeff[ix][iy]); + } + fprintf (stdout, "\n"); + } + + psVector *P = psPolynomial2DEvalVector (polyCheb, x, y); + + FILE *g1 = fopen ("cheb.dat", "w"); + for (int i = 0; i < x->n; i++) { + fprintf (g1, "%lf %lf %lf %lf\n", x->data.F64[i], y->data.F64[i], p->data.F64[i], P->data.F64[i]); + } + fclose (g1); + } +} +/// gcc -o tap_psPolyFit2DCheb tap_psPolyFit2DCheb.c -lm -lpslib -I /home/real/eugene/src/psconfig/ipp-trunk-20210518.lin64/include/pslib -L /home/real/eugene/src/psconfig/ipp-trunk-20210518.lin64/lib