From 2c277e6b466b6d90df27ba3e1846c81d15965829 Mon Sep 17 00:00:00 2001 From: Ian Sullivan Date: Sat, 6 Sep 2025 11:58:58 -0700 Subject: [PATCH 1/6] Add log message with final numbers of diaSources and sky sources --- python/lsst/ip/diffim/detectAndMeasure.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/python/lsst/ip/diffim/detectAndMeasure.py b/python/lsst/ip/diffim/detectAndMeasure.py index 2e895fea2..16b3195bb 100644 --- a/python/lsst/ip/diffim/detectAndMeasure.py +++ b/python/lsst/ip/diffim/detectAndMeasure.py @@ -808,7 +808,10 @@ def processResults(self, science, matchedTemplate, difference, sources, idFactor # This option allows returning sky sources, # even if there are no diaSources measurementResults.diaSources = diaSources - + self.log.info("Measured %d diaSources and %d sky sources", + np.count_nonzero(~diaSources["sky_source"]), + np.count_nonzero(diaSources["sky_source"]) + ) return measurementResults def _deblend(self, difference, positiveFootprints, negativeFootprints): From f96c290bde5f64172d476a6a77efb269d7162605 Mon Sep 17 00:00:00 2001 From: Ian Sullivan Date: Sat, 6 Sep 2025 12:11:34 -0700 Subject: [PATCH 2/6] Add checkMask utility to reject PSF candidates based on template mask --- python/lsst/ip/diffim/subtractImages.py | 16 +++++---- python/lsst/ip/diffim/utils.py | 43 ++++++++++++++++++++++++- tests/test_subtractTask.py | 6 ++-- 3 files changed, 55 insertions(+), 10 deletions(-) diff --git a/python/lsst/ip/diffim/subtractImages.py b/python/lsst/ip/diffim/subtractImages.py index b0fa9b85c..732741b4b 100644 --- a/python/lsst/ip/diffim/subtractImages.py +++ b/python/lsst/ip/diffim/subtractImages.py @@ -26,7 +26,7 @@ import lsst.afw.image import lsst.afw.math import lsst.geom -from lsst.ip.diffim.utils import evaluateMeanPsfFwhm, getPsfFwhm, computeDifferenceImageMetrics +from lsst.ip.diffim.utils import evaluateMeanPsfFwhm, getPsfFwhm, computeDifferenceImageMetrics, checkMask from lsst.meas.algorithms import ScaleVarianceTask, ScienceSourceSelectorTask import lsst.pex.config import lsst.pipe.base @@ -270,7 +270,6 @@ class AlardLuptonSubtractBaseConfig(lsst.pex.config.Config): dtype=str, default=("NO_DATA", "BAD", "SAT", "EDGE", "FAKE"), doc="Mask planes to exclude when selecting sources for PSF matching.", - deprecated="No longer used. Will be removed after v30" ) badMaskPlanes = lsst.pex.config.ListField( dtype=str, @@ -443,7 +442,7 @@ def run(self, template, science, sources, visitSummary=None): convolveTemplate = self.chooseConvolutionMethod(template, science) - selectSources = self._sourceSelector(sources, science.getBBox()) + selectSources = self._sourceSelector(sources, science.getBBox(), template.mask) kernelResult = self.runMakeKernel(template, science, selectSources, convolveTemplate=convolveTemplate) @@ -627,7 +626,7 @@ def runKernelSourceDetection(self, template, science): scienceFwhmPix=self.sciencePsfSize) # return sources - return self._sourceSelector(sources, science.getBBox(), fallback=True) + return self._sourceSelector(sources, science.getBBox(), template.mask, fallback=True) def runConvolveTemplate(self, template, science, psfMatchingKernel, backgroundModel=None): """Convolve the template image with a PSF-matching kernel and subtract @@ -908,7 +907,7 @@ def _convolveExposure(self, exposure, kernel, convolutionControl, else: return convolvedExposure[bbox] - def _sourceSelector(self, sources, bbox, fallback=False): + def _sourceSelector(self, sources, bbox, mask, 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 @@ -920,6 +919,8 @@ def _sourceSelector(self, sources, bbox, fallback=False): Input source catalog to select sources from. bbox : `lsst.geom.Box2I` Bounding box of the science image. + mask : `lsst.afw.image.Mask` + The mask plane of the template to use to reject kernel candidates. Returns ------- @@ -945,6 +946,9 @@ def _sourceSelector(self, sources, bbox, fallback=False): np.count_nonzero(~bboxSelected), rejectRadius) selected &= bboxSelected selectSources = sources[selected].copy(deep=True) + # Optionally remove sources that land on masked pixels + maskSelected = checkMask(mask, selectSources, self.config.excludeMaskPlanes, checkAdjacent=True) + selectSources = selectSources[maskSelected].copy(deep=True) # Trim selectSources if they exceed ``maxKernelSources``. # Keep the highest signal-to-noise sources of those selected. if (len(selectSources) > self.config.maxKernelSources) & (self.config.maxKernelSources > 0): @@ -1211,7 +1215,7 @@ def run(self, template, science, sources, visitSummary=None): interpolateBadMaskPlanes=True) self.metadata["convolvedExposure"] = "Preconvolution" try: - selectSources = self._sourceSelector(sources, science.getBBox()) + selectSources = self._sourceSelector(sources, science.getBBox(), template.mask) subtractResults = self.runPreconvolve(template, science, matchedScience, selectSources, scienceKernel) diff --git a/python/lsst/ip/diffim/utils.py b/python/lsst/ip/diffim/utils.py index fc0ade02e..f43103c4e 100644 --- a/python/lsst/ip/diffim/utils.py +++ b/python/lsst/ip/diffim/utils.py @@ -23,7 +23,7 @@ __all__ = ["evaluateMeanPsfFwhm", "getPsfFwhm", "getKernelCenterDisplacement", - "computeDifferenceImageMetrics", + "computeDifferenceImageMetrics", "checkMask" ] import itertools @@ -454,3 +454,44 @@ def populate_sattle_visit_cache(visit_info, historical=False): "historical": historical}) r.raise_for_status() + + +def checkMask(mask, sources, excludeMaskPlanes, checkAdjacent=True): + """Exclude sources that are located on or adjacent to masked pixels. + + Parameters + ---------- + mask : `lsst.afw.image.Mask` + The image mask plane to use to reject sources + based on the location of their centroid on the ccd. + sources : `lsst.afw.table.SourceCatalog` + The source catalog to evaluate. + excludeMaskPlanes : `list` of `str` + List of the names of the mask planes to exclude. + + Returns + ------- + flags : `numpy.ndarray` of `bool` + Array indicating whether each source in the catalog should be + kept (True) or rejected (False) based on the value of the + mask plane at its location. + """ + setExcludeMaskPlanes = [ + maskPlane for maskPlane in excludeMaskPlanes if maskPlane in mask.getMaskPlaneDict() + ] + + excludePixelMask = mask.getPlaneBitMask(setExcludeMaskPlanes) + + xv = (np.rint(sources.getX() - mask.getX0())).astype(int) + yv = (np.rint(sources.getY() - mask.getY0())).astype(int) + + flags = np.ones(len(sources), dtype=bool) + if checkAdjacent: + pixRange = (0, -1, 1) + else: + pixRange = (0,) + for j in pixRange: + for i in pixRange: + mv = mask.array[yv + j, xv + i] + flags *= np.bitwise_and(mv, excludePixelMask) == 0 + return flags diff --git a/tests/test_subtractTask.py b/tests/test_subtractTask.py index 6be9d7cf7..15f19eb2f 100644 --- a/tests/test_subtractTask.py +++ b/tests/test_subtractTask.py @@ -117,9 +117,9 @@ def test_incomplete_template_coverage(self): border = 20 xSize = 400 ySize = 400 - science, sources = makeTestImage(psfSize=3.0, noiseLevel=noiseLevel, noiseSeed=6, nSrc=50, + science, sources = makeTestImage(psfSize=3.0, noiseLevel=noiseLevel, noiseSeed=6, nSrc=100, xSize=xSize, ySize=ySize) - template, _ = makeTestImage(psfSize=2.0, noiseLevel=noiseLevel, noiseSeed=7, nSrc=50, + template, _ = makeTestImage(psfSize=2.0, noiseLevel=noiseLevel, noiseSeed=7, nSrc=100, templateBorderSize=border, doApplyCalibration=True, xSize=xSize, ySize=ySize) @@ -485,7 +485,7 @@ def _run_and_check_sources(sourcesIn, maxKernelSources=1000, minKernelSources=3, signalToNoise = sources.getPsfInstFlux()/sources.getPsfInstFluxErr() signalToNoise = signalToNoise[~sources[badSourceFlag]] signalToNoise.sort() - selectSources = task._sourceSelector(sources, science.getBBox()) + selectSources = task._sourceSelector(sources, science.getBBox(), template.mask) self.assertEqual(nGoodSources, len(selectSources)) signalToNoiseOut = selectSources.getPsfInstFlux()/selectSources.getPsfInstFluxErr() signalToNoiseOut.sort() From a2edd10d10b5e627fb7e26683c0fe897c455f4c4 Mon Sep 17 00:00:00 2001 From: Ian Sullivan Date: Sat, 6 Sep 2025 12:22:03 -0700 Subject: [PATCH 3/6] Add missing docstrings --- python/lsst/ip/diffim/subtractImages.py | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/python/lsst/ip/diffim/subtractImages.py b/python/lsst/ip/diffim/subtractImages.py index 732741b4b..13ba9cfd1 100644 --- a/python/lsst/ip/diffim/subtractImages.py +++ b/python/lsst/ip/diffim/subtractImages.py @@ -921,6 +921,11 @@ def _sourceSelector(self, sources, bbox, mask, fallback=False): Bounding box of the science image. mask : `lsst.afw.image.Mask` The mask plane of the template to use to reject kernel candidates. + fallback : `bool`, optional + Switch indicating the source selector is being called after + running the fallback source detection subtask, which does not run a + full set of measurement plugins and can't use the same settings for + the source selector. Returns ------- @@ -930,6 +935,9 @@ def _sourceSelector(self, sources, bbox, mask, fallback=False): Raises ------ + InsufficientKernelSourcesError + An AlgorithmError that is raised if there are not enough PSF + candidates to construct the PSF matching kernel. RuntimeError If there are too few sources to compute the PSF matching kernel remaining after source selection. From f46a218b658b4689b5ecc7160750c4d9877808d4 Mon Sep 17 00:00:00 2001 From: Ian Sullivan Date: Sat, 6 Sep 2025 12:23:43 -0700 Subject: [PATCH 4/6] Add a mask plane for high variance regions of the template --- python/lsst/ip/diffim/getTemplate.py | 46 +++++++++++++++++++++++++++- 1 file changed, 45 insertions(+), 1 deletion(-) diff --git a/python/lsst/ip/diffim/getTemplate.py b/python/lsst/ip/diffim/getTemplate.py index 639f67026..b1b8c1a2c 100644 --- a/python/lsst/ip/diffim/getTemplate.py +++ b/python/lsst/ip/diffim/getTemplate.py @@ -34,7 +34,7 @@ from lsst.skymap import BaseSkyMap from lsst.ip.diffim.dcrModel import DcrModel -from lsst.meas.algorithms import CoaddPsf, CoaddPsfConfig +from lsst.meas.algorithms import CoaddPsf, CoaddPsfConfig, SubtractBackgroundTask from lsst.utils.timer import timeMethod __all__ = [ @@ -103,6 +103,17 @@ class GetTemplateConfig( doc="Configuration for CoaddPsf", dtype=CoaddPsfConfig, ) + varianceBackground = pexConfig.ConfigurableField( + target=SubtractBackgroundTask, + doc="Task to estimate the background variance.", + ) + highVarianceThreshold = pexConfig.RangeField( + dtype=float, + default=4, + min=1, + doc="Set the HIGH_VARIANCE mask plane for regions with variance" + " greater than the median by this factor.", + ) def setDefaults(self): # Use a smaller cache: per SeparableKernel.computeCache, this should @@ -115,6 +126,19 @@ def setDefaults(self): self.warp.warpingKernelName = "lanczos3" self.coaddPsf.warpingKernelName = self.warp.warpingKernelName + # Background subtraction of the variance plane + self.varianceBackground.algorithm = "LINEAR" + self.varianceBackground.binSize = 32 + self.varianceBackground.useApprox = False + self.varianceBackground.statisticsProperty = "MEDIAN" + self.varianceBackground.doFilterSuperPixels = True + self.varianceBackground.ignoredPixelMask = ["BAD", + "EDGE", + "DETECTED", + "DETECTED_NEGATIVE", + "NO_DATA", + ] + class GetTemplateTask(pipeBase.PipelineTask): ConfigClass = GetTemplateConfig @@ -137,6 +161,7 @@ def __init__(self, *args, **kwargs): type=float, doc="Weight for each exposure, used to make the CoaddPsf; should always be 1.", ) + self.makeSubtask("varianceBackground") def runQuantum(self, butlerQC, inputRefs, outputRefs): inputs = butlerQC.get(inputRefs) @@ -357,11 +382,30 @@ def run(self, *, coaddExposureHandles, bbox, wcs, dataIds, physical_filter): for c in catalogs: catalog.extend(c) + # Set a mask plane for any regions with exceptionally high variance. + self.checkHighVariance(template) template.setPsf(self._makePsf(template, catalog, wcs)) template.setFilter(afwImage.FilterLabel(band, physical_filter)) template.setPhotoCalib(photoCalib) return pipeBase.Struct(template=template) + def checkHighVariance(self, template): + """Set a mask plane for regions with unusually high variance. + + Parameters + ---------- + template : `lsst.afw.image.Exposure` + The warped template exposure, which will be modified in place. + """ + highVarianceMaskPlaneBit = template.mask.addMaskPlane("HIGH_VARIANCE") + template.mask.getPlaneBitMask("HIGH_VARIANCE") + varianceExposure = template.clone() + varianceExposure.image.array = varianceExposure.variance.array + varianceBackground = self.varianceBackground.run(varianceExposure).background.getImage().array + threshold = self.config.highVarianceThreshold*np.nanmedian(varianceBackground) + highVariancePix = varianceBackground > threshold + template.mask.array[highVariancePix] |= 2**highVarianceMaskPlaneBit + @staticmethod def _checkInputs(dataIds, coaddExposures): """Check that the all the dataIds are from the same band and that From b008426f479b8bcb9088cdd222bcfb797724224d Mon Sep 17 00:00:00 2001 From: Ian Sullivan Date: Sat, 6 Sep 2025 12:24:10 -0700 Subject: [PATCH 5/6] Exclude PSF candidates in high variance regions of the template --- python/lsst/ip/diffim/subtractImages.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/python/lsst/ip/diffim/subtractImages.py b/python/lsst/ip/diffim/subtractImages.py index 13ba9cfd1..e0d998af6 100644 --- a/python/lsst/ip/diffim/subtractImages.py +++ b/python/lsst/ip/diffim/subtractImages.py @@ -268,8 +268,8 @@ class AlardLuptonSubtractBaseConfig(lsst.pex.config.Config): ) excludeMaskPlanes = lsst.pex.config.ListField( dtype=str, - default=("NO_DATA", "BAD", "SAT", "EDGE", "FAKE"), - doc="Mask planes to exclude when selecting sources for PSF matching.", + default=("NO_DATA", "BAD", "SAT", "EDGE", "FAKE", "HIGH_VARIANCE"), + doc="Template mask planes to exclude when selecting sources for PSF matching.", ) badMaskPlanes = lsst.pex.config.ListField( dtype=str, @@ -278,7 +278,7 @@ class AlardLuptonSubtractBaseConfig(lsst.pex.config.Config): ) preserveTemplateMask = lsst.pex.config.ListField( dtype=str, - default=("NO_DATA", "BAD",), + default=("NO_DATA", "BAD", "HIGH_VARIANCE"), doc="Mask planes from the template to propagate to the image difference." ) renameTemplateMask = lsst.pex.config.ListField( From 5475126d9ffef2e36af995d941330365c55eccc5 Mon Sep 17 00:00:00 2001 From: Ian Sullivan Date: Sat, 6 Sep 2025 12:24:41 -0700 Subject: [PATCH 6/6] Add pixelFlags for diaSources measured in high variance regions of the difference image --- python/lsst/ip/diffim/detectAndMeasure.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/python/lsst/ip/diffim/detectAndMeasure.py b/python/lsst/ip/diffim/detectAndMeasure.py index 16b3195bb..0b90d4778 100644 --- a/python/lsst/ip/diffim/detectAndMeasure.py +++ b/python/lsst/ip/diffim/detectAndMeasure.py @@ -416,9 +416,9 @@ def setDefaults(self): # Keep track of which footprints contain streaks self.measurement.plugins["base_PixelFlags"].masksFpAnywhere = [ - "STREAK", "INJECTED", "INJECTED_TEMPLATE"] + "STREAK", "INJECTED", "INJECTED_TEMPLATE", "HIGH_VARIANCE"] self.measurement.plugins["base_PixelFlags"].masksFpCenter = [ - "STREAK", "INJECTED", "INJECTED_TEMPLATE"] + "STREAK", "INJECTED", "INJECTED_TEMPLATE", "HIGH_VARIANCE"] self.skySources.avoidMask = ["DETECTED", "DETECTED_NEGATIVE", "BAD", "NO_DATA", "EDGE"] def validate(self):