dipoleMeasurement.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 numpy as np
 
import lsst.afw.image as afwImage
import lsst.geom as geom
import lsst.pex.config as pexConfig
import lsst.meas.deblender.baseline as deblendBaseline
from lsst.meas.base.pluginRegistry import register
from lsst.meas.base import SingleFrameMeasurementTask, SingleFrameMeasurementConfig, \
    SingleFramePluginConfig, SingleFramePlugin
import lsst.afw.display as afwDisplay
from lsst.utils.logging import getLogger
 
__all__ = ("DipoleMeasurementConfig", "DipoleMeasurementTask", "DipoleAnalysis", "DipoleDeblender",
           "SourceFlagChecker", "ClassificationDipoleConfig", "ClassificationDipolePlugin")
 
 
class ClassificationDipoleConfig(SingleFramePluginConfig):
    """Configuration for classification of detected diaSources as dipole or not"""
    minSn = pexConfig.Field(
        doc="Minimum quadrature sum of positive+negative lobe S/N to be considered a dipole",
        dtype=float, default=np.sqrt(2) * 5.0,
    )
    maxFluxRatio = pexConfig.Field(
        doc="Maximum flux ratio in either lobe to be considered a dipole",
        dtype=float, default=0.65
    )
 
 
@register("ip_diffim_ClassificationDipole")
class ClassificationDipolePlugin(SingleFramePlugin):
    """A plugin to classify whether a diaSource is a dipole.
    """
 
    ConfigClass = ClassificationDipoleConfig
 
    @classmethod
    def getExecutionOrder(cls):
        """
        Returns
        -------
        result : `callable`
        """
        return cls.APCORR_ORDER
 
    def __init__(self, config, name, schema, metadata):
        SingleFramePlugin.__init__(self, config, name, schema, metadata)
        self.dipoleAnalysis = DipoleAnalysis()
        self.keyProbability = schema.addField(name + "_value", type="D",
                                              doc="Set to 1 for dipoles, else 0.")
        self.keyFlag = schema.addField(name + "_flag", type="Flag", doc="Set to 1 for any fatal failure.")
 
    def measure(self, measRecord, exposure):
        passesSn = self.dipoleAnalysis.getSn(measRecord) > self.config.minSn
        negFlux = np.abs(measRecord.get("ip_diffim_PsfDipoleFlux_neg_instFlux"))
        negFluxFlag = measRecord.get("ip_diffim_PsfDipoleFlux_neg_flag")
        posFlux = np.abs(measRecord.get("ip_diffim_PsfDipoleFlux_pos_instFlux"))
        posFluxFlag = measRecord.get("ip_diffim_PsfDipoleFlux_pos_flag")
 
        if negFluxFlag or posFluxFlag:
            self.failfail(measRecord)
            # continue on to classify
 
        totalFlux = negFlux + posFlux
 
        # If negFlux or posFlux are NaN, these evaluate to False
        passesFluxNeg = (negFlux / totalFlux) < self.config.maxFluxRatio
        passesFluxPos = (posFlux / totalFlux) < self.config.maxFluxRatio
        if (passesSn and passesFluxPos and passesFluxNeg):
            val = 1.0
        else:
            val = 0.0
 
        measRecord.set(self.keyProbability, val)
 
    def fail(self, measRecord, error=None):
        measRecord.set(self.keyFlag, True)
 
 
class DipoleMeasurementConfig(SingleFrameMeasurementConfig):
    """Measurement of detected diaSources as dipoles"""
 
    def setDefaults(self):
        SingleFrameMeasurementConfig.setDefaults(self)
        self.pluginsplugins = ["base_CircularApertureFlux",
                        "base_PixelFlags",
                        "base_SkyCoord",
                        "base_PsfFlux",
                        "ip_diffim_NaiveDipoleCentroid",
                        "ip_diffim_NaiveDipoleFlux",
                        "ip_diffim_PsfDipoleFlux",
                        "ip_diffim_ClassificationDipole",
                        ]
 
        self.slots.calibFlux = None
        self.slots.modelFlux = None
        self.slots.gaussianFlux = None
        self.slots.shape = None
        self.slots.centroid = "ip_diffim_NaiveDipoleCentroid"
        self.doReplaceWithNoisedoReplaceWithNoise = False
 
 
class DipoleMeasurementTask(SingleFrameMeasurementTask):
    """Measurement of Sources, specifically ones from difference images, for characterization as dipoles
 
    Parameters
    ----------
    sources : 'lsst.afw.table.SourceCatalog'
        Sources that will be measured
    badFlags : `list` of `dict`
        A list of flags that will be used to determine if there was a measurement problem
 
    Notes
    -----
    The list of badFlags will be used to make a list of keys to check for measurement flags on.  By
    default the centroid keys are added to this list
 
    Description
 
    This class provides a default configuration for running Source measurement on image differences.
 
    .. code-block:: py
 
        class DipoleMeasurementConfig(SingleFrameMeasurementConfig):
            "Measurement of detected diaSources as dipoles"
            def setDefaults(self):
                SingleFrameMeasurementConfig.setDefaults(self)
                self.plugins = ["base_CircularApertureFlux",
                                "base_PixelFlags",
                                "base_SkyCoord",
                                "base_PsfFlux",
                                "ip_diffim_NaiveDipoleCentroid",
                                "ip_diffim_NaiveDipoleFlux",
                                "ip_diffim_PsfDipoleFlux",
                                "ip_diffim_ClassificationDipole",
                                ]
                self.slots.calibFlux = None
                self.slots.modelFlux = None
                self.slots.instFlux = None
                self.slots.shape = None
                self.slots.centroid = "ip_diffim_NaiveDipoleCentroid"
                self.doReplaceWithNoise = False
 
    These plugins enabled by default allow the user to test the hypothesis that the Source is a dipole.
    This includes a set of measurements derived from intermediate base classes
    DipoleCentroidAlgorithm and DipoleFluxAlgorithm.
    Their respective algorithm control classes are defined in
    DipoleCentroidControl and DipoleFluxControl.
    Each centroid and flux measurement will have _neg (negative)
    and _pos (positive lobe) fields.
 
    The first set of measurements uses a "naive" alrogithm
    for centroid and flux measurements, implemented in
    NaiveDipoleCentroidControl and NaiveDipoleFluxControl.
    The algorithm uses a naive 3x3 weighted moment around
    the nominal centroids of each peak in the Source Footprint.  These algorithms fill the table fields
    ip_diffim_NaiveDipoleCentroid* and ip_diffim_NaiveDipoleFlux*
 
    The second set of measurements undertakes a joint-Psf model on the negative
    and positive lobe simultaneously. This fit simultaneously solves for the negative and positive
    lobe centroids and fluxes using non-linear least squares minimization.
    The fields are stored in table elements ip_diffim_PsfDipoleFlux*.
 
    Because this Task is just a config for SingleFrameMeasurementTask, the same result may be acheived by
    manually editing the config and running SingleFrameMeasurementTask. For example:
 
    .. code-block:: py
 
        config = SingleFrameMeasurementConfig()
        config.plugins.names = ["base_PsfFlux",
                                "ip_diffim_PsfDipoleFlux",
                                "ip_diffim_NaiveDipoleFlux",
                                "ip_diffim_NaiveDipoleCentroid",
                                "ip_diffim_ClassificationDipole",
                                "base_CircularApertureFlux",
                                "base_SkyCoord"]
 
        config.slots.calibFlux = None
        config.slots.modelFlux = None
        config.slots.instFlux = None
        config.slots.shape = None
        config.slots.centroid = "ip_diffim_NaiveDipoleCentroid"
        config.doReplaceWithNoise = False
 
        schema = afwTable.SourceTable.makeMinimalSchema()
        task = SingleFrameMeasurementTask(schema, config=config)-
 
    Debug variables
 
    The ``pipetask`` command line interface supports a
    flag --debug to import @b debug.py from your PYTHONPATH.  The relevant contents of debug.py
    for this Task include:
 
    .. code-block:: py
 
        import sys
        import lsstDebug
        def DebugInfo(name):
            di = lsstDebug.getInfo(name)
            if name == "lsst.ip.diffim.dipoleMeasurement":
                di.display = True                 # enable debug output
                di.maskTransparency = 90          # display mask transparency
                di.displayDiaSources = True       # show exposure with dipole results
            return di
        lsstDebug.Info = DebugInfo
        lsstDebug.frame = 1
 
        config.slots.calibFlux = None
        config.slots.modelFlux = None
        config.slots.gaussianFlux = None
        config.slots.shape = None
        config.slots.centroid = "ip_diffim_NaiveDipoleCentroid"
        config.doReplaceWithNoise = False
 
    This code is dipoleMeasTask.py in the examples directory, and can be run as e.g.
 
    .. code-block:: none
 
        examples/dipoleMeasTask.py
        examples/dipoleMeasTask.py --debug
        examples/dipoleMeasTask.py --debug --image /path/to/image.fits
 
 
 
    Start the processing by parsing the command line, where the user has the option of
    enabling debugging output and/or sending their own image for demonstration
    (in case they have not downloaded the afwdata package).
 
    .. code-block:: py
 
        if __name__ == "__main__":
            import argparse
            parser = argparse.ArgumentParser(
                description="Demonstrate the use of SourceDetectionTask and DipoleMeasurementTask")
            parser.add_argument('--debug', '-d', action="store_true", help="Load debug.py?", default=False)
            parser.add_argument("--image", "-i", help="User defined image", default=None)
            args = parser.parse_args()
            if args.debug:
                try:
                    import debug
                    debug.lsstDebug.frame = 2
                except ImportError as e:
                    print(e, file=sys.stderr)
            run(args)
 
    The processing occurs in the run function.  We first extract an exposure from disk or afwdata, displaying
    it if requested:
 
    .. code-block:: py
 
        def run(args):
            exposure = loadData(args.image)
            if args.debug:
                afwDisplay.Display(frame=1).mtv(exposure)
 
    Create a default source schema that we will append fields to as we add more algorithms:
 
    .. code-block:: py
 
        schema = afwTable.SourceTable.makeMinimalSchema()
 
    Create the detection and measurement Tasks, with some minor tweaking of their configs:
 
    .. code-block:: py
 
            # Create the detection task
        config = SourceDetectionTask.ConfigClass()
        config.thresholdPolarity = "both"
        config.background.isNanSafe = True
        config.thresholdValue = 3
        detectionTask = SourceDetectionTask(config=config, schema=schema)
        # And the measurement Task
        config = DipoleMeasurementTask.ConfigClass()
        config.plugins.names.remove('base_SkyCoord')
        algMetadata = dafBase.PropertyList()
        measureTask = DipoleMeasurementTask(schema, algMetadata, config=config)
 
    Having fully initialied the schema, we create a Source table from it:
 
    .. code-block:: py
 
        # Create the output table
        tab = afwTable.SourceTable.make(schema)
 
    Run detection:
 
    .. code-block:: py
 
        # Process the data
        results = detectionTask.run(tab, exposure)
 
    Because we are looking for dipoles, we need to merge the positive and negative detections:
 
    .. code-block:: py
 
        # Merge the positve and negative sources
        fpSet = results.fpSets.positive
        growFootprint = 2
        fpSet.merge(results.fpSets.negative, growFootprint, growFootprint, False)
        diaSources = afwTable.SourceCatalog(tab)
        fpSet.makeSources(diaSources)
        print("Merged %s Sources into %d diaSources (from %d +ve, %d -ve)" % (len(results.sources),
                                                                          len(diaSources),
                                                                          results.fpSets.numPos,
                                                                          results.fpSets.numNeg))
 
    Finally, perform measurement (both standard and dipole-specialized) on the merged sources:
 
    .. code-block:: py
 
        measureTask.run(diaSources, exposure)
 
    Optionally display debugging information:
 
    .. code-block:: py
 
        # Display dipoles if debug enabled
        if args.debug:
            dpa = DipoleAnalysis()
            dpa.displayDipoles(exposure, diaSources)
 
    """
    ConfigClass = DipoleMeasurementConfig
    _DefaultName = "dipoleMeasurement"
 
 
 
 
class SourceFlagChecker(object):
    """Functor class to check whether a diaSource has flags set that should cause it to be labeled bad."""
 
    def __init__(self, sources, badFlags=None):
        self.badFlags = ['base_PixelFlags_flag_edge', 'base_PixelFlags_flag_interpolatedCenter',
                         'base_PixelFlags_flag_saturatedCenter']
        if badFlags is not None:
            for flag in badFlags:
                self.badFlags.append(flag)
        self.keys = [sources.getSchema().find(name).key for name in self.badFlags]
        self.keys.append(sources.table.getCentroidFlagKey())
 
    def __call__(self, source):
        """Call the source flag checker on a single Source
 
        Parameters
        ----------
        source :
            Source that will be examined
        """
        for k in self.keys:
            if source.get(k):
                return False
        return True
 
 
class DipoleAnalysis(object):
    """Functor class that provides (S/N, position, orientation) of measured dipoles"""
 
    def __init__(self):
        pass
 
    def __call__(self, source):
        """Parse information returned from dipole measurement
 
        Parameters
        ----------
        source : `lsst.afw.table.SourceRecord`
            The source that will be examined"""
        return self.getSn(source), self.getCentroid(source), self.getOrientation(source)
 
    def getSn(self, source):
        """Get the total signal-to-noise of the dipole; total S/N is from positive and negative lobe
 
        Parameters
        ----------
        source : `lsst.afw.table.SourceRecord`
            The source that will be examined"""
 
        posflux = source.get("ip_diffim_PsfDipoleFlux_pos_instFlux")
        posfluxErr = source.get("ip_diffim_PsfDipoleFlux_pos_instFluxErr")
        negflux = source.get("ip_diffim_PsfDipoleFlux_neg_instFlux")
        negfluxErr = source.get("ip_diffim_PsfDipoleFlux_neg_instFluxErr")
 
        # Not a dipole!
        if (posflux < 0) is (negflux < 0):
            return 0
 
        return np.sqrt((posflux/posfluxErr)**2 + (negflux/negfluxErr)**2)
 
    def getCentroid(self, source):
        """Get the centroid of the dipole; average of positive and negative lobe
 
        Parameters
        ----------
        source : `lsst.afw.table.SourceRecord`
            The source that will be examined"""
 
        negCenX = source.get("ip_diffim_PsfDipoleFlux_neg_centroid_x")
        negCenY = source.get("ip_diffim_PsfDipoleFlux_neg_centroid_y")
        posCenX = source.get("ip_diffim_PsfDipoleFlux_pos_centroid_x")
        posCenY = source.get("ip_diffim_PsfDipoleFlux_pos_centroid_y")
        if (np.isinf(negCenX) or np.isinf(negCenY) or np.isinf(posCenX) or np.isinf(posCenY)):
            return None
 
        center = geom.Point2D(0.5*(negCenX+posCenX),
                              0.5*(negCenY+posCenY))
        return center
 
    def getOrientation(self, source):
        """Calculate the orientation of dipole; vector from negative to positive lobe
 
        Parameters
        ----------
        source : `lsst.afw.table.SourceRecord`
            The source that will be examined"""
 
        negCenX = source.get("ip_diffim_PsfDipoleFlux_neg_centroid_x")
        negCenY = source.get("ip_diffim_PsfDipoleFlux_neg_centroid_y")
        posCenX = source.get("ip_diffim_PsfDipoleFlux_pos_centroid_x")
        posCenY = source.get("ip_diffim_PsfDipoleFlux_pos_centroid_y")
        if (np.isinf(negCenX) or np.isinf(negCenY) or np.isinf(posCenX) or np.isinf(posCenY)):
            return None
 
        dx, dy = posCenX-negCenX, posCenY-negCenY
        angle = geom.Angle(np.arctan2(dx, dy), geom.radians)
        return angle
 
    def displayDipoles(self, exposure, sources):
        """Display debugging information on the detected dipoles
 
        Parameters
        ----------
        exposure : `lsst.afw.image.Exposure`
            Image the dipoles were measured on
        sources : `lsst.afw.table.SourceCatalog`
            The set of diaSources that were measured"""
 
        import lsstDebug
        display = lsstDebug.Info(__name__).display
        displayDiaSources = lsstDebug.Info(__name__).displayDiaSources
        maskTransparency = lsstDebug.Info(__name__).maskTransparency
        if not maskTransparency:
            maskTransparency = 90
        disp = afwDisplay.Display(frame=lsstDebug.frame)
        disp.setMaskTransparency(maskTransparency)
        disp.mtv(exposure)
 
        if display and displayDiaSources:
            with disp.Buffering():
                for source in sources:
                    cenX, cenY = source.get("ipdiffim_DipolePsfFlux_centroid")
                    if np.isinf(cenX) or np.isinf(cenY):
                        cenX, cenY = source.getCentroid()
 
                    isdipole = source.get("ip_diffim_ClassificationDipole_value")
                    if isdipole and np.isfinite(isdipole):
                        # Dipole
                        ctype = afwDisplay.GREEN
                    else:
                        # Not dipole
                        ctype = afwDisplay.RED
 
                    disp.dot("o", cenX, cenY, size=2, ctype=ctype)
 
                    negCenX = source.get("ip_diffim_PsfDipoleFlux_neg_centroid_x")
                    negCenY = source.get("ip_diffim_PsfDipoleFlux_neg_centroid_y")
                    posCenX = source.get("ip_diffim_PsfDipoleFlux_pos_centroid_x")
                    posCenY = source.get("ip_diffim_PsfDipoleFlux_pos_centroid_y")
                    if (np.isinf(negCenX) or np.isinf(negCenY) or np.isinf(posCenX) or np.isinf(posCenY)):
                        continue
 
                    disp.line([(negCenX, negCenY), (posCenX, posCenY)], ctype=afwDisplay.YELLOW)
 
            lsstDebug.frame += 1
 
 
class DipoleDeblender(object):
    """Functor to deblend a source as a dipole, and return a new source with deblended footprints.
 
       This necessarily overrides some of the functionality from
       meas_algorithms/python/lsst/meas/algorithms/deblend.py since we
       need a single source that contains the blended peaks, not
       multiple children sources.  This directly calls the core
       deblending code deblendBaseline.deblend (optionally _fitPsf for
       debugging).
 
       Not actively being used, but there is a unit test for it in
       dipoleAlgorithm.py.
    """
 
    def __init__(self):
        # Set up defaults to send to deblender
 
        # Always deblend as Psf
        self.psfChisqCut1 = self.psfChisqCut2 = self.psfChisqCut2b = np.inf
        self.log = getLogger('lsst.ip.diffim.DipoleDeblender')
        self.sigma2fwhm = 2. * np.sqrt(2. * np.log(2.))
 
    def __call__(self, source, exposure):
        fp = source.getFootprint()
        peaks = fp.getPeaks()
        peaksF = [pk.getF() for pk in peaks]
        fbb = fp.getBBox()
        fmask = afwImage.Mask(fbb)
        fmask.setXY0(fbb.getMinX(), fbb.getMinY())
        fp.spans.setMask(fmask, 1)
 
        psf = exposure.getPsf()
        psfSigPix = psf.computeShape(psf.getAveragePosition()).getDeterminantRadius()
        psfFwhmPix = psfSigPix * self.sigma2fwhm
        subimage = afwImage.ExposureF(exposure, bbox=fbb, deep=True)
        cpsf = deblendBaseline.CachingPsf(psf)
 
        # if fewer than 2 peaks, just return a copy of the source
        if len(peaks) < 2:
            return source.getTable().copyRecord(source)
 
        # make sure you only deblend 2 peaks; take the brighest and faintest
        speaks = [(p.getPeakValue(), p) for p in peaks]
        speaks.sort()
        dpeaks = [speaks[0][1], speaks[-1][1]]
 
        # and only set these peaks in the footprint (peaks is mutable)
        peaks.clear()
        for peak in dpeaks:
            peaks.append(peak)
 
        if True:
            # Call top-level deblend task
            fpres = deblendBaseline.deblend(fp, exposure.getMaskedImage(), psf, psfFwhmPix,
                                            log=self.log,
                                            psfChisqCut1=self.psfChisqCut1,
                                            psfChisqCut2=self.psfChisqCut2,
                                            psfChisqCut2b=self.psfChisqCut2b)
        else:
            # Call lower-level _fit_psf task
 
            # Prepare results structure
            fpres = deblendBaseline.DeblenderResult(fp, exposure.getMaskedImage(), psf, psfFwhmPix, self.log)
 
            for pki, (pk, pkres, pkF) in enumerate(zip(dpeaks, fpres.deblendedParents[0].peaks, peaksF)):
                self.log.debug('Peak %i', pki)
                deblendBaseline._fitPsf(fp, fmask, pk, pkF, pkres, fbb, dpeaks, peaksF, self.log,
                                        cpsf, psfFwhmPix,
                                        subimage.getMaskedImage().getImage(),
                                        subimage.getMaskedImage().getVariance(),
                                        self.psfChisqCut1, self.psfChisqCut2, self.psfChisqCut2b)
 
        deblendedSource = source.getTable().copyRecord(source)
        deblendedSource.setParent(source.getId())
        peakList = deblendedSource.getFootprint().getPeaks()
        peakList.clear()
 
        for i, peak in enumerate(fpres.deblendedParents[0].peaks):
            if peak.psfFitFlux > 0:
                suffix = "pos"
            else:
                suffix = "neg"
            c = peak.psfFitCenter
            self.log.info("deblended.centroid.dipole.psf.%s %f %f",
                          suffix, c[0], c[1])
            self.log.info("deblended.chi2dof.dipole.%s %f",
                          suffix, peak.psfFitChisq / peak.psfFitDof)
            self.log.info("deblended.flux.dipole.psf.%s %f",
                          suffix, peak.psfFitFlux * np.sum(peak.templateImage.getArray()))
            peakList.append(peak.peak)
        return deblendedSource