Changeset 35455 for trunk/psModules

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.

Location:

trunk/psModules/src

Files:

• camera/pmReadoutFake.c (modified) (1 diff)
• imcombine/pmPSFEnvelope.c (modified) (2 diffs)
• imcombine/pmStack.c (modified) (2 diffs)
• imcombine/pmSubtraction.c (modified) (3 diffs)
• imcombine/pmSubtractionEquation.c (modified) (2 diffs)
• objects/pmPSF.c (modified) (3 diffs)

trunk/psModules/src/camera/pmReadoutFake.c

@@ -152 +152 @@
             return false;
         }
-
+        
         psTrace("psModules.camera", 10, "Adding source at %f,%f with flux %f\n",
                 fakeModel->params->data.F32[PM_PAR_XPOS], fakeModel->params->data.F32[PM_PAR_YPOS],

trunk/psModules/src/imcombine/pmPSFEnvelope.c

@@ -37 +37 @@
 
 
-//#define TESTING                         // Enable test output
+// #define TESTING                         // Enable test output
 // #define PEAK_NORM                       // Normalise peaks?
 #define PEAK_FLUX 1.0e4                 // Peak flux for each source
@@ -190 +190 @@
         psf->residuals = NULL;
         pmModelClassSetLimits(PM_MODEL_LIMITS_MODERATE);
+        psLogMsg("psModules",PS_LOG_INFO,"Matching Input %d",i);
+#define CIRCULARIZE true
         if (!pmReadoutFakeFromSources(fakeRO, fakeSize, fakeSize, fakes, 0, xOffset, yOffset, psf,
-                                      NAN, radius, true, false)) {
+                                      NAN, radius, CIRCULARIZE, false)) {
             psError(PS_ERR_UNKNOWN, false, "Unable to generate fake readout.");
             psFree(envelope);

trunk/psModules/src/imcombine/pmStack.c

@@ -41 +41 @@
 /* #define TEST_X 5745                       // x coordinate to examine */
 /* #define TEST_Y 5331                       // y coordinate to examine */
-#define TEST_X 25
-#define TEST_Y 25
+// #define TEST_X 972
+// #define TEST_Y 3213
+#define TEST_X 3289
+#define TEST_Y 4810
 #define TEST_RADIUS 0.5                 // Radius to examine
 # endif
@@ -752 +754 @@
     if (PS_SQR(x - TEST_X) + PS_SQR(y - TEST_Y) <= PS_SQR(TEST_RADIUS)) {
         for (int i = 0; i < numGood; i++) {
-            fprintf(stderr,"Input %d, pixel %d,%d (%" PRIu16 "): %f %f (%f) %f %f %d %x %x -> %x %x\n",
-                    i, x, y, pixelSources->data.U16[i], pixelData->data.F32[i], pixelVariances->data.F32[i],
-                    addVariance->data.F32[i], pixelWeights->data.F32[i], pixelExps->data.F32[i],
-                    pixelSuspects->data.U8[i], badMaskBits, suspectMaskBits, *badMask, *goodMask);
+            fprintf(stderr,"Input %d, pixel %d,%d (%" PRIu16 "): %f %f (%g) %g %f %d %x %x -> %x %x\n",
+                    i, x, y, pixelSources->data.U16[i],
+                    pixelData->data.F32[i], pixelVariances->data.F32[i],
+                    addVariance->data.F32[i],
+                    pixelWeights->data.F32[i], pixelExps->data.F32[i],
+                    pixelSuspects->data.U8[i],
+                    badMaskBits, suspectMaskBits,
+                    *badMask, *goodMask);
         }
     }

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

trunk/psModules/src/imcombine/pmSubtractionEquation.c

@@ -751 +751 @@
     stamps->normValue2 = stats->robustMedian;
 
-    psLogMsg ("psModules.imcombine", PS_LOG_INFO, "norm (1): %f (%f)\n", stamps->normValue, stamps->normValue2);
+    psLogMsg ("psModules.imcombine", PS_LOG_INFO, "norm (1): %f (%f) %ld\n", stamps->normValue, stamps->normValue2,norms->n);
 
     psFree(stats);
@@ -790 +790 @@
         }
     }
-
+    //    psLogMsg("psModules.imcombine",PS_LOG_INFO, "stampNorm: %f %f %d %d %d %f %f",
+    //       stamp->x,stamp->y,footprint,normWindow1,normWindow2,
+    //       normI1,normI2);
     stamp->norm = normI2 / normI1;
     stamp->normI1 = normI1;

trunk/psModules/src/objects/pmPSF.c

@@ -342 +342 @@
     PS_ASSERT_PTR_NON_NULL(modelPar, axes);
 
-    bool useReff = true;
+    bool useReff = false;
     useReff |= (type == pmModelClassGetType ("PS_MODEL_SERSIC"));
     useReff |= (type == pmModelClassGetType ("PS_MODEL_DEV"));
@@ -350 +350 @@
         shape.sx  = modelPar[PM_PAR_SXX];
         shape.sy  = modelPar[PM_PAR_SYY];
-        shape.sxy = modelPar[PM_PAR_SXY];
+        shape.sxy = modelPar[PM_PAR_SXY]; // This is potentially wrong?
     } else {
         shape.sx  = modelPar[PM_PAR_SXX] / M_SQRT2;
@@ -384 +384 @@
     psEllipseShape shape = psEllipseAxesToShape (axes);
 
-    bool useReff = true;
-    useReff |= pmModelClassGetType ("PS_MODEL_SERSIC");
-    useReff |= pmModelClassGetType ("PS_MODEL_DEV");
-    useReff |= pmModelClassGetType ("PS_MODEL_EXP");
+    bool useReff = false;
+    useReff |= ( type == pmModelClassGetType ("PS_MODEL_SERSIC"));
+    useReff |= ( type == pmModelClassGetType ("PS_MODEL_DEV"));
+    useReff |= ( type == pmModelClassGetType ("PS_MODEL_EXP"));
 
     if (useReff) {
         modelPar[PM_PAR_SXX] = shape.sx;
         modelPar[PM_PAR_SYY] = shape.sy;
-        modelPar[PM_PAR_SXY] = shape.sxy;
+        modelPar[PM_PAR_SXY] = shape.sxy; // This is potentially wrong?
     } else {
         modelPar[PM_PAR_SXX] = shape.sx * M_SQRT2;