isrFunctions.py

Go to the documentation of this file.

#
# LSST Data Management System
# Copyright 2008, 2009, 2010 LSST Corporation.
#
# This product includes software developed by the
# LSST Project (http://www.lsst.org/).
#
# 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 LSST License Statement and
# the GNU General Public License along with this program.  If not,
# see <http://www.lsstcorp.org/LegalNotices/>.
#
 
__all__ = [
    "applyGains",
    "attachTransmissionCurve",
    "biasCorrection",
    "checkFilter",
    "compareCameraKeywords",
    "countMaskedPixels",
    "createPsf",
    "darkCorrection",
    "flatCorrection",
    "gainContext",
    "getPhysicalFilter",
    "growMasks",
    "maskDECamEdgeBleed",
    "maskE2VEdgeBleed",
    "maskITLEdgeBleed",
    "maskITLSatSag",
    "maskITLDip",
    "illuminationCorrection",
    "interpolateDefectList",
    "interpolateFromMask",
    "makeThresholdMask",
    "saturationCorrection",
    "setBadRegions",
    "transposeMaskedImage",
    "trimToMatchCalibBBox",
    "updateVariance",
    "widenSaturationTrails",
    "getExposureGains",
    "getExposureReadNoises",
]
 
import logging
import math
import numpy
 
import lsst.geom
import lsst.afw.image as afwImage
import lsst.afw.detection as afwDetection
import lsst.afw.math as afwMath
import lsst.meas.algorithms as measAlg
import lsst.afw.cameraGeom as camGeom
 
from lsst.afw.geom import SpanSet, Stencil
from lsst.meas.algorithms.detection import SourceDetectionTask
 
from contextlib import contextmanager
 
from .defects import Defects
 
 
def createPsf(fwhm):
    """Make a double Gaussian PSF.
 
    Parameters
    ----------
    fwhm : scalar
        FWHM of double Gaussian smoothing kernel.
 
    Returns
    -------
    psf : `lsst.meas.algorithms.DoubleGaussianPsf`
        The created smoothing kernel.
    """
    ksize = 4*int(fwhm) + 1
    return measAlg.DoubleGaussianPsf(ksize, ksize, fwhm/(2*math.sqrt(2*math.log(2))))
 
 
def transposeMaskedImage(maskedImage):
    """Make a transposed copy of a masked image.
 
    Parameters
    ----------
    maskedImage : `lsst.afw.image.MaskedImage`
        Image to process.
 
    Returns
    -------
    transposed : `lsst.afw.image.MaskedImage`
        The transposed copy of the input image.
    """
    transposed = maskedImage.Factory(lsst.geom.Extent2I(maskedImage.getHeight(), maskedImage.getWidth()))
    transposed.getImage().getArray()[:] = maskedImage.getImage().getArray().T
    transposed.getMask().getArray()[:] = maskedImage.getMask().getArray().T
    transposed.getVariance().getArray()[:] = maskedImage.getVariance().getArray().T
    return transposed
 
 
def interpolateDefectList(maskedImage, defectList, fwhm, fallbackValue=None,
                          maskNameList=None, useLegacyInterp=True):
    """Interpolate over defects specified in a defect list.
 
    Parameters
    ----------
    maskedImage : `lsst.afw.image.MaskedImage`
        Image to process.
    defectList : `lsst.meas.algorithms.Defects`
        List of defects to interpolate over.
    fwhm : `float`
        FWHM of double Gaussian smoothing kernel.
    fallbackValue : scalar, optional
        Fallback value if an interpolated value cannot be determined.
        If None, then the clipped mean of the image is used.
    maskNameList : `list [string]`
        List of the defects to interpolate over (used for GP interpolator).
    useLegacyInterp : `bool`
        Use the legacy interpolation (polynomial interpolation) if True. Use
        Gaussian Process interpolation if False.
 
    Notes
    -----
    The ``fwhm`` parameter is used to create a PSF, but the underlying
    interpolation code (`lsst.meas.algorithms.interpolateOverDefects`) does
    not currently make use of this information in legacy Interpolation, but use
    if for the Gaussian Process as an estimation of the correlation lenght.
    """
    psf = createPsf(fwhm)
    if fallbackValue is None:
        fallbackValue = afwMath.makeStatistics(maskedImage.getImage(), afwMath.MEANCLIP).getValue()
    if 'INTRP' not in maskedImage.getMask().getMaskPlaneDict():
        maskedImage.getMask().addMaskPlane('INTRP')
 
    # Hardcoded fwhm value. PSF estimated latter in step1,
    # not in ISR.
    if useLegacyInterp:
        kwargs = {}
        fwhm = fwhm
    else:
        # tested on a dozens of images and looks a good set of
        # hyperparameters, but cannot guarrenty this is optimal,
        # need further testing.
        kwargs = {"bin_spacing": 20,
                  "threshold_dynamic_binning": 2000,
                  "threshold_subdivide": 20000}
        fwhm = 15
 
    measAlg.interpolateOverDefects(maskedImage, psf, defectList,
                                   fallbackValue=fallbackValue,
                                   useFallbackValueAtEdge=True,
                                   fwhm=fwhm,
                                   useLegacyInterp=useLegacyInterp,
                                   maskNameList=maskNameList, **kwargs)
    return maskedImage
 
 
def makeThresholdMask(maskedImage, threshold, growFootprints=1, maskName='SAT'):
    """Mask pixels based on threshold detection.
 
    Parameters
    ----------
    maskedImage : `lsst.afw.image.MaskedImage`
        Image to process.  Only the mask plane is updated.
    threshold : scalar
        Detection threshold.
    growFootprints : scalar, optional
        Number of pixels to grow footprints of detected regions.
    maskName : str, optional
        Mask plane name, or list of names to convert
 
    Returns
    -------
    defectList : `lsst.meas.algorithms.Defects`
        Defect list constructed from pixels above the threshold.
    """
    # find saturated regions
    thresh = afwDetection.Threshold(threshold)
    fs = afwDetection.FootprintSet(maskedImage, thresh)
 
    if growFootprints > 0:
        fs = afwDetection.FootprintSet(fs, rGrow=growFootprints, isotropic=False)
    fpList = fs.getFootprints()
 
    # set mask
    mask = maskedImage.getMask()
    bitmask = mask.getPlaneBitMask(maskName)
    afwDetection.setMaskFromFootprintList(mask, fpList, bitmask)
 
    return Defects.fromFootprintList(fpList)
 
 
def growMasks(mask, radius=0, maskNameList=['BAD'], maskValue="BAD"):
    """Grow a mask by an amount and add to the requested plane.
 
    Parameters
    ----------
    mask : `lsst.afw.image.Mask`
        Mask image to process.
    radius : scalar
        Amount to grow the mask.
    maskNameList : `str` or `list` [`str`]
        Mask names that should be grown.
    maskValue : `str`
        Mask plane to assign the newly masked pixels to.
    """
    if radius > 0:
        spans = SpanSet.fromMask(mask, mask.getPlaneBitMask(maskNameList))
        # Use MANHATTAN for equivalence with 'isotropic=False` footprint grows,
        # but CIRCLE is probably better and might be just as fast.
        spans = spans.dilated(radius, Stencil.MANHATTAN)
        spans = spans.clippedTo(mask.getBBox())
        spans.setMask(mask, mask.getPlaneBitMask(maskValue))
 
 
def maskE2VEdgeBleed(exposure, e2vEdgeBleedSatMinArea=10000,
                     e2vEdgeBleedSatMaxArea=100000,
                     e2vEdgeBleedYMax=350,
                     saturatedMaskName="SAT", log=None):
    """Mask edge bleeds in E2V detectors.
 
    Parameters
    ----------
    exposure : `lsst.afw.image.Exposure`
        Exposure to apply masking to.
    e2vEdgeBleedSatMinArea : `int`, optional
        Minimum limit of saturated cores footprint area.
    e2vEdgeBleedSatMaxArea : `int`, optional
        Maximum limit of saturated cores footprint area.
    e2vEdgeBleedYMax: `float`, optional
        Height of edge bleed masking.
    saturatedMaskName : `str`, optional
        Mask name for saturation.
    log : `logging.Logger`, optional
        Logger to handle messages.
    """
 
    log = log if log else logging.getLogger(__name__)
 
    maskedImage = exposure.maskedImage
    saturatedBit = maskedImage.mask.getPlaneBitMask(saturatedMaskName)
 
    thresh = afwDetection.Threshold(saturatedBit, afwDetection.Threshold.BITMASK)
 
    fpList = afwDetection.FootprintSet(exposure.mask, thresh).getFootprints()
 
    satAreas = numpy.asarray([fp.getArea() for fp in fpList])
    largeAreas, = numpy.where((satAreas >= e2vEdgeBleedSatMinArea)
                              & (satAreas < e2vEdgeBleedSatMaxArea))
    for largeAreasIndex in largeAreas:
        fpCore = fpList[largeAreasIndex]
        xCore, yCore = fpCore.getCentroid()
        xCore = int(xCore)
        yCore = int(yCore)
 
        for amp in exposure.getDetector():
            if amp.getBBox().contains(xCore, yCore):
                ampName = amp.getName()
                if ampName[:2] == 'C0':
                    # Check that the footprint reaches the bottom of the
                    # amplifier.
                    if fpCore.getBBox().getMinY() == 0:
                        # This is a large saturation footprint that hits the
                        # edge, and is thus classified as an edge bleed.
 
                        # TODO DM-50587: Optimize number of rows to mask by
                        # looking at the median signal level as a function of
                        # row number on the right side of the saturation trail.
 
                        log.info("Found E2V edge bleed in amp %s, column %d.", ampName, xCore)
                        maskedImage.mask[amp.getBBox()].array[:e2vEdgeBleedYMax, :] |= saturatedBit
 
 
def maskDECamEdgeBleed(exposure, satMinArea=10000, satMaxArea=100000, approachRows=20,
                       nSigma=5.0, nRowsCheck=20, minLowPixelsPerRow=30, minLowPixelsExtent=10,
                       marginFraction=0.125, saturatedMaskName="SAT", log=None):
    """Mask DECam-style edge bleeds with low rows next to the read register.
 
    For each amplifier that contains a large saturated footprint reaching
    within ``approachRows`` of its read edge, this function confirms the dip
    from the per-row count of pixels more than ``nSigma`` below the clipped
    sky, measures how far the dip extends inward, and marks those rows
    (plus a margin) with the saturation mask plane.
 
    Parameters
    ----------
    exposure : `lsst.afw.image.Exposure`
        Assembled exposure to mask.
    satMinArea : `int`, optional
        Minimum area (pixels) of a saturated footprint to be considered.
    satMaxArea : `int`, optional
        Maximum area (pixels) of a saturated footprint to be considered.
    approachRows : `int`, optional
        A footprint must come within this many rows of the read edge.
    nSigma : `float`, optional
        A pixel is "low" if it is more than this many sigma below sky.
    nRowsCheck : `int`, optional
        Number of rows from the read edge used to confirm the bleed.
    minLowPixelsPerRow : `int`, optional
        Mean number of low pixels per row over the check rows required
        to confirm the bleed.
    minLowPixelsExtent : `int`, optional
        Number of low pixels a row must exceed to count toward the bleed
        height; the scan stops after five consecutive rows that do not.
        Should be smaller than ``minLowPixelsPerRow``.
    marginFraction : `float`, optional
        Extra rows masked beyond the measured height, as a fraction of
        that height (plus one row).
    saturatedMaskName : `str`, optional
        Mask plane name for saturation.
    log : `logging.Logger`, optional
        Logger to handle messages.
 
    Raises
    ------
    RuntimeError
        Raised if the exposure has no detector with amplifiers.
    """
    log = log if log else logging.getLogger(__name__)
 
    detector = exposure.getDetector()
    if detector is None or len(detector) == 0:
        raise RuntimeError("Cannot mask DECam edge bleeds: exposure has no detector with amplifiers.")
 
    maskedImage = exposure.maskedImage
    saturatedBit = maskedImage.mask.getPlaneBitMask(saturatedMaskName)
    unusableBits = maskedImage.mask.getPlaneBitMask([saturatedMaskName, "BAD", "NO_DATA"])
    statsBits = unusableBits | maskedImage.mask.getPlaneBitMask("SUSPECT")
 
    thresh = afwDetection.Threshold(saturatedBit, afwDetection.Threshold.BITMASK)
    fpList = afwDetection.FootprintSet(maskedImage.mask, thresh).getFootprints()
    candidates = []
    for fp in fpList:
        if satMinArea <= fp.getArea() < satMaxArea:
            xCore, yCore = fp.getCentroid()
            candidates.append((int(xCore), int(yCore), fp.getBBox()))
    if not candidates:
        log.debug("No saturated footprints in the DECam edge bleed area range.")
        return
 
    statsCtrl = afwMath.StatisticsControl()
    statsCtrl.setAndMask(statsBits)
    # Number of consecutive rows without a dip that ends the height scan.
    nRowsStop = 5
 
    for amp in detector:
        ampBBox = amp.getBBox()
        if not exposure.getBBox().contains(ampBBox):
            continue
        readTop = amp.getReadoutCorner() in (camGeom.ReadoutCorner.UL, camGeom.ReadoutCorner.UR)
 
        if readTop:
            reachesReadEdge = any(ampBBox.contains(x, y)
                                  and fpBBox.getMaxY() >= ampBBox.getMaxY() - approachRows
                                  for x, y, fpBBox in candidates)
        else:
            reachesReadEdge = any(ampBBox.contains(x, y)
                                  and fpBBox.getMinY() <= ampBBox.getMinY() + approachRows
                                  for x, y, fpBBox in candidates)
        if not reachesReadEdge:
            continue
 
        ampImage = maskedImage[ampBBox]
        stats = afwMath.makeStatistics(ampImage, afwMath.MEANCLIP | afwMath.STDEVCLIP | afwMath.NPOINT,
                                       statsCtrl)
        if stats.getValue(afwMath.NPOINT) < 1000:
            log.debug("Skipping DECam edge bleed check in amp %s: too few unmasked pixels.",
                      amp.getName())
            continue
        sky = stats.getValue(afwMath.MEANCLIP)
        sigma = stats.getValue(afwMath.STDEVCLIP)
 
        usable = (ampImage.mask.array & unusableBits) == 0
        low = (ampImage.image.array < sky - nSigma*sigma) & usable
        # Row profiles with index 0 at the read edge, increasing inward.
        profile = low.sum(axis=1)
        nUsable = usable.sum(axis=1)
        if readTop:
            profile = profile[::-1]
            nUsable = nUsable[::-1]
 
        # Skip edge rows that cannot hold enough low pixels to confirm a dip.
        enoughUsable = numpy.nonzero(nUsable >= minLowPixelsPerRow)[0]
        if len(enoughUsable) == 0:
            log.debug("Skipping DECam edge bleed check in amp %s: no rows with usable pixels.",
                      amp.getName())
            continue
        start = enoughUsable[0]
 
        if profile[start:start + nRowsCheck].mean() <= minLowPixelsPerRow:
            log.debug("Saturated footprint reaches read edge of amp %s but no dip found.",
                      amp.getName())
            continue
 
        height = 0
        nBelow = 0
        for i in range(start, len(profile)):
            if profile[i] > minLowPixelsExtent:
                height = i + 1
                nBelow = 0
            else:
                nBelow += 1
                if nBelow >= nRowsStop:
                    break
        height += int(height*marginFraction) + 1
        height = min(height, ampBBox.getHeight())
 
        if readTop:
            maskedImage.mask[ampBBox].array[-height:, :] |= saturatedBit
        else:
            maskedImage.mask[ampBBox].array[:height, :] |= saturatedBit
        log.info("Found DECam edge bleed in amp %s at the %s read edge; masked %d rows.",
                 amp.getName(), "top" if readTop else "bottom", height)
 
 
def maskITLEdgeBleed(ccdExposure, badAmpDict,
                     fpCore, itlEdgeBleedSatMinArea=10000,
                     itlEdgeBleedSatMaxArea=100000,
                     itlEdgeBleedThreshold=5000.,
                     itlEdgeBleedModelConstant=0.02,
                     saturatedMaskName="SAT", log=None):
    """Mask edge bleeds in ITL detectors.
 
    Parameters
    ----------
    ccdExposure : `lsst.afw.image.Exposure`
        Exposure to apply masking to.
    badAmpDict : `dict` [`str`, `bool`]
        Dictionary of amplifiers, keyed by name, value is True if
        amplifier is fully masked.
    fpCore : `lsst.afw.detection._detection.Footprint`
        Footprint of saturated core.
    itlEdgeBleedThreshold : `float`, optional
        Threshold above median sky background for edge bleed detection
        (electron units).
    itlEdgeBleedModelConstant : `float`, optional
        Constant in the decaying exponential in the edge bleed masking.
    saturatedMaskName : `str`, optional
        Mask name for saturation.
    log : `logging.Logger`, optional
        Logger to handle messages.
    """
 
    log = log if log else logging.getLogger(__name__)
 
    # Get median of amplifier saturation level
    satLevel = numpy.nanmedian([ccdExposure.metadata[f"LSST ISR SATURATION LEVEL {amp.getName()}"]
                                for amp in ccdExposure.getDetector() if not badAmpDict[amp.getName()]])
 
    # 1. we check if there are several cores in the footprint:
    # Get centroid of saturated core
    xCore, yCore = fpCore.getCentroid()
    # Turn the Y detector coordinate into Y footprint coordinate
    yCoreFP = int(yCore) - fpCore.getBBox().getMinY()
    # Now test if there is one or more cores by checking if the slice at the
    # center is full of saturated pixels or has several segments of saturated
    # columns (i.e. several cores with trails)
    checkCoreNbRow = fpCore.getSpans().asArray()[yCoreFP, :]
    nbCore = 0
    indexSwitchTrue = []
    indexSwitchFalse = []
    if checkCoreNbRow[0]:
        # If the slice starts with saturated pixels
        inSatSegment = True
        nbCore = 1
        indexSwitchTrue.append(0)
    else:
        # If the slice starts with non saturated pixels
        inSatSegment = False
 
    for i, value in enumerate(checkCoreNbRow):
        if value:
            if not inSatSegment:
                indexSwitchTrue.append(i)
                # nbCore is the number of detected cores.
                nbCore += 1
                inSatSegment = True
        elif inSatSegment:
            indexSwitchFalse.append(i)
            inSatSegment = False
 
    # 1. we look for edge bleed in saturated cores in the footprint
    if nbCore == 2:
        # we now estimate the x coordinates of the edges of the subfootprint
        # for each core
        xEdgesCores = [0]
        xEdgesCores.append(int((indexSwitchTrue[1] + indexSwitchFalse[0])/2))
        xEdgesCores.append(fpCore.getSpans().asArray().shape[1])
        # Get the X and Y footprint coordinates of the cores
        for i in range(nbCore):
            subfp = fpCore.getSpans().asArray()[:, xEdgesCores[i]:xEdgesCores[i+1]]
            xCoreFP = int(xEdgesCores[i] + numpy.argmax(numpy.sum(subfp, axis=0)))
            # turn into X coordinate in detector space
            xCore = xCoreFP + fpCore.getBBox().getMinX()
            # get Y footprint coordinate of the core
            # by trimming the edges where edge bleeds are potentially dominant
            if subfp.shape[0] <= 200:
                yCoreFP = int(numpy.argmax(numpy.sum(subfp, axis=1)))
            else:
                yCoreFP = int(numpy.argmax(numpy.sum(subfp[100:-100, :],
                                                     axis=1)))
                yCoreFP = 100+yCoreFP
 
            # Estimate the width of the saturated core
            widthSat = numpy.sum(subfp[int(yCoreFP), :])
 
            subfpArea = numpy.sum(subfp)
            if subfpArea > itlEdgeBleedSatMinArea and subfpArea < itlEdgeBleedSatMaxArea:
                _applyMaskITLEdgeBleed(ccdExposure, xCore,
                                       satLevel, widthSat,
                                       itlEdgeBleedThreshold,
                                       itlEdgeBleedModelConstant,
                                       saturatedMaskName, log)
    elif nbCore > 2:
        # TODO DM-49736: support N cores in saturated footprint
        log.warning(
            "Too many (%d) cores in saturated footprint to mask edge bleeds.",
            nbCore,
        )
    else:
        # Get centroid of saturated core
        xCore, yCore = fpCore.getCentroid()
        # Turn the Y detector coordinate into Y footprint coordinate
        yCoreFP = yCore - fpCore.getBBox().getMinY()
        # Get the number of saturated columns around the centroid
        widthSat = numpy.sum(fpCore.getSpans().asArray()[int(yCoreFP), :])
        _applyMaskITLEdgeBleed(ccdExposure, xCore,
                               satLevel, widthSat, itlEdgeBleedThreshold,
                               itlEdgeBleedModelConstant, saturatedMaskName, log)
 
 
def _applyMaskITLEdgeBleed(ccdExposure, xCore,
                           satLevel, widthSat,
                           itlEdgeBleedThreshold=5000.,
                           itlEdgeBleedModelConstant=0.03,
                           saturatedMaskName="SAT", log=None):
    """Apply ITL edge bleed masking model.
 
    Parameters
    ----------
    ccdExposure : `lsst.afw.image.Exposure`
        Exposure to apply masking to.
    xCore: `int`
        X coordinate of the saturated core.
    satLevel: `float`
        Minimum saturation level of the detector.
    widthSat: `float`
        Width of the saturated core.
    itlEdgeBleedThreshold : `float`, optional
        Threshold above median sky background for edge bleed detection
        (electron units).
    itlEdgeBleedModelConstant : `float`, optional
        Constant in the decaying exponential in the edge bleed masking.
    saturatedMaskName : `str`, optional
        Mask name for saturation.
    log : `logging.Logger`, optional
        Logger to handle messages.
    """
    log = log if log else logging.getLogger(__name__)
 
    maskedImage = ccdExposure.maskedImage
    xmax = maskedImage.image.array.shape[1]
    saturatedBit = maskedImage.mask.getPlaneBitMask(saturatedMaskName)
 
    for amp in ccdExposure.getDetector():
        # Select the 2 top and bottom amplifiers around the saturated
        # core with a potential edge bleed by selecting the amplifiers
        # that have the same X coordinate as the saturated core.
        # As we don't care about the Y coordinate, we set it to the
        # center of the BBox.
        yBox = amp.getBBox().getCenter()[1]
        if amp.getBBox().contains(xCore, yBox):
 
            # Get the amp name
            ampName = amp.getName()
 
            # Because in ITLs the edge bleed happens on both edges
            # of the detector, we make a cutout around
            # both the top and bottom
            # edge bleed candidates around the saturated core.
            # We flip the cutout of the top amplifier
            # to then work with the same coordinates for both.
            # The way of selecting top vs bottom amp
            # is very specific to ITL.
            if ampName[:2] == 'C1':
                sliceImage = maskedImage.image.array[:200, :]
                sliceMask = maskedImage.mask.array[:200, :]
            elif ampName[:2] == 'C0':
                sliceImage = numpy.flipud(maskedImage.image.array[-200:, :])
                sliceMask = numpy.flipud(maskedImage.mask.array[-200:, :])
 
            # The middle columns of edge bleeds often have
            # high counts, so  we check there is an edge bleed
            # by looking at a small image up to 50 pixels from the edge
            # and around the saturated columns
            # of the saturated core, and checking its median is
            # above the sky background by itlEdgeBleedThreshold
 
            # If the centroid is too close to the edge of the detector
            # (within 5 pixels), we set the limit to the mean check
            # to the edge of the detector
            lowerRangeSmall = int(xCore)-5
            upperRangeSmall = int(xCore)+5
            if lowerRangeSmall < 0:
                lowerRangeSmall = 0
            if upperRangeSmall > xmax:
                upperRangeSmall = xmax
            ampImageBG = numpy.median(maskedImage[amp.getBBox()].image.array)
            edgeMedian = numpy.median(sliceImage[:50, lowerRangeSmall:upperRangeSmall])
            if edgeMedian > (ampImageBG + itlEdgeBleedThreshold):
 
                log.info("Found ITL edge bleed in amp %s, column %d.", ampName, xCore)
 
                # We need an estimate of the maximum width
                # of the edge bleed for our masking model
                # so we now estimate it by measuring the width of
                # areas above 60 percent of the saturation level
                # close to the edge,
                # in a cutout up to 100 pixels from the edge,
                # with a width of around the width of an amplifier.
                subImageXMin = int(xCore)-250
                subImageXMax = int(xCore)+250
                if subImageXMin < 0:
                    subImageXMin = 0
                elif subImageXMax > xmax:
                    subImageXMax = xmax
 
                subImage = sliceImage[:100, subImageXMin:subImageXMax]
                maxWidthEdgeBleed = numpy.max(numpy.sum(subImage > 0.45*satLevel,
                                                        axis=1))
 
                # Mask edge bleed with a decaying exponential model
                for y in range(200):
                    edgeBleedHalfWidth = \
                        int(((maxWidthEdgeBleed)*numpy.exp(-itlEdgeBleedModelConstant*y)
                             + widthSat)/2.)
                    lowerRange = int(xCore)-edgeBleedHalfWidth
                    upperRange = int(xCore)+edgeBleedHalfWidth
                    # If the edge bleed model goes outside the detector
                    # we set the limit for the masking
                    # to the edge of the detector
                    if lowerRange < 0:
                        lowerRange = 0
                    if upperRange > xmax:
                        upperRange = xmax
                    sliceMask[y, lowerRange:upperRange] |= saturatedBit
 
 
def maskITLSatSag(ccdExposure, fpCore, saturatedMaskName="SAT"):
    """Mask columns presenting saturation sag in saturated footprints in
    ITL detectors.
 
    Parameters
    ----------
    ccdExposure : `lsst.afw.image.Exposure`
        Exposure to apply masking to.
    fpCore : `lsst.afw.detection._detection.Footprint`
        Footprint of saturated core.
    saturatedMaskName : `str`, optional
        Mask name for saturation.
    """
 
    # TODO DM-49736: add a flux level check to apply masking
 
    maskedImage = ccdExposure.maskedImage
    saturatedBit = maskedImage.mask.getPlaneBitMask(saturatedMaskName)
 
    cc = numpy.sum(fpCore.getSpans().asArray(), axis=0)
    # Mask full columns that have 20 percent of the height of the footprint
    # saturated
    columnsToMaskFP = numpy.where(cc > fpCore.getSpans().asArray().shape[0]/5.)
 
    columnsToMask = [x + int(fpCore.getBBox().getMinX()) for x in columnsToMaskFP]
    maskedImage.mask.array[:, columnsToMask] |= saturatedBit
 
 
def maskITLDip(exposure, detectorConfig, maskPlaneNames=["SUSPECT", "ITL_DIP"], log=None):
    """Add mask bits according to the ITL dip model.
 
    Parameters
    ----------
    exposure : `lsst.afw.image.Exposure`
        Exposure to do ITL dip masking.
    detectorConfig : `lsst.ip.isr.overscanAmpConfig.OverscanDetectorConfig`
        Configuration for this detector.
    maskPlaneNames : `list [`str`], optional
        Name of the ITL Dip mask planes.
    log : `logging.Logger`, optional
        If not set, a default logger will be used.
    """
    if detectorConfig.itlDipBackgroundFraction == 0.0:
        # Nothing to do.
        return
 
    if log is None:
        log = logging.getLogger(__name__)
 
    thresh = afwDetection.Threshold(
        exposure.mask.getPlaneBitMask("SAT"),
        afwDetection.Threshold.BITMASK,
    )
    fpList = afwDetection.FootprintSet(exposure.mask, thresh).getFootprints()
 
    heights = numpy.asarray([fp.getBBox().getHeight() for fp in fpList])
 
    largeHeights, = numpy.where(heights >= detectorConfig.itlDipMinHeight)
 
    if len(largeHeights) == 0:
        return
 
    # Get the approximate image background.
    approxBackground = numpy.median(exposure.image.array)
    maskValue = exposure.mask.getPlaneBitMask(maskPlaneNames)
 
    maskBak = exposure.mask.array.copy()
    nMaskedCols = 0
 
    for index in largeHeights:
        fp = fpList[index]
        center = fp.getCentroid()
 
        nSat = numpy.sum(fp.getSpans().asArray(), axis=0)
        width = numpy.sum(nSat > detectorConfig.itlDipMinHeight)
 
        if width < detectorConfig.itlDipMinWidth:
            continue
 
        width = numpy.clip(width, None, detectorConfig.itlDipMaxWidth)
 
        dipMax = detectorConfig.itlDipBackgroundFraction * approxBackground * width
 
        # Assume sky-noise dominated; we could add in read noise here.
        if dipMax < detectorConfig.itlDipMinBackgroundNoiseFraction * numpy.sqrt(approxBackground):
            continue
 
        minCol = int(center.getX() - (detectorConfig.itlDipWidthScale * width) / 2.)
        maxCol = int(center.getX() + (detectorConfig.itlDipWidthScale * width) / 2.)
        minCol = numpy.clip(minCol, 0, None)
        maxCol = numpy.clip(maxCol, None, exposure.mask.array.shape[1] - 1)
 
        log.info(
            "Found ITL dip (width %d; bkg %.2f); masking column %d to %d.",
            width,
            approxBackground,
            minCol,
            maxCol,
        )
 
        exposure.mask.array[:, minCol: maxCol + 1] |= maskValue
 
        nMaskedCols += (maxCol - minCol + 1)
 
    if nMaskedCols > detectorConfig.itlDipMaxColsPerImage:
        log.warning(
            "Too many (%d) columns would be masked on this image from dip masking; restoring original mask.",
            nMaskedCols,
        )
        exposure.mask.array[:, :] = maskBak
 
 
def interpolateFromMask(maskedImage, fwhm, growSaturatedFootprints=1,
                        maskNameList=['SAT'], fallbackValue=None, useLegacyInterp=True):
    """Interpolate over defects identified by a particular set of mask planes.
 
    Parameters
    ----------
    maskedImage : `lsst.afw.image.MaskedImage`
        Image to process.
    fwhm : `float`
        FWHM of double Gaussian smoothing kernel.
    growSaturatedFootprints : scalar, optional
        Number of pixels to grow footprints for saturated pixels.
    maskNameList : `List` of `str`, optional
        Mask plane name.
    fallbackValue : scalar, optional
        Value of last resort for interpolation.
 
    Notes
    -----
    The ``fwhm`` parameter is used to create a PSF, but the underlying
    interpolation code (`lsst.meas.algorithms.interpolateOverDefects`) does
    not currently make use of this information.
    """
    mask = maskedImage.getMask()
 
    if growSaturatedFootprints > 0 and "SAT" in maskNameList:
        # If we are interpolating over an area larger than the original masked
        # region, we need to expand the original mask bit to the full area to
        # explain why we interpolated there.
        growMasks(mask, radius=growSaturatedFootprints, maskNameList=['SAT'], maskValue="SAT")
 
    thresh = afwDetection.Threshold(mask.getPlaneBitMask(maskNameList), afwDetection.Threshold.BITMASK)
    fpSet = afwDetection.FootprintSet(mask, thresh)
    defectList = Defects.fromFootprintList(fpSet.getFootprints())
 
    interpolateDefectList(maskedImage, defectList, fwhm, fallbackValue=fallbackValue,
                          maskNameList=maskNameList, useLegacyInterp=useLegacyInterp)
 
    return maskedImage
 
 
def saturationCorrection(maskedImage, saturation, fwhm, growFootprints=1, interpolate=True, maskName='SAT',
                         fallbackValue=None, useLegacyInterp=True):
    """Mark saturated pixels and optionally interpolate over them
 
    Parameters
    ----------
    maskedImage : `lsst.afw.image.MaskedImage`
        Image to process.
    saturation  : scalar
        Saturation level used as the detection threshold.
    fwhm : `float`
        FWHM of double Gaussian smoothing kernel.
    growFootprints : scalar, optional
        Number of pixels to grow footprints of detected regions.
    interpolate : Bool, optional
        If True, saturated pixels are interpolated over.
    maskName : str, optional
        Mask plane name.
    fallbackValue : scalar, optional
        Value of last resort for interpolation.
 
    Notes
    -----
    The ``fwhm`` parameter is used to create a PSF, but the underlying
    interpolation code (`lsst.meas.algorithms.interpolateOverDefects`) does
    not currently make use of this information.
    """
    defectList = makeThresholdMask(
        maskedImage=maskedImage,
        threshold=saturation,
        growFootprints=growFootprints,
        maskName=maskName,
    )
    if interpolate:
        interpolateDefectList(maskedImage, defectList, fwhm, fallbackValue=fallbackValue,
                              maskNameList=[maskName], useLegacyInterp=useLegacyInterp)
 
    return maskedImage
 
 
def trimToMatchCalibBBox(rawMaskedImage, calibMaskedImage):
    """Compute number of edge trim pixels to match the calibration data.
 
    Use the dimension difference between the raw exposure and the
    calibration exposure to compute the edge trim pixels.  This trim
    is applied symmetrically, with the same number of pixels masked on
    each side.
 
    Parameters
    ----------
    rawMaskedImage : `lsst.afw.image.MaskedImage`
        Image to trim.
    calibMaskedImage : `lsst.afw.image.MaskedImage`
        Calibration image to draw new bounding box from.
 
    Returns
    -------
    replacementMaskedImage : `lsst.afw.image.MaskedImage`
        ``rawMaskedImage`` trimmed to the appropriate size.
 
    Raises
    ------
    RuntimeError
       Raised if ``rawMaskedImage`` cannot be symmetrically trimmed to
       match ``calibMaskedImage``.
    """
    nx, ny = rawMaskedImage.getBBox().getDimensions() - calibMaskedImage.getBBox().getDimensions()
    if nx != ny:
        raise RuntimeError("Raw and calib maskedImages are trimmed differently in X and Y.")
    if nx % 2 != 0:
        raise RuntimeError("Calibration maskedImage is trimmed unevenly in X.")
    if nx < 0:
        raise RuntimeError("Calibration maskedImage is larger than raw data.")
 
    nEdge = nx//2
    if nEdge > 0:
        replacementMaskedImage = rawMaskedImage[nEdge:-nEdge, nEdge:-nEdge, afwImage.LOCAL]
        SourceDetectionTask.setEdgeBits(
            rawMaskedImage,
            replacementMaskedImage.getBBox(),
            rawMaskedImage.getMask().getPlaneBitMask("EDGE")
        )
    else:
        replacementMaskedImage = rawMaskedImage
 
    return replacementMaskedImage
 
 
def biasCorrection(maskedImage, biasMaskedImage, trimToFit=False):
    """Apply bias correction in place.
 
    Parameters
    ----------
    maskedImage : `lsst.afw.image.MaskedImage`
       Image to process.  The image is modified by this method.
    biasMaskedImage : `lsst.afw.image.MaskedImage`
        Bias image of the same size as ``maskedImage``
    trimToFit : `Bool`, optional
        If True, raw data is symmetrically trimmed to match
        calibration size.
 
    Raises
    ------
    RuntimeError
        Raised if ``maskedImage`` and ``biasMaskedImage`` do not have
        the same size.
 
    """
    if trimToFit:
        maskedImage = trimToMatchCalibBBox(maskedImage, biasMaskedImage)
 
    if maskedImage.getBBox(afwImage.LOCAL) != biasMaskedImage.getBBox(afwImage.LOCAL):
        raise RuntimeError("maskedImage bbox %s != biasMaskedImage bbox %s" %
                           (maskedImage.getBBox(afwImage.LOCAL), biasMaskedImage.getBBox(afwImage.LOCAL)))
    maskedImage -= biasMaskedImage
 
 
def darkCorrection(maskedImage, darkMaskedImage, expScale, darkScale, invert=False, trimToFit=False):
    """Apply dark correction in place.
 
    Parameters
    ----------
    maskedImage : `lsst.afw.image.MaskedImage`
       Image to process.  The image is modified by this method.
    darkMaskedImage : `lsst.afw.image.MaskedImage`
        Dark image of the same size as ``maskedImage``.
    expScale : scalar
        Dark exposure time for ``maskedImage``.
    darkScale : scalar
        Dark exposure time for ``darkMaskedImage``.
    invert : `Bool`, optional
        If True, re-add the dark to an already corrected image.
    trimToFit : `Bool`, optional
        If True, raw data is symmetrically trimmed to match
        calibration size.
 
    Raises
    ------
    RuntimeError
        Raised if ``maskedImage`` and ``darkMaskedImage`` do not have
        the same size.
 
    Notes
    -----
    The dark correction is applied by calculating:
        maskedImage -= dark * expScaling / darkScaling
    """
    if trimToFit:
        maskedImage = trimToMatchCalibBBox(maskedImage, darkMaskedImage)
 
    if maskedImage.getBBox(afwImage.LOCAL) != darkMaskedImage.getBBox(afwImage.LOCAL):
        raise RuntimeError("maskedImage bbox %s != darkMaskedImage bbox %s" %
                           (maskedImage.getBBox(afwImage.LOCAL), darkMaskedImage.getBBox(afwImage.LOCAL)))
 
    scale = expScale / darkScale
    if not invert:
        maskedImage.scaledMinus(scale, darkMaskedImage)
    else:
        maskedImage.scaledPlus(scale, darkMaskedImage)
 
 
def updateVariance(maskedImage, gain, readNoise, replace=True):
    """Set the variance plane based on the image plane.
 
    The maskedImage must have units of `adu` (if gain != 1.0) or
    electron (if gain == 1.0). This routine will always produce a
    variance plane in the same units as the image.
 
    Parameters
    ----------
    maskedImage : `lsst.afw.image.MaskedImage`
        Image to process.  The variance plane is modified.
    gain : scalar
        The amplifier gain in electron/adu.
    readNoise : scalar
        The amplifier read noise in electron/pixel.
    replace : `bool`, optional
        Replace the current variance?  If False, the image
        variance will be added to the current variance plane.
    """
    var = maskedImage.variance
    if replace:
        var[:, :] = maskedImage.image
    else:
        var[:, :] += maskedImage.image
    var /= gain
    var += (readNoise/gain)**2
 
 
def flatCorrection(maskedImage, flatMaskedImage, scalingType, userScale=1.0, invert=False, trimToFit=False):
    """Apply flat correction in place.
 
    Parameters
    ----------
    maskedImage : `lsst.afw.image.MaskedImage`
        Image to process.  The image is modified.
    flatMaskedImage : `lsst.afw.image.MaskedImage`
        Flat image of the same size as ``maskedImage``
    scalingType : str
        Flat scale computation method.  Allowed values are 'MEAN',
        'MEDIAN', or 'USER'.
    userScale : scalar, optional
        Scale to use if ``scalingType='USER'``.
    invert : `Bool`, optional
        If True, unflatten an already flattened image.
    trimToFit : `Bool`, optional
        If True, raw data is symmetrically trimmed to match
        calibration size.
 
    Raises
    ------
    RuntimeError
        Raised if ``maskedImage`` and ``flatMaskedImage`` do not have
        the same size or if ``scalingType`` is not an allowed value.
    """
    if trimToFit:
        maskedImage = trimToMatchCalibBBox(maskedImage, flatMaskedImage)
 
    if maskedImage.getBBox(afwImage.LOCAL) != flatMaskedImage.getBBox(afwImage.LOCAL):
        raise RuntimeError("maskedImage bbox %s != flatMaskedImage bbox %s" %
                           (maskedImage.getBBox(afwImage.LOCAL), flatMaskedImage.getBBox(afwImage.LOCAL)))
 
    # Figure out scale from the data
    # Ideally the flats are normalized by the calibration product pipeline,
    # but this allows some flexibility in the case that the flat is created by
    # some other mechanism.
    if scalingType in ('MEAN', 'MEDIAN'):
        scalingType = afwMath.stringToStatisticsProperty(scalingType)
        flatScale = afwMath.makeStatistics(flatMaskedImage.image, scalingType).getValue()
    elif scalingType == 'USER':
        flatScale = userScale
    else:
        raise RuntimeError('%s : %s not implemented' % ("flatCorrection", scalingType))
 
    if not invert:
        maskedImage.scaledDivides(1.0/flatScale, flatMaskedImage)
    else:
        maskedImage.scaledMultiplies(1.0/flatScale, flatMaskedImage)
 
 
def illuminationCorrection(maskedImage, illumMaskedImage, illumScale, trimToFit=True):
    """Apply illumination correction in place.
 
    Parameters
    ----------
    maskedImage : `lsst.afw.image.MaskedImage`
        Image to process.  The image is modified.
    illumMaskedImage : `lsst.afw.image.MaskedImage`
        Illumination correction image of the same size as ``maskedImage``.
    illumScale : scalar
        Scale factor for the illumination correction.
    trimToFit : `Bool`, optional
        If True, raw data is symmetrically trimmed to match
        calibration size.
 
    Raises
    ------
    RuntimeError
        Raised if ``maskedImage`` and ``illumMaskedImage`` do not have
        the same size.
    """
    if trimToFit:
        maskedImage = trimToMatchCalibBBox(maskedImage, illumMaskedImage)
 
    if maskedImage.getBBox(afwImage.LOCAL) != illumMaskedImage.getBBox(afwImage.LOCAL):
        raise RuntimeError("maskedImage bbox %s != illumMaskedImage bbox %s" %
                           (maskedImage.getBBox(afwImage.LOCAL), illumMaskedImage.getBBox(afwImage.LOCAL)))
 
    maskedImage.scaledDivides(1.0/illumScale, illumMaskedImage)
 
 
@contextmanager
def gainContext(exp, image, apply, gains=None, invert=False, isTrimmed=True):
    """Context manager that applies and removes gain.
 
    Parameters
    ----------
    exp : `lsst.afw.image.Exposure`
        Exposure to apply/remove gain.
    image : `lsst.afw.image.Image`
        Image to apply/remove gain.
    apply : `bool`
        If True, apply and remove the amplifier gain.
    gains : `dict` [`str`, `float`], optional
        A dictionary, keyed by amplifier name, of the gains to use.
        If gains is None, the nominal gains in the amplifier object are used.
    invert : `bool`, optional
        Invert the gains (e.g. convert electrons to adu temporarily)?
    isTrimmed : `bool`, optional
        Is this a trimmed exposure?
 
    Yields
    ------
    exp : `lsst.afw.image.Exposure`
        Exposure with the gain applied.
    """
    # check we have all of them if provided because mixing and matching would
    # be a real mess
    if gains and apply is True:
        ampNames = [amp.getName() for amp in exp.getDetector()]
        for ampName in ampNames:
            if ampName not in gains.keys():
                raise RuntimeError(f"Gains provided to gain context, but no entry found for amp {ampName}")
 
    if apply:
        ccd = exp.getDetector()
        for amp in ccd:
            sim = image.Factory(image, amp.getBBox() if isTrimmed else amp.getRawBBox())
            if gains:
                gain = gains[amp.getName()]
            else:
                gain = amp.getGain()
            if invert:
                sim /= gain
            else:
                sim *= gain
 
    try:
        yield exp
    finally:
        if apply:
            ccd = exp.getDetector()
            for amp in ccd:
                sim = image.Factory(image, amp.getBBox() if isTrimmed else amp.getRawBBox())
                if gains:
                    gain = gains[amp.getName()]
                else:
                    gain = amp.getGain()
                if invert:
                    sim *= gain
                else:
                    sim /= gain
 
 
def attachTransmissionCurve(exposure, opticsTransmission=None, filterTransmission=None,
                            sensorTransmission=None, atmosphereTransmission=None):
    """Attach a TransmissionCurve to an Exposure, given separate curves for
    different components.
 
    Parameters
    ----------
    exposure : `lsst.afw.image.Exposure`
        Exposure object to modify by attaching the product of all given
        ``TransmissionCurves`` in post-assembly trimmed detector coordinates.
        Must have a valid ``Detector`` attached that matches the detector
        associated with sensorTransmission.
    opticsTransmission : `lsst.afw.image.TransmissionCurve`
        A ``TransmissionCurve`` that represents the throughput of the optics,
        to be evaluated in focal-plane coordinates.
    filterTransmission : `lsst.afw.image.TransmissionCurve`
        A ``TransmissionCurve`` that represents the throughput of the filter
        itself, to be evaluated in focal-plane coordinates.
    sensorTransmission : `lsst.afw.image.TransmissionCurve`
        A ``TransmissionCurve`` that represents the throughput of the sensor
        itself, to be evaluated in post-assembly trimmed detector coordinates.
    atmosphereTransmission : `lsst.afw.image.TransmissionCurve`
        A ``TransmissionCurve`` that represents the throughput of the
        atmosphere, assumed to be spatially constant.
 
    Returns
    -------
    combined : `lsst.afw.image.TransmissionCurve`
        The TransmissionCurve attached to the exposure.
 
    Notes
    -----
    All ``TransmissionCurve`` arguments are optional; if none are provided, the
    attached ``TransmissionCurve`` will have unit transmission everywhere.
    """
    combined = afwImage.TransmissionCurve.makeIdentity()
    if atmosphereTransmission is not None:
        combined *= atmosphereTransmission
    if opticsTransmission is not None:
        combined *= opticsTransmission
    if filterTransmission is not None:
        combined *= filterTransmission
    detector = exposure.getDetector()
    fpToPix = detector.getTransform(fromSys=camGeom.FOCAL_PLANE,
                                    toSys=camGeom.PIXELS)
    combined = combined.transformedBy(fpToPix)
    if sensorTransmission is not None:
        combined *= sensorTransmission
    exposure.getInfo().setTransmissionCurve(combined)
    return combined
 
 
def applyGains(exposure, normalizeGains=False, ptcGains=None, isTrimmed=True):
    """Scale an exposure by the amplifier gains.
 
    Parameters
    ----------
    exposure : `lsst.afw.image.Exposure`
        Exposure to process.  The image is modified.
    normalizeGains : `Bool`, optional
        If True, then amplifiers are scaled to force the median of
        each amplifier to equal the median of those medians.
    ptcGains : `dict`[`str`], optional
        Dictionary keyed by amp name containing the PTC gains.
    isTrimmed : `bool`, optional
        Is the input image trimmed?
    """
    ccd = exposure.getDetector()
    ccdImage = exposure.getMaskedImage()
 
    medians = []
    for amp in ccd:
        if isTrimmed:
            sim = ccdImage.Factory(ccdImage, amp.getBBox())
        else:
            sim = ccdImage.Factory(ccdImage, amp.getRawBBox())
        if ptcGains:
            sim *= ptcGains[amp.getName()]
        else:
            sim *= amp.getGain()
 
        if normalizeGains:
            medians.append(numpy.median(sim.getImage().getArray()))
 
    if normalizeGains:
        median = numpy.median(numpy.array(medians))
        for index, amp in enumerate(ccd):
            if isTrimmed:
                sim = ccdImage.Factory(ccdImage, amp.getBBox())
            else:
                sim = ccdImage.Factory(ccdImage, amp.getRawBBox())
            if medians[index] != 0.0:
                sim *= median/medians[index]
 
 
def widenSaturationTrails(mask):
    """Grow the saturation trails by an amount dependent on the width of the
    trail.
 
    Parameters
    ----------
    mask : `lsst.afw.image.Mask`
        Mask which will have the saturated areas grown.
    """
 
    extraGrowDict = {}
    for i in range(1, 6):
        extraGrowDict[i] = 0
    for i in range(6, 8):
        extraGrowDict[i] = 1
    for i in range(8, 10):
        extraGrowDict[i] = 3
    extraGrowMax = 4
 
    if extraGrowMax <= 0:
        return
 
    saturatedBit = mask.getPlaneBitMask("SAT")
 
    xmin, ymin = mask.getBBox().getMin()
    width = mask.getWidth()
 
    thresh = afwDetection.Threshold(saturatedBit, afwDetection.Threshold.BITMASK)
    fpList = afwDetection.FootprintSet(mask, thresh).getFootprints()
 
    for fp in fpList:
        for s in fp.getSpans():
            x0, x1 = s.getX0(), s.getX1()
 
            extraGrow = extraGrowDict.get(x1 - x0 + 1, extraGrowMax)
            if extraGrow > 0:
                y = s.getY() - ymin
                x0 -= xmin + extraGrow
                x1 -= xmin - extraGrow
 
                if x0 < 0:
                    x0 = 0
                if x1 >= width - 1:
                    x1 = width - 1
 
                mask.array[y, x0:x1+1] |= saturatedBit
 
 
def setBadRegions(exposure, badStatistic="MEDIAN"):
    """Set all BAD areas of the chip to the average of the rest of the exposure
 
    Parameters
    ----------
    exposure : `lsst.afw.image.Exposure`
        Exposure to mask.  The exposure mask is modified.
    badStatistic : `str`, optional
        Statistic to use to generate the replacement value from the
        image data.  Allowed values are 'MEDIAN' or 'MEANCLIP'.
 
    Returns
    -------
    badPixelCount : scalar
        Number of bad pixels masked.
    badPixelValue : scalar
        Value substituted for bad pixels.
 
    Raises
    ------
    RuntimeError
        Raised if `badStatistic` is not an allowed value.
    """
    if badStatistic == "MEDIAN":
        statistic = afwMath.MEDIAN
    elif badStatistic == "MEANCLIP":
        statistic = afwMath.MEANCLIP
    else:
        raise RuntimeError("Impossible method %s of bad region correction" % badStatistic)
 
    mi = exposure.getMaskedImage()
    mask = mi.getMask()
    BAD = mask.getPlaneBitMask("BAD")
    INTRP = mask.getPlaneBitMask("INTRP")
 
    sctrl = afwMath.StatisticsControl()
    sctrl.setAndMask(BAD)
    value = afwMath.makeStatistics(mi, statistic, sctrl).getValue()
 
    maskArray = mask.getArray()
    imageArray = mi.getImage().getArray()
    badPixels = numpy.logical_and((maskArray & BAD) > 0, (maskArray & INTRP) == 0)
    imageArray[:] = numpy.where(badPixels, value, imageArray)
 
    return badPixels.sum(), value
 
 
def checkFilter(exposure, filterList, log):
    """Check to see if an exposure is in a filter specified by a list.
 
    The goal of this is to provide a unified filter checking interface
    for all filter dependent stages.
 
    Parameters
    ----------
    exposure : `lsst.afw.image.Exposure`
        Exposure to examine.
    filterList : `list` [`str`]
        List of physical_filter names to check.
    log : `logging.Logger`
        Logger to handle messages.
 
    Returns
    -------
    result : `bool`
        True if the exposure's filter is contained in the list.
    """
    if len(filterList) == 0:
        return False
    thisFilter = exposure.getFilter()
    if thisFilter is None:
        log.warning("No FilterLabel attached to this exposure!")
        return False
 
    thisPhysicalFilter = getPhysicalFilter(thisFilter, log)
    if thisPhysicalFilter in filterList:
        return True
    elif thisFilter.bandLabel in filterList:
        if log:
            log.warning("Physical filter (%s) should be used instead of band %s for filter configurations"
                        " (%s)", thisPhysicalFilter, thisFilter.bandLabel, filterList)
        return True
    else:
        return False
 
 
def getPhysicalFilter(filterLabel, log):
    """Get the physical filter label associated with the given filterLabel.
 
    If ``filterLabel`` is `None` or there is no physicalLabel attribute
    associated with the given ``filterLabel``, the returned label will be
    "Unknown".
 
    Parameters
    ----------
    filterLabel : `lsst.afw.image.FilterLabel`
        The `lsst.afw.image.FilterLabel` object from which to derive the
        physical filter label.
    log : `logging.Logger`
        Logger to handle messages.
 
    Returns
    -------
    physicalFilter : `str`
        The value returned by the physicalLabel attribute of ``filterLabel`` if
        it exists, otherwise set to \"Unknown\".
    """
    if filterLabel is None:
        physicalFilter = "Unknown"
        log.warning("filterLabel is None.  Setting physicalFilter to \"Unknown\".")
    else:
        try:
            physicalFilter = filterLabel.physicalLabel
        except RuntimeError:
            log.warning("filterLabel has no physicalLabel attribute.  Setting physicalFilter to \"Unknown\".")
            physicalFilter = "Unknown"
    return physicalFilter
 
 
def countMaskedPixels(maskedIm, maskPlane):
    """Count the number of pixels in a given mask plane.
 
    Parameters
    ----------
    maskedIm : `~lsst.afw.image.MaskedImage`
        Masked image to examine.
    maskPlane : `str`
        Name of the mask plane to examine.
 
    Returns
    -------
    nPix : `int`
        Number of pixels in the requested mask plane.
    """
    maskBit = maskedIm.mask.getPlaneBitMask(maskPlane)
    nPix = numpy.where(numpy.bitwise_and(maskedIm.mask.array, maskBit))[0].flatten().size
    return nPix
 
 
def getExposureGains(exposure):
    """Get the per-amplifier gains used for this exposure.
 
    Parameters
    ----------
    exposure : `lsst.afw.image.Exposure`
        The exposure to find gains for.
 
    Returns
    -------
    gains : `dict` [`str` `float`]
        Dictionary of gain values, keyed by amplifier name.
        Returns empty dict when detector is None.
    """
    det = exposure.getDetector()
    if det is None:
        return dict()
 
    metadata = exposure.getMetadata()
    gains = {}
    for amp in det:
        ampName = amp.getName()
        # The key may use the new LSST ISR or the old LSST prefix
        if (key1 := f"LSST ISR GAIN {ampName}") in metadata:
            gains[ampName] = metadata[key1]
        elif (key2 := f"LSST GAIN {ampName}") in metadata:
            gains[ampName] = metadata[key2]
        else:
            gains[ampName] = amp.getGain()
    return gains
 
 
def getExposureReadNoises(exposure):
    """Get the per-amplifier read noise used for this exposure.
 
    Parameters
    ----------
    exposure : `lsst.afw.image.Exposure`
        The exposure to find read noise for.
 
    Returns
    -------
    readnoises : `dict` [`str` `float`]
        Dictionary of read noise values, keyed by amplifier name.
        Returns empty dict when detector is None.
    """
    det = exposure.getDetector()
    if det is None:
        return dict()
 
    metadata = exposure.getMetadata()
    readnoises = {}
    for amp in det:
        ampName = amp.getName()
        # The key may use the new LSST ISR or the old LSST prefix
        if (key1 := f"LSST ISR READNOISE {ampName}") in metadata:
            readnoises[ampName] = metadata[key1]
        elif (key2 := f"LSST READNOISE {ampName}") in metadata:
            readnoises[ampName] = metadata[key2]
        else:
            readnoises[ampName] = amp.getReadNoise()
    return readnoises
 
 
def isTrimmedExposure(exposure):
    """Check if the unused pixels (pre-/over-scan pixels) have
    been trimmed from an exposure.
 
    Parameters
    ----------
    exposure : `lsst.afw.image.Exposure`
        The exposure to check.
 
    Returns
    -------
    result : `bool`
        True if the image is trimmed, else False.
    """
    return exposure.getDetector().getBBox() == exposure.getBBox()
 
 
def isTrimmedImage(image, detector):
    """Check if the unused pixels (pre-/over-scan pixels) have
    been trimmed from an image
 
    Parameters
    ----------
    image : `lsst.afw.image.Image`
        The image to check.
    detector : `lsst.afw.cameraGeom.Detector`
        The detector associated with the image.
 
    Returns
    -------
    result : `bool`
        True if the image is trimmed, else False.
    """
    return detector.getBBox() == image.getBBox()
 
 
def compareCameraKeywords(
    doRaiseOnCalibMismatch,
    cameraKeywordsToCompare,
    exposureMetadata,
    calib,
    calibName,
    log=None,
):
    """Compare header keywords to confirm camera states match.
 
    Parameters
    ----------
    doRaiseOnCalibMismatch : `bool`
        Raise on calibration mismatch?  Otherwise, log a warning.
    cameraKeywordsToCompare : `list` [`str`]
        List of camera keywords to compare.
    exposureMetadata : `lsst.daf.base.PropertyList`
        Header for the exposure being processed.
    calib : `lsst.afw.image.Exposure` or `lsst.ip.isr.IsrCalib`
        Calibration to be applied.
    calibName : `str`
        Calib type for log message.
    log : `logging.Logger`, optional
        Logger to handle messages.
    """
    try:
        calibMetadata = calib.metadata
    except AttributeError:
        return
 
    log = log if log else logging.getLogger(__name__)
 
    missingKeywords = []
    for keyword in cameraKeywordsToCompare:
        exposureValue = exposureMetadata.get(keyword, None)
        if exposureValue is None:
            log.debug("Sequencer keyword %s not found in exposure metadata.", keyword)
            continue
 
        calibValue = calibMetadata.get(keyword, None)
 
        # We don't log here if there is a missing keyword.
        if calibValue is None:
            missingKeywords.append(keyword)
            continue
 
        if exposureValue != calibValue:
            if doRaiseOnCalibMismatch:
                raise RuntimeError(
                    "Sequencer mismatch for %s [%s]: exposure: %s calib: %s",
                    calibName,
                    keyword,
                    exposureValue,
                    calibValue,
                )
            else:
                log.warning(
                    "Sequencer mismatch for %s [%s]: exposure: %s calib: %s",
                    calibName,
                    keyword,
                    exposureValue,
                    calibValue,
                )
                exposureMetadata[f"ISR {calibName.upper()} SEQUENCER MISMATCH"] = True
 
    if missingKeywords:
        log.info(
            "Calibration %s missing keywords %s, which were not checked.",
            calibName,
            ",".join(missingKeywords),
        )
 
 
def symmetrize(inputArray):
    """ Copy array over 4 quadrants prior to convolution.
 
    Parameters
    ----------
    inputarray : `numpy.array`
        Input array to symmetrize.
 
    Returns
    -------
    aSym : `numpy.array`
        Symmetrized array.
    """
    targetShape = list(inputArray.shape)
    r1, r2 = inputArray.shape[-1], inputArray.shape[-2]
    targetShape[-1] = 2*r1-1
    targetShape[-2] = 2*r2-1
    aSym = numpy.ndarray(tuple(targetShape))
    aSym[..., r2-1:, r1-1:] = inputArray
    aSym[..., r2-1:, r1-1::-1] = inputArray
    aSym[..., r2-1::-1, r1-1::-1] = inputArray
    aSym[..., r2-1::-1, r1-1:] = inputArray
 
    return aSym