Changeset 23960 for trunk/psModules/src/imcombine/pmPSFEnvelope.c

Timestamp:

Apr 23, 2009, 5:48:09 PM (17 years ago)

Author:

Paul Price

Message:

Doing our best to cope when working with skycells on the edge of the FoV --- e.g., the PSF is not necessarily defined in the centre of the image. If the model or moments aren't good, then don't use the source to fit the PSF.

File:

• trunk/psModules/src/imcombine/pmPSFEnvelope.c (modified) (9 diffs)

trunk/psModules/src/imcombine/pmPSFEnvelope.c

@@ -40 +40 @@
 #define PSF_STATS PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV // Statistics options for measuring PSF
 #define SOURCE_FIT_ITERATIONS 100       // Number of iterations for source fitting
+#define MODEL_MASK (PM_MODEL_STATUS_NONCONVERGE | PM_MODEL_STATUS_OFFIMAGE | \
+                    PM_MODEL_STATUS_BADARGS | PM_MODEL_STATUS_LIMITS) // Mask to apply to models
 
 
@@ -112 +114 @@
     psImageInit(envelope, SKY_VALUE);
     pmReadout *fakeRO = pmReadoutAlloc(NULL); // Fake readout
-    float maxRadius = 0.0;              // Maximum radius of sources
     psVector *numbers = psVectorAlloc(numFakes, PS_TYPE_S32); // Number of detections for each source
     psVectorInit(numbers, 0);
@@ -143 +144 @@
 
             double flux = fakeRO->image->data.F32[(int)y][(int)x];
-            if (flux > 0) {
-                // The source is present
-                numbers->data.S32[j]++;
+            if (!isfinite(flux) || flux < 0) {
+                continue;
             }
             float norm = PEAK_FLUX / flux; // Normalisation for source
@@ -156 +156 @@
             psFree(subEnv);
 
-            // Get the radius
-            pmModel *model = pmModelFromPSFforXY(psf, x, y, PEAK_FLUX); // Model for source
-            psAssert (model, "failed to generate model: should this be an error or not?");
-            float srcRadius = model->modelRadius(model->params, PS_SQR(VARIANCE_VAL)); // Radius for source
-            if (srcRadius > maxRadius) {
-                maxRadius = srcRadius;
-            }
-            psFree(model);
+            // If we got this far, the source is decent
+            numbers->data.S32[j]++;
         }
 
@@ -182 +176 @@
     psFree(fakeRO);
 
-    if (maxRadius > radius) {
-        maxRadius = radius;
-    }
-
 #ifdef TESTING
     {
@@ -261 +251 @@
 
     // Reset the sources to point to the new pixels, and measure the moments in preparation for PSF fitting
+    int numMoments = 0;                 // Number of moments measured
     for (int i = 0; i < numFakes; i++) {
         pmSource *source = fakes->data[i]; // Fake source
@@ -275 +266 @@
         source->maskObj = NULL;
 
-        if (!pmSourceDefinePixels(source, readout, x, y, maxRadius)) {
+        if (!pmSourceDefinePixels(source, readout, x, y, radius)) {
             psError(PS_ERR_UNKNOWN, false, "Unable to define pixels for source.");
             psFree(readout);
@@ -282 +273 @@
         }
 
-        if (!pmSourceMoments(source, maxRadius)) {
-            psError(PS_ERR_UNKNOWN, false, "Unable to measure moments for source.");
-            psFree(readout);
-            psFree(fakes);
-            return NULL;
-        }
+        if (!pmSourceMoments(source, radius)) {
+            // Can't do anything about it; limp along as best we can
+            psErrorClear();
+            continue;
+        }
+        numMoments++;
+    }
+
+    if (numMoments == 0) {
+        psError(PS_ERR_UNKNOWN, true, "Unable to measure moments for sources.");
+        psFree(fakes);
+        psFree(readout);
+        return NULL;
     }
 
@@ -296 +294 @@
     options->poissonErrorsParams = true;
     options->stats = psStatsAlloc(PSF_STATS);
-    options->radius = maxRadius;
+    options->radius = radius;
     options->psfTrendMode = PM_TREND_MAP;
     options->psfTrendNx = xOrder;