Changeset 20755 for trunk/ppStack/src/ppStackMatch.c

Timestamp:

Nov 14, 2008, 12:31:34 AM (17 years ago)

Author:

Paul Price

Message:

Filter out sources for stamps before making the fake image. Correct
the fake image with the background around the stamp stars. This seems
to be important, to prevent the convolved image becoming noisy and
around zero.

File:

• trunk/ppStack/src/ppStackMatch.c (modified) (6 diffs)

trunk/ppStack/src/ppStackMatch.c

@@ -80 +80 @@
             continue;
         }
-        coords->data.F32[0] = source->peak->xf;
-        coords->data.F32[1] = source->peak->yf;
+        coords->data.F64[0] = source->peak->xf;
+        coords->data.F64[1] = source->peak->yf;
 
         long numWithin = psTreeWithin(tree, coords, exclusion); // Number within exclusion zone
-        psTrace("ppStack", 5, "Source at %.0f,%.0f has %ld sources in exclusion zone",
-                coords->data.F32[0], coords->data.F32[1], numWithin);
+        psTrace("ppStack", 5, "Source at %.0lf,%.0lf has %ld sources in exclusion zone",
+                coords->data.F64[0], coords->data.F64[1], numWithin);
         if (numWithin == 1) {
             // Only itself inside the exclusion zone
@@ -101 +101 @@
 }
 
+#define BG_SIZE 64                      // Large box half-size in which to measure background
+
+// Generate a background model of the readout we're matching
+psImage *stackBackgroundModel(pmReadout *ro, psMaskType maskVal, const psArray *sources, int size)
+{
+    psImage *image = ro->image, *mask = ro->mask; // Image and mask of readout
+    int numCols = image->numCols, numRows = image->numRows; // Size of image
+
+    psImage *bgImage = psImageAlloc(numCols, numRows, PS_TYPE_F32); // Background image
+    psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN); // Statistics for measuring background
+    psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS, 0); // Random number generator
+
+    for (int i = 0; i < sources->n; i++) {
+        pmSource *source = sources->data[i]; // Source of interest
+        if (!source) {
+            continue;
+        }
+
+        int x = source->peak->xf + 0.5, y = source->peak->yf + 0.5; // Coordinates of source
+
+        int xMin = PS_MAX(x - BG_SIZE, 0), xMax = PS_MIN(x + BG_SIZE, numCols);
+        int yMin = PS_MAX(y - BG_SIZE, 0), yMax = PS_MIN(y + BG_SIZE, numCols);
+
+        psRegion region = psRegionSet(xMin, xMax, yMin, yMax); // Region of interest
+        psImage *subImage = psImageSubset(image, region); // Subimage
+        psImage *subMask = psImageSubset(mask, region); // Subimage mask
+        if (!psImageBackground(stats, NULL, subImage, subMask, maskVal, rng)) {
+            // Can't do anything about it
+            psErrorClear();
+            continue;
+        }
+        psFree(subImage);
+        psFree(subMask);
+
+        float value = stats->robustMedian; // Value to set background to
+
+        int uMin = PS_MAX(x - size, 0), uMax = PS_MIN(x + size, numCols);
+        int vMin = PS_MAX(y - size, 0), vMax = PS_MIN(y + size, numCols);
+
+        for (int v = vMin; v < vMax; v++) {
+            for (int u = uMin; u < uMax; u++) {
+                bgImage->data.F32[v][u] = value;
+            }
+        }
+    }
+
+    psFree(stats);
+    psFree(rng);
+
+    return bgImage;
+}
 
 bool ppStackMatch(pmReadout *readout, psArray **regions, psArray **kernels, float *chi2, float *weighting,
@@ -242 +293 @@
                     psErrorClear();
                 } else {
-                    minFlux = 0.1 * bg->robustStdev;
+                    minFlux = bg->robustStdev;
                 }
                 psFree(bg);
             }
 
+#if 0
             float maxMag = -INFINITY;   // Maximum magnitude of sources
             for (int i = 0; i < sources->n; i++) {
@@ -256 +308 @@
             }
             minFlux = PS_MIN(FAINT_SOURCE_FRAC * powf(10.0, -0.4 * maxMag), minFlux);
+#endif
 
             pmReadout *fake = pmReadoutAlloc(NULL); // Fake readout with target PSF
 
-            if (!pmReadoutFakeFromSources(fake, readout->image->numCols, readout->image->numRows, sources,
-                                          NULL, NULL, psf, minFlux, 0, false)) {
+            // For the sake of stamps, remove nearby sources
+            psArray *stampSources = stackSourcesFilter(sources, footprint); // Filtered list of sources
+
+            if (!pmReadoutFakeFromSources(fake, readout->image->numCols, readout->image->numRows,
+                                          stampSources, NULL, NULL, psf, minFlux, 0, false)) {
                 psError(PS_ERR_UNKNOWN, false, "Unable to generate fake image with target PSF.");
                 psFree(fake);
@@ -267 +323 @@
                 return false;
             }
+
+            // Add the background into the target image
+            psImage *bgImage = stackBackgroundModel(readout, maskVal, stampSources, footprint); // Image of background
+            psBinaryOp(fake->image, fake->image, "+", bgImage);
+            psFree(bgImage);
+
 
 #ifdef TESTING
@@ -301 +363 @@
                 }
             }
-
-            // For the sake of stamps, remove nearby sources
-            psArray *stampSources = stackSourcesFilter(sources, footprint); // Filtered list of sources
 
             if (threads > 0) {