Changeset 35455 for trunk/psModules/src/imcombine/pmSubtraction.c

Timestamp:

Apr 30, 2013, 6:10:41 PM (13 years ago)

Author:

watersc1

Message:

Changes to stacking code:

• New input FWHM clipping algorithm using mixture models
• New recipes/reductions that further reduce FWHM clipping
• pmSubtractionRejectStamps now fits flux/chi^{2 values in log-log space to prevent outliers biasing the fit (and causing images to be rejected due to a small number of bad stamps)}
• pmPSF now correctly chooses between useReff forms as appropriate
• pmPSF change fixes circularization issues in pmReadoutFakeFromSources
• pmPSFEnvelope somewhat more clear about circularization
• Changed debugging in pmSubtractionEquation to add information.

I do not see any conflicts with the ppSub code.

File:

• trunk/psModules/src/imcombine/pmSubtraction.c (modified) (3 diffs)

trunk/psModules/src/imcombine/pmSubtraction.c

@@ -37 +37 @@
 #define USE_KERNEL_ERR                  // Use kernel error image?
 #define NUM_COVAR_POS 5                 // Number of positions for covariance calculation
+#define USE_LOGFIT_REJECT
 
 // XXX we need to pass these fwhm values elsewhere.  These should go on one of the structure, but 
@@ -966 +967 @@
     psPolynomial1D *model = psPolynomial1DAlloc (PS_POLYNOMIAL_ORD, 2);
 
+#ifdef USE_LOGFIT_REJECT
+    // CZW: Since flux and chisq can range over many orders of magnitude, use a log-log fit to calculate the
+    // RMS scatter.  This should be more robust against outliers that can impact the fit quality.
+    psVector *logchisq = psVectorAlloc(match->chisq->n,PS_TYPE_F32);
+    psVector *logflux  = psVectorAlloc(match->fluxes->n,PS_TYPE_F32);
+    for (int qq = 0; qq < match->chisq->n; qq++) {
+      if ((match->chisq->data.F32[qq] > 0)&&
+          (match->fluxes->data.F32[qq] > 0)) {
+        logchisq->data.F32[qq] = log10(match->chisq->data.F32[qq]);
+        logflux->data.F32[qq]  = log10(match->fluxes->data.F32[qq]);
+      }
+      else {
+        // Ignore negative values.  Maybe we could do an offset or something, but not today.
+        match->stampMask->data.PS_TYPE_VECTOR_MASK_DATA[qq] |= PM_SUBTRACTION_STAMP_REJECTED;
+      } 
+      //      if (!(match->stampMask->data.PS_TYPE_VECTOR_MASK_DATA[qq] & 0xff)) {
+      //        psLogMsg("psModules.imcombine", PS_LOG_INFO, "RRRRRR: %d %g %g %g %g\n",
+      //                 qq,match->chisq->data.F32[qq],match->fluxes->data.F32[qq],
+      //                 0.0,0.0);
+      // }
+
+    }
+    bool result = psVectorClipFitPolynomial1D(model, stats, match->stampMask, 0xff, logchisq, NULL, logflux);
+    if (!result) {
+        psError(PM_ERR_DATA, false, "Unable to measure statistics for deviations.");
+        psFree(model);
+        psFree(stats);
+        return -1;
+    }
+
+    // However, since we've done a log-log fit, we can't rely on the statistics in the stats object
+    // To accurately represent the distribution.  I've tried to massage them (as dlog10(X) = abs(1/(X log(10))) dX),
+    // but ended up deciding to just manually calculate the residual, and then pass that to a robust stats object.
+    psStats *residStats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
+    psVector *residData = psVectorAlloc(match->chisq->n,PS_TYPE_F32);
+    int counter = 0;
+    for (int qq = 0; qq < match->chisq->n; qq++) {
+      if (!(match->stampMask->data.PS_TYPE_VECTOR_MASK_DATA[qq] & 0xff)) {
+        double vv = pow(10,
+                        model->coeff[0] +
+                        model->coeff[1] * log10(match->fluxes->data.F32[qq]) +
+                        model->coeff[2] * pow(log10(match->fluxes->data.F32[qq]),2));
+        residData->data.F32[qq] = match->chisq->data.F32[qq] - vv;
+        counter++;
+        //                      psLogMsg("psModules.imcombine", PS_LOG_INFO, "SSSSS: %d %g %g %g %g\n",
+        //                       qq,match->chisq->data.F32[qq],match->fluxes->data.F32[qq],
+        //                       vv,residData->data.F32[qq]);
+      }
+    }
+    psVectorStats(residStats,residData,NULL,match->stampMask,0xff);
+    if (isnan(residStats->robustMedian)) {
+        psError(PM_ERR_DATA, false, "Unable to measure statistics for deviations.");
+        psFree(model);
+        psFree(stats);
+        psFree(logchisq);
+        psFree(logflux);
+        psFree(residStats);
+        psFree(residData);
+
+        return -1;
+    }
+
+    kernels->mean = residStats->robustMedian;
+    kernels->rms =  residStats->robustStdev;
+    kernels->numStamps = counter;
+
+    psLogMsg ("pmPSFtry", 4, "chisq vs flux resid: %f +/- %f\n", residStats->robustMedian,residStats->robustStdev);
+    psFree(logchisq);
+    psFree(logflux);
+    psFree(residStats);
+    psFree(residData);
+
+#else
+    // CZW: Otherwise, use the original linear fit code.
     bool result = psVectorClipFitPolynomial1D(model, stats, match->stampMask, 0xff, match->chisq, NULL, match->fluxes);
     if (!result) {
@@ -979 +1054 @@
         return -1;
     }
-
     kernels->mean = stats->sampleMean;
-    kernels->rms = stats->sampleStdev;
+    kernels->rms =  stats->sampleStdev;
     kernels->numStamps = stats->clippedNvalues;
+    psLogMsg ("pmPSFtry", 4, "chisq vs flux resid: %f +/- %f\n", stats->sampleMean, stats->sampleStdev);
+#endif
 
     psLogMsg ("pmPSFtry", 4, "chisq vs flux model: %e + %e flux + %e flux^2\n", model->coeff[0], model->coeff[1], model->coeff[2]);
-    psLogMsg ("pmPSFtry", 4, "chisq vs flux resid: %f +/- %f\n", stats->sampleMean, stats->sampleStdev);
     psLogMsg("psModules.imcombine", PS_LOG_INFO, "Mean deviation from %d stamps: %lf +/- %lf",  kernels->numStamps, kernels->mean, kernels->rms);