Changeset 25027 for branches/pap/psModules/src/imcombine/pmPSFEnvelope.c

Timestamp:

Aug 7, 2009, 4:08:25 PM (17 years ago)

Author:

Paul Price

Message:

Merging trunk (r25026) to get up-to-date on old branch.

Location:

branches/pap

Files:

• . (modified) (1 prop)
• psModules (modified) (1 prop)
• psModules/src/imcombine/pmPSFEnvelope.c (modified) (16 diffs)

branches/pap/psModules/src/imcombine/pmPSFEnvelope.c

@@ -33 +33 @@
 
 
-//#define TESTING                         // Enable test output
+// #define TESTING                         // Enable test output
 #define PEAK_FLUX 1.0e4                 // Peak flux for each source
 #define SKY_VALUE 0.0e0                 // Sky value for fake image
@@ -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
+    float maxRadius = 0.0;              // Maximum radius for sources
+    psVector *numbers = psVectorAlloc(numFakes, PS_TYPE_S32); // Number of detections for each source
+    psVectorInit(numbers, 0);
     for (int i = 0; i < inputs->n; i++) {
         pmPSF *psf = inputs->data[i];   // PSF of interest
@@ -127 +131 @@
             psFree(xOffset);
             psFree(fakes);
+            psFree(numbers);
             psf->residuals = resid;
             return NULL;
@@ -140 +145 @@
 
             double flux = fakeRO->image->data.F32[(int)y][(int)x];
+            if (!isfinite(flux) || flux < 0) {
+                continue;
+            }
             float norm = PEAK_FLUX / flux; // Normalisation for source
             psRegion region = psRegionSet(x - radius, x + radius, y - radius, y + radius); // PSF region
@@ -151 +159 @@
             // 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?");
+            if (!model || (model->flags & MODEL_MASK)) {
+                continue;
+            }
             float srcRadius = model->modelRadius(model->params, PS_SQR(VARIANCE_VAL)); // Radius for source
+            if (srcRadius == 0) {
+                continue;
+            }
             if (srcRadius > maxRadius) {
                 maxRadius = srcRadius;
             }
-            psFree(model);
+
+            // If we got this far, the source is decent
+            numbers->data.S32[j]++;
         }
 
@@ -175 +190 @@
     psFree(fakeRO);
 
-    if (maxRadius > radius) {
-        maxRadius = radius;
-    }
-
 #ifdef TESTING
     {
@@ -190 +201 @@
 
     // Put the fake sources onto a full-size image
+    psArray *goodFakes = psArrayAllocEmpty(numFakes); // Good fake sources
     pmReadout *readout = pmReadoutAlloc(NULL); // Readout to contain envelope pixels
     readout->image = psImageAlloc(numCols, numRows, PS_TYPE_F32);
@@ -195 +207 @@
     for (int i = 0; i < numFakes; i++) {
         pmSource *source = fakes->data[i]; // Fake source
+        if (numbers->data.S32[i] > 0) {
+            psArrayAdd(goodFakes, goodFakes->n, source);
+        }
+
         // Position of source on fake image
         int xFake = source->peak->x + xOffset->data.S32[i];
@@ -213 +229 @@
             psFree(yOffset);
             psFree(fakes);
+            psFree(numbers);
             return NULL;
         }
@@ -220 +237 @@
     psFree(yOffset);
     psFree(envelope);
+    psFree(numbers);
+
+    psFree(fakes);
+    fakes = goodFakes;
+    numFakes = fakes->n;
+
+    if (numFakes == 0.0) {
+        psError(PS_ERR_UNKNOWN, false, "No fake sources are suitable for PSF fitting.");
+        psFree(fakes);
+        psFree(readout);
+        return false;
+    }
 
     // XXX Setting the variance seems to be an art
@@ -232 +261 @@
     psImageInit(readout->mask, 0);
 
+    if (maxRadius > radius) {
+        maxRadius = radius;
+    }
+
 #ifdef TESTING
     {
@@ -243 +276 @@
 
     // 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
@@ -257 +291 @@
         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);
@@ -264 +298 @@
         }
 
-        if (!pmSourceMoments(source, maxRadius)) {
-            psError(PS_ERR_UNKNOWN, false, "Unable to measure moments for source.");
-            psFree(readout);
-            psFree(fakes);
-            return NULL;
-        }
+        // measure the source moments: tophat windowing, no pixel S/N cutoff
+        if (!pmSourceMoments(source, maxRadius, 0.0, 1.0)) {
+            // 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;
     }
 
@@ -286 +328 @@
     options->psfFieldXo = 0;
     options->psfFieldYo = 0;
+    options->chiFluxTrend = false;      // All sources have similar flux, so fitting a trend often fails
 
     pmSourceFitModelInit(SOURCE_FIT_ITERATIONS, 0.01, VARIANCE_VAL, true);