Changeset 25406 for branches/eam_branches/20090715/ppStack/src/ppStackMatch.c

Timestamp:

Sep 15, 2009, 4:02:42 PM (17 years ago)

Author:

eugene

Message:

updates from trunk

File:

• branches/eam_branches/20090715/ppStack/src/ppStackMatch.c (modified) (8 diffs)

branches/eam_branches/20090715/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
+                     PM_SOURCE_MODE_CR_LIMIT | PM_SOURCE_MODE_EXT_LIMIT) // Mask to apply to input sources
 #define FAINT_SOURCE_FRAC 1.0e-4         // Set minimum flux to this fraction of faintest source flux
 #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)
@@ -222 +256 @@
             }
             psFitsClose(fits);
-
-            // Add in variance factor
-            pmSubtractionKernels *kernels = psMetadataLookupPtr(NULL, conv->analysis,
-                                                                PM_SUBTRACTION_ANALYSIS_KERNEL); // Kernels
-            float vf = pmSubtractionVarianceFactor(kernels, 0.0, 0.0, false); // Variance factor
-            psMetadataItem *vfItem = psMetadataLookup(readout->parent->concepts, "CELL.VARFACTOR");
-            if (!isfinite(vf)) {
-                vf = 1.0;
-            }
-            if (isfinite(vfItem->data.F32)) {
-                vfItem->data.F32 *= vf;
-            } else {
-                vfItem->data.F32 = vf;
-            }
 
             if (!readImage(&readout->image, options->imageNames->data[index], config) ||
@@ -302 +322 @@
 
             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,
+                                          NAN, footprint + size, false, true)) {
                 psError(PS_ERR_UNKNOWN, false, "Unable to generate fake image with target PSF.");
                 psFree(fake);
@@ -413 +433 @@
             }
 
-            // Set the variance factor
-            psMetadataItem *vfItem = psMetadataLookup(readout->parent->concepts, "CELL.VARFACTOR");
-            float vf = psMetadataLookupF32(NULL, conv->analysis, PM_SUBTRACTION_ANALYSIS_VARFACTOR_1);
-            if (!isfinite(vf)) {
-                vf = 1.0;
-            }
-            if (isfinite(vfItem->data.F32)) {
-                vfItem->data.F32 *= vf;
-            } else {
-                vfItem->data.F32 = vf;
-            }
-
             // Replace original images with convolved
             psFree(readout->image);
@@ -512 +520 @@
         psBinaryOp(readout->variance, readout->variance, "*", psScalarAlloc(PS_SQR(norm), PS_TYPE_F32));
     }
-   
+
     // Ensure the background value is zero
     psStats *bg = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); // Statistics for background
@@ -521 +529 @@
     } else {
       if (!psMetadataLookupBool(NULL, config->arguments, "PPSTACK.SKIP.BG.SUB")) {
-        psLogMsg("ppStack", PS_LOG_INFO, "Correcting convolved image background by %lf (+/- %lf)",
-                 psStatsGetValue(bg, PS_STAT_ROBUST_MEDIAN), psStatsGetValue(bg, PS_STAT_ROBUST_STDEV));
-        (void)psBinaryOp(readout->image, readout->image, "-",
-                         psScalarAlloc(psStatsGetValue(bg, PS_STAT_ROBUST_MEDIAN), PS_TYPE_F32));
+        psLogMsg("ppStack", PS_LOG_INFO, "Correcting convolved image background by %lf (+/- %lf)",
+                 psStatsGetValue(bg, PS_STAT_ROBUST_MEDIAN), psStatsGetValue(bg, PS_STAT_ROBUST_STDEV));
+        (void)psBinaryOp(readout->image, readout->image, "-",
+                         psScalarAlloc(psStatsGetValue(bg, PS_STAT_ROBUST_MEDIAN), PS_TYPE_F32));
       }
     }
-    
+
+    if (!stackRenormaliseReadout(config, readout)) {
+        psFree(rng);
+        psFree(bg);
+        return false;
+    }
 
     // Measure the variance level for the weighting