Changeset 42510

Timestamp:

Aug 17, 2023, 3:24:45 PM (3 years ago)

Author:

eugene

Message:

convert 2D fitting to IRLS-based fitting; use sample median instead of robust median for initial pass

File:

• branches/eam_branches/ipp-20230313/psphot/src/psphotMakeResiduals.c (modified) (13 diffs)

branches/eam_branches/ipp-20230313/psphot/src/psphotMakeResiduals.c

@@ -95 +95 @@
     psVector *yC = psVectorAllocEmpty (100, PS_TYPE_F32);
 
-# define TESTRESID 1
+# define TESTRESID 0
 # if (TESTRESID)
     psArray *sourceRawCube = psArrayAllocEmpty (1);
@@ -101 +101 @@
     psArray *sourceVarCube = psArrayAllocEmpty (1);
     psArray *sourceMskCube = psArrayAllocEmpty (1);
+
+    psImage *refImage = NULL;
 # endif
 
@@ -121 +123 @@
 # if (TESTRESID)
         psImage *tmpImage = psImageCopy(NULL, image, PS_TYPE_F32);
-        psArrayAdd (sourceRawCube, 1, tmpImage);
-        psFree (tmpImage);
+        if (refImage == NULL) {
+          refImage = tmpImage;
+          psArrayAdd (sourceRawCube, 1, tmpImage);
+        } else {
+          if (tmpImage->numCols != refImage->numCols) {
+            fprintf (stderr, "mismatched image sizes (%d : %f, %f)\n", i, model->params->data.F32[PM_PAR_XPOS], model->params->data.F32[PM_PAR_YPOS]);
+            continue;
+          }
+          if (tmpImage->numRows != refImage->numRows) {
+            fprintf (stderr, "mismatched image sizes (%d : %f, %f)\n", i, model->params->data.F32[PM_PAR_XPOS], model->params->data.F32[PM_PAR_YPOS]);
+            continue;
+          }
+          psArrayAdd (sourceRawCube, 1, tmpImage);
+          psFree (tmpImage);
+        }
 # endif
 
@@ -199 +214 @@
     psVector *fmasks  = psVectorAlloc (input->n, PS_TYPE_VECTOR_MASK);
 
-    // statistic to use to determine baseline for clipping
-    psStats *fluxClip     = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
-    psStats *fluxClipDef  = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
     // statistic to use to determine output flux
     // XXX make API to convert statOption for MEAN/MEDIAN in to corresponding STDEV?
@@ -207 +219 @@
     psStats *fluxStatsDef = psStatsAlloc (statOption | PS_STAT_SAMPLE_STDEV);
 
+    // use this for the IRLS fitting below (defines niter)
+    psStats *statsIRLS = psStatsAlloc(PS_STAT_CLIPPED_MEAN); statsIRLS->clipIter = 10; // max number of iterations
+    psPolynomial2D *polyIRLS = psPolynomial2DAlloc(PS_POLYNOMIAL_ORD, 1, 1);
+
     // Use psF64 to minimize overflow problems?
     psImage *A = psImageAlloc(3, 3, PS_TYPE_F64); // Least-squares matrix
@@ -217 +233 @@
 
             int nGoodPixel = 0;              // pixel is off the image
+
+            // skip pixels outside of the requested radius
+            float radius = hypot((ox - 0.5*resid->Ro->numCols), (oy - 0.5*resid->Ro->numRows));
+            if (radius > radiusMax) {
+              resid->mask->data.PM_TYPE_RESID_MASK_DATA[oy][ox] = 1;
+              continue;
+            }
 
             // build the vector of data values for this output pixel
@@ -268 +291 @@
 
             // measure the robust median to determine a baseline reference value
-            *fluxClip = *fluxClipDef;
+            psStats *fluxClip = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
             if (!psVectorStats (fluxClip, fluxes, NULL, fmasks, fmaskVal)) {
                 psError(PSPHOT_ERR_CONFIG, false, "Error calculating residual stats");
+                psFree (fluxClip);
                 return false;
             }
-            if (isnan(fluxClip->robustMedian)) {
+            if (isnan(fluxClip->sampleMedian)) {
                 resid->Ro->data.F32[oy][ox] = 0.0;
                 resid->Rx->data.F32[oy][ox] = 0.0;
                 resid->Ry->data.F32[oy][ox] = 0.0;
                 resid->mask->data.PM_TYPE_RESID_MASK_DATA[oy][ox] = badMask;
+                psFree (fluxClip);
                 continue;
             }
 
-            if (fabs(fluxClip->robustMedian) > 1.5) {
-              fprintf (stderr, "warning: %d, %d : residual pixel has funny initial median value: %f\n", ox, oy, fluxClip->robustMedian);
+            if (fabs(fluxClip->sampleMedian) > 1.5) {
+              fprintf (stderr, "warning: %d, %d : residual pixel has funny initial median value: %f\n", ox, oy, fluxClip->sampleMedian);
             }
 
@@ -288 +313 @@
             int nKeep = 0;
             for (int i = 0; i < fluxes->n; i++) {
-                float delta = fluxes->data.F32[i] - fluxClip->robustMedian;
+                float delta = fluxes->data.F32[i] - fluxClip->sampleMedian;
                 float sigma = sqrt (dfluxes->data.F32[i]);
                 float swing = fabs(delta) / sigma;
@@ -298 +323 @@
                 if (!fmasks->data.PS_TYPE_VECTOR_MASK_DATA[i]) nKeep++;
             }
+            psFree (fluxClip);
 
             if (nKeep < 5) {
@@ -311 +337 @@
                 }
 
-                // XXX this test should go at the top of the loop
-                float radius = hypot((ox - 0.5*resid->Ro->numCols), (oy - 0.5*resid->Ro->numRows));
-                if (radius > radiusMax) {
-                  resid->mask->data.PM_TYPE_RESID_MASK_DATA[oy][ox] = 1;
-                  continue;
-                }
+                resid->Rx->data.F32[oy][ox] = resid->Ry->data.F32[oy][ox] = 0.0;
 
                 resid->Ro->data.F32[oy][ox] = psStatsGetValue(fluxStats, statOption);
-                resid->Rx->data.F32[oy][ox] = resid->Ry->data.F32[oy][ox] = 0.0;
-
                 if (isnan(resid->Ro->data.F32[oy][ox])) {
                     resid->Ro->data.F32[oy][ox] = 0.0;
-                    resid->Rx->data.F32[oy][ox] = 0.0;
-                    resid->Ry->data.F32[oy][ox] = 0.0;
                     resid->mask->data.PM_TYPE_RESID_MASK_DATA[oy][ox] = badMask;
                     continue;
@@ -335 +352 @@
                 assert (SPATIAL_ORDER == 1);
 
-                float radius = hypot((ox - 0.5*resid->Ro->numCols), (oy - 0.5*resid->Ro->numRows));
-                if (radius > radiusMax) {
-                  resid->mask->data.PM_TYPE_RESID_MASK_DATA[oy][ox] = 1;
-                  continue;
-                }
-
-                psImageInit(A, 0.0);
-                psVectorInit(B, 0.0);
-                for (int i = 0; i < fluxes->n; i++) {
-                    if (fmasks->data.PS_TYPE_VECTOR_MASK_DATA[i]) continue;
-                    B->data.F64[0] += fluxes->data.F32[i]/dfluxes->data.F32[i];
-                    B->data.F64[1] += fluxes->data.F32[i]*xC->data.F32[i]/dfluxes->data.F32[i];
-                    B->data.F64[2] += fluxes->data.F32[i]*yC->data.F32[i]/dfluxes->data.F32[i];
-
-                    A->data.F64[0][0] += 1.0/dfluxes->data.F32[i];
-                    A->data.F64[1][0] += xC->data.F32[i]/dfluxes->data.F32[i];
-                    A->data.F64[2][0] += yC->data.F32[i]/dfluxes->data.F32[i];
-
-                    A->data.F64[1][1] += PS_SQR(xC->data.F32[i])/dfluxes->data.F32[i];
-                    A->data.F64[2][2] += PS_SQR(yC->data.F32[i])/dfluxes->data.F32[i];
-                    A->data.F64[1][2] += xC->data.F32[i]*yC->data.F32[i]/dfluxes->data.F32[i];
-                }
-
-                A->data.F64[0][1] = A->data.F64[1][0];
-                A->data.F64[0][2] = A->data.F64[2][0];
-                A->data.F64[2][1] = A->data.F64[1][2];
-
-                if (!psMatrixGJSolve(A, B)) {
-                    resid->mask->data.PM_TYPE_RESID_MASK_DATA[oy][ox] = 1;
-                    psWarning("Singular matrix solving for (y,x) = (%d,%d)'s residuals, masking", oy, ox);
-                    continue;
-                }
-
-                resid->Ro->data.F32[oy][ox] = B->data.F64[0];
-                resid->Rx->data.F32[oy][ox] = B->data.F64[1];
-                resid->Ry->data.F32[oy][ox] = B->data.F64[2];
-
-                float dRo = sqrt(A->data.F32[0][0]);
-
-                if (fabs(B->data.F64[0]) > 1.5) {
-                  fprintf (stderr, "warning: %d, %d : residual pixel has funny fit value: %f\n", ox, oy, B->data.F64[0]);
+                // we have the following vectors to describe this pixel:
+                // fluxes, dfluxes, fmasks, xC, yC
+
+                psVectorIRLSFitPolynomial2D (polyIRLS, statsIRLS, fmasks, fmaskVal, fluxes, dfluxes, xC, yC);
+
+                resid->Ro->data.F32[oy][ox] = polyIRLS->coeff[0][0];
+                resid->Rx->data.F32[oy][ox] = polyIRLS->coeff[1][0];
+                resid->Ry->data.F32[oy][ox] = polyIRLS->coeff[0][1];
+
+                // is the error well-defined?
+                float dRo = polyIRLS->coeffErr[0][0];
+
+                if (fabs(resid->Ro->data.F32[oy][ox]) > 1.5) {
+                  fprintf (stderr, "warning: %d, %d : residual pixel has funny fit value: %f\n", ox, oy, resid->Ro->data.F32[oy][ox]);
                 }
 
@@ -442 +432 @@
     psFree (fluxStats);
     psFree (fluxStatsDef);
-    psFree (fluxClip);
-    psFree (fluxClipDef);
 
     if (resid != NULL && psTraceGetLevel("psphot") > 5) {
@@ -453 +441 @@
     }
 
+    psFree (statsIRLS);
+    psFree (polyIRLS);
+
     psf->residuals = resid;
     return (resid != NULL) ? true : false;