lsst.meas.algorithms.pcaPsfDeterminer.PcaPsfDeterminerTask Class Reference

Inheritance diagram for lsst.meas.algorithms.pcaPsfDeterminer.PcaPsfDeterminerTask:

Public Member Functions
def	determinePsf (self, exposure, psfCandidateList, metadata=None, flagKey=None)
def	determinePsf (self, exposure, psfCandidateList, metadata=None)

Static Public Attributes
	ConfigClass = PcaPsfDeterminerConfig
bool	usesMatches = False

Detailed Description

A measurePsfTask psf estimator.

Definition at line 162 of file pcaPsfDeterminer.py.

Member Function Documentation

determinePsf() [1/2]

def lsst.meas.algorithms.psfDeterminer.BasePsfDeterminerTask.determinePsf (   self,   exposure,   psfCandidateList,   metadata = None  )

inherited

Determine a PSF model.

Parameters
----------
exposure : `lsst.afw.Exposure`
    Exposure containing the psf candidates.
psdCandidateList : `list` [`lsst.meas.algorithms.PsfCandidate`]
    A sequence of PSF candidates; typically obtained by
    detecting sources and then running them through a star
    selector.
metadata : `str`
    A place to save interesting items.

Returns
-------
psf : `lsst.afw.detection.Psf`
    The fit PSF.
cellSet : `lsst.afw.math.SpatialCellSet`
    The spatial cell set used to determine the PSF

Definition at line 78 of file psfDeterminer.py.

    def determinePsf(self, exposure, psfCandidateList, metadata=None):
        """Determine a PSF model.
 
        Parameters
        ----------
        exposure : `lsst.afw.Exposure`
            Exposure containing the psf candidates.
        psdCandidateList : `list` [`lsst.meas.algorithms.PsfCandidate`]
            A sequence of PSF candidates; typically obtained by
            detecting sources and then running them through a star
            selector.
        metadata : `str`
            A place to save interesting items.
 
        Returns
        -------
        psf : `lsst.afw.detection.Psf`
            The fit PSF.
        cellSet : `lsst.afw.math.SpatialCellSet`
            The spatial cell set used to determine the PSF
        """
        raise NotImplementedError("BasePsfDeterminerTask is abstract, subclasses must override this method")
 
 

determinePsf() [2/2]

def lsst.meas.algorithms.pcaPsfDeterminer.PcaPsfDeterminerTask.determinePsf	(		self,
			exposure,
			psfCandidateList,
			metadata = None,
			flagKey = None
	)

Determine a PCA PSF model for an exposure given a list of PSF candidates.

Parameters
----------
exposure : `lsst.afw.image.Exposure`
   Exposure containing the psf candidates.
psfCandidateList : `list` of `lsst.meas.algorithms.PsfCandidate`
   A sequence of PSF candidates typically obtained by detecting sources
   and then running them through a star selector.
metadata : `lsst.daf.base import PropertyList` or `None`, optional
   A home for interesting tidbits of information.
flagKey : `str`, optional
   Schema key used to mark sources actually used in PSF determination.

Returns
-------
psf : `lsst.meas.algorithms.PcaPsf`
   The measured PSF.
psfCellSet : `lsst.afw.math.SpatialCellSet`
   The PSF candidates.

Definition at line 205 of file pcaPsfDeterminer.py.

    def determinePsf(self, exposure, psfCandidateList, metadata=None, flagKey=None):
        """Determine a PCA PSF model for an exposure given a list of PSF candidates.
 
        Parameters
        ----------
        exposure : `lsst.afw.image.Exposure`
           Exposure containing the psf candidates.
        psfCandidateList : `list` of `lsst.meas.algorithms.PsfCandidate`
           A sequence of PSF candidates typically obtained by detecting sources
           and then running them through a star selector.
        metadata : `lsst.daf.base import PropertyList` or `None`, optional
           A home for interesting tidbits of information.
        flagKey : `str`, optional
           Schema key used to mark sources actually used in PSF determination.
 
        Returns
        -------
        psf : `lsst.meas.algorithms.PcaPsf`
           The measured PSF.
        psfCellSet : `lsst.afw.math.SpatialCellSet`
           The PSF candidates.
        """
        import lsstDebug
        display = lsstDebug.Info(__name__).display
        displayExposure = lsstDebug.Info(__name__).displayExposure     # display the Exposure + spatialCells
        displayPsfCandidates = lsstDebug.Info(__name__).displayPsfCandidates  # show the viable candidates
        displayIterations = lsstDebug.Info(__name__).displayIterations  # display on each PSF iteration
        displayPsfComponents = lsstDebug.Info(__name__).displayPsfComponents  # show the PCA components
        displayResiduals = lsstDebug.Info(__name__).displayResiduals         # show residuals
        displayPsfMosaic = lsstDebug.Info(__name__).displayPsfMosaic   # show mosaic of reconstructed PSF(x,y)
        # match Kernel amplitudes for spatial plots
        matchKernelAmplitudes = lsstDebug.Info(__name__).matchKernelAmplitudes
        # Keep matplotlib alive post mortem
        keepMatplotlibPlots = lsstDebug.Info(__name__).keepMatplotlibPlots
        displayPsfSpatialModel = lsstDebug.Info(__name__).displayPsfSpatialModel  # Plot spatial model?
        showBadCandidates = lsstDebug.Info(__name__).showBadCandidates  # Include bad candidates
        # Normalize residuals by object amplitude
        normalizeResiduals = lsstDebug.Info(__name__).normalizeResiduals
        pause = lsstDebug.Info(__name__).pause                         # Prompt user after each iteration?
 
        if display:
            afwDisplay.setDefaultMaskTransparency(75)
        if display > 1:
            pause = True
 
        mi = exposure.getMaskedImage()
 
        if len(psfCandidateList) == 0:
            raise RuntimeError("No PSF candidates supplied.")
 
        # construct and populate a spatial cell set
        bbox = mi.getBBox()
        psfCellSet = afwMath.SpatialCellSet(bbox, self.config.sizeCellX, self.config.sizeCellY)
        sizes = []
        for i, psfCandidate in enumerate(psfCandidateList):
            if psfCandidate.getSource().getPsfFluxFlag():  # bad measurement
                continue
 
            try:
                psfCellSet.insertCandidate(psfCandidate)
            except Exception as e:
                self.log.debug("Skipping PSF candidate %d of %d: %s", i, len(psfCandidateList), e)
                continue
            source = psfCandidate.getSource()
 
            quad = afwGeom.Quadrupole(source.getIxx(), source.getIyy(), source.getIxy())
            axes = afwEll.Axes(quad)
            sizes.append(axes.getA())
        if len(sizes) == 0:
            raise RuntimeError("No usable PSF candidates supplied")
        nEigenComponents = self.config.nEigenComponents  # initial version
 
        if self.config.kernelSize >= 15:
            self.log.warning("WARNING: NOT scaling kernelSize by stellar quadrupole moment "
                             "because config.kernelSize=%s >= 15; "
                             "using config.kernelSize as the width, instead",
                             self.config.kernelSize)
            actualKernelSize = int(self.config.kernelSize)
        else:
            medSize = numpy.median(sizes)
            actualKernelSize = 2*int(self.config.kernelSize*math.sqrt(medSize) + 0.5) + 1
            if actualKernelSize < self.config.kernelSizeMin:
                actualKernelSize = self.config.kernelSizeMin
            if actualKernelSize > self.config.kernelSizeMax:
                actualKernelSize = self.config.kernelSizeMax
 
            if display:
                print("Median size=%s" % (medSize,))
        self.log.trace("Kernel size=%s", actualKernelSize)
 
        # Set size of image returned around candidate
        psfCandidateList[0].setHeight(actualKernelSize)
        psfCandidateList[0].setWidth(actualKernelSize)
 
        if self.config.doRejectBlends:
            # Remove blended candidates completely
            blendedCandidates = []  # Candidates to remove; can't do it while iterating
            for cell, cand in candidatesIter(psfCellSet, False):
                if len(cand.getSource().getFootprint().getPeaks()) > 1:
                    blendedCandidates.append((cell, cand))
                    continue
            if display:
                print("Removing %d blended Psf candidates" % len(blendedCandidates))
            for cell, cand in blendedCandidates:
                cell.removeCandidate(cand)
            if sum(1 for cand in candidatesIter(psfCellSet, False)) == 0:
                raise RuntimeError("All PSF candidates removed as blends")
 
        if display:
            if displayExposure:
                disp = afwDisplay.Display(frame=0)
                disp.mtv(exposure, title="psf determination")
                utils.showPsfSpatialCells(exposure, psfCellSet, self.config.nStarPerCell, symb="o",
                                          ctype=afwDisplay.CYAN, ctypeUnused=afwDisplay.YELLOW,
                                          size=4, display=disp)
 
        #
        # Do a PCA decomposition of those PSF candidates
        #
        reply = "y"                         # used in interactive mode
        for iterNum in range(self.config.nIterForPsf):
            if display and displayPsfCandidates:  # Show a mosaic of usable PSF candidates
 
                stamps = []
                for cell in psfCellSet.getCellList():
                    for cand in cell.begin(not showBadCandidates):  # maybe include bad candidates
                        try:
                            im = cand.getMaskedImage()
 
                            chi2 = cand.getChi2()
                            if chi2 > 1e100:
                                chi2 = numpy.nan
 
                            stamps.append((im, "%d%s" %
                                           (utils.splitId(cand.getSource().getId(), True)["objId"], chi2),
                                           cand.getStatus()))
                        except Exception:
                            continue
 
                if len(stamps) == 0:
                    print("WARNING: No PSF candidates to show; try setting showBadCandidates=True")
                else:
                    mos = afwDisplay.utils.Mosaic()
                    for im, label, status in stamps:
                        im = type(im)(im, True)
                        try:
                            im /= afwMath.makeStatistics(im, afwMath.MAX).getValue()
                        except NotImplementedError:
                            pass
 
                        mos.append(im, label,
                                   (afwDisplay.GREEN if status == afwMath.SpatialCellCandidate.GOOD else
                                    afwDisplay.YELLOW if status == afwMath.SpatialCellCandidate.UNKNOWN else
                                    afwDisplay.RED))
 
                    disp8 = afwDisplay.Display(frame=8)
                    mos.makeMosaic(display=disp8, title="Psf Candidates")
 
            # Re-fit until we don't have any candidates with naughty chi^2 values influencing the fit
            cleanChi2 = False  # Any naughty (negative/NAN) chi^2 values?
            while not cleanChi2:
                cleanChi2 = True
                #
                # First, estimate the PSF
                #
                psf, eigenValues, nEigenComponents, fitChi2 = \
                    self._fitPsf(exposure, psfCellSet, actualKernelSize, nEigenComponents)
                #
                # In clipping, allow all candidates to be innocent until proven guilty on this iteration.
                # Throw out any prima facie guilty candidates (naughty chi^2 values)
                #
                for cell in psfCellSet.getCellList():
                    awfulCandidates = []
                    for cand in cell.begin(False):  # include bad candidates
                        cand.setStatus(afwMath.SpatialCellCandidate.UNKNOWN)  # until proven guilty
                        rchi2 = cand.getChi2()
                        if not numpy.isfinite(rchi2) or rchi2 <= 0:
                            # Guilty prima facie
                            awfulCandidates.append(cand)
                            cleanChi2 = False
                            self.log.debug("chi^2=%s; id=%s",
                                           cand.getChi2(), cand.getSource().getId())
                    for cand in awfulCandidates:
                        if display:
                            print("Removing bad candidate: id=%d, chi^2=%f" %
                                  (cand.getSource().getId(), cand.getChi2()))
                        cell.removeCandidate(cand)
 
            #
            # Clip out bad fits based on reduced chi^2
            #
            badCandidates = list()
            for cell in psfCellSet.getCellList():
                for cand in cell.begin(False):  # include bad candidates
                    rchi2 = cand.getChi2()  # reduced chi^2 when fitting PSF to candidate
                    assert rchi2 > 0
                    if rchi2 > self.config.reducedChi2ForPsfCandidates:
                        badCandidates.append(cand)
 
            badCandidates.sort(key=lambda x: x.getChi2(), reverse=True)
            numBad = numCandidatesToReject(len(badCandidates), iterNum,
                                           self.config.nIterForPsf)
            for i, c in zip(range(numBad), badCandidates):
                if display:
                    chi2 = c.getChi2()
                    if chi2 > 1e100:
                        chi2 = numpy.nan
 
                    print("Chi^2 clipping %-4d  %.2g" % (c.getSource().getId(), chi2))
                c.setStatus(afwMath.SpatialCellCandidate.BAD)
 
            #
            # Clip out bad fits based on spatial fitting.
            #
            # This appears to be better at getting rid of sources that have a single dominant kernel component
            # (other than the zeroth; e.g., a nearby contaminant) because the surrounding sources (which help
            # set the spatial model) don't contain that kernel component, and so the spatial modeling
            # downweights the component.
            #
 
            residuals = list()
            candidates = list()
            kernel = psf.getKernel()
            noSpatialKernel = psf.getKernel()
            for cell in psfCellSet.getCellList():
                for cand in cell.begin(False):
                    candCenter = lsst.geom.PointD(cand.getXCenter(), cand.getYCenter())
                    try:
                        im = cand.getMaskedImage(kernel.getWidth(), kernel.getHeight())
                    except Exception:
                        continue
 
                    fit = fitKernelParamsToImage(noSpatialKernel, im, candCenter)
                    params = fit[0]
                    kernels = fit[1]
                    amp = 0.0
                    for p, k in zip(params, kernels):
                        amp += p*k.getSum()
 
                    predict = [kernel.getSpatialFunction(k)(candCenter.getX(), candCenter.getY()) for
                               k in range(kernel.getNKernelParameters())]
 
                    residuals.append([a/amp - p for a, p in zip(params, predict)])
                    candidates.append(cand)
 
            residuals = numpy.array(residuals)
 
            for k in range(kernel.getNKernelParameters()):
                if False:
                    # Straight standard deviation
                    mean = residuals[:, k].mean()
                    rms = residuals[:, k].std()
                elif False:
                    # Using interquartile range
                    sr = numpy.sort(residuals[:, k])
                    mean = (sr[int(0.5*len(sr))] if len(sr)%2 else
                            0.5*(sr[int(0.5*len(sr))] + sr[int(0.5*len(sr)) + 1]))
                    rms = 0.74*(sr[int(0.75*len(sr))] - sr[int(0.25*len(sr))])
                else:
                    stats = afwMath.makeStatistics(residuals[:, k], afwMath.MEANCLIP | afwMath.STDEVCLIP)
                    mean = stats.getValue(afwMath.MEANCLIP)
                    rms = stats.getValue(afwMath.STDEVCLIP)
 
                rms = max(1.0e-4, rms)  # Don't trust RMS below this due to numerical issues
 
                if display:
                    print("Mean for component %d is %f" % (k, mean))
                    print("RMS for component %d is %f" % (k, rms))
                badCandidates = list()
                for i, cand in enumerate(candidates):
                    if numpy.fabs(residuals[i, k] - mean) > self.config.spatialReject*rms:
                        badCandidates.append(i)
 
                badCandidates.sort(key=lambda x: numpy.fabs(residuals[x, k] - mean), reverse=True)
 
                numBad = numCandidatesToReject(len(badCandidates), iterNum,
                                               self.config.nIterForPsf)
 
                for i, c in zip(range(min(len(badCandidates), numBad)), badCandidates):
                    cand = candidates[c]
                    if display:
                        print("Spatial clipping %d (%f,%f) based on %d: %f vs %f" %
                              (cand.getSource().getId(), cand.getXCenter(), cand.getYCenter(), k,
                               residuals[badCandidates[i], k], self.config.spatialReject*rms))
                    cand.setStatus(afwMath.SpatialCellCandidate.BAD)
 
            #
            # Display results
            #
            if display and displayIterations:
                if displayExposure:
                    if iterNum > 0:
                        disp.erase()
                    utils.showPsfSpatialCells(exposure, psfCellSet, self.config.nStarPerCell, showChi2=True,
                                              symb="o", size=8, display=disp, ctype=afwDisplay.YELLOW,
                                              ctypeBad=afwDisplay.RED, ctypeUnused=afwDisplay.MAGENTA)
                    if self.config.nStarPerCellSpatialFit != self.config.nStarPerCell:
                        utils.showPsfSpatialCells(exposure, psfCellSet, self.config.nStarPerCellSpatialFit,
                                                  symb="o", size=10, display=disp,
                                                  ctype=afwDisplay.YELLOW, ctypeBad=afwDisplay.RED)
                if displayResiduals:
                    while True:
                        try:
                            disp4 = afwDisplay.Display(frame=4)
                            utils.showPsfCandidates(exposure, psfCellSet, psf=psf, display=disp4,
                                                    normalize=normalizeResiduals,
                                                    showBadCandidates=showBadCandidates)
                            disp5 = afwDisplay.Display(frame=5)
                            utils.showPsfCandidates(exposure, psfCellSet, psf=psf, display=disp5,
                                                    normalize=normalizeResiduals,
                                                    showBadCandidates=showBadCandidates,
                                                    variance=True)
                        except Exception:
                            if not showBadCandidates:
                                showBadCandidates = True
                                continue
                        break
 
                if displayPsfComponents:
                    disp6 = afwDisplay.Display(frame=6)
                    utils.showPsf(psf, eigenValues, display=disp6)
                if displayPsfMosaic:
                    disp7 = afwDisplay.Display(frame=7)
                    utils.showPsfMosaic(exposure, psf, display=disp7, showFwhm=True)
                    disp7.scale('linear', 0, 1)
                if displayPsfSpatialModel:
                    utils.plotPsfSpatialModel(exposure, psf, psfCellSet, showBadCandidates=True,
                                              matchKernelAmplitudes=matchKernelAmplitudes,
                                              keepPlots=keepMatplotlibPlots)
 
                if pause:
                    while True:
                        try:
                            reply = input("Next iteration? [ynchpqQs] ").strip()
                        except EOFError:
                            reply = "n"
 
                        reply = reply.split()
                        if reply:
                            reply, args = reply[0], reply[1:]
                        else:
                            reply = ""
 
                        if reply in ("", "c", "h", "n", "p", "q", "Q", "s", "y"):
                            if reply == "c":
                                pause = False
                            elif reply == "h":
                                print("c[ontinue without prompting] h[elp] n[o] p[db] q[uit displaying] "
                                      "s[ave fileName] y[es]")
                                continue
                            elif reply == "p":
                                import pdb
                                pdb.set_trace()
                            elif reply == "q":
                                display = False
                            elif reply == "Q":
                                sys.exit(1)
                            elif reply == "s":
                                fileName = args.pop(0)
                                if not fileName:
                                    print("Please provide a filename")
                                    continue
 
                                print("Saving to %s" % fileName)
                                utils.saveSpatialCellSet(psfCellSet, fileName=fileName)
                                continue
                            break
                        else:
                            print("Unrecognised response: %s" % reply, file=sys.stderr)
 
                    if reply == "n":
                        break
 
        # One last time, to take advantage of the last iteration
        psf, eigenValues, nEigenComponents, fitChi2 = \
            self._fitPsf(exposure, psfCellSet, actualKernelSize, nEigenComponents)
 
        #
        # Display code for debugging
        #
        if display and reply != "n":
            disp = afwDisplay.Display(frame=0)
            if displayExposure:
                utils.showPsfSpatialCells(exposure, psfCellSet, self.config.nStarPerCell, showChi2=True,
                                          symb="o", ctype=afwDisplay.YELLOW, ctypeBad=afwDisplay.RED,
                                          size=8, display=disp)
                if self.config.nStarPerCellSpatialFit != self.config.nStarPerCell:
                    utils.showPsfSpatialCells(exposure, psfCellSet, self.config.nStarPerCellSpatialFit,
                                              symb="o", ctype=afwDisplay.YELLOW, ctypeBad=afwDisplay.RED,
                                              size=10, display=disp)
                if displayResiduals:
                    disp4 = afwDisplay.Display(frame=4)
                    utils.showPsfCandidates(exposure, psfCellSet, psf=psf, display=disp4,
                                            normalize=normalizeResiduals,
                                            showBadCandidates=showBadCandidates)
 
            if displayPsfComponents:
                disp6 = afwDisplay.Display(frame=6)
                utils.showPsf(psf, eigenValues, display=disp6)
 
            if displayPsfMosaic:
                disp7 = afwDisplay.Display(frame=7)
                utils.showPsfMosaic(exposure, psf, display=disp7, showFwhm=True)
                disp7.scale("linear", 0, 1)
            if displayPsfSpatialModel:
                utils.plotPsfSpatialModel(exposure, psf, psfCellSet, showBadCandidates=True,
                                          matchKernelAmplitudes=matchKernelAmplitudes,
                                          keepPlots=keepMatplotlibPlots)
        #
        # Generate some QA information
        #
        # Count PSF stars
        #
        numGoodStars = 0
        numAvailStars = 0
 
        avgX = 0.0
        avgY = 0.0
 
        for cell in psfCellSet.getCellList():
            for cand in cell.begin(False):  # don't ignore BAD stars
                numAvailStars += 1
 
            for cand in cell.begin(True):  # do ignore BAD stars
                src = cand.getSource()
                if flagKey is not None:
                    src.set(flagKey, True)
                avgX += src.getX()
                avgY += src.getY()
                numGoodStars += 1
 
        avgX /= numGoodStars
        avgY /= numGoodStars
 
        if metadata is not None:
            metadata["spatialFitChi2"] = fitChi2
            metadata["numGoodStars"] = numGoodStars
            metadata["numAvailStars"] = numAvailStars
            metadata["avgX"] = avgX
            metadata["avgY"] = avgY
 
        psf = PcaPsf(psf.getKernel(), lsst.geom.Point2D(avgX, avgY))
 
        return psf, psfCellSet
 
 

Member Data Documentation

ConfigClass

lsst.meas.algorithms.pcaPsfDeterminer.PcaPsfDeterminerTask.ConfigClass = PcaPsfDeterminerConfig

static

Definition at line 165 of file pcaPsfDeterminer.py.

usesMatches

bool lsst.meas.algorithms.psfDeterminer.BasePsfDeterminerTask.usesMatches = False

staticinherited

Definition at line 70 of file psfDeterminer.py.

---

The documentation for this class was generated from the following file:

• /j/snowflake/release/lsstsw/stack/lsst-scipipe-0.7.0/Linux64/meas_algorithms/22.0.1-32-gf62ce7b1+aa4237961e/python/lsst/meas/algorithms/pcaPsfDeterminer.py