From 93e2e56562cdeb80d4b455bca074614a4e85532f Mon Sep 17 00:00:00 2001 From: Ian Sullivan Date: Fri, 15 Aug 2025 15:31:09 -0700 Subject: [PATCH 1/3] Move fallback kernel source detection to new method Also apply the source_selector to the kernel candidates detected in the fallback method --- python/lsst/ip/diffim/makeKernel.py | 3 +- python/lsst/ip/diffim/subtractImages.py | 45 ++++++++++++++++++------- 2 files changed, 35 insertions(+), 13 deletions(-) diff --git a/python/lsst/ip/diffim/makeKernel.py b/python/lsst/ip/diffim/makeKernel.py index fee5b1183..b20cd0efd 100644 --- a/python/lsst/ip/diffim/makeKernel.py +++ b/python/lsst/ip/diffim/makeKernel.py @@ -80,7 +80,8 @@ def setDefaults(self): # Minimal set of measurments for star selection self.selectMeasurement.algorithms.names.clear() self.selectMeasurement.algorithms.names = ('base_SdssCentroid', 'base_PsfFlux', 'base_PixelFlags', - 'base_SdssShape', 'base_GaussianFlux', 'base_SkyCoord') + 'base_SdssShape', 'base_GaussianFlux', 'base_SkyCoord', + 'base_ClassificationSizeExtendedness') self.selectMeasurement.slots.modelFlux = None self.selectMeasurement.slots.apFlux = None self.selectMeasurement.slots.calibFlux = None diff --git a/python/lsst/ip/diffim/subtractImages.py b/python/lsst/ip/diffim/subtractImages.py index 8813d7a37..1bfacbb8a 100644 --- a/python/lsst/ip/diffim/subtractImages.py +++ b/python/lsst/ip/diffim/subtractImages.py @@ -541,17 +541,7 @@ def runMakeKernel(self, template, science, sources, convolveTemplate=True): if self.config.allowKernelSourceDetection and convolveTemplate: self.log.warning("Error encountered trying to construct the matching kernel" f" Running source detection and retrying. {e}") - kernelSize = self.makeKernel.makeKernelBasisList( - referenceFwhmPix, targetFwhmPix)[0].getWidth() - sigmaToFwhm = 2*np.log(2*np.sqrt(2)) - candidateList = self.makeKernel.makeCandidateList(reference, target, kernelSize, - candidateList=None, - sigma=targetFwhmPix/sigmaToFwhm) - kernelSources = self.makeKernel.selectKernelSources(reference, target, - candidateList=candidateList, - preconvolved=False, - templateFwhmPix=referenceFwhmPix, - scienceFwhmPix=targetFwhmPix) + kernelSources = self.runKernelSourceDetection(template, science) kernelResult = self.makeKernel.run(reference, target, kernelSources, preconvolved=False, templateFwhmPix=referenceFwhmPix, @@ -571,6 +561,38 @@ def runMakeKernel(self, template, science, sources, convolveTemplate=True): psfMatchingKernel=kernelResult.psfMatchingKernel, kernelSources=kernelSources) + def runKernelSourceDetection(self, template, science): + """Run detection on the science image and use the template mask plane + to reject candidate sources. + + Parameters + ---------- + template : `lsst.afw.image.ExposureF` + Template exposure, warped to match the science exposure. + science : `lsst.afw.image.ExposureF` + Science exposure to subtract from the template. + + Returns + ------- + kernelSources : `lsst.afw.table.SourceCatalog` + Sources from the input catalog to use to construct the + PSF-matching kernel. + """ + kernelSize = self.makeKernel.makeKernelBasisList( + self.templatePsfSize, self.sciencePsfSize)[0].getWidth() + sigmaToFwhm = 2*np.log(2*np.sqrt(2)) + candidateList = self.makeKernel.makeCandidateList(template, science, kernelSize, + candidateList=None, + sigma=self.sciencePsfSize/sigmaToFwhm) + sources = self.makeKernel.selectKernelSources(template, science, + candidateList=candidateList, + preconvolved=False, + templateFwhmPix=self.templatePsfSize, + scienceFwhmPix=self.sciencePsfSize) + + # return sources + return self._sourceSelector(sources, science.getBBox()) + def runConvolveTemplate(self, template, science, psfMatchingKernel, backgroundModel=None): """Convolve the template image with a PSF-matching kernel and subtract from the science image. @@ -875,7 +897,6 @@ def _sourceSelector(self, sources, bbox): If there are too few sources to compute the PSF matching kernel remaining after source selection. """ - selected = self.sourceSelector.selectSources(sources).selected if self.config.restrictKernelEdgeSources: rejectRadius = 2*self.config.makeKernel.kernel.active.kernelSize From 4eaf81a58b7f7234d35b942f6627c8cffe194eec Mon Sep 17 00:00:00 2001 From: Ian Sullivan Date: Fri, 15 Aug 2025 15:32:06 -0700 Subject: [PATCH 2/3] Add simplified image subtraction task --- python/lsst/ip/diffim/subtractImages.py | 172 ++++++++++++++++++++++++ tests/test_subtractTask.py | 59 +++++++- 2 files changed, 228 insertions(+), 3 deletions(-) diff --git a/python/lsst/ip/diffim/subtractImages.py b/python/lsst/ip/diffim/subtractImages.py index 1bfacbb8a..589c72dcf 100644 --- a/python/lsst/ip/diffim/subtractImages.py +++ b/python/lsst/ip/diffim/subtractImages.py @@ -37,6 +37,7 @@ __all__ = ["AlardLuptonSubtractConfig", "AlardLuptonSubtractTask", "AlardLuptonPreconvolveSubtractConfig", "AlardLuptonPreconvolveSubtractTask", + "SimplifiedSubtractConfig", "SimplifiedSubtractTask", "InsufficientKernelSourcesError"] _dimensions = ("instrument", "visit", "detector") @@ -154,6 +155,24 @@ class AlardLuptonSubtractConnections(SubtractInputConnections, SubtractImageOutp pass +class SimplifiedSubtractConnections(SubtractInputConnections, SubtractImageOutputConnections): + inputPsfMatchingKernel = connectionTypes.Input( + doc="Kernel used to PSF match the science and template images.", + dimensions=("instrument", "visit", "detector"), + storageClass="MatchingKernel", + name="{fakesType}{coaddName}Diff_psfMatchKernel", + ) + + def __init__(self, *, config=None): + super().__init__(config=config) + del self.sources + if config.useExistingKernel: + del self.psfMatchingKernel + del self.kernelSources + else: + del self.inputPsfMatchingKernel + + class AlardLuptonSubtractBaseConfig(lsst.pex.config.Config): makeKernel = lsst.pex.config.ConfigurableField( target=MakeKernelTask, @@ -1377,6 +1396,159 @@ def _shapeTest(exp1, exp2, fwhmExposureBuffer, fwhmExposureGrid): return xTest | yTest +class SimplifiedSubtractConfig(AlardLuptonSubtractBaseConfig, lsst.pipe.base.PipelineTaskConfig, + pipelineConnections=SimplifiedSubtractConnections): + mode = lsst.pex.config.ChoiceField( + dtype=str, + default="convolveTemplate", + allowed={"auto": "Choose which image to convolve at runtime.", + "convolveScience": "Only convolve the science image.", + "convolveTemplate": "Only convolve the template image."}, + doc="Choose which image to convolve at runtime, or require that a specific image is convolved." + ) + useExistingKernel = lsst.pex.config.Field( + dtype=bool, + default=True, + doc="Use a pre-existing PSF matching kernel?" + "If False, source detection and measurement will be run." + ) + + +class SimplifiedSubtractTask(AlardLuptonSubtractTask): + """Compute the image difference of a science and template image using + the Alard & Lupton (1998) algorithm. + """ + ConfigClass = SimplifiedSubtractConfig + _DefaultName = "simplifiedSubtract" + + @timeMethod + def run(self, template, science, visitSummary=None, inputPsfMatchingKernel=None): + """PSF match, subtract, and decorrelate two images. + + Parameters + ---------- + template : `lsst.afw.image.ExposureF` + Template exposure, warped to match the science exposure. + science : `lsst.afw.image.ExposureF` + Science exposure to subtract from the template. + visitSummary : `lsst.afw.table.ExposureCatalog`, optional + Exposure catalog with external calibrations to be applied. Catalog + uses the detector id for the catalog id, sorted on id for fast + lookup. + + Returns + ------- + results : `lsst.pipe.base.Struct` + ``difference`` : `lsst.afw.image.ExposureF` + Result of subtracting template and science. + ``matchedTemplate`` : `lsst.afw.image.ExposureF` + Warped and PSF-matched template exposure. + ``backgroundModel`` : `lsst.afw.math.Function2D` + Background model that was fit while solving for the + PSF-matching kernel + ``psfMatchingKernel`` : `lsst.afw.math.Kernel` + Kernel used to PSF-match the convolved image. + ``kernelSources` : `lsst.afw.table.SourceCatalog` + Sources detected on the science image that were used to + construct the PSF-matching kernel. + + Raises + ------ + RuntimeError + If an unsupported convolution mode is supplied. + RuntimeError + If there are too few sources to calculate the PSF matching kernel. + lsst.pipe.base.NoWorkFound + Raised if fraction of good pixels, defined as not having NO_DATA + set, is less then the configured requiredTemplateFraction + """ + self._prepareInputs(template, science, visitSummary=visitSummary) + + convolveTemplate = self.chooseConvolutionMethod(template, science) + + if self.config.useExistingKernel: + psfMatchingKernel = inputPsfMatchingKernel + backgroundModel = None + kernelSources = None + else: + kernelResult = self.runMakeKernel(template, science, convolveTemplate=convolveTemplate) + psfMatchingKernel = kernelResult.psfMatchingKernel + kernelSources = kernelResult.kernelSources + if self.config.doSubtractBackground: + backgroundModel = kernelResult.backgroundModel + else: + backgroundModel = None + if convolveTemplate: + subtractResults = self.runConvolveTemplate(template, science, psfMatchingKernel, + backgroundModel=backgroundModel) + else: + subtractResults = self.runConvolveScience(template, science, psfMatchingKernel, + backgroundModel=backgroundModel) + if kernelSources is not None: + subtractResults.kernelSources = kernelSources + return subtractResults + + def runMakeKernel(self, template, science, convolveTemplate=True): + """Construct the PSF-matching kernel. + + Parameters + ---------- + template : `lsst.afw.image.ExposureF` + Template exposure, warped to match the science exposure. + science : `lsst.afw.image.ExposureF` + Science exposure to subtract from the template. + sources : `lsst.afw.table.SourceCatalog` + Identified sources on the science exposure. This catalog is used to + select sources in order to perform the AL PSF matching on stamp + images around them. + convolveTemplate : `bool`, optional + Construct the matching kernel to convolve the template? + + Returns + ------- + results : `lsst.pipe.base.Struct` + ``backgroundModel`` : `lsst.afw.math.Function2D` + Background model that was fit while solving for the + PSF-matching kernel + ``psfMatchingKernel`` : `lsst.afw.math.Kernel` + Kernel used to PSF-match the convolved image. + ``kernelSources` : `lsst.afw.table.SourceCatalog` + Sources from the input catalog that were used to construct the + PSF-matching kernel. + """ + if convolveTemplate: + reference = template + target = science + referenceFwhmPix = self.templatePsfSize + targetFwhmPix = self.sciencePsfSize + else: + reference = science + target = template + referenceFwhmPix = self.sciencePsfSize + targetFwhmPix = self.templatePsfSize + try: + # The try..except block catches any error, and raises + # NoWorkFound if the template coverage is insufficient. Otherwise, + # the original error is raised. + kernelSources = self.runKernelSourceDetection(template, science) + kernelResult = self.makeKernel.run(reference, target, kernelSources, + preconvolved=False, + templateFwhmPix=referenceFwhmPix, + scienceFwhmPix=targetFwhmPix) + except (RuntimeError, lsst.pex.exceptions.Exception) as e: + self.log.warning("Failed to match template. Checking coverage") + # Raise NoWorkFound if template fraction is insufficient + checkTemplateIsSufficient(template[science.getBBox()], science, self.log, + self.config.minTemplateFractionForExpectedSuccess, + exceptionMessage="Template coverage lower than expected to succeed." + f" Failure is tolerable: {e}") + # checkTemplateIsSufficient did not raise NoWorkFound, so raise original exception + raise e + return lsst.pipe.base.Struct(backgroundModel=kernelResult.backgroundModel, + psfMatchingKernel=kernelResult.psfMatchingKernel, + kernelSources=kernelSources) + + def _interpolateImage(maskedImage, badMaskPlanes, fallbackValue=None): """Replace masked image pixels with interpolated values. diff --git a/tests/test_subtractTask.py b/tests/test_subtractTask.py index 907e2c50c..6be9d7cf7 100644 --- a/tests/test_subtractTask.py +++ b/tests/test_subtractTask.py @@ -40,11 +40,15 @@ class AlardLuptonSubtractTestBase: - def _setup_subtraction(self, **kwargs): + def _setup_subtraction(self, fluxField="truth_instFlux", errField="truth_instFluxErr", **kwargs): """Setup and configure the image subtraction PipelineTask. Parameters ---------- + fluxField : `str`, optional + Name of the flux field in the source catalog. + errField : `str`, optional + Name of the flux error field in the source catalog. **kwargs Any additional config parameters to set. @@ -55,8 +59,8 @@ def _setup_subtraction(self, **kwargs): """ config = self.subtractTask.ConfigClass() config.doSubtractBackground = False - config.sourceSelector.signalToNoise.fluxField = "truth_instFlux" - config.sourceSelector.signalToNoise.errField = "truth_instFluxErr" + config.sourceSelector.signalToNoise.fluxField = fluxField + config.sourceSelector.signalToNoise.errField = errField config.sourceSelector.doUnresolved = True config.sourceSelector.doIsolated = True config.sourceSelector.doRequirePrimary = True @@ -1331,6 +1335,55 @@ def _run_and_check_images(doDecorrelation): _run_and_check_images(doDecorrelation=False) +class SimplifiedSubtractTest(AlardLuptonSubtractTestBase, lsst.utils.tests.TestCase): + subtractTask = subtractImages.SimplifiedSubtractTask + + def test_runSimplifiedTaskWithExistingKernel(self): + """Test that the simplified task produces the same output as + `AlardLuptonSubtractTask` if it uses the AL kernel. + """ + noiseLevel = 1. + science, sources = makeTestImage(psfSize=3.0, noiseLevel=noiseLevel, noiseSeed=6) + template, _ = makeTestImage(psfSize=2.0, noiseLevel=noiseLevel, noiseSeed=7, + templateBorderSize=20, doApplyCalibration=True) + alTask = AlardLuptonSubtractTest._setup_subtraction(AlardLuptonSubtractTest()) + task = self._setup_subtraction(useExistingKernel=True) + + alResults = alTask.run(template.clone(), science.clone(), sources) + results = task.run(template.clone(), science.clone(), + inputPsfMatchingKernel=alResults.psfMatchingKernel) + + self.assertMaskedImagesEqual(alResults.difference, results.difference) + + def test_runSimplifiedTaskWithSourceDetection(self): + """Test that the simplified task with source detection produces + reasonable output. + """ + noiseLevel = 1. + science, sources = makeTestImage(psfSize=3.0, noiseLevel=noiseLevel, noiseSeed=6) + template, _ = makeTestImage(psfSize=2.0, noiseLevel=noiseLevel, noiseSeed=7, + templateBorderSize=20, doApplyCalibration=True) + task = self._setup_subtraction(useExistingKernel=False, + fluxField="base_PsfFlux_instFlux", + errField="base_PsfFlux_instFluxErr", + ) + + output = task.run(template, science) + + # There shoud be no NaN values in the difference image + self.assertTrue(np.all(np.isfinite(output.difference.image.array))) + # Mean of difference image should be close to zero. + meanError = noiseLevel/np.sqrt(output.difference.image.array.size) + # Make sure to include pixels with the DETECTED mask bit set. + statsCtrl = makeStats(badMaskPlanes=("EDGE", "BAD", "NO_DATA", "DETECTED", "DETECTED_NEGATIVE")) + differenceMean = computeRobustStatistics(output.difference.image, output.difference.mask, statsCtrl) + self.assertFloatsAlmostEqual(differenceMean, 0, atol=5*meanError) + # stddev of difference image should be close to expected value. + differenceStd = computeRobustStatistics(output.difference.image, output.difference.mask, + makeStats(), statistic=afwMath.STDEV) + self.assertFloatsAlmostEqual(differenceStd, np.sqrt(2)*noiseLevel, rtol=0.1) + + def setup_module(module): lsst.utils.tests.init() From 830e0b4dc9199164fe5bd92ab6a77696ccbe6ad2 Mon Sep 17 00:00:00 2001 From: Ian Sullivan Date: Sat, 23 Aug 2025 20:19:21 -0700 Subject: [PATCH 3/3] Use fallback source selector if running kernel source detection --- python/lsst/ip/diffim/subtractImages.py | 26 ++++++++++++++++++++++--- 1 file changed, 23 insertions(+), 3 deletions(-) diff --git a/python/lsst/ip/diffim/subtractImages.py b/python/lsst/ip/diffim/subtractImages.py index 589c72dcf..b0fa9b85c 100644 --- a/python/lsst/ip/diffim/subtractImages.py +++ b/python/lsst/ip/diffim/subtractImages.py @@ -228,6 +228,14 @@ class AlardLuptonSubtractBaseConfig(lsst.pex.config.Config): target=ScienceSourceSelectorTask, doc="Task to select sources to be used for PSF matching.", ) + fallbackSourceSelector = lsst.pex.config.ConfigurableField( + target=ScienceSourceSelectorTask, + doc="Task to select sources to be used for PSF matching." + "Used only if the kernel calculation fails and" + "`allowKernelSourceDetection` is set. The fallback source detection" + " will not include all of the same plugins as the original source " + " detection, so not all of the same flags can be used.", + ) detectionThreshold = lsst.pex.config.Field( dtype=float, default=10, @@ -301,6 +309,14 @@ def setDefaults(self): self.sourceSelector.doSignalToNoise = True # apply signal to noise filter self.sourceSelector.signalToNoise.minimum = 10 self.sourceSelector.signalToNoise.maximum = 500 + self.fallbackSourceSelector.doSkySources = False # Do not include sky sources + self.fallbackSourceSelector.doSignalToNoise = True # apply signal to noise filter + self.fallbackSourceSelector.signalToNoise.minimum = 10 + # The following two configs should not be necessary to be turned on for + # PSF-matching, and the fallback kernel source selection will fail if + # they are set since it does not run deblending. + self.fallbackSourceSelector.doIsolated = False # Do not apply isolated star selection + self.fallbackSourceSelector.doRequirePrimary = False # Do not apply primary flag selection class AlardLuptonSubtractConfig(AlardLuptonSubtractBaseConfig, lsst.pipe.base.PipelineTaskConfig, @@ -327,6 +343,7 @@ def __init__(self, **kwargs): self.makeSubtask("decorrelate") self.makeSubtask("makeKernel") self.makeSubtask("sourceSelector") + self.makeSubtask("fallbackSourceSelector") if self.config.doScaleVariance: self.makeSubtask("scaleVariance") @@ -610,7 +627,7 @@ def runKernelSourceDetection(self, template, science): scienceFwhmPix=self.sciencePsfSize) # return sources - return self._sourceSelector(sources, science.getBBox()) + return self._sourceSelector(sources, science.getBBox(), fallback=True) def runConvolveTemplate(self, template, science, psfMatchingKernel, backgroundModel=None): """Convolve the template image with a PSF-matching kernel and subtract @@ -891,7 +908,7 @@ def _convolveExposure(self, exposure, kernel, convolutionControl, else: return convolvedExposure[bbox] - def _sourceSelector(self, sources, bbox): + def _sourceSelector(self, sources, bbox, fallback=False): """Select sources from a catalog that meet the selection criteria. The selection criteria include any configured parameters of the `sourceSelector` subtask, as well as distance from the edge if @@ -916,7 +933,10 @@ def _sourceSelector(self, sources, bbox): If there are too few sources to compute the PSF matching kernel remaining after source selection. """ - selected = self.sourceSelector.selectSources(sources).selected + if fallback: + selected = self.fallbackSourceSelector.selectSources(sources).selected + else: + selected = self.sourceSelector.selectSources(sources).selected if self.config.restrictKernelEdgeSources: rejectRadius = 2*self.config.makeKernel.kernel.active.kernelSize bbox.grow(-rejectRadius)