Changeset 25766 for branches/eam_branches/20090820/ppStack/src/ppStackMatch.c

Timestamp:

Oct 2, 2009, 5:10:19 PM (17 years ago)

Author:

eugene

Message:

merge changes from trunk

Location:

branches/eam_branches/20090820

Files:

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

branches/eam_branches/20090820/ppStack/src/ppStackMatch.c

@@ -14 +14 @@
 #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) // Mask to apply to input sources
-#define FAINT_SOURCE_FRAC 1.0e-4         // Set minimum flux to this fraction of faintest source flux
+                     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
 
@@ -87 +87 @@
     x->n = y->n = numGood;
 
-    psTree *tree = psTreePlant(2, 2, x, y); // kd tree
+    psTree *tree = psTreePlant(2, 2, PS_TREE_EUCLIDEAN, x, y); // kd tree
 
     psArray *filtered = psArrayAllocEmpty(numGood); // Filtered list of sources
@@ -162 +162 @@
 }
 
+// Renormalise a readout's variance map
+bool stackRenormaliseReadout(const pmConfig *config, // Configuration
+                             pmReadout *readout      // Readout to renormalise
+    )
+{
+    bool mdok; // Status of metadata lookups
+
+    psMetadata *recipe = psMetadataLookupPtr(NULL, config->recipes, PPSTACK_RECIPE); // Recipe for ppStack
+    psAssert(recipe, "Need PPSTACK recipe");
+
+    if (!psMetadataLookupBool(&mdok, recipe, "RENORM")) return true;
+
+    int num = psMetadataLookupS32(&mdok, recipe, "RENORM.NUM");
+    if (!mdok) {
+        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "RENORM.NUM is not set in the recipe");
+        return false;
+    }
+    float minValid = psMetadataLookupF32(&mdok, recipe, "RENORM.MIN");
+    if (!mdok) {
+        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "RENORM.MIN is not set in the recipe");
+        return false;
+    }
+    float maxValid = psMetadataLookupF32(&mdok, recipe, "RENORM.MAX");
+    if (!mdok) {
+        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "RENORM.MAX is not set in the recipe");
+        return false;
+    }
+
+    psImageMaskType maskBad = pmConfigMaskGet("BLANK", config); // Bits to mask
+
+    return pmReadoutVarianceRenormalise(readout, maskBad, num, minValid, maxValid);
+}
+
+
 
 bool ppStackMatch(pmReadout *readout, ppStackOptions *options, int index, const pmConfig *config)
@@ -232 +266 @@
             psRegion *region = psMetadataLookupPtr(NULL, conv->analysis,
                                                    PM_SUBTRACTION_ANALYSIS_REGION); // Convolution region
-
-            pmSubtractionAnalysis(readout->analysis, kernels, region,
+            pmSubtractionKernels *kernels = psMetadataLookupPtr(NULL, conv->analysis,
+                                                                PM_SUBTRACTION_ANALYSIS_KERNEL);
+
+            pmSubtractionAnalysis(readout->analysis, NULL, kernels, region,
                                   readout->image->numCols, readout->image->numRows);
 
@@ -283 +319 @@
             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(PS_ERR_UNKNOWN, 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],
@@ -288 +339 @@
 
             if (!pmReadoutFakeFromSources(fake, readout->image->numCols, readout->image->numRows,
-                                          stampSources, NULL, NULL, options->psf, NAN, footprint + size,
-                                          false, true)) {
+                                          stampSources, SOURCE_MASK, NULL, NULL, options->psf,
+                                          minFlux, footprint + size, false, true)) {
                 psError(PS_ERR_UNKNOWN, false, "Unable to generate fake image with target PSF.");
                 psFree(fake);
@@ -335 +386 @@
                                                                PM_SUBTRACTION_ANALYSIS_KERNEL); // Conv kernel
             if (kernel) {
-                if (!pmSubtractionMatchPrecalc(conv, NULL, readout, fake, readout->analysis,
+                if (!pmSubtractionMatchPrecalc(NULL, conv, fake, readout, readout->analysis,
                                                stride, sysError, maskVal, maskBad, maskPoor,
                                                poorFrac, badFrac)) {
@@ -349 +400 @@
                 }
             } else {
-                if (!pmSubtractionMatch(conv, NULL, readout, fake, footprint, stride, regionSize, spacing,
+                if (!pmSubtractionMatch(NULL, conv, fake, readout, footprint, stride, regionSize, spacing,
                                         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)) {
+                                        PM_SUBTRACTION_MODE_2)) {
                     psError(PS_ERR_UNKNOWN, false, "Unable to match images.");
                     psFree(fake);
@@ -502 +553 @@
     }
 
+    if (!stackRenormaliseReadout(config, readout)) {
+        psFree(rng);
+        psFree(bg);
+        return false;
+    }
 
     // Measure the variance level for the weighting