Inheritance diagram for lsst.meas.algorithms.pcaPsfDeterminer.PcaPsfDeterminerTask:
Public Member Functions | |
| determinePsf (self, exposure, psfCandidateList, metadata=None, flagKey=None) | |
Static Public Attributes | |
| ConfigClass = PcaPsfDeterminerConfig | |
Protected Member Functions | |
| _fitPsf (self, exposure, psfCellSet, kernelSize, nEigenComponents) | |
Detailed Description
A measurePsfTask psf estimator.
Definition at line 149 of file pcaPsfDeterminer.py.
Member Function Documentation
◆ _fitPsf()
|
protected |
Definition at line 154 of file pcaPsfDeterminer.py.
154 def _fitPsf(self, exposure, psfCellSet, kernelSize, nEigenComponents):
155 PsfCandidateF.setPixelThreshold(self.config.pixelThreshold)
156 PsfCandidateF.setMaskBlends(self.config.doMaskBlends)
157 #
158 # Loop trying to use nEigenComponents, but allowing smaller numbers if necessary
159 #
160 for nEigen in range(nEigenComponents, 0, -1):
161 # Determine KL components
162 try:
163 kernel, eigenValues = createKernelFromPsfCandidates(
164 psfCellSet, exposure.getDimensions(), exposure.getXY0(), nEigen,
165 self.config.spatialOrder, kernelSize, self.config.nStarPerCell,
166 bool(self.config.constantWeight))
167
168 break # OK, we can get nEigen components
170 if nEigen == 1: # can't go any lower
171 raise IndexError("No viable PSF candidates survive")
172
173 self.log.warning("%s: reducing number of eigen components", e.what())
174 #
175 # We got our eigen decomposition so let's use it
176 #
177 # Express eigenValues in units of reduced chi^2 per star
178 size = kernelSize + 2*self.config.borderWidth
179 nu = size*size - 1 # number of degrees of freedom/star for chi^2
180 eigenValues = [val/float(countPsfCandidates(psfCellSet, self.config.nStarPerCell)*nu)
181 for val in eigenValues]
182
183 # Fit spatial model
184 status, chi2 = fitSpatialKernelFromPsfCandidates(
185 kernel, psfCellSet, bool(self.config.nonLinearSpatialFit),
186 self.config.nStarPerCellSpatialFit, self.config.tolerance, self.config.lam)
187
188 psf = PcaPsf(kernel)
189
190 return psf, eigenValues, nEigen, chi2
191
Reports attempts to exceed implementation-defined length limits for some classes.
Definition Runtime.h:76
◆ determinePsf()
| 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.
Reimplemented from lsst.meas.algorithms.psfDeterminer.BasePsfDeterminerTask.
Definition at line 192 of file pcaPsfDeterminer.py.
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
217 display = lsstDebug.Info(__name__).display
218 displayExposure = lsstDebug.Info(__name__).displayExposure # display the Exposure + spatialCells
219 displayPsfCandidates = lsstDebug.Info(__name__).displayPsfCandidates # show the viable candidates
223 displayPsfMosaic = lsstDebug.Info(__name__).displayPsfMosaic # show mosaic of reconstructed PSF(x,y)
224 # match Kernel amplitudes for spatial plots
225 matchKernelAmplitudes = lsstDebug.Info(__name__).matchKernelAmplitudes
226 # Keep matplotlib alive post mortem
227 keepMatplotlibPlots = lsstDebug.Info(__name__).keepMatplotlibPlots
230 # Normalize residuals by object amplitude
231 normalizeResiduals = lsstDebug.Info(__name__).normalizeResiduals
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 # construct and populate a spatial cell set
245 bbox = mi.getBBox()
246 psfCellSet = afwMath.SpatialCellSet(bbox, self.config.sizeCellX, self.config.sizeCellY)
247 sizes = []
248 for i, psfCandidate in enumerate(psfCandidateList):
249 if psfCandidate.getSource().getPsfFluxFlag(): # bad measurement
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
259 quad = afwGeom.Quadrupole(source.getIxx(), source.getIyy(), source.getIxy())
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 # initial version
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 # Remove blended candidates completely
281 blendedCandidates = [] # Candidates to remove; can't do it while iterating
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 # Do a PCA decomposition of those PSF candidates
303 #
304 reply = "y" # used in interactive mode
305 for iterNum in range(self.config.nIterForPsf):
306 if display and displayPsfCandidates: # Show a mosaic of usable PSF candidates
307
308 stamps = []
309 for cell in psfCellSet.getCellList():
310 for cand in cell.begin(not showBadCandidates): # maybe include bad candidates
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:
331 im /= afwMath.makeStatistics(im, afwMath.MAX).getValue()
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 # Re-fit until we don't have any candidates with naughty chi^2 values influencing the fit
344 cleanChi2 = False # Any naughty (negative/NAN) chi^2 values?
345 while not cleanChi2:
346 cleanChi2 = True
347 #
348 # First, estimate the PSF
349 #
350 psf, eigenValues, nEigenComponents, fitChi2 = \
351 self._fitPsf(exposure, psfCellSet, actualKernelSize, nEigenComponents)
352 #
353 # In clipping, allow all candidates to be innocent until proven guilty on this iteration.
354 # Throw out any prima facie guilty candidates (naughty chi^2 values)
355 #
356 for cell in psfCellSet.getCellList():
357 awfulCandidates = []
358 for cand in cell.begin(False): # include bad candidates
359 cand.setStatus(afwMath.SpatialCellCandidate.UNKNOWN) # until proven guilty
360 rchi2 = cand.getChi2()
361 if not numpy.isfinite(rchi2) or rchi2 <= 0:
362 # Guilty prima facie
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 # Clip out bad fits based on reduced chi^2
375 #
376 badCandidates = list()
377 for cell in psfCellSet.getCellList():
378 for cand in cell.begin(False): # include bad candidates
379 rchi2 = cand.getChi2() # reduced chi^2 when fitting PSF to candidate
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 # Clip out bad fits based on spatial fitting.
398 #
399 # This appears to be better at getting rid of sources that have a single dominant kernel component
400 # (other than the zeroth; e.g., a nearby contaminant) because the surrounding sources (which help
401 # set the spatial model) don't contain that kernel component, and so the spatial modeling
402 # downweights the component.
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):
411 candCenter = lsst.geom.PointD(cand.getXCenter(), cand.getYCenter())
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 # Straight standard deviation
435 mean = residuals[:, k].mean()
436 rms = residuals[:, k].std()
437 elif False:
438 # Using interquartile range
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:
444 stats = afwMath.makeStatistics(residuals[:, k], afwMath.MEANCLIP | afwMath.STDEVCLIP)
445 mean = stats.getValue(afwMath.MEANCLIP)
446 rms = stats.getValue(afwMath.STDEVCLIP)
447
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
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 # Display results
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:
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 # One last time, to take advantage of the last iteration
559 psf, eigenValues, nEigenComponents, fitChi2 = \
560 self._fitPsf(exposure, psfCellSet, actualKernelSize, nEigenComponents)
561
562 #
563 # Display code for debugging
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 # Generate some QA information
595 #
596 # Count PSF stars
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): # don't ignore BAD stars
606 numAvailStars += 1
607
608 for cand in cell.begin(True): # do ignore BAD stars
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
626 psf = PcaPsf(psf.getKernel(), lsst.geom.Point2D(avgX, avgY))
627
628 return psf, psfCellSet
629
630
An ellipse core with quadrupole moments as parameters.
Definition Quadrupole.h:47
A collection of SpatialCells covering an entire image.
Definition SpatialCell.h:383
Definition lsstDebug.py:28
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)
Definition Statistics.h:361
STL namespace.
Member Data Documentation
◆ ConfigClass
|
static |
Definition at line 152 of file pcaPsfDeterminer.py.
The documentation for this class was generated from the following file:
- /j/snowflake/release/lsstsw/stack/lsst-scipipe-8.0.0/Linux64/meas_algorithms/ge6526c86ff+0433e6603d/python/lsst/meas/algorithms/pcaPsfDeterminer.py
