Index: trunk/psphot/src/psphotEfficiency.c
===================================================================
--- trunk/psphot/src/psphotEfficiency.c	(revision 34352)
+++ trunk/psphot/src/psphotEfficiency.c	(revision 34353)
@@ -57,7 +57,37 @@
     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,4 +97,5 @@
     *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,5 +219,8 @@
 	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,4 +410,34 @@
         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,4 +448,5 @@
         float sourceRadius = PS_MAX(radius, model->modelRadius(model->params, minFlux)); // Radius for source
         psFree(model);
+#endif
 
         psArray *sources = psArrayAllocEmpty(numSources); // Sources in this bin