22__all__ = [
"PcaPsfDeterminerConfig",
"PcaPsfDeterminerTask"]
35from .psfDeterminer
import BasePsfDeterminerTask, psfDeterminerRegistry
36from ._algorithmsLib
import PsfCandidateF
37from ._algorithmsLib
import createKernelFromPsfCandidates, countPsfCandidates, \
38 fitSpatialKernelFromPsfCandidates, fitKernelParamsToImage
39from ._algorithmsLib
import PcaPsf
44 """Return the number of PSF candidates to be rejected.
46 The number of candidates being rejected on each iteration gradually
47 increases, so that on the Nth of M iterations we reject N/M of the bad
52 numBadCandidates : `int`
53 Number of bad candidates under consideration.
56 The number of the current PSF iteration.
59 The total number of PSF iterations.
64 Number of candidates to reject.
66 return int(numBadCandidates*(numIter + 1)//totalIter + 0.5)
70 nonLinearSpatialFit = pexConfig.Field[bool](
71 doc=
"Use non-linear fitter for spatial variation of Kernel",
74 nEigenComponents = pexConfig.Field[int](
75 doc=
"number of eigen components for PSF kernel creation",
77 check=
lambda x: x >= 1
79 spatialOrder = pexConfig.Field[int](
80 doc=
"specify spatial order for PSF kernel creation",
83 sizeCellX = pexConfig.Field[int](
84 doc=
"size of cell used to determine PSF (pixels, column direction)",
87 check=
lambda x: x >= 10,
89 sizeCellY = pexConfig.Field[int](
90 doc=
"size of cell used to determine PSF (pixels, row direction)",
91 default=sizeCellX.default,
93 check=
lambda x: x >= 10,
95 nStarPerCell = pexConfig.Field[int](
96 doc=
"number of stars per psf cell for PSF kernel creation",
99 borderWidth = pexConfig.Field[int](
100 doc=
"Number of pixels to ignore around the edge of PSF candidate postage stamps",
103 nStarPerCellSpatialFit = pexConfig.Field[int](
104 doc=
"number of stars per psf Cell for spatial fitting",
107 constantWeight = pexConfig.Field[bool](
108 doc=
"Should each PSF candidate be given the same weight, independent of magnitude?",
111 nIterForPsf = pexConfig.Field[int](
112 doc=
"number of iterations of PSF candidate star list",
115 tolerance = pexConfig.Field[float](
116 doc=
"tolerance of spatial fitting",
119 lam = pexConfig.Field[float](
120 doc=
"floor for variance is lam*data",
123 reducedChi2ForPsfCandidates = pexConfig.Field[float](
124 doc=
"for psf candidate evaluation",
127 spatialReject = pexConfig.Field[float](
128 doc=
"Rejection threshold (stdev) for candidates based on spatial fit",
131 pixelThreshold = pexConfig.Field[float](
132 doc=
"Threshold (stdev) for rejecting extraneous pixels around candidate; applied if positive",
135 doRejectBlends = pexConfig.Field[bool](
136 doc=
"Reject candidates that are blended?",
139 doMaskBlends = pexConfig.Field[bool](
140 doc=
"Mask blends in image?",
150 """A measurePsfTask psf estimator.
152 ConfigClass = PcaPsfDeterminerConfig
154 def _fitPsf(self, exposure, psfCellSet, kernelSize, nEigenComponents):
155 PsfCandidateF.setPixelThreshold(self.config.pixelThreshold)
156 PsfCandidateF.setMaskBlends(self.config.doMaskBlends)
160 for nEigen
in range(nEigenComponents, 0, -1):
163 kernel, eigenValues = createKernelFromPsfCandidates(
164 psfCellSet, exposure.getDimensions(), exposure.getXY0(), nEigen,
165 self.config.spatialOrder, kernelSize, self.config.nStarPerCell,
166 bool(self.config.constantWeight))
171 raise IndexError(
"No viable PSF candidates survive")
173 self.log.warning(
"%s: reducing number of eigen components", e.what())
178 size = kernelSize + 2*self.config.borderWidth
180 eigenValues = [val/float(countPsfCandidates(psfCellSet, self.config.nStarPerCell)*nu)
181 for val
in eigenValues]
184 status, chi2 = fitSpatialKernelFromPsfCandidates(
185 kernel, psfCellSet, bool(self.config.nonLinearSpatialFit),
186 self.config.nStarPerCellSpatialFit, self.config.tolerance, self.config.lam)
190 return psf, eigenValues, nEigen, chi2
192 def determinePsf(self, exposure, psfCandidateList, metadata=None, flagKey=None):
193 """Determine a PCA PSF model for an exposure given a list of PSF candidates.
198 Exposure containing the psf candidates.
200 A sequence of PSF candidates typically obtained by detecting sources
201 and then running them through a star selector.
202 metadata : `
lsst.daf.base import PropertyList`
or `
None`, optional
203 A home
for interesting tidbits of information.
204 flagKey : `str`, optional
205 Schema key used to mark sources actually used
in PSF determination.
217 displayPsfCandidates =
lsstDebug.Info(__name__).displayPsfCandidates
219 displayPsfComponents =
lsstDebug.Info(__name__).displayPsfComponents
223 matchKernelAmplitudes =
lsstDebug.Info(__name__).matchKernelAmplitudes
225 keepMatplotlibPlots =
lsstDebug.Info(__name__).keepMatplotlibPlots
226 displayPsfSpatialModel =
lsstDebug.Info(__name__).displayPsfSpatialModel
233 afwDisplay.setDefaultMaskTransparency(75)
237 mi = exposure.getMaskedImage()
239 if len(psfCandidateList) == 0:
240 raise RuntimeError(
"No PSF candidates supplied.")
246 for i, psfCandidate
in enumerate(psfCandidateList):
247 if psfCandidate.getSource().getPsfFluxFlag():
251 psfCellSet.insertCandidate(psfCandidate)
252 except Exception
as e:
253 self.log.debug(
"Skipping PSF candidate %d of %d: %s", i, len(psfCandidateList), e)
255 source = psfCandidate.getSource()
258 axes = afwEll.Axes(quad)
259 sizes.append(axes.getA())
261 raise RuntimeError(
"No usable PSF candidates supplied")
262 nEigenComponents = self.config.nEigenComponents
264 actualKernelSize = int(self.config.stampSize)
267 print(
"Median size=%s" % (numpy.median(sizes),))
269 self.log.trace(
"Kernel size=%s", actualKernelSize)
271 if actualKernelSize > psfCandidateList[0].getWidth():
272 self.log.warning(
"Using a region (%d x %d) larger than kernelSize (%d) set while making PSF "
273 "candidates. Consider setting a larger value for kernelSize for "
274 "`makePsfCandidates` to avoid this warning.",
275 actualKernelSize, actualKernelSize, psfCandidateList[0].getWidth())
277 if self.config.doRejectBlends:
279 blendedCandidates = []
281 if len(cand.getSource().getFootprint().getPeaks()) > 1:
282 blendedCandidates.append((cell, cand))
285 print(
"Removing %d blended Psf candidates" % len(blendedCandidates))
286 for cell, cand
in blendedCandidates:
287 cell.removeCandidate(cand)
289 raise RuntimeError(
"All PSF candidates removed as blends")
293 disp = afwDisplay.Display(frame=0)
294 disp.mtv(exposure, title=
"psf determination")
295 utils.showPsfSpatialCells(exposure, psfCellSet, self.config.nStarPerCell, symb=
"o",
296 ctype=afwDisplay.CYAN, ctypeUnused=afwDisplay.YELLOW,
297 size=4, display=disp)
303 for iterNum
in range(self.config.nIterForPsf):
304 if display
and displayPsfCandidates:
307 for cell
in psfCellSet.getCellList():
308 for cand
in cell.begin(
not showBadCandidates):
310 im = cand.getMaskedImage()
312 chi2 = cand.getChi2()
316 stamps.append((im,
"%d%s" %
317 (utils.splitId(cand.getSource().getId(),
True)[
"objId"], chi2),
323 print(
"WARNING: No PSF candidates to show; try setting showBadCandidates=True")
325 mos = afwDisplay.utils.Mosaic()
326 for im, label, status
in stamps:
327 im =
type(im)(im,
True)
330 except NotImplementedError:
333 mos.append(im, label,
334 (afwDisplay.GREEN
if status == afwMath.SpatialCellCandidate.GOOD
else
335 afwDisplay.YELLOW
if status == afwMath.SpatialCellCandidate.UNKNOWN
else
338 disp8 = afwDisplay.Display(frame=8)
339 mos.makeMosaic(display=disp8, title=
"Psf Candidates")
348 psf, eigenValues, nEigenComponents, fitChi2 = \
349 self.
_fitPsf(exposure, psfCellSet, actualKernelSize, nEigenComponents)
354 for cell
in psfCellSet.getCellList():
356 for cand
in cell.begin(
False):
357 cand.setStatus(afwMath.SpatialCellCandidate.UNKNOWN)
358 rchi2 = cand.getChi2()
359 if not numpy.isfinite(rchi2)
or rchi2 <= 0:
361 awfulCandidates.append(cand)
363 self.log.debug(
"chi^2=%s; id=%s",
364 cand.getChi2(), cand.getSource().getId())
365 for cand
in awfulCandidates:
367 print(
"Removing bad candidate: id=%d, chi^2=%f" %
368 (cand.getSource().getId(), cand.getChi2()))
369 cell.removeCandidate(cand)
374 badCandidates =
list()
375 for cell
in psfCellSet.getCellList():
376 for cand
in cell.begin(
False):
377 rchi2 = cand.getChi2()
379 if rchi2 > self.config.reducedChi2ForPsfCandidates:
380 badCandidates.append(cand)
382 badCandidates.sort(key=
lambda x: x.getChi2(), reverse=
True)
384 self.config.nIterForPsf)
385 for i, c
in zip(range(numBad), badCandidates):
391 print(
"Chi^2 clipping %-4d %.2g" % (c.getSource().getId(), chi2))
392 c.setStatus(afwMath.SpatialCellCandidate.BAD)
405 kernel = psf.getKernel()
406 noSpatialKernel = psf.getKernel()
407 for cell
in psfCellSet.getCellList():
408 for cand
in cell.begin(
False):
411 im = cand.getMaskedImage(actualKernelSize, actualKernelSize)
415 fit = fitKernelParamsToImage(noSpatialKernel, im, candCenter)
419 for p, k
in zip(params, kernels):
422 predict = [kernel.getSpatialFunction(k)(candCenter.getX(), candCenter.getY())
for
423 k
in range(kernel.getNKernelParameters())]
425 residuals.append([a/amp - p
for a, p
in zip(params, predict)])
426 candidates.append(cand)
428 residuals = numpy.array(residuals)
430 for k
in range(kernel.getNKernelParameters()):
433 mean = residuals[:, k].mean()
434 rms = residuals[:, k].
std()
437 sr = numpy.sort(residuals[:, k])
438 mean = (sr[int(0.5*len(sr))]
if len(sr)%2
else
439 0.5*(sr[int(0.5*len(sr))] + sr[int(0.5*len(sr)) + 1]))
440 rms = 0.74*(sr[int(0.75*len(sr))] - sr[int(0.25*len(sr))])
443 mean = stats.getValue(afwMath.MEANCLIP)
444 rms = stats.getValue(afwMath.STDEVCLIP)
446 rms =
max(1.0e-4, rms)
449 print(
"Mean for component %d is %f" % (k, mean))
450 print(
"RMS for component %d is %f" % (k, rms))
451 badCandidates =
list()
452 for i, cand
in enumerate(candidates):
453 if numpy.fabs(residuals[i, k] - mean) > self.config.spatialReject*rms:
454 badCandidates.append(i)
456 badCandidates.sort(key=
lambda x: numpy.fabs(residuals[x, k] - mean), reverse=
True)
459 self.config.nIterForPsf)
461 for i, c
in zip(range(
min(len(badCandidates), numBad)), badCandidates):
464 print(
"Spatial clipping %d (%f,%f) based on %d: %f vs %f" %
465 (cand.getSource().getId(), cand.getXCenter(), cand.getYCenter(), k,
466 residuals[badCandidates[i], k], self.config.spatialReject*rms))
467 cand.setStatus(afwMath.SpatialCellCandidate.BAD)
472 if display
and displayIterations:
476 utils.showPsfSpatialCells(exposure, psfCellSet, self.config.nStarPerCell, showChi2=
True,
477 symb=
"o", size=8, display=disp, ctype=afwDisplay.YELLOW,
478 ctypeBad=afwDisplay.RED, ctypeUnused=afwDisplay.MAGENTA)
479 if self.config.nStarPerCellSpatialFit != self.config.nStarPerCell:
480 utils.showPsfSpatialCells(exposure, psfCellSet, self.config.nStarPerCellSpatialFit,
481 symb=
"o", size=10, display=disp,
482 ctype=afwDisplay.YELLOW, ctypeBad=afwDisplay.RED)
486 disp4 = afwDisplay.Display(frame=4)
487 utils.showPsfCandidates(exposure, psfCellSet, psf=psf, display=disp4,
488 normalize=normalizeResiduals,
489 showBadCandidates=showBadCandidates)
490 disp5 = afwDisplay.Display(frame=5)
491 utils.showPsfCandidates(exposure, psfCellSet, psf=psf, display=disp5,
492 normalize=normalizeResiduals,
493 showBadCandidates=showBadCandidates,
496 if not showBadCandidates:
497 showBadCandidates =
True
501 if displayPsfComponents:
502 disp6 = afwDisplay.Display(frame=6)
503 utils.showPsf(psf, eigenValues, display=disp6)
505 disp7 = afwDisplay.Display(frame=7)
506 utils.showPsfMosaic(exposure, psf, display=disp7, showFwhm=
True)
507 disp7.scale(
'linear', 0, 1)
508 if displayPsfSpatialModel:
509 utils.plotPsfSpatialModel(exposure, psf, psfCellSet, showBadCandidates=
True,
510 matchKernelAmplitudes=matchKernelAmplitudes,
511 keepPlots=keepMatplotlibPlots)
516 reply = input(
"Next iteration? [ynchpqQs] ").
strip()
520 reply = reply.split()
522 reply, args = reply[0], reply[1:]
526 if reply
in (
"",
"c",
"h",
"n",
"p",
"q",
"Q",
"s",
"y"):
530 print(
"c[ontinue without prompting] h[elp] n[o] p[db] q[uit displaying] "
531 "s[ave fileName] y[es]")
541 fileName = args.pop(0)
543 print(
"Please provide a filename")
546 print(
"Saving to %s" % fileName)
547 utils.saveSpatialCellSet(psfCellSet, fileName=fileName)
551 print(
"Unrecognised response: %s" % reply, file=sys.stderr)
557 psf, eigenValues, nEigenComponents, fitChi2 = \
558 self.
_fitPsf(exposure, psfCellSet, actualKernelSize, nEigenComponents)
563 if display
and reply !=
"n":
564 disp = afwDisplay.Display(frame=0)
566 utils.showPsfSpatialCells(exposure, psfCellSet, self.config.nStarPerCell, showChi2=
True,
567 symb=
"o", ctype=afwDisplay.YELLOW, ctypeBad=afwDisplay.RED,
568 size=8, display=disp)
569 if self.config.nStarPerCellSpatialFit != self.config.nStarPerCell:
570 utils.showPsfSpatialCells(exposure, psfCellSet, self.config.nStarPerCellSpatialFit,
571 symb=
"o", ctype=afwDisplay.YELLOW, ctypeBad=afwDisplay.RED,
572 size=10, display=disp)
574 disp4 = afwDisplay.Display(frame=4)
575 utils.showPsfCandidates(exposure, psfCellSet, psf=psf, display=disp4,
576 normalize=normalizeResiduals,
577 showBadCandidates=showBadCandidates)
579 if displayPsfComponents:
580 disp6 = afwDisplay.Display(frame=6)
581 utils.showPsf(psf, eigenValues, display=disp6)
584 disp7 = afwDisplay.Display(frame=7)
585 utils.showPsfMosaic(exposure, psf, display=disp7, showFwhm=
True)
586 disp7.scale(
"linear", 0, 1)
587 if displayPsfSpatialModel:
588 utils.plotPsfSpatialModel(exposure, psf, psfCellSet, showBadCandidates=
True,
589 matchKernelAmplitudes=matchKernelAmplitudes,
590 keepPlots=keepMatplotlibPlots)
602 for cell
in psfCellSet.getCellList():
603 for cand
in cell.begin(
False):
606 for cand
in cell.begin(
True):
607 src = cand.getSource()
608 if flagKey
is not None:
609 src.set(flagKey,
True)
617 if metadata
is not None:
618 metadata[
"spatialFitChi2"] = fitChi2
619 metadata[
"numGoodStars"] = numGoodStars
620 metadata[
"numAvailStars"] = numAvailStars
621 metadata[
"avgX"] = avgX
622 metadata[
"avgY"] = avgY
626 return psf, psfCellSet
630 """Generator for Psf candidates.
632 This allows two 'for' loops to be reduced to one.
637 SpatialCellSet of PSF candidates.
638 ignoreBad : `bool`, optional
639 Ignore candidates flagged
as BAD?
648 for cell
in psfCellSet.getCellList():
649 for cand
in cell.begin(ignoreBad):
653psfDeterminerRegistry.register(
"pca", PcaPsfDeterminerTask)
An ellipse core with quadrupole moments as parameters.
A class to contain the data, WCS, and other information needed to describe an image of the sky.
Class to ensure constraints for spatial modeling.
A collection of SpatialCells covering an entire image.
Represent a PSF as a linear combination of PCA (== Karhunen-Loeve) basis functions.
Class stored in SpatialCells for spatial Psf fitting.
determinePsf(self, exposure, psfCandidateList, metadata=None, flagKey=None)
_fitPsf(self, exposure, psfCellSet, kernelSize, nEigenComponents)
Reports attempts to exceed implementation-defined length limits for some classes.
daf::base::PropertyList * list
Statistics makeStatistics(lsst::afw::image::Image< Pixel > const &img, lsst::afw::image::Mask< image::MaskPixel > const &msk, int const flags, StatisticsControl const &sctrl=StatisticsControl())
Handle a watered-down front-end to the constructor (no variance)
numCandidatesToReject(numBadCandidates, numIter, totalIter)
candidatesIter(psfCellSet, ignoreBad=True)