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

Timestamp:

Nov 13, 2008, 10:58:27 PM (17 years ago)

Author:

Paul Price

Message:

Filter source list used for stamps to remove sources that have nearby neighbours. We don't want to choose as stamps something that isn't a real source, but more like a ton of unresolved sources.

File:

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

trunk/ppStack/src/ppStackMatch.c

@@ -47 +47 @@
 #endif
 
+static psArray *stackSourcesFilter(psArray *sources, // Source list to filter
+                                   int exclusion // Exclusion zone, pixels
+    )
+{
+    psAssert(sources && sources->n > 0, "Require array of sources");
+    if (exclusion <= 0) {
+        return psMemIncrRefCounter(sources);
+    }
+
+    int num = sources->n;               // Number of sources
+    psVector *x = psVectorAlloc(num, PS_TYPE_F32), *y = psVectorAlloc(num, PS_TYPE_F32); // Coordinates
+    int numGood = 0;                    // Number of good sources
+    for (int i = 0; i < num; i++) {
+        pmSource *source = sources->data[i]; // Source of interest
+        if (!source) {
+            continue;
+        }
+        x->data.F32[numGood] = source->peak->xf;
+        y->data.F32[numGood] = source->peak->yf;
+        numGood++;
+    }
+    x->n = y->n = numGood;
+
+    psTree *tree = psTreePlant(2, 2, x, y); // kd tree
+
+    psArray *filtered = psArrayAllocEmpty(numGood); // Filtered list of sources
+    psVector *coords = psVectorAlloc(2, PS_TYPE_F64); // Coordinates of source
+    int numFiltered = 0;                // Number of filtered sources
+    for (int i = 0; i < num; i++) {
+        pmSource *source = sources->data[i]; // Source of interest
+        if (!source) {
+            continue;
+        }
+        coords->data.F32[0] = source->peak->xf;
+        coords->data.F32[1] = source->peak->yf;
+        if (psTreeWithin(tree, coords, exclusion) == 1) {
+            // Only itself inside the exclusion zone
+            filtered = psArrayAdd(filtered, filtered->n, source);
+        } else {
+            numFiltered++;
+        }
+    }
+    psFree(coords);
+    psFree(tree);
+
+    psLogMsg("ppStack", PS_LOG_INFO, "Filtered out %d of %d sources", numFiltered, numGood);
+
+    return filtered;
+}
+
+
 bool ppStackMatch(pmReadout *readout, psArray **regions, psArray **kernels, float *chi2, float *weighting,
-                  const psArray *sources, const pmPSF *psf, psRandom *rng, const pmConfig *config)
+                  psArray *sources, const pmPSF *psf, psRandom *rng, const pmConfig *config)
 {
     assert(readout);
@@ -247 +298 @@
             }
 
+            // For the sake of stamps, remove nearby sources
+            psArray *stampSources = stackSourcesFilter(sources, footprint); // Filtered list of sources
+
             if (threads > 0) {
                 pmSubtractionThreadsInit(readout, fake);
@@ -253 +307 @@
             // Do the image matching
             if (!pmSubtractionMatch(output, NULL, readout, fake, footprint, stride, regionSize, spacing,
-                                    threshold, sources, stampsName, type, size, order, widths, orders, inner,
-                                    ringsOrder, binning, penalty, optimum, optWidths, optOrder, optThresh,
-                                    iter, rej, sysError, maskVal, maskBad, maskPoor, poorFrac, badFrac,
-                                    PM_SUBTRACTION_MODE_1)) {
+                                    threshold, stampSources, stampsName, type, size, order, widths, orders,
+                                    inner, ringsOrder, binning, penalty, optimum, optWidths, optOrder,
+                                    optThresh, iter, rej, sysError, maskVal, maskBad, maskPoor, poorFrac,
+                                    badFrac, PM_SUBTRACTION_MODE_1)) {
                 psError(PS_ERR_UNKNOWN, false, "Unable to match images.");
                 psFree(fake);
                 psFree(optWidths);
+                psFree(stampSources);
                 psFree(output);
                 return false;
@@ -286 +341 @@
             psFree(fake);
             psFree(optWidths);
+            psFree(stampSources);
 
             if (threads > 0) {