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.
197 exposure : `lsst.afw.image.Exposure`
198 Exposure containing the psf candidates.
199 psfCandidateList : `list` of `lsst.meas.algorithms.PsfCandidate`
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.
209 psf : `lsst.meas.algorithms.PcaPsf`
211 psfCellSet : `lsst.afw.math.SpatialCellSet`
219 displayPsfCandidates =
lsstDebug.Info(__name__).displayPsfCandidates
221 displayPsfComponents =
lsstDebug.Info(__name__).displayPsfComponents
225 matchKernelAmplitudes =
lsstDebug.Info(__name__).matchKernelAmplitudes
227 keepMatplotlibPlots =
lsstDebug.Info(__name__).keepMatplotlibPlots
228 displayPsfSpatialModel =
lsstDebug.Info(__name__).displayPsfSpatialModel
235 afwDisplay.setDefaultMaskTransparency(75)
239 mi = exposure.getMaskedImage()
241 if len(psfCandidateList) == 0:
242 raise RuntimeError(
"No PSF candidates supplied.")
248 for i, psfCandidate
in enumerate(psfCandidateList):
249 if psfCandidate.getSource().getPsfFluxFlag():
253 psfCellSet.insertCandidate(psfCandidate)
254 except Exception
as e:
255 self.log.debug(
"Skipping PSF candidate %d of %d: %s", i, len(psfCandidateList), e)
257 source = psfCandidate.getSource()
260 axes = afwEll.Axes(quad)
261 sizes.append(axes.getA())
263 raise RuntimeError(
"No usable PSF candidates supplied")
264 nEigenComponents = self.config.nEigenComponents
266 actualKernelSize = int(self.config.stampSize)
269 print(
"Median size=%s" % (numpy.median(sizes),))
271 self.log.trace(
"Kernel size=%s", actualKernelSize)
273 if actualKernelSize > psfCandidateList[0].getWidth():
274 self.log.warning(
"Using a region (%d x %d) larger than kernelSize (%d) set while making PSF "
275 "candidates. Consider setting a larger value for kernelSize for "
276 "`makePsfCandidates` to avoid this warning.",
277 actualKernelSize, actualKernelSize, psfCandidateList[0].getWidth())
279 if self.config.doRejectBlends:
281 blendedCandidates = []
283 if len(cand.getSource().getFootprint().getPeaks()) > 1:
284 blendedCandidates.append((cell, cand))
287 print(
"Removing %d blended Psf candidates" % len(blendedCandidates))
288 for cell, cand
in blendedCandidates:
289 cell.removeCandidate(cand)
291 raise RuntimeError(
"All PSF candidates removed as blends")
295 disp = afwDisplay.Display(frame=0)
296 disp.mtv(exposure, title=
"psf determination")
297 utils.showPsfSpatialCells(exposure, psfCellSet, self.config.nStarPerCell, symb=
"o",
298 ctype=afwDisplay.CYAN, ctypeUnused=afwDisplay.YELLOW,
299 size=4, display=disp)
305 for iterNum
in range(self.config.nIterForPsf):
306 if display
and displayPsfCandidates:
309 for cell
in psfCellSet.getCellList():
310 for cand
in cell.begin(
not showBadCandidates):
312 im = cand.getMaskedImage()
314 chi2 = cand.getChi2()
318 stamps.append((im,
"%d%s" %
319 (utils.splitId(cand.getSource().getId(),
True)[
"objId"], chi2),
325 print(
"WARNING: No PSF candidates to show; try setting showBadCandidates=True")
327 mos = afwDisplay.utils.Mosaic()
328 for im, label, status
in stamps:
329 im = type(im)(im,
True)
332 except NotImplementedError:
335 mos.append(im, label,
336 (afwDisplay.GREEN
if status == afwMath.SpatialCellCandidate.GOOD
else
337 afwDisplay.YELLOW
if status == afwMath.SpatialCellCandidate.UNKNOWN
else
340 disp8 = afwDisplay.Display(frame=8)
341 mos.makeMosaic(display=disp8, title=
"Psf Candidates")
350 psf, eigenValues, nEigenComponents, fitChi2 = \
351 self.
_fitPsf(exposure, psfCellSet, actualKernelSize, nEigenComponents)
356 for cell
in psfCellSet.getCellList():
358 for cand
in cell.begin(
False):
359 cand.setStatus(afwMath.SpatialCellCandidate.UNKNOWN)
360 rchi2 = cand.getChi2()
361 if not numpy.isfinite(rchi2)
or rchi2 <= 0:
363 awfulCandidates.append(cand)
365 self.log.debug(
"chi^2=%s; id=%s",
366 cand.getChi2(), cand.getSource().getId())
367 for cand
in awfulCandidates:
369 print(
"Removing bad candidate: id=%d, chi^2=%f" %
370 (cand.getSource().getId(), cand.getChi2()))
371 cell.removeCandidate(cand)
376 badCandidates = list()
377 for cell
in psfCellSet.getCellList():
378 for cand
in cell.begin(
False):
379 rchi2 = cand.getChi2()
381 if rchi2 > self.config.reducedChi2ForPsfCandidates:
382 badCandidates.append(cand)
384 badCandidates.sort(key=
lambda x: x.getChi2(), reverse=
True)
386 self.config.nIterForPsf)
387 for i, c
in zip(range(numBad), badCandidates):
393 print(
"Chi^2 clipping %-4d %.2g" % (c.getSource().getId(), chi2))
394 c.setStatus(afwMath.SpatialCellCandidate.BAD)
407 kernel = psf.getKernel()
408 noSpatialKernel = psf.getKernel()
409 for cell
in psfCellSet.getCellList():
410 for cand
in cell.begin(
False):
413 im = cand.getMaskedImage(actualKernelSize, actualKernelSize)
417 fit = fitKernelParamsToImage(noSpatialKernel, im, candCenter)
421 for p, k
in zip(params, kernels):
424 predict = [kernel.getSpatialFunction(k)(candCenter.getX(), candCenter.getY())
for
425 k
in range(kernel.getNKernelParameters())]
427 residuals.append([a/amp - p
for a, p
in zip(params, predict)])
428 candidates.append(cand)
430 residuals = numpy.array(residuals)
432 for k
in range(kernel.getNKernelParameters()):
435 mean = residuals[:, k].mean()
436 rms = residuals[:, k].
std()
439 sr = numpy.sort(residuals[:, k])
440 mean = (sr[int(0.5*len(sr))]
if len(sr)%2
else
441 0.5*(sr[int(0.5*len(sr))] + sr[int(0.5*len(sr)) + 1]))
442 rms = 0.74*(sr[int(0.75*len(sr))] - sr[int(0.25*len(sr))])
445 mean = stats.getValue(afwMath.MEANCLIP)
446 rms = stats.getValue(afwMath.STDEVCLIP)
448 rms =
max(1.0e-4, rms)
451 print(
"Mean for component %d is %f" % (k, mean))
452 print(
"RMS for component %d is %f" % (k, rms))
453 badCandidates = list()
454 for i, cand
in enumerate(candidates):
455 if numpy.fabs(residuals[i, k] - mean) > self.config.spatialReject*rms:
456 badCandidates.append(i)
458 badCandidates.sort(key=
lambda x: numpy.fabs(residuals[x, k] - mean), reverse=
True)
461 self.config.nIterForPsf)
463 for i, c
in zip(range(
min(len(badCandidates), numBad)), badCandidates):
466 print(
"Spatial clipping %d (%f,%f) based on %d: %f vs %f" %
467 (cand.getSource().getId(), cand.getXCenter(), cand.getYCenter(), k,
468 residuals[badCandidates[i], k], self.config.spatialReject*rms))
469 cand.setStatus(afwMath.SpatialCellCandidate.BAD)
474 if display
and displayIterations:
478 utils.showPsfSpatialCells(exposure, psfCellSet, self.config.nStarPerCell, showChi2=
True,
479 symb=
"o", size=8, display=disp, ctype=afwDisplay.YELLOW,
480 ctypeBad=afwDisplay.RED, ctypeUnused=afwDisplay.MAGENTA)
481 if self.config.nStarPerCellSpatialFit != self.config.nStarPerCell:
482 utils.showPsfSpatialCells(exposure, psfCellSet, self.config.nStarPerCellSpatialFit,
483 symb=
"o", size=10, display=disp,
484 ctype=afwDisplay.YELLOW, ctypeBad=afwDisplay.RED)
488 disp4 = afwDisplay.Display(frame=4)
489 utils.showPsfCandidates(exposure, psfCellSet, psf=psf, display=disp4,
490 normalize=normalizeResiduals,
491 showBadCandidates=showBadCandidates)
492 disp5 = afwDisplay.Display(frame=5)
493 utils.showPsfCandidates(exposure, psfCellSet, psf=psf, display=disp5,
494 normalize=normalizeResiduals,
495 showBadCandidates=showBadCandidates,
498 if not showBadCandidates:
499 showBadCandidates =
True
503 if displayPsfComponents:
504 disp6 = afwDisplay.Display(frame=6)
505 utils.showPsf(psf, eigenValues, display=disp6)
507 disp7 = afwDisplay.Display(frame=7)
508 utils.showPsfMosaic(exposure, psf, display=disp7, showFwhm=
True)
509 disp7.scale(
'linear', 0, 1)
510 if displayPsfSpatialModel:
511 utils.plotPsfSpatialModel(exposure, psf, psfCellSet, showBadCandidates=
True,
512 matchKernelAmplitudes=matchKernelAmplitudes,
513 keepPlots=keepMatplotlibPlots)
518 reply = input(
"Next iteration? [ynchpqQs] ").
strip()
522 reply = reply.split()
524 reply, args = reply[0], reply[1:]
528 if reply
in (
"",
"c",
"h",
"n",
"p",
"q",
"Q",
"s",
"y"):
532 print(
"c[ontinue without prompting] h[elp] n[o] p[db] q[uit displaying] "
533 "s[ave fileName] y[es]")
543 fileName = args.pop(0)
545 print(
"Please provide a filename")
548 print(
"Saving to %s" % fileName)
549 utils.saveSpatialCellSet(psfCellSet, fileName=fileName)
553 print(
"Unrecognised response: %s" % reply, file=sys.stderr)
559 psf, eigenValues, nEigenComponents, fitChi2 = \
560 self.
_fitPsf(exposure, psfCellSet, actualKernelSize, nEigenComponents)
565 if display
and reply !=
"n":
566 disp = afwDisplay.Display(frame=0)
568 utils.showPsfSpatialCells(exposure, psfCellSet, self.config.nStarPerCell, showChi2=
True,
569 symb=
"o", ctype=afwDisplay.YELLOW, ctypeBad=afwDisplay.RED,
570 size=8, display=disp)
571 if self.config.nStarPerCellSpatialFit != self.config.nStarPerCell:
572 utils.showPsfSpatialCells(exposure, psfCellSet, self.config.nStarPerCellSpatialFit,
573 symb=
"o", ctype=afwDisplay.YELLOW, ctypeBad=afwDisplay.RED,
574 size=10, display=disp)
576 disp4 = afwDisplay.Display(frame=4)
577 utils.showPsfCandidates(exposure, psfCellSet, psf=psf, display=disp4,
578 normalize=normalizeResiduals,
579 showBadCandidates=showBadCandidates)
581 if displayPsfComponents:
582 disp6 = afwDisplay.Display(frame=6)
583 utils.showPsf(psf, eigenValues, display=disp6)
586 disp7 = afwDisplay.Display(frame=7)
587 utils.showPsfMosaic(exposure, psf, display=disp7, showFwhm=
True)
588 disp7.scale(
"linear", 0, 1)
589 if displayPsfSpatialModel:
590 utils.plotPsfSpatialModel(exposure, psf, psfCellSet, showBadCandidates=
True,
591 matchKernelAmplitudes=matchKernelAmplitudes,
592 keepPlots=keepMatplotlibPlots)
604 for cell
in psfCellSet.getCellList():
605 for cand
in cell.begin(
False):
608 for cand
in cell.begin(
True):
609 src = cand.getSource()
610 if flagKey
is not None:
611 src.set(flagKey,
True)
619 if metadata
is not None:
620 metadata[
"spatialFitChi2"] = fitChi2
621 metadata[
"numGoodStars"] = numGoodStars
622 metadata[
"numAvailStars"] = numAvailStars
623 metadata[
"avgX"] = avgX
624 metadata[
"avgY"] = avgY
628 return psf, psfCellSet