subtractImages.py

Go to the documentation of this file.

# 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 dataclasses
 
from astropy import units as u
from astropy.stats import gaussian_fwhm_to_sigma
import numpy as np
 
import lsst.afw.detection as afwDetection
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, record_from_obs_info, setSourceFootprints,
                                  get_difference_image_provenance)
from lsst.meas.algorithms import ScaleVarianceTask, ScienceSourceSelectorTask
import lsst.pex.config
import lsst.pipe.base
import lsst.pex.exceptions
from lsst.pipe.base import connectionTypes
from . import MakeKernelTask, DecorrelateALKernelTask
from lsst.utils.timer import timeMethod
 
__all__ = ["AlardLuptonSubtractConfig", "AlardLuptonSubtractTask",
           "AlardLuptonPreconvolveSubtractConfig", "AlardLuptonPreconvolveSubtractTask",
           "SimplifiedSubtractConfig", "SimplifiedSubtractTask",
           "InsufficientKernelSourcesError"]
 
_dimensions = ("instrument", "visit", "detector")
_defaultTemplates = {"coaddName": "deep", "fakesType": ""}
 
 
class InsufficientKernelSourcesError(lsst.pipe.base.AlgorithmError):
    """Raised when there are too few sources to calculate the PSF matching
    kernel.
    """
    def __init__(self, *, nSources, nRequired):
        msg = (f"Only {nSources} sources were selected for PSF matching,"
               f" but {nRequired} are required.")
        super().__init__(msg)
        self.nSources = nSources
        self.nRequired = nRequired
 
    @property
    def metadata(self):
        return {"nSources": self.nSources,
                "nRequired": self.nRequired
                }
 
 
class SubtractInputConnections(lsst.pipe.base.PipelineTaskConnections,
                               dimensions=_dimensions,
                               defaultTemplates=_defaultTemplates):
    template = connectionTypes.Input(
        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."
        " 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"
    )
    sources = connectionTypes.Input(
        doc="Sources measured on the science exposure; "
            "used to select sources for making the matching kernel.",
        dimensions=("instrument", "visit", "detector"),
        storageClass="SourceCatalog",
        name="{fakesType}src"
    )
    visitSummary = connectionTypes.Input(
        doc=("Per-visit catalog with final calibration objects. "
             "These catalogs use the detector id for the catalog id, "
             "sorted on id for fast lookup."),
        dimensions=("instrument", "visit"),
        storageClass="ExposureCatalog",
        name="finalVisitSummary",
    )
 
    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,
                                     dimensions=_dimensions,
                                     defaultTemplates=_defaultTemplates):
    difference = connectionTypes.Output(
        doc="Result of subtracting convolved template from science image.",
        dimensions=("instrument", "visit", "detector"),
        storageClass="ExposureF",
        name="{fakesType}{coaddName}Diff_differenceTempExp",
    )
    matchedTemplate = connectionTypes.Output(
        doc="Warped and PSF-matched template used to create `subtractedExposure`.",
        dimensions=("instrument", "visit", "detector"),
        storageClass="ExposureF",
        name="{fakesType}{coaddName}Diff_matchedExp",
    )
    psfMatchingKernel = connectionTypes.Output(
        doc="Kernel used to PSF match the science and template images.",
        dimensions=("instrument", "visit", "detector"),
        storageClass="MatchingKernel",
        name="{fakesType}{coaddName}Diff_psfMatchKernel",
    )
    kernelSources = connectionTypes.Output(
        doc="Final selection of sources used for psf matching.",
        dimensions=("instrument", "visit", "detector"),
        storageClass="SourceCatalog",
        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,
                                     defaultTemplates=_defaultTemplates):
    scoreExposure = connectionTypes.Output(
        doc="The maximum likelihood image, used for the detection of diaSources.",
        dimensions=("instrument", "visit", "detector"),
        storageClass="ExposureF",
        name="{fakesType}{coaddName}Diff_scoreTempExp",
    )
    psfMatchingKernel = connectionTypes.Output(
        doc="Kernel used to PSF match the science and template images.",
        dimensions=("instrument", "visit", "detector"),
        storageClass="MatchingKernel",
        name="{fakesType}{coaddName}Diff_psfScoreMatchKernel",
    )
    kernelSources = connectionTypes.Output(
        doc="Final selection of sources used for psf matching.",
        dimensions=("instrument", "visit", "detector"),
        storageClass="SourceCatalog",
        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
 
 
class SimplifiedSubtractConnections(SubtractInputConnections, SubtractImageOutputConnections):
    inputPsfMatchingKernel = connectionTypes.Input(
        doc="Kernel used to PSF match the science and template images.",
        dimensions=("instrument", "visit", "detector"),
        storageClass="MatchingKernel",
        name="{fakesType}{coaddName}Diff_psfMatchKernel",
    )
 
    def __init__(self, *, config=None):
        super().__init__(config=config)
        del self.sources
        if config.useExistingKernel:
            del self.psfMatchingKernel
            del self.kernelSources
            if config.image_type == "future":
                self.inputPsfMatchingKernel = dataclasses.replace(self.inputPsfMatchingKernel,
                                                                  storageClass="ConvolutionKernel")
        else:
            del self.inputPsfMatchingKernel
 
 
class AlardLuptonSubtractBaseConfig(lsst.pex.config.Config):
    makeKernel = lsst.pex.config.ConfigurableField(
        target=MakeKernelTask,
        doc="Task to construct a matching kernel for convolution.",
    )
    doDecorrelation = lsst.pex.config.Field(
        dtype=bool,
        default=True,
        doc="Perform diffim decorrelation to undo pixel correlation due to A&L "
        "kernel convolution? If True, also update the diffim PSF."
    )
    decorrelate = lsst.pex.config.ConfigurableField(
        target=DecorrelateALKernelTask,
        doc="Task to decorrelate the image difference.",
    )
    requiredTemplateFraction = lsst.pex.config.Field(
        dtype=float,
        default=0.1,
        doc="Raise NoWorkFound and do not attempt image subtraction if template covers less than this "
        " fraction of pixels. Setting to 0 will always attempt image subtraction."
    )
    minTemplateFractionForExpectedSuccess = lsst.pex.config.Field(
        dtype=float,
        default=0.2,
        doc="Raise NoWorkFound if PSF-matching fails and template covers less than this fraction of pixels."
        " If the fraction of pixels covered by the template is less than this value (and greater than"
        " requiredTemplateFraction) this task is attempted but failure is anticipated and tolerated."
    )
    doScaleVariance = lsst.pex.config.Field(
        dtype=bool,
        default=True,
        doc="Scale variance of the science image? Note that the template variance is NOT scaled"
        " here. The template variance may be scaled independently in ``GetTemplateTask``."
    )
    scaleVariance = lsst.pex.config.ConfigurableField(
        target=ScaleVarianceTask,
        doc="Subtask to rescale the variance of the template to the statistically expected level."
    )
    doSubtractBackground = lsst.pex.config.Field(
        doc="Subtract the background fit when solving the kernel? "
        "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=(
            "Replace science Exposure's calibration objects with those"
            " in visitSummary.  Ignored if `doApplyFinalizedPsf is True."
        ),
        dtype=bool,
        default=False,
    )
    sourceSelector = lsst.pex.config.ConfigurableField(
        target=ScienceSourceSelectorTask,
        doc="Task to select sources to be used for PSF matching.",
    )
    fallbackSourceSelector = lsst.pex.config.ConfigurableField(
        target=ScienceSourceSelectorTask,
        doc="Task to select sources to be used for PSF matching."
        "Used only if the kernel calculation fails and"
        "`allowKernelSourceDetection` is set. The fallback source detection"
        " will not include all of the same plugins as the original source "
        " detection, so not all of the same flags can be used.",
    )
    detectionThreshold = lsst.pex.config.Field(
        dtype=float,
        default=10,
        doc="Minimum signal to noise ratio of detected sources "
        "to use for calculating the PSF matching kernel.",
        deprecated="No longer used. Will be removed after v30"
    )
    detectionThresholdMax = lsst.pex.config.Field(
        dtype=float,
        default=500,
        doc="Maximum signal to noise ratio of detected sources "
        "to use for calculating the PSF matching kernel.",
        deprecated="No longer used. Will be removed after v30"
    )
    restrictKernelEdgeSources = lsst.pex.config.Field(
        dtype=bool,
        default=True,
        doc="Exclude sources close to the edge from the kernel calculation?"
    )
    maxKernelSources = lsst.pex.config.Field(
        dtype=int,
        default=1000,
        doc="Maximum number of sources to use for calculating the PSF matching kernel."
        "Set to -1 to disable."
    )
    minKernelSources = lsst.pex.config.Field(
        dtype=int,
        default=3,
        doc="Minimum number of sources needed for calculating the PSF matching kernel."
    )
    excludeMaskPlanes = lsst.pex.config.ListField(
        dtype=str,
        default=("NO_DATA", "BAD", "SAT", "EDGE", "FAKE", "HIGH_VARIANCE"),
        doc="Template mask planes to exclude when selecting sources for PSF matching.",
    )
    badMaskPlanes = lsst.pex.config.ListField(
        dtype=str,
        default=("NO_DATA", "BAD", "SAT", "EDGE"),
        doc="Mask planes to interpolate over."
    )
    preserveTemplateMask = lsst.pex.config.ListField(
        dtype=str,
        default=("NO_DATA", "BAD", "HIGH_VARIANCE"),
        doc="Mask planes from the template to propagate to the image difference."
    )
    renameTemplateMask = lsst.pex.config.ListField(
        dtype=str,
        default=("SAT", "INJECTED", "INJECTED_CORE",),
        doc="Mask planes from the template to propagate to the image difference"
        "with '_TEMPLATE' appended to the name."
    )
    preserveMaskPlanes = lsst.pex.config.ListField(
        dtype=str,
        default=("INJECTED", "INJECTED_CORE", "INJECTED_TEMPLATE", "INJECTED_CORE_TEMPLATE"),
        doc="Mask planes to preserve without dilation when convolving the image.",
    )
    allowKernelSourceDetection = lsst.pex.config.Field(
        dtype=bool,
        default=False,
        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"
        # Always include background fitting in the kernel fit,
        # even if it is not subtracted
        self.makeKernel.kernel.active.fitForBackground = True
        self.makeKernel.kernel.active.spatialKernelOrder = 1
        self.makeKernel.kernel.active.spatialBgOrder = 2
        # Shared source selector settings
        doSkySources = False  # Do not include sky sources
        doSignalToNoise = True  # apply signal to noise filter
        doUnresolved = True  # apply star-galaxy separation
        signalToNoiseMinimum = 10
        signalToNoiseMaximum = 500
        self.sourceSelector.doIsolated = True  # apply isolated star selection
        self.sourceSelector.doRequirePrimary = True  # apply primary flag selection
        self.sourceSelector.doUnresolved = doUnresolved
        self.sourceSelector.doSkySources = doSkySources
        self.sourceSelector.doSignalToNoise = doSignalToNoise
        self.sourceSelector.signalToNoise.minimum = signalToNoiseMinimum
        self.sourceSelector.signalToNoise.maximum = signalToNoiseMaximum
        # The following two configs should not be necessary to be turned on for
        # PSF-matching, and the fallback kernel source selection will fail if
        # they are set since it does not run deblending.
        self.fallbackSourceSelector.doIsolated = False  # Do not apply isolated star selection
        self.fallbackSourceSelector.doRequirePrimary = False  # Do not apply primary flag selection
        self.fallbackSourceSelector.doUnresolved = doUnresolved
        self.fallbackSourceSelector.doSkySources = doSkySources
        self.fallbackSourceSelector.doSignalToNoise = doSignalToNoise
        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):
    mode = lsst.pex.config.ChoiceField(
        dtype=str,
        default="convolveTemplate",
        allowed={"auto": "Choose which image to convolve at runtime.",
                 "convolveScience": "Only convolve the science image.",
                 "convolveTemplate": "Only convolve the template image."},
        doc="Choose which image to convolve at runtime, or require that a specific image is convolved."
    )
 
 
class AlardLuptonSubtractTask(lsst.pipe.base.PipelineTask):
    """Compute the image difference of a science and template image using
    the Alard & Lupton (1998) algorithm.
    """
    ConfigClass = AlardLuptonSubtractConfig
    _DefaultName = "alardLuptonSubtract"
    usePreconvolution = False
    """Whether this task preconvolves the science image with its own PSF
    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")
        self.makeSubtask("makeKernel")
        self.makeSubtask("sourceSelector")
        self.makeSubtask("fallbackSourceSelector")
        if self.config.doScaleVariance:
            self.makeSubtask("scaleVariance")
 
        self.convolutionControl = lsst.afw.math.ConvolutionControl()
        # Normalization is an extra, unnecessary, calculation and will result
        #  in mis-subtraction of the images if there are calibration errors.
        self.convolutionControl.setDoNormalize(False)
        self.convolutionControl.setDoCopyEdge(True)
 
    def _applyExternalCalibrations(self, exposure, visitSummary):
        """Replace calibrations (psf, and ApCorrMap) on this exposure with
        external ones.".
 
        Parameters
        ----------
        exposure : `lsst.afw.image.exposure.Exposure`
            Input exposure to adjust calibrations.
        visitSummary : `lsst.afw.table.ExposureCatalog`
            Exposure catalog with external calibrations to be applied. Catalog
            uses the detector id for the catalog id, sorted on id for fast
            lookup.
 
        Returns
        -------
        exposure : `lsst.afw.image.exposure.Exposure`
            Exposure with adjusted calibrations.
        """
        detectorId = exposure.info.getDetector().getId()
 
        row = visitSummary.find(detectorId)
        if row is None:
            self.log.warning("Detector id %s not found in external calibrations catalog; "
                             "Using original calibrations.", detectorId)
        else:
            psf = row.getPsf()
            apCorrMap = row.getApCorrMap()
            if psf is None:
                self.log.warning("Detector id %s has None for psf in "
                                 "external calibrations catalog; Using original psf and aperture correction.",
                                 detectorId)
            elif apCorrMap is None:
                self.log.warning("Detector id %s has None for apCorrMap in "
                                 "external calibrations catalog; Using original psf and aperture correction.",
                                 detectorId)
            else:
                exposure.setPsf(psf)
                exposure.info.setApCorrMap(apCorrMap)
 
        return exposure
 
    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, 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.
 
        Parameters
        ----------
        template : `lsst.afw.image.ExposureF`
            Template exposure, warped to match the science exposure.
        science : `lsst.afw.image.ExposureF`
            Science exposure to subtract from the template.
        sources : `lsst.afw.table.SourceCatalog`
            Identified sources on the science exposure. This catalog is used to
            select sources in order to perform the AL PSF matching on stamp
            images around them.
        visitSummary : `lsst.afw.table.ExposureCatalog`, optional
            Exposure catalog with external calibrations to be applied. Catalog
            uses the detector id for the catalog id, sorted on id for fast
            lookup.
 
        Returns
        -------
        results : `lsst.pipe.base.Struct`
            ``difference`` : `lsst.afw.image.ExposureF`
                Result of subtracting template and science.
            ``matchedTemplate`` : `lsst.afw.image.ExposureF`
                Warped and PSF-matched template exposure.
            ``backgroundModel`` : `lsst.afw.math.Function2D`
                Background model that was fit while solving for the
                PSF-matching kernel
            ``psfMatchingKernel`` : `lsst.afw.math.Kernel`
                Kernel used to PSF-match the convolved image.
            ``kernelSources` : `lsst.afw.table.SourceCatalog`
                Sources from the input catalog that were used to construct the
                PSF-matching kernel.
        """
        self._prepareInputs(template, science, visitSummary=visitSummary)
 
        convolveTemplate = self.chooseConvolutionMethod(template, science)
        self.matchedPsfSize = self.sciencePsfSize if convolveTemplate else self.templatePsfSize
 
        kernelResult = self.runMakeKernel(template, science, sources=sources,
                                          convolveTemplate=convolveTemplate,
                                          runSourceDetection=False)
 
        if self.config.doSubtractBackground:
            backgroundModel = kernelResult.backgroundModel
        else:
            backgroundModel = None
        if convolveTemplate:
            subtractResults = self.runConvolveTemplate(template, science, kernelResult.psfMatchingKernel,
                                                       backgroundModel=backgroundModel)
        else:
            subtractResults = self.runConvolveScience(template, science, kernelResult.psfMatchingKernel,
                                                      backgroundModel=backgroundModel)
        subtractResults.kernelSources = kernelResult.kernelSources
 
        metrics = computeDifferenceImageMetrics(science, subtractResults.difference, sources)
 
        self.metadata["differenceFootprintRatioMean"] = metrics.differenceFootprintRatioMean
        self.metadata["differenceFootprintRatioStdev"] = metrics.differenceFootprintRatioStdev
        self.metadata["differenceFootprintSkyRatioMean"] = metrics.differenceFootprintSkyRatioMean
        self.metadata["differenceFootprintSkyRatioStdev"] = metrics.differenceFootprintSkyRatioStdev
        self.log.info("Mean, stdev of ratio of difference to science "
                      "pixels in star footprints: %5.4f, %5.4f",
                      self.metadata["differenceFootprintRatioMean"],
                      self.metadata["differenceFootprintRatioStdev"])
 
        return subtractResults
 
    def chooseConvolutionMethod(self, template, science):
        """Determine whether the template should be convolved with the PSF
        matching kernel.
 
        Parameters
        ----------
        template : `lsst.afw.image.ExposureF`
            Template exposure, warped to match the science exposure.
        science : `lsst.afw.image.ExposureF`
            Science exposure to subtract from the template.
 
        Returns
        -------
        convolveTemplate : `bool`
            Convolve the template to match the two images?
 
        Raises
        ------
        RuntimeError
            If an unsupported convolution mode is supplied.
        """
        if self.usePreconvolution:
            raise RuntimeError("Choosing a convolution method is incompatible with preconvolution!")
        if self.config.mode == "auto":
            convolveTemplate = _shapeTest(template,
                                          science,
                                          fwhmExposureBuffer=self.config.makeKernel.fwhmExposureBuffer,
                                          fwhmExposureGrid=self.config.makeKernel.fwhmExposureGrid)
            if convolveTemplate:
                if self.sciencePsfSize < self.templatePsfSize:
                    self.log.info("Average template PSF size is greater, "
                                  "but science PSF greater in one dimension: convolving template image.")
                else:
                    self.log.info("Science PSF size is greater: convolving template image.")
            else:
                self.log.info("Template PSF size is greater: convolving science image.")
        elif self.config.mode == "convolveTemplate":
            self.log.info("`convolveTemplate` is set: convolving template image.")
            convolveTemplate = True
        elif self.config.mode == "convolveScience":
            self.log.info("`convolveScience` is set: convolving science image.")
            convolveTemplate = False
        else:
            raise RuntimeError(f"Cannot handle AlardLuptonSubtract mode: {self.config.mode}")
        return convolveTemplate
 
    def runMakeKernel(self, template, science, sources=None, convolveTemplate=True, runSourceDetection=False):
        """Construct the PSF-matching kernel. Not used for preconvolution.
 
        Parameters
        ----------
        template : `lsst.afw.image.ExposureF`
            Template exposure, warped to match the science exposure.
        science : `lsst.afw.image.ExposureF`
            Science exposure to subtract from the template.
        sources : `lsst.afw.table.SourceCatalog`
            Identified sources on the science exposure. This catalog is used to
            select sources in order to perform the AL PSF matching on stamp
            images around them.
            Not used if ``runSourceDetection`` is set.
        convolveTemplate : `bool`, optional
            Construct the matching kernel to convolve the template?
        runSourceDetection : `bool`, optional
            Run a minimal version of source detection to determine kernel
            candidates? If False, a source list to select kernel candidates
            from must be supplied.
 
        Returns
        -------
        results : `lsst.pipe.base.Struct`
            ``backgroundModel`` : `lsst.afw.math.Function2D`
                Background model that was fit while solving for the
                PSF-matching kernel
            ``psfMatchingKernel`` : `lsst.afw.math.Kernel`
                Kernel used to PSF-match the convolved image.
            ``kernelSources` : `lsst.afw.table.SourceCatalog`
                Sources from the input catalog that were used to construct the
                PSF-matching kernel.
        """
        if self.usePreconvolution:
            raise RuntimeError("Incorrect matching kernel calculation configured. "
                               "`runMakeKernel` can't be called if `usePreconvolution` is set.")
        if convolveTemplate:
            reference = template
            target = science
            referenceFwhmPix = self.templatePsfSize
            targetFwhmPix = self.sciencePsfSize
        else:
            reference = science
            target = template
            referenceFwhmPix = self.sciencePsfSize
            targetFwhmPix = self.templatePsfSize
        try:
            if runSourceDetection:
                kernelSources = self.runKernelSourceDetection(template, science)
            else:
                kernelSources = self._sourceSelector(template, science, sources)
            kernelResult = self.makeKernel.run(reference, target, kernelSources,
                                               preconvolved=False,
                                               templateFwhmPix=referenceFwhmPix,
                                               scienceFwhmPix=targetFwhmPix)
        except (RuntimeError, lsst.pex.exceptions.Exception) as e:
            self.log.warning("Failed to match template. Checking coverage")
            #  Raise NoWorkFound if template fraction is insufficient
            checkTemplateIsSufficient(template[science.getBBox()], science, self.log,
                                      self.config.minTemplateFractionForExpectedSuccess,
                                      exceptionMessage="Template coverage lower than expected to succeed."
                                      f" Failure is tolerable: {e}")
            #  checkTemplateIsSufficient did not raise NoWorkFound, so raise original exception
            raise e
 
        return lsst.pipe.base.Struct(backgroundModel=kernelResult.backgroundModel,
                                     psfMatchingKernel=kernelResult.psfMatchingKernel,
                                     kernelSources=kernelSources)
 
    def runKernelSourceDetection(self, template, science):
        """Run detection on the science image and use the template mask plane
        to reject candidate sources.
 
        Parameters
        ----------
        template : `lsst.afw.image.ExposureF`
            Template exposure, warped to match the science exposure.
        science : `lsst.afw.image.ExposureF`
            Science exposure to subtract from the template.
 
        Returns
        -------
        kernelSources : `lsst.afw.table.SourceCatalog`
                Sources from the input catalog to use to construct the
                PSF-matching kernel.
        """
        kernelSize = self.makeKernel.makeKernelBasisList(
            self.templatePsfSize, self.matchedPsfSize)[0].getWidth()
        sources = self.makeKernel.makeCandidateList(template, science, kernelSize,
                                                    candidateList=None,
                                                    sigma=gaussian_fwhm_to_sigma*self.sciencePsfSize)
        return self._sourceSelector(template, science, sources, fallback=True)
 
    def runConvolveTemplate(self, template, science, psfMatchingKernel, backgroundModel=None):
        """Convolve the template image with a PSF-matching kernel and subtract
        from the science image.
 
        Parameters
        ----------
        template : `lsst.afw.image.ExposureF`
            Template exposure, warped to match the science exposure.
        science : `lsst.afw.image.ExposureF`
            Science exposure to subtract from the template.
        psfMatchingKernel : `lsst.afw.math.Kernel`
            Kernel to be used to PSF-match the science image to the template.
        backgroundModel : `lsst.afw.math.Function2D`, optional
            Background model that was fit while solving for the PSF-matching
            kernel.
 
        Returns
        -------
        results : `lsst.pipe.base.Struct`
 
            ``difference`` : `lsst.afw.image.ExposureF`
                Result of subtracting template and science.
            ``matchedTemplate`` : `lsst.afw.image.ExposureF`
                Warped and PSF-matched template exposure.
            ``backgroundModel`` : `lsst.afw.math.Function2D`
                Background model that was fit while solving for the PSF-matching kernel
            ``psfMatchingKernel`` : `lsst.afw.math.Kernel`
                Kernel used to PSF-match the template to the science image.
        """
        self.metadata["convolvedExposure"] = "Template"
 
        matchedTemplate = self._convolveExposure(template, psfMatchingKernel,
                                                 self.convolutionControl,
                                                 bbox=science.getBBox(),
                                                 psf=science.psf,
                                                 photoCalib=science.photoCalib)
 
        difference = _subtractImages(science, matchedTemplate, backgroundModel=backgroundModel)
        correctedExposure = self.finalize(template, science, difference,
                                          psfMatchingKernel,
                                          templateMatched=True)
 
        return lsst.pipe.base.Struct(difference=correctedExposure,
                                     matchedTemplate=matchedTemplate,
                                     matchedScience=science,
                                     backgroundModel=backgroundModel,
                                     psfMatchingKernel=psfMatchingKernel)
 
    def runConvolveScience(self, template, science, psfMatchingKernel, backgroundModel=None):
        """Convolve the science image with a PSF-matching kernel and subtract
        the template image.
 
        Parameters
        ----------
        template : `lsst.afw.image.ExposureF`
            Template exposure, warped to match the science exposure.
        science : `lsst.afw.image.ExposureF`
            Science exposure to subtract from the template.
        psfMatchingKernel : `lsst.afw.math.Kernel`
            Kernel to be used to PSF-match the science image to the template.
        backgroundModel : `lsst.afw.math.Function2D`, optional
            Background model that was fit while solving for the PSF-matching
            kernel.
 
        Returns
        -------
        results : `lsst.pipe.base.Struct`
 
            ``difference`` : `lsst.afw.image.ExposureF`
                Result of subtracting template and science.
            ``matchedTemplate`` : `lsst.afw.image.ExposureF`
                Warped template exposure. Note that in this case, the template
                is not PSF-matched to the science image.
            ``backgroundModel`` : `lsst.afw.math.Function2D`
                Background model that was fit while solving for the PSF-matching kernel
            ``psfMatchingKernel`` : `lsst.afw.math.Kernel`
               Kernel used to PSF-match the science image to the template.
        """
        self.metadata["convolvedExposure"] = "Science"
        bbox = science.getBBox()
 
        kernelImage = lsst.afw.image.ImageD(psfMatchingKernel.getDimensions())
        xcen, ycen = bbox.getCenter()
        norm = psfMatchingKernel.computeImage(kernelImage, doNormalize=False, x=xcen, y=ycen)
 
        matchedScience = self._convolveExposure(science, psfMatchingKernel,
                                                self.convolutionControl,
                                                psf=template.psf)
 
        # Place back on native photometric scale
        matchedScience.maskedImage /= norm
        matchedTemplate = template.clone()[bbox]
        matchedTemplate.setPhotoCalib(science.photoCalib)
 
        if backgroundModel is not None:
            # We must invert the background model if the matching kernel is solved for the science image.
            invertedBackground = backgroundModel.clone()
            invertedBackground.setParameters([-p for p in backgroundModel.getParameters()])
            backgroundModel = invertedBackground
 
        difference = _subtractImages(matchedScience, matchedTemplate, backgroundModel=backgroundModel)
 
        correctedExposure = self.finalize(template, science, difference,
                                          psfMatchingKernel,
                                          templateMatched=False)
 
        return lsst.pipe.base.Struct(difference=correctedExposure,
                                     matchedTemplate=matchedTemplate,
                                     matchedScience=matchedScience,
                                     backgroundModel=backgroundModel,
                                     psfMatchingKernel=psfMatchingKernel)
 
    def finalize(self, template, science, difference, kernel,
                 templateMatched=True,
                 preConvMode=False,
                 preConvKernel=None,
                 spatiallyVarying=False):
        """Decorrelate the difference image to undo the noise correlations
        caused by convolution.
 
        Parameters
        ----------
        template : `lsst.afw.image.ExposureF`
            Template exposure, warped to match the science exposure.
        science : `lsst.afw.image.ExposureF`
            Science exposure to subtract from the template.
        difference : `lsst.afw.image.ExposureF`
            Result of subtracting template and science.
        kernel : `lsst.afw.math.Kernel`
            An (optionally spatially-varying) PSF matching kernel
        templateMatched : `bool`, optional
            Was the template PSF-matched to the science image?
        preConvMode : `bool`, optional
            Was the science image preconvolved with its own PSF
            before PSF matching the template?
        preConvKernel : `lsst.afw.detection.Psf`, optional
            If not `None`, then the science image was pre-convolved with
            (the reflection of) this kernel. Must be normalized to sum to 1.
        spatiallyVarying : `bool`, optional
            Compute the decorrelation kernel spatially varying across the image?
 
        Returns
        -------
        correctedExposure : `lsst.afw.image.ExposureF`
            The decorrelated image difference.
        """
        if self.config.doDecorrelation:
            self.log.info("Decorrelating image difference.")
            # We have cleared the template mask plane, so copy the mask plane of
            # the image difference so that we can calculate correct statistics
            # during decorrelation
            correctedExposure = self.decorrelate.run(science, template[science.getBBox()], difference, kernel,
                                                     templateMatched=templateMatched,
                                                     preConvMode=preConvMode,
                                                     preConvKernel=preConvKernel,
                                                     spatiallyVarying=spatiallyVarying).correctedExposure
        else:
            self.log.info("NOT decorrelating image difference.")
            correctedExposure = difference
        return correctedExposure
 
    def _calculateMagLim(self, exposure, nsigma=5.0, fallbackPsfSize=None):
        """Calculate the point-source limiting magnitude of an exposure.
 
        The limit is the flux of a point source detected at ``nsigma`` in
        the median pixel noise, integrated over the effective area of the
        PSF, so it does not depend on the units of the pixels.
 
        Parameters
        ----------
        exposure : `lsst.afw.image.Exposure`
            The exposure to calculate the limiting magnitude for.
        nsigma : `float`, optional
            The detection threshold in sigma.
        fallbackPsfSize : `float`, optional
            PSF FWHM in pixels to use if the exposure PSF cannot be
            evaluated.
 
        Returns
        -------
        maglim : `float`
            The limiting AB magnitude, or `numpy.nan` if the exposure has no
            photometric calibration, no usable pixels, or no PSF.
        """
        if exposure.photoCalib is None:
            return np.nan
        try:
            psfArea = computePSFNoiseEquivalentArea(exposure.getPsf())
        except (lsst.pex.exceptions.InvalidParameterError,
                afwDetection.InvalidPsfError,
                lsst.pex.exceptions.RangeError):
            if fallbackPsfSize is None:
                self.log.info("Unable to evaluate PSF, setting maglim to nan")
                return np.nan
            self.log.info("Unable to evaluate PSF, using fallback FWHM %f", fallbackPsfSize)
            # Effective area of a Gaussian PSF with this FWHM.
            psfArea = 4*np.pi*(fallbackPsfSize/np.sqrt(8*np.log(2)))**2
        mask = exposure.mask
        good = (mask.array & mask.getPlaneBitMask(self.config.badMaskPlanes)) == 0
        variance = exposure.variance.array[good]
        variance = variance[np.isfinite(variance) & (variance > 0)]
        if variance.size == 0:
            return np.nan
        fluxLim = nsigma*np.sqrt(np.median(variance)*psfArea)
        return exposure.photoCalib.instFluxToMagnitude(fluxLim)
 
    @staticmethod
    def _validateExposures(template, science):
        """Check that the WCS of the two Exposures match, the template bbox
        contains the science bbox, and that the bands match.
 
        Parameters
        ----------
        template : `lsst.afw.image.ExposureF`
            Template exposure, warped to match the science exposure.
        science : `lsst.afw.image.ExposureF`
            Science exposure to subtract from the template.
 
        Raises
        ------
        AssertionError
            Raised if the WCS of the template is not equal to the science WCS,
            if the science image is not fully contained in the template
            bounding box, or if the bands do not match.
        """
        assert template.wcs == science.wcs, \
            "Template and science exposure WCS are not identical."
        templateBBox = template.getBBox()
        scienceBBox = science.getBBox()
        assert science.filter.bandLabel == template.filter.bandLabel, \
            "Science and template exposures have different bands: %s, %s" % \
            (science.filter, template.filter)
 
        assert templateBBox.contains(scienceBBox), \
            "Template bbox does not contain all of the science image."
 
    def _convolveExposure(self, exposure, kernel, convolutionControl,
                          bbox=None,
                          psf=None,
                          photoCalib=None,
                          interpolateBadMaskPlanes=False,
                          ):
        """Convolve an exposure with the given kernel.
 
        Parameters
        ----------
        exposure : `lsst.afw.Exposure`
            exposure to convolve.
        kernel : `lsst.afw.math.LinearCombinationKernel`
            PSF matching kernel computed in the ``makeKernel`` subtask.
        convolutionControl : `lsst.afw.math.ConvolutionControl`
            Configuration for convolve algorithm.
        bbox : `lsst.geom.Box2I`, optional
            Bounding box to trim the convolved exposure to.
        psf : `lsst.afw.detection.Psf`, optional
            Point spread function (PSF) to set for the convolved exposure.
        photoCalib : `lsst.afw.image.PhotoCalib`, optional
            Photometric calibration of the convolved exposure.
        interpolateBadMaskPlanes : `bool`, optional
            If set, interpolate over mask planes specified in
            ``config.badMaskPlanes`` before convolving the image.
 
        Returns
        -------
        convolvedExp : `lsst.afw.Exposure`
            The convolved image.
        """
        convolvedExposure = exposure.clone()
        if psf is not None:
            convolvedExposure.setPsf(psf)
        if photoCalib is not None:
            convolvedExposure.setPhotoCalib(photoCalib)
        if interpolateBadMaskPlanes and self.config.badMaskPlanes is not None:
            nInterp = _interpolateImage(convolvedExposure.maskedImage,
                                        self.config.badMaskPlanes)
            self.metadata["nInterpolated"] = nInterp
 
        # Snapshot the footprints of mask planes that must not be dilated by
        # the convolution.
        preservePlanes = [mp for mp in self.config.preserveMaskPlanes
                          if mp in convolvedExposure.mask.getMaskPlaneDict()]
        maskResetDict = {
            mp: (convolvedExposure.mask.array
                 & convolvedExposure.mask.getPlaneBitMask(mp)) > 0
            for mp in preservePlanes
        }
 
        convolvedImage = lsst.afw.image.MaskedImageF(convolvedExposure.getBBox())
        lsst.afw.math.convolve(convolvedImage, convolvedExposure.maskedImage, kernel, convolutionControl)
        convolvedExposure.setMaskedImage(convolvedImage)
 
        # Undo the convolution's dilation of the preserved planes: clear the
        # dilated bits, then restore each mask plane for the pixels that were
        # previously set.
        self._clearMask(convolvedExposure.mask, clearMaskPlanes=preservePlanes)
        for maskPlane, maskSetPixels in maskResetDict.items():
            bit = convolvedExposure.mask.getPlaneBitMask(maskPlane)
            convolvedExposure.mask.array[maskSetPixels] |= bit
 
        if bbox is None:
            return convolvedExposure
        else:
            return convolvedExposure[bbox]
 
    def _sourceSelector(self, template, science, sources, fallback=False):
        """Select sources from a catalog that meet the selection criteria.
        The selection criteria include any configured parameters of the
        `sourceSelector` subtask, as well as checking the science and template
        mask planes.
 
        Parameters
        ----------
        template : `lsst.afw.image.ExposureF`
            Template exposure, warped to match the science exposure.
        science : `lsst.afw.image.ExposureF`
            Science exposure to subtract from the template.
        sources : `lsst.afw.table.SourceCatalog`
            Input source catalog to select sources from.
        fallback : `bool`, optional
            Switch indicating the source selector is being called after
            running the fallback source detection subtask, which does not run a
            full set of measurement plugins and can't use the same settings for
            the source selector.
 
        Returns
        -------
        kernelSources : `lsst.afw.table.SourceCatalog`
            The input source catalog, with flagged and low signal-to-noise
            sources removed and footprints added.
 
        Raises
        ------
        InsufficientKernelSourcesError
            An AlgorithmError that is raised if there are not enough PSF
            candidates to construct the PSF matching kernel.
        """
        if fallback:
            selected = self.fallbackSourceSelector.selectSources(sources).selected
        else:
            selected = self.sourceSelector.selectSources(sources).selected
        # It is OK to use just self.matchedPsfSize for the science PSF here,
        # since we are just using it to calculate the size of the matching
        # kernel.
        kSize = self.makeKernel.makeKernelBasisList(self.templatePsfSize, self.matchedPsfSize)[0].getWidth()
        selectSources = sources[selected].copy(deep=True)
        # Set the footprints, to be used in `makeKernel` and `checkMask`.
        kernelSources = setSourceFootprints(selectSources, kernelSize=kSize)
        bbox = science.getBBox()
        if self.usePreconvolution:
            # Exclude a wider buffer around the edge of the image to
            # account for an extra convolution.
            bbox.grow(-kSize)
        if self.config.restrictKernelEdgeSources:
            bbox.grow(-kSize)
        # Remove sources that land on masked pixels
        scienceSelected = checkMask(science.mask[bbox], kernelSources, self.config.excludeMaskPlanes)
        templateSelected = checkMask(template.mask[bbox], kernelSources, self.config.excludeMaskPlanes)
        maskSelected = scienceSelected & templateSelected
        kernelSources = kernelSources[maskSelected].copy(deep=True)
        # Trim kernelSources if they exceed ``maxKernelSources``.
        # Keep the highest signal-to-noise sources of those selected.
        if (len(kernelSources) > self.config.maxKernelSources) & (self.config.maxKernelSources > 0):
            signalToNoise = kernelSources.getPsfInstFlux()/kernelSources.getPsfInstFluxErr()
            indices = np.argsort(signalToNoise)
            indices = indices[-self.config.maxKernelSources:]
            selected = np.zeros(len(kernelSources), dtype=bool)
            selected[indices] = True
            kernelSources = kernelSources[selected].copy(deep=True)
 
        self.log.info("%i/%i=%.1f%% of sources selected for PSF matching from the input catalog",
                      len(kernelSources), len(sources), 100*len(kernelSources)/len(sources))
        if len(kernelSources) < self.config.minKernelSources:
            self.log.error("Too few sources to calculate the PSF matching kernel: "
                           "%i selected but %i needed for the calculation.",
                           len(kernelSources), self.config.minKernelSources)
            if self.config.allowKernelSourceDetection and not fallback:
                # The fallback source detection pipeline calls this method, so
                # allowing source detection in that case would create an endless
                # loop
                kernelSources = self.runKernelSourceDetection(template, science)
            else:
                raise InsufficientKernelSourcesError(nSources=len(kernelSources),
                                                     nRequired=self.config.minKernelSources)
 
        self.metadata["nPsfSources"] = len(kernelSources)
 
        return kernelSources
 
    def _prepareInputs(self, template, science, visitSummary=None):
        """Perform preparatory calculations common to all Alard&Lupton Tasks.
 
        Parameters
        ----------
        template : `lsst.afw.image.ExposureF`
            Template exposure, warped to match the science exposure. The
            variance plane of the template image is modified in place.
        science : `lsst.afw.image.ExposureF`
            Science exposure to subtract from the template. The variance plane
            of the science image is modified in place.
        visitSummary : `lsst.afw.table.ExposureCatalog`, optional
            Exposure catalog with external calibrations to be applied.  Catalog
            uses the detector id for the catalog id, sorted on id for fast
            lookup.
        """
        self._validateExposures(template, science)
        if visitSummary is not None:
            self._applyExternalCalibrations(science, visitSummary=visitSummary)
        templateCoverageFraction = checkTemplateIsSufficient(
            template[science.getBBox()], science, self.log,
            requiredTemplateFraction=self.config.requiredTemplateFraction,
            exceptionMessage="Not attempting subtraction. To force subtraction,"
            " set config requiredTemplateFraction=0"
        )
        self.metadata["templateCoveragePercent"] = 100*templateCoverageFraction
 
        if self.config.doScaleVariance:
            # Scale the variance of the science image before
            # convolution, subtraction, or decorrelation so that it has the
            # correct ratio. Note that the template variance is scaled
            # independently in ``GetTemplateTask``.
            sciVarFactor = self.scaleVariance.run(science.maskedImage)
            self.log.info("Science variance scaling factor: %.2f", sciVarFactor)
            self.metadata["scaleScienceVarianceFactor"] = sciVarFactor
 
        # Erase existing detection mask planes.
        #  We don't want the detection mask from the science image
        self.updateMasks(template, science)
 
        # Calling getPsfFwhm on template.psf fails on some rare occasions when
        # the template has no input exposures at the average position of the
        # stars. So we try getPsfFwhm first on template, and if that fails we
        # evaluate the PSF on a grid specified by fwhmExposure* fields.
        # To keep consistent definitions for PSF size on the template and
        # science images, we use the same method for both.
        # In the try block below, we catch two exceptions:
        # 1. InvalidParameterError, in case the point where we are evaluating
        #    the PSF lands in a gap in the template.
        # 2. RangeError, in case the template coverage is so poor that we end
        #    up near a region with no data.
        try:
            self.templatePsfSize = getPsfFwhm(template.psf)
            self.sciencePsfSize = getPsfFwhm(science.psf)
        except lsst.pex.exceptions.Exception:
            # Catch a broad range of exceptions, since some are C++ only
            # Catching:
            #  - lsst::geom::SingularTransformException
            #  - lsst.pex.exceptions.InvalidParameterError
            #  - lsst.pex.exceptions.RangeError
            self.log.info("Unable to evaluate PSF at the average position. "
                          "Evaluting PSF on a grid of points."
                          )
            self.templatePsfSize = evaluateMeanPsfFwhm(
                template,
                fwhmExposureBuffer=self.config.makeKernel.fwhmExposureBuffer,
                fwhmExposureGrid=self.config.makeKernel.fwhmExposureGrid
            )
            self.sciencePsfSize = evaluateMeanPsfFwhm(
                science,
                fwhmExposureBuffer=self.config.makeKernel.fwhmExposureBuffer,
                fwhmExposureGrid=self.config.makeKernel.fwhmExposureGrid
            )
        self.log.info("Science PSF FWHM: %f pixels", self.sciencePsfSize)
        self.log.info("Template PSF FWHM: %f pixels", self.templatePsfSize)
        self.metadata["sciencePsfSize"] = self.sciencePsfSize
        self.metadata["templatePsfSize"] = self.templatePsfSize
 
        #  Calculate estimated image depths, i.e., limiting magnitudes
        maglim_science = self._calculateMagLim(science, fallbackPsfSize=self.sciencePsfSize)
        if np.isnan(maglim_science):
            self.log.warning("Limiting magnitude of the science image is NaN!")
        fluxlim_science = (maglim_science*u.ABmag).to_value(u.nJy)
        maglim_template = self._calculateMagLim(template, fallbackPsfSize=self.templatePsfSize)
        if np.isnan(maglim_template):
            self.log.info("Cannot evaluate template limiting mag; adopting science limiting mag for diffim")
            maglim_diffim = maglim_science
        else:
            fluxlim_template = (maglim_template*u.ABmag).to_value(u.nJy)
            maglim_diffim = (np.sqrt(fluxlim_science**2 + fluxlim_template**2)*u.nJy).to(u.ABmag).value
        self.metadata["scienceLimitingMagnitude"] = maglim_science
        self.metadata["templateLimitingMagnitude"] = maglim_template
        self.metadata["diffimLimitingMagnitude"] = maglim_diffim
 
    def updateMasks(self, template, science):
        """Update the science and template mask planes before differencing.
 
        Parameters
        ----------
        template : `lsst.afw.image.Exposure`
            Template exposure, warped to match the science exposure.
            The template mask planes will be erased, except for a few specified
            in the task config.
        science : `lsst.afw.image.Exposure`
            Science exposure to subtract from the template.
            The DETECTED and DETECTED_NEGATIVE mask planes of the science image
            will be erased.
        """
        self._clearMask(science.mask, clearMaskPlanes=["DETECTED", "DETECTED_NEGATIVE"])
 
        # We will clear ALL template mask planes, except for those specified
        # via the `preserveTemplateMask` config. Mask planes specified via
        # the `renameTemplateMask` config will be copied to new planes with
        # "_TEMPLATE" appended to their names, and the original mask plane will
        # be cleared.
        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
        # NOTE: This can be removed in DM-40796
        if "FAKE" in science.mask.getMaskPlaneDict().keys():
            self.log.info("Adding injected mask plane to science image")
            self._renameMaskPlanes(science.mask, "FAKE", "INJECTED")
        if "FAKE" in template.mask.getMaskPlaneDict().keys():
            self.log.info("Adding injected mask plane to template image")
            self._renameMaskPlanes(template.mask, "FAKE", "INJECTED_TEMPLATE")
            if "INJECTED" in renameMaskPlanes:
                renameMaskPlanes.remove("INJECTED")
            if "INJECTED_TEMPLATE" in clearMaskPlanes:
                clearMaskPlanes.remove("INJECTED_TEMPLATE")
 
        for maskPlane in renameMaskPlanes:
            self._renameMaskPlanes(template.mask, maskPlane, maskPlane + "_TEMPLATE")
        self._clearMask(template.mask, clearMaskPlanes=clearMaskPlanes)
 
    @staticmethod
    def _renameMaskPlanes(mask, maskPlane, newMaskPlane):
        """Rename a mask plane by adding the new name and copying the data.
 
        Parameters
        ----------
        mask : `lsst.afw.image.Mask`
            The mask image to update in place.
        maskPlane : `str`
            The name of the existing mask plane to copy.
        newMaskPlane : `str`
            The new name of the mask plane that will be added.
            If the mask plane already exists, it will be updated in place.
        """
        mask.addMaskPlane(newMaskPlane)
        originBitMask = mask.getPlaneBitMask(maskPlane)
        destinationBitMask = mask.getPlaneBitMask(newMaskPlane)
        mask.array |= ((mask.array & originBitMask) > 0)*destinationBitMask
 
    def _clearMask(self, mask, clearMaskPlanes=None):
        """Clear the mask plane of an exposure.
 
        Parameters
        ----------
        mask : `lsst.afw.image.Mask`
            The mask plane to erase, which will be modified in place.
        clearMaskPlanes : `list` of `str`, optional
            Erase the specified mask planes.
            If not supplied, the entire mask will be erased.
        """
        if clearMaskPlanes is None:
            clearMaskPlanes = list(mask.getMaskPlaneDict().keys())
 
        bitMaskToClear = mask.getPlaneBitMask(clearMaskPlanes)
        mask &= ~bitMaskToClear
 
 
class AlardLuptonPreconvolveSubtractConnections(SubtractInputConnections,
                                                SubtractScoreOutputConnections):
    pass
 
 
class AlardLuptonPreconvolveSubtractConfig(AlardLuptonSubtractBaseConfig, lsst.pipe.base.PipelineTaskConfig,
                                           pipelineConnections=AlardLuptonPreconvolveSubtractConnections):
    pass
 
 
class AlardLuptonPreconvolveSubtractTask(AlardLuptonSubtractTask):
    """Subtract a template from a science image, convolving the science image
    before computing the kernel, and also convolving the template before
    subtraction.
    """
    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,
        convolve the template image with a PSF-matching kernel and subtract
        from the preconvolved science image.
 
        Parameters
        ----------
        template : `lsst.afw.image.ExposureF`
            The template image, which has previously been warped to the science
            image. The template bbox will be padded by a few pixels compared to
            the science bbox.
        science : `lsst.afw.image.ExposureF`
            The science exposure.
        sources : `lsst.afw.table.SourceCatalog`
            Identified sources on the science exposure. This catalog is used to
            select sources in order to perform the AL PSF matching on stamp
            images around them.
        visitSummary : `lsst.afw.table.ExposureCatalog`, optional
            Exposure catalog with complete external calibrations. Catalog uses
            the detector id for the catalog id, sorted on id for fast lookup.
 
        Returns
        -------
        results : `lsst.pipe.base.Struct`
            ``scoreExposure`` : `lsst.afw.image.ExposureF`
                Result of subtracting the convolved template and science
                images. Attached PSF is that of the original science image.
            ``matchedTemplate`` : `lsst.afw.image.ExposureF`
                Warped and PSF-matched template exposure. Attached PSF is that
                of the original science image.
            ``matchedScience`` : `lsst.afw.image.ExposureF`
                The science exposure after convolving with its own PSF.
                Attached PSF is that of the original science image.
            ``backgroundModel`` : `lsst.afw.math.Function2D`
                Background model that was fit while solving for the
                PSF-matching kernel
            ``psfMatchingKernel`` : `lsst.afw.math.Kernel`
                Final kernel used to PSF-match the template to the science
                image.
        """
        self._prepareInputs(template, science, visitSummary=visitSummary)
 
        convolutionKernel = self._makePreconvolutionKernel(science.psf)
        matchedScience = self._convolveExposure(science, convolutionKernel, self.convolutionControl,
                                                interpolateBadMaskPlanes=True)
        self.metadata["convolvedExposure"] = "Preconvolution"
 
        self.matchedPsfSize = self.sciencePsfSize*np.sqrt(2)
        self.log.info("Preconvolved science PSF FWHM: %f pixels", self.matchedPsfSize)
        self.metadata["preconvolvedSciencePsfSize"] = self.matchedPsfSize
        try:
            kernelSources = self._sourceSelector(template, matchedScience, sources)
            subtractResults = self.runPreconvolve(template, science, matchedScience,
                                                  kernelSources, convolutionKernel)
 
        except (RuntimeError, lsst.pex.exceptions.Exception) as e:
            self.log.warning("Failed to match template. Checking coverage")
            #  Raise NoWorkFound if template fraction is insufficient
            checkTemplateIsSufficient(template[science.getBBox()], science, self.log,
                                      self.config.minTemplateFractionForExpectedSuccess,
                                      exceptionMessage="Template coverage lower than expected to succeed."
                                      f" Failure is tolerable: {e}")
            #  checkTemplateIsSufficient did not raise NoWorkFound, so raise original exception
            raise e
 
        return subtractResults
 
    @staticmethod
    def _flagScoreEdge(mask, innerBBox):
        """Set the EDGE mask bit on pixels outside a known-valid region.
 
        Parameters
        ----------
        mask : `~lsst.afw.image.Mask`
            Exposure mask that will be modified in place. Must have
            an ``EDGE`` mask plane.
        innerBBox : `~lsst.geom.Box2I`
            The valid inner region. Pixels
            outside this bbox will have their ``EDGE`` bit set.
        """
        bbox = mask.getBBox()
        edgeBit = mask.getPlaneBitMask("EDGE")
        dx0 = innerBBox.getMinX() - bbox.getMinX()
        dx1 = bbox.getMaxX() - innerBBox.getMaxX()
        dy0 = innerBBox.getMinY() - bbox.getMinY()
        dy1 = bbox.getMaxY() - innerBBox.getMaxY()
        if dy0 > 0:
            mask.array[:dy0, :] |= edgeBit
        if dy1 > 0:
            mask.array[-dy1:, :] |= edgeBit
        if dx0 > 0:
            mask.array[:, :dx0] |= edgeBit
        if dx1 > 0:
            mask.array[:, -dx1:] |= edgeBit
 
    @staticmethod
    def _makePreconvolutionKernel(psf):
        """Build a normalized, reflected matched-filter kernel from a PSF.
 
        Convolving an image with this kernel is equivalent to correlating
        the image with the PSF, so peaks in the output align with the PSF's
        centroid — even for asymmetric PSFs. The kernel is evaluated at the
        PSF's average position and returned as a constant
        `~lsst.afw.math.Kernel`.
 
        Parameters
        ----------
        psf : `~lsst.afw.detection.Psf`
            The PSF to derive the preconvolution kernel from.
 
        Returns
        -------
        kernel : `~lsst.afw.math.Kernel`
            The PSF reflected about both axes, normalized to sum to one.
 
        Raises
        ------
        ValueError
            Raised if the PSF kernel has an even size along either axis.
            It's not possible to center an even-sized kernel.
        """
        avgPos = psf.getAveragePosition()
        localKernel = psf.getLocalKernel(avgPos)
        dims = localKernel.getDimensions()
        if dims.x % 2 == 0 or dims.y % 2 == 0:
            raise ValueError(
                f"Preconvolution requires an odd-sized PSF kernel, got {dims.x}x{dims.y}. "
            )
        kimg = lsst.afw.image.ImageD(dims)
        localKernel.computeImage(kimg, doNormalize=True)  # normalize to unit sum
        # Reflect about the kernel center. PSF kernels are odd-sized,
        # so ``[::-1, ::-1]`` places the peak at the same pixel.
        kimg.array[...] = kimg.array[::-1, ::-1]
        return lsst.afw.math.FixedKernel(kimg)
 
    def runPreconvolve(self, template, science, matchedScience, kernelSources, preConvKernel):
        """Convolve the science image with its own PSF, then convolve the
        template with a matching kernel and subtract to form the Score
        exposure.
 
        Parameters
        ----------
        template : `lsst.afw.image.ExposureF`
            Template exposure, warped to match the science exposure.
        science : `lsst.afw.image.ExposureF`
            Science exposure to subtract from the template.
        matchedScience : `lsst.afw.image.ExposureF`
            The science exposure, convolved with the reflection of its own PSF.
        kernelSources : `lsst.afw.table.SourceCatalog`
            Identified sources on the science exposure. This catalog is used to
            select sources in order to perform the AL PSF matching on stamp
            images around them.
        preConvKernel : `lsst.afw.math.Kernel`
            The kernel that was used to preconvolve the ``science``
            exposure. Must be normalized to sum to 1.
 
        Returns
        -------
        results : `lsst.pipe.base.Struct`
 
            ``scoreExposure`` : `lsst.afw.image.ExposureF`
                Result of subtracting the convolved template and science
                images. Attached PSF is that of the original science image.
            ``matchedTemplate`` : `lsst.afw.image.ExposureF`
                Warped and PSF-matched template exposure. Attached PSF is that
                of the original science image.
            ``matchedScience`` : `lsst.afw.image.ExposureF`
                The science exposure after convolving with its own PSF.
                Attached PSF is that of the original science image.
            ``backgroundModel`` : `lsst.afw.math.Function2D`
                Background model that was fit while solving for the
                PSF-matching kernel
            ``psfMatchingKernel`` : `lsst.afw.math.Kernel`
                Final kernel used to PSF-match the template to the science
                image.
        """
        bbox = science.getBBox()
        innerBBox = preConvKernel.shrinkBBox(bbox)
 
        kernelResult = self.makeKernel.run(template[innerBBox], matchedScience[innerBBox], kernelSources,
                                           preconvolved=True,
                                           templateFwhmPix=self.templatePsfSize,
                                           scienceFwhmPix=self.matchedPsfSize)
 
        matchedTemplate = self._convolveExposure(template, kernelResult.psfMatchingKernel,
                                                 self.convolutionControl,
                                                 bbox=bbox,
                                                 psf=science.psf,
                                                 interpolateBadMaskPlanes=True,
                                                 photoCalib=science.photoCalib)
        score = _subtractImages(matchedScience, matchedTemplate,
                                backgroundModel=(kernelResult.backgroundModel
                                                 if self.config.doSubtractBackground else None))
        correctedScore = self.finalize(template[bbox], science, score,
                                       kernelResult.psfMatchingKernel,
                                       templateMatched=True, preConvMode=True,
                                       preConvKernel=preConvKernel)
 
        # Flag the outer ``preConvKernel/2``-wide border as EDGE.
        self._flagScoreEdge(correctedScore.mask, innerBBox)
 
        return lsst.pipe.base.Struct(scoreExposure=correctedScore,
                                     matchedTemplate=matchedTemplate,
                                     matchedScience=matchedScience,
                                     backgroundModel=kernelResult.backgroundModel,
                                     psfMatchingKernel=kernelResult.psfMatchingKernel,
                                     kernelSources=kernelSources)
 
 
def checkTemplateIsSufficient(templateExposure, scienceExposure, logger, requiredTemplateFraction=0.,
                              exceptionMessage=""):
    """Raise NoWorkFound if template coverage < requiredTemplateFraction
 
    Parameters
    ----------
    templateExposure : `lsst.afw.image.ExposureF`
        The template exposure to check
    logger : `logging.Logger`
        Logger for printing output.
    requiredTemplateFraction : `float`, optional
        Fraction of pixels of the science image required to have coverage
        in the template.
    exceptionMessage : `str`, optional
        Message to include in the exception raised if the template coverage
        is insufficient.
 
    Returns
    -------
    templateCoverageFraction: `float`
        Fraction of pixels in the template with data.
 
    Raises
    ------
    lsst.pipe.base.NoWorkFound
        Raised if fraction of good pixels, defined as not having NO_DATA
        set, is less than the requiredTemplateFraction
    """
    # Count the number of pixels with the NO_DATA mask bit set
    # counting NaN pixels is insufficient because pixels without data are often intepolated over)
    noTemplate = templateExposure.mask.array & templateExposure.mask.getPlaneBitMask('NO_DATA')
    # Also need to account for missing data in the science image,
    # because template coverage there doesn't help
    noScience = scienceExposure.mask.array & scienceExposure.mask.getPlaneBitMask('NO_DATA')
    pixNoData = np.count_nonzero(noTemplate | noScience)
    pixGood = templateExposure.getBBox().getArea() - pixNoData
    templateCoverageFraction = pixGood/templateExposure.getBBox().getArea()
    logger.info("template has %d good pixels (%.1f%%)", pixGood, 100*templateCoverageFraction)
 
    if templateCoverageFraction < requiredTemplateFraction:
        message = ("Insufficient Template Coverage. (%.1f%% < %.1f%%)" % (
                   100*templateCoverageFraction,
                   100*requiredTemplateFraction))
        raise lsst.pipe.base.NoWorkFound(message + " " + exceptionMessage)
    return templateCoverageFraction
 
 
def _subtractImages(science, template, backgroundModel=None):
    """Subtract template from science, propagating relevant metadata.
 
    Parameters
    ----------
    science : `lsst.afw.Exposure`
        The input science image.
    template : `lsst.afw.Exposure`
        The template to subtract from the science image.
    backgroundModel : `lsst.afw.MaskedImage`, optional
        Differential background model
 
    Returns
    -------
    difference : `lsst.afw.Exposure`
        The subtracted image.
    """
    difference = science.clone()
    if backgroundModel is not None:
        difference.maskedImage -= backgroundModel
    difference.maskedImage -= template.maskedImage
    return difference
 
 
def _shapeTest(exp1, exp2, fwhmExposureBuffer, fwhmExposureGrid):
    """Determine that the PSF of ``exp1`` is not wider than that of ``exp2``.
 
    Parameters
    ----------
    exp1 : `~lsst.afw.image.Exposure`
        Exposure with the reference point spread function (PSF) to evaluate.
    exp2 : `~lsst.afw.image.Exposure`
        Exposure with a candidate point spread function (PSF) to evaluate.
    fwhmExposureBuffer : `float`
        Fractional buffer margin to be left out of all sides of the image
        during the construction of the grid to compute mean PSF FWHM in an
        exposure, if the PSF is not available at its average position.
    fwhmExposureGrid : `int`
        Grid size to compute the mean FWHM in an exposure, if the PSF is not
        available at its average position.
    Returns
    -------
    result : `bool`
        True if ``exp1`` has a PSF that is not wider than that of ``exp2`` in
        either dimension.
    """
    try:
        shape1 = getPsfFwhm(exp1.psf, average=False)
        shape2 = getPsfFwhm(exp2.psf, average=False)
    except (lsst.pex.exceptions.InvalidParameterError, lsst.pex.exceptions.RangeError):
        shape1 = evaluateMeanPsfFwhm(exp1,
                                     fwhmExposureBuffer=fwhmExposureBuffer,
                                     fwhmExposureGrid=fwhmExposureGrid
                                     )
        shape2 = evaluateMeanPsfFwhm(exp2,
                                     fwhmExposureBuffer=fwhmExposureBuffer,
                                     fwhmExposureGrid=fwhmExposureGrid
                                     )
        return shape1 <= shape2
 
    # Results from getPsfFwhm is a tuple of two values, one for each dimension.
    xTest = shape1[0] <= shape2[0]
    yTest = shape1[1] <= shape2[1]
    return xTest | yTest
 
 
class SimplifiedSubtractConfig(AlardLuptonSubtractBaseConfig, lsst.pipe.base.PipelineTaskConfig,
                               pipelineConnections=SimplifiedSubtractConnections):
    mode = lsst.pex.config.ChoiceField(
        dtype=str,
        default="convolveTemplate",
        allowed={"auto": "Choose which image to convolve at runtime.",
                 "convolveScience": "Only convolve the science image.",
                 "convolveTemplate": "Only convolve the template image."},
        doc="Choose which image to convolve at runtime, or require that a specific image is convolved."
    )
    useExistingKernel = lsst.pex.config.Field(
        dtype=bool,
        default=True,
        doc="Use a pre-existing PSF matching kernel?"
        "If False, source detection and measurement will be run."
    )
 
 
class SimplifiedSubtractTask(AlardLuptonSubtractTask):
    """Compute the image difference of a science and template image using
    the Alard & Lupton (1998) algorithm.
    """
    ConfigClass = SimplifiedSubtractConfig
    _DefaultName = "simplifiedSubtract"
 
    @timeMethod
    def run(self, template, science, visitSummary=None, inputPsfMatchingKernel=None):
        """PSF match, subtract, and decorrelate two images.
 
        Parameters
        ----------
        template : `lsst.afw.image.ExposureF`
            Template exposure, warped to match the science exposure.
        science : `lsst.afw.image.ExposureF`
            Science exposure to subtract from the template.
        visitSummary : `lsst.afw.table.ExposureCatalog`, optional
            Exposure catalog with external calibrations to be applied. Catalog
            uses the detector id for the catalog id, sorted on id for fast
            lookup.
        inputPsfMatchingKernel : `lsst.afw.math.Kernel`, optional
            Pre-existing PSF matching kernel to use for convolution.
            Required, and only used, if ``config.useExistingKernel`` is set.
 
        Returns
        -------
        results : `lsst.pipe.base.Struct`
            ``difference`` : `lsst.afw.image.ExposureF`
                Result of subtracting template and science.
            ``matchedTemplate`` : `lsst.afw.image.ExposureF`
                Warped and PSF-matched template exposure.
            ``backgroundModel`` : `lsst.afw.math.Function2D`
                Background model that was fit while solving for the
                PSF-matching kernel
            ``psfMatchingKernel`` : `lsst.afw.math.Kernel`
                Kernel used to PSF-match the convolved image.
            ``kernelSources` : `lsst.afw.table.SourceCatalog`
                Sources detected on the science image that were used to
                construct the PSF-matching kernel.
 
        Raises
        ------
        lsst.pipe.base.NoWorkFound
            Raised if fraction of good pixels, defined as not having NO_DATA
            set, is less then the configured requiredTemplateFraction
        """
        self._prepareInputs(template, science, visitSummary=visitSummary)
 
        convolveTemplate = self.chooseConvolutionMethod(template, science)
        self.matchedPsfSize = self.sciencePsfSize if convolveTemplate else self.templatePsfSize
 
        if self.config.useExistingKernel:
            psfMatchingKernel = inputPsfMatchingKernel
            backgroundModel = None
            kernelSources = None
        else:
            kernelResult = self.runMakeKernel(template, science, convolveTemplate=convolveTemplate,
                                              runSourceDetection=True)
            psfMatchingKernel = kernelResult.psfMatchingKernel
            kernelSources = kernelResult.kernelSources
            if self.config.doSubtractBackground:
                backgroundModel = kernelResult.backgroundModel
            else:
                backgroundModel = None
        if convolveTemplate:
            subtractResults = self.runConvolveTemplate(template, science, psfMatchingKernel,
                                                       backgroundModel=backgroundModel)
        else:
            subtractResults = self.runConvolveScience(template, science, psfMatchingKernel,
                                                      backgroundModel=backgroundModel)
        if kernelSources is not None:
            subtractResults.kernelSources = kernelSources
        return subtractResults
 
 
def _interpolateImage(maskedImage, badMaskPlanes, fallbackValue=None):
    """Replace masked image pixels with interpolated values.
 
    Parameters
    ----------
    maskedImage : `lsst.afw.image.MaskedImage`
        Image on which to perform interpolation.
    badMaskPlanes : `list` of `str`
        List of mask planes to interpolate over.
    fallbackValue : `float`, optional
        Value to set when interpolation fails.
 
    Returns
    -------
    result: `float`
        The number of masked pixels that were replaced.
    """
    imgBadMaskPlanes = [
        maskPlane for maskPlane in badMaskPlanes if maskPlane in maskedImage.mask.getMaskPlaneDict()
    ]
 
    image = maskedImage.image.array
    badPixels = (maskedImage.mask.array & maskedImage.mask.getPlaneBitMask(imgBadMaskPlanes)) > 0
    image[badPixels] = np.nan
    if fallbackValue is None:
        fallbackValue = np.nanmedian(image)
    # For this initial implementation, skip the interpolation and just fill with
    # the median value.
    image[badPixels] = fallbackValue
    return np.sum(badPixels)