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

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

File:

• branches/eam_branches/ipp-dev-20210817/psLib/src/astro/psCoord.c (modified) (9 diffs)

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