DM-50717: Run background subtraction in detectAndMeasure

python/lsst/ip/diffim/detectAndMeasure.py

@@ -85,6 +85,12 @@ class DetectAndMeasureConnections(pipeBase.PipelineTaskConnections,
         storageClass="ExposureF",
         name="{fakesType}{coaddName}Diff_differenceExp",
     )
+    differenceBackground = pipeBase.connectionTypes.Output(
+        doc="Background model that was subtracted from the difference image.",
+        dimensions=("instrument", "visit", "detector"),
+        storageClass="Background",
+        name="difference_background",
+    )
     maskedStreaks = pipeBase.connectionTypes.Output(
         doc='Catalog of streak fit parameters for the difference image.',
         storageClass="ArrowNumpyDict",
@@ -96,6 +102,8 @@ def __init__(self, *, config):
         super().__init__(config=config)
         if not (self.config.writeStreakInfo and self.config.doMaskStreaks):
             self.outputs.remove("maskedStreaks")
+        if not (self.config.doSubtractBackground and self.config.doWriteBackground):
+            self.outputs.remove("differenceBackground")
 
 
 class DetectAndMeasureConfig(pipeBase.PipelineTaskConfig,
@@ -117,6 +125,24 @@ class DetectAndMeasureConfig(pipeBase.PipelineTaskConfig,
         default=False,
         doc="Add columns to the source table to hold analysis metrics?"
     )
+    doSubtractBackground = pexConfig.Field(
+        dtype=bool,
+        doc="Subtract a background model from the image before detection?",
+        default=True,
+    )
+    doWriteBackground = pexConfig.Field(
+        dtype=bool,
+        doc="Persist the fitted background model?",
+        default=False,
+    )
+    subtractInitialBackground = pexConfig.ConfigurableField(
+        target=lsst.meas.algorithms.SubtractBackgroundTask,
+        doc="Task to perform intial background subtraction, before first detection pass.",
+    )
+    subtractFinalBackground = pexConfig.ConfigurableField(
+        target=lsst.meas.algorithms.SubtractBackgroundTask,
+        doc="Task to perform final background subtraction, after first detection pass.",
+    )
     detection = pexConfig.ConfigurableField(
         target=SourceDetectionTask,
         doc="Final source detection for diaSource measurement",
@@ -215,6 +241,20 @@ class DetectAndMeasureConfig(pipeBase.PipelineTaskConfig,
     idGenerator = DetectorVisitIdGeneratorConfig.make_field()
 
     def setDefaults(self):
+        # Background subtraction
+        # Use a small binsize for the first pass to reduce detections on glints
+        #  and extended structures. Should not affect the detectability of
+        #  faint diaSources
+        self.subtractInitialBackground.binSize = 8
+        self.subtractInitialBackground.useApprox = False
+        self.subtractInitialBackground.statisticsProperty = "MEDIAN"
+        self.subtractInitialBackground.doFilterSuperPixels = True
+        # Use a larger binsize for the final background subtraction, to reduce
+        #  over-subtraction of bright objects.
+        self.subtractFinalBackground.binSize = 32
+        self.subtractFinalBackground.useApprox = False
+        self.subtractFinalBackground.statisticsProperty = "MEDIAN"
+        self.subtractFinalBackground.doFilterSuperPixels = True
         # DiaSource Detection
         self.detection.thresholdPolarity = "both"
         self.detection.thresholdValue = 5.0
@@ -226,7 +266,7 @@ def setDefaults(self):
 
         # Copy configs for binned streak detection from the base detection task
         self.streakDetection.thresholdType = self.detection.thresholdType
-        self.streakDetection.reEstimateBackground = self.detection.reEstimateBackground
+        self.streakDetection.reEstimateBackground = False
         self.streakDetection.excludeMaskPlanes = self.detection.excludeMaskPlanes
         self.streakDetection.thresholdValue = self.detection.thresholdValue
         # Only detect positive streaks
@@ -284,6 +324,9 @@ def __init__(self, **kwargs):
         afwTable.CoordKey.addErrorFields(self.schema)
 
         self.algMetadata = dafBase.PropertyList()
+        if self.config.doSubtractBackground:
+            self.makeSubtask("subtractInitialBackground")
+            self.makeSubtask("subtractFinalBackground")
         self.makeSubtask("detection", schema=self.schema)
         if self.config.doDeblend:
             self.makeSubtask("deblend", schema=self.schema)
@@ -381,41 +424,70 @@ def run(self, science, matchedTemplate, difference,
                 Subtracted exposure with detection mask applied.
             ``diaSources``  : `lsst.afw.table.SourceCatalog`
                 The catalog of detected sources.
+            ``differenceBackground`` : `lsst.afw.math.BackgroundList`
+                Background that was subtracted from the difference image.
         """
         if idFactory is None:
             idFactory = lsst.meas.base.IdGenerator().make_table_id_factory()
 
-        self._prepareInputs(difference)
+        if self.config.doSubtractBackground:
+            # Run background subtraction before clearing the mask planes
+            detectionExposure = difference.clone()
+            background = self.subtractInitialBackground.run(detectionExposure).background
+        else:
+            detectionExposure = difference
+            background = afwMath.BackgroundList()
+
+        self._prepareInputs(detectionExposure)
 
         # Don't use the idFactory until after deblend+merge, so that we aren't
         # generating ids that just get thrown away (footprint merge doesn't
         # know about past ids).
         table = afwTable.SourceTable.make(self.schema)
         results = self.detection.run(
             table=table,
-            exposure=difference,
+            exposure=detectionExposure,
             doSmooth=True,
+            background=background
         )
 
+        if self.config.doSubtractBackground:
+            # Run background subtraction again after detecting peaks
+            # but before measurement
+            # First update the mask using the detection image
+            difference.setMask(detectionExposure.mask)
+            background = self.subtractFinalBackground.run(difference).background
+
+            # Re-run detection to get final footprints
+            table = afwTable.SourceTable.make(self.schema)
+            results = self.detection.run(
+                table=table,
+                exposure=difference,
+                doSmooth=True,
+                background=background
+            )
+
         if self.config.doDeblend:
             sources, positives, negatives = self._deblend(difference,
                                                           results.positive,
                                                           results.negative)
 
-            return self.processResults(science, matchedTemplate, difference,
-                                       sources, idFactory,
-                                       positives=positives,
-                                       negatives=negatives)
-
+            result = self.processResults(science, matchedTemplate, difference,
+                                         sources, idFactory,
+                                         positives=positives,
+                                         negatives=negatives)
+            result.differenceBackground = background
         else:
             positives = afwTable.SourceCatalog(self.schema)
             results.positive.makeSources(positives)
             negatives = afwTable.SourceCatalog(self.schema)
             results.negative.makeSources(negatives)
-            return self.processResults(science, matchedTemplate, difference,
-                                       results.sources, idFactory,
-                                       positives=positives,
-                                       negatives=negatives)
+            result = self.processResults(science, matchedTemplate, difference,
+                                         results.sources, idFactory,
+                                         positives=positives,
+                                         negatives=negatives)
+            result.differenceBackground = background
+        return result
 
     def _prepareInputs(self, difference):
         """Ensure that we start with an empty detection and deblended mask.

python/lsst/ip/diffim/subtractImages.py

@@ -175,7 +175,7 @@ class AlardLuptonSubtractBaseConfig(lsst.pex.config.Config):
     doSubtractBackground = lsst.pex.config.Field(
         doc="Subtract the background fit when solving the kernel?",
         dtype=bool,
-        default=True,
+        default=False,
     )
     doApplyExternalCalibrations = lsst.pex.config.Field(
         doc=(

tests/test_detectAndMeasure.py

@@ -99,7 +99,7 @@ def _check_values(self, values, minValue=None, maxValue=None):
         if maxValue is not None:
             self.assertTrue(np.all(values <= maxValue))
 
-    def _setup_detection(self, doSkySources=False, nSkySources=5, **kwargs):
+    def _setup_detection(self, doSkySources=False, nSkySources=5, doSubtractBackground=False, **kwargs):
         """Setup and configure the detection and measurement PipelineTask.
 
         Parameters
@@ -134,6 +134,7 @@ def _setup_detection(self, doSkySources=False, nSkySources=5, **kwargs):
         config.idGenerator.packer["observation"].n_detectors = 99
         config.idGenerator.n_releases = 8
         config.idGenerator.release_id = 2
+        config.doSubtractBackground = doSubtractBackground
         self.idGenerator = config.idGenerator.apply(dataId)
 
         return self.detectionTask(config=config)
@@ -335,6 +336,55 @@ def _detection_wrapper(positive=True):
         _detection_wrapper(positive=True)
         _detection_wrapper(positive=False)
 
+    def test_detect_transients_with_background(self):
+        """Run detection on a difference image containing transients and a background.
+        """
+        # Set up the simulated images
+        noiseLevel = 1.
+        staticSeed = 1
+        transientSeed = 6
+        fluxLevel = 500
+        xSize = 512
+        ySize = 512
+        x0 = 123
+        y0 = 456
+        kwargs = {"seed": staticSeed, "psfSize": 2.4, "fluxLevel": fluxLevel,
+                  "xSize": xSize, "ySize": ySize, "x0": x0, "y0": y0}
+        params = [2.2, 2.1, 2.0, 1.2, 1.1, 1.0]
+
+        bbox2D = lsst.geom.Box2D(lsst.geom.Point2D(x0, y0), lsst.geom.Extent2D(xSize, ySize))
+        background_model = afwMath.Chebyshev1Function2D(params, bbox2D)
+        science, sources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=6,
+                                         background=background_model, **kwargs)
+        matchedTemplate, _ = makeTestImage(noiseLevel=noiseLevel/4, noiseSeed=7, **kwargs)
+
+        # Configure the detection Task
+        detectionTask = self._setup_detection(doMerge=False, doSubtractBackground=True)
+        kwargs["seed"] = transientSeed
+        kwargs["nSrc"] = 10
+        kwargs["fluxLevel"] = 1000
+
+        # Run detection and check the results
+        def _detection_wrapper(positive=True):
+            transients, transientSources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=8, **kwargs)
+            difference = science.clone()
+            difference.maskedImage -= matchedTemplate.maskedImage
+            if positive:
+                difference.maskedImage += transients.maskedImage
+            else:
+                difference.maskedImage -= transients.maskedImage
+            # NOTE: NoiseReplacer (run by forcedMeasurement) can modify the
+            # science image if we've e.g. removed parents post-deblending.
+            # Pass a clone of the science image, so that it doesn't disrupt
+            # later tests.
+            output = detectionTask.run(science.clone(), matchedTemplate, difference)
+            refIds = []
+            scale = 1. if positive else -1.
+            for transient in transientSources:
+                self._check_diaSource(output.diaSources, transient, refIds=refIds, scale=scale)
+        _detection_wrapper(positive=True)
+        _detection_wrapper(positive=False)
+
     def test_missing_mask_planes(self):
         """Check that detection runs with missing mask planes.
         """

tests/test_subtractTask.py

@@ -528,6 +528,11 @@ def test_order_equal_images(self):
     def test_background_subtraction(self):
         """Check that we can recover the background,
         and that it is subtracted correctly in the difference image.
+
+        NOTE: Background subtraction is now turned off by default in
+        subtractImages. It is now run in detectAndMeasure instead, but since the
+        code to run background subtraction is not being removed this test should
+        stay to make sure it continues functioning as intended.
         """
         noiseLevel = 1.
         xSize = 512
@@ -971,7 +976,7 @@ def test_metadata_metrics(self):
 
         # The mean ratio metric should be much worse on the "bad" subtraction.
         self.assertLess(subtractTask_good.metadata['differenceFootprintRatioMean'], 0.02)
-        self.assertGreater(subtractTask_bad.metadata['differenceFootprintRatioMean'], 0.17)
+        self.assertGreater(subtractTask_bad.metadata['differenceFootprintRatioMean'], 0.12)
 
 
 class AlardLuptonPreconvolveSubtractTest(AlardLuptonSubtractTestBase, lsst.utils.tests.TestCase):