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/computeSpatiallySampledMetrics.py b/python/lsst/ip/diffim/computeSpatiallySampledMetrics.py
index dd8949f6e..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 <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,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):
@@ -88,7 +96,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',
@@ -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(
diff --git a/python/lsst/ip/diffim/detectAndMeasure.py b/python/lsst/ip/diffim/detectAndMeasure.py
index 78ea68652..719b80a92 100644
--- a/python/lsst/ip/diffim/detectAndMeasure.py
+++ b/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,10 +616,28 @@ 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,
@@ -603,18 +645,75 @@ 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, 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):
diff --git a/python/lsst/ip/diffim/getTemplate.py b/python/lsst/ip/diffim/getTemplate.py
index 74e7a7588..92d18dd77 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,9 +37,15 @@
 import lsst.pex.exceptions
 import lsst.pipe.base as pipeBase
 
+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
+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
 
@@ -63,6 +70,22 @@ 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"),
+    )
+    obs_info = pipeBase.connectionTypes.Input(
+        doc="Observation metadata of the exposure the template is built for."
+        " 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"),
     )
     wcs = pipeBase.connectionTypes.Input(
         doc="WCS of the exposure that we will construct the template for.",
@@ -96,11 +119,17 @@ 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,
                 deferGraphConstraint=False,
             )
+        if config.output_image_type == "future":
+            self.template = dataclasses.replace(self.template, storageClass="DifferenceImage")
+        else:
+            del self.obs_info
 
 
 class GetTemplateConfig(
@@ -172,6 +201,23 @@ 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",
+    )
+    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
@@ -225,8 +271,10 @@ 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")
+        obsInfo = inputs.pop("obs_info", None)
         coaddExposures = inputs.pop("coaddExposures")
         skymap = inputs.pop("skyMap")
 
@@ -243,8 +291,105 @@ def runQuantum(self, butlerQC, inputRefs, outputRefs):
             physical_filter=physical_filter,
             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,
+                                                              inputRefs.coaddExposures,
+                                                              detector=detector,
+                                                              exposureRecord=exposureRecord)
         butlerQC.put(outputs, outputRefs)
 
+    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`
+        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`
+            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.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
+            for, which supplies the observation metadata recorded on the
+            output.
+
+        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 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.
+        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}
+        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=planeMap,
+        )
+        del template
+        convertedTemplate.templates = DifferenceImageTemplateInfo.from_legacy_psf(
+            convertedTemplate.sky_projection.pixel_frame,
+            coaddPsf,
+            skymaps.pop(),
+            butler_info,
+            log=self.log,
+        )
+        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
@@ -350,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`]]
@@ -382,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 = {}
@@ -864,6 +1010,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,
@@ -895,11 +1047,15 @@ 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")
         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 4c39a924b..3ca166dea 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 <https://www.gnu.org/licenses/>.
 
+import dataclasses
+
 from astropy import units as u
 from astropy.stats import gaussian_fwhm_to_sigma
 import numpy as np
@@ -27,9 +29,12 @@
 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.ip.diffim.utils import (evaluateMeanPsfFwhm, getPsfFwhm,
                                   computeDifferenceImageMetrics, computePSFNoiseEquivalentArea,
-                                  checkMask, 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
@@ -69,13 +74,17 @@ 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"
     )
     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"
@@ -100,6 +109,9 @@ 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")
+            self.science = dataclasses.replace(self.science, storageClass="VisitImage")
 
 
 class SubtractImageOutputConnections(lsst.pipe.base.PipelineTaskConnections,
@@ -130,6 +142,15 @@ 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")
+            self.psfMatchingKernel = dataclasses.replace(self.psfMatchingKernel,
+                                                         storageClass="ConvolutionKernel")
+
 
 class SubtractScoreOutputConnections(lsst.pipe.base.PipelineTaskConnections,
                                      dimensions=_dimensions,
@@ -153,6 +174,13 @@ 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")
+            self.psfMatchingKernel = dataclasses.replace(self.psfMatchingKernel,
+                                                         storageClass="ConvolutionKernel")
+
 
 class AlardLuptonSubtractConnections(SubtractInputConnections, SubtractImageOutputConnections):
     pass
@@ -172,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
 
@@ -219,6 +250,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=(
@@ -302,6 +334,19 @@ 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="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 and"
+                      " lsst.afw.math.Kernel.",
+            "future": "Read and write lsst.images.DifferenceImage and"
+                      " lsst.images.convolution_kernels.ConvolutionKernel.",
+        },
+        optional=False,
+        default="legacy",
+    )
 
     def setDefaults(self):
         self.makeKernel.kernel.name = "AL"
@@ -334,6 +379,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):
@@ -358,6 +412,11 @@ 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 are converted in future
+    mode (`tuple` [`str`]). Subclasses with different image outputs override
+    this."""
+
     def __init__(self, **kwargs):
         super().__init__(**kwargs)
         self.makeSubtask("decorrelate")
@@ -416,16 +475,107 @@ def _applyExternalCalibrations(self, exposure, visitSummary):
 
     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
+            # Read the `TemplateInfo` and the observation metadata from the
+            # template while it is still in the `DifferenceImage` format.
+            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()
+            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)
         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, exposureRecord, templateInfo=templateInfo,
+                                           photometricScaling=photometricScaling)
+
         butlerQC.put(results, outputRefs)
 
+    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
+        ``results`` is replaced by an `lsst.images.DifferenceImage`, with the
+        PSF matching kernel attached to it. ``psfMatchingKernel`` is replaced
+        by the same kernel as an
+        `lsst.images.convolution_kernels.ImageBasisConvolutionKernel`.
+
+        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.
+        templateInfo : `list` [`lsst.images.DifferenceImageTemplateInfo`], \
+                optional
+            Record of the coadds that went into the template, taken from the
+            template this task subtracted.
+        photometricScaling : `lsst.images.fields.BaseField`, optional
+            Photometric scaling of the science image.
+
+        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.
+
+        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
+        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,
+                exposure_record=exposureRecord,
+            )
+            if kernel is not None:
+                image.kernel = kernel
+            if templateInfo:
+                image.templates = templateInfo
+            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):
         """PSF match, subtract, and decorrelate two images.
@@ -1137,10 +1287,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
@@ -1215,6 +1370,9 @@ 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",)
 
     def run(self, template, science, sources, visitSummary=None):
         """Preconvolve the science image with its own PSF,
diff --git a/python/lsst/ip/diffim/utils.py b/python/lsst/ip/diffim/utils.py
index 6e3823d82..9b0940749 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", "get_difference_image_provenance",
            ]
 
 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,106 @@
 _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, usually from the quantum data ID.
+    visit : `int`
+        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``.
+
+    Returns
+    -------
+    record : `lsst.daf.butler.DimensionRecord`
+        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
+    ------
+    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.
+    """
+    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 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_computeSpatiallySampledMetrics.py b/tests/test_computeSpatiallySampledMetrics.py
new file mode 100644
index 000000000..8bb99c1e3
--- /dev/null
+++ b/tests/test_computeSpatiallySampledMetrics.py
@@ -0,0 +1,153 @@
+# 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 <https://www.gnu.org/licenses/>.
+
+import types
+import unittest
+
+import numpy as np
+
+import lsst.afw.image
+import lsst.afw.table
+import lsst.ip.diffim
+import lsst.utils.tests
+from lsst.images.convolution_kernels import ImageBasisConvolutionKernel
+
+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.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 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()
+
+
+class MemoryTestCase(lsst.utils.tests.MemoryTestCase):
+    pass
+
+
+if __name__ == "__main__":
+    lsst.utils.tests.init()
+    unittest.main()
diff --git a/tests/test_detectAndMeasure.py b/tests/test_detectAndMeasure.py
index 58315d35d..a2f26858f 100644
--- a/tests/test_detectAndMeasure.py
+++ b/tests/test_detectAndMeasure.py
@@ -22,19 +22,31 @@
 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,
+                         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 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
@@ -43,7 +55,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:
@@ -1307,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
@@ -1557,6 +1599,478 @@ 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(),
+        )
+        # 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):
+        """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)
+
+    @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
+        ----------
+        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.
+        photometricScaling : `lsst.images.fields.BaseField`, optional
+            Photometric scaling of the science image.
+
+        Returns
+        -------
+        image : `lsst.images.DifferenceImage`
+            The converted image.
+        """
+        image = DifferenceImage.from_legacy(
+            exposure, unit=astropy.units.nJy, exposure_record=self.exposureRecord)
+        if templateInfo is not 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):
+        """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)
+        # 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
+        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, legacy):
+        """Check the properties common to every converted output image.
+
+        Parameters
+        ----------
+        image : `lsst.images.DifferenceImage`
+            The converted image.
+        legacy : `lsst.afw.image.ExposureF`
+            The legacy output exposure it was converted from.
+        """
+        self.assertIsInstance(image, DifferenceImage)
+        self.assertEqual(image.unit, astropy.units.nJy)
+        self.assertIsNotNone(image.detector)
+        # `_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``,
+        and that the background is a separate output only in legacy mode.
+        """
+        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
+                    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)
+                    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."""
+        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 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)
+        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")
+        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, {"subtractedMeasuredExposure": inputImage})
+
+        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.
+
+        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")
+        inputImage = self._makeFutureInput(difference)
+        output = task.run(science, matchedTemplate, inputImage.to_legacy(), sources)
+        legacy = output.subtractedMeasuredExposure
+        for n, plane in enumerate(planes):
+            legacy.mask.addMaskPlane(plane)
+            legacy.mask.array[0, n] |= legacy.mask.getPlaneBitMask(plane)
+
+        task.convert_outputs_to_future(output, {"subtractedMeasuredExposure": inputImage})
+
+        image = output.subtractedMeasuredExposure
+        self.assertIsInstance(image, DifferenceImage)
+        # 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")
+        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.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)
+
+    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")
+        inputImage = self._makeFutureInput(difference)
+        output = task.run(science, matchedTemplate, inputImage.to_legacy(), sources)
+        self.assertEqual(len(output.differenceBackground), 0)
+        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):
+        """``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()
+                photometricScaling = self._makePhotometricScaling(difference)
+                if imageType == "future":
+                    difference = self._makeFutureInput(difference, templateInfo, kernel,
+                                                       photometricScaling)
+                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)
+                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.
+                    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
+        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.
+        """
+        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)
+        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, photometricScaling),
+                                  scoreExposure=self._makeFutureInput(score, templateInfo, scoreKernel,
+                                                                      photometricScaling),
+                                  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)
+        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
+        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..1bb0a81be 100644
--- a/tests/test_getTemplate.py
+++ b/tests/test_getTemplate.py
@@ -24,25 +24,31 @@
 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
+from lsst.images.tests import compare_masked_image_to_legacy
 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
 
-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.
@@ -77,6 +83,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))
 
@@ -260,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")
@@ -278,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")
@@ -294,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")
@@ -316,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")
@@ -347,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")
@@ -376,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")
@@ -391,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")
@@ -424,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")
@@ -593,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")
@@ -638,6 +639,287 @@ def _configureAndRunTask(doScaleVariance, varianceScale=1.):
                                      factor*scaleFactor, rtol=1e-5)
         self.assertImagesAlmostEqual(templateLow.variance, templateOn.variance, rtol=1e-5)
 
+    def _runLegacyForFuture(self, raiseOnUndefinedMaskMap=True):
+        """Build a template from tract 0 with a task configured for the
+        future output type, without converting it.
+
+        The task is kept as ``self.futureTask``.
+
+        Parameters
+        ----------
+        raiseOnUndefinedMaskMap : `bool`, optional
+            Value of the task config field of the same name.
+
+        Returns
+        -------
+        result : `lsst.pipe.base.Struct`
+            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.
+        """
+        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))
+        result = task.run(coaddExposureHandles={0: self.patches[0]},
+                          bbox=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, 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.
+
+        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 pixels and the detector it was built on.
+        """
+        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())
+        # 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)
+        legacyTemplate = result.template.clone()
+
+        result.template = self.futureTask.convert_outputs_to_future(
+            result.template, self._coaddRefs(0), detector=detector,
+            exposureRecord=self._exposureRecord())
+
+        template = result.template
+        self.assertIsInstance(template, lsst.images.DifferenceImage)
+        # 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.
+        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 _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):
+        """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()
+        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`.
+        """
+        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.
+        """
+        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, self._coaddRefs(0), detector=detector,
+            exposureRecord=self._exposureRecord())
+
+        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")
+        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.
+        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.
+        """
+        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..ab318d9b1 100644
--- a/tests/test_subtractTask.py
+++ b/tests/test_subtractTask.py
@@ -20,25 +20,34 @@
 # along with this program.  If not, see <https://www.gnu.org/licenses/>.
 
 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.daf.butler import DataCoordinate, DimensionUniverse
+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 NoWorkFound
+from lsst.pipe.base import InvalidQuantumError, NoWorkFound
 import lsst.utils.tests
 import numpy as np
 from lsst.ip.diffim.utils import (computeRobustStatistics, computePSFNoiseEquivalentArea,
                                   evaluateMeanPsfFwhm, getPsfFwhm)
 from lsst.pex.exceptions import InvalidParameterError
 
-from utils import makeStats, makeTestImage, CustomCoaddPsf
+from utils import (makeStats, makeTestExposureRecord, makeTestImage, makeTestVisitInfo,
+                   CustomCoaddPsf)
 
 
 class AlardLuptonSubtractTestBase:
@@ -64,7 +73,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
@@ -198,6 +206,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 +1469,494 @@ 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.
+    """
+    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
+        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)
+        for exposure in (science, template):
+            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.
+
+        Parameters
+        ----------
+        image : `lsst.images.DifferenceImage`
+            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)
+        # 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"])
+        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
+        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
+        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")
+                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 and kernel output.
+        """
+        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 science image is read in the new type to get its
+                # photometric scaling.
+                self.assertEqual(connections.science.storageClass, "VisitImage")
+                if "psfMatchingKernel" in connections.outputs:
+                    self.assertEqual(connections.psfMatchingKernel.storageClass, "ConvolutionKernel")
+
+    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, "ConvolutionKernel")
+        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())
+        legacy = {name: getattr(results, name).clone() for name in ("difference", "matchedTemplate")}
+        legacyKernel = results.psfMatchingKernel
+
+        task.convert_outputs_to_future(results, self.exposureRecord)
+
+        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())
+        self._check_converted_kernel(results.psfMatchingKernel, legacyKernel, science.getBBox())
+        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)
+        legacyKernel = results.psfMatchingKernel
+
+        task.convert_outputs_to_future(results, self.exposureRecord)
+
+        self._check_converted(results.scoreExposure)
+        self._check_converted_kernel(results.psfMatchingKernel, legacyKernel, science.getBBox())
+
+    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.exposureRecord)
+
+        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.
+        """
+        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.exposureRecord)
+
+        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)
+
+    @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.
+
+        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, exposure_record=self.exposureRecord)
+        if templateInfo is not None:
+            image.templates = templateInfo
+        return image
+
+    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.exposureRecord, 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.exposureRecord)
+
+        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()]
+        photometricScaling = self._makePhotometricScaling(science)
+        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)
+                    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=inputScience,
+                                          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,
+                                      ImageBasisConvolutionKernel if imageType == "future"
+                                      else afwMath.LinearCombinationKernel)
+                if imageType == "future":
+                    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
+        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=self._make_future_science(science), 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.
+        """
+        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):
+            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": self._make_future_science(science, photometricScaling)}
+                if useExistingKernel:
+                    refs["inputPsfMatchingKernel"] = ImageBasisConvolutionKernel.from_legacy(
+                        alResults.psfMatchingKernel
+                    )
+                butlerQC = _RecordingQuantumContext(self.dataId)
+
+                task.runQuantum(butlerQC, _FakeRefs(**refs), _FakeRefs())
+
+                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)
+                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:
+    """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/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..bfc28ff99 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,77 @@ 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 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,
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)