crosstalk.py

Go to the documentation of this file.

#
# LSST Data Management System
# Copyright 2008-2017 AURA/LSST.
#
# 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 <https://www.lsstcorp.org/LegalNotices/>.
#
"""
Apply intra-detector crosstalk corrections
"""
import numpy as np
 
import lsst.afw.math
import lsst.afw.table
import lsst.afw.detection
from lsst.pex.config import Config, Field, ChoiceField, ListField
from lsst.pipe.base import Task
 
__all__ = ["CrosstalkConfig", "CrosstalkTask", "subtractCrosstalk", "writeCrosstalkCoeffs",
           "NullCrosstalkTask"]
 
 
class CrosstalkConfig(Config):
    """Configuration for intra-detector crosstalk removal."""
    minPixelToMask = Field(
        dtype=float,
        doc="Set crosstalk mask plane for pixels over this value.",
        default=45000
    )
    crosstalkMaskPlane = Field(
        dtype=str,
        doc="Name for crosstalk mask plane.",
        default="CROSSTALK"
    )
    crosstalkBackgroundMethod = ChoiceField(
        dtype=str,
        doc="Type of background subtraction to use when applying correction.",
        default="None",
        allowed={
            "None": "Do no background subtraction.",
            "AMP": "Subtract amplifier-by-amplifier background levels.",
            "DETECTOR": "Subtract detector level background."
        },
    )
    useConfigCoefficients = Field(
        dtype=bool,
        doc="Ignore the detector crosstalk information in favor of CrosstalkConfig values?",
        default=False,
    )
    crosstalkValues = ListField(
        dtype=float,
        doc=("Amplifier-indexed crosstalk coefficients to use.  This should be arranged as a 1 x nAmp**2 "
             "list of coefficients, such that when reshaped by crosstalkShape, the result is nAmp x nAmp. "
             "This matrix should be structured so CT * [amp0 amp1 amp2 ...]^T returns the column "
             "vector [corr0 corr1 corr2 ...]^T."),
        default=[0.0],
    )
    crosstalkShape = ListField(
        dtype=int,
        doc="Shape of the coefficient array.  This should be equal to [nAmp, nAmp].",
        default=[1],
    )
 
    def getCrosstalk(self, detector=None):
        """Return a 2-D numpy array of crosstalk coefficients in the proper shape.
 
        Parameters
        ----------
        detector : `lsst.afw.cameraGeom.detector`
            Detector that is to be crosstalk corrected.
 
        Returns
        -------
        coeffs : `numpy.ndarray`
            Crosstalk coefficients that can be used to correct the detector.
 
        Raises
        ------
        RuntimeError
            Raised if no coefficients could be generated from this detector/configuration.
        """
        if self.useConfigCoefficients is True:
            coeffs = np.array(self.crosstalkValues).reshape(self.crosstalkShape)
            if detector is not None:
                nAmp = len(detector)
                if coeffs.shape != (nAmp, nAmp):
                    raise RuntimeError("Constructed crosstalk coeffients do not match detector shape. " +
                                       f"{coeffs.shape} {nAmp}")
            return coeffs
        elif detector is not None and detector.hasCrosstalk() is True:
            # Assume the detector defines itself consistently.
            return detector.getCrosstalk()
        else:
            raise RuntimeError("Attempted to correct crosstalk without crosstalk coefficients")
 
    def hasCrosstalk(self, detector=None):
        """Return a boolean indicating if crosstalk coefficients exist.
 
        Parameters
        ----------
        detector : `lsst.afw.cameraGeom.detector`
            Detector that is to be crosstalk corrected.
 
        Returns
        -------
        hasCrosstalk : `bool`
            True if this detector/configuration has crosstalk coefficients defined.
        """
        if self.useConfigCoefficients is True and self.crosstalkValues is not None:
            return True
        elif detector is not None and detector.hasCrosstalk() is True:
            return True
        else:
            return False
 
 
class CrosstalkTask(Task):
    """Apply intra-detector crosstalk correction."""
    ConfigClass = CrosstalkConfig
    _DefaultName = 'isrCrosstalk'
 
    def prepCrosstalk(self, dataRef):
        """Placeholder for crosstalk preparation method, e.g., for inter-detector crosstalk.
 
        Parameters
        ----------
        dataRef : `daf.persistence.butlerSubset.ButlerDataRef`
            Butler reference of the detector data to be processed.
 
        See also
        --------
        lsst.obs.decam.crosstalk.DecamCrosstalkTask.prepCrosstalk
        """
        return
 
    def run(self, exposure, crosstalkSources=None, isTrimmed=False):
        """Apply intra-detector crosstalk correction
 
        Parameters
        ----------
        exposure : `lsst.afw.image.Exposure`
            Exposure for which to remove crosstalk.
        crosstalkSources : `defaultdict`, optional
            Image data and crosstalk coefficients from other detectors/amps that are
            sources of crosstalk in exposure.
            The default for intra-detector crosstalk here is None.
        isTrimmed : `bool`
            The image is already trimmed.
            This should no longer be needed once DM-15409 is resolved.
 
        Raises
        ------
        RuntimeError
            Raised if called for a detector that does not have a
            crosstalk correction.
        TypeError
            Raised if crosstalkSources is not None
            and not a numpy array or a dictionary.
        """
        if crosstalkSources is not None:
            if isinstance(crosstalkSources, np.ndarray):
                coeffs = crosstalkSources
            elif isinstance(crosstalkSources, dict):
                # Nested dictionary produced by `measureCrosstalk.py`
                # There are two keys first in the nested dictionary
                for fKey, fValue in crosstalkSources.items():
                    for sKey, sValue in fValue.items():
                        firstKey = fKey
                        secondKey = sKey
                tempDict = crosstalkSources[firstKey][secondKey]
                coeffs = []
                for thirdKey in tempDict:
                    tempList = []
                    for fourthKey in tempDict[thirdKey]:
                        value = tempDict[thirdKey][fourthKey]
                        tempList.append(value)
                    coeffs.append(tempList)
                coeffs = np.array(coeffs)
            else:
                raise TypeError("crosstalkSources not of the correct type: `np.array` or `dict`")
        else:
            detector = exposure.getDetector()
            if not self.config.hasCrosstalk(detector=detector):
                raise RuntimeError("Attempted to correct crosstalk without crosstalk coefficients")
            coeffs = self.config.getCrosstalk(detector=detector)
 
        self.log.info("Applying crosstalk correction.")
        subtractCrosstalk(exposure, crosstalkCoeffs=coeffs,
                          minPixelToMask=self.config.minPixelToMask,
                          crosstalkStr=self.config.crosstalkMaskPlane, isTrimmed=isTrimmed,
                          backgroundMethod=self.config.crosstalkBackgroundMethod)
 
 
# Flips required to get the corner to the lower-left
# (an arbitrary choice; flips are relative, so the choice of reference here is not important)
X_FLIP = {lsst.afw.cameraGeom.ReadoutCorner.LL: False, lsst.afw.cameraGeom.ReadoutCorner.LR: True,
          lsst.afw.cameraGeom.ReadoutCorner.UL: False, lsst.afw.cameraGeom.ReadoutCorner.UR: True}
Y_FLIP = {lsst.afw.cameraGeom.ReadoutCorner.LL: False, lsst.afw.cameraGeom.ReadoutCorner.LR: False,
          lsst.afw.cameraGeom.ReadoutCorner.UL: True, lsst.afw.cameraGeom.ReadoutCorner.UR: True}
 
 
class NullCrosstalkTask(CrosstalkTask):
    def run(self, exposure, crosstalkSources=None):
        self.log.info("Not performing any crosstalk correction")
 
 
def extractAmp(image, amp, corner, isTrimmed=False):
    """Return an image of the amp
 
    The returned image will have the amp's readout corner in the
    nominated `corner`.
 
    Parameters
    ----------
    image : `lsst.afw.image.Image` or `lsst.afw.image.MaskedImage`
        Image containing the amplifier of interest.
    amp : `lsst.afw.table.AmpInfoRecord`
        Amplifier information.
    corner : `lsst.afw.table.ReadoutCorner` or `None`
        Corner in which to put the amp's readout corner, or `None` for
        no flipping.
    isTrimmed : `bool`
        The image is already trimmed.
        This should no longer be needed once DM-15409 is resolved.
 
    Returns
    -------
    output : `lsst.afw.image.Image`
        Image of the amplifier in the standard configuration.
    """
    output = image[amp.getBBox() if isTrimmed else amp.getRawDataBBox()]
    ampCorner = amp.getReadoutCorner()
    # Flipping is necessary only if the desired configuration doesn't match what we currently have
    xFlip = X_FLIP[corner] ^ X_FLIP[ampCorner]
    yFlip = Y_FLIP[corner] ^ Y_FLIP[ampCorner]
    return lsst.afw.math.flipImage(output, xFlip, yFlip)
 
 
def calculateBackground(mi, badPixels=["BAD"]):
    """Calculate median background in image
 
    Getting a great background model isn't important for crosstalk correction,
    since the crosstalk is at a low level. The median should be sufficient.
 
    Parameters
    ----------
    mi : `lsst.afw.image.MaskedImage`
        MaskedImage for which to measure background.
    badPixels : `list` of `str`
        Mask planes to ignore.
 
    Returns
    -------
    bg : `float`
        Median background level.
    """
    mask = mi.getMask()
    stats = lsst.afw.math.StatisticsControl()
    stats.setAndMask(mask.getPlaneBitMask(badPixels))
    return lsst.afw.math.makeStatistics(mi, lsst.afw.math.MEDIAN, stats).getValue()
 
 
def subtractCrosstalk(exposure, crosstalkCoeffs=None,
                      badPixels=["BAD"], minPixelToMask=45000,
                      crosstalkStr="CROSSTALK", isTrimmed=False,
                      backgroundMethod="None"):
    """Subtract the intra-detector crosstalk from an exposure
 
    We set the mask plane indicated by ``crosstalkStr`` in a target amplifier
    for pixels in a source amplifier that exceed `minPixelToMask`. Note that
    the correction is applied to all pixels in the amplifier, but only those
    that have a substantial crosstalk are masked with ``crosstalkStr``.
 
    The uncorrected image is used as a template for correction. This is good
    enough if the crosstalk is small (e.g., coefficients < ~ 1e-3), but if it's
    larger you may want to iterate.
 
    This method needs unittests (DM-18876), but such testing requires
    DM-18610 to allow the test detector to have the crosstalk
    parameters set.
 
    Parameters
    ----------
    exposure : `lsst.afw.image.Exposure`
        Exposure for which to subtract crosstalk.
    crosstalkCoeffs : `numpy.ndarray`
        Coefficients to use to correct crosstalk.
    badPixels : `list` of `str`
        Mask planes to ignore.
    minPixelToMask : `float`
        Minimum pixel value (relative to the background level) in
        source amplifier for which to set ``crosstalkStr`` mask plane
        in target amplifier.
    crosstalkStr : `str`
        Mask plane name for pixels greatly modified by crosstalk.
    isTrimmed : `bool`
        The image is already trimmed.
        This should no longer be needed once DM-15409 is resolved.
    backgroundMethod : `str`
        Method used to subtract the background.  "AMP" uses
        amplifier-by-amplifier background levels, "DETECTOR" uses full
        exposure/maskedImage levels.  Any other value results in no
        background subtraction.
    """
    mi = exposure.getMaskedImage()
    mask = mi.getMask()
 
    ccd = exposure.getDetector()
    numAmps = len(ccd)
    if crosstalkCoeffs is None:
        coeffs = ccd.getCrosstalk()
    else:
        coeffs = crosstalkCoeffs
    assert coeffs.shape == (numAmps, numAmps)
 
    # Set background level based on the requested method.  The
    # thresholdBackground holds the offset needed so that we only mask
    # pixels high relative to the background, not in an absolute
    # sense.
    thresholdBackground = calculateBackground(mi, badPixels)
 
    backgrounds = [0.0 for amp in ccd]
    if backgroundMethod is None:
        pass
    elif backgroundMethod == "AMP":
        backgrounds = [calculateBackground(mi[amp.getBBox()], badPixels) for amp in ccd]
    elif backgroundMethod == "DETECTOR":
        backgrounds = [calculateBackground(mi, badPixels) for amp in ccd]
 
    # Set the crosstalkStr bit for the bright pixels (those which will have significant crosstalk correction)
    crosstalkPlane = mask.addMaskPlane(crosstalkStr)
    footprints = lsst.afw.detection.FootprintSet(mi, lsst.afw.detection.Threshold(minPixelToMask +
                                                                                  thresholdBackground))
    footprints.setMask(mask, crosstalkStr)
    crosstalk = mask.getPlaneBitMask(crosstalkStr)
 
    # Do pixel level crosstalk correction.
    subtrahend = mi.Factory(mi.getBBox())
    subtrahend.set((0, 0, 0))
    for ii, iAmp in enumerate(ccd):
        iImage = subtrahend[iAmp.getBBox() if isTrimmed else iAmp.getRawDataBBox()]
        for jj, jAmp in enumerate(ccd):
            if ii == jj:
                assert coeffs[ii, jj] == 0.0
            if coeffs[ii, jj] == 0.0:
                continue
 
            jImage = extractAmp(mi, jAmp, iAmp.getReadoutCorner(), isTrimmed)
            jImage.getMask().getArray()[:] &= crosstalk  # Remove all other masks
            jImage -= backgrounds[jj]
 
            iImage.scaledPlus(coeffs[ii, jj], jImage)
 
    # Set crosstalkStr bit only for those pixels that have been significantly modified (i.e., those
    # masked as such in 'subtrahend'), not necessarily those that are bright originally.
    mask.clearMaskPlane(crosstalkPlane)
    mi -= subtrahend  # also sets crosstalkStr bit for bright pixels
 
 
def writeCrosstalkCoeffs(outputFileName, coeff, det=None, crosstalkName="Unknown", indent=2):
    """Write a yaml file containing the crosstalk coefficients
 
    The coeff array is indexed by [i, j] where i and j are amplifiers
    corresponding to the amplifiers in det
 
    Parameters
    ----------
    outputFileName : `str`
        Name of output yaml file
    coeff : `numpy.array(namp, namp)`
        numpy array of coefficients
    det : `lsst.afw.cameraGeom.Detector`
        Used to provide the list of amplifier names;
        if None use ['0', '1', ...]
    ccdType : `str`
        Name of detector, used to index the yaml file
        If all detectors are identical could be the type (e.g. ITL)
    indent : `int`
        Indent width to use when writing the yaml file
    """
 
    if det is None:
        ampNames = [str(i) for i in range(coeff.shape[0])]
    else:
        ampNames = [a.getName() for a in det]
 
    assert coeff.shape == (len(ampNames), len(ampNames))
 
    dIndent = indent
    indent = 0
    with open(outputFileName, "w") as fd:
        print(indent*" " + "crosstalk :", file=fd)
        indent += dIndent
        print(indent*" " + "%s :" % crosstalkName, file=fd)
        indent += dIndent
 
        for i, ampNameI in enumerate(ampNames):
            print(indent*" " + "%s : {" % ampNameI, file=fd)
            indent += dIndent
            print(indent*" ", file=fd, end='')
 
            for j, ampNameJ in enumerate(ampNames):
                print("%s : %11.4e, " % (ampNameJ, coeff[i, j]), file=fd,
                      end='\n' + indent*" " if j%4 == 3 else '')
            print("}", file=fd)
 
            indent -= dIndent