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.
194
195 Parameters
196 ----------
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.
206
207 Returns
208 -------
209 psf : `lsst.meas.algorithms.PcaPsf`
210 The measured PSF.
211 psfCellSet : `lsst.afw.math.SpatialCellSet`
212 The PSF candidates.
213 """
214 psfCandidateList = self.downsampleCandidates(psfCandidateList)
215
216 import lsstDebug
219 displayPsfCandidates =
lsstDebug.Info(__name__).displayPsfCandidates
221 displayPsfComponents =
lsstDebug.Info(__name__).displayPsfComponents
224
225 matchKernelAmplitudes =
lsstDebug.Info(__name__).matchKernelAmplitudes
226
227 keepMatplotlibPlots =
lsstDebug.Info(__name__).keepMatplotlibPlots
228 displayPsfSpatialModel =
lsstDebug.Info(__name__).displayPsfSpatialModel
230
233
234 if display:
235 afwDisplay.setDefaultMaskTransparency(75)
236 if display > 1:
237 pause = True
238
239 mi = exposure.getMaskedImage()
240
241 if len(psfCandidateList) == 0:
242 raise RuntimeError("No PSF candidates supplied.")
243
244
245 bbox = mi.getBBox()
247 sizes = []
248 for i, psfCandidate in enumerate(psfCandidateList):
249 if psfCandidate.getSource().getPsfFluxFlag():
250 continue
251
252 try:
253 psfCellSet.insertCandidate(psfCandidate)
254 except Exception as e:
255 self.log.debug("Skipping PSF candidate %d of %d: %s", i, len(psfCandidateList), e)
256 continue
257 source = psfCandidate.getSource()
258
260 axes = afwEll.Axes(quad)
261 sizes.append(axes.getA())
262 if len(sizes) == 0:
263 raise RuntimeError("No usable PSF candidates supplied")
264 nEigenComponents = self.config.nEigenComponents
265
266 actualKernelSize = int(self.config.stampSize)
267
268 if display:
269 print("Median size=%s" % (numpy.median(sizes),))
270
271 self.log.trace("Kernel size=%s", actualKernelSize)
272
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())
278
279 if self.config.doRejectBlends:
280
281 blendedCandidates = []
282 for cell, cand in candidatesIter(psfCellSet, False):
283 if len(cand.getSource().getFootprint().getPeaks()) > 1:
284 blendedCandidates.append((cell, cand))
285 continue
286 if display:
287 print("Removing %d blended Psf candidates" % len(blendedCandidates))
288 for cell, cand in blendedCandidates:
289 cell.removeCandidate(cand)
290 if sum(1 for cand in candidatesIter(psfCellSet, False)) == 0:
291 raise RuntimeError("All PSF candidates removed as blends")
292
293 if display:
294 if displayExposure:
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)
300
301
302
303
304 reply = "y"
305 for iterNum in range(self.config.nIterForPsf):
306 if display and displayPsfCandidates:
307
308 stamps = []
309 for cell in psfCellSet.getCellList():
310 for cand in cell.begin(not showBadCandidates):
311 try:
312 im = cand.getMaskedImage()
313
314 chi2 = cand.getChi2()
315 if chi2 > 1e100:
316 chi2 = numpy.nan
317
318 stamps.append((im, "%d%s" %
319 (utils.splitId(cand.getSource().getId(), True)["objId"], chi2),
320 cand.getStatus()))
321 except Exception:
322 continue
323
324 if len(stamps) == 0:
325 print("WARNING: No PSF candidates to show; try setting showBadCandidates=True")
326 else:
327 mos = afwDisplay.utils.Mosaic()
328 for im, label, status in stamps:
329 im = type(im)(im, True)
330 try:
332 except NotImplementedError:
333 pass
334
335 mos.append(im, label,
336 (afwDisplay.GREEN if status == afwMath.SpatialCellCandidate.GOOD else
337 afwDisplay.YELLOW if status == afwMath.SpatialCellCandidate.UNKNOWN else
338 afwDisplay.RED))
339
340 disp8 = afwDisplay.Display(frame=8)
341 mos.makeMosaic(display=disp8, title="Psf Candidates")
342
343
344 cleanChi2 = False
345 while not cleanChi2:
346 cleanChi2 = True
347
348
349
350 psf, eigenValues, nEigenComponents, fitChi2 = \
351 self._fitPsf(exposure, psfCellSet, actualKernelSize, nEigenComponents)
352
353
354
355
356 for cell in psfCellSet.getCellList():
357 awfulCandidates = []
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:
362
363 awfulCandidates.append(cand)
364 cleanChi2 = False
365 self.log.debug("chi^2=%s; id=%s",
366 cand.getChi2(), cand.getSource().getId())
367 for cand in awfulCandidates:
368 if display:
369 print("Removing bad candidate: id=%d, chi^2=%f" %
370 (cand.getSource().getId(), cand.getChi2()))
371 cell.removeCandidate(cand)
372
373
374
375
376 badCandidates = list()
377 for cell in psfCellSet.getCellList():
378 for cand in cell.begin(False):
379 rchi2 = cand.getChi2()
380 assert rchi2 > 0
381 if rchi2 > self.config.reducedChi2ForPsfCandidates:
382 badCandidates.append(cand)
383
384 badCandidates.sort(key=lambda x: x.getChi2(), reverse=True)
385 numBad = numCandidatesToReject(len(badCandidates), iterNum,
386 self.config.nIterForPsf)
387 for i, c in zip(range(numBad), badCandidates):
388 if display:
389 chi2 = c.getChi2()
390 if chi2 > 1e100:
391 chi2 = numpy.nan
392
393 print("Chi^2 clipping %-4d %.2g" % (c.getSource().getId(), chi2))
394 c.setStatus(afwMath.SpatialCellCandidate.BAD)
395
396
397
398
399
400
401
402
403
404
405 residuals = list()
406 candidates = list()
407 kernel = psf.getKernel()
408 noSpatialKernel = psf.getKernel()
409 for cell in psfCellSet.getCellList():
410 for cand in cell.begin(False):
412 try:
413 im = cand.getMaskedImage(actualKernelSize, actualKernelSize)
414 except Exception:
415 continue
416
417 fit = fitKernelParamsToImage(noSpatialKernel, im, candCenter)
418 params = fit[0]
419 kernels = fit[1]
420 amp = 0.0
421 for p, k in zip(params, kernels):
422 amp += p*k.getSum()
423
424 predict = [kernel.getSpatialFunction(k)(candCenter.getX(), candCenter.getY()) for
425 k in range(kernel.getNKernelParameters())]
426
427 residuals.append([a/amp - p for a, p in zip(params, predict)])
428 candidates.append(cand)
429
430 residuals = numpy.array(residuals)
431
432 for k in range(kernel.getNKernelParameters()):
433 if False:
434
435 mean = residuals[:, k].mean()
436 rms = residuals[:, k].
std()
437 elif False:
438
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))])
443 else:
445 mean = stats.getValue(afwMath.MEANCLIP)
446 rms = stats.getValue(afwMath.STDEVCLIP)
447
448 rms = max(1.0e-4, rms)
449
450 if display:
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)
457
458 badCandidates.sort(key=lambda x: numpy.fabs(residuals[x, k] - mean), reverse=True)
459
460 numBad = numCandidatesToReject(len(badCandidates), iterNum,
461 self.config.nIterForPsf)
462
463 for i, c in zip(range(min(len(badCandidates), numBad)), badCandidates):
464 cand = candidates[c]
465 if display:
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)
470
471
472
473
474 if display and displayIterations:
475 if displayExposure:
476 if iterNum > 0:
477 disp.erase()
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)
485 if displayResiduals:
486 while True:
487 try:
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,
496 variance=True)
497 except Exception:
498 if not showBadCandidates:
499 showBadCandidates = True
500 continue
501 break
502
503 if displayPsfComponents:
504 disp6 = afwDisplay.Display(frame=6)
505 utils.showPsf(psf, eigenValues, display=disp6)
506 if displayPsfMosaic:
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)
514
515 if pause:
516 while True:
517 try:
518 reply = input("Next iteration? [ynchpqQs] ").strip()
519 except EOFError:
520 reply = "n"
521
522 reply = reply.split()
523 if reply:
524 reply, args = reply[0], reply[1:]
525 else:
526 reply = ""
527
528 if reply in ("", "c", "h", "n", "p", "q", "Q", "s", "y"):
529 if reply == "c":
530 pause = False
531 elif reply == "h":
532 print("c[ontinue without prompting] h[elp] n[o] p[db] q[uit displaying] "
533 "s[ave fileName] y[es]")
534 continue
535 elif reply == "p":
536 import pdb
537 pdb.set_trace()
538 elif reply == "q":
539 display = False
540 elif reply == "Q":
541 sys.exit(1)
542 elif reply == "s":
543 fileName = args.pop(0)
544 if not fileName:
545 print("Please provide a filename")
546 continue
547
548 print("Saving to %s" % fileName)
549 utils.saveSpatialCellSet(psfCellSet, fileName=fileName)
550 continue
551 break
552 else:
553 print("Unrecognised response: %s" % reply, file=sys.stderr)
554
555 if reply == "n":
556 break
557
558
559 psf, eigenValues, nEigenComponents, fitChi2 = \
560 self._fitPsf(exposure, psfCellSet, actualKernelSize, nEigenComponents)
561
562
563
564
565 if display and reply != "n":
566 disp = afwDisplay.Display(frame=0)
567 if displayExposure:
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)
575 if displayResiduals:
576 disp4 = afwDisplay.Display(frame=4)
577 utils.showPsfCandidates(exposure, psfCellSet, psf=psf, display=disp4,
578 normalize=normalizeResiduals,
579 showBadCandidates=showBadCandidates)
580
581 if displayPsfComponents:
582 disp6 = afwDisplay.Display(frame=6)
583 utils.showPsf(psf, eigenValues, display=disp6)
584
585 if displayPsfMosaic:
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)
593
594
595
596
597
598 numGoodStars = 0
599 numAvailStars = 0
600
601 avgX = 0.0
602 avgY = 0.0
603
604 for cell in psfCellSet.getCellList():
605 for cand in cell.begin(False):
606 numAvailStars += 1
607
608 for cand in cell.begin(True):
609 src = cand.getSource()
610 if flagKey is not None:
611 src.set(flagKey, True)
612 avgX += src.getX()
613 avgY += src.getY()
614 numGoodStars += 1
615
616 avgX /= numGoodStars
617 avgY /= numGoodStars
618
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
625
627
628 return psf, psfCellSet
629
630
An ellipse core with quadrupole moments as parameters.
A collection of SpatialCells covering an entire image.
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)