From 4f4dc42790fdccb8bc26fb441dc2dba4fd8ff43b Mon Sep 17 00:00:00 2001 From: John Parejko Date: Thu, 8 May 2025 15:11:19 -0700 Subject: [PATCH] Move background subtraction to detectAndMeasure --- python/lsst/ip/diffim/detectAndMeasure.py | 96 ++++++++++++++++++++--- python/lsst/ip/diffim/subtractImages.py | 2 +- tests/test_detectAndMeasure.py | 52 +++++++++++- tests/test_subtractTask.py | 7 +- 4 files changed, 142 insertions(+), 15 deletions(-) diff --git a/python/lsst/ip/diffim/detectAndMeasure.py b/python/lsst/ip/diffim/detectAndMeasure.py index c7d59150f..9d059a4b5 100644 --- a/python/lsst/ip/diffim/detectAndMeasure.py +++ b/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,11 +424,21 @@ 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 @@ -393,29 +446,48 @@ def run(self, science, matchedTemplate, difference, 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. diff --git a/python/lsst/ip/diffim/subtractImages.py b/python/lsst/ip/diffim/subtractImages.py index f1ab71e54..2e8a8acff 100644 --- a/python/lsst/ip/diffim/subtractImages.py +++ b/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=( diff --git a/tests/test_detectAndMeasure.py b/tests/test_detectAndMeasure.py index 81c82e470..86b8c315d 100644 --- a/tests/test_detectAndMeasure.py +++ b/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. """ diff --git a/tests/test_subtractTask.py b/tests/test_subtractTask.py index ca6e1c47e..845c318fc 100644 --- a/tests/test_subtractTask.py +++ b/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):