DM-53153: Add safeguards on image statistics based scalings in detectCoaddPeaks

python/lsst/meas/algorithms/dynamicDetection.py

@@ -7,7 +7,7 @@
 
 import numpy as np
 
-from lsst.pex.config import Field, ConfigurableField
+from lsst.pex.config import Field, ConfigurableField, FieldValidationError
 
 from .detection import SourceDetectionConfig, SourceDetectionTask
 from .skyObjects import SkyObjectsTask
@@ -16,7 +16,7 @@
 from lsst.afw.geom import makeCdMatrix, makeSkyWcs, SpanSet
 from lsst.afw.table import SourceCatalog, SourceTable
 from lsst.meas.base import ForcedMeasurementTask
-from lsst.pipe.base import Struct
+from lsst.pipe.base import NoWorkFound, Struct
 
 import lsst.afw.image
 import lsst.afw.math
@@ -67,8 +67,25 @@ class DynamicDetectionConfig(SourceDetectionConfig):
                                 doc="Multiplier for the negative (relative to positive) polarity "
                                 "detections threshold to use for first pass (to find sky objects).")
     skyObjects = ConfigurableField(target=SkyObjectsTask, doc="Generate sky objects.")
+    minGoodPixelFraction = Field(dtype=float, default=0.005,
+                                 doc="Minimum fraction of 'good' pixels required to be deemed "
+                                 "worthwhile for an attempt at further processing.")
+    doThresholdScaling = Field(dtype=bool, default=True,
+                               doc="Scale the threshold level to get empirically measured S/N requested?")
+    minThresholdScaleFactor = Field(dtype=float, default=0.1, optional=True,
+                                    doc="Mininum threshold scaling allowed (i.e. it will be set to this "
+                                    "if the computed value is smaller than it). Set to None for no limit.")
+    maxThresholdScaleFactor = Field(dtype=float, default=10.0, optional=True,
+                                    doc="Maximum threshold scaling allowed (i.e. it will be set to this "
+                                    "if the computed value is greater than it). Set to None for no limit.")
     doBackgroundTweak = Field(dtype=bool, default=True,
                               doc="Tweak background level so median PSF flux of sky objects is zero?")
+    minBackgroundTweak = Field(dtype=float, default=-100.0, optional=True,
+                               doc="Mininum background tweak allowed (i.e. it will be set to this "
+                               "if the computed value is smaller than it). Set to None for no limit.")
+    maxBackgroundTweak = Field(dtype=float, default=100.0, optional=True,
+                               doc="Maximum background tweak allowed (i.e. it will be set to this "
+                               "if the computed value is greater than it). Set to None for no limit.")
     minFractionSources = Field(dtype=float, default=0.02,
                                doc="Minimum fraction of the requested number of sky sources for dynamic "
                                "detection to be considered a success. If the number of good sky sources "
@@ -126,6 +143,22 @@ def validate(self):
         if self.doApplyFlatBackgroundRatio:
             raise ValueError("DynamicDetectionTask does not support doApplyFlatBackgroundRatio.")
 
+        if self.doThresholdScaling:
+            if self.minThresholdScaleFactor and self.maxThresholdScaleFactor:
+                if self.minThresholdScaleFactor > self.maxThresholdScaleFactor:
+                    msg = "minThresholdScaleFactor must be <= maxThresholdScaleFactor"
+                    raise FieldValidationError(
+                        DynamicDetectionConfig.doThresholdScaling, self, msg,
+                    )
+
+        if self.doBackgroundTweak:
+            if self.minBackgroundTweak and self.maxBackgroundTweak:
+                if self.minBackgroundTweak > self.maxBackgroundTweak:
+                    msg = "minBackgroundTweak must be <= maxBackgroundTweak"
+                    raise FieldValidationError(
+                        DynamicDetectionConfig.doBackgroundTweak, self, msg,
+                    )
+
 
 class DynamicDetectionTask(SourceDetectionTask):
     """Detection of sources on an image with a dynamic threshold
@@ -399,8 +432,12 @@ def detectFootprints(self, exposure, doSmooth=True, sigma=None, clearMask=True,
         nPix = maskedImage.mask.array.size
         badPixelMask = lsst.afw.image.Mask.getPlaneBitMask(["NO_DATA", "BAD"])
         nGoodPix = np.sum(maskedImage.mask.array & badPixelMask == 0)
-        self.log.info("Number of good data pixels (i.e. not NO_DATA or BAD): {} ({:.1f}% of total)".
+        self.log.info("Number of good data pixels (i.e. not NO_DATA or BAD): {} ({:.2f}% of total)".
                       format(nGoodPix, 100*nGoodPix/nPix))
+        if nGoodPix/nPix < self.config.minGoodPixelFraction:
+            msg = (f"Image has a very low good pixel fraction ({nGoodPix} of {nPix}), so not worth further "
+                   "consideration")
+            raise NoWorkFound(msg)
 
         with self.tempWideBackgroundContext(exposure):
             # Could potentially smooth with a wider kernel than the PSF in
@@ -409,36 +446,52 @@ def detectFootprints(self, exposure, doSmooth=True, sigma=None, clearMask=True,
             psf = self.getPsf(exposure, sigma=sigma)
             convolveResults = self.convolveImage(maskedImage, psf, doSmooth=doSmooth)
 
-            if self.config.doBrightPrelimDetection:
-                brightDetectedMask = self._computeBrightDetectionMask(maskedImage, convolveResults)
-
-            middle = convolveResults.middle
-            sigma = convolveResults.sigma
-            prelim = self.applyThreshold(
-                middle, maskedImage.getBBox(), factor=self.config.prelimThresholdFactor,
-                factorNeg=self.config.prelimNegMultiplier*self.config.prelimThresholdFactor
-            )
-            self.finalizeFootprints(
-                maskedImage.mask, prelim, sigma, factor=self.config.prelimThresholdFactor,
-                factorNeg=self.config.prelimNegMultiplier*self.config.prelimThresholdFactor
-            )
-            if self.config.doBrightPrelimDetection:
-                # Combine prelim and bright detection masks for multiplier
-                # determination.
-                maskedImage.mask.array |= brightDetectedMask
+            if self.config.doThresholdScaling:
+                if self.config.doBrightPrelimDetection:
+                    brightDetectedMask = self._computeBrightDetectionMask(maskedImage, convolveResults)
+            else:
+                prelim = None
+                factor = 1.0
 
-            # Calculate the proper threshold
             # seed needs to fit in a C++ 'int' so pybind doesn't choke on it
             seed = (expId if expId is not None else int(maskedImage.image.array.sum())) % (2**31 - 1)
-            threshResults = self.calculateThreshold(exposure, seed, sigma=sigma)
-            minMultiplicative = 0.5
-            if threshResults.multiplicative < minMultiplicative:
-                self.log.warning("threshResults.multiplicative = %.2f is less than minimum value (%.2f). "
-                                 "Setting to %.2f.", threshResults.multiplicative, minMultiplicative,
-                                 minMultiplicative)
-            factor = max(minMultiplicative, threshResults.multiplicative)
-            self.log.info("Modifying configured detection threshold by factor %.2f to %.2f",
-                          factor, factor*self.config.thresholdValue)
+
+            middle = convolveResults.middle
+            sigma = convolveResults.sigma
+            if self.config.doThresholdScaling:
+                prelim = self.applyThreshold(
+                    middle, maskedImage.getBBox(), factor=self.config.prelimThresholdFactor,
+                    factorNeg=self.config.prelimNegMultiplier*self.config.prelimThresholdFactor
+                )
+                self.finalizeFootprints(
+                    maskedImage.mask, prelim, sigma, factor=self.config.prelimThresholdFactor,
+                    factorNeg=self.config.prelimNegMultiplier*self.config.prelimThresholdFactor
+                )
+                if self.config.doBrightPrelimDetection:
+                    # Combine prelim and bright detection masks for multiplier
+                    # determination.
+                    maskedImage.mask.array |= brightDetectedMask
+
+                # Calculate the proper threshold
+                threshResults = self.calculateThreshold(exposure, seed, sigma=sigma)
+                if (self.config.minThresholdScaleFactor
+                        and threshResults.multiplicative < self.config.minThresholdScaleFactor):
+                    self.log.warning("Measured threshold scaling factor (%.2f) is outside [min, max] "
+                                     "bounds [%.2f, %.2f].  Setting factor to lower limit: %.2f.",
+                                     threshResults.multiplicative, self.config.minThresholdScaleFactor,
+                                     self.config.maxThresholdScaleFactor, self.config.minThresholdScaleFactor)
+                    factor = self.config.minThresholdScaleFactor
+                elif (self.config.maxThresholdScaleFactor
+                        and threshResults.multiplicative > self.config.maxThresholdScaleFactor):
+                    self.log.warning("Measured threshold scaling factor (%.2f) is outside [min, max] "
+                                     "bounds [%.2f, %.2f].  Setting factor to upper limit: %.2f.",
+                                     threshResults.multiplicative, self.config.minThresholdScaleFactor,
+                                     self.config.maxThresholdScaleFactor, self.config.maxThresholdScaleFactor)
+                    factor = self.config.maxThresholdScaleFactor
+                else:
+                    factor = threshResults.multiplicative
+                self.log.info("Modifying configured detection threshold by factor %.2f to %.2f",
+                              factor, factor*self.config.thresholdValue)
 
             growOverride = None
             inFinalize = True
@@ -516,6 +569,18 @@ def detectFootprints(self, exposure, doSmooth=True, sigma=None, clearMask=True,
                 self.finalizeFootprints(maskedImage.mask, tweakDetResults, sigma, factor=factor)
                 bgLevel = self.calculateThreshold(exposure, seed, sigma=sigma, minFractionSourcesFactor=0.5,
                                                   isBgTweak=True).additive
+                if self.config.minBackgroundTweak and bgLevel < self.config.minBackgroundTweak:
+                    self.log.warning("Measured background tweak (%.2f) is outside [min, max] bounds "
+                                     "[%.2f, %.2f].  Setting tweak to lower limit: %.2f.", bgLevel,
+                                     self.config.minBackgroundTweak, self.config.maxBackgroundTweak,
+                                     self.config.minBackgroundTweak)
+                    bgLevel = self.config.minBackgroundTweak
+                if self.config.maxBackgroundTweak and bgLevel > self.config.maxBackgroundTweak:
+                    self.log.warning("Measured background tweak (%.2f) is outside [min, max] bounds "
+                                     "[%.2f, %.2f].  Setting tweak to upper limit: %.2f.", bgLevel,
+                                     self.config.minBackgroundTweak, self.config.maxBackgroundTweak,
+                                     self.config.maxBackgroundTweak)
+                    bgLevel = self.config.maxBackgroundTweak
             finally:
                 maskedImage.mask.array[:] = originalMask
         else:

tests/test_dynamicDetection.py

@@ -7,6 +7,7 @@
 from lsst.geom import Box2I, Extent2I, Point2D, Point2I, SpherePoint, degrees
 from lsst.meas.algorithms import DynamicDetectionTask
 from lsst.meas.algorithms.testUtils import plantSources
+from lsst.pex.config import FieldValidationError
 
 
 class DynamicDetectionTest(lsst.utils.tests.TestCase):
@@ -42,37 +43,108 @@ def setUp(self):
     def tearDown(self):
         del self.exposure
 
-    def check(self, expectFactor):
+    def check(self, expectFactor, config):
         schema = SourceTable.makeMinimalSchema()
-        task = DynamicDetectionTask(config=self.config, schema=schema)
+        task = DynamicDetectionTask(config=config, schema=schema)
         table = SourceTable.make(schema)
 
         results = task.run(table, self.exposure, expId=12345)
         self.assertFloatsAlmostEqual(results.factor, expectFactor, rtol=self.rtol)
 
     def testVanilla(self):
         """Dynamic detection used as normal detection."""
-        self.check(1.0)
+        self.check(1.0, self.config)
 
     def testDynamic(self):
         """Modify the variance plane, and see if the task is able to determine
         what we did.
         """
         factor = 2.0
         self.exposure.maskedImage.variance /= factor
-        self.check(1.0/np.sqrt(factor))
+        self.check(1.0/np.sqrt(factor), self.config)
 
     def testNoWcs(self):
         """Check that dynamic detection runs when the exposure wcs is None."""
         self.exposure.setWcs(None)
-        self.check(1.0)
+        self.check(1.0, self.config)
 
     def testMinimalSkyObjects(self):
         """Check that dynamic detection runs when there are a relatively small
         number of sky objects.
         """
-        self.config.skyObjects.nSources = int(0.1 * self.config.skyObjects.nSources)
-        self.check(1.0)
+        config = DynamicDetectionTask.ConfigClass()
+        config.skyObjects.nSources = int(0.1 * self.config.skyObjects.nSources)
+        self.check(1.0, config)
+
+    def testDynamicNoLimits(self):
+        """Test setting the threshold scaling and background tweak limits to
+        None (i.e. no limits imposed).
+        """
+        config = DynamicDetectionTask.ConfigClass()
+        config.minThresholdScaleFactor = None
+        config.maxThresholdScaleFactor = None
+        config.minBackgroundTweak = None
+        config.maxBackgroundTweak = None
+        self.check(1.0, config)
+
+    def testNoThresholdScaling(self):
+        """Check that dynamic detection runs when doThresholdScaling is False.
+        """
+        config = DynamicDetectionTask.ConfigClass()
+        config.doThresholdScaling = False
+        self.check(1.0, config)
+
+    def testNoBackgroundTweak(self):
+        """Check that dynamic detection runs when doBackgroundTweak is False.
+        """
+        config = DynamicDetectionTask.ConfigClass()
+        config.doBackgroundTweak = False
+        self.check(1.0, config)
+
+    def testThresholdScalingAndNoBackgroundTweak(self):
+        """Check that dynamic detection runs when both doThresholdScaling and
+        doBackgroundTweak are False.
+        """
+        config = DynamicDetectionTask.ConfigClass()
+        config.doThresholdScaling = False
+        config.doBackgroundTweak = False
+        self.check(1.0, config)
+
+    def testThresholdsOutsideBounds(self):
+        """Check that dynamic detection properly sets threshold limits.
+        """
+        schema = SourceTable.makeMinimalSchema()
+        config = DynamicDetectionTask.ConfigClass()
+        config.minThresholdScaleFactor = 1.05
+        config.maxBackgroundTweak = -0.1
+        table = SourceTable.make(schema)
+        task = DynamicDetectionTask(config=config, schema=schema)
+        task.run(table, self.exposure, expId=12345)
+
+        config = DynamicDetectionTask.ConfigClass()
+        config.maxThresholdScaleFactor = 0.99
+        config.minBackgroundTweak = 0.1
+        table = SourceTable.make(schema)
+        task = DynamicDetectionTask(config=config, schema=schema)
+        task.run(table, self.exposure, expId=12345)
+
+    def testConfigValidation(self):
+        """Check that the field validation is working correctly.
+        """
+        schema = SourceTable.makeMinimalSchema()
+        config = DynamicDetectionTask.ConfigClass()
+        config.minThresholdScaleFactor = 1.05
+        config.maxThresholdScaleFactor = 1.01
+        with self.assertRaisesRegex(
+                FieldValidationError, "minThresholdScaleFactor must be <= maxThresholdScaleFactor"):
+            DynamicDetectionTask(config=config, schema=schema)
+
+        config = DynamicDetectionTask.ConfigClass()
+        config.minBackgroundTweak = 2.0
+        config.maxBackgroundTweak = 1.0
+        with self.assertRaisesRegex(
+                FieldValidationError, "minBackgroundTweak must be <= maxBackgroundTweak"):
+            DynamicDetectionTask(config=config, schema=schema)
 
     def testNoSkyObjects(self):
         """Check that dynamic detection runs when there are no sky objects.