Changeset 34353 for trunk/psphot

Timestamp:

Aug 24, 2012, 12:37:46 PM (14 years ago)

Author:

bills

Message:

Change psphotEfficiency to not fault if one of the readouts fails. Print message
clear the error and move on to next readout (if any).
When calculating the normilization and sourceRadius sample a number of points to avoid
failures if the input happens to not have pixels near the center of the image.

File:

• trunk/psphot/src/psphotEfficiency.c (modified) (5 diffs)

trunk/psphot/src/psphotEfficiency.c

@@ -57 +57 @@
     psFree(stats);
 
+#ifndef USE_SINGLE_POINT_FOR_MODEL
     // Need to normalise out difference between Gaussian and real PSF
+    int sizeX = ro->variance->numCols / 16;
+    int sizeY = ro->variance->numRows / 16;
+    int numPoints = 0;
+    float sum =  0;
+    for (int y = sizeY * 0.5 ; y < ro->variance->numRows; y += sizeY) {
+        for (int x = sizeX * 0.5 ; x < ro->variance->numCols; x += sizeX) {
+            pmModel *normModel = pmModelFromPSFforXY(psf, (float) x, (float) y, 1.0); // model for normalization
+
+            if (!normModel || (normModel->flags & MODEL_MASK)) {
+                psFree(normModel);
+                continue;
+            }
+            float flux = normModel->modelFlux(normModel->params); // Total flux for peak of 1.0
+            psFree(normModel);
+            if (!isfinite(flux)) {
+                continue;
+            }
+            numPoints++;
+            sum += flux;
+        }
+    }
+    if (!isfinite(sum) || numPoints == 0) {
+        psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to generate PSF model for any of %d points.", sizeX * sizeY);
+        return false;
+    }
+    *norm = sum / numPoints;
+#else
+    // This can fail for readout that has few good pixels. It's better to sample many points.
     pmModel *normModel = pmModelFromPSFforXY(psf, 0.5 * ro->variance->numCols,
-                                             0.5 * ro->variance->numRows, 1.0); // Model for normalisation
+                                         0.5 * ro->variance->numRows, 1.0); // Model for normalisation
+    psFree(normModel);
     if (!normModel || (normModel->flags & MODEL_MASK)) {
         psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to generate PSF model.");
@@ -67 +97 @@
     *norm = normModel->modelFlux(normModel->params); // Total flux for peak of 1.0
     psFree(normModel);
+#endif
 
     // The signal-to-noise of the smoothed image is: S/N ~ I_smooth / sqrt(variance * factor)
@@ -188 +219 @@
         if (!psphotEfficiencyReadout (config, view, filerule, i, recipe)) {
             psError (PSPHOT_ERR_CONFIG, false, "failed to measure detection efficiency for %s entry %d", filerule, i);
-            return false;
+            psErrorStackPrint(stderr, " ");
+            psErrorClear();
+            // send message to log as well
+            psLogMsg ("psphot", PS_LOG_WARN, "failed to measure detection efficiency for %s entry %d", filerule, i);
         }
     }
@@ -376 +410 @@
         float mag = magLim + magOffsets->data.F32[i]; // Magnitude for bin
         float peak = powf(10.0, -0.4 * mag) / norm;   // Peak flux
+#ifndef USE_SINGLE_POINT_FOR_MODEL
+        int sizeX = numCols / 16;
+        int sizeY = numRows / 16;
+        int numPoints = 0;
+        float sum =  0;
+        for (int y = sizeY * 0.5 ; y < numRows; y += sizeY) {
+            for (int x = sizeX * 0.5 ; x < numCols; x += sizeX) {
+                // Need to normalise out difference between Gaussian and real PSF
+                pmModel *model = pmModelFromPSFforXY(psf, (float) x, (float) y, peak); 
+
+                if (!model || (model->flags & MODEL_MASK)) {
+                    psFree(model);
+                    continue;
+                }
+                float sourceRadius = PS_MAX(radius, model->modelRadius(model->params, minFlux)); // Radius for source
+                psFree(model);
+                if (!isfinite(sourceRadius)) {
+                    continue;
+                }
+                numPoints++;
+                sum += sourceRadius;
+            }
+        }
+        if (!isfinite(sum) || numPoints == 0) {
+            psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to generate PSF model for any of %d points.", sizeX * sizeY);
+            return false;
+        }
+        float sourceRadius = sum / numPoints;
+#else
+        // This old way can fail for heavily masked images
         pmModel *model = pmModelFromPSFforXY(psf, numCols / 2.0, numRows / 2.0, peak); // Model for source
         if (!model || (model->flags & MODEL_MASK)) {
@@ -384 +448 @@
         float sourceRadius = PS_MAX(radius, model->modelRadius(model->params, minFlux)); // Radius for source
         psFree(model);
+#endif
 
         psArray *sources = psArrayAllocEmpty(numSources); // Sources in this bin