Changeset 28484 for branches/pap/ppStack/src/ppStackMatch.c

Timestamp:

Jun 24, 2010, 2:59:09 PM (16 years ago)

Author:

Paul Price

Message:

Merging trunk in advance of reintegrating into trunk.

Location:

branches/pap

Files:

• . (modified) (1 prop)
• ppStack/src/ppStackMatch.c (modified) (10 diffs)

branches/pap/ppStack/src/ppStackMatch.c

@@ -13 +13 @@
 #define MAG_IGNORE 50                   // Ignore magnitudes fainter than this --- they're not real!
 #define FAKE_SIZE 1                     // Size of fake convolution kernel
-#define SOURCE_MASK (PM_SOURCE_MODE_FAIL | PM_SOURCE_MODE_DEFECT | PM_SOURCE_MODE_SATURATED | \
-                     PM_SOURCE_MODE_CR_LIMIT | PM_SOURCE_MODE_EXT_LIMIT) // Mask to apply to input sources
-#define NOISE_FRACTION 0.01             // Set minimum flux to this fraction of noise
 #define COVAR_FRAC 0.01                 // Truncation fraction for covariance matrix
 
@@ -49 +46 @@
 #endif
 
-// Get coordinates from a source
-static void coordsFromSource(float *x, float *y, const pmSource *source)
-{
-    assert(x && y);
-    assert(source);
-
-    if (source->modelPSF) {
-        *x = source->modelPSF->params->data.F32[PM_PAR_XPOS];
-        *y = source->modelPSF->params->data.F32[PM_PAR_YPOS];
-    } else {
-        *x = source->peak->xf;
-        *y = source->peak->yf;
-    }
-    return;
-}
-
-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;
-        }
-        coordsFromSource(&x->data.F32[numGood], &y->data.F32[numGood], source);
-        numGood++;
-    }
-    x->n = y->n = numGood;
-
-    psTree *tree = psTreePlant(2, 2, PS_TREE_EUCLIDEAN, 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;
-        }
-        float xSource, ySource;         // Coordinates of source
-        coordsFromSource(&xSource, &ySource, source);
-
-        coords->data.F64[0] = xSource;
-        coords->data.F64[1] = ySource;
-
-        long numWithin = psTreeWithin(tree, coords, exclusion); // Number within exclusion zone
-        psTrace("ppStack", 9, "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
-            filtered = psArrayAdd(filtered, filtered->n, source);
-        } else {
-            numFiltered++;
-        }
-    }
-    psFree(coords);
-    psFree(tree);
-    psFree(x);
-    psFree(y);
-
-    psLogMsg("ppStack", PS_LOG_INFO, "Filtered out %d of %d sources", numFiltered, numGood);
-
-    return filtered;
-}
 
 // Add background into the fake image
@@ -202 +126 @@
 
 
-bool ppStackMatch(pmReadout *readout, ppStackOptions *options, int index, const pmConfig *config)
+bool ppStackMatch(pmReadout *readout, const psImage *target, ppStackOptions *options,
+                  int index, const pmConfig *config)
 {
     assert(readout);
@@ -286 +211 @@
         } else {
 #endif
-
             // Normal operations here
-            psAssert(options->psf, "Require target PSF");
-            psAssert(options->sourceLists && options->sourceLists->data[index], "Require source list");
+            psAssert(target, "Require target PSF image");
 
             int order = psMetadataLookupS32(NULL, ppsub, "SPATIAL.ORDER"); // Spatial polynomial order
@@ -341 +264 @@
             }
 
-            pmReadout *fake = pmReadoutAlloc(NULL); // Fake readout with target PSF
-
-            psStats *bg = psStatsAlloc(PS_STAT_ROBUST_STDEV); // Statistics for background
-            psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); // Random number generator
-            if (!psImageBackground(bg, NULL, readout->image, readout->mask, maskVal | maskBad, rng)) {
-                psError(PPSTACK_ERR_DATA, false, "Can't measure background for image.");
-                psFree(fake);
-                psFree(optWidths);
-                psFree(conv);
-                psFree(bg);
-                psFree(rng);
-                return false;
-            }
-            float minFlux = NOISE_FRACTION * bg->robustStdev; // Minimum flux level for fake image
-            psFree(rng);
-            psFree(bg);
-
-            // For the sake of stamps, remove nearby sources
-            psArray *stampSources = stackSourcesFilter(options->sourceLists->data[index],
-                                                       footprint); // Filtered list of sources
-
-            bool oldThreads = pmReadoutFakeThreads(true); // Old threading state
-            if (!pmReadoutFakeFromSources(fake, readout->image->numCols, readout->image->numRows,
-                                          stampSources, SOURCE_MASK, NULL, NULL, options->psf,
-                                          minFlux, footprint + size, false, true)) {
-                psError(PPSTACK_ERR_DATA, false, "Unable to generate fake image with target PSF.");
-                psFree(fake);
-                psFree(optWidths);
-                psFree(conv);
-                return false;
-            }
-            pmReadoutFakeThreads(oldThreads);
-
+            pmReadout *fake = pmReadoutAlloc(NULL);
+            fake->image = psImageCopy(NULL, target, PS_TYPE_F32);
             fake->mask = psImageCopy(NULL, readout->mask, PS_TYPE_IMAGE_MASK);
 
@@ -420 +312 @@
                     psFree(fake);
                     psFree(optWidths);
-                    psFree(stampSources);
                     psFree(conv);
                     if (threads > 0) {
@@ -436 +327 @@
                     psFree(fake);
                     psFree(optWidths);
-                    psFree(stampSources);
                     psFree(conv);
                     psFree(widthsCopy);
@@ -446 +336 @@
 
                 if (!pmSubtractionMatch(NULL, conv, fake, readout, footprint, stride, regionSize, spacing,
-                                        threshold, stampSources, stampsName, type, size, order, widthsCopy,
+                                        threshold, options->sources, stampsName, type, size, order, widthsCopy,
                                         orders, inner, ringsOrder, binning, penalty,
                                         optimum, optWidths, optOrder, optThresh, iter, rej, normFrac,
@@ -454 +344 @@
                     psFree(fake);
                     psFree(optWidths);
-                    psFree(stampSources);
                     psFree(conv);
                     psFree(widthsCopy);
@@ -492 +381 @@
             psFree(fake);
             psFree(optWidths);
-            psFree(stampSources);
 
             if (threads > 0) {