From e3b02dad9bee88aae29a5874689425b28bfdaa22 Mon Sep 17 00:00:00 2001 From: Ian Sullivan Date: Thu, 17 Sep 2026 11:06:16 -0700 Subject: [PATCH 01/17] Use full detector not just bbox in getTemplate component load --- .../lsst.ip.diffim.SimplifiedSubtractTask.rst | 4 ++-- python/lsst/ip/diffim/getTemplate.py | 16 ++++++++++++++-- 2 files changed, 16 insertions(+), 4 deletions(-) diff --git a/doc/lsst.ip.diffim/tasks/lsst.ip.diffim.SimplifiedSubtractTask.rst b/doc/lsst.ip.diffim/tasks/lsst.ip.diffim.SimplifiedSubtractTask.rst index 8ad68e8d3..dd0812546 100644 --- a/doc/lsst.ip.diffim/tasks/lsst.ip.diffim.SimplifiedSubtractTask.rst +++ b/doc/lsst.ip.diffim/tasks/lsst.ip.diffim.SimplifiedSubtractTask.rst @@ -36,7 +36,7 @@ Use a local version of the AP pipeline configuration to run the pipeline by reus rewarpTemplate: class: lsst.ip.diffim.getTemplate.GetTemplateTask config: - connections.bbox: preliminary_visit_image.bbox + connections.detector: preliminary_visit_image.detector connections.wcs: preliminary_visit_image.wcs connections.coaddExposures: template_coadd connections.template: template_detector @@ -77,7 +77,7 @@ Alternatively, the pipeline can compute the PSF matching kernel internally by pe rewarpTemplate: class: lsst.ip.diffim.getTemplate.GetTemplateTask config: - connections.bbox: preliminary_visit_image.bbox + connections.detector: preliminary_visit_image.detector connections.wcs: preliminary_visit_image.wcs connections.coaddExposures: template_coadd connections.template: template_detector diff --git a/python/lsst/ip/diffim/getTemplate.py b/python/lsst/ip/diffim/getTemplate.py index 74e7a7588..a7f3c66e3 100644 --- a/python/lsst/ip/diffim/getTemplate.py +++ b/python/lsst/ip/diffim/getTemplate.py @@ -63,6 +63,14 @@ class GetTemplateConnections( name="{fakesType}calexp.bbox", storageClass="Box2I", dimensions=("instrument", "visit", "detector"), + deprecated="Replaced by the `detector` connection. Will be removed after v31.", + ) + detector = pipeBase.connectionTypes.Input( + doc="Detector of the exposure that we will construct the template for." + " Its bounding box sets the geometry of the output template.", + name="{fakesType}calexp.detector", + storageClass="Detector", + dimensions=("instrument", "visit", "detector"), ) wcs = pipeBase.connectionTypes.Input( doc="WCS of the exposure that we will construct the template for.", @@ -96,6 +104,8 @@ class GetTemplateConnections( def __init__(self, *, config=None): super().__init__(config=config) + # Kept only so that older config files still load. + del self.bbox if config.requireCoaddAtGraphBuild: self.coaddExposures = dataclasses.replace( self.coaddExposures, @@ -225,7 +235,8 @@ def __init__(self, *args, **kwargs): def runQuantum(self, butlerQC, inputRefs, outputRefs): inputs = butlerQC.get(inputRefs) - bbox = inputs.pop("bbox") + detector = inputs.pop("detector") + bbox = detector.getBBox() wcs = inputs.pop("wcs") coaddExposures = inputs.pop("coaddExposures") skymap = inputs.pop("skyMap") @@ -895,7 +906,8 @@ class GetDcrTemplateTask(GetTemplateTask): def runQuantum(self, butlerQC, inputRefs, outputRefs): inputs = butlerQC.get(inputRefs) - bbox = inputs.pop("bbox") + detector = inputs.pop("detector") + bbox = detector.getBBox() wcs = inputs.pop("wcs") dcrCoaddExposureHandles = inputs.pop("dcrCoadds") skymap = inputs.pop("skyMap") From 78673ee0f3b77c3d432bfd97819c5624a3890a66 Mon Sep 17 00:00:00 2001 From: Ian Sullivan Date: Fri, 25 Sep 2026 14:28:29 -0700 Subject: [PATCH 02/17] Make sure renamed mask planes are not cleared --- python/lsst/ip/diffim/subtractImages.py | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/python/lsst/ip/diffim/subtractImages.py b/python/lsst/ip/diffim/subtractImages.py index 4c39a924b..25baf102a 100644 --- a/python/lsst/ip/diffim/subtractImages.py +++ b/python/lsst/ip/diffim/subtractImages.py @@ -1137,10 +1137,15 @@ def updateMasks(self, template, science): # the `renameTemplateMask` config will be copied to new planes with # "_TEMPLATE" appended to their names, and the original mask plane will # be cleared. - clearMaskPlanes = [mp for mp in template.mask.getMaskPlaneDict().keys() - if mp not in self.config.preserveTemplateMask] renameMaskPlanes = [mp for mp in self.config.renameTemplateMask if mp in template.mask.getMaskPlaneDict().keys()] + # afw's mask planes are global, so if the "*_TEMPLATE" mask planes were + # defined anywhere previously they would now also be cleared unless we + # exclude them. + renamedMaskPlanes = {maskPlane + "_TEMPLATE" for maskPlane in renameMaskPlanes} + clearMaskPlanes = [mp for mp in template.mask.getMaskPlaneDict().keys() + if mp not in self.config.preserveTemplateMask + and mp not in renamedMaskPlanes] # propagate the mask plane related to Fake source injection # NOTE: the fake source injection sets FAKE plane, but it should be INJECTED From d4a95c7d95308ee851a92dfa5c81c91feaeb57b3 Mon Sep 17 00:00:00 2001 From: Ian Sullivan Date: Fri, 25 Sep 2026 17:55:37 -0700 Subject: [PATCH 03/17] Support missing mask planes for spatially sampled metrics --- .../diffim/computeSpatiallySampledMetrics.py | 25 +++-- tests/test_computeSpatiallySampledMetrics.py | 104 ++++++++++++++++++ 2 files changed, 120 insertions(+), 9 deletions(-) create mode 100644 tests/test_computeSpatiallySampledMetrics.py diff --git a/python/lsst/ip/diffim/computeSpatiallySampledMetrics.py b/python/lsst/ip/diffim/computeSpatiallySampledMetrics.py index dd8949f6e..2f698e514 100644 --- a/python/lsst/ip/diffim/computeSpatiallySampledMetrics.py +++ b/python/lsst/ip/diffim/computeSpatiallySampledMetrics.py @@ -29,7 +29,7 @@ import lsst.pipe.base as pipeBase import lsst.pex.config as pexConfig -from lsst.ip.diffim.utils import getPsfFwhm, angleMean, evaluateMaskFraction, getKernelCenterDisplacement +from lsst.ip.diffim.utils import getPsfFwhm, angleMean, getKernelCenterDisplacement from lsst.meas.algorithms import SkyObjectsTask from lsst.pex.exceptions import InvalidParameterError, RangeError from lsst.utils.timer import timeMethod @@ -88,7 +88,8 @@ class SpatiallySampledMetricsConfig(pipeBase.PipelineTaskConfig, """ metricsMaskPlanes = lsst.pex.config.ListField( dtype=str, - doc="List of mask planes to include in metrics", + doc="List of mask planes to include in metrics. " + "A mask plane that is not registered reports 0 pixels set.", default=('BAD', 'CLIPPED', 'CR', 'DETECTED', 'DETECTED_NEGATIVE', 'EDGE', 'INEXACT_PSF', 'INJECTED', 'INJECTED_TEMPLATE', 'INTRP', 'NOT_DEBLENDED', 'NO_DATA', 'REJECTED', 'SAT', 'SAT_TEMPLATE', 'SENSOR_EDGE', 'STREAK', 'SUSPECT', @@ -247,12 +248,16 @@ def run(self, science, template, difference, diaSources, psfMatchingKernel): self.metricSources.run(mask=science.mask, seed=difference.info.id, catalog=spatiallySampledMetrics) - metricsMaskPlanes = [] + # Mask plane bits are looked up in the global afw dict, not the + # image's own. A plane that is not registered in this process has no + # pixels set in any image. + metricsMaskPlanes = {} for maskPlane in self.config.metricsMaskPlanes: try: - metricsMaskPlanes.append(maskPlane) + metricsMaskPlanes[maskPlane] = afwImage.Mask.getPlaneBitMask(maskPlane) except InvalidParameterError: - self.log.info("Unable to calculate metrics for mask plane %s: not in image"%maskPlane) + self.log.info("Mask plane %s is not registered; its metrics will be zero.", maskPlane) + metricsMaskPlanes[maskPlane] = 0 for src in spatiallySampledMetrics: self._evaluateLocalMetric(src, science, template, difference, diaSources, @@ -275,8 +280,9 @@ def _evaluateLocalMetric(self, src, science, template, difference, diaSources, The science image. difference : `lsst.afw.image.Exposure` Result of subtracting template from the science image. - metricsMaskPlanes : `list` of `str` - Mask planes to calculate metrics from. + metricsMaskPlanes : `dict` [`str`, `int`] + Mask planes to calculate metrics from, mapped to their bitmask. + The bitmask is 0 for a plane that is not registered. psfMatchingKernel : `~lsst.afw.math.LinearCombinationKernel` The PSF matching kernel of the subtraction to evaluate. """ @@ -324,9 +330,10 @@ def _evaluateLocalMetric(self, src, science, template, difference, diaSources, src.set('diffim_value', diffimVal) src.set('diffim_variance', diffimVar) src.set('diffim_chi2PerPix', self._diffimChi2PerPix(difference[bbox])) - for maskPlane in metricsMaskPlanes: + differenceMask = difference.mask[bbox].array + for maskPlane, bitMask in metricsMaskPlanes.items(): src.set("%s_mask_fraction"%maskPlane.lower(), - evaluateMaskFraction(difference.mask[bbox], maskPlane) + np.count_nonzero(differenceMask & bitMask)/differenceMask.size ) krnlSum, dx, dy, direction, length = getKernelCenterDisplacement( diff --git a/tests/test_computeSpatiallySampledMetrics.py b/tests/test_computeSpatiallySampledMetrics.py new file mode 100644 index 000000000..4f12d1d13 --- /dev/null +++ b/tests/test_computeSpatiallySampledMetrics.py @@ -0,0 +1,104 @@ +# This file is part of ip_diffim. +# +# Developed for the LSST Data Management System. +# This product includes software developed by the LSST Project +# (https://www.lsst.org). +# See the COPYRIGHT file at the top-level directory of this distribution +# for details of code ownership. +# +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . + +import unittest + +import numpy as np + +import lsst.afw.image +import lsst.afw.table +import lsst.ip.diffim +import lsst.utils.tests + +from utils import makeTestImage + + +class ComputeSpatiallySampledMetricsTest(lsst.utils.tests.TestCase): + """Test the metrics that are sampled across the difference image.""" + + def setUp(self): + noiseLevel = 1. + science, sources = makeTestImage(psfSize=2.4, noiseLevel=noiseLevel, noiseSeed=6) + template, _ = makeTestImage(psfSize=2.0, noiseLevel=noiseLevel, noiseSeed=7, + templateBorderSize=20, doApplyCalibration=True) + config = lsst.ip.diffim.AlardLuptonSubtractTask.ConfigClass() + config.doSubtractBackground = False + config.sourceSelector.signalToNoise.fluxField = "truth_instFlux" + config.sourceSelector.signalToNoise.errField = "truth_instFluxErr" + subtraction = lsst.ip.diffim.AlardLuptonSubtractTask(config=config).run(template, science, sources) + self.science = subtraction.matchedScience + self.template = template + self.difference = subtraction.difference + self.psfMatchingKernel = subtraction.psfMatchingKernel + self.diaSources = self._makeDiaSources() + + @staticmethod + def _makeDiaSources(): + """Return an empty catalog with the fields the metrics read from the + detected sources. + """ + schema = lsst.afw.table.SourceTable.makeMinimalSchema() + schema.addField("x", "D", "Centroid x.") + schema.addField("y", "D", "Centroid y.") + schema.addField("isDipole", "Flag", "Is this source a dipole?") + schema.addField("dipoleAngle", "D", "Dipole orientation.") + schema.addField("dipoleLength", "D", "Dipole separation.") + return lsst.afw.table.SourceCatalog(schema) + + def _run(self, **kwargs): + """Run the metrics task on the images made in `setUp`.""" + config = lsst.ip.diffim.SpatiallySampledMetricsTask.ConfigClass() + config.update(**kwargs) + task = lsst.ip.diffim.SpatiallySampledMetricsTask(config=config) + return task.run(self.science, self.template, self.difference, self.diaSources, + self.psfMatchingKernel).spatiallySampledMetrics + + def testMaskPlaneNotInImage(self): + """A configured mask plane that is not registered has a fraction of + zero, rather than raising. + + An image converted from `lsst.images.DifferenceImage` only defines the + planes that have pixels set, so INJECTED is missing from an ordinary + AP run. + """ + present = "DETECTED" + missing = "NOT_A_MASK_PLANE" + self.assertIn(present, self.difference.mask.getMaskPlaneDict()) + self.assertNotIn(missing, self.difference.mask.getMaskPlaneDict()) + + metrics = self._run(metricsMaskPlanes=[present, missing]) + + self.assertGreater(len(metrics), 0) + self.assertTrue(np.all(metrics[f"{missing.lower()}_mask_fraction"] == 0)) + self.assertTrue(np.all(np.isfinite(metrics[f"{present.lower()}_mask_fraction"]))) + + +def setup_module(module): + lsst.utils.tests.init() + + +class MemoryTestCase(lsst.utils.tests.MemoryTestCase): + pass + + +if __name__ == "__main__": + lsst.utils.tests.init() + unittest.main() From 41e7364a07100716f679f745f8a8a42473225c8e Mon Sep 17 00:00:00 2001 From: Ian Sullivan Date: Mon, 5 Oct 2026 12:20:34 -0700 Subject: [PATCH 04/17] Convert the difference imaging tasks to optionally use `lsst.images` storage formats This adds a dependency on `lsst.images` --- python/lsst/ip/diffim/detectAndMeasure.py | 126 +++++- python/lsst/ip/diffim/getTemplate.py | 82 ++++ python/lsst/ip/diffim/subtractImages.py | 135 +++++- tests/test_detectAndMeasure.py | 396 +++++++++++++++++- tests/test_getTemplate.py | 215 +++++++++- tests/test_subtractTask.py | 481 ++++++++++++++++++++++ ups/ip_diffim.table | 1 + 7 files changed, 1430 insertions(+), 6 deletions(-) diff --git a/python/lsst/ip/diffim/detectAndMeasure.py b/python/lsst/ip/diffim/detectAndMeasure.py index 78ea68652..7fb6c8768 100644 --- a/python/lsst/ip/diffim/detectAndMeasure.py +++ b/python/lsst/ip/diffim/detectAndMeasure.py @@ -19,16 +19,21 @@ # You should have received a copy of the GNU General Public License # along with this program. If not, see . +import dataclasses import numpy as np import requests import os +import astropy.units + import lsst.afw.detection as afwDetection import lsst.afw.image as afwImage import lsst.afw.math as afwMath import lsst.afw.table as afwTable import lsst.daf.base as dafBase import lsst.geom +from lsst.images import DifferenceImage +from lsst.images.fields import field_from_legacy_background from lsst.ip.diffim.utils import (evaluateMaskFraction, computeDifferenceImageMetrics, populate_sattle_visit_cache) from lsst.meas.algorithms import SkyObjectsTask, SourceDetectionTask, SetPrimaryFlagsTask, MaskStreaksTask @@ -183,6 +188,10 @@ def __init__(self, *, config): self.outputs.remove("differenceBackground") if not (self.config.writeGlintInfo): self.outputs.remove("glintTrailInfo") + if self.config.image_type == "future": + self.difference = dataclasses.replace(self.difference, storageClass="DifferenceImage") + self.subtractedMeasuredExposure = dataclasses.replace(self.subtractedMeasuredExposure, + storageClass="DifferenceImage") class DetectAndMeasureConfig(pipeBase.PipelineTaskConfig, @@ -399,6 +408,16 @@ class DetectAndMeasureConfig(pipeBase.PipelineTaskConfig, "closest in time to the exposure." ) idGenerator = DetectorVisitIdGeneratorConfig.make_field() + image_type = pexConfig.ChoiceField[str]( + "Which image type to use for the measured difference image (and, for" + " DetectAndMeasureScoreTask, the measured score image).", + allowed={ + "legacy": "Write as a lsst.afw.image.ExposureF.", + "future": "Write as a lsst.images.DifferenceImage.", + }, + optional=False, + default="legacy", + ) def setDefaults(self): # Background subtraction @@ -592,10 +611,31 @@ def runQuantum(self, butlerQC: pipeBase.QuantumContext, inputs = butlerQC.get(inputRefs) idGenerator = self.config.idGenerator.apply(butlerQC.quantum.dataId) idFactory = idGenerator.make_table_id_factory() - # Specify the fields that `annotate` needs below, to ensure they - # exist, even as None. + # Assign the science image Detector to the output DifferenceImage + detector = inputs["science"].getDetector() + + templateInfo = None + matchingKernels = {} + if self.config.image_type == "future": + difference = inputs["difference"] + # Read the `TemplateInfo` from the difference while it is still in + # the `DifferenceImage` format, since legacy exposures don't carry + # it. + templateInfo = difference.templates + # Assign the kernel from the input difference image to the + # corresponding output image + matchingKernels["subtractedMeasuredExposure"] = difference.kernel + inputs["difference"] = difference.to_legacy() + if "scoreExposure" in inputs: + scoreExposure = inputs["scoreExposure"] + matchingKernels["scoreMeasuredExposure"] = scoreExposure.kernel + inputs["scoreExposure"] = scoreExposure.to_legacy() + + # Specify the fields that `annotate` and `convert_outputs_to_future` + # need below, to ensure they exist, even as None. measurementResults = pipeBase.Struct( subtractedMeasuredExposure=None, + scoreMeasuredExposure=None, diaSources=None, maskedStreaks=None, differenceBackground=None, @@ -603,18 +643,93 @@ def runQuantum(self, butlerQC: pipeBase.QuantumContext, try: self.run(**inputs, idFactory=idFactory, measurementResults=measurementResults) except pipeBase.AlgorithmError as e: + # Annotate before converting error = pipeBase.AnnotatedPartialOutputsError.annotate( e, self, measurementResults.subtractedMeasuredExposure, + measurementResults.scoreMeasuredExposure, measurementResults.diaSources, measurementResults.maskedStreaks, log=self.log ) + if self.config.image_type == "future": + self.convert_outputs_to_future(measurementResults, butlerQC.quantum.dataId, + detector=detector, matchingKernels=matchingKernels, + templateInfo=templateInfo) butlerQC.put(measurementResults, outputRefs) raise error from e + if self.config.image_type == "future": + self.convert_outputs_to_future(measurementResults, butlerQC.quantum.dataId, + detector=detector, matchingKernels=matchingKernels, + templateInfo=templateInfo) butlerQC.put(measurementResults, outputRefs) + def convert_outputs_to_future(self, results, data_id, detector=None, matchingKernels=None, + templateInfo=None): + """Convert the output images in a result struct to `lsst.images` types. + + This replaces ``results.subtractedMeasuredExposure`` and + ``results.scoreMeasuredExposure`` with `lsst.images.DifferenceImage` + instances, skipping either one that is absent or `None`. + + Parameters + ---------- + results : `lsst.pipe.base.Struct` + Output struct to read and modify in place. + data_id : `lsst.daf.butler.DataCoordinate` + The data ID of the images, which must have ``instrument`` and + ``visit`` keys. + detector : `lsst.afw.cameraGeom.Detector`, optional + Detector to set on an output image that has none; all of these + images are on the science image's pixel grid, so this should be + the science image's detector. + matchingKernels : `lsst.afw.math.Kernel` + Kernel used to PSF-match the subtracted template image. + templateInfo : `list` [`lsst.images.DifferenceImageTemplateInfo`] + Records of the coadds that went into the template, each holding + the second moments of that coadd's PSF and the region where it + overlapped the science image. + + Notes + ----- + This uses the default legacy mask plane map of + `lsst.images.get_legacy_difference_image_mask_planes` and adds any + planes from `~lsst.images.get_legacy_optional_mask_planes` that are set + for any pixels in the image. + """ + background = getattr(results, "differenceBackground", None) + for name in ("subtractedMeasuredExposure", "scoreMeasuredExposure"): + exposure = getattr(results, name, None) + if exposure is None: + continue + if exposure.getDetector() is None and detector is not None: + exposure.setDetector(detector) + image = DifferenceImage.from_legacy( + exposure, + unit=astropy.units.nJy, + instrument=data_id["instrument"], + visit=data_id["visit"], + ) + if templateInfo: + image.templates = templateInfo + if matchingKernels is not None: + image.kernel = matchingKernels[name] + if background is not None and len(background) > 0: + # The same background is subtracted from the difference image + # and the score image, so it applies to both. + image.backgrounds.add( + "subtracted", + field_from_legacy_background(background, bounds=image.bbox, unit=image.unit), + description="Background subtracted from this image before detection.", + is_subtracted=True, + ) + setattr(results, name, image) + if background is not None and len(background) == 0: + # An empty BackgroundList cannot be converted to a field, and means + # that no background was subtracted. + self.log.debug("No background model to attach to the output images.") + @timeMethod def run(self, science, matchedTemplate, difference, kernelSources=None, idFactory=None, measurementResults=None): @@ -1483,6 +1598,13 @@ class DetectAndMeasureScoreConnections(DetectAndMeasureConnections): name="{fakesType}{coaddName}Diff_scoreExp", ) + def __init__(self, *, config): + super().__init__(config=config) + if self.config.image_type == "future": + self.scoreExposure = dataclasses.replace(self.scoreExposure, storageClass="DifferenceImage") + self.scoreMeasuredExposure = dataclasses.replace(self.scoreMeasuredExposure, + storageClass="DifferenceImage") + class DetectAndMeasureScoreConfig(DetectAndMeasureConfig, pipelineConnections=DetectAndMeasureScoreConnections): diff --git a/python/lsst/ip/diffim/getTemplate.py b/python/lsst/ip/diffim/getTemplate.py index a7f3c66e3..e052610a5 100644 --- a/python/lsst/ip/diffim/getTemplate.py +++ b/python/lsst/ip/diffim/getTemplate.py @@ -21,6 +21,7 @@ import collections import dataclasses +import astropy.units from astropy.stats import gaussian_sigma_to_fwhm import numpy as np @@ -36,6 +37,7 @@ import lsst.pex.exceptions import lsst.pipe.base as pipeBase +from lsst.images import DifferenceImage, DifferenceImageTemplateInfo, get_legacy_template_mask_planes from lsst.skymap import BaseSkyMap from lsst.ip.diffim.dcrModel import DcrModel from lsst.ip.diffim.utils import evaluateMeanPsfFwhm, getPsfFwhm @@ -111,6 +113,8 @@ def __init__(self, *, config=None): self.coaddExposures, deferGraphConstraint=False, ) + if config.output_image_type == "future": + self.template = dataclasses.replace(self.template, storageClass="DifferenceImage") class GetTemplateConfig( @@ -182,6 +186,15 @@ class GetTemplateConfig( doc="If True, include the coadd dataset existence in the" " initial butler query during QuantumGraph generation.", ) + output_image_type = pexConfig.ChoiceField[str]( + "Which image type to use for the output template.", + allowed={ + "legacy": "Write as a lsst.afw.image.ExposureF.", + "future": "Write as a lsst.images.DifferenceImage.", + }, + optional=False, + default="legacy", + ) def setDefaults(self): # Use a smaller cache: per SeparableKernel.computeCache, this should @@ -254,8 +267,77 @@ def runQuantum(self, butlerQC, inputRefs, outputRefs): physical_filter=physical_filter, visit=outputRefs.template.dataId["visit"], ) + if self.config.output_image_type == "future": + outputs.template = self.convert_outputs_to_future(outputs.template, + butlerQC.quantum.dataId, + inputRefs.coaddExposures, + detector=detector) butlerQC.put(outputs, outputRefs) + def convert_outputs_to_future(self, template, data_id, coadd_refs, detector): + """Convert a template image to DifferenceImage format. + + This replaces ``template`` with an `lsst.images.DifferenceImage` + instance and records the coadds contributed to it. The Gaussian PSF + set by `run` is converted to an + `lsst.images.psfs.GaussianPointSpreadFunction`. + + Parameters + ---------- + template : `lsst.afw.image.Exposure` + Output struct to read and modify in place. Its ``template`` must + have a detector with per-amplifier raw geometry and a field angle + transform, which is why `runQuantum` sets the science image's + detector on it first. + data_id : `lsst.daf.butler.DataCoordinate` + Data ID of the science image the template was built for; supplies + the instrument name and visit id recorded on the output. + coadd_refs : `list` [`lsst.daf.butler.DatasetRef`] + References to the coadds that may have gone into the template. + May be a superset of the coadds that were used; supplies the + dataset id and RUN collection recorded for each coadd that was. + detector : `lsst.afw.image.Detector` + Detector of the science image that the template was built for. + + Returns + ------- + convertedTemplate : `lsst.images.DifferenceImage` + The converted image. + + Raises + ------ + RuntimeError + Raised if the template has pixels set in a mask plane that + `lsst.images.get_legacy_template_mask_planes` does not map, or if + ``coadd_refs`` spans more than one skymap. + """ + # The template PSF is a Gaussian, so rebuild the CoaddPsf from the + # template inputs to record which coadds contributed. + coaddPsf = self._makePsf(template, template.getInfo().getCoaddInputs().ccds, template.getWcs()) + template.setDetector(detector) + skymaps = {ref.dataId["skymap"] for ref in coadd_refs} + if len(skymaps) != 1: + raise RuntimeError("Expected the template inputs to come from one skymap;" + f" got {sorted(skymaps)}.") + butler_info = {(ref.dataId["tract"], ref.dataId["patch"]): (ref.id, ref.run) + for ref in coadd_refs} + convertedTemplate = DifferenceImage.from_legacy( + template, + unit=astropy.units.nJy, + plane_map=get_legacy_template_mask_planes(), + instrument=data_id["instrument"], + visit=data_id["visit"], + ) + del template + convertedTemplate.templates = DifferenceImageTemplateInfo.from_legacy_psf( + convertedTemplate.sky_projection.pixel_frame, + coaddPsf, + skymaps.pop(), + butler_info, + log=self.log, + ) + return convertedTemplate + def getExposures(self, coaddExposureHandles, bbox, skymap, wcs): """Return a data structure containing the coadds that overlap the specified bbox projected onto the sky, and a corresponding data diff --git a/python/lsst/ip/diffim/subtractImages.py b/python/lsst/ip/diffim/subtractImages.py index 25baf102a..de2a338dc 100644 --- a/python/lsst/ip/diffim/subtractImages.py +++ b/python/lsst/ip/diffim/subtractImages.py @@ -19,6 +19,8 @@ # You should have received a copy of the GNU General Public License # along with this program. If not, see . +import dataclasses + from astropy import units as u from astropy.stats import gaussian_fwhm_to_sigma import numpy as np @@ -27,6 +29,9 @@ import lsst.afw.image import lsst.afw.math import lsst.geom +from lsst.images import DifferenceImage +from lsst.images.convolution_kernels import ImageBasisConvolutionKernel +from lsst.images.fields import ChebyshevField from lsst.ip.diffim.utils import (evaluateMeanPsfFwhm, getPsfFwhm, computeDifferenceImageMetrics, computePSFNoiseEquivalentArea, checkMask, setSourceFootprints) @@ -69,7 +74,9 @@ class SubtractInputConnections(lsst.pipe.base.PipelineTaskConnections, dimensions=_dimensions, defaultTemplates=_defaultTemplates): template = connectionTypes.Input( - doc="Input warped template to subtract.", + doc="Input warped template to subtract." + " With image_type='future' this is read as a lsst.images.DifferenceImage, so that the" + " record of the coadds that went into it can be passed on to the outputs.", dimensions=("instrument", "visit", "detector"), storageClass="ExposureF", name="{fakesType}{coaddName}Diff_templateExp" @@ -100,6 +107,8 @@ def __init__(self, *, config=None): super().__init__(config=config) if not config.doApplyExternalCalibrations: del self.visitSummary + if config.image_type == "future": + self.template = dataclasses.replace(self.template, storageClass="DifferenceImage") class SubtractImageOutputConnections(lsst.pipe.base.PipelineTaskConnections, @@ -130,6 +139,13 @@ class SubtractImageOutputConnections(lsst.pipe.base.PipelineTaskConnections, name="{fakesType}{coaddName}Diff_psfMatchSources" ) + def __init__(self, *, config=None): + super().__init__(config=config) + if config.image_type == "future": + self.difference = dataclasses.replace(self.difference, storageClass="DifferenceImage") + self.matchedTemplate = dataclasses.replace(self.matchedTemplate, + storageClass="DifferenceImage") + class SubtractScoreOutputConnections(lsst.pipe.base.PipelineTaskConnections, dimensions=_dimensions, @@ -153,6 +169,11 @@ class SubtractScoreOutputConnections(lsst.pipe.base.PipelineTaskConnections, name="{fakesType}{coaddName}Diff_psfScoreMatchSources" ) + def __init__(self, *, config=None): + super().__init__(config=config) + if config.image_type == "future": + self.scoreExposure = dataclasses.replace(self.scoreExposure, storageClass="DifferenceImage") + class AlardLuptonSubtractConnections(SubtractInputConnections, SubtractImageOutputConnections): pass @@ -302,6 +323,16 @@ class AlardLuptonSubtractBaseConfig(lsst.pex.config.Config): doc="Re-run source detection for kernel candidates if an error is" " encountered while calculating the matching kernel." ) + image_type = lsst.pex.config.ChoiceField[str]( + doc="Which image type to use for the image outputs of this task " + "(difference and matchedTemplate, or scoreExposure).", + allowed={ + "legacy": "Write as a lsst.afw.image.ExposureF.", + "future": "Write as a lsst.images.DifferenceImage.", + }, + optional=False, + default="legacy", + ) def setDefaults(self): self.makeKernel.kernel.name = "AL" @@ -358,6 +389,16 @@ class AlardLuptonSubtractTask(lsst.pipe.base.PipelineTask): before kernel-matching. Subclasses that preconvolve override this to `True`.""" + futureImageOutputs = ("difference", "matchedTemplate") + """Names of the images in the results struct that will be converted, in the + same order as the output connections that carry them(`tuple` [`str`]). + Subclasses with different image outputs override this.""" + + futureBackgroundOutput = "difference" + """Name of the difference image in the results struct that the background + should be attached to (`str`). Subclasses that subtract it from + a different image override this.""" + def __init__(self, **kwargs): super().__init__(**kwargs) self.makeSubtask("decorrelate") @@ -416,16 +457,102 @@ def _applyExternalCalibrations(self, exposure, visitSummary): def runQuantum(self, butlerQC, inputRefs, outputRefs): inputs = butlerQC.get(inputRefs) + templateInfo = None + if self.config.image_type == "future": + template = inputs["template"] + # Read the `TemplateInfo` from the template while it is still in the + # `DifferenceImage` format, since legacy exposures don't carry it. + templateInfo = template.templates + inputs["template"] = template.to_legacy() try: results = self.run(**inputs) except lsst.pipe.base.AlgorithmError as e: error = lsst.pipe.base.AnnotatedPartialOutputsError.annotate(e, self, log=self.log) - # No partial outputs for butler to put + # No partial outputs for butler to put, so nothing to convert. raise error from e + if self.config.image_type == "future": + self.convert_outputs_to_future(results, butlerQC.quantum.dataId, templateInfo=templateInfo) + butlerQC.put(results, outputRefs) + def convert_outputs_to_future(self, results, data_id, templateInfo=None): + """Convert the image outputs to `lsst.images` types. + + Each image named in `futureImageOutputs` that is present on + ``results`` is replaced by an `lsst.images.DifferenceImage`, with the + PSF matching kernel attached to it, and the differential background + attached to the one named by `futureBackgroundOutput`. + + Parameters + ---------- + results : `lsst.pipe.base.Struct` + Output struct to read and modify in place. + data_id : `lsst.daf.butler.DataCoordinate` + The data ID of the science image. + templates : `list` [`lsst.images.DifferenceImageTemplateInfo`], + optional (for unit tests) + Record of the coadds that went into the template, taken from the + template this task subtracted. Attached to every converted + output. + + Raises + ------ + TypeError + Raised if the PSF matching kernel cannot be stored on an + `lsst.images.DifferenceImage`. + RuntimeError + Raised if an output has pixels set in a mask plane that + `lsst.images.get_legacy_difference_image_mask_planes` does not + map, or if an output has no detector and there is no science + image detector to take one from. + + Notes + ----- + This uses the default legacy mask plane map of + `lsst.images.get_legacy_difference_image_mask_planes` and adds any + planes from `~lsst.images.get_legacy_optional_mask_planes` that are set + for any pixels in the image. + """ + # The science image and every output built from it share the science + # pixel grid, so its detector applies to all of them. + science = results.matchedScience + try: + kernel = ImageBasisConvolutionKernel.from_legacy(results.psfMatchingKernel) + except (TypeError, ValueError, AttributeError) as e: + raise TypeError("The supplied PSF matching kernel cannot be attached to an" + f" lsst.images.DifferenceImage: {e}") from e + background = None + if self.config.doSubtractBackground: + background = ChebyshevField.from_legacy_function2(results.backgroundModel, unit=u.nJy) + for name in self.futureImageOutputs: + exposure = getattr(results, name, None) + if exposure is None: + continue + if exposure.getDetector() is None: + # Fall back on the science detector for images without a + # detector of their own (e.g. matchedTemplate) + exposure.setDetector(science.getDetector()) + image = DifferenceImage.from_legacy( + exposure, + unit=u.nJy, + instrument=data_id["instrument"], + visit=data_id["visit"], + ) + if kernel is not None: + image.kernel = kernel + if templateInfo: + image.templates = templateInfo + if background is not None and name == self.futureBackgroundOutput: + image.backgrounds.add( + "subtracted", background, + description="Differential background subtracted from this image when solving for" + " the PSF matching kernel.", + is_subtracted=True, + ) + setattr(results, name, image) + @timeMethod def run(self, template, science, sources, visitSummary=None): """PSF match, subtract, and decorrelate two images. @@ -1220,6 +1347,10 @@ class AlardLuptonPreconvolveSubtractTask(AlardLuptonSubtractTask): ConfigClass = AlardLuptonPreconvolveSubtractConfig _DefaultName = "alardLuptonPreconvolveSubtract" usePreconvolution = True + # This task just writes the score image instead of the difference image and + # matchedTemplate + futureImageOutputs = ("scoreExposure",) + futureBackgroundOutput = "scoreExposure" def run(self, template, science, sources, visitSummary=None): """Preconvolve the science image with its own PSF, diff --git a/tests/test_detectAndMeasure.py b/tests/test_detectAndMeasure.py index 58315d35d..feabd61c0 100644 --- a/tests/test_detectAndMeasure.py +++ b/tests/test_detectAndMeasure.py @@ -22,19 +22,28 @@ import numpy as np import os import unittest +import uuid from unittest import mock import requests +import astropy.units + import lsst.afw.detection as afwDetection import lsst.afw.geom as afwGeom import lsst.afw.image as afwImage import lsst.afw.table as afwTable import lsst.afw.math as afwMath import lsst.geom +from lsst.images import Box, DifferenceImage, DifferenceImageTemplateInfo, Polygon +from lsst.images.convolution_kernels import ImageBasisConvolutionKernel +import lsst.images.psfs +import lsst.images.serialization from lsst.ip.diffim import detectAndMeasure, subtractImages from lsst.afw.table import IdFactory -from lsst.afw.cameraGeom.testUtils import DetectorWrapper +from lsst.afw.cameraGeom.testUtils import CameraWrapper, DetectorWrapper import lsst.meas.algorithms as measAlg +import lsst.pipe.base +import lsst.pipe.base.testUtils from lsst.pipe.base import InvalidQuantumError, UpstreamFailureNoWorkFound, AlgorithmError import lsst.utils.tests import lsst.meas.base.tests @@ -1557,6 +1566,391 @@ def makeVisitInfo(id=1): ) +class DetectAndMeasureOutputFormatTest(DetectAndMeasureTestBase, lsst.utils.tests.TestCase): + """Tests of the conversion of the output images to `lsst.images` types. + """ + detectionTask = detectAndMeasure.DetectAndMeasureTask + + @classmethod + def setUpClass(cls): + super().setUpClass() + # `DifferenceImage.from_legacy` needs a detector with per-amplifier + # raw geometry and a FIELD_ANGLE transform; the simpler + # `DetectorWrapper` detector used elsewhere in this file does not + # have either, and fails to convert. + cls.detector = list(CameraWrapper().camera)[0] + + def setUp(self): + self.dataId = lsst.daf.butler.DataCoordinate.standardize( + instrument="I", + visit=42, + detector=12, + universe=lsst.daf.butler.DimensionUniverse(), + ) + + @staticmethod + def _makeTemplateInfo(patch=42): + """Return one record of a coadd that went into a template. + + Parameters + ---------- + patch : `int`, optional + Identifier of the patch the record describes. + + Returns + ------- + info : `lsst.images.DifferenceImageTemplateInfo` + The record. + """ + return DifferenceImageTemplateInfo( + skymap="skymap", tract=9813, patch=patch, dataset_id=uuid.uuid4(), + dataset_run="a/template/run", bounds=Polygon.from_box(Box.factory[0:50, 0:50]), + psf_shape_xx=4.0, psf_shape_yy=4.0, psf_shape_xy=0.0, psf_shape_flag=False) + + @staticmethod + def _makeKernel(scale=1.0): + """Return a PSF matching kernel of the one form that can be attached + to a difference image. + + Parameters + ---------- + scale : `float`, optional + Coefficient of the first basis kernel, to tell two kernels apart. + + Returns + ------- + kernel : `lsst.images.convolution_kernels.ImageBasisConvolutionKernel` + The kernel. + """ + basis = [] + for sigma in (1.0, 2.0): + image = afwImage.ImageD(5, 5) + y, x = np.mgrid[-2:3, -2:3] + image.array[:] = np.exp(-(x**2 + y**2)/(2*sigma**2)) + basis.append(afwMath.FixedKernel(image)) + bbox = lsst.geom.Box2D(lsst.geom.Box2I(lsst.geom.Point2I(0, 0), lsst.geom.Extent2I(100, 100))) + legacy = afwMath.LinearCombinationKernel(basis, afwMath.Chebyshev1Function2D(1, bbox)) + legacy.setSpatialParameters([[scale, 0.0, 0.0], [0.5, 0.0, 0.0]]) + return ImageBasisConvolutionKernel.from_legacy(legacy) + + def _makeFutureInput(self, exposure, templateInfo=None, kernel=None): + """Convert a legacy image to a `lsst.images.DifferenceImage`. + + Parameters + ---------- + exposure : `lsst.afw.image.ExposureF` + Image from `_make_images`. + templateInfo : `list` [`lsst.images.DifferenceImageTemplateInfo`], optional + Record of the coadds that went into the template. + kernel : `lsst.images.convolution_kernels.ImageBasisConvolutionKernel`, optional + PSF matching kernel. + + Returns + ------- + image : `lsst.images.DifferenceImage` + The converted image. + """ + image = DifferenceImage.from_legacy( + exposure, unit=astropy.units.nJy, + instrument=self.dataId["instrument"], visit=self.dataId["visit"]) + if templateInfo is not None: + image.templates = templateInfo + if kernel is not None: + image.kernel = kernel + return image + + def _make_images(self, withScore=False): + """Make a science/template/difference image set for detection. + + The template is not actually subtracted from the science image, so + that there are plenty of detections; configure the detection task + with ``badSubtractionRatioThreshold=1.`` to match. + + Parameters + ---------- + withScore : `bool`, optional + Also make and return a score (maximum likelihood) image? + + Returns + ------- + images : `tuple` + The science exposure, the source catalog, the matched template, + the difference image, and (if ``withScore``) the score image. + """ + kernelSize = 31 + kwargs = {"seed": 1, "psfSize": 2.4, "fluxLevel": 500, "x0": 12345, "y0": 67890, + "kernelSize": kernelSize, "templateBorderSize": kernelSize//2} + science, sources = makeTestImage(noiseLevel=1., noiseSeed=6, **kwargs) + matchedTemplate, _ = makeTestImage(noiseLevel=0.25, noiseSeed=7, **kwargs) + # `makeTestImage` does not set a detector, but conversion needs one. + science.setDetector(self.detector) + matchedTemplate.setDetector(self.detector) + difference = science.clone() + if not withScore: + return science, sources, matchedTemplate, difference + subtractTask = subtractImages.AlardLuptonPreconvolveSubtractTask() + score = subtractTask._convolveExposure(difference, science.psf.getKernel(), + subtractTask.convolutionControl) + score.setDetector(self.detector) + return science, sources, matchedTemplate, difference, score + + def _setup_task(self, **kwargs): + """Configure a detection task to go with `_make_images`. + + Parameters + ---------- + **kwargs + Any additional config parameters to set. + + Returns + ------- + task : `lsst.pipe.base.PipelineTask` + The configured task to use for detection and measurement. + """ + # The template is not actually subtracted from the science image in + # `_make_images`, so the subtraction residual metric threshold has to + # be set high. + kwargs.setdefault("badSubtractionRatioThreshold", 1.) + return self._setup_detection(doSkySources=False, **kwargs) + + def _check_converted(self, image, legacyBbox): + """Check the properties common to every converted output image. + + Parameters + ---------- + image : `lsst.images.DifferenceImage` + The converted image. + legacyBbox : `lsst.geom.Box2I` + Bounding box of the legacy exposure it was converted from. + """ + self.assertIsInstance(image, DifferenceImage) + self.assertEqual(image.bbox.x.start, legacyBbox.getMinX()) + self.assertEqual(image.bbox.x.stop, legacyBbox.getMaxX() + 1) + self.assertEqual(image.bbox.y.start, legacyBbox.getMinY()) + self.assertEqual(image.bbox.y.stop, legacyBbox.getMaxY() + 1) + self.assertEqual(image.unit, astropy.units.nJy) + self.assertIsNotNone(image.detector) + # Populating this needs the visit summary PhotoCalib, which this task + # does not have as an input, so it is deliberately left unset. + self.assertIsNone(image.photometric_scaling) + + def test_connections(self): + """Check that the output storage classes follow ``image_type``, + and that the background output is kept in both modes. + """ + for configClass, connectionsClass in ( + (detectAndMeasure.DetectAndMeasureConfig, detectAndMeasure.DetectAndMeasureConnections), + (detectAndMeasure.DetectAndMeasureScoreConfig, + detectAndMeasure.DetectAndMeasureScoreConnections), + ): + for imageType, storageClass in (("legacy", "ExposureF"), ("future", "DifferenceImage")): + with self.subTest(configClass=configClass.__name__, imageType=imageType): + config = configClass() + config.image_type = imageType + config.doSubtractBackground = True + config.doWriteBackground = True + connections = connectionsClass(config=config) + self.assertEqual(connections.subtractedMeasuredExposure.storageClass, storageClass) + if connectionsClass is detectAndMeasure.DetectAndMeasureScoreConnections: + self.assertEqual(connections.scoreMeasuredExposure.storageClass, storageClass) + # The background is a separate dataset in both modes. + self.assertIn("differenceBackground", connections.outputs) + self.assertEqual(connections.differenceBackground.storageClass, "Background") + + def test_run_returns_legacy_types(self): + """``run`` always returns legacy types, even in future mode.""" + science, sources, matchedTemplate, difference = self._make_images() + task = self._setup_task(image_type="future") + output = task.run(science, matchedTemplate, difference, sources) + self.assertIsInstance(output.subtractedMeasuredExposure, afwImage.ExposureF) + + def test_legacy_mode_does_not_convert(self): + """Legacy mode is unchanged, with no conversion.. + """ + science, sources, matchedTemplate, difference = self._make_images() + task = self._setup_task(doSubtractBackground=True, image_type="legacy") + output = task.run(science, matchedTemplate, difference, sources) + self.assertIsInstance(output.subtractedMeasuredExposure, afwImage.ExposureF) + self.assertIsInstance(output.differenceBackground, afwMath.BackgroundList) + lsst.pipe.base.testUtils.assertValidOutput(task, output) + + def test_convert_outputs_to_future(self): + """Convert the outputs of ``DetectAndMeasureTask``.""" + science, sources, matchedTemplate, difference = self._make_images() + task = self._setup_task(doSubtractBackground=True, image_type="future") + output = task.run(science, matchedTemplate, difference, sources) + legacyBbox = output.subtractedMeasuredExposure.getBBox() + background = output.differenceBackground + task.convert_outputs_to_future(output, self.dataId) + + self._check_converted(output.subtractedMeasuredExposure, legacyBbox) + self.assertIn("subtracted", output.subtractedMeasuredExposure.backgrounds) + # The separate background dataset is still written in future mode. + self.assertIs(output.differenceBackground, background) + self.assertIsInstance(output.differenceBackground, afwMath.BackgroundList) + + def test_convert_outputs_to_future_source_injection(self): + """Convert a difference image from a pipeline that injects sources. + + The injection planes are optional, so the conversion adds them for + an image that has pixels set in them instead of raising. + """ + planes = ("INJECTED", "INJECTED_CORE", "INJECTED_TEMPLATE", "INJECTED_CORE_TEMPLATE") + science, sources, matchedTemplate, difference = self._make_images() + task = self._setup_task(image_type="future") + output = task.run(science, matchedTemplate, difference, sources) + mask = output.subtractedMeasuredExposure.mask + for n, plane in enumerate(planes): + mask.addMaskPlane(plane) + mask.array[0, n] |= mask.getPlaneBitMask(plane) + + task.convert_outputs_to_future(output, self.dataId) + + image = output.subtractedMeasuredExposure + self.assertIsInstance(image, DifferenceImage) + for plane in planes: + with self.subTest(plane=plane): + self.assertEqual(np.count_nonzero(image.mask.get(plane)), 1) + + def test_convert_outputs_to_future_score(self): + """Convert the outputs of ``DetectAndMeasureScoreTask``.""" + self.detectionTask = detectAndMeasure.DetectAndMeasureScoreTask + science, sources, matchedTemplate, difference, score = self._make_images(withScore=True) + task = self._setup_task(doSubtractBackground=True, image_type="future") + output = task.run(science, matchedTemplate, difference, score, sources) + legacyBboxes = {name: getattr(output, name).getBBox() + for name in ("subtractedMeasuredExposure", "scoreMeasuredExposure")} + task.convert_outputs_to_future(output, self.dataId) + + for name, legacyBbox in legacyBboxes.items(): + with self.subTest(name=name): + self._check_converted(getattr(output, name), legacyBbox) + # The same background is subtracted from both images. + self.assertIn("subtracted", getattr(output, name).backgrounds) + + def test_convert_outputs_to_future_empty_background(self): + """An empty background list cannot be converted to a field, and must + not break the conversion of the images. + """ + science, sources, matchedTemplate, difference = self._make_images() + task = self._setup_task(doSubtractBackground=False, image_type="future") + output = task.run(science, matchedTemplate, difference, sources) + self.assertEqual(len(output.differenceBackground), 0) + legacyBbox = output.subtractedMeasuredExposure.getBBox() + task.convert_outputs_to_future(output, self.dataId) + self._check_converted(output.subtractedMeasuredExposure, legacyBbox) + self.assertEqual(len(output.subtractedMeasuredExposure.backgrounds), 0) + + def test_stray_struct_field_is_ignored(self): + """``DetectAndMeasureTask`` seeds ``scoreMeasuredExposure`` in its + result struct even though it has no such output connection; check + that such a stray field is ignored. + """ + science, sources, matchedTemplate, difference = self._make_images() + task = self._setup_task(image_type="legacy") + output = task.run(science, matchedTemplate, difference, sources) + output.scoreMeasuredExposure = difference + connections = task.config.ConnectionsClass(config=task.config) + self.assertNotIn("scoreMeasuredExposure", connections.outputs) + lsst.pipe.base.testUtils.assertValidOutput(task, output) + + def test_run_quantum(self): + """``runQuantum`` converts before putting in future mode, and leaves + the legacy types alone otherwise. + """ + templateInfo = [self._makeTemplateInfo()] + kernel = self._makeKernel() + for imageType, expectedType in (("legacy", afwImage.ExposureF), + ("future", DifferenceImage)): + with self.subTest(imageType=imageType): + science, sources, matchedTemplate, difference = self._make_images() + if imageType == "future": + difference = self._makeFutureInput(difference, templateInfo, kernel) + task = self._setup_task(doSubtractBackground=True, image_type=imageType) + butlerQC = _RecordingQuantumContext(self.dataId) + task.runQuantum(butlerQC, + _FakeRefs(science=science, matchedTemplate=matchedTemplate, + difference=difference, kernelSources=sources), + _FakeRefs()) + out = butlerQC.put_values.subtractedMeasuredExposure + self.assertIsInstance(out, expectedType) + # The background is a separate dataset in both modes. + self.assertIsInstance(butlerQC.put_values.differenceBackground, afwMath.BackgroundList) + if imageType == "future": + # The provenance of the image this was measured from + # survives, rather than being worked out again. + self.assertEqual(out.templates, templateInfo) + self.assertEqual(out.kernel, kernel) + + def test_run_quantum_score_task_provenance(self): + """Each output takes the kernel of the image it is measured from, + which for the score task are two different kernels. + """ + self.detectionTask = detectAndMeasure.DetectAndMeasureScoreTask + science, sources, matchedTemplate, difference, score = self._make_images(withScore=True) + templateInfo = [self._makeTemplateInfo()] + differenceKernel = self._makeKernel(scale=1.0) + scoreKernel = self._makeKernel(scale=3.0) + self.assertNotEqual(differenceKernel, scoreKernel) + task = self._setup_task(image_type="future") + butlerQC = _RecordingQuantumContext(self.dataId) + + task.runQuantum(butlerQC, + _FakeRefs(science=science, matchedTemplate=matchedTemplate, + difference=self._makeFutureInput(difference, templateInfo, + differenceKernel), + scoreExposure=self._makeFutureInput(score, templateInfo, scoreKernel), + kernelSources=sources), + _FakeRefs()) + + out = butlerQC.put_values.subtractedMeasuredExposure + scoreOut = butlerQC.put_values.scoreMeasuredExposure + self.assertEqual(out.templates, templateInfo) + self.assertEqual(scoreOut.templates, templateInfo) + self.assertEqual(out.kernel, differenceKernel) + self.assertEqual(scoreOut.kernel, scoreKernel) + + def test_connections_future_reads_difference_image(self): + """The inputs whose provenance is copied ask for the type that + carries it, and only in future mode. + """ + for imageType, expected in (("legacy", "ExposureF"), ("future", "DifferenceImage")): + with self.subTest(imageType=imageType): + config = detectAndMeasure.DetectAndMeasureScoreTask.ConfigClass() + config.image_type = imageType + connections = config.connections.ConnectionsClass(config=config) + self.assertEqual(connections.difference.storageClass, expected) + self.assertEqual(connections.scoreExposure.storageClass, expected) + # Nothing is copied off these, so the butler converts them. + self.assertEqual(connections.science.storageClass, "ExposureF") + self.assertEqual(connections.matchedTemplate.storageClass, "ExposureF") + + +class _FakeRefs: + """Stand-in for an input or output ``QuantizedConnection``.""" + def __init__(self, **kwargs): + self.__dict__.update(kwargs) + + +class _RecordingQuantumContext: + """Minimal `~lsst.pipe.base.QuantumContext` that records what was put. + + Parameters + ---------- + dataId : `lsst.daf.butler.DataCoordinate` + Data ID to report as the quantum data ID. + """ + def __init__(self, dataId): + self.quantum = _FakeRefs(dataId=dataId) + self.put_values = None + + def get(self, refs): + return dict(refs.__dict__) + + def put(self, values, refs): + self.put_values = values + + class MockResponse: """Provide a mock for requests.put calls""" def __init__(self, json_data, status_code, text): diff --git a/tests/test_getTemplate.py b/tests/test_getTemplate.py index 3b407eec5..ede0b58b0 100644 --- a/tests/test_getTemplate.py +++ b/tests/test_getTemplate.py @@ -24,21 +24,26 @@ import unittest import unittest.mock +import astropy.units as u from astropy.stats import gaussian_sigma_to_fwhm import numpy as np +import lsst.afw.cameraGeom.testUtils import lsst.afw.detection import lsst.afw.geom import lsst.afw.image import lsst.afw.math import lsst.afw.table -from lsst.daf.butler import DataCoordinate, DimensionUniverse +from lsst.daf.butler import DataCoordinate, DatasetRef, DatasetType, DimensionUniverse import lsst.geom +import lsst.images +import lsst.images.psfs import lsst.ip.diffim import lsst.meas.algorithms import lsst.meas.base.tests import lsst.pex.exceptions import lsst.pipe.base as pipeBase +import lsst.pipe.base.testUtils import lsst.skymap import lsst.utils.tests @@ -77,6 +82,7 @@ def setUp(self): self.skymap = self._makeSkymap() self.patches = collections.defaultdict(list) self.dataIds = collections.defaultdict(list) + self.coaddRefs = {} self.exposure = self._makeExposure() self.varianceBox = lsst.geom.Box2I(lsst.geom.Point2I(0, 0), lsst.geom.Point2I(180, 180)) @@ -638,6 +644,213 @@ def _configureAndRunTask(doScaleVariance, varianceScale=1.): factor*scaleFactor, rtol=1e-5) self.assertImagesAlmostEqual(templateLow.variance, templateOn.variance, rtol=1e-5) + def _runLegacyForFuture(self): + """Build a template from tract 0 with a task configured for the + future output type, and attach the detector that ``runQuantum`` reads + from the science image, but do not convert it. + + Returns + ------- + result : `lsst.pipe.base.Struct` + The legacy output struct, with a detector on the template. + box : `lsst.geom.Box2I` + The bounding box the template was requested on, before the task + grew it by the template border. + """ + config = lsst.ip.diffim.GetTemplateTask.ConfigClass() + config.output_image_type = "future" + task = lsst.ip.diffim.GetTemplateTask(config=config) + self.visit = 9876 + box = lsst.geom.Box2I(lsst.geom.Point2I(0, 0), lsst.geom.Point2I(180, 180)) + # Task modifies the input bbox, so pass a copy. + result = task.run(coaddExposureHandles={0: self.patches[0]}, + bbox=lsst.geom.Box2I(box), + wcs=self.exposure.wcs, + dataIds={0: self.dataIds[0]}, + physical_filter="a_test", + visit=self.visit) + # run() always produces legacy types, whatever output_image_type is. + self.assertIsInstance(result.template, lsst.afw.image.ExposureF) + self.futureTask = task + return result, box + + def _runFuture(self): + """Build a template and convert it, the way ``runQuantum`` does. + + Returns + ------- + result : `lsst.pipe.base.Struct` + The converted output struct. + box : `lsst.geom.Box2I` + The bounding box the template was requested on, before the task + grew it by the template border. + """ + # lsst.images needs per-amplifier raw geometry and a field angle + # transform, which the trivial test detectors do not have. + detector = list(lsst.afw.cameraGeom.testUtils.CameraWrapper().camera)[0] + result, box = self._runLegacyForFuture() + result.template = self.futureTask.convert_outputs_to_future( + result.template, {"instrument": "testCam", "visit": self.visit}, self._coaddRefs(0), + detector=detector) + return result, box, detector + + def _coaddRefs(self, tract): + """Return butler references for the coadds of one tract, like the + ones ``runQuantum`` passes on from its ``coaddExposures`` input. + + Parameters + ---------- + tract : `int` + Id of the tract whose coadds to make references for. + + Returns + ------- + refs : `list` [`lsst.daf.butler.DatasetRef`] + One reference per patch of that tract. The same references are + returned on every call, because each new one would get a new + dataset id. + """ + if tract not in self.coaddRefs: + datasetType = DatasetType("template_coadd", ("tract", "patch", "band", "skymap"), + "ExposureF", universe=DimensionUniverse()) + self.coaddRefs[tract] = [DatasetRef(datasetType, dataId, run="test_run") + for dataId in self.dataIds[tract]] + return self.coaddRefs[tract] + + def testConvertOutputsToFuture(self): + """Test that the output template is converted to a DifferenceImage + that keeps the pixel geometry and the detector it was built on. + """ + result, box, detector = self._runFuture() + + self.assertIsInstance(result.template, lsst.images.DifferenceImage) + self.assertEqual(result.template.unit, u.nJy) + self.assertIsNotNone(result.template.detector) + self.assertEqual(result.template.detector.name, detector.getName()) + # The template is grown by the border, so it is larger than both the + # requested box and the detector. + expectedBox = lsst.geom.Box2I(box) + expectedBox.grow(lsst.ip.diffim.GetTemplateTask.ConfigClass().templateBorderSize) + self.assertEqual(result.template.bbox.to_legacy(), expectedBox) + self.assertEqual(result.template.obs_info.visit_id, self.visit) + + def testConvertOutputsToFutureCoaddMaskPlanes(self): + """Test that the coadd mask planes a warped template carries are + converted, and come back when the image is converted to legacy. + + The template is converted with the template plane map, so the coadd + planes are named there; the source injection planes are optional, and + are added because this template has pixels set in them. + """ + detector = list(lsst.afw.cameraGeom.testUtils.CameraWrapper().camera)[0] + planes = ("CLIPPED", "REJECTED", "INEXACT_PSF", "SENSOR_EDGE", "HIGH_VARIANCE", + "INJECTED", "INJECTED_CORE") + result, _ = self._runLegacyForFuture() + mask = result.template.mask + for n, plane in enumerate(planes): + mask.addMaskPlane(plane) + mask.array[0, n] |= mask.getPlaneBitMask(plane) + + result.template = self.futureTask.convert_outputs_to_future( + result.template, {"instrument": "testCam", "visit": self.visit}, self._coaddRefs(0), + detector=detector) + + template = result.template + self.assertIsInstance(template, lsst.images.DifferenceImage) + for n, plane in enumerate(planes): + with self.subTest(plane=plane): + self.assertIn(plane, template.mask.schema.names) + self.assertEqual(np.count_nonzero(template.mask.get(plane)), 1) + # The butler converts a DifferenceImage to an ExposureF with the + # difference image plane map, which does not name the coadd planes; + # they keep their own names instead of being dropped. + legacy = template.to_legacy() + for n, plane in enumerate(planes): + with self.subTest(plane=plane): + self.assertIn(plane, legacy.mask.getMaskPlaneDict()) + self.assertEqual( + np.count_nonzero(legacy.mask.array & legacy.mask.getPlaneBitMask(plane)), 1) + + def testConvertOutputsToFutureTemplates(self): + """The coadds that went into the template are recorded on the + converted image. + """ + result, _, _ = self._runFuture() + refs = {(ref.dataId["tract"], ref.dataId["patch"]): ref for ref in self._coaddRefs(0)} + + templates = result.template.templates + self.assertGreater(len(templates), 0) + for template in templates: + with self.subTest(tract=template.tract, patch=template.patch): + ref = refs[template.tract, template.patch] + self.assertEqual(template.skymap, "skymap") + self.assertEqual(template.dataset_id, ref.id) + self.assertEqual(template.dataset_run, ref.run) + self.assertFalse(template.psf_shape_flag) + self.assertGreater(template.psf_shape_xx, 0) + self.assertGreater(template.bounds.area, 0) + + def testConvertOutputsToFuturePsf(self): + """The Gaussian PSF set by `run` is kept, defined over the detector. + """ + result, _ = self._runLegacyForFuture() + detector = list(lsst.afw.cameraGeom.testUtils.CameraWrapper().camera)[0] + legacyPsf = result.template.getPsf() + + result.template = self.futureTask.convert_outputs_to_future( + result.template, {"instrument": "testCam", "visit": self.visit}, self._coaddRefs(0), + detector=detector) + + psf = result.template.psf + self.assertIsInstance(psf, lsst.images.psfs.GaussianPointSpreadFunction) + self.assertFloatsAlmostEqual(psf.sigma, legacyPsf.getSigma()) + self.assertEqual(psf.bounds.to_legacy(), detector.getBBox()) + self.assertEqual(psf.kernel_bbox.to_legacy().getDimensions(), legacyPsf.getDimensions()) + + def testConvertOutputsToFuturePsfToLegacy(self): + """The Gaussian survives conversion back to an Exposure, which is how + a task that has not been converted reads the template. + """ + result, _, _ = self._runFuture() + + psf = result.template.to_legacy().getPsf() + self.assertIsInstance(psf, lsst.afw.detection.GaussianPsf) + self.assertFloatsAlmostEqual(psf.getSigma(), result.template.psf.sigma) + + +class GetTemplateConnectionsTestCase(lsst.utils.tests.TestCase): + """Test the connections that ``output_image_type`` switches on. + + These only need the config classes, not a built template. + """ + + def testOutputImageType(self): + """Test that the future output type changes the template storage + class, and leaves the input connections alone. + """ + Connections = lsst.ip.diffim.GetTemplateTask.ConfigClass.ConnectionsClass + config = lsst.ip.diffim.GetTemplateTask.ConfigClass() + + connections = Connections(config=config) + self.assertEqual(connections.template.storageClass, "ExposureF") + + config.output_image_type = "future" + connections = Connections(config=config) + self.assertEqual(connections.template.storageClass, "DifferenceImage") + # The dataset name is unchanged; only its storage class differs. + self.assertEqual(connections.template.name, "goodSeeingDiff_templateExp") + # The detector the conversion needs comes from the science image, so + # the future mode adds no input of its own. + self.assertEqual(connections.detector.name, "calexp.detector") + self.assertEqual(connections.detector.storageClass, "Detector") + + def testLintConnections(self): + """Check that the connections are self-consistent in both modes. + """ + for task in (lsst.ip.diffim.GetTemplateTask, lsst.ip.diffim.GetDcrTemplateTask): + with self.subTest(task=task.__name__): + lsst.pipe.base.testUtils.lintConnections(task.ConfigClass.ConnectionsClass) + def setup_module(module): lsst.utils.tests.init() diff --git a/tests/test_subtractTask.py b/tests/test_subtractTask.py index cf46b6794..18a67d99f 100644 --- a/tests/test_subtractTask.py +++ b/tests/test_subtractTask.py @@ -20,15 +20,18 @@ # along with this program. If not, see . import unittest +import uuid from astropy import units as u import lsst.afw.detection as afwDetection import lsst.afw.image as afwImage +import lsst.afw.cameraGeom.testUtils import lsst.afw.math as afwMath import lsst.afw.table as afwTable import lsst.geom import lsst.meas.algorithms as measAlg +from lsst.images import Box, DifferenceImage, DifferenceImageTemplateInfo, Polygon from lsst.ip.diffim import subtractImages, InsufficientKernelSourcesError from lsst.pex.config import FieldValidationError from lsst.pipe.base import NoWorkFound @@ -198,6 +201,36 @@ def test_clear_template_mask(self): else: self.assertTrue(diffimSum >= scienceSum) + def test_renamed_template_mask_survives_registered_plane(self): + """A renamed template plane keeps its pixels whether or not the plane + it is renamed to is already registered. + + afw's mask plane dictionary is shared by every mask in the process, so + reading an image that defines SAT_TEMPLATE puts that plane in the + template's dictionary as well, with no pixels set. + """ + for preRegistered in (False, True): + with self.subTest(preRegistered=preRegistered): + science, _ = makeTestImage(psfSize=self.midPsfSize, noiseLevel=1., noiseSeed=6) + template, _ = makeTestImage(psfSize=self.goodPsfSize, noiseLevel=1., noiseSeed=7, + templateBorderSize=20, doApplyCalibration=True) + template.mask.addMaskPlane("SAT") + satBitMask = template.mask.getPlaneBitMask("SAT") + template.mask.array[10:20, 10:20] |= satBitMask + expected = np.count_nonzero(template.mask.array & satBitMask) + self.assertGreater(expected, 0) + + if preRegistered: + afwImage.Mask(1, 1).addMaskPlane("SAT_TEMPLATE") + + task = self._setup_subtraction(mode="convolveTemplate") + task.updateMasks(template, science) + + satTemplateBitMask = template.mask.getPlaneBitMask("SAT_TEMPLATE") + self.assertEqual(np.count_nonzero(template.mask.array & satTemplateBitMask), expected) + # The plane the pixels were renamed from is still cleared. + self.assertEqual(np.count_nonzero(template.mask.array & satBitMask), 0) + def test_equal_images(self): """Test that running with enough sources produces reasonable output, with the same size psf in the template and science. @@ -1431,6 +1464,454 @@ def test_runSimplifiedTaskWithSourceDetection(self): self.assertFloatsAlmostEqual(differenceStd, np.sqrt(2)*noiseLevel, rtol=0.1) +class OutputImageTypeTest(lsst.utils.tests.TestCase): + """Tests of the optional conversion of the image outputs to + `lsst.images` types. + """ + # `makeTestImage` does not attach a detector or a VisitInfo, so supply a + # detector with real amplifier geometry and a simple data id. + dataId = {"instrument": "testCam", "visit": 12345} + + def setUp(self): + self.detector = list(lsst.afw.cameraGeom.testUtils.CameraWrapper().camera)[0] + + def _setup_subtraction(self, taskClass, **kwargs): + """Configure one of the subtraction tasks with the shared test + settings. + + Parameters + ---------- + taskClass : `lsst.pipe.base.PipelineTask` + The subtraction task class to configure. + **kwargs + Any additional config parameters to set. + + Returns + ------- + task : `lsst.pipe.base.PipelineTask` + The configured task. + """ + helper = AlardLuptonSubtractTestBase() + helper.subtractTask = taskClass + return helper._setup_subtraction(**kwargs) + + def _make_images(self): + """Make a science and template image pair to subtract. + + Returns + ------- + science : `lsst.afw.image.ExposureF` + The science image, with a detector attached. + template : `lsst.afw.image.ExposureF` + The template image, with no detector attached. + sources : `lsst.afw.table.SourceCatalog` + Sources detected on the science image. + """ + noiseLevel = 1. + science, sources = makeTestImage(psfSize=2.4, noiseLevel=noiseLevel, noiseSeed=6) + template, _ = makeTestImage(psfSize=2.0, noiseLevel=noiseLevel, noiseSeed=7, + templateBorderSize=20, doApplyCalibration=True) + science.setDetector(self.detector) + return science, template, sources + + def _check_converted(self, image, kernelExpected=True): + """Check that one converted image is a fully-formed DifferenceImage. + + Parameters + ---------- + image : `lsst.images.DifferenceImage` + The converted output to check. + kernelExpected : `bool`, optional + Should the PSF matching kernel be attached? + """ + self.assertIsInstance(image, DifferenceImage) + self.assertEqual(image.unit, u.nJy) + # `makeTestImage` does not attach a VisitInfo, so the data id supplied + # to the conversion is the only source of these. + frame = image.sky_projection.pixel_frame + self.assertEqual(frame.instrument, self.dataId["instrument"]) + self.assertEqual(frame.visit, self.dataId["visit"]) + self.assertEqual(image.detector.id, self.detector.getId()) + if kernelExpected: + self.assertIsNotNone(image.kernel) + else: + # An unset kernel raises, rather than returning None. + with self.assertRaises(AttributeError): + image.kernel + # Populating this needs the visitSummary PhotoCalib; it is deferred. + self.assertIsNone(image.photometric_scaling) + + def test_connections_legacy(self): + """Test that the default configuration leaves the output connections + with the legacy storage classes. + """ + for configClass, connectionsClass in ( + (subtractImages.AlardLuptonSubtractConfig, + subtractImages.AlardLuptonSubtractConnections), + (subtractImages.AlardLuptonPreconvolveSubtractConfig, + subtractImages.AlardLuptonPreconvolveSubtractConnections), + (subtractImages.SimplifiedSubtractConfig, + subtractImages.SimplifiedSubtractConnections), + ): + with self.subTest(configClass=configClass.__name__): + config = configClass() + self.assertEqual(config.image_type, "legacy") + connections = connectionsClass(config=config) + for name in connections.outputs: + if name in ("difference", "matchedTemplate", "scoreExposure"): + self.assertEqual(getattr(connections, name).storageClass, "ExposureF") + + def test_connections_future(self): + """Test that ``image_type="future"`` swaps the storage class of + every image output, and only those. + """ + for configClass, connectionsClass, expected in ( + (subtractImages.AlardLuptonSubtractConfig, + subtractImages.AlardLuptonSubtractConnections, + ("difference", "matchedTemplate")), + (subtractImages.AlardLuptonPreconvolveSubtractConfig, + subtractImages.AlardLuptonPreconvolveSubtractConnections, + ("scoreExposure",)), + (subtractImages.SimplifiedSubtractConfig, + subtractImages.SimplifiedSubtractConnections, + ("difference", "matchedTemplate")), + ): + with self.subTest(configClass=configClass.__name__): + config = configClass() + config.image_type = "future" + connections = connectionsClass(config=config) + for name in expected: + self.assertEqual(getattr(connections, name).storageClass, "DifferenceImage") + # The kernel has no converter to or from the new types, so it + # must remain a separate output in both modes. + if "psfMatchingKernel" in connections.outputs: + self.assertEqual(connections.psfMatchingKernel.storageClass, "MatchingKernel") + + def test_connections_future_simplified_existing_kernel(self): + """Test the storage class swap for the connection variant that + consumes an existing kernel. + """ + config = subtractImages.SimplifiedSubtractConfig() + config.image_type = "future" + config.useExistingKernel = True + connections = subtractImages.SimplifiedSubtractConnections(config=config) + self.assertEqual(connections.difference.storageClass, "DifferenceImage") + self.assertEqual(connections.matchedTemplate.storageClass, "DifferenceImage") + self.assertEqual(connections.inputPsfMatchingKernel.storageClass, "MatchingKernel") + self.assertNotIn("psfMatchingKernel", connections.outputs) + + def test_convert_outputs_to_future(self): + """Test converting the outputs of `AlardLuptonSubtractTask`. + """ + science, template, sources = self._make_images() + task = self._setup_subtraction(subtractImages.AlardLuptonSubtractTask, + image_type="future") + results = task.run(template, science, sources) + # The template has no detector, so neither does the matched template. + self.assertIsNone(results.matchedTemplate.getDetector()) + legacyDifference = results.difference.image.array.copy() + + task.convert_outputs_to_future(results, self.dataId) + + self._check_converted(results.difference) + self._check_converted(results.matchedTemplate) + self.assertEqual(results.difference.bbox.to_legacy(), science.getBBox()) + np.testing.assert_array_equal(results.difference.image.array, legacyDifference) + # The kernel outputs are untouched. + self.assertIsInstance(results.psfMatchingKernel, afwMath.LinearCombinationKernel) + self.assertIsInstance(results.kernelSources, afwTable.SourceCatalog) + + def test_convert_outputs_to_future_score(self): + """Test converting the outputs of + `AlardLuptonPreconvolveSubtractTask`, which writes a score image. + """ + science, template, sources = self._make_images() + task = self._setup_subtraction(subtractImages.AlardLuptonPreconvolveSubtractTask, + image_type="future") + results = task.run(template, science, sources) + + task.convert_outputs_to_future(results, self.dataId) + + self._check_converted(results.scoreExposure) + self.assertIsInstance(results.psfMatchingKernel, afwMath.LinearCombinationKernel) + + def test_convert_outputs_to_future_simplified(self): + """Test converting the outputs of `SimplifiedSubtractTask` running + with a pre-existing kernel. + """ + science, template, sources = self._make_images() + alTask = self._setup_subtraction(subtractImages.AlardLuptonSubtractTask) + alResults = alTask.run(template.clone(), science.clone(), sources) + task = self._setup_subtraction(subtractImages.SimplifiedSubtractTask, + useExistingKernel=True, + image_type="future") + results = task.run(template.clone(), science.clone(), + inputPsfMatchingKernel=alResults.psfMatchingKernel) + + task.convert_outputs_to_future(results, self.dataId) + + self._check_converted(results.difference) + self._check_converted(results.matchedTemplate) + + def test_convert_outputs_to_future_missing_output(self): + """Test that conversion skips outputs that are absent from the result + struct, instead of failing. + """ + science, template, sources = self._make_images() + task = self._setup_subtraction(subtractImages.AlardLuptonSubtractTask, + image_type="future") + results = task.run(template, science, sources) + del results.matchedTemplate + + task.convert_outputs_to_future(results, self.dataId) + + self._check_converted(results.difference) + self.assertFalse(hasattr(results, "matchedTemplate")) + + @staticmethod + def _makeTemplateInfo(tract=9813, patch=42): + """Return one record of a coadd that went into a template. + + Parameters + ---------- + tract, patch : `int`, optional + Identifiers of the coadd the record describes. + + Returns + ------- + info : `lsst.images.DifferenceImageTemplateInfo` + The record. + """ + return DifferenceImageTemplateInfo( + skymap="skymap", tract=tract, patch=patch, dataset_id=uuid.uuid4(), + dataset_run="a/template/run", bounds=Polygon.from_box(Box.factory[0:50, 0:50]), + psf_shape_xx=4.0, psf_shape_yy=4.0, psf_shape_xy=0.0, psf_shape_flag=False) + + def _make_future_template(self, template, templateInfo=None): + """Convert a legacy template to the type the ``future`` mode reads. + + Parameters + ---------- + template : `lsst.afw.image.ExposureF` + Template from `_make_images`, which has no detector. + templateInfo : `list` [`lsst.images.DifferenceImageTemplateInfo`], optional + Record of the coadds that went into it, left off if `None`. + + Returns + ------- + template : `lsst.images.DifferenceImage` + The converted template. + """ + template = template.clone() + template.setDetector(self.detector) + image = DifferenceImage.from_legacy(template, unit=u.nJy, instrument=self.dataId["instrument"], + visit=self.dataId["visit"]) + if templateInfo is not None: + image.templates = templateInfo + return image + + def _check_background(self, image, legacyBackground): + """Check that a converted image carries the background + that was subtracted from it. + + Parameters + ---------- + image : `lsst.images.DifferenceImage` + The converted image to check. + legacyBackground : `lsst.afw.math.Chebyshev1Function2D` + Background model from the result struct. + """ + self.assertEqual(list(image.backgrounds), ["subtracted"]) + field = image.backgrounds.subtracted.field + self.assertEqual(field.unit, u.nJy) + bbox = field.bounds.bbox.to_legacy() + rendered = field.render().array + for y in np.linspace(bbox.getMinY(), bbox.getMaxY(), 5, dtype=int): + for x in np.linspace(bbox.getMinX(), bbox.getMaxX(), 5, dtype=int): + self.assertFloatsAlmostEqual(rendered[y - bbox.getMinY(), x - bbox.getMinX()], + legacyBackground(float(x), float(y)), rtol=1e-12) + + def test_convert_outputs_to_future_background(self): + """The background goes on the difference image, which is + the one it was subtracted from. + """ + science, template, sources = self._make_images() + task = self._setup_subtraction(subtractImages.AlardLuptonSubtractTask, + doSubtractBackground=True, + image_type="future") + results = task.run(template, science, sources) + legacyBackground = results.backgroundModel + self.assertNotEqual(list(legacyBackground.getParameters()), [0.0]*len( + legacyBackground.getParameters())) + + task.convert_outputs_to_future(results, self.dataId) + + self._check_background(results.difference, legacyBackground) + self.assertEqual(results.difference.backgrounds.subtracted.field.bounds.bbox.to_legacy(), + science.getBBox()) + # The background was only taken out of the difference image. + self.assertEqual(len(results.matchedTemplate.backgrounds), 0) + + def test_convert_outputs_to_future_background_score(self): + """The score image gets the background, over the inner region that + the kernel was fit on. + """ + science, template, sources = self._make_images() + task = self._setup_subtraction(subtractImages.AlardLuptonPreconvolveSubtractTask, + doSubtractBackground=True, + image_type="future") + results = task.run(template, science, sources) + legacyBackground = results.backgroundModel + + task.convert_outputs_to_future(results, self.dataId) + + self._check_background(results.scoreExposure, legacyBackground) + bbox = results.scoreExposure.backgrounds.subtracted.field.bounds.bbox.to_legacy() + self.assertTrue(results.scoreExposure.bbox.to_legacy().contains(bbox)) + self.assertNotEqual(bbox, results.scoreExposure.bbox.to_legacy()) + + def test_convert_outputs_to_future_background_not_subtracted(self): + """A background the task fit but did not subtract is left off the + outputs. + """ + science, template, sources = self._make_images() + for taskClass, output in ((subtractImages.AlardLuptonSubtractTask, "difference"), + (subtractImages.AlardLuptonPreconvolveSubtractTask, "scoreExposure")): + with self.subTest(taskClass=taskClass.__name__): + task = self._setup_subtraction(taskClass, doSubtractBackground=False, + image_type="future") + results = task.run(template.clone(), science.clone(), sources) + + task.convert_outputs_to_future(results, self.dataId) + + self.assertEqual(len(getattr(results, output).backgrounds), 0) + + def test_convert_outputs_to_future_background_convolve_science(self): + """In ``convolveScience`` mode the actual background used is inverted. + """ + science, template, sources = self._make_images() + task = self._setup_subtraction(subtractImages.AlardLuptonSubtractTask, + mode="convolveScience", + doSubtractBackground=True, + image_type="future") + results = task.run(template, science, sources) + self.assertEqual(task.metadata["convolvedExposure"], "Science") + legacyBackground = results.backgroundModel + + task.convert_outputs_to_future(results, self.dataId) + + self._check_background(results.difference, legacyBackground) + + def test_convert_outputs_to_future_templates(self): + """The record of the coadds that went into the template is copied to + every converted output. + """ + science, template, sources = self._make_images() + task = self._setup_subtraction(subtractImages.AlardLuptonSubtractTask, + image_type="future") + results = task.run(template, science, sources) + templateInfo = [self._makeTemplateInfo(patch=41), self._makeTemplateInfo(patch=42)] + + task.convert_outputs_to_future(results, self.dataId, templateInfo=templateInfo) + + for image in (results.difference, results.matchedTemplate): + self.assertEqual([(t.tract, t.patch) for t in image.templates], + [(9813, 41), (9813, 42)]) + + def test_convert_outputs_to_future_no_templates(self): + """A template with no record of its coadds leaves the outputs without + one, rather than with an empty list. + """ + science, template, sources = self._make_images() + task = self._setup_subtraction(subtractImages.AlardLuptonSubtractTask, + image_type="future") + results = task.run(template, science, sources) + + task.convert_outputs_to_future(results, self.dataId) + + with self.assertRaises(AttributeError): + results.difference.templates + + def test_run_quantum_puts_future_types(self): + """``runQuantum`` converts before putting in future mode, and does not + convert in legacy mode. + """ + science, template, sources = self._make_images() + templateInfo = [self._makeTemplateInfo()] + for imageType, expectedType in (("legacy", afwImage.ExposureF), + ("future", DifferenceImage)): + with self.subTest(imageType=imageType): + task = self._setup_subtraction(subtractImages.AlardLuptonSubtractTask, + image_type=imageType) + if imageType == "future": + inputTemplate = self._make_future_template(template, templateInfo) + else: + inputTemplate = template.clone() + butlerQC = _RecordingQuantumContext(self.dataId) + task.runQuantum(butlerQC, + _FakeRefs(template=inputTemplate, science=science.clone(), + sources=sources), + _FakeRefs()) + self.assertIsInstance(butlerQC.put_values.difference, expectedType) + self.assertIsInstance(butlerQC.put_values.matchedTemplate, expectedType) + # The kernel is a separate output in both modes. + self.assertIsInstance(butlerQC.put_values.psfMatchingKernel, + afwMath.LinearCombinationKernel) + if imageType == "future": + self.assertEqual(butlerQC.put_values.difference.templates, templateInfo) + self.assertEqual(butlerQC.put_values.matchedTemplate.templates, templateInfo) + + def test_run_quantum_simplified_puts_future_types(self): + """``runQuantum`` converts the outputs of `SimplifiedSubtractTask`, + which loads an existing kernel rather than fitting one. + """ + science, template, sources = self._make_images() + templateInfo = [self._makeTemplateInfo()] + alTask = self._setup_subtraction(subtractImages.AlardLuptonSubtractTask) + alResults = alTask.run(template.clone(), science.clone(), sources) + for useExistingKernel in (True, False): + with self.subTest(useExistingKernel=useExistingKernel): + task = self._setup_subtraction(subtractImages.SimplifiedSubtractTask, + useExistingKernel=useExistingKernel, + image_type="future") + refs = {"template": self._make_future_template(template, templateInfo), + "science": science.clone()} + if useExistingKernel: + refs["inputPsfMatchingKernel"] = alResults.psfMatchingKernel + butlerQC = _RecordingQuantumContext(self.dataId) + + task.runQuantum(butlerQC, _FakeRefs(**refs), _FakeRefs()) + + self._check_converted(butlerQC.put_values.difference) + self._check_converted(butlerQC.put_values.matchedTemplate) + self.assertEqual(butlerQC.put_values.difference.templates, templateInfo) + + +class _FakeRefs: + """Stand-in for an input or output ``QuantizedConnection``.""" + def __init__(self, **kwargs): + self.__dict__.update(kwargs) + + +class _RecordingQuantumContext: + """Minimal `~lsst.pipe.base.QuantumContext` that records what was put. + + Parameters + ---------- + dataId : `dict` + Data ID to report as the quantum data ID. + """ + def __init__(self, dataId): + self.quantum = _FakeRefs(dataId=dataId) + self.put_values = None + + def get(self, refs): + return dict(refs.__dict__) + + def put(self, values, refs): + self.put_values = values + + def setup_module(module): lsst.utils.tests.init() diff --git a/ups/ip_diffim.table b/ups/ip_diffim.table index 12a20e86a..3802bd9d0 100644 --- a/ups/ip_diffim.table +++ b/ups/ip_diffim.table @@ -1,5 +1,6 @@ setupRequired(afw) setupRequired(daf_base) +setupRequired(images) setupRequired(log) setupRequired(meas_astrom) setupRequired(meas_algorithms) From 29ddb4f475af1225ecafe327d37170e343e1a043 Mon Sep 17 00:00:00 2001 From: Ian Sullivan Date: Sun, 27 Sep 2026 20:37:44 +0200 Subject: [PATCH 05/17] Use full exposure record for lsst.images conversion --- python/lsst/ip/diffim/getTemplate.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/python/lsst/ip/diffim/getTemplate.py b/python/lsst/ip/diffim/getTemplate.py index e052610a5..7af65d1be 100644 --- a/python/lsst/ip/diffim/getTemplate.py +++ b/python/lsst/ip/diffim/getTemplate.py @@ -271,10 +271,11 @@ def runQuantum(self, butlerQC, inputRefs, outputRefs): outputs.template = self.convert_outputs_to_future(outputs.template, butlerQC.quantum.dataId, inputRefs.coaddExposures, - detector=detector) + detector=detector, + exposureRecord=outputRefs.template) butlerQC.put(outputs, outputRefs) - def convert_outputs_to_future(self, template, data_id, coadd_refs, detector): + def convert_outputs_to_future(self, template, data_id, coadd_refs, detector, exposureRecord): """Convert a template image to DifferenceImage format. This replaces ``template`` with an `lsst.images.DifferenceImage` @@ -323,10 +324,9 @@ def convert_outputs_to_future(self, template, data_id, coadd_refs, detector): for ref in coadd_refs} convertedTemplate = DifferenceImage.from_legacy( template, + exposure_record=exposureRecord, unit=astropy.units.nJy, plane_map=get_legacy_template_mask_planes(), - instrument=data_id["instrument"], - visit=data_id["visit"], ) del template convertedTemplate.templates = DifferenceImageTemplateInfo.from_legacy_psf( From 4e0c0baf358abb4922fc0f0900e6c5bac4d17d63 Mon Sep 17 00:00:00 2001 From: Ian Sullivan Date: Sun, 27 Sep 2026 22:54:06 +0200 Subject: [PATCH 06/17] Read the obs_info to provide an exposureRecord for lsst.images conversions --- python/lsst/ip/diffim/detectAndMeasure.py | 26 ++++---- python/lsst/ip/diffim/getTemplate.py | 30 ++++++--- python/lsst/ip/diffim/subtractImages.py | 22 ++++--- python/lsst/ip/diffim/utils.py | 67 ++++++++++++++++++++ tests/test_detectAndMeasure.py | 22 ++++--- tests/test_getTemplate.py | 52 ++++++++++----- tests/test_subtractTask.py | 45 ++++++++----- tests/test_utils.py | 54 +++++++++++++++- tests/utils.py | 77 ++++++++++++++++++++++- 9 files changed, 326 insertions(+), 69 deletions(-) diff --git a/python/lsst/ip/diffim/detectAndMeasure.py b/python/lsst/ip/diffim/detectAndMeasure.py index 7fb6c8768..91e1f3e90 100644 --- a/python/lsst/ip/diffim/detectAndMeasure.py +++ b/python/lsst/ip/diffim/detectAndMeasure.py @@ -35,7 +35,7 @@ from lsst.images import DifferenceImage from lsst.images.fields import field_from_legacy_background from lsst.ip.diffim.utils import (evaluateMaskFraction, computeDifferenceImageMetrics, - populate_sattle_visit_cache) + populate_sattle_visit_cache, record_from_obs_info) from lsst.meas.algorithms import SkyObjectsTask, SourceDetectionTask, SetPrimaryFlagsTask, MaskStreaksTask from lsst.meas.algorithms import FindGlintTrailsTask, FindCosmicRaysConfig, findCosmicRays from lsst.meas.base import ForcedMeasurementTask, ApplyApCorrTask, DetectorVisitIdGeneratorConfig @@ -615,13 +615,16 @@ def runQuantum(self, butlerQC: pipeBase.QuantumContext, detector = inputs["science"].getDetector() templateInfo = None + exposureRecord = None matchingKernels = {} if self.config.image_type == "future": difference = inputs["difference"] - # Read the `TemplateInfo` from the difference while it is still in - # the `DifferenceImage` format, since legacy exposures don't carry - # it. + # Read the `TemplateInfo` and the observation metadata from the + # difference while it is still in the `DifferenceImage` format. templateInfo = difference.templates + dataId = butlerQC.quantum.dataId + exposureRecord = record_from_obs_info(difference.obs_info, dataId["instrument"], + dataId["visit"], dataId.universe) # Assign the kernel from the input difference image to the # corresponding output image matchingKernels["subtractedMeasuredExposure"] = difference.kernel @@ -654,18 +657,18 @@ def runQuantum(self, butlerQC: pipeBase.QuantumContext, log=self.log ) if self.config.image_type == "future": - self.convert_outputs_to_future(measurementResults, butlerQC.quantum.dataId, + self.convert_outputs_to_future(measurementResults, exposureRecord, detector=detector, matchingKernels=matchingKernels, templateInfo=templateInfo) butlerQC.put(measurementResults, outputRefs) raise error from e if self.config.image_type == "future": - self.convert_outputs_to_future(measurementResults, butlerQC.quantum.dataId, + self.convert_outputs_to_future(measurementResults, exposureRecord, detector=detector, matchingKernels=matchingKernels, templateInfo=templateInfo) butlerQC.put(measurementResults, outputRefs) - def convert_outputs_to_future(self, results, data_id, detector=None, matchingKernels=None, + def convert_outputs_to_future(self, results, exposureRecord, detector=None, matchingKernels=None, templateInfo=None): """Convert the output images in a result struct to `lsst.images` types. @@ -677,9 +680,9 @@ def convert_outputs_to_future(self, results, data_id, detector=None, matchingKer ---------- results : `lsst.pipe.base.Struct` Output struct to read and modify in place. - data_id : `lsst.daf.butler.DataCoordinate` - The data ID of the images, which must have ``instrument`` and - ``visit`` keys. + exposureRecord : `lsst.daf.butler.DimensionRecord` + The ``exposure`` record of the observation, which supplies the + observation metadata recorded on each converted output. detector : `lsst.afw.cameraGeom.Detector`, optional Detector to set on an output image that has none; all of these images are on the science image's pixel grid, so this should be @@ -708,8 +711,7 @@ def convert_outputs_to_future(self, results, data_id, detector=None, matchingKer image = DifferenceImage.from_legacy( exposure, unit=astropy.units.nJy, - instrument=data_id["instrument"], - visit=data_id["visit"], + exposure_record=exposureRecord, ) if templateInfo: image.templates = templateInfo diff --git a/python/lsst/ip/diffim/getTemplate.py b/python/lsst/ip/diffim/getTemplate.py index 7af65d1be..f21691efe 100644 --- a/python/lsst/ip/diffim/getTemplate.py +++ b/python/lsst/ip/diffim/getTemplate.py @@ -40,7 +40,7 @@ from lsst.images import DifferenceImage, DifferenceImageTemplateInfo, get_legacy_template_mask_planes from lsst.skymap import BaseSkyMap from lsst.ip.diffim.dcrModel import DcrModel -from lsst.ip.diffim.utils import evaluateMeanPsfFwhm, getPsfFwhm +from lsst.ip.diffim.utils import evaluateMeanPsfFwhm, getPsfFwhm, record_from_obs_info from lsst.meas.algorithms import CoaddPsf, CoaddPsfConfig, SubtractBackgroundTask, ScaleVarianceTask from lsst.utils.timer import timeMethod @@ -74,6 +74,13 @@ class GetTemplateConnections( storageClass="Detector", dimensions=("instrument", "visit", "detector"), ) + obs_info = pipeBase.connectionTypes.Input( + doc="Observation metadata of the exposure the template is built for." + " Only read with output_image_type='future'.", + name="{fakesType}calexp.obs_info", + storageClass="ObservationInfo", + dimensions=("instrument", "visit", "detector"), + ) wcs = pipeBase.connectionTypes.Input( doc="WCS of the exposure that we will construct the template for.", name="{fakesType}calexp.wcs", @@ -115,6 +122,8 @@ def __init__(self, *, config=None): ) if config.output_image_type == "future": self.template = dataclasses.replace(self.template, storageClass="DifferenceImage") + else: + del self.obs_info class GetTemplateConfig( @@ -251,6 +260,7 @@ def runQuantum(self, butlerQC, inputRefs, outputRefs): detector = inputs.pop("detector") bbox = detector.getBBox() wcs = inputs.pop("wcs") + obsInfo = inputs.pop("obs_info", None) coaddExposures = inputs.pop("coaddExposures") skymap = inputs.pop("skyMap") @@ -268,14 +278,16 @@ def runQuantum(self, butlerQC, inputRefs, outputRefs): visit=outputRefs.template.dataId["visit"], ) if self.config.output_image_type == "future": + dataId = butlerQC.quantum.dataId + exposureRecord = record_from_obs_info(obsInfo, dataId["instrument"], dataId["visit"], + dataId.universe) outputs.template = self.convert_outputs_to_future(outputs.template, - butlerQC.quantum.dataId, inputRefs.coaddExposures, detector=detector, - exposureRecord=outputRefs.template) + exposureRecord=exposureRecord) butlerQC.put(outputs, outputRefs) - def convert_outputs_to_future(self, template, data_id, coadd_refs, detector, exposureRecord): + def convert_outputs_to_future(self, template, coadd_refs, detector, exposureRecord): """Convert a template image to DifferenceImage format. This replaces ``template`` with an `lsst.images.DifferenceImage` @@ -290,15 +302,16 @@ def convert_outputs_to_future(self, template, data_id, coadd_refs, detector, exp have a detector with per-amplifier raw geometry and a field angle transform, which is why `runQuantum` sets the science image's detector on it first. - data_id : `lsst.daf.butler.DataCoordinate` - Data ID of the science image the template was built for; supplies - the instrument name and visit id recorded on the output. coadd_refs : `list` [`lsst.daf.butler.DatasetRef`] References to the coadds that may have gone into the template. May be a superset of the coadds that were used; supplies the dataset id and RUN collection recorded for each coadd that was. detector : `lsst.afw.image.Detector` Detector of the science image that the template was built for. + exposureRecord : `lsst.daf.butler.DimensionRecord` + The ``exposure`` record of the observation the template was built + for, which supplies the observation metadata recorded on the + output. Returns ------- @@ -994,6 +1007,9 @@ def runQuantum(self, butlerQC, inputRefs, outputRefs): dcrCoaddExposureHandles = inputs.pop("dcrCoadds") skymap = inputs.pop("skyMap") visitInfo = inputs.pop("visitInfo") + # This task never converts its output, so the observation metadata + # that the conversion needs is dropped here. + inputs.pop("obs_info", None) # This should not happen with a properly configured execution context. assert not inputs, "runQuantum got more inputs than expected" diff --git a/python/lsst/ip/diffim/subtractImages.py b/python/lsst/ip/diffim/subtractImages.py index de2a338dc..d559fcd90 100644 --- a/python/lsst/ip/diffim/subtractImages.py +++ b/python/lsst/ip/diffim/subtractImages.py @@ -34,7 +34,7 @@ from lsst.images.fields import ChebyshevField from lsst.ip.diffim.utils import (evaluateMeanPsfFwhm, getPsfFwhm, computeDifferenceImageMetrics, computePSFNoiseEquivalentArea, - checkMask, setSourceFootprints) + checkMask, record_from_obs_info, setSourceFootprints) from lsst.meas.algorithms import ScaleVarianceTask, ScienceSourceSelectorTask import lsst.pex.config import lsst.pipe.base @@ -458,11 +458,15 @@ def _applyExternalCalibrations(self, exposure, visitSummary): def runQuantum(self, butlerQC, inputRefs, outputRefs): inputs = butlerQC.get(inputRefs) templateInfo = None + exposureRecord = None if self.config.image_type == "future": template = inputs["template"] - # Read the `TemplateInfo` from the template while it is still in the - # `DifferenceImage` format, since legacy exposures don't carry it. + # Read the `TemplateInfo` and the observation metadata from the + # template while it is still in the `DifferenceImage` format. templateInfo = template.templates + dataId = butlerQC.quantum.dataId + exposureRecord = record_from_obs_info(template.obs_info, dataId["instrument"], + dataId["visit"], dataId.universe) inputs["template"] = template.to_legacy() try: @@ -473,11 +477,11 @@ def runQuantum(self, butlerQC, inputRefs, outputRefs): raise error from e if self.config.image_type == "future": - self.convert_outputs_to_future(results, butlerQC.quantum.dataId, templateInfo=templateInfo) + self.convert_outputs_to_future(results, exposureRecord, templateInfo=templateInfo) butlerQC.put(results, outputRefs) - def convert_outputs_to_future(self, results, data_id, templateInfo=None): + def convert_outputs_to_future(self, results, exposureRecord, templateInfo=None): """Convert the image outputs to `lsst.images` types. Each image named in `futureImageOutputs` that is present on @@ -489,8 +493,9 @@ def convert_outputs_to_future(self, results, data_id, templateInfo=None): ---------- results : `lsst.pipe.base.Struct` Output struct to read and modify in place. - data_id : `lsst.daf.butler.DataCoordinate` - The data ID of the science image. + exposureRecord : `lsst.daf.butler.DimensionRecord` + The ``exposure`` record of the observation, which supplies the + observation metadata recorded on each converted output. templates : `list` [`lsst.images.DifferenceImageTemplateInfo`], optional (for unit tests) Record of the coadds that went into the template, taken from the @@ -537,8 +542,7 @@ def convert_outputs_to_future(self, results, data_id, templateInfo=None): image = DifferenceImage.from_legacy( exposure, unit=u.nJy, - instrument=data_id["instrument"], - visit=data_id["visit"], + exposure_record=exposureRecord, ) if kernel is not None: image.kernel = kernel diff --git a/python/lsst/ip/diffim/utils.py b/python/lsst/ip/diffim/utils.py index 6e3823d82..de2296083 100644 --- a/python/lsst/ip/diffim/utils.py +++ b/python/lsst/ip/diffim/utils.py @@ -24,6 +24,7 @@ __all__ = ["evaluateMeanPsfFwhm", "getPsfFwhm", "getKernelCenterDisplacement", "computeDifferenceImageMetrics", "checkMask", "setSourceFootprints", + "record_from_obs_info", ] import itertools @@ -31,11 +32,13 @@ import requests from astropy.stats import gaussian_sigma_to_fwhm +import astropy.units as u import numpy as np import lsst.afw.detection as afwDetection import lsst.afw.image as afwImage import lsst.afw.math as afwMath +from lsst.daf.butler import Timespan import lsst.geom as geom from lsst.pex.exceptions import InvalidParameterError, RangeError import lsst.pipe.base @@ -45,6 +48,70 @@ _LOG = getLogger(__name__) +def record_from_obs_info(obs_info, instrument, visit, universe): + """Build an ``exposure`` dimension record from an observation info. + + Parameters + ---------- + obs_info : `astro_metadata_translator.ObservationInfo` + Observation info of an image that the image to be converted is + derived from. + instrument : `str` + Name of the instrument, from the pixel frame of the image + ``obs_info`` came from. + visit : `int` + Id of the visit, from the pixel frame of the image ``obs_info`` + came from. + universe : `lsst.daf.butler.DimensionUniverse` + Dimension universe that defines the record's schema, usually + ``butlerQC.quantum.dataId.universe``. + + Returns + ------- + record : `lsst.daf.butler.DimensionRecord` + The ``exposure`` record of the observation. + + Raises + ------ + ValueError + Raised if ``obs_info`` is `None`, which means the image it came from + was not in an `lsst.images` format. + + Notes + ----- + The instrument and the visit are passed separately because an + observation info does not carry them in the form the conversion needs: + its ``instrument`` is read from the legacy FITS header rather than from + the record, and its ``exposure_id`` is the id that the visit info held. + Taking both from the pixel frame of the image that supplied + ``obs_info`` puts the converted image on the same instrument and visit + as the image it came from. + """ + if obs_info is None: + raise ValueError("Cannot rebuild an exposure record without an observation info.") + azimuth = zenithAngle = None + if obs_info.altaz_begin is not None: + azimuth = obs_info.altaz_begin.az.deg + zenithAngle = 90.0 - obs_info.altaz_begin.alt.deg + exposureTime = obs_info.exposure_time_requested + return universe["exposure"].RecordClass( + instrument=instrument, + id=visit, + obs_id=obs_info.observation_id, + group=obs_info.exposure_group, + day_obs=obs_info.observing_day, + physical_filter=obs_info.physical_filter, + exposure_time=None if exposureTime is None else exposureTime.to_value(u.s), + seq_num=obs_info.observation_counter, + seq_start=obs_info.group_counter_start, + seq_end=obs_info.group_counter_end, + can_see_sky=obs_info.can_see_sky, + azimuth=azimuth, + zenith_angle=zenithAngle, + timespan=Timespan(begin=obs_info.datetime_begin, end=obs_info.datetime_end), + ) + + def getKernelCenterDisplacement(kernel, x, y, image=None): """Calculate the PSF matching kernel peak offset from the nominal position. diff --git a/tests/test_detectAndMeasure.py b/tests/test_detectAndMeasure.py index feabd61c0..f06c99b36 100644 --- a/tests/test_detectAndMeasure.py +++ b/tests/test_detectAndMeasure.py @@ -52,7 +52,7 @@ import lsst.geom as geom import lsst.pex.config as pexConfig -from utils import makeTestImage, checkMask +from utils import makeTestExposureRecord, makeTestImage, makeTestVisitInfo, checkMask class DetectAndMeasureTestBase: @@ -1587,6 +1587,11 @@ def setUp(self): detector=12, universe=lsst.daf.butler.DimensionUniverse(), ) + # The exposure record that `runQuantum` rebuilds from the + # observation metadata of the input difference image. + self.exposureRecord = makeTestExposureRecord(self.dataId.universe, + instrument=self.dataId["instrument"], + visit=self.dataId["visit"]) @staticmethod def _makeTemplateInfo(patch=42): @@ -1651,8 +1656,7 @@ def _makeFutureInput(self, exposure, templateInfo=None, kernel=None): The converted image. """ image = DifferenceImage.from_legacy( - exposure, unit=astropy.units.nJy, - instrument=self.dataId["instrument"], visit=self.dataId["visit"]) + exposure, unit=astropy.units.nJy, exposure_record=self.exposureRecord) if templateInfo is not None: image.templates = templateInfo if kernel is not None: @@ -1682,9 +1686,11 @@ def _make_images(self, withScore=False): "kernelSize": kernelSize, "templateBorderSize": kernelSize//2} science, sources = makeTestImage(noiseLevel=1., noiseSeed=6, **kwargs) matchedTemplate, _ = makeTestImage(noiseLevel=0.25, noiseSeed=7, **kwargs) - # `makeTestImage` does not set a detector, but conversion needs one. + # Add components to the basic test image science.setDetector(self.detector) matchedTemplate.setDetector(self.detector) + for exposure in (science, matchedTemplate): + exposure.info.setVisitInfo(makeTestVisitInfo(self.dataId["visit"])) difference = science.clone() if not withScore: return science, sources, matchedTemplate, difference @@ -1781,7 +1787,7 @@ def test_convert_outputs_to_future(self): output = task.run(science, matchedTemplate, difference, sources) legacyBbox = output.subtractedMeasuredExposure.getBBox() background = output.differenceBackground - task.convert_outputs_to_future(output, self.dataId) + task.convert_outputs_to_future(output, self.exposureRecord) self._check_converted(output.subtractedMeasuredExposure, legacyBbox) self.assertIn("subtracted", output.subtractedMeasuredExposure.backgrounds) @@ -1804,7 +1810,7 @@ def test_convert_outputs_to_future_source_injection(self): mask.addMaskPlane(plane) mask.array[0, n] |= mask.getPlaneBitMask(plane) - task.convert_outputs_to_future(output, self.dataId) + task.convert_outputs_to_future(output, self.exposureRecord) image = output.subtractedMeasuredExposure self.assertIsInstance(image, DifferenceImage) @@ -1820,7 +1826,7 @@ def test_convert_outputs_to_future_score(self): output = task.run(science, matchedTemplate, difference, score, sources) legacyBboxes = {name: getattr(output, name).getBBox() for name in ("subtractedMeasuredExposure", "scoreMeasuredExposure")} - task.convert_outputs_to_future(output, self.dataId) + task.convert_outputs_to_future(output, self.exposureRecord) for name, legacyBbox in legacyBboxes.items(): with self.subTest(name=name): @@ -1837,7 +1843,7 @@ def test_convert_outputs_to_future_empty_background(self): output = task.run(science, matchedTemplate, difference, sources) self.assertEqual(len(output.differenceBackground), 0) legacyBbox = output.subtractedMeasuredExposure.getBBox() - task.convert_outputs_to_future(output, self.dataId) + task.convert_outputs_to_future(output, self.exposureRecord) self._check_converted(output.subtractedMeasuredExposure, legacyBbox) self.assertEqual(len(output.subtractedMeasuredExposure.backgrounds), 0) diff --git a/tests/test_getTemplate.py b/tests/test_getTemplate.py index ede0b58b0..731069d5e 100644 --- a/tests/test_getTemplate.py +++ b/tests/test_getTemplate.py @@ -38,6 +38,7 @@ import lsst.geom import lsst.images import lsst.images.psfs +from lsst.images.tests import compare_masked_image_to_legacy import lsst.ip.diffim import lsst.meas.algorithms import lsst.meas.base.tests @@ -47,6 +48,8 @@ import lsst.skymap import lsst.utils.tests +from utils import makeTestExposureRecord + from utils import generate_data_id # Change this to True, `setup display_ds9`, and open ds9 (or use another afw @@ -690,10 +693,16 @@ def _runFuture(self): detector = list(lsst.afw.cameraGeom.testUtils.CameraWrapper().camera)[0] result, box = self._runLegacyForFuture() result.template = self.futureTask.convert_outputs_to_future( - result.template, {"instrument": "testCam", "visit": self.visit}, self._coaddRefs(0), - detector=detector) + result.template, self._coaddRefs(0), detector=detector, + exposureRecord=self._exposureRecord()) return result, box, detector + def _exposureRecord(self): + """Return the exposure record that ``runQuantum`` rebuilds from the + observation metadata of the science image. + """ + return makeTestExposureRecord(DimensionUniverse(), visit=self.visit) + def _coaddRefs(self, tract): """Return butler references for the coadds of one tract, like the ones ``runQuantum`` passes on from its ``coaddExposures`` input. @@ -719,11 +728,18 @@ def _coaddRefs(self, tract): def testConvertOutputsToFuture(self): """Test that the output template is converted to a DifferenceImage - that keeps the pixel geometry and the detector it was built on. + that keeps the pixels and the detector it was built on. """ - result, box, detector = self._runFuture() + detector = list(lsst.afw.cameraGeom.testUtils.CameraWrapper().camera)[0] + result, box = self._runLegacyForFuture() + legacy = result.template.clone() + result.template = self.futureTask.convert_outputs_to_future( + result.template, self._coaddRefs(0), detector=detector, + exposureRecord=self._exposureRecord()) self.assertIsInstance(result.template, lsst.images.DifferenceImage) + compare_masked_image_to_legacy(result.template, legacy.maskedImage, + plane_map=lsst.images.get_legacy_template_mask_planes()) self.assertEqual(result.template.unit, u.nJy) self.assertIsNotNone(result.template.detector) self.assertEqual(result.template.detector.name, detector.getName()) @@ -750,17 +766,21 @@ def testConvertOutputsToFutureCoaddMaskPlanes(self): for n, plane in enumerate(planes): mask.addMaskPlane(plane) mask.array[0, n] |= mask.getPlaneBitMask(plane) + legacyTemplate = result.template.clone() result.template = self.futureTask.convert_outputs_to_future( - result.template, {"instrument": "testCam", "visit": self.visit}, self._coaddRefs(0), - detector=detector) + result.template, self._coaddRefs(0), detector=detector, + exposureRecord=self._exposureRecord()) template = result.template self.assertIsInstance(template, lsst.images.DifferenceImage) - for n, plane in enumerate(planes): - with self.subTest(plane=plane): - self.assertIn(plane, template.mask.schema.names) - self.assertEqual(np.count_nonzero(template.mask.get(plane)), 1) + # Name the optional planes that have pixels set in the map, so that + # dropping one fails the comparison instead of skipping it. The others + # are left out because they may be registered by unrelated tests. + optional = lsst.images.get_legacy_optional_mask_planes() + planeMap = lsst.images.get_legacy_template_mask_planes() + planeMap.update({plane: optional[plane] for plane in planes if plane in optional}) + compare_masked_image_to_legacy(template, legacyTemplate.maskedImage, plane_map=planeMap) # The butler converts a DifferenceImage to an ExposureF with the # difference image plane map, which does not name the coadd planes; # they keep their own names instead of being dropped. @@ -798,8 +818,8 @@ def testConvertOutputsToFuturePsf(self): legacyPsf = result.template.getPsf() result.template = self.futureTask.convert_outputs_to_future( - result.template, {"instrument": "testCam", "visit": self.visit}, self._coaddRefs(0), - detector=detector) + result.template, self._coaddRefs(0), detector=detector, + exposureRecord=self._exposureRecord()) psf = result.template.psf self.assertIsInstance(psf, lsst.images.psfs.GaussianPointSpreadFunction) @@ -833,16 +853,20 @@ def testOutputImageType(self): connections = Connections(config=config) self.assertEqual(connections.template.storageClass, "ExposureF") + self.assertNotIn("obs_info", connections.inputs) config.output_image_type = "future" connections = Connections(config=config) self.assertEqual(connections.template.storageClass, "DifferenceImage") # The dataset name is unchanged; only its storage class differs. self.assertEqual(connections.template.name, "goodSeeingDiff_templateExp") - # The detector the conversion needs comes from the science image, so - # the future mode adds no input of its own. + # The detector the conversion needs comes from the science image. self.assertEqual(connections.detector.name, "calexp.detector") self.assertEqual(connections.detector.storageClass, "Detector") + # The observation metadata is only read in future mode, where it is + # the one input the conversion adds. + self.assertEqual(connections.obs_info.name, "calexp.obs_info") + self.assertEqual(connections.obs_info.storageClass, "ObservationInfo") def testLintConnections(self): """Check that the connections are self-consistent in both modes. diff --git a/tests/test_subtractTask.py b/tests/test_subtractTask.py index 18a67d99f..ab9262df9 100644 --- a/tests/test_subtractTask.py +++ b/tests/test_subtractTask.py @@ -31,6 +31,7 @@ import lsst.afw.table as afwTable import lsst.geom import lsst.meas.algorithms as measAlg +from lsst.daf.butler import DataCoordinate, DimensionUniverse from lsst.images import Box, DifferenceImage, DifferenceImageTemplateInfo, Polygon from lsst.ip.diffim import subtractImages, InsufficientKernelSourcesError from lsst.pex.config import FieldValidationError @@ -41,7 +42,8 @@ evaluateMeanPsfFwhm, getPsfFwhm) from lsst.pex.exceptions import InvalidParameterError -from utils import makeStats, makeTestImage, CustomCoaddPsf +from utils import (makeStats, makeTestExposureRecord, makeTestImage, makeTestVisitInfo, + CustomCoaddPsf) class AlardLuptonSubtractTestBase: @@ -1468,12 +1470,20 @@ class OutputImageTypeTest(lsst.utils.tests.TestCase): """Tests of the optional conversion of the image outputs to `lsst.images` types. """ - # `makeTestImage` does not attach a detector or a VisitInfo, so supply a - # detector with real amplifier geometry and a simple data id. - dataId = {"instrument": "testCam", "visit": 12345} - def setUp(self): + # `makeTestImage` does not attach a detector or a VisitInfo, so + # supply a detector with real amplifier geometry, a data id, and the + # exposure record the conversion takes its metadata from. self.detector = list(lsst.afw.cameraGeom.testUtils.CameraWrapper().camera)[0] + self.dataId = DataCoordinate.standardize( + instrument="testCam", + visit=12345, + detector=self.detector.getId(), + universe=DimensionUniverse(), + ) + self.exposureRecord = makeTestExposureRecord(self.dataId.universe, + instrument=self.dataId["instrument"], + visit=self.dataId["visit"]) def _setup_subtraction(self, taskClass, **kwargs): """Configure one of the subtraction tasks with the shared test @@ -1512,6 +1522,8 @@ def _make_images(self): template, _ = makeTestImage(psfSize=2.0, noiseLevel=noiseLevel, noiseSeed=7, templateBorderSize=20, doApplyCalibration=True) science.setDetector(self.detector) + for exposure in (science, template): + exposure.info.setVisitInfo(makeTestVisitInfo(self.dataId["visit"])) return science, template, sources def _check_converted(self, image, kernelExpected=True): @@ -1611,7 +1623,7 @@ def test_convert_outputs_to_future(self): self.assertIsNone(results.matchedTemplate.getDetector()) legacyDifference = results.difference.image.array.copy() - task.convert_outputs_to_future(results, self.dataId) + task.convert_outputs_to_future(results, self.exposureRecord) self._check_converted(results.difference) self._check_converted(results.matchedTemplate) @@ -1630,7 +1642,7 @@ def test_convert_outputs_to_future_score(self): image_type="future") results = task.run(template, science, sources) - task.convert_outputs_to_future(results, self.dataId) + task.convert_outputs_to_future(results, self.exposureRecord) self._check_converted(results.scoreExposure) self.assertIsInstance(results.psfMatchingKernel, afwMath.LinearCombinationKernel) @@ -1648,7 +1660,7 @@ def test_convert_outputs_to_future_simplified(self): results = task.run(template.clone(), science.clone(), inputPsfMatchingKernel=alResults.psfMatchingKernel) - task.convert_outputs_to_future(results, self.dataId) + task.convert_outputs_to_future(results, self.exposureRecord) self._check_converted(results.difference) self._check_converted(results.matchedTemplate) @@ -1663,7 +1675,7 @@ def test_convert_outputs_to_future_missing_output(self): results = task.run(template, science, sources) del results.matchedTemplate - task.convert_outputs_to_future(results, self.dataId) + task.convert_outputs_to_future(results, self.exposureRecord) self._check_converted(results.difference) self.assertFalse(hasattr(results, "matchedTemplate")) @@ -1704,8 +1716,7 @@ def _make_future_template(self, template, templateInfo=None): """ template = template.clone() template.setDetector(self.detector) - image = DifferenceImage.from_legacy(template, unit=u.nJy, instrument=self.dataId["instrument"], - visit=self.dataId["visit"]) + image = DifferenceImage.from_legacy(template, unit=u.nJy, exposure_record=self.exposureRecord) if templateInfo is not None: image.templates = templateInfo return image @@ -1744,7 +1755,7 @@ def test_convert_outputs_to_future_background(self): self.assertNotEqual(list(legacyBackground.getParameters()), [0.0]*len( legacyBackground.getParameters())) - task.convert_outputs_to_future(results, self.dataId) + task.convert_outputs_to_future(results, self.exposureRecord) self._check_background(results.difference, legacyBackground) self.assertEqual(results.difference.backgrounds.subtracted.field.bounds.bbox.to_legacy(), @@ -1763,7 +1774,7 @@ def test_convert_outputs_to_future_background_score(self): results = task.run(template, science, sources) legacyBackground = results.backgroundModel - task.convert_outputs_to_future(results, self.dataId) + task.convert_outputs_to_future(results, self.exposureRecord) self._check_background(results.scoreExposure, legacyBackground) bbox = results.scoreExposure.backgrounds.subtracted.field.bounds.bbox.to_legacy() @@ -1782,7 +1793,7 @@ def test_convert_outputs_to_future_background_not_subtracted(self): image_type="future") results = task.run(template.clone(), science.clone(), sources) - task.convert_outputs_to_future(results, self.dataId) + task.convert_outputs_to_future(results, self.exposureRecord) self.assertEqual(len(getattr(results, output).backgrounds), 0) @@ -1798,7 +1809,7 @@ def test_convert_outputs_to_future_background_convolve_science(self): self.assertEqual(task.metadata["convolvedExposure"], "Science") legacyBackground = results.backgroundModel - task.convert_outputs_to_future(results, self.dataId) + task.convert_outputs_to_future(results, self.exposureRecord) self._check_background(results.difference, legacyBackground) @@ -1812,7 +1823,7 @@ def test_convert_outputs_to_future_templates(self): results = task.run(template, science, sources) templateInfo = [self._makeTemplateInfo(patch=41), self._makeTemplateInfo(patch=42)] - task.convert_outputs_to_future(results, self.dataId, templateInfo=templateInfo) + task.convert_outputs_to_future(results, self.exposureRecord, templateInfo=templateInfo) for image in (results.difference, results.matchedTemplate): self.assertEqual([(t.tract, t.patch) for t in image.templates], @@ -1827,7 +1838,7 @@ def test_convert_outputs_to_future_no_templates(self): image_type="future") results = task.run(template, science, sources) - task.convert_outputs_to_future(results, self.dataId) + task.convert_outputs_to_future(results, self.exposureRecord) with self.assertRaises(AttributeError): results.difference.templates diff --git a/tests/test_utils.py b/tests/test_utils.py index 34a54b7d9..2d04fc5a2 100644 --- a/tests/test_utils.py +++ b/tests/test_utils.py @@ -23,13 +23,65 @@ import requests from unittest import mock -from lsst.ip.diffim.utils import populate_sattle_visit_cache +import lsst.afw.cameraGeom.testUtils +import lsst.daf.butler +from lsst.images import obs_info_from_legacy +from lsst.images.tests import get_dp2_exposure_record +from lsst.ip.diffim.utils import populate_sattle_visit_cache, record_from_obs_info import lsst.utils.tests from test_detectAndMeasure import makeVisitInfo, MockResponse from utils import makeTestImage +class ExposureRecordTest(lsst.utils.tests.TestCase): + """Tests of rebuilding an exposure record from an observation info.""" + + def setUp(self): + self.universe = lsst.daf.butler.DimensionUniverse() + self.record = get_dp2_exposure_record(self.universe) + # lsst.images needs per-amplifier raw geometry and a field angle + # transform, which the trivial test detectors do not have. + self.detector = list(lsst.afw.cameraGeom.testUtils.CameraWrapper().camera)[0] + self.obsInfo = self._makeObsInfo(self.record) + + def _makeObsInfo(self, record): + """Return the observation info that a converted image built from + ``record`` would carry. + """ + # The conversion requires the visit info and the record to agree + # on the exposure id. + return obs_info_from_legacy(makeVisitInfo(id=record.id), record, self.detector, + detector_exposure_id=12345) + + def _rebuild(self): + """Return the record rebuilt from the observation info, with additional + identifiers added. + """ + return record_from_obs_info(self.obsInfo, self.record.instrument, self.record.id, self.universe) + + def test_fields(self): + """Every field the conversion reads comes back with the value the + original record had. + """ + record = self._rebuild() + for field in ("instrument", "id", "obs_id", "group", "day_obs", "physical_filter", + "exposure_time", "seq_num", "seq_start", "seq_end", "can_see_sky", "timespan"): + self.assertEqual(getattr(record, field), getattr(self.record, field), msg=field) + self.assertFloatsAlmostEqual(record.azimuth, self.record.azimuth, rtol=1e-14) + self.assertFloatsAlmostEqual(record.zenith_angle, self.record.zenith_angle, rtol=1e-14) + + def test_conversion_unchanged(self): + """Converting with the rebuilt record gives the same observation + info as converting with the original one. + """ + self.assertEqual(self._makeObsInfo(self._rebuild()), self.obsInfo) + + def test_missing_obs_info(self): + with self.assertRaises(ValueError): + record_from_obs_info(None, self.record.instrument, self.record.id, self.universe) + + class UtilsTest(lsst.utils.tests.TestCase): def test_populate_sattle(self): diff --git a/tests/utils.py b/tests/utils.py index 47f5fce8b..eaa4500a0 100644 --- a/tests/utils.py +++ b/tests/utils.py @@ -24,6 +24,7 @@ __all__ = ["DipoleTestImage"] +import astropy.time import numpy as np import lsst.geom as geom import lsst.afw.detection as afwDet @@ -31,9 +32,11 @@ import lsst.afw.detection as afwDetection import lsst.afw.geom as afwGeom import lsst.afw.image as afwImage +from lsst.afw.coord import Observatory +import lsst.daf.base as dafBase import lsst.afw.math as afwMath import lsst.afw.table as afwTable -from lsst.daf.butler import DataCoordinate, DimensionUniverse +from lsst.daf.butler import DataCoordinate, DimensionUniverse, Timespan import lsst.meas.algorithms as measAlg import lsst.meas.base as measBase from lsst.meas.algorithms.testUtils import plantSources @@ -1152,6 +1155,78 @@ def getAveragePosition(self): return geom.Point2D(-10000, -10000) +def makeTestVisitInfo(visit=12345): + """Return a visit info for the test images. + + Converting an exposure to an `lsst.images` type reads the observation + metadata it still carries from its visit info, so an exposure that will + be converted needs one. The date matches the timespan of + `makeTestExposureRecord`. + + Parameters + ---------- + visit : `int`, optional + Id of the visit to record. + + Returns + ------- + visitInfo : `lsst.afw.image.VisitInfo` + The visit info. + """ + return afwImage.VisitInfo( + id=visit, + exposureTime=30.0, + date=dafBase.DateTime(60815.000174, dafBase.DateTime.MJD, dafBase.DateTime.TAI), + boresightRaDec=geom.SpherePoint(23.1, 73.2, geom.degrees), + boresightAzAlt=geom.SpherePoint(120.0, 60.0, geom.degrees), + boresightRotAngle=73.2*geom.degrees, + rotType=afwImage.RotType.SKY, + observatory=Observatory(11.1*geom.degrees, 22.2*geom.degrees, 0.333), + ) + + +def makeTestExposureRecord(universe, instrument="testCam", visit=12345): + """Return an ``exposure`` dimension record for the test images. + + The record stands in for the one a task reads from the butler, and holds + the observation metadata that converting an image to an `lsst.images` + type needs. + + Parameters + ---------- + universe : `lsst.daf.butler.DimensionUniverse` + Dimension universe that defines the record's schema. + instrument : `str`, optional + Name of the instrument to record. + visit : `int`, optional + Id of the visit to record. The conversion takes the visit id of an + image from the id of its exposure record. + + Returns + ------- + record : `lsst.daf.butler.DimensionRecord` + The ``exposure`` record. + """ + begin = astropy.time.Time(60815.0, format="mjd", scale="tai") + exposureTime = 30.0 + return universe["exposure"].RecordClass( + instrument=instrument, + id=visit, + obs_id=f"test_{visit}", + group=str(visit), + day_obs=20250520, + physical_filter="test filter", + exposure_time=exposureTime, + seq_num=1, + seq_start=1, + seq_end=1, + can_see_sky=True, + azimuth=120.0, + zenith_angle=30.0, + timespan=Timespan(begin=begin, end=begin + astropy.time.TimeDelta(exposureTime, format="sec")), + ) + + def generate_data_id(*, tract: int = 9813, patch: int = 42, From 40c986fb9a9685fee1c9a3b32f90d4b4f70552f2 Mon Sep 17 00:00:00 2001 From: Ian Sullivan Date: Sun, 27 Sep 2026 23:39:49 +0200 Subject: [PATCH 07/17] Warn when a coadd for the template has unmapped mask planes Only relevant with the `lsst.images` conversion, and the original error can be raised if `raiseOnUndefinedMaskMap` is set. --- python/lsst/ip/diffim/getTemplate.py | 48 +++++++++++++++++++++-- tests/test_getTemplate.py | 57 +++++++++++++++++++++++++++- 2 files changed, 100 insertions(+), 5 deletions(-) diff --git a/python/lsst/ip/diffim/getTemplate.py b/python/lsst/ip/diffim/getTemplate.py index f21691efe..8ff7c9873 100644 --- a/python/lsst/ip/diffim/getTemplate.py +++ b/python/lsst/ip/diffim/getTemplate.py @@ -37,7 +37,12 @@ import lsst.pex.exceptions import lsst.pipe.base as pipeBase -from lsst.images import DifferenceImage, DifferenceImageTemplateInfo, get_legacy_template_mask_planes +from lsst.images import ( + DifferenceImage, + DifferenceImageTemplateInfo, + get_legacy_optional_mask_planes, + get_legacy_template_mask_planes, +) from lsst.skymap import BaseSkyMap from lsst.ip.diffim.dcrModel import DcrModel from lsst.ip.diffim.utils import evaluateMeanPsfFwhm, getPsfFwhm, record_from_obs_info @@ -204,6 +209,14 @@ class GetTemplateConfig( optional=False, default="legacy", ) + raiseOnUndefinedMaskMap = pexConfig.Field( + dtype=bool, + default=False, + doc="Raise if the template has pixels set in a mask plane that" + " `lsst.images` does not map? If False, log a warning and drop that" + " plane from the output instead." + " Only used with output_image_type='future'.", + ) def setDefaults(self): # Use a smaller cache: per SeparableKernel.computeCache, this should @@ -322,8 +335,9 @@ def convert_outputs_to_future(self, template, coadd_refs, detector, exposureReco ------ RuntimeError Raised if the template has pixels set in a mask plane that - `lsst.images.get_legacy_template_mask_planes` does not map, or if - ``coadd_refs`` spans more than one skymap. + `lsst.images.get_legacy_template_mask_planes` does not map and + ``config.raiseOnUndefinedMaskMap`` is True, or if ``coadd_refs`` + spans more than one skymap. """ # The template PSF is a Gaussian, so rebuild the CoaddPsf from the # template inputs to record which coadds contributed. @@ -335,11 +349,14 @@ def convert_outputs_to_future(self, template, coadd_refs, detector, exposureReco f" got {sorted(skymaps)}.") butler_info = {(ref.dataId["tract"], ref.dataId["patch"]): (ref.id, ref.run) for ref in coadd_refs} + planeMap = get_legacy_template_mask_planes() + if not self.config.raiseOnUndefinedMaskMap: + self._clearUnmappedMaskPlanes(template.mask, planeMap) convertedTemplate = DifferenceImage.from_legacy( template, exposure_record=exposureRecord, unit=astropy.units.nJy, - plane_map=get_legacy_template_mask_planes(), + plane_map=planeMap, ) del template convertedTemplate.templates = DifferenceImageTemplateInfo.from_legacy_psf( @@ -351,6 +368,29 @@ def convert_outputs_to_future(self, template, coadd_refs, detector, exposureReco ) return convertedTemplate + def _clearUnmappedMaskPlanes(self, mask, planeMap): + """Clear the pixels of every mask plane that conversion would reject. + + Parameters + ---------- + mask : `lsst.afw.image.Mask` + Template mask to modify in place. + planeMap : `dict` [`str`, `lsst.images.MaskPlane`] + Plane map the template will be converted with. The optional + planes of `lsst.images.get_legacy_optional_mask_planes` are kept + as well, since conversion adds them when they have pixels set. + """ + keep = planeMap.keys() | get_legacy_optional_mask_planes().keys() + for name in mask.getMaskPlaneDict(): + if name in keep: + continue + bitmask = mask.getPlaneBitMask(name) + nSet = np.count_nonzero(mask.array & bitmask) + if nSet: + self.log.warning("Dropping mask plane %s from the template: %d pixels are set, but" + " the template plane map does not include it.", name, nSet) + mask.array &= ~bitmask + def getExposures(self, coaddExposureHandles, bbox, skymap, wcs): """Return a data structure containing the coadds that overlap the specified bbox projected onto the sky, and a corresponding data diff --git a/tests/test_getTemplate.py b/tests/test_getTemplate.py index 731069d5e..9e3357ac3 100644 --- a/tests/test_getTemplate.py +++ b/tests/test_getTemplate.py @@ -647,11 +647,16 @@ def _configureAndRunTask(doScaleVariance, varianceScale=1.): factor*scaleFactor, rtol=1e-5) self.assertImagesAlmostEqual(templateLow.variance, templateOn.variance, rtol=1e-5) - def _runLegacyForFuture(self): + def _runLegacyForFuture(self, raiseOnUndefinedMaskMap=True): """Build a template from tract 0 with a task configured for the future output type, and attach the detector that ``runQuantum`` reads from the science image, but do not convert it. + Parameters + ---------- + raiseOnUndefinedMaskMap : `bool`, optional + Value of the task config field of the same name. + Returns ------- result : `lsst.pipe.base.Struct` @@ -662,6 +667,7 @@ def _runLegacyForFuture(self): """ config = lsst.ip.diffim.GetTemplateTask.ConfigClass() config.output_image_type = "future" + config.raiseOnUndefinedMaskMap = raiseOnUndefinedMaskMap task = lsst.ip.diffim.GetTemplateTask(config=config) self.visit = 9876 box = lsst.geom.Box2I(lsst.geom.Point2I(0, 0), lsst.geom.Point2I(180, 180)) @@ -791,6 +797,55 @@ def testConvertOutputsToFutureCoaddMaskPlanes(self): self.assertEqual( np.count_nonzero(legacy.mask.array & legacy.mask.getPlaneBitMask(plane)), 1) + def _addUnmappedMaskPlane(self, template): + """Set one pixel in the template's BRIGHT_OBJECT plane, which the + template plane map does not include, and one in its CLIPPED plane, + which it does. + """ + mask = template.mask + for n, plane in enumerate(("BRIGHT_OBJECT", "CLIPPED")): + mask.addMaskPlane(plane) + mask.array[0, n] |= mask.getPlaneBitMask(plane) + + def testConvertOutputsToFutureUnmappedMaskPlaneRaises(self): + """By default, a template with pixels set in an unmapped mask plane + cannot be converted. + """ + detector = list(lsst.afw.cameraGeom.testUtils.CameraWrapper().camera)[0] + result, _ = self._runLegacyForFuture() + self._addUnmappedMaskPlane(result.template) + with self.assertRaisesRegex(RuntimeError, "BRIGHT_OBJECT"): + self.futureTask.convert_outputs_to_future( + result.template, self._coaddRefs(0), detector=detector, + exposureRecord=self._exposureRecord()) + + def testConvertOutputsToFutureUnmappedMaskPlaneWarns(self): + """With raiseOnUndefinedMaskMap=False, an unmapped mask plane is + dropped with a warning and the mapped planes are kept. + """ + detector = list(lsst.afw.cameraGeom.testUtils.CameraWrapper().camera)[0] + result, _ = self._runLegacyForFuture(raiseOnUndefinedMaskMap=False) + self._addUnmappedMaskPlane(result.template) + with self.assertLogs(self.futureTask.log.name, level="WARNING") as cm: + template = self.futureTask.convert_outputs_to_future( + result.template, self._coaddRefs(0), detector=detector, + exposureRecord=self._exposureRecord()) + self.assertIn("BRIGHT_OBJECT", "\n".join(cm.output)) + self.assertIsInstance(template, lsst.images.DifferenceImage) + self.assertNotIn("BRIGHT_OBJECT", template.mask.schema.names) + self.assertEqual(np.count_nonzero(template.mask.get("CLIPPED")), 1) + + def testConvertOutputsToFutureLosesProvenance(self): + """The coadd inputs that `run` attaches are not carried by + `lsst.images.DifferenceImage`. + + This pins limitation (2) of ``output_image_type``. + """ + legacy, _ = self._runLegacyForFuture() + self.assertTrue(legacy.template.getInfo().hasCoaddInputs()) + result, _, _ = self._runFuture() + self.assertFalse(result.template.to_legacy().getInfo().hasCoaddInputs()) + def testConvertOutputsToFutureTemplates(self): """The coadds that went into the template are recorded on the converted image. From 3b1e1b59f6b83ef8866b4cd23f0e30dd4b5ede11 Mon Sep 17 00:00:00 2001 From: Ian Sullivan Date: Mon, 28 Sep 2026 14:07:19 +0200 Subject: [PATCH 08/17] Fix stale docs for the lsst.images conversion --- python/lsst/ip/diffim/detectAndMeasure.py | 20 +++++++++------- python/lsst/ip/diffim/getTemplate.py | 17 ++++++++----- python/lsst/ip/diffim/subtractImages.py | 29 +++++++++++------------ python/lsst/ip/diffim/utils.py | 14 +++++------ tests/test_detectAndMeasure.py | 9 +++---- tests/test_getTemplate.py | 17 ++++++------- tests/test_subtractTask.py | 7 +++--- tests/utils.py | 5 ++-- 8 files changed, 61 insertions(+), 57 deletions(-) diff --git a/python/lsst/ip/diffim/detectAndMeasure.py b/python/lsst/ip/diffim/detectAndMeasure.py index 91e1f3e90..fc2232788 100644 --- a/python/lsst/ip/diffim/detectAndMeasure.py +++ b/python/lsst/ip/diffim/detectAndMeasure.py @@ -409,11 +409,11 @@ class DetectAndMeasureConfig(pipeBase.PipelineTaskConfig, ) idGenerator = DetectorVisitIdGeneratorConfig.make_field() image_type = pexConfig.ChoiceField[str]( - "Which image type to use for the measured difference image (and, for" - " DetectAndMeasureScoreTask, the measured score image).", + "Image type of the input and measured difference images (and, for" + " DetectAndMeasureScoreTask, the input and measured score images).", allowed={ - "legacy": "Write as a lsst.afw.image.ExposureF.", - "future": "Write as a lsst.images.DifferenceImage.", + "legacy": "Read and write lsst.afw.image.ExposureF.", + "future": "Read and write lsst.images.DifferenceImage.", }, optional=False, default="legacy", @@ -687,12 +687,16 @@ def convert_outputs_to_future(self, results, exposureRecord, detector=None, matc Detector to set on an output image that has none; all of these images are on the science image's pixel grid, so this should be the science image's detector. - matchingKernels : `lsst.afw.math.Kernel` - Kernel used to PSF-match the subtracted template image. - templateInfo : `list` [`lsst.images.DifferenceImageTemplateInfo`] + matchingKernels : `dict` [`str`, \ + `lsst.images.convolution_kernels.ConvolutionKernel`], optional + Kernel used to PSF-match the template, keyed by the name of the + output image it is attached to. An image with no entry, or an + entry of `None`, gets no kernel. + templateInfo : `list` [`lsst.images.DifferenceImageTemplateInfo`], \ + optional Records of the coadds that went into the template, each holding the second moments of that coadd's PSF and the region where it - overlapped the science image. + overlapped the science image. Not attached if `None` or empty. Notes ----- diff --git a/python/lsst/ip/diffim/getTemplate.py b/python/lsst/ip/diffim/getTemplate.py index 8ff7c9873..9758caf06 100644 --- a/python/lsst/ip/diffim/getTemplate.py +++ b/python/lsst/ip/diffim/getTemplate.py @@ -81,7 +81,8 @@ class GetTemplateConnections( ) obs_info = pipeBase.connectionTypes.Input( doc="Observation metadata of the exposure the template is built for." - " Only read with output_image_type='future'.", + " Only read with output_image_type='future', which needs this to set" + " the obs_info component of a lsst.images.VisitImage.", name="{fakesType}calexp.obs_info", storageClass="ObservationInfo", dimensions=("instrument", "visit", "detector"), @@ -311,15 +312,13 @@ def convert_outputs_to_future(self, template, coadd_refs, detector, exposureReco Parameters ---------- template : `lsst.afw.image.Exposure` - Output struct to read and modify in place. Its ``template`` must - have a detector with per-amplifier raw geometry and a field angle - transform, which is why `runQuantum` sets the science image's - detector on it first. + Template to convert. Its detector is replaced with ``detector`` + and its mask is modified in place. coadd_refs : `list` [`lsst.daf.butler.DatasetRef`] References to the coadds that may have gone into the template. May be a superset of the coadds that were used; supplies the dataset id and RUN collection recorded for each coadd that was. - detector : `lsst.afw.image.Detector` + detector : `lsst.afw.cameraGeom.Detector` Detector of the science image that the template was built for. exposureRecord : `lsst.daf.butler.DimensionRecord` The ``exposure`` record of the observation the template was built @@ -1010,6 +1009,12 @@ def __init__(self, *, config=None): class GetDcrTemplateConfig( GetTemplateConfig, pipelineConnections=GetDcrTemplateConnections ): + """Config for GetDcrTemplateTask. + + Only ``output_image_type="legacy"`` is supported: the task always writes + an `lsst.afw.image.ExposureF`. + """ + numSubfilters = pexConfig.Field( doc="Number of subfilters in the DcrCoadd.", dtype=int, diff --git a/python/lsst/ip/diffim/subtractImages.py b/python/lsst/ip/diffim/subtractImages.py index d559fcd90..418ccd218 100644 --- a/python/lsst/ip/diffim/subtractImages.py +++ b/python/lsst/ip/diffim/subtractImages.py @@ -324,11 +324,11 @@ class AlardLuptonSubtractBaseConfig(lsst.pex.config.Config): " encountered while calculating the matching kernel." ) image_type = lsst.pex.config.ChoiceField[str]( - doc="Which image type to use for the image outputs of this task " - "(difference and matchedTemplate, or scoreExposure).", + doc="Image type of the input template and the image outputs of this" + " task (difference and matchedTemplate, or scoreExposure).", allowed={ - "legacy": "Write as a lsst.afw.image.ExposureF.", - "future": "Write as a lsst.images.DifferenceImage.", + "legacy": "Read and write lsst.afw.image.ExposureF.", + "future": "Read and write lsst.images.DifferenceImage.", }, optional=False, default="legacy", @@ -390,14 +390,14 @@ class AlardLuptonSubtractTask(lsst.pipe.base.PipelineTask): `True`.""" futureImageOutputs = ("difference", "matchedTemplate") - """Names of the images in the results struct that will be converted, in the - same order as the output connections that carry them(`tuple` [`str`]). - Subclasses with different image outputs override this.""" + """Names of the images in the results struct that are converted in future + mode (`tuple` [`str`]). Subclasses with different image outputs override + this.""" futureBackgroundOutput = "difference" - """Name of the difference image in the results struct that the background - should be attached to (`str`). Subclasses that subtract it from - a different image override this.""" + """Name of the image in the results struct that the background is attached + to (`str`). Subclasses that subtract it from a different image override + this.""" def __init__(self, **kwargs): super().__init__(**kwargs) @@ -496,11 +496,11 @@ def convert_outputs_to_future(self, results, exposureRecord, templateInfo=None): exposureRecord : `lsst.daf.butler.DimensionRecord` The ``exposure`` record of the observation, which supplies the observation metadata recorded on each converted output. - templates : `list` [`lsst.images.DifferenceImageTemplateInfo`], - optional (for unit tests) + templateInfo : `list` [`lsst.images.DifferenceImageTemplateInfo`], \ + optional Record of the coadds that went into the template, taken from the template this task subtracted. Attached to every converted - output. + output unless `None` or empty. Raises ------ @@ -510,8 +510,7 @@ def convert_outputs_to_future(self, results, exposureRecord, templateInfo=None): RuntimeError Raised if an output has pixels set in a mask plane that `lsst.images.get_legacy_difference_image_mask_planes` does not - map, or if an output has no detector and there is no science - image detector to take one from. + map. Notes ----- diff --git a/python/lsst/ip/diffim/utils.py b/python/lsst/ip/diffim/utils.py index de2296083..b37e0b5a5 100644 --- a/python/lsst/ip/diffim/utils.py +++ b/python/lsst/ip/diffim/utils.py @@ -57,11 +57,10 @@ def record_from_obs_info(obs_info, instrument, visit, universe): Observation info of an image that the image to be converted is derived from. instrument : `str` - Name of the instrument, from the pixel frame of the image - ``obs_info`` came from. + Name of the instrument, usually from the quantum data ID. visit : `int` - Id of the visit, from the pixel frame of the image ``obs_info`` - came from. + Id of the visit, usually from the quantum data ID. It is also used + as the id of the exposure. universe : `lsst.daf.butler.DimensionUniverse` Dimension universe that defines the record's schema, usually ``butlerQC.quantum.dataId.universe``. @@ -69,7 +68,9 @@ def record_from_obs_info(obs_info, instrument, visit, universe): Returns ------- record : `lsst.daf.butler.DimensionRecord` - The ``exposure`` record of the observation. + An ``exposure`` record holding the fields that `lsst.images` reads + when it converts a legacy image. Fields it does not read, such as + ``dark_time`` and ``target_name``, are `None`. Raises ------ @@ -83,9 +84,6 @@ def record_from_obs_info(obs_info, instrument, visit, universe): observation info does not carry them in the form the conversion needs: its ``instrument`` is read from the legacy FITS header rather than from the record, and its ``exposure_id`` is the id that the visit info held. - Taking both from the pixel frame of the image that supplied - ``obs_info`` puts the converted image on the same instrument and visit - as the image it came from. """ if obs_info is None: raise ValueError("Cannot rebuild an exposure record without an observation info.") diff --git a/tests/test_detectAndMeasure.py b/tests/test_detectAndMeasure.py index f06c99b36..7b08b72c5 100644 --- a/tests/test_detectAndMeasure.py +++ b/tests/test_detectAndMeasure.py @@ -1645,9 +1645,11 @@ def _makeFutureInput(self, exposure, templateInfo=None, kernel=None): ---------- exposure : `lsst.afw.image.ExposureF` Image from `_make_images`. - templateInfo : `list` [`lsst.images.DifferenceImageTemplateInfo`], optional + templateInfo : `list` [`lsst.images.DifferenceImageTemplateInfo`], + optional Record of the coadds that went into the template. - kernel : `lsst.images.convolution_kernels.ImageBasisConvolutionKernel`, optional + kernel : `lsst.images.convolution_kernels.ImageBasisConvolutionKernel`, + optional PSF matching kernel. Returns @@ -1771,8 +1773,7 @@ def test_run_returns_legacy_types(self): self.assertIsInstance(output.subtractedMeasuredExposure, afwImage.ExposureF) def test_legacy_mode_does_not_convert(self): - """Legacy mode is unchanged, with no conversion.. - """ + """Legacy mode leaves the outputs unconverted.""" science, sources, matchedTemplate, difference = self._make_images() task = self._setup_task(doSubtractBackground=True, image_type="legacy") output = task.run(science, matchedTemplate, difference, sources) diff --git a/tests/test_getTemplate.py b/tests/test_getTemplate.py index 9e3357ac3..a60f1afa9 100644 --- a/tests/test_getTemplate.py +++ b/tests/test_getTemplate.py @@ -48,9 +48,7 @@ import lsst.skymap import lsst.utils.tests -from utils import makeTestExposureRecord - -from utils import generate_data_id +from utils import generate_data_id, makeTestExposureRecord # Change this to True, `setup display_ds9`, and open ds9 (or use another afw # display backend) to show the tract/patch layouts on the image. @@ -649,8 +647,9 @@ def _configureAndRunTask(doScaleVariance, varianceScale=1.): def _runLegacyForFuture(self, raiseOnUndefinedMaskMap=True): """Build a template from tract 0 with a task configured for the - future output type, and attach the detector that ``runQuantum`` reads - from the science image, but do not convert it. + future output type, without converting it. + + The task is kept as ``self.futureTask``. Parameters ---------- @@ -660,7 +659,7 @@ def _runLegacyForFuture(self, raiseOnUndefinedMaskMap=True): Returns ------- result : `lsst.pipe.base.Struct` - The legacy output struct, with a detector on the template. + The legacy output struct. box : `lsst.geom.Box2I` The bounding box the template was requested on, before the task grew it by the template border. @@ -808,8 +807,8 @@ def _addUnmappedMaskPlane(self, template): mask.array[0, n] |= mask.getPlaneBitMask(plane) def testConvertOutputsToFutureUnmappedMaskPlaneRaises(self): - """By default, a template with pixels set in an unmapped mask plane - cannot be converted. + """With raiseOnUndefinedMaskMap=True, a template with pixels set in + an unmapped mask plane cannot be converted. """ detector = list(lsst.afw.cameraGeom.testUtils.CameraWrapper().camera)[0] result, _ = self._runLegacyForFuture() @@ -838,8 +837,6 @@ def testConvertOutputsToFutureUnmappedMaskPlaneWarns(self): def testConvertOutputsToFutureLosesProvenance(self): """The coadd inputs that `run` attaches are not carried by `lsst.images.DifferenceImage`. - - This pins limitation (2) of ``output_image_type``. """ legacy, _ = self._runLegacyForFuture() self.assertTrue(legacy.template.getInfo().hasCoaddInputs()) diff --git a/tests/test_subtractTask.py b/tests/test_subtractTask.py index ab9262df9..43a9d7909 100644 --- a/tests/test_subtractTask.py +++ b/tests/test_subtractTask.py @@ -1538,8 +1538,8 @@ def _check_converted(self, image, kernelExpected=True): """ self.assertIsInstance(image, DifferenceImage) self.assertEqual(image.unit, u.nJy) - # `makeTestImage` does not attach a VisitInfo, so the data id supplied - # to the conversion is the only source of these. + # The instrument and visit come from the exposure record passed to + # the conversion. frame = image.sky_projection.pixel_frame self.assertEqual(frame.instrument, self.dataId["instrument"]) self.assertEqual(frame.visit, self.dataId["visit"]) @@ -1706,7 +1706,8 @@ def _make_future_template(self, template, templateInfo=None): ---------- template : `lsst.afw.image.ExposureF` Template from `_make_images`, which has no detector. - templateInfo : `list` [`lsst.images.DifferenceImageTemplateInfo`], optional + templateInfo : `list` [`lsst.images.DifferenceImageTemplateInfo`], + optional Record of the coadds that went into it, left off if `None`. Returns diff --git a/tests/utils.py b/tests/utils.py index eaa4500a0..bfc28ff99 100644 --- a/tests/utils.py +++ b/tests/utils.py @@ -1188,9 +1188,8 @@ def makeTestVisitInfo(visit=12345): def makeTestExposureRecord(universe, instrument="testCam", visit=12345): """Return an ``exposure`` dimension record for the test images. - The record stands in for the one a task reads from the butler, and holds - the observation metadata that converting an image to an `lsst.images` - type needs. + The record holds the observation metadata that converting an image to an + `lsst.images` type needs. Parameters ---------- From 462c852fdccdbde51ba090d79f7c2036b7e525b8 Mon Sep 17 00:00:00 2001 From: Ian Sullivan Date: Mon, 28 Sep 2026 13:43:38 +0200 Subject: [PATCH 09/17] Raise an error if required input components are missing when using lsst.images format images --- python/lsst/ip/diffim/detectAndMeasure.py | 16 +++++---- python/lsst/ip/diffim/subtractImages.py | 7 ++-- python/lsst/ip/diffim/utils.py | 40 ++++++++++++++++++++++- tests/test_detectAndMeasure.py | 20 ++++++++++++ tests/test_subtractTask.py | 16 ++++++++- 5 files changed, 87 insertions(+), 12 deletions(-) diff --git a/python/lsst/ip/diffim/detectAndMeasure.py b/python/lsst/ip/diffim/detectAndMeasure.py index fc2232788..be0aee2e8 100644 --- a/python/lsst/ip/diffim/detectAndMeasure.py +++ b/python/lsst/ip/diffim/detectAndMeasure.py @@ -35,7 +35,8 @@ from lsst.images import DifferenceImage from lsst.images.fields import field_from_legacy_background from lsst.ip.diffim.utils import (evaluateMaskFraction, computeDifferenceImageMetrics, - populate_sattle_visit_cache, record_from_obs_info) + populate_sattle_visit_cache, record_from_obs_info, + get_difference_image_provenance) from lsst.meas.algorithms import SkyObjectsTask, SourceDetectionTask, SetPrimaryFlagsTask, MaskStreaksTask from lsst.meas.algorithms import FindGlintTrailsTask, FindCosmicRaysConfig, findCosmicRays from lsst.meas.base import ForcedMeasurementTask, ApplyApCorrTask, DetectorVisitIdGeneratorConfig @@ -619,19 +620,21 @@ def runQuantum(self, butlerQC: pipeBase.QuantumContext, matchingKernels = {} if self.config.image_type == "future": difference = inputs["difference"] + dataId = butlerQC.quantum.dataId # Read the `TemplateInfo` and the observation metadata from the # difference while it is still in the `DifferenceImage` format. - templateInfo = difference.templates - dataId = butlerQC.quantum.dataId + templateInfo = get_difference_image_provenance(difference, "templates", "difference", dataId) exposureRecord = record_from_obs_info(difference.obs_info, dataId["instrument"], dataId["visit"], dataId.universe) # Assign the kernel from the input difference image to the # corresponding output image - matchingKernels["subtractedMeasuredExposure"] = difference.kernel + matchingKernels["subtractedMeasuredExposure"] = get_difference_image_provenance( + difference, "kernel", "difference", dataId) inputs["difference"] = difference.to_legacy() if "scoreExposure" in inputs: scoreExposure = inputs["scoreExposure"] - matchingKernels["scoreMeasuredExposure"] = scoreExposure.kernel + matchingKernels["scoreMeasuredExposure"] = get_difference_image_provenance( + scoreExposure, "kernel", "scoreExposure", dataId) inputs["scoreExposure"] = scoreExposure.to_legacy() # Specify the fields that `annotate` and `convert_outputs_to_future` @@ -690,8 +693,7 @@ def convert_outputs_to_future(self, results, exposureRecord, detector=None, matc matchingKernels : `dict` [`str`, \ `lsst.images.convolution_kernels.ConvolutionKernel`], optional Kernel used to PSF-match the template, keyed by the name of the - output image it is attached to. An image with no entry, or an - entry of `None`, gets no kernel. + output image it is attached to. If `None`, no kernel is attached. templateInfo : `list` [`lsst.images.DifferenceImageTemplateInfo`], \ optional Records of the coadds that went into the template, each holding diff --git a/python/lsst/ip/diffim/subtractImages.py b/python/lsst/ip/diffim/subtractImages.py index 418ccd218..6c3816900 100644 --- a/python/lsst/ip/diffim/subtractImages.py +++ b/python/lsst/ip/diffim/subtractImages.py @@ -34,7 +34,8 @@ from lsst.images.fields import ChebyshevField from lsst.ip.diffim.utils import (evaluateMeanPsfFwhm, getPsfFwhm, computeDifferenceImageMetrics, computePSFNoiseEquivalentArea, - checkMask, record_from_obs_info, setSourceFootprints) + checkMask, record_from_obs_info, setSourceFootprints, + get_difference_image_provenance) from lsst.meas.algorithms import ScaleVarianceTask, ScienceSourceSelectorTask import lsst.pex.config import lsst.pipe.base @@ -461,10 +462,10 @@ def runQuantum(self, butlerQC, inputRefs, outputRefs): exposureRecord = None if self.config.image_type == "future": template = inputs["template"] + dataId = butlerQC.quantum.dataId # Read the `TemplateInfo` and the observation metadata from the # template while it is still in the `DifferenceImage` format. - templateInfo = template.templates - dataId = butlerQC.quantum.dataId + templateInfo = get_difference_image_provenance(template, "templates", "template", dataId) exposureRecord = record_from_obs_info(template.obs_info, dataId["instrument"], dataId["visit"], dataId.universe) inputs["template"] = template.to_legacy() diff --git a/python/lsst/ip/diffim/utils.py b/python/lsst/ip/diffim/utils.py index b37e0b5a5..9b0940749 100644 --- a/python/lsst/ip/diffim/utils.py +++ b/python/lsst/ip/diffim/utils.py @@ -24,7 +24,7 @@ __all__ = ["evaluateMeanPsfFwhm", "getPsfFwhm", "getKernelCenterDisplacement", "computeDifferenceImageMetrics", "checkMask", "setSourceFootprints", - "record_from_obs_info", + "record_from_obs_info", "get_difference_image_provenance", ] import itertools @@ -110,6 +110,44 @@ def record_from_obs_info(obs_info, instrument, visit, universe): ) +def get_difference_image_provenance(image, attribute, connection, dataId): + """Get the PSF matching kernel or the template records of an input + difference image. + + Parameters + ---------- + image : `lsst.images.DifferenceImage` + Input image read by a task in future mode. + attribute : {"kernel", "templates"} + Name of the attribute to get. + connection : `str` + Name of the input connection the image was read from. + dataId : `lsst.daf.butler.DataCoordinate` + Data ID of the quantum, for the error message. + + Returns + ------- + value : `lsst.images.convolution_kernels.ConvolutionKernel` or \ + `list` [`lsst.images.DifferenceImageTemplateInfo`] + Value of the attribute. + + Raises + ------ + lsst.pipe.base.InvalidQuantumError + Raised if the image does not have the attribute, which happens when + the input was written as a legacy Exposure and converted on read. + """ + try: + return getattr(image, attribute) + except AttributeError as e: + what = "PSF matching kernel" if attribute == "kernel" else "template coadd records" + raise lsst.pipe.base.InvalidQuantumError( + f"Input {connection!r} for {dataId} has no {what}, which is required for image_type='future'." + " This is likely due to reading a legacy Exposure, which lacks this component. To fix," + " reprocess this image starting from `getTemplate`." + ) from e + + def getKernelCenterDisplacement(kernel, x, y, image=None): """Calculate the PSF matching kernel peak offset from the nominal position. diff --git a/tests/test_detectAndMeasure.py b/tests/test_detectAndMeasure.py index 7b08b72c5..567fb22ae 100644 --- a/tests/test_detectAndMeasure.py +++ b/tests/test_detectAndMeasure.py @@ -1889,6 +1889,26 @@ def test_run_quantum(self): self.assertEqual(out.templates, templateInfo) self.assertEqual(out.kernel, kernel) + def test_run_quantum_without_provenance(self): + """An input difference image with no kernel or template records + raises an error that names what is missing, and writes nothing. + """ + templateInfo = [self._makeTemplateInfo()] + kernel = self._makeKernel() + for missing, templates, inputKernel in (("PSF matching kernel", templateInfo, None), + ("template coadd records", None, kernel)): + with self.subTest(missing=missing): + science, sources, matchedTemplate, difference = self._make_images() + task = self._setup_task(image_type="future") + butlerQC = _RecordingQuantumContext(self.dataId) + refs = _FakeRefs(science=science, matchedTemplate=matchedTemplate, + difference=self._makeFutureInput(difference, templates, inputKernel), + kernelSources=sources) + + with self.assertRaisesRegex(lsst.pipe.base.InvalidQuantumError, f"'difference'.*{missing}"): + task.runQuantum(butlerQC, refs, _FakeRefs()) + self.assertIsNone(butlerQC.put_values) + def test_run_quantum_score_task_provenance(self): """Each output takes the kernel of the image it is measured from, which for the score task are two different kernels. diff --git a/tests/test_subtractTask.py b/tests/test_subtractTask.py index 43a9d7909..8a2c7b107 100644 --- a/tests/test_subtractTask.py +++ b/tests/test_subtractTask.py @@ -35,7 +35,7 @@ from lsst.images import Box, DifferenceImage, DifferenceImageTemplateInfo, Polygon from lsst.ip.diffim import subtractImages, InsufficientKernelSourcesError from lsst.pex.config import FieldValidationError -from lsst.pipe.base import NoWorkFound +from lsst.pipe.base import InvalidQuantumError, NoWorkFound import lsst.utils.tests import numpy as np from lsst.ip.diffim.utils import (computeRobustStatistics, computePSFNoiseEquivalentArea, @@ -1873,6 +1873,20 @@ def test_run_quantum_puts_future_types(self): self.assertEqual(butlerQC.put_values.difference.templates, templateInfo) self.assertEqual(butlerQC.put_values.matchedTemplate.templates, templateInfo) + def test_run_quantum_template_without_provenance(self): + """A template with no record of its coadds raises an error that + names what is missing, and writes nothing. + """ + science, template, sources = self._make_images() + task = self._setup_subtraction(subtractImages.AlardLuptonSubtractTask, image_type="future") + butlerQC = _RecordingQuantumContext(self.dataId) + refs = _FakeRefs(template=self._make_future_template(template), science=science.clone(), + sources=sources) + + with self.assertRaisesRegex(InvalidQuantumError, "'template'.*template coadd records"): + task.runQuantum(butlerQC, refs, _FakeRefs()) + self.assertIsNone(butlerQC.put_values) + def test_run_quantum_simplified_puts_future_types(self): """``runQuantum`` converts the outputs of `SimplifiedSubtractTask`, which loads an existing kernel rather than fitting one. From 8eca35b4c57d663c77308a4ce6392ae79bb2a70e Mon Sep 17 00:00:00 2001 From: Ian Sullivan Date: Mon, 28 Sep 2026 21:38:37 +0200 Subject: [PATCH 10/17] Set the photometric_scaling of the DifferenceImage from the science image in future mode --- python/lsst/ip/diffim/detectAndMeasure.py | 15 +++- python/lsst/ip/diffim/subtractImages.py | 21 +++-- tests/test_detectAndMeasure.py | 39 +++++++-- tests/test_subtractTask.py | 99 ++++++++++++++++++++--- 4 files changed, 149 insertions(+), 25 deletions(-) diff --git a/python/lsst/ip/diffim/detectAndMeasure.py b/python/lsst/ip/diffim/detectAndMeasure.py index be0aee2e8..512caad17 100644 --- a/python/lsst/ip/diffim/detectAndMeasure.py +++ b/python/lsst/ip/diffim/detectAndMeasure.py @@ -618,6 +618,7 @@ def runQuantum(self, butlerQC: pipeBase.QuantumContext, templateInfo = None exposureRecord = None matchingKernels = {} + photometricScaling = None if self.config.image_type == "future": difference = inputs["difference"] dataId = butlerQC.quantum.dataId @@ -630,6 +631,7 @@ def runQuantum(self, butlerQC: pipeBase.QuantumContext, # corresponding output image matchingKernels["subtractedMeasuredExposure"] = get_difference_image_provenance( difference, "kernel", "difference", dataId) + photometricScaling = difference.photometric_scaling inputs["difference"] = difference.to_legacy() if "scoreExposure" in inputs: scoreExposure = inputs["scoreExposure"] @@ -662,17 +664,19 @@ def runQuantum(self, butlerQC: pipeBase.QuantumContext, if self.config.image_type == "future": self.convert_outputs_to_future(measurementResults, exposureRecord, detector=detector, matchingKernels=matchingKernels, - templateInfo=templateInfo) + templateInfo=templateInfo, + photometricScaling=photometricScaling) butlerQC.put(measurementResults, outputRefs) raise error from e if self.config.image_type == "future": self.convert_outputs_to_future(measurementResults, exposureRecord, detector=detector, matchingKernels=matchingKernels, - templateInfo=templateInfo) + templateInfo=templateInfo, + photometricScaling=photometricScaling) butlerQC.put(measurementResults, outputRefs) def convert_outputs_to_future(self, results, exposureRecord, detector=None, matchingKernels=None, - templateInfo=None): + templateInfo=None, photometricScaling=None): """Convert the output images in a result struct to `lsst.images` types. This replaces ``results.subtractedMeasuredExposure`` and @@ -699,6 +703,9 @@ def convert_outputs_to_future(self, results, exposureRecord, detector=None, matc Records of the coadds that went into the template, each holding the second moments of that coadd's PSF and the region where it overlapped the science image. Not attached if `None` or empty. + photometricScaling : `lsst.images.fields.BaseField`, optional + Photometric scaling of the science image, taken from the input + difference image. Not attached if `None`. Notes ----- @@ -723,6 +730,8 @@ def convert_outputs_to_future(self, results, exposureRecord, detector=None, matc image.templates = templateInfo if matchingKernels is not None: image.kernel = matchingKernels[name] + if photometricScaling is not None: + image.photometric_scaling = photometricScaling if background is not None and len(background) > 0: # The same background is subtracted from the difference image # and the score image, so it applies to both. diff --git a/python/lsst/ip/diffim/subtractImages.py b/python/lsst/ip/diffim/subtractImages.py index 6c3816900..dc1797cf3 100644 --- a/python/lsst/ip/diffim/subtractImages.py +++ b/python/lsst/ip/diffim/subtractImages.py @@ -83,7 +83,9 @@ class SubtractInputConnections(lsst.pipe.base.PipelineTaskConnections, name="{fakesType}{coaddName}Diff_templateExp" ) science = connectionTypes.Input( - doc="Input science exposure to subtract from.", + doc="Input science exposure to subtract from." + " With image_type='future' this is read as a lsst.images.VisitImage, so that its" + " photometric scaling can be passed on to the outputs.", dimensions=("instrument", "visit", "detector"), storageClass="ExposureF", name="{fakesType}calexp" @@ -110,6 +112,7 @@ def __init__(self, *, config=None): del self.visitSummary if config.image_type == "future": self.template = dataclasses.replace(self.template, storageClass="DifferenceImage") + self.science = dataclasses.replace(self.science, storageClass="VisitImage") class SubtractImageOutputConnections(lsst.pipe.base.PipelineTaskConnections, @@ -460,6 +463,7 @@ def runQuantum(self, butlerQC, inputRefs, outputRefs): inputs = butlerQC.get(inputRefs) templateInfo = None exposureRecord = None + photometricScaling = None if self.config.image_type == "future": template = inputs["template"] dataId = butlerQC.quantum.dataId @@ -469,6 +473,8 @@ def runQuantum(self, butlerQC, inputRefs, outputRefs): exposureRecord = record_from_obs_info(template.obs_info, dataId["instrument"], dataId["visit"], dataId.universe) inputs["template"] = template.to_legacy() + photometricScaling = inputs["science"].photometric_scaling + inputs["science"] = inputs["science"].to_legacy() try: results = self.run(**inputs) @@ -478,11 +484,13 @@ def runQuantum(self, butlerQC, inputRefs, outputRefs): raise error from e if self.config.image_type == "future": - self.convert_outputs_to_future(results, exposureRecord, templateInfo=templateInfo) + self.convert_outputs_to_future(results, exposureRecord, templateInfo=templateInfo, + photometricScaling=photometricScaling) butlerQC.put(results, outputRefs) - def convert_outputs_to_future(self, results, exposureRecord, templateInfo=None): + def convert_outputs_to_future(self, results, exposureRecord, templateInfo=None, + photometricScaling=None): """Convert the image outputs to `lsst.images` types. Each image named in `futureImageOutputs` that is present on @@ -500,8 +508,9 @@ def convert_outputs_to_future(self, results, exposureRecord, templateInfo=None): templateInfo : `list` [`lsst.images.DifferenceImageTemplateInfo`], \ optional Record of the coadds that went into the template, taken from the - template this task subtracted. Attached to every converted - output unless `None` or empty. + template this task subtracted. + photometricScaling : `lsst.images.fields.BaseField`, optional + Photometric scaling of the science image. Raises ------ @@ -548,6 +557,8 @@ def convert_outputs_to_future(self, results, exposureRecord, templateInfo=None): image.kernel = kernel if templateInfo: image.templates = templateInfo + if photometricScaling is not None: + image.photometric_scaling = photometricScaling if background is not None and name == self.futureBackgroundOutput: image.backgrounds.add( "subtracted", background, diff --git a/tests/test_detectAndMeasure.py b/tests/test_detectAndMeasure.py index 567fb22ae..a6d688134 100644 --- a/tests/test_detectAndMeasure.py +++ b/tests/test_detectAndMeasure.py @@ -36,6 +36,7 @@ import lsst.geom from lsst.images import Box, DifferenceImage, DifferenceImageTemplateInfo, Polygon from lsst.images.convolution_kernels import ImageBasisConvolutionKernel +from lsst.images.fields import ChebyshevField import lsst.images.psfs import lsst.images.serialization from lsst.ip.diffim import detectAndMeasure, subtractImages @@ -1638,7 +1639,24 @@ def _makeKernel(scale=1.0): legacy.setSpatialParameters([[scale, 0.0, 0.0], [0.5, 0.0, 0.0]]) return ImageBasisConvolutionKernel.from_legacy(legacy) - def _makeFutureInput(self, exposure, templateInfo=None, kernel=None): + @staticmethod + def _makePhotometricScaling(exposure): + """Return a constant photometric scaling for an input image. + + Parameters + ---------- + exposure : `lsst.afw.image.ExposureF` + Image from `_make_images`. + + Returns + ------- + photometricScaling : `lsst.images.fields.ChebyshevField` + The scaling. + """ + return ChebyshevField(Box.from_legacy(exposure.getBBox()), np.array([[0.25]]), + unit=astropy.units.nJy/astropy.units.electron) + + def _makeFutureInput(self, exposure, templateInfo=None, kernel=None, photometricScaling=None): """Convert a legacy image to a `lsst.images.DifferenceImage`. Parameters @@ -1651,6 +1669,8 @@ def _makeFutureInput(self, exposure, templateInfo=None, kernel=None): kernel : `lsst.images.convolution_kernels.ImageBasisConvolutionKernel`, optional PSF matching kernel. + photometricScaling : `lsst.images.fields.BaseField`, optional + Photometric scaling of the science image. Returns ------- @@ -1663,6 +1683,8 @@ def _makeFutureInput(self, exposure, templateInfo=None, kernel=None): image.templates = templateInfo if kernel is not None: image.kernel = kernel + if photometricScaling is not None: + image.photometric_scaling = photometricScaling return image def _make_images(self, withScore=False): @@ -1738,8 +1760,7 @@ def _check_converted(self, image, legacyBbox): self.assertEqual(image.bbox.y.stop, legacyBbox.getMaxY() + 1) self.assertEqual(image.unit, astropy.units.nJy) self.assertIsNotNone(image.detector) - # Populating this needs the visit summary PhotoCalib, which this task - # does not have as an input, so it is deliberately left unset. + # `convert_outputs_to_future` was not given a photometric scaling. self.assertIsNone(image.photometric_scaling) def test_connections(self): @@ -1871,8 +1892,10 @@ def test_run_quantum(self): ("future", DifferenceImage)): with self.subTest(imageType=imageType): science, sources, matchedTemplate, difference = self._make_images() + photometricScaling = self._makePhotometricScaling(difference) if imageType == "future": - difference = self._makeFutureInput(difference, templateInfo, kernel) + difference = self._makeFutureInput(difference, templateInfo, kernel, + photometricScaling) task = self._setup_task(doSubtractBackground=True, image_type=imageType) butlerQC = _RecordingQuantumContext(self.dataId) task.runQuantum(butlerQC, @@ -1888,6 +1911,7 @@ def test_run_quantum(self): # survives, rather than being worked out again. self.assertEqual(out.templates, templateInfo) self.assertEqual(out.kernel, kernel) + self.assertIs(out.photometric_scaling, photometricScaling) def test_run_quantum_without_provenance(self): """An input difference image with no kernel or template records @@ -1919,13 +1943,14 @@ def test_run_quantum_score_task_provenance(self): differenceKernel = self._makeKernel(scale=1.0) scoreKernel = self._makeKernel(scale=3.0) self.assertNotEqual(differenceKernel, scoreKernel) + photometricScaling = self._makePhotometricScaling(difference) task = self._setup_task(image_type="future") butlerQC = _RecordingQuantumContext(self.dataId) task.runQuantum(butlerQC, _FakeRefs(science=science, matchedTemplate=matchedTemplate, difference=self._makeFutureInput(difference, templateInfo, - differenceKernel), + differenceKernel, photometricScaling), scoreExposure=self._makeFutureInput(score, templateInfo, scoreKernel), kernelSources=sources), _FakeRefs()) @@ -1936,6 +1961,10 @@ def test_run_quantum_score_task_provenance(self): self.assertEqual(scoreOut.templates, templateInfo) self.assertEqual(out.kernel, differenceKernel) self.assertEqual(scoreOut.kernel, scoreKernel) + # Both outputs are in the science image's units, so both take the + # scaling of the input difference image. + self.assertIs(out.photometric_scaling, photometricScaling) + self.assertIs(scoreOut.photometric_scaling, photometricScaling) def test_connections_future_reads_difference_image(self): """The inputs whose provenance is copied ask for the type that diff --git a/tests/test_subtractTask.py b/tests/test_subtractTask.py index 8a2c7b107..84d04ee1c 100644 --- a/tests/test_subtractTask.py +++ b/tests/test_subtractTask.py @@ -32,7 +32,8 @@ import lsst.geom import lsst.meas.algorithms as measAlg from lsst.daf.butler import DataCoordinate, DimensionUniverse -from lsst.images import Box, DifferenceImage, DifferenceImageTemplateInfo, Polygon +from lsst.images import Box, DifferenceImage, DifferenceImageTemplateInfo, Polygon, VisitImage +from lsst.images.fields import ChebyshevField from lsst.ip.diffim import subtractImages, InsufficientKernelSourcesError from lsst.pex.config import FieldValidationError from lsst.pipe.base import InvalidQuantumError, NoWorkFound @@ -1526,7 +1527,7 @@ def _make_images(self): exposure.info.setVisitInfo(makeTestVisitInfo(self.dataId["visit"])) return science, template, sources - def _check_converted(self, image, kernelExpected=True): + def _check_converted(self, image, kernelExpected=True, photometricScaling=None): """Check that one converted image is a fully-formed DifferenceImage. Parameters @@ -1535,6 +1536,9 @@ def _check_converted(self, image, kernelExpected=True): The converted output to check. kernelExpected : `bool`, optional Should the PSF matching kernel be attached? + photometricScaling : `lsst.images.fields.BaseField`, optional + Photometric scaling the image should carry, or `None` if it should + have none. """ self.assertIsInstance(image, DifferenceImage) self.assertEqual(image.unit, u.nJy) @@ -1550,8 +1554,10 @@ def _check_converted(self, image, kernelExpected=True): # An unset kernel raises, rather than returning None. with self.assertRaises(AttributeError): image.kernel - # Populating this needs the visitSummary PhotoCalib; it is deferred. - self.assertIsNone(image.photometric_scaling) + if photometricScaling is None: + self.assertIsNone(image.photometric_scaling) + else: + self.assertIs(image.photometric_scaling, photometricScaling) def test_connections_legacy(self): """Test that the default configuration leaves the output connections @@ -1572,6 +1578,7 @@ def test_connections_legacy(self): for name in connections.outputs: if name in ("difference", "matchedTemplate", "scoreExposure"): self.assertEqual(getattr(connections, name).storageClass, "ExposureF") + self.assertEqual(connections.science.storageClass, "ExposureF") def test_connections_future(self): """Test that ``image_type="future"`` swaps the storage class of @@ -1594,6 +1601,9 @@ def test_connections_future(self): connections = connectionsClass(config=config) for name in expected: self.assertEqual(getattr(connections, name).storageClass, "DifferenceImage") + # The science image is read in the new type to get its + # photometric scaling. + self.assertEqual(connections.science.storageClass, "VisitImage") # The kernel has no converter to or from the new types, so it # must remain a separate output in both modes. if "psfMatchingKernel" in connections.outputs: @@ -1665,6 +1675,24 @@ def test_convert_outputs_to_future_simplified(self): self._check_converted(results.difference) self._check_converted(results.matchedTemplate) + def test_convert_outputs_to_future_photometric_scaling(self): + """Test that every converted output carries the photometric scaling + of the science image. + """ + science, template, sources = self._make_images() + photometricScaling = self._makePhotometricScaling(science) + for taskClass, names in ((subtractImages.AlardLuptonSubtractTask, ("difference", "matchedTemplate")), + (subtractImages.AlardLuptonPreconvolveSubtractTask, ("scoreExposure",))): + with self.subTest(taskClass=taskClass.__name__): + task = self._setup_subtraction(taskClass, image_type="future") + results = task.run(template.clone(), science.clone(), sources) + + task.convert_outputs_to_future(results, self.exposureRecord, + photometricScaling=photometricScaling) + + for name in names: + self._check_converted(getattr(results, name), photometricScaling=photometricScaling) + def test_convert_outputs_to_future_missing_output(self): """Test that conversion skips outputs that are absent from the result struct, instead of failing. @@ -1699,6 +1727,45 @@ def _makeTemplateInfo(tract=9813, patch=42): dataset_run="a/template/run", bounds=Polygon.from_box(Box.factory[0:50, 0:50]), psf_shape_xx=4.0, psf_shape_yy=4.0, psf_shape_xy=0.0, psf_shape_flag=False) + @staticmethod + def _makePhotometricScaling(science): + """Return a photometric scaling that is easy to tell apart from the + one converted from the science image's own calibration. + + Parameters + ---------- + science : `lsst.afw.image.ExposureF` + Science image from `_make_images`. + + Returns + ------- + photometricScaling : `lsst.images.fields.ChebyshevField` + A constant scaling over the science image. + """ + return ChebyshevField(Box.from_legacy(science.getBBox()), np.array([[0.25]]), + unit=u.nJy/u.electron) + + def _make_future_science(self, science, photometricScaling=None): + """Convert a legacy science image to the type the ``future`` mode + reads. + + Parameters + ---------- + science : `lsst.afw.image.ExposureF` + Science image from `_make_images`. + photometricScaling : `lsst.images.fields.BaseField`, optional + Photometric scaling to attach in place of the converted one. + + Returns + ------- + science : `lsst.images.VisitImage` + The converted science image. + """ + image = VisitImage.from_legacy(science.clone(), unit=u.nJy, exposure_record=self.exposureRecord) + if photometricScaling is not None: + image.photometric_scaling = photometricScaling + return image + def _make_future_template(self, template, templateInfo=None): """Convert a legacy template to the type the ``future`` mode reads. @@ -1850,6 +1917,7 @@ def test_run_quantum_puts_future_types(self): """ science, template, sources = self._make_images() templateInfo = [self._makeTemplateInfo()] + photometricScaling = self._makePhotometricScaling(science) for imageType, expectedType in (("legacy", afwImage.ExposureF), ("future", DifferenceImage)): with self.subTest(imageType=imageType): @@ -1857,11 +1925,13 @@ def test_run_quantum_puts_future_types(self): image_type=imageType) if imageType == "future": inputTemplate = self._make_future_template(template, templateInfo) + inputScience = self._make_future_science(science, photometricScaling) else: inputTemplate = template.clone() + inputScience = science.clone() butlerQC = _RecordingQuantumContext(self.dataId) task.runQuantum(butlerQC, - _FakeRefs(template=inputTemplate, science=science.clone(), + _FakeRefs(template=inputTemplate, science=inputScience, sources=sources), _FakeRefs()) self.assertIsInstance(butlerQC.put_values.difference, expectedType) @@ -1870,8 +1940,10 @@ def test_run_quantum_puts_future_types(self): self.assertIsInstance(butlerQC.put_values.psfMatchingKernel, afwMath.LinearCombinationKernel) if imageType == "future": - self.assertEqual(butlerQC.put_values.difference.templates, templateInfo) - self.assertEqual(butlerQC.put_values.matchedTemplate.templates, templateInfo) + for name in ("difference", "matchedTemplate"): + image = getattr(butlerQC.put_values, name) + self.assertEqual(image.templates, templateInfo) + self.assertIs(image.photometric_scaling, photometricScaling) def test_run_quantum_template_without_provenance(self): """A template with no record of its coadds raises an error that @@ -1880,8 +1952,8 @@ def test_run_quantum_template_without_provenance(self): science, template, sources = self._make_images() task = self._setup_subtraction(subtractImages.AlardLuptonSubtractTask, image_type="future") butlerQC = _RecordingQuantumContext(self.dataId) - refs = _FakeRefs(template=self._make_future_template(template), science=science.clone(), - sources=sources) + refs = _FakeRefs(template=self._make_future_template(template), + science=self._make_future_science(science), sources=sources) with self.assertRaisesRegex(InvalidQuantumError, "'template'.*template coadd records"): task.runQuantum(butlerQC, refs, _FakeRefs()) @@ -1893,6 +1965,7 @@ def test_run_quantum_simplified_puts_future_types(self): """ science, template, sources = self._make_images() templateInfo = [self._makeTemplateInfo()] + photometricScaling = self._makePhotometricScaling(science) alTask = self._setup_subtraction(subtractImages.AlardLuptonSubtractTask) alResults = alTask.run(template.clone(), science.clone(), sources) for useExistingKernel in (True, False): @@ -1901,15 +1974,17 @@ def test_run_quantum_simplified_puts_future_types(self): useExistingKernel=useExistingKernel, image_type="future") refs = {"template": self._make_future_template(template, templateInfo), - "science": science.clone()} + "science": self._make_future_science(science, photometricScaling)} if useExistingKernel: refs["inputPsfMatchingKernel"] = alResults.psfMatchingKernel butlerQC = _RecordingQuantumContext(self.dataId) task.runQuantum(butlerQC, _FakeRefs(**refs), _FakeRefs()) - self._check_converted(butlerQC.put_values.difference) - self._check_converted(butlerQC.put_values.matchedTemplate) + self._check_converted(butlerQC.put_values.difference, + photometricScaling=photometricScaling) + self._check_converted(butlerQC.put_values.matchedTemplate, + photometricScaling=photometricScaling) self.assertEqual(butlerQC.put_values.difference.templates, templateInfo) From a106c7e0ac4170568bed198527de8fc1f840a031 Mon Sep 17 00:00:00 2001 From: Ian Sullivan Date: Tue, 29 Sep 2026 00:24:53 +0200 Subject: [PATCH 11/17] Add a deprecation warning to doSubtractBackground in subtractImages --- python/lsst/ip/diffim/subtractImages.py | 1 + tests/test_computeSpatiallySampledMetrics.py | 1 - tests/test_subtractTask.py | 1 - 3 files changed, 1 insertion(+), 2 deletions(-) diff --git a/python/lsst/ip/diffim/subtractImages.py b/python/lsst/ip/diffim/subtractImages.py index dc1797cf3..7bf0140b6 100644 --- a/python/lsst/ip/diffim/subtractImages.py +++ b/python/lsst/ip/diffim/subtractImages.py @@ -244,6 +244,7 @@ class AlardLuptonSubtractBaseConfig(lsst.pex.config.Config): "It is generally better to instead subtract the background in detectAndMeasure.", dtype=bool, default=False, + deprecated="Subtract the background in detectAndMeasure instead. Will be removed after v31.", ) doApplyExternalCalibrations = lsst.pex.config.Field( doc=( diff --git a/tests/test_computeSpatiallySampledMetrics.py b/tests/test_computeSpatiallySampledMetrics.py index 4f12d1d13..8b1abf786 100644 --- a/tests/test_computeSpatiallySampledMetrics.py +++ b/tests/test_computeSpatiallySampledMetrics.py @@ -40,7 +40,6 @@ def setUp(self): template, _ = makeTestImage(psfSize=2.0, noiseLevel=noiseLevel, noiseSeed=7, templateBorderSize=20, doApplyCalibration=True) config = lsst.ip.diffim.AlardLuptonSubtractTask.ConfigClass() - config.doSubtractBackground = False config.sourceSelector.signalToNoise.fluxField = "truth_instFlux" config.sourceSelector.signalToNoise.errField = "truth_instFluxErr" subtraction = lsst.ip.diffim.AlardLuptonSubtractTask(config=config).run(template, science, sources) diff --git a/tests/test_subtractTask.py b/tests/test_subtractTask.py index 84d04ee1c..167b6dabb 100644 --- a/tests/test_subtractTask.py +++ b/tests/test_subtractTask.py @@ -70,7 +70,6 @@ def _setup_subtraction(self, fluxField="truth_instFlux", errField="truth_instFlu The configured Task to use for detection and measurement. """ config = self.subtractTask.ConfigClass() - config.doSubtractBackground = False config.restrictKernelEdgeSources = False config.sourceSelector.signalToNoise.fluxField = fluxField config.sourceSelector.signalToNoise.errField = errField From ab86ad04cd4fdc159cacc4b8e3429cb55239cca2 Mon Sep 17 00:00:00 2001 From: Ian Sullivan Date: Tue, 29 Sep 2026 00:30:39 +0200 Subject: [PATCH 12/17] Reuse input DifferenceImage in detectAndMeasure The input image already carries the kernel, template records, observation metadata, PSF, and photometric scaling, so only the pixels, mask, and background need to be written back. The subtractImages conversion no longer attaches the background fit with the kernel, which is unused, so the detectAndMeasure background is the only one recorded. --- python/lsst/ip/diffim/detectAndMeasure.py | 108 +++++++--------------- python/lsst/ip/diffim/subtractImages.py | 20 +--- tests/test_detectAndMeasure.py | 99 ++++++++++++-------- tests/test_subtractTask.py | 92 ------------------ 4 files changed, 99 insertions(+), 220 deletions(-) diff --git a/python/lsst/ip/diffim/detectAndMeasure.py b/python/lsst/ip/diffim/detectAndMeasure.py index 512caad17..4f063028d 100644 --- a/python/lsst/ip/diffim/detectAndMeasure.py +++ b/python/lsst/ip/diffim/detectAndMeasure.py @@ -24,7 +24,6 @@ import requests import os -import astropy.units import lsst.afw.detection as afwDetection import lsst.afw.image as afwImage @@ -32,10 +31,10 @@ import lsst.afw.table as afwTable import lsst.daf.base as dafBase import lsst.geom -from lsst.images import DifferenceImage +from lsst.images import Mask, get_legacy_difference_image_mask_planes from lsst.images.fields import field_from_legacy_background from lsst.ip.diffim.utils import (evaluateMaskFraction, computeDifferenceImageMetrics, - populate_sattle_visit_cache, record_from_obs_info, + populate_sattle_visit_cache, get_difference_image_provenance) from lsst.meas.algorithms import SkyObjectsTask, SourceDetectionTask, SetPrimaryFlagsTask, MaskStreaksTask from lsst.meas.algorithms import FindGlintTrailsTask, FindCosmicRaysConfig, findCosmicRays @@ -612,32 +611,22 @@ def runQuantum(self, butlerQC: pipeBase.QuantumContext, inputs = butlerQC.get(inputRefs) idGenerator = self.config.idGenerator.apply(butlerQC.quantum.dataId) idFactory = idGenerator.make_table_id_factory() - # Assign the science image Detector to the output DifferenceImage - detector = inputs["science"].getDetector() - templateInfo = None - exposureRecord = None - matchingKernels = {} - photometricScaling = None + # Input image that each output is measured from, keyed by output name. + inputImages = {} if self.config.image_type == "future": - difference = inputs["difference"] dataId = butlerQC.quantum.dataId - # Read the `TemplateInfo` and the observation metadata from the - # difference while it is still in the `DifferenceImage` format. - templateInfo = get_difference_image_provenance(difference, "templates", "difference", dataId) - exposureRecord = record_from_obs_info(difference.obs_info, dataId["instrument"], - dataId["visit"], dataId.universe) - # Assign the kernel from the input difference image to the - # corresponding output image - matchingKernels["subtractedMeasuredExposure"] = get_difference_image_provenance( - difference, "kernel", "difference", dataId) - photometricScaling = difference.photometric_scaling - inputs["difference"] = difference.to_legacy() - if "scoreExposure" in inputs: - scoreExposure = inputs["scoreExposure"] - matchingKernels["scoreMeasuredExposure"] = get_difference_image_provenance( - scoreExposure, "kernel", "scoreExposure", dataId) - inputs["scoreExposure"] = scoreExposure.to_legacy() + for output, connection in (("subtractedMeasuredExposure", "difference"), + ("scoreMeasuredExposure", "scoreExposure")): + if connection not in inputs: + continue + image = inputs[connection] + # The output is written back into this image, so it must + # already carry the kernel and the template records. + get_difference_image_provenance(image, "kernel", connection, dataId) + get_difference_image_provenance(image, "templates", connection, dataId) + inputImages[output] = image + inputs[connection] = image.to_legacy() # Specify the fields that `annotate` and `convert_outputs_to_future` # need below, to ensure they exist, even as None. @@ -662,50 +651,30 @@ def runQuantum(self, butlerQC: pipeBase.QuantumContext, log=self.log ) if self.config.image_type == "future": - self.convert_outputs_to_future(measurementResults, exposureRecord, - detector=detector, matchingKernels=matchingKernels, - templateInfo=templateInfo, - photometricScaling=photometricScaling) + self.convert_outputs_to_future(measurementResults, inputImages) butlerQC.put(measurementResults, outputRefs) raise error from e if self.config.image_type == "future": - self.convert_outputs_to_future(measurementResults, exposureRecord, - detector=detector, matchingKernels=matchingKernels, - templateInfo=templateInfo, - photometricScaling=photometricScaling) + self.convert_outputs_to_future(measurementResults, inputImages) butlerQC.put(measurementResults, outputRefs) - def convert_outputs_to_future(self, results, exposureRecord, detector=None, matchingKernels=None, - templateInfo=None, photometricScaling=None): - """Convert the output images in a result struct to `lsst.images` types. + def convert_outputs_to_future(self, results, inputImages): + """Write the output images in a result struct back into the + `lsst.images` images they were measured from. - This replaces ``results.subtractedMeasuredExposure`` and - ``results.scoreMeasuredExposure`` with `lsst.images.DifferenceImage` - instances, skipping either one that is absent or `None`. + The image, variance, and mask planes of each input image are replaced + by those of the corresponding legacy output, and the background + subtracted by this task is attached. The kernel, template records, + and observation metadata of the input are kept. Outputs that are + absent or `None` are skipped. Parameters ---------- results : `lsst.pipe.base.Struct` Output struct to read and modify in place. - exposureRecord : `lsst.daf.butler.DimensionRecord` - The ``exposure`` record of the observation, which supplies the - observation metadata recorded on each converted output. - detector : `lsst.afw.cameraGeom.Detector`, optional - Detector to set on an output image that has none; all of these - images are on the science image's pixel grid, so this should be - the science image's detector. - matchingKernels : `dict` [`str`, \ - `lsst.images.convolution_kernels.ConvolutionKernel`], optional - Kernel used to PSF-match the template, keyed by the name of the - output image it is attached to. If `None`, no kernel is attached. - templateInfo : `list` [`lsst.images.DifferenceImageTemplateInfo`], \ - optional - Records of the coadds that went into the template, each holding - the second moments of that coadd's PSF and the region where it - overlapped the science image. Not attached if `None` or empty. - photometricScaling : `lsst.images.fields.BaseField`, optional - Photometric scaling of the science image, taken from the input - difference image. Not attached if `None`. + inputImages : `dict` [`str`, `lsst.images.DifferenceImage`] + Input image each output was measured from, keyed by the name of + the output. Each is modified in place. Notes ----- @@ -715,23 +684,16 @@ def convert_outputs_to_future(self, results, exposureRecord, detector=None, matc for any pixels in the image. """ background = getattr(results, "differenceBackground", None) - for name in ("subtractedMeasuredExposure", "scoreMeasuredExposure"): + for name, image in inputImages.items(): exposure = getattr(results, name, None) if exposure is None: continue - if exposure.getDetector() is None and detector is not None: - exposure.setDetector(detector) - image = DifferenceImage.from_legacy( - exposure, - unit=astropy.units.nJy, - exposure_record=exposureRecord, - ) - if templateInfo: - image.templates = templateInfo - if matchingKernels is not None: - image.kernel = matchingKernels[name] - if photometricScaling is not None: - image.photometric_scaling = photometricScaling + # These are no-ops when the legacy output is a view of the input, + # as returned by `~lsst.images.DifferenceImage.to_legacy`. + image.image.array[...] = exposure.image.array + image.variance.array[...] = exposure.variance.array + image.mask = Mask.from_legacy(exposure.mask, get_legacy_difference_image_mask_planes(), + sky_projection=image.sky_projection) if background is not None and len(background) > 0: # The same background is subtracted from the difference image # and the score image, so it applies to both. diff --git a/python/lsst/ip/diffim/subtractImages.py b/python/lsst/ip/diffim/subtractImages.py index 7bf0140b6..546dee839 100644 --- a/python/lsst/ip/diffim/subtractImages.py +++ b/python/lsst/ip/diffim/subtractImages.py @@ -31,7 +31,6 @@ import lsst.geom from lsst.images import DifferenceImage from lsst.images.convolution_kernels import ImageBasisConvolutionKernel -from lsst.images.fields import ChebyshevField from lsst.ip.diffim.utils import (evaluateMeanPsfFwhm, getPsfFwhm, computeDifferenceImageMetrics, computePSFNoiseEquivalentArea, checkMask, record_from_obs_info, setSourceFootprints, @@ -399,11 +398,6 @@ class AlardLuptonSubtractTask(lsst.pipe.base.PipelineTask): mode (`tuple` [`str`]). Subclasses with different image outputs override this.""" - futureBackgroundOutput = "difference" - """Name of the image in the results struct that the background is attached - to (`str`). Subclasses that subtract it from a different image override - this.""" - def __init__(self, **kwargs): super().__init__(**kwargs) self.makeSubtask("decorrelate") @@ -496,8 +490,7 @@ def convert_outputs_to_future(self, results, exposureRecord, templateInfo=None, Each image named in `futureImageOutputs` that is present on ``results`` is replaced by an `lsst.images.DifferenceImage`, with the - PSF matching kernel attached to it, and the differential background - attached to the one named by `futureBackgroundOutput`. + PSF matching kernel attached to it. Parameters ---------- @@ -538,9 +531,6 @@ def convert_outputs_to_future(self, results, exposureRecord, templateInfo=None, except (TypeError, ValueError, AttributeError) as e: raise TypeError("The supplied PSF matching kernel cannot be attached to an" f" lsst.images.DifferenceImage: {e}") from e - background = None - if self.config.doSubtractBackground: - background = ChebyshevField.from_legacy_function2(results.backgroundModel, unit=u.nJy) for name in self.futureImageOutputs: exposure = getattr(results, name, None) if exposure is None: @@ -560,13 +550,6 @@ def convert_outputs_to_future(self, results, exposureRecord, templateInfo=None, image.templates = templateInfo if photometricScaling is not None: image.photometric_scaling = photometricScaling - if background is not None and name == self.futureBackgroundOutput: - image.backgrounds.add( - "subtracted", background, - description="Differential background subtracted from this image when solving for" - " the PSF matching kernel.", - is_subtracted=True, - ) setattr(results, name, image) @timeMethod @@ -1366,7 +1349,6 @@ class AlardLuptonPreconvolveSubtractTask(AlardLuptonSubtractTask): # This task just writes the score image instead of the difference image and # matchedTemplate futureImageOutputs = ("scoreExposure",) - futureBackgroundOutput = "scoreExposure" def run(self, template, science, sources, visitSummary=None): """Preconvolve the science image with its own PSF, diff --git a/tests/test_detectAndMeasure.py b/tests/test_detectAndMeasure.py index a6d688134..67fdd2ec5 100644 --- a/tests/test_detectAndMeasure.py +++ b/tests/test_detectAndMeasure.py @@ -34,11 +34,13 @@ import lsst.afw.table as afwTable import lsst.afw.math as afwMath import lsst.geom -from lsst.images import Box, DifferenceImage, DifferenceImageTemplateInfo, Polygon +from lsst.images import (Box, DifferenceImage, DifferenceImageTemplateInfo, Polygon, + get_legacy_difference_image_mask_planes, get_legacy_optional_mask_planes) from lsst.images.convolution_kernels import ImageBasisConvolutionKernel from lsst.images.fields import ChebyshevField import lsst.images.psfs import lsst.images.serialization +from lsst.images.tests import compare_masked_image_to_legacy from lsst.ip.diffim import detectAndMeasure, subtractImages from lsst.afw.table import IdFactory from lsst.afw.cameraGeom.testUtils import CameraWrapper, DetectorWrapper @@ -1743,25 +1745,24 @@ def _setup_task(self, **kwargs): kwargs.setdefault("badSubtractionRatioThreshold", 1.) return self._setup_detection(doSkySources=False, **kwargs) - def _check_converted(self, image, legacyBbox): + def _check_converted(self, image, legacy): """Check the properties common to every converted output image. Parameters ---------- image : `lsst.images.DifferenceImage` The converted image. - legacyBbox : `lsst.geom.Box2I` - Bounding box of the legacy exposure it was converted from. + legacy : `lsst.afw.image.ExposureF` + The legacy output exposure it was converted from. """ self.assertIsInstance(image, DifferenceImage) - self.assertEqual(image.bbox.x.start, legacyBbox.getMinX()) - self.assertEqual(image.bbox.x.stop, legacyBbox.getMaxX() + 1) - self.assertEqual(image.bbox.y.start, legacyBbox.getMinY()) - self.assertEqual(image.bbox.y.stop, legacyBbox.getMaxY() + 1) self.assertEqual(image.unit, astropy.units.nJy) self.assertIsNotNone(image.detector) - # `convert_outputs_to_future` was not given a photometric scaling. + # `_makeFutureInput` was not given a photometric scaling. self.assertIsNone(image.photometric_scaling) + compare_masked_image_to_legacy(image, legacy.maskedImage, + plane_map=get_legacy_difference_image_mask_planes()) + self.assertTrue(np.any(image.mask.get("DETECTED"))) def test_connections(self): """Check that the output storage classes follow ``image_type``, @@ -1806,17 +1807,37 @@ def test_convert_outputs_to_future(self): """Convert the outputs of ``DetectAndMeasureTask``.""" science, sources, matchedTemplate, difference = self._make_images() task = self._setup_task(doSubtractBackground=True, image_type="future") - output = task.run(science, matchedTemplate, difference, sources) - legacyBbox = output.subtractedMeasuredExposure.getBBox() + inputImage = self._makeFutureInput(difference) + output = task.run(science, matchedTemplate, inputImage.to_legacy(), sources) + legacy = output.subtractedMeasuredExposure background = output.differenceBackground - task.convert_outputs_to_future(output, self.exposureRecord) + task.convert_outputs_to_future(output, {"subtractedMeasuredExposure": inputImage}) - self._check_converted(output.subtractedMeasuredExposure, legacyBbox) - self.assertIn("subtracted", output.subtractedMeasuredExposure.backgrounds) + self.assertIs(output.subtractedMeasuredExposure, inputImage) + self._check_converted(inputImage, legacy) + self.assertEqual(list(inputImage.backgrounds), ["subtracted"]) # The separate background dataset is still written in future mode. self.assertIs(output.differenceBackground, background) self.assertIsInstance(output.differenceBackground, afwMath.BackgroundList) + def test_convert_outputs_to_future_copied_pixels(self): + """The output pixels are copied into the input image when the legacy + image that was measured does not share them. + """ + science, sources, matchedTemplate, difference = self._make_images() + task = self._setup_task(doSubtractBackground=True, image_type="future") + inputImage = self._makeFutureInput(difference) + original = inputImage.image.array.copy() + output = task.run(science, matchedTemplate, inputImage.to_legacy(copy=True), sources) + legacy = output.subtractedMeasuredExposure + # Background subtraction changed the pixels of the copy only. + self.assertFalse(np.array_equal(legacy.image.array, original)) + np.testing.assert_array_equal(inputImage.image.array, original) + + task.convert_outputs_to_future(output, {"subtractedMeasuredExposure": inputImage}) + + self._check_converted(inputImage, legacy) + def test_convert_outputs_to_future_source_injection(self): """Convert a difference image from a pipeline that injects sources. @@ -1826,33 +1847,39 @@ def test_convert_outputs_to_future_source_injection(self): planes = ("INJECTED", "INJECTED_CORE", "INJECTED_TEMPLATE", "INJECTED_CORE_TEMPLATE") science, sources, matchedTemplate, difference = self._make_images() task = self._setup_task(image_type="future") - output = task.run(science, matchedTemplate, difference, sources) - mask = output.subtractedMeasuredExposure.mask + inputImage = self._makeFutureInput(difference) + output = task.run(science, matchedTemplate, inputImage.to_legacy(), sources) + legacy = output.subtractedMeasuredExposure for n, plane in enumerate(planes): - mask.addMaskPlane(plane) - mask.array[0, n] |= mask.getPlaneBitMask(plane) + legacy.mask.addMaskPlane(plane) + legacy.mask.array[0, n] |= legacy.mask.getPlaneBitMask(plane) - task.convert_outputs_to_future(output, self.exposureRecord) + task.convert_outputs_to_future(output, {"subtractedMeasuredExposure": inputImage}) image = output.subtractedMeasuredExposure self.assertIsInstance(image, DifferenceImage) - for plane in planes: - with self.subTest(plane=plane): - self.assertEqual(np.count_nonzero(image.mask.get(plane)), 1) + # Name the optional planes in the map, so that dropping one fails + # the comparison instead of skipping it. + planeMap = get_legacy_difference_image_mask_planes() | get_legacy_optional_mask_planes() + compare_masked_image_to_legacy(image, legacy.maskedImage, plane_map=planeMap) def test_convert_outputs_to_future_score(self): """Convert the outputs of ``DetectAndMeasureScoreTask``.""" self.detectionTask = detectAndMeasure.DetectAndMeasureScoreTask science, sources, matchedTemplate, difference, score = self._make_images(withScore=True) task = self._setup_task(doSubtractBackground=True, image_type="future") - output = task.run(science, matchedTemplate, difference, score, sources) - legacyBboxes = {name: getattr(output, name).getBBox() - for name in ("subtractedMeasuredExposure", "scoreMeasuredExposure")} - task.convert_outputs_to_future(output, self.exposureRecord) - - for name, legacyBbox in legacyBboxes.items(): + inputImages = {"subtractedMeasuredExposure": self._makeFutureInput(difference), + "scoreMeasuredExposure": self._makeFutureInput(score)} + output = task.run(science, matchedTemplate, + inputImages["subtractedMeasuredExposure"].to_legacy(), + inputImages["scoreMeasuredExposure"].to_legacy(), sources) + legacyOutputs = {name: getattr(output, name) for name in inputImages} + task.convert_outputs_to_future(output, inputImages) + + for name, legacy in legacyOutputs.items(): with self.subTest(name=name): - self._check_converted(getattr(output, name), legacyBbox) + self.assertIs(getattr(output, name), inputImages[name]) + self._check_converted(getattr(output, name), legacy) # The same background is subtracted from both images. self.assertIn("subtracted", getattr(output, name).backgrounds) @@ -1862,11 +1889,12 @@ def test_convert_outputs_to_future_empty_background(self): """ science, sources, matchedTemplate, difference = self._make_images() task = self._setup_task(doSubtractBackground=False, image_type="future") - output = task.run(science, matchedTemplate, difference, sources) + inputImage = self._makeFutureInput(difference) + output = task.run(science, matchedTemplate, inputImage.to_legacy(), sources) self.assertEqual(len(output.differenceBackground), 0) - legacyBbox = output.subtractedMeasuredExposure.getBBox() - task.convert_outputs_to_future(output, self.exposureRecord) - self._check_converted(output.subtractedMeasuredExposure, legacyBbox) + legacy = output.subtractedMeasuredExposure + task.convert_outputs_to_future(output, {"subtractedMeasuredExposure": inputImage}) + self._check_converted(output.subtractedMeasuredExposure, legacy) self.assertEqual(len(output.subtractedMeasuredExposure.backgrounds), 0) def test_stray_struct_field_is_ignored(self): @@ -1951,7 +1979,8 @@ def test_run_quantum_score_task_provenance(self): _FakeRefs(science=science, matchedTemplate=matchedTemplate, difference=self._makeFutureInput(difference, templateInfo, differenceKernel, photometricScaling), - scoreExposure=self._makeFutureInput(score, templateInfo, scoreKernel), + scoreExposure=self._makeFutureInput(score, templateInfo, scoreKernel, + photometricScaling), kernelSources=sources), _FakeRefs()) @@ -1961,8 +1990,6 @@ def test_run_quantum_score_task_provenance(self): self.assertEqual(scoreOut.templates, templateInfo) self.assertEqual(out.kernel, differenceKernel) self.assertEqual(scoreOut.kernel, scoreKernel) - # Both outputs are in the science image's units, so both take the - # scaling of the input difference image. self.assertIs(out.photometric_scaling, photometricScaling) self.assertIs(scoreOut.photometric_scaling, photometricScaling) diff --git a/tests/test_subtractTask.py b/tests/test_subtractTask.py index 167b6dabb..ccb3221e6 100644 --- a/tests/test_subtractTask.py +++ b/tests/test_subtractTask.py @@ -1788,98 +1788,6 @@ def _make_future_template(self, template, templateInfo=None): image.templates = templateInfo return image - def _check_background(self, image, legacyBackground): - """Check that a converted image carries the background - that was subtracted from it. - - Parameters - ---------- - image : `lsst.images.DifferenceImage` - The converted image to check. - legacyBackground : `lsst.afw.math.Chebyshev1Function2D` - Background model from the result struct. - """ - self.assertEqual(list(image.backgrounds), ["subtracted"]) - field = image.backgrounds.subtracted.field - self.assertEqual(field.unit, u.nJy) - bbox = field.bounds.bbox.to_legacy() - rendered = field.render().array - for y in np.linspace(bbox.getMinY(), bbox.getMaxY(), 5, dtype=int): - for x in np.linspace(bbox.getMinX(), bbox.getMaxX(), 5, dtype=int): - self.assertFloatsAlmostEqual(rendered[y - bbox.getMinY(), x - bbox.getMinX()], - legacyBackground(float(x), float(y)), rtol=1e-12) - - def test_convert_outputs_to_future_background(self): - """The background goes on the difference image, which is - the one it was subtracted from. - """ - science, template, sources = self._make_images() - task = self._setup_subtraction(subtractImages.AlardLuptonSubtractTask, - doSubtractBackground=True, - image_type="future") - results = task.run(template, science, sources) - legacyBackground = results.backgroundModel - self.assertNotEqual(list(legacyBackground.getParameters()), [0.0]*len( - legacyBackground.getParameters())) - - task.convert_outputs_to_future(results, self.exposureRecord) - - self._check_background(results.difference, legacyBackground) - self.assertEqual(results.difference.backgrounds.subtracted.field.bounds.bbox.to_legacy(), - science.getBBox()) - # The background was only taken out of the difference image. - self.assertEqual(len(results.matchedTemplate.backgrounds), 0) - - def test_convert_outputs_to_future_background_score(self): - """The score image gets the background, over the inner region that - the kernel was fit on. - """ - science, template, sources = self._make_images() - task = self._setup_subtraction(subtractImages.AlardLuptonPreconvolveSubtractTask, - doSubtractBackground=True, - image_type="future") - results = task.run(template, science, sources) - legacyBackground = results.backgroundModel - - task.convert_outputs_to_future(results, self.exposureRecord) - - self._check_background(results.scoreExposure, legacyBackground) - bbox = results.scoreExposure.backgrounds.subtracted.field.bounds.bbox.to_legacy() - self.assertTrue(results.scoreExposure.bbox.to_legacy().contains(bbox)) - self.assertNotEqual(bbox, results.scoreExposure.bbox.to_legacy()) - - def test_convert_outputs_to_future_background_not_subtracted(self): - """A background the task fit but did not subtract is left off the - outputs. - """ - science, template, sources = self._make_images() - for taskClass, output in ((subtractImages.AlardLuptonSubtractTask, "difference"), - (subtractImages.AlardLuptonPreconvolveSubtractTask, "scoreExposure")): - with self.subTest(taskClass=taskClass.__name__): - task = self._setup_subtraction(taskClass, doSubtractBackground=False, - image_type="future") - results = task.run(template.clone(), science.clone(), sources) - - task.convert_outputs_to_future(results, self.exposureRecord) - - self.assertEqual(len(getattr(results, output).backgrounds), 0) - - def test_convert_outputs_to_future_background_convolve_science(self): - """In ``convolveScience`` mode the actual background used is inverted. - """ - science, template, sources = self._make_images() - task = self._setup_subtraction(subtractImages.AlardLuptonSubtractTask, - mode="convolveScience", - doSubtractBackground=True, - image_type="future") - results = task.run(template, science, sources) - self.assertEqual(task.metadata["convolvedExposure"], "Science") - legacyBackground = results.backgroundModel - - task.convert_outputs_to_future(results, self.exposureRecord) - - self._check_background(results.difference, legacyBackground) - def test_convert_outputs_to_future_templates(self): """The record of the coadds that went into the template is copied to every converted output. From 168dc6fbb9a535cb77121bb1f56a1b21663437ad Mon Sep 17 00:00:00 2001 From: Ian Sullivan Date: Tue, 29 Sep 2026 00:31:31 +0200 Subject: [PATCH 13/17] Reject doSubtractBackground in future mode --- python/lsst/ip/diffim/subtractImages.py | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/python/lsst/ip/diffim/subtractImages.py b/python/lsst/ip/diffim/subtractImages.py index 546dee839..8d86a2089 100644 --- a/python/lsst/ip/diffim/subtractImages.py +++ b/python/lsst/ip/diffim/subtractImages.py @@ -369,6 +369,15 @@ def setDefaults(self): self.fallbackSourceSelector.signalToNoise.minimum = signalToNoiseMinimum self.fallbackSourceSelector.signalToNoise.maximum = signalToNoiseMaximum + def validate(self): + super().validate() + if self.doSubtractBackground and self.image_type == "future": + raise lsst.pex.config.FieldValidationError( + AlardLuptonSubtractBaseConfig.doSubtractBackground, self, + "The background subtracted when solving the kernel is not recorded on" + " image_type='future' outputs; subtract the background in detectAndMeasure instead." + ) + class AlardLuptonSubtractConfig(AlardLuptonSubtractBaseConfig, lsst.pipe.base.PipelineTaskConfig, pipelineConnections=AlardLuptonSubtractConnections): From e55396b4198bfcf10990855cc3e47de48a454734 Mon Sep 17 00:00:00 2001 From: Ian Sullivan Date: Wed, 7 Oct 2026 17:16:32 -0700 Subject: [PATCH 14/17] Restore the difference mask after background fit The score image variant runs its final detection on the score image, so the first-pass detections used to fit the background stayed in the DETECTED planes of the output difference image. --- python/lsst/ip/diffim/detectAndMeasure.py | 17 +++++++------ tests/test_detectAndMeasure.py | 30 +++++++++++++++++++++++ 2 files changed, 40 insertions(+), 7 deletions(-) diff --git a/python/lsst/ip/diffim/detectAndMeasure.py b/python/lsst/ip/diffim/detectAndMeasure.py index 4f063028d..cf717cc9c 100644 --- a/python/lsst/ip/diffim/detectAndMeasure.py +++ b/python/lsst/ip/diffim/detectAndMeasure.py @@ -876,16 +876,19 @@ def _fitAndSubtractBackground(self, differenceExposure, scoreExposure=None): background=background, clearMask=True, ) - # Use the temporary detection mask for the final background subtraction. - # The detection mask planes will be cleared before the final detection - # step, so it is OK if they get set for differenceExposure. - detectedBit = differenceExposure.mask.getPlaneBitMask(["DETECTED"]) + # Use the temporary detection mask for the final background subtraction, + # then restore the original mask. The score image variant never + # clears the detection planes of the difference image. + mask = differenceExposure.mask + originalMask = mask.array.copy() + detectedBit = mask.getPlaneBitMask(["DETECTED"]) detectedPix = detectionExposure.mask.array & detectedBit > 0 - detectedNegativeBit = differenceExposure.mask.getPlaneBitMask(["DETECTED_NEGATIVE"]) + detectedNegativeBit = mask.getPlaneBitMask(["DETECTED_NEGATIVE"]) detectedNegativePix = detectionExposure.mask.array & detectedNegativeBit > 0 - differenceExposure.mask.array[detectedPix] |= detectedBit - differenceExposure.mask.array[detectedNegativePix] |= detectedNegativeBit + mask.array[detectedPix] |= detectedBit + mask.array[detectedNegativePix] |= detectedNegativeBit background = self.subtractFinalBackground.run(differenceExposure).background + mask.array[...] = originalMask if scoreExposure is not None: # The preconvolution kernel is normalized to 1, so the same # background level applies to the difference and score images. diff --git a/tests/test_detectAndMeasure.py b/tests/test_detectAndMeasure.py index 67fdd2ec5..c56ab2a27 100644 --- a/tests/test_detectAndMeasure.py +++ b/tests/test_detectAndMeasure.py @@ -1319,6 +1319,36 @@ def _run_and_check_detections(positive=True): _run_and_check_detections(positive=True) _run_and_check_detections(positive=False) + def test_background_detection_mask_not_kept(self): + """The detection masks of the difference image should match those of + the score image, without the detections used to fit the background. + """ + kwargs = {"seed": 1, "psfSize": 2.4, "fluxLevel": 500, "kernelSize": 31, + "templateBorderSize": 15} + science, sources = makeTestImage(noiseLevel=1., noiseSeed=6, **kwargs) + matchedTemplate, _ = makeTestImage(noiseLevel=0.25, noiseSeed=7, **kwargs) + difference = science.clone() + # Small-scale structure that the first-pass background fit removes + # but the final fit does not, so the two passes detect different + # pixels. + difference.image.array += 0.3*np.sin(np.arange(difference.getWidth())/20.) + subtractTask = subtractImages.AlardLuptonPreconvolveSubtractTask() + score = subtractTask._convolveExposure(difference, science.psf.getKernel(), + subtractTask.convolutionControl) + detectionTask = self._setup_detection(doSubtractBackground=True, doSkySources=False, + badSubtractionRatioThreshold=1.) + + output = detectionTask.run(science, matchedTemplate, difference, score, sources) + + differenceMask = output.subtractedMeasuredExposure.mask + scoreMask = output.scoreMeasuredExposure.mask + for plane in ("DETECTED", "DETECTED_NEGATIVE"): + with self.subTest(plane=plane): + differencePix = differenceMask.array & differenceMask.getPlaneBitMask(plane) > 0 + scorePix = scoreMask.array & scoreMask.getPlaneBitMask(plane) > 0 + self.assertGreater(scorePix.sum(), 0) + np.testing.assert_array_equal(differencePix, scorePix) + def test_mask_cosmic_rays(self): """Cosmic rays detected on the difference image should propagate to the mask of the returned (measured) exposure. This version creates From ea6fe1b7a27bea9a33837221e466f7ca9794456e Mon Sep 17 00:00:00 2001 From: Ian Sullivan Date: Wed, 7 Oct 2026 17:24:28 -0700 Subject: [PATCH 15/17] Deprecate doWriteBackground in detectAndMeasure In future mode the background is recorded on the output images, so the separate difference_background dataset is no longer written. --- python/lsst/ip/diffim/detectAndMeasure.py | 9 +++++++-- tests/test_detectAndMeasure.py | 21 ++++++++++++++------- 2 files changed, 21 insertions(+), 9 deletions(-) diff --git a/python/lsst/ip/diffim/detectAndMeasure.py b/python/lsst/ip/diffim/detectAndMeasure.py index cf717cc9c..719b80a92 100644 --- a/python/lsst/ip/diffim/detectAndMeasure.py +++ b/python/lsst/ip/diffim/detectAndMeasure.py @@ -184,7 +184,9 @@ def __init__(self, *, config): self.inputs.remove("kernelSources") if not (self.config.writeStreakInfo and self.config.doMaskStreaks): self.outputs.remove("maskedStreaks") - if not (self.config.doSubtractBackground and self.config.doWriteBackground): + # In future mode the background is recorded on the output images. + if (self.config.image_type == "future" + or not (self.config.doSubtractBackground and self.config.doWriteBackground)): self.outputs.remove("differenceBackground") if not (self.config.writeGlintInfo): self.outputs.remove("glintTrailInfo") @@ -221,8 +223,11 @@ class DetectAndMeasureConfig(pipeBase.PipelineTaskConfig, ) doWriteBackground = pexConfig.Field( dtype=bool, - doc="Persist the fitted background model?", + doc="Persist the fitted background model? Ignored when" + " image_type='future', which records it on the output images.", default=False, + deprecated="Use image_type='future', which records the background on the output images." + " Will be removed after v31.", ) doCalculateResidualMetics = pexConfig.Field( dtype=bool, diff --git a/tests/test_detectAndMeasure.py b/tests/test_detectAndMeasure.py index c56ab2a27..a2f26858f 100644 --- a/tests/test_detectAndMeasure.py +++ b/tests/test_detectAndMeasure.py @@ -1796,7 +1796,7 @@ def _check_converted(self, image, legacy): def test_connections(self): """Check that the output storage classes follow ``image_type``, - and that the background output is kept in both modes. + and that the background is a separate output only in legacy mode. """ for configClass, connectionsClass in ( (detectAndMeasure.DetectAndMeasureConfig, detectAndMeasure.DetectAndMeasureConnections), @@ -1808,14 +1808,18 @@ def test_connections(self): config = configClass() config.image_type = imageType config.doSubtractBackground = True - config.doWriteBackground = True + with self.assertWarns(FutureWarning): + config.doWriteBackground = True connections = connectionsClass(config=config) self.assertEqual(connections.subtractedMeasuredExposure.storageClass, storageClass) if connectionsClass is detectAndMeasure.DetectAndMeasureScoreConnections: self.assertEqual(connections.scoreMeasuredExposure.storageClass, storageClass) - # The background is a separate dataset in both modes. - self.assertIn("differenceBackground", connections.outputs) - self.assertEqual(connections.differenceBackground.storageClass, "Background") + if imageType == "legacy": + self.assertIn("differenceBackground", connections.outputs) + self.assertEqual(connections.differenceBackground.storageClass, "Background") + else: + # The background is recorded on the output images. + self.assertNotIn("differenceBackground", connections.outputs) def test_run_returns_legacy_types(self): """``run`` always returns legacy types, even in future mode.""" @@ -1962,8 +1966,11 @@ def test_run_quantum(self): _FakeRefs()) out = butlerQC.put_values.subtractedMeasuredExposure self.assertIsInstance(out, expectedType) - # The background is a separate dataset in both modes. - self.assertIsInstance(butlerQC.put_values.differenceBackground, afwMath.BackgroundList) + if imageType == "legacy": + self.assertIsInstance(butlerQC.put_values.differenceBackground, + afwMath.BackgroundList) + else: + self.assertIsNotNone(out.backgrounds.subtracted) if imageType == "future": # The provenance of the image this was measured from # survives, rather than being worked out again. From 791cb4261bb994db2fde670693a5d968060b7381 Mon Sep 17 00:00:00 2001 From: Ian Sullivan Date: Wed, 7 Oct 2026 17:31:39 -0700 Subject: [PATCH 16/17] Use ConvolutionKernel for kernels in future mode Needs the daf_butler formatter for ConvolutionKernel datasets. --- .../diffim/computeSpatiallySampledMetrics.py | 25 +++++++ python/lsst/ip/diffim/subtractImages.py | 25 +++++-- tests/test_computeSpatiallySampledMetrics.py | 50 ++++++++++++++ tests/test_subtractTask.py | 67 ++++++++++++++----- 4 files changed, 147 insertions(+), 20 deletions(-) diff --git a/python/lsst/ip/diffim/computeSpatiallySampledMetrics.py b/python/lsst/ip/diffim/computeSpatiallySampledMetrics.py index 2f698e514..3006991e6 100644 --- a/python/lsst/ip/diffim/computeSpatiallySampledMetrics.py +++ b/python/lsst/ip/diffim/computeSpatiallySampledMetrics.py @@ -19,6 +19,8 @@ # You should have received a copy of the GNU General Public License # along with this program. If not, see . +import dataclasses + import numpy as np import scipy.signal @@ -81,6 +83,12 @@ class SpatiallySampledMetricsConnections(pipeBase.PipelineTaskConnections, name="{fakesType}{coaddName}Diff_spatiallySampledMetrics", ) + def __init__(self, *, config=None): + super().__init__(config=config) + if config.image_type == "future": + self.psfMatchingKernel = dataclasses.replace(self.psfMatchingKernel, + storageClass="ConvolutionKernel") + class SpatiallySampledMetricsConfig(pipeBase.PipelineTaskConfig, pipelineConnections=SpatiallySampledMetricsConnections): @@ -100,6 +108,16 @@ class SpatiallySampledMetricsConfig(pipeBase.PipelineTaskConfig, target=SkyObjectsTask, doc="Generate QA metric sources", ) + image_type = pexConfig.ChoiceField[str]( + doc="Type of the input PSF matching kernel. The input images are read" + " as lsst.afw.image.ExposureF either way.", + allowed={ + "legacy": "Read lsst.afw.math.Kernel.", + "future": "Read lsst.images.convolution_kernels.ConvolutionKernel.", + }, + optional=False, + default="legacy", + ) def setDefaults(self): self.metricSources.avoidMask = ["NO_DATA", "EDGE"] @@ -216,6 +234,13 @@ def __init__(self, **kwargs): " relative to the science PSF L2 norm. Larger values indicate worse" " PSF matching. Assumes the kernel was solved to convolve the template.") + def runQuantum(self, butlerQC, inputRefs, outputRefs): + inputs = butlerQC.get(inputRefs) + if self.config.image_type == "future": + inputs["psfMatchingKernel"] = inputs["psfMatchingKernel"].to_legacy() + outputs = self.run(**inputs) + butlerQC.put(outputs, outputRefs) + @timeMethod def run(self, science, template, difference, diaSources, psfMatchingKernel): """Calculate difference image metrics on specific locations across the images diff --git a/python/lsst/ip/diffim/subtractImages.py b/python/lsst/ip/diffim/subtractImages.py index 8d86a2089..3ca166dea 100644 --- a/python/lsst/ip/diffim/subtractImages.py +++ b/python/lsst/ip/diffim/subtractImages.py @@ -148,6 +148,8 @@ def __init__(self, *, config=None): self.difference = dataclasses.replace(self.difference, storageClass="DifferenceImage") self.matchedTemplate = dataclasses.replace(self.matchedTemplate, storageClass="DifferenceImage") + self.psfMatchingKernel = dataclasses.replace(self.psfMatchingKernel, + storageClass="ConvolutionKernel") class SubtractScoreOutputConnections(lsst.pipe.base.PipelineTaskConnections, @@ -176,6 +178,8 @@ def __init__(self, *, config=None): super().__init__(config=config) if config.image_type == "future": self.scoreExposure = dataclasses.replace(self.scoreExposure, storageClass="DifferenceImage") + self.psfMatchingKernel = dataclasses.replace(self.psfMatchingKernel, + storageClass="ConvolutionKernel") class AlardLuptonSubtractConnections(SubtractInputConnections, SubtractImageOutputConnections): @@ -196,6 +200,9 @@ def __init__(self, *, config=None): if config.useExistingKernel: del self.psfMatchingKernel del self.kernelSources + if config.image_type == "future": + self.inputPsfMatchingKernel = dataclasses.replace(self.inputPsfMatchingKernel, + storageClass="ConvolutionKernel") else: del self.inputPsfMatchingKernel @@ -328,11 +335,14 @@ class AlardLuptonSubtractBaseConfig(lsst.pex.config.Config): " encountered while calculating the matching kernel." ) image_type = lsst.pex.config.ChoiceField[str]( - doc="Image type of the input template and the image outputs of this" - " task (difference and matchedTemplate, or scoreExposure).", + doc="Image type of the input template, the image outputs of this" + " task (difference and matchedTemplate, or scoreExposure), and the PSF" + " matching kernel.", allowed={ - "legacy": "Read and write lsst.afw.image.ExposureF.", - "future": "Read and write lsst.images.DifferenceImage.", + "legacy": "Read and write lsst.afw.image.ExposureF and" + " lsst.afw.math.Kernel.", + "future": "Read and write lsst.images.DifferenceImage and" + " lsst.images.convolution_kernels.ConvolutionKernel.", }, optional=False, default="legacy", @@ -479,6 +489,8 @@ def runQuantum(self, butlerQC, inputRefs, outputRefs): inputs["template"] = template.to_legacy() photometricScaling = inputs["science"].photometric_scaling inputs["science"] = inputs["science"].to_legacy() + if "inputPsfMatchingKernel" in inputs: + inputs["inputPsfMatchingKernel"] = inputs["inputPsfMatchingKernel"].to_legacy() try: results = self.run(**inputs) @@ -499,7 +511,9 @@ def convert_outputs_to_future(self, results, exposureRecord, templateInfo=None, Each image named in `futureImageOutputs` that is present on ``results`` is replaced by an `lsst.images.DifferenceImage`, with the - PSF matching kernel attached to it. + PSF matching kernel attached to it. ``psfMatchingKernel`` is replaced + by the same kernel as an + `lsst.images.convolution_kernels.ImageBasisConvolutionKernel`. Parameters ---------- @@ -560,6 +574,7 @@ def convert_outputs_to_future(self, results, exposureRecord, templateInfo=None, if photometricScaling is not None: image.photometric_scaling = photometricScaling setattr(results, name, image) + results.psfMatchingKernel = kernel @timeMethod def run(self, template, science, sources, visitSummary=None): diff --git a/tests/test_computeSpatiallySampledMetrics.py b/tests/test_computeSpatiallySampledMetrics.py index 8b1abf786..8bb99c1e3 100644 --- a/tests/test_computeSpatiallySampledMetrics.py +++ b/tests/test_computeSpatiallySampledMetrics.py @@ -19,6 +19,7 @@ # You should have received a copy of the GNU General Public License # along with this program. If not, see . +import types import unittest import numpy as np @@ -27,6 +28,7 @@ import lsst.afw.table import lsst.ip.diffim import lsst.utils.tests +from lsst.images.convolution_kernels import ImageBasisConvolutionKernel from utils import makeTestImage @@ -89,6 +91,54 @@ def testMaskPlaneNotInImage(self): self.assertTrue(np.all(metrics[f"{missing.lower()}_mask_fraction"] == 0)) self.assertTrue(np.all(np.isfinite(metrics[f"{present.lower()}_mask_fraction"]))) + def testConnections(self): + """``image_type`` sets the storage class of the kernel input only.""" + for imageType, kernelStorageClass in (("legacy", "MatchingKernel"), + ("future", "ConvolutionKernel")): + with self.subTest(imageType=imageType): + config = lsst.ip.diffim.SpatiallySampledMetricsTask.ConfigClass() + config.image_type = imageType + connections = config.ConnectionsClass(config=config) + self.assertEqual(connections.psfMatchingKernel.storageClass, kernelStorageClass) + for name in ("science", "template", "difference"): + self.assertEqual(getattr(connections, name).storageClass, "ExposureF") + + def testRunQuantumFuture(self): + """``runQuantum`` converts a future kernel back to the legacy type, + giving the same metrics as a legacy run. + """ + inputs = {"science": self.science, "template": self.template, "difference": self.difference, + "diaSources": self.diaSources} + metrics = {} + for imageType, kernel in ( + ("legacy", self.psfMatchingKernel), + ("future", ImageBasisConvolutionKernel.from_legacy(self.psfMatchingKernel)), + ): + config = lsst.ip.diffim.SpatiallySampledMetricsTask.ConfigClass() + config.image_type = imageType + task = lsst.ip.diffim.SpatiallySampledMetricsTask(config=config) + butlerQC = _RecordingQuantumContext() + task.runQuantum(butlerQC, types.SimpleNamespace(**inputs, psfMatchingKernel=kernel), None) + metrics[imageType] = butlerQC.put_values.spatiallySampledMetrics + + self.assertGreater(len(metrics["legacy"]), 0) + for name in ("psfMatchingKernel_sum", "psfMatchingKernel_dx", "psfMatchingKernel_dy", + "psfMatchingKernel_residualNorm"): + np.testing.assert_allclose(metrics["future"][name], metrics["legacy"][name], + rtol=1e-6, err_msg=name) + + +class _RecordingQuantumContext: + """Minimal `~lsst.pipe.base.QuantumContext` that records what was put.""" + def __init__(self): + self.put_values = None + + def get(self, refs): + return dict(refs.__dict__) + + def put(self, values, refs): + self.put_values = values + def setup_module(module): lsst.utils.tests.init() diff --git a/tests/test_subtractTask.py b/tests/test_subtractTask.py index ccb3221e6..ab318d9b1 100644 --- a/tests/test_subtractTask.py +++ b/tests/test_subtractTask.py @@ -32,8 +32,11 @@ import lsst.geom import lsst.meas.algorithms as measAlg from lsst.daf.butler import DataCoordinate, DimensionUniverse -from lsst.images import Box, DifferenceImage, DifferenceImageTemplateInfo, Polygon, VisitImage +from lsst.images import (Box, DifferenceImage, DifferenceImageTemplateInfo, Polygon, VisitImage, + get_legacy_difference_image_mask_planes) +from lsst.images.convolution_kernels import ImageBasisConvolutionKernel from lsst.images.fields import ChebyshevField +from lsst.images.tests import compare_masked_image_to_legacy from lsst.ip.diffim import subtractImages, InsufficientKernelSourcesError from lsst.pex.config import FieldValidationError from lsst.pipe.base import InvalidQuantumError, NoWorkFound @@ -1526,6 +1529,28 @@ def _make_images(self): exposure.info.setVisitInfo(makeTestVisitInfo(self.dataId["visit"])) return science, template, sources + def _check_converted_kernel(self, kernel, legacyKernel, bbox): + """Check that a converted kernel evaluates the same as the legacy + kernel it came from. + + Parameters + ---------- + kernel : `lsst.images.convolution_kernels.ImageBasisConvolutionKernel` + The converted kernel to check. + legacyKernel : `lsst.afw.math.LinearCombinationKernel` + The kernel before conversion. + bbox : `lsst.geom.Box2I` + Region of the image the kernel was fit on. + """ + self.assertIsInstance(kernel, ImageBasisConvolutionKernel) + roundTrip = kernel.to_legacy() + for point in (bbox.getMin(), bbox.getCenter(), bbox.getMax()): + expected = afwImage.ImageD(legacyKernel.getDimensions()) + legacyKernel.computeImage(expected, False, point.x, point.y) + actual = afwImage.ImageD(roundTrip.getDimensions()) + roundTrip.computeImage(actual, False, point.x, point.y) + np.testing.assert_allclose(actual.array, expected.array, rtol=0, atol=1e-12) + def _check_converted(self, image, kernelExpected=True, photometricScaling=None): """Check that one converted image is a fully-formed DifferenceImage. @@ -1578,10 +1603,13 @@ def test_connections_legacy(self): if name in ("difference", "matchedTemplate", "scoreExposure"): self.assertEqual(getattr(connections, name).storageClass, "ExposureF") self.assertEqual(connections.science.storageClass, "ExposureF") + for name in ("psfMatchingKernel", "inputPsfMatchingKernel"): + if name in connections.outputs or name in connections.inputs: + self.assertEqual(getattr(connections, name).storageClass, "MatchingKernel") def test_connections_future(self): """Test that ``image_type="future"`` swaps the storage class of - every image output, and only those. + every image and kernel output. """ for configClass, connectionsClass, expected in ( (subtractImages.AlardLuptonSubtractConfig, @@ -1603,10 +1631,8 @@ def test_connections_future(self): # The science image is read in the new type to get its # photometric scaling. self.assertEqual(connections.science.storageClass, "VisitImage") - # The kernel has no converter to or from the new types, so it - # must remain a separate output in both modes. if "psfMatchingKernel" in connections.outputs: - self.assertEqual(connections.psfMatchingKernel.storageClass, "MatchingKernel") + self.assertEqual(connections.psfMatchingKernel.storageClass, "ConvolutionKernel") def test_connections_future_simplified_existing_kernel(self): """Test the storage class swap for the connection variant that @@ -1618,7 +1644,7 @@ def test_connections_future_simplified_existing_kernel(self): connections = subtractImages.SimplifiedSubtractConnections(config=config) self.assertEqual(connections.difference.storageClass, "DifferenceImage") self.assertEqual(connections.matchedTemplate.storageClass, "DifferenceImage") - self.assertEqual(connections.inputPsfMatchingKernel.storageClass, "MatchingKernel") + self.assertEqual(connections.inputPsfMatchingKernel.storageClass, "ConvolutionKernel") self.assertNotIn("psfMatchingKernel", connections.outputs) def test_convert_outputs_to_future(self): @@ -1630,16 +1656,18 @@ def test_convert_outputs_to_future(self): results = task.run(template, science, sources) # The template has no detector, so neither does the matched template. self.assertIsNone(results.matchedTemplate.getDetector()) - legacyDifference = results.difference.image.array.copy() + legacy = {name: getattr(results, name).clone() for name in ("difference", "matchedTemplate")} + legacyKernel = results.psfMatchingKernel task.convert_outputs_to_future(results, self.exposureRecord) - self._check_converted(results.difference) - self._check_converted(results.matchedTemplate) + for name, exposure in legacy.items(): + with self.subTest(name=name): + self._check_converted(getattr(results, name)) + compare_masked_image_to_legacy(getattr(results, name), exposure.maskedImage, + plane_map=get_legacy_difference_image_mask_planes()) self.assertEqual(results.difference.bbox.to_legacy(), science.getBBox()) - np.testing.assert_array_equal(results.difference.image.array, legacyDifference) - # The kernel outputs are untouched. - self.assertIsInstance(results.psfMatchingKernel, afwMath.LinearCombinationKernel) + self._check_converted_kernel(results.psfMatchingKernel, legacyKernel, science.getBBox()) self.assertIsInstance(results.kernelSources, afwTable.SourceCatalog) def test_convert_outputs_to_future_score(self): @@ -1650,11 +1678,12 @@ def test_convert_outputs_to_future_score(self): task = self._setup_subtraction(subtractImages.AlardLuptonPreconvolveSubtractTask, image_type="future") results = task.run(template, science, sources) + legacyKernel = results.psfMatchingKernel task.convert_outputs_to_future(results, self.exposureRecord) self._check_converted(results.scoreExposure) - self.assertIsInstance(results.psfMatchingKernel, afwMath.LinearCombinationKernel) + self._check_converted_kernel(results.psfMatchingKernel, legacyKernel, science.getBBox()) def test_convert_outputs_to_future_simplified(self): """Test converting the outputs of `SimplifiedSubtractTask` running @@ -1845,7 +1874,8 @@ def test_run_quantum_puts_future_types(self): self.assertIsInstance(butlerQC.put_values.matchedTemplate, expectedType) # The kernel is a separate output in both modes. self.assertIsInstance(butlerQC.put_values.psfMatchingKernel, - afwMath.LinearCombinationKernel) + ImageBasisConvolutionKernel if imageType == "future" + else afwMath.LinearCombinationKernel) if imageType == "future": for name in ("difference", "matchedTemplate"): image = getattr(butlerQC.put_values, name) @@ -1883,7 +1913,9 @@ def test_run_quantum_simplified_puts_future_types(self): refs = {"template": self._make_future_template(template, templateInfo), "science": self._make_future_science(science, photometricScaling)} if useExistingKernel: - refs["inputPsfMatchingKernel"] = alResults.psfMatchingKernel + refs["inputPsfMatchingKernel"] = ImageBasisConvolutionKernel.from_legacy( + alResults.psfMatchingKernel + ) butlerQC = _RecordingQuantumContext(self.dataId) task.runQuantum(butlerQC, _FakeRefs(**refs), _FakeRefs()) @@ -1893,6 +1925,11 @@ def test_run_quantum_simplified_puts_future_types(self): self._check_converted(butlerQC.put_values.matchedTemplate, photometricScaling=photometricScaling) self.assertEqual(butlerQC.put_values.difference.templates, templateInfo) + self.assertIsInstance(butlerQC.put_values.psfMatchingKernel, ImageBasisConvolutionKernel) + if useExistingKernel: + # The input kernel is the one the subtraction used. + self._check_converted_kernel(butlerQC.put_values.difference.kernel, + alResults.psfMatchingKernel, science.getBBox()) class _FakeRefs: From 83a5f3de3e91732dee88fec07aefa8bd0c5d099d Mon Sep 17 00:00:00 2001 From: Ian Sullivan Date: Wed, 7 Oct 2026 19:28:11 -0700 Subject: [PATCH 17/17] Copy bbox in GetTemplateTask.run instead of growing it --- python/lsst/ip/diffim/getTemplate.py | 7 ++++--- tests/test_getTemplate.py | 27 ++++++++++----------------- 2 files changed, 14 insertions(+), 20 deletions(-) diff --git a/python/lsst/ip/diffim/getTemplate.py b/python/lsst/ip/diffim/getTemplate.py index 9758caf06..92d18dd77 100644 --- a/python/lsst/ip/diffim/getTemplate.py +++ b/python/lsst/ip/diffim/getTemplate.py @@ -495,9 +495,9 @@ def run(self, *, coaddExposureHandles, bbox, wcs, dataIds, physical_filter, visi `lsst.afw.image.Exposure`]] Coadds to be mosaicked, indexed on tract id. bbox : `lsst.geom.Box2I` - Template Bounding box of the detector geometry onto which to - resample the ``coaddExposureHandles``. Modified in-place to include the - template border. + Input Bounding box of the detector geometry onto which to + resample the ``coaddExposureHandles``. The output template bbox + will be grown by ``templateBorderSize``. wcs : `lsst.afw.geom.SkyWcs` Template WCS onto which to resample the ``coaddExposureHandles``. dataIds : `dict` [`int`, `list` [`lsst.daf.butler.DataCoordinate`]] @@ -527,6 +527,7 @@ def run(self, *, coaddExposureHandles, bbox, wcs, dataIds, physical_filter, visi """ band, photoCalib = self._checkInputs(dataIds, coaddExposureHandles) + bbox = geom.Box2I(bbox) bbox.grow(self.config.templateBorderSize) warped = {} diff --git a/tests/test_getTemplate.py b/tests/test_getTemplate.py index a60f1afa9..1bb0a81be 100644 --- a/tests/test_getTemplate.py +++ b/tests/test_getTemplate.py @@ -267,9 +267,8 @@ def testRunOneTractInput(self): box = lsst.geom.Box2I(lsst.geom.Point2I(0, 0), lsst.geom.Point2I(180, 180)) task = lsst.ip.diffim.GetTemplateTask() # Restrict to tract 0, since the box fits in just that tract. - # Task modifies the input bbox, so pass a copy. result = task.run(coaddExposureHandles={0: self.patches[0]}, - bbox=lsst.geom.Box2I(box), + bbox=box, wcs=self.exposure.wcs, dataIds={0: self.dataIds[0]}, physical_filter="a_test") @@ -285,9 +284,8 @@ def testRunOneTractMultipleInputs(self): """ box = lsst.geom.Box2I(lsst.geom.Point2I(0, 0), lsst.geom.Point2I(180, 180)) task = lsst.ip.diffim.GetTemplateTask() - # Task modifies the input bbox, so pass a copy. result = task.run(coaddExposureHandles=self.patches, - bbox=lsst.geom.Box2I(box), + bbox=box, wcs=self.exposure.wcs, dataIds=self.dataIds, physical_filter="a_test") @@ -301,9 +299,8 @@ def testRunTwoTracts(self): """ box = lsst.geom.Box2I(lsst.geom.Point2I(200, 200), lsst.geom.Point2I(600, 600)) task = lsst.ip.diffim.GetTemplateTask() - # Task modifies the input bbox, so pass a copy. result = task.run(coaddExposureHandles=self.patches, - bbox=lsst.geom.Box2I(box), + bbox=box, wcs=self.exposure.wcs, dataIds=self.dataIds, physical_filter="a_test") @@ -323,9 +320,8 @@ def testRunCalibratesToNanojansky(self): templates = {} for photoCalib in (self.exposure.photoCalib, lsst.afw.image.PhotoCalib(1.0)): self._remakePatches(photoCalib) - # Task modifies the input bbox, so pass a copy. result = task.run(coaddExposureHandles={0: self.patches[0]}, - bbox=lsst.geom.Box2I(box), + bbox=box, wcs=self.exposure.wcs, dataIds={0: self.dataIds[0]}, physical_filter="a_test") @@ -354,7 +350,7 @@ def testGaussianPsfGridFallback(self): box = lsst.geom.Box2I(lsst.geom.Point2I(0, 0), lsst.geom.Point2I(180, 180)) task = lsst.ip.diffim.GetTemplateTask() result = task.run(coaddExposureHandles={0: self.patches[0]}, - bbox=lsst.geom.Box2I(box), + bbox=box, wcs=self.exposure.wcs, dataIds={0: self.dataIds[0]}, physical_filter="a_test") @@ -383,7 +379,7 @@ def testRunNoTemplate(self): task = lsst.ip.diffim.GetTemplateTask() with self.assertRaisesRegex(lsst.pipe.base.NoWorkFound, "No patches found"): task.run(coaddExposureHandles=self.patches, - bbox=lsst.geom.Box2I(box), + bbox=box, wcs=self.exposure.wcs, dataIds=self.dataIds, physical_filter="a_test") @@ -398,9 +394,8 @@ def testMissingPatches(self): self.patches[0].pop(1) box = lsst.geom.Box2I(lsst.geom.Point2I(0, 0), lsst.geom.Point2I(180, 180)) task = lsst.ip.diffim.GetTemplateTask() - # Task modifies the input bbox, so pass a copy. result = task.run(coaddExposureHandles=self.patches, - bbox=lsst.geom.Box2I(box), + bbox=box, wcs=self.exposure.wcs, dataIds=self.dataIds, physical_filter="a_test") @@ -431,7 +426,7 @@ def testNanInputs(self, box=None, nInput=None): box = lsst.geom.Box2I(lsst.geom.Point2I(200, 200), lsst.geom.Point2I(600, 600)) task = lsst.ip.diffim.GetTemplateTask() result = task.run(coaddExposureHandles=self.patches, - bbox=lsst.geom.Box2I(box), + bbox=box, wcs=self.exposure.wcs, dataIds=self.dataIds, physical_filter="a_test") @@ -600,9 +595,8 @@ def _configureAndRunTask(doScaleVariance, varianceScale=1.): config = lsst.ip.diffim.GetTemplateTask.ConfigClass() config.doScaleVariance = doScaleVariance task = lsst.ip.diffim.GetTemplateTask(config=config) - # Task modifies the input bbox, so pass a copy. result = task.run(coaddExposureHandles={0: self._scaleInputVariance(0, varianceScale)}, - bbox=lsst.geom.Box2I(box), + bbox=box, wcs=self.exposure.wcs, dataIds={0: self.dataIds[0]}, physical_filter="a_test") @@ -670,9 +664,8 @@ def _runLegacyForFuture(self, raiseOnUndefinedMaskMap=True): task = lsst.ip.diffim.GetTemplateTask(config=config) self.visit = 9876 box = lsst.geom.Box2I(lsst.geom.Point2I(0, 0), lsst.geom.Point2I(180, 180)) - # Task modifies the input bbox, so pass a copy. result = task.run(coaddExposureHandles={0: self.patches[0]}, - bbox=lsst.geom.Box2I(box), + bbox=box, wcs=self.exposure.wcs, dataIds={0: self.dataIds[0]}, physical_filter="a_test",