DM-56114: Add optional lsst.image format for output images

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

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 <https://www.gnu.org/licenses/>.
 
+import dataclasses
+
 import numpy as np
 import scipy.signal
 
@@ -29,7 +31,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
@@ -81,14 +83,21 @@ 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):
     """Config for SpatiallySampledMetricsTask
     """
     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',
@@ -99,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"]
@@ -215,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
@@ -247,12 +273,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 +305,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 +355,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(

python/lsst/ip/diffim/detectAndMeasure.py

@@ -19,18 +19,23 @@
 # You should have received a copy of the GNU General Public License
 # along with this program.  If not, see <https://www.gnu.org/licenses/>.
 
+import dataclasses
 import numpy as np
 import requests
 import os
 
+
 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 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)
+                                  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
 from lsst.meas.base import ForcedMeasurementTask, ApplyApCorrTask, DetectorVisitIdGeneratorConfig
@@ -179,10 +184,16 @@ 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")
+        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,
@@ -212,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,
@@ -399,6 +413,16 @@ class DetectAndMeasureConfig(pipeBase.PipelineTaskConfig,
             "closest in time to the exposure."
     )
     idGenerator = DetectorVisitIdGeneratorConfig.make_field()
+    image_type = pexConfig.ChoiceField[str](
+        "Image type of the input and measured difference images (and, for"
+        " DetectAndMeasureScoreTask, the input and measured score images).",
+        allowed={
+            "legacy": "Read and write lsst.afw.image.ExposureF.",
+            "future": "Read and write lsst.images.DifferenceImage.",
+        },
+        optional=False,
+        default="legacy",
+    )
 
     def setDefaults(self):
         # Background subtraction
@@ -592,29 +616,104 @@ 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.
+
+        # Input image that each output is measured from, keyed by output name.
+        inputImages = {}
+        if self.config.image_type == "future":
+            dataId = butlerQC.quantum.dataId
+            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.
         measurementResults = pipeBase.Struct(
             subtractedMeasuredExposure=None,
+            scoreMeasuredExposure=None,
             diaSources=None,
             maskedStreaks=None,
             differenceBackground=None,
         )
         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, inputImages)
             butlerQC.put(measurementResults, outputRefs)
             raise error from e
+        if self.config.image_type == "future":
+            self.convert_outputs_to_future(measurementResults, inputImages)
         butlerQC.put(measurementResults, outputRefs)
 
+    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.
+
+        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.
+        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
+        -----
+        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, image in inputImages.items():
+            exposure = getattr(results, name, None)
+            if exposure is None:
+                continue
+            # 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.
+                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):
@@ -782,16 +881,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.
@@ -1483,6 +1585,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):