LSST Applications 27.0.0,g0265f82a02+469cd937ee,g02d81e74bb+21ad69e7e1,g1470d8bcf6+cbe83ee85a,g2079a07aa2+e67c6346a6,g212a7c68fe+04a9158687,g2305ad1205+94392ce272,g295015adf3+81dd352a9d,g2bbee38e9b+469cd937ee,g337abbeb29+469cd937ee,g3939d97d7f+72a9f7b576,g487adcacf7+71499e7cba,g50ff169b8f+5929b3527e,g52b1c1532d+a6fc98d2e7,g591dd9f2cf+df404f777f,g5a732f18d5+be83d3ecdb,g64a986408d+21ad69e7e1,g858d7b2824+21ad69e7e1,g8a8a8dda67+a6fc98d2e7,g99cad8db69+f62e5b0af5,g9ddcbc5298+d4bad12328,ga1e77700b3+9c366c4306,ga8c6da7877+71e4819109,gb0e22166c9+25ba2f69a1,gb6a65358fc+469cd937ee,gbb8dafda3b+69d3c0e320,gc07e1c2157+a98bf949bb,gc120e1dc64+615ec43309,gc28159a63d+469cd937ee,gcf0d15dbbd+72a9f7b576,gdaeeff99f8+a38ce5ea23,ge6526c86ff+3a7c1ac5f1,ge79ae78c31+469cd937ee,gee10cc3b42+a6fc98d2e7,gf1cff7945b+21ad69e7e1,gfbcc870c63+9a11dc8c8f
LSST Data Management Base Package
Loading...
Searching...
No Matches
pcaPsfDeterminer.py
Go to the documentation of this file.
1# This file is part of meas_algorithms.
2#
3# Developed for the LSST Data Management System.
4# This product includes software developed by the LSST Project
5# (https://www.lsst.org).
6# See the COPYRIGHT file at the top-level directory of this distribution
7# for details of code ownership.
8#
9# This program is free software: you can redistribute it and/or modify
10# it under the terms of the GNU General Public License as published by
11# the Free Software Foundation, either version 3 of the License, or
12# (at your option) any later version.
13#
14# This program is distributed in the hope that it will be useful,
15# but WITHOUT ANY WARRANTY; without even the implied warranty of
16# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17# GNU General Public License for more details.
18#
19# You should have received a copy of the GNU General Public License
20# along with this program. If not, see <https://www.gnu.org/licenses/>.
21
22__all__ = ["PcaPsfDeterminerConfig", "PcaPsfDeterminerTask"]
23
24import sys
25
26import numpy
27
28import lsst.pex.config as pexConfig
29import lsst.pex.exceptions as pexExceptions
30import lsst.geom
31import lsst.afw.geom as afwGeom
32import lsst.afw.geom.ellipses as afwEll
33import lsst.afw.display as afwDisplay
34import lsst.afw.math as afwMath
35from .psfDeterminer import BasePsfDeterminerTask, psfDeterminerRegistry
36from ._algorithmsLib import PsfCandidateF
37from ._algorithmsLib import createKernelFromPsfCandidates, countPsfCandidates, \
38 fitSpatialKernelFromPsfCandidates, fitKernelParamsToImage
39from ._algorithmsLib import PcaPsf
40from . import utils
41
42
43def numCandidatesToReject(numBadCandidates, numIter, totalIter):
44 """Return the number of PSF candidates to be rejected.
45
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
48 candidates.
49
50 Parameters
51 ----------
52 numBadCandidates : `int`
53 Number of bad candidates under consideration.
54
55 numIter : `int`
56 The number of the current PSF iteration.
57
58 totalIter : `int`
59 The total number of PSF iterations.
60
61 Returns
62 -------
63 return : `int`
64 Number of candidates to reject.
65 """
66 return int(numBadCandidates*(numIter + 1)//totalIter + 0.5)
67
68
69class PcaPsfDeterminerConfig(BasePsfDeterminerTask.ConfigClass):
70 nonLinearSpatialFit = pexConfig.Field[bool](
71 doc="Use non-linear fitter for spatial variation of Kernel",
72 default=False,
73 )
74 nEigenComponents = pexConfig.Field[int](
75 doc="number of eigen components for PSF kernel creation",
76 default=4,
77 check=lambda x: x >= 1
78 )
79 spatialOrder = pexConfig.Field[int](
80 doc="specify spatial order for PSF kernel creation",
81 default=2,
82 )
83 sizeCellX = pexConfig.Field[int](
84 doc="size of cell used to determine PSF (pixels, column direction)",
85 default=256,
86 # minValue = 10,
87 check=lambda x: x >= 10,
88 )
89 sizeCellY = pexConfig.Field[int](
90 doc="size of cell used to determine PSF (pixels, row direction)",
91 default=sizeCellX.default,
92 # minValue = 10,
93 check=lambda x: x >= 10,
94 )
95 nStarPerCell = pexConfig.Field[int](
96 doc="number of stars per psf cell for PSF kernel creation",
97 default=3,
98 )
99 borderWidth = pexConfig.Field[int](
100 doc="Number of pixels to ignore around the edge of PSF candidate postage stamps",
101 default=0,
102 )
103 nStarPerCellSpatialFit = pexConfig.Field[int](
104 doc="number of stars per psf Cell for spatial fitting",
105 default=5,
106 )
107 constantWeight = pexConfig.Field[bool](
108 doc="Should each PSF candidate be given the same weight, independent of magnitude?",
109 default=True,
110 )
111 nIterForPsf = pexConfig.Field[int](
112 doc="number of iterations of PSF candidate star list",
113 default=3,
114 )
115 tolerance = pexConfig.Field[float](
116 doc="tolerance of spatial fitting",
117 default=1e-2,
118 )
119 lam = pexConfig.Field[float](
120 doc="floor for variance is lam*data",
121 default=0.05,
122 )
123 reducedChi2ForPsfCandidates = pexConfig.Field[float](
124 doc="for psf candidate evaluation",
125 default=2.0,
126 )
127 spatialReject = pexConfig.Field[float](
128 doc="Rejection threshold (stdev) for candidates based on spatial fit",
129 default=3.0,
130 )
131 pixelThreshold = pexConfig.Field[float](
132 doc="Threshold (stdev) for rejecting extraneous pixels around candidate; applied if positive",
133 default=0.0,
134 )
135 doRejectBlends = pexConfig.Field[bool](
136 doc="Reject candidates that are blended?",
137 default=False,
138 )
139 doMaskBlends = pexConfig.Field[bool](
140 doc="Mask blends in image?",
141 default=True,
142 )
143
144 def setDefaults(self):
145 super().setDefaults()
146 self.stampSize = 41
147
148
150 """A measurePsfTask psf estimator.
151 """
152 ConfigClass = PcaPsfDeterminerConfig
153
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
169 except pexExceptions.LengthError as e:
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
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
220 displayIterations = lsstDebug.Info(__name__).displayIterations # display on each PSF iteration
221 displayPsfComponents = lsstDebug.Info(__name__).displayPsfComponents # show the PCA components
222 displayResiduals = lsstDebug.Info(__name__).displayResiduals # show residuals
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
228 displayPsfSpatialModel = lsstDebug.Info(__name__).displayPsfSpatialModel # Plot spatial model?
229 showBadCandidates = lsstDebug.Info(__name__).showBadCandidates # Include bad candidates
230 # Normalize residuals by object amplitude
231 normalizeResiduals = lsstDebug.Info(__name__).normalizeResiduals
232 pause = lsstDebug.Info(__name__).pause # Prompt user after each iteration?
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
448 rms = max(1.0e-4, rms) # Don't trust RMS below this due to numerical issues
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 # 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:
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 # 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
631def candidatesIter(psfCellSet, ignoreBad=True):
632 """Generator for Psf candidates.
633
634 This allows two 'for' loops to be reduced to one.
635
636 Parameters
637 ----------
638 psfCellSet : `lsst.afw.math.SpatialCellSet`
639 SpatialCellSet of PSF candidates.
640 ignoreBad : `bool`, optional
641 Ignore candidates flagged as BAD?
642
643 Yields
644 -------
645 cell : `lsst.afw.math.SpatialCell`
646 A SpatialCell.
647 cand : `lsst.meas.algorithms.PsfCandidate`
648 A PsfCandidate.
649 """
650 for cell in psfCellSet.getCellList():
651 for cand in cell.begin(ignoreBad):
652 yield (cell, cand)
653
654
655psfDeterminerRegistry.register("pca", PcaPsfDeterminerTask)
int min
int max
An ellipse core with quadrupole moments as parameters.
Definition Quadrupole.h:47
A collection of SpatialCells covering an entire image.
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.
Definition Runtime.h:76
bool strip
Definition fits.cc:930
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
numCandidatesToReject(numBadCandidates, numIter, totalIter)
candidatesIter(psfCellSet, ignoreBad=True)
STL namespace.