LSST Applications 26.0.0,g0265f82a02+6660c170cc,g07994bdeae+30b05a742e,g0a0026dc87+17526d298f,g0a60f58ba1+17526d298f,g0e4bf8285c+96dd2c2ea9,g0ecae5effc+c266a536c8,g1e7d6db67d+6f7cb1f4bb,g26482f50c6+6346c0633c,g2bbee38e9b+6660c170cc,g2cc88a2952+0a4e78cd49,g3273194fdb+f6908454ef,g337abbeb29+6660c170cc,g337c41fc51+9a8f8f0815,g37c6e7c3d5+7bbafe9d37,g44018dc512+6660c170cc,g4a941329ef+4f7594a38e,g4c90b7bd52+5145c320d2,g58be5f913a+bea990ba40,g635b316a6c+8d6b3a3e56,g67924a670a+bfead8c487,g6ae5381d9b+81bc2a20b4,g93c4d6e787+26b17396bd,g98cecbdb62+ed2cb6d659,g98ffbb4407+81bc2a20b4,g9ddcbc5298+7f7571301f,ga1e77700b3+99e9273977,gae46bcf261+6660c170cc,gb2715bf1a1+17526d298f,gc86a011abf+17526d298f,gcf0d15dbbd+96dd2c2ea9,gdaeeff99f8+0d8dbea60f,gdb4ec4c597+6660c170cc,ge23793e450+96dd2c2ea9,gf041782ebf+171108ac67
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 -------
210 The measured PSF.
211 psfCellSet : `lsst.afw.math.SpatialCellSet`
212 The PSF candidates.
213 """
214 import lsstDebug
215 display = lsstDebug.Info(__name__).display
216 displayExposure = lsstDebug.Info(__name__).displayExposure # display the Exposure + spatialCells
217 displayPsfCandidates = lsstDebug.Info(__name__).displayPsfCandidates # show the viable candidates
218 displayIterations = lsstDebug.Info(__name__).displayIterations # display on each PSF iteration
219 displayPsfComponents = lsstDebug.Info(__name__).displayPsfComponents # show the PCA components
220 displayResiduals = lsstDebug.Info(__name__).displayResiduals # show residuals
221 displayPsfMosaic = lsstDebug.Info(__name__).displayPsfMosaic # show mosaic of reconstructed PSF(x,y)
222 # match Kernel amplitudes for spatial plots
223 matchKernelAmplitudes = lsstDebug.Info(__name__).matchKernelAmplitudes
224 # Keep matplotlib alive post mortem
225 keepMatplotlibPlots = lsstDebug.Info(__name__).keepMatplotlibPlots
226 displayPsfSpatialModel = lsstDebug.Info(__name__).displayPsfSpatialModel # Plot spatial model?
227 showBadCandidates = lsstDebug.Info(__name__).showBadCandidates # Include bad candidates
228 # Normalize residuals by object amplitude
229 normalizeResiduals = lsstDebug.Info(__name__).normalizeResiduals
230 pause = lsstDebug.Info(__name__).pause # Prompt user after each iteration?
231
232 if display:
233 afwDisplay.setDefaultMaskTransparency(75)
234 if display > 1:
235 pause = True
236
237 mi = exposure.getMaskedImage()
238
239 if len(psfCandidateList) == 0:
240 raise RuntimeError("No PSF candidates supplied.")
241
242 # construct and populate a spatial cell set
243 bbox = mi.getBBox()
244 psfCellSet = afwMath.SpatialCellSet(bbox, self.config.sizeCellX, self.config.sizeCellY)
245 sizes = []
246 for i, psfCandidate in enumerate(psfCandidateList):
247 if psfCandidate.getSource().getPsfFluxFlag(): # bad measurement
248 continue
249
250 try:
251 psfCellSet.insertCandidate(psfCandidate)
252 except Exception as e:
253 self.log.debug("Skipping PSF candidate %d of %d: %s", i, len(psfCandidateList), e)
254 continue
255 source = psfCandidate.getSource()
256
257 quad = afwGeom.Quadrupole(source.getIxx(), source.getIyy(), source.getIxy())
258 axes = afwEll.Axes(quad)
259 sizes.append(axes.getA())
260 if len(sizes) == 0:
261 raise RuntimeError("No usable PSF candidates supplied")
262 nEigenComponents = self.config.nEigenComponents # initial version
263
264 actualKernelSize = int(self.config.stampSize)
265
266 if display:
267 print("Median size=%s" % (numpy.median(sizes),))
268
269 self.log.trace("Kernel size=%s", actualKernelSize)
270
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())
276
277 if self.config.doRejectBlends:
278 # Remove blended candidates completely
279 blendedCandidates = [] # Candidates to remove; can't do it while iterating
280 for cell, cand in candidatesIter(psfCellSet, False):
281 if len(cand.getSource().getFootprint().getPeaks()) > 1:
282 blendedCandidates.append((cell, cand))
283 continue
284 if display:
285 print("Removing %d blended Psf candidates" % len(blendedCandidates))
286 for cell, cand in blendedCandidates:
287 cell.removeCandidate(cand)
288 if sum(1 for cand in candidatesIter(psfCellSet, False)) == 0:
289 raise RuntimeError("All PSF candidates removed as blends")
290
291 if display:
292 if displayExposure:
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)
298
299 #
300 # Do a PCA decomposition of those PSF candidates
301 #
302 reply = "y" # used in interactive mode
303 for iterNum in range(self.config.nIterForPsf):
304 if display and displayPsfCandidates: # Show a mosaic of usable PSF candidates
305
306 stamps = []
307 for cell in psfCellSet.getCellList():
308 for cand in cell.begin(not showBadCandidates): # maybe include bad candidates
309 try:
310 im = cand.getMaskedImage()
311
312 chi2 = cand.getChi2()
313 if chi2 > 1e100:
314 chi2 = numpy.nan
315
316 stamps.append((im, "%d%s" %
317 (utils.splitId(cand.getSource().getId(), True)["objId"], chi2),
318 cand.getStatus()))
319 except Exception:
320 continue
321
322 if len(stamps) == 0:
323 print("WARNING: No PSF candidates to show; try setting showBadCandidates=True")
324 else:
325 mos = afwDisplay.utils.Mosaic()
326 for im, label, status in stamps:
327 im = type(im)(im, True)
328 try:
329 im /= afwMath.makeStatistics(im, afwMath.MAX).getValue()
330 except NotImplementedError:
331 pass
332
333 mos.append(im, label,
334 (afwDisplay.GREEN if status == afwMath.SpatialCellCandidate.GOOD else
335 afwDisplay.YELLOW if status == afwMath.SpatialCellCandidate.UNKNOWN else
336 afwDisplay.RED))
337
338 disp8 = afwDisplay.Display(frame=8)
339 mos.makeMosaic(display=disp8, title="Psf Candidates")
340
341 # Re-fit until we don't have any candidates with naughty chi^2 values influencing the fit
342 cleanChi2 = False # Any naughty (negative/NAN) chi^2 values?
343 while not cleanChi2:
344 cleanChi2 = True
345 #
346 # First, estimate the PSF
347 #
348 psf, eigenValues, nEigenComponents, fitChi2 = \
349 self._fitPsf(exposure, psfCellSet, actualKernelSize, nEigenComponents)
350 #
351 # In clipping, allow all candidates to be innocent until proven guilty on this iteration.
352 # Throw out any prima facie guilty candidates (naughty chi^2 values)
353 #
354 for cell in psfCellSet.getCellList():
355 awfulCandidates = []
356 for cand in cell.begin(False): # include bad candidates
357 cand.setStatus(afwMath.SpatialCellCandidate.UNKNOWN) # until proven guilty
358 rchi2 = cand.getChi2()
359 if not numpy.isfinite(rchi2) or rchi2 <= 0:
360 # Guilty prima facie
361 awfulCandidates.append(cand)
362 cleanChi2 = False
363 self.log.debug("chi^2=%s; id=%s",
364 cand.getChi2(), cand.getSource().getId())
365 for cand in awfulCandidates:
366 if display:
367 print("Removing bad candidate: id=%d, chi^2=%f" %
368 (cand.getSource().getId(), cand.getChi2()))
369 cell.removeCandidate(cand)
370
371 #
372 # Clip out bad fits based on reduced chi^2
373 #
374 badCandidates = list()
375 for cell in psfCellSet.getCellList():
376 for cand in cell.begin(False): # include bad candidates
377 rchi2 = cand.getChi2() # reduced chi^2 when fitting PSF to candidate
378 assert rchi2 > 0
379 if rchi2 > self.config.reducedChi2ForPsfCandidates:
380 badCandidates.append(cand)
381
382 badCandidates.sort(key=lambda x: x.getChi2(), reverse=True)
383 numBad = numCandidatesToReject(len(badCandidates), iterNum,
384 self.config.nIterForPsf)
385 for i, c in zip(range(numBad), badCandidates):
386 if display:
387 chi2 = c.getChi2()
388 if chi2 > 1e100:
389 chi2 = numpy.nan
390
391 print("Chi^2 clipping %-4d %.2g" % (c.getSource().getId(), chi2))
392 c.setStatus(afwMath.SpatialCellCandidate.BAD)
393
394 #
395 # Clip out bad fits based on spatial fitting.
396 #
397 # This appears to be better at getting rid of sources that have a single dominant kernel component
398 # (other than the zeroth; e.g., a nearby contaminant) because the surrounding sources (which help
399 # set the spatial model) don't contain that kernel component, and so the spatial modeling
400 # downweights the component.
401 #
402
403 residuals = list()
404 candidates = list()
405 kernel = psf.getKernel()
406 noSpatialKernel = psf.getKernel()
407 for cell in psfCellSet.getCellList():
408 for cand in cell.begin(False):
409 candCenter = lsst.geom.PointD(cand.getXCenter(), cand.getYCenter())
410 try:
411 im = cand.getMaskedImage(actualKernelSize, actualKernelSize)
412 except Exception:
413 continue
414
415 fit = fitKernelParamsToImage(noSpatialKernel, im, candCenter)
416 params = fit[0]
417 kernels = fit[1]
418 amp = 0.0
419 for p, k in zip(params, kernels):
420 amp += p*k.getSum()
421
422 predict = [kernel.getSpatialFunction(k)(candCenter.getX(), candCenter.getY()) for
423 k in range(kernel.getNKernelParameters())]
424
425 residuals.append([a/amp - p for a, p in zip(params, predict)])
426 candidates.append(cand)
427
428 residuals = numpy.array(residuals)
429
430 for k in range(kernel.getNKernelParameters()):
431 if False:
432 # Straight standard deviation
433 mean = residuals[:, k].mean()
434 rms = residuals[:, k].std()
435 elif False:
436 # Using interquartile range
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))])
441 else:
442 stats = afwMath.makeStatistics(residuals[:, k], afwMath.MEANCLIP | afwMath.STDEVCLIP)
443 mean = stats.getValue(afwMath.MEANCLIP)
444 rms = stats.getValue(afwMath.STDEVCLIP)
445
446 rms = max(1.0e-4, rms) # Don't trust RMS below this due to numerical issues
447
448 if display:
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)
455
456 badCandidates.sort(key=lambda x: numpy.fabs(residuals[x, k] - mean), reverse=True)
457
458 numBad = numCandidatesToReject(len(badCandidates), iterNum,
459 self.config.nIterForPsf)
460
461 for i, c in zip(range(min(len(badCandidates), numBad)), badCandidates):
462 cand = candidates[c]
463 if display:
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)
468
469 #
470 # Display results
471 #
472 if display and displayIterations:
473 if displayExposure:
474 if iterNum > 0:
475 disp.erase()
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)
483 if displayResiduals:
484 while True:
485 try:
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,
494 variance=True)
495 except Exception:
496 if not showBadCandidates:
497 showBadCandidates = True
498 continue
499 break
500
501 if displayPsfComponents:
502 disp6 = afwDisplay.Display(frame=6)
503 utils.showPsf(psf, eigenValues, display=disp6)
504 if displayPsfMosaic:
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)
512
513 if pause:
514 while True:
515 try:
516 reply = input("Next iteration? [ynchpqQs] ").strip()
517 except EOFError:
518 reply = "n"
519
520 reply = reply.split()
521 if reply:
522 reply, args = reply[0], reply[1:]
523 else:
524 reply = ""
525
526 if reply in ("", "c", "h", "n", "p", "q", "Q", "s", "y"):
527 if reply == "c":
528 pause = False
529 elif reply == "h":
530 print("c[ontinue without prompting] h[elp] n[o] p[db] q[uit displaying] "
531 "s[ave fileName] y[es]")
532 continue
533 elif reply == "p":
534 import pdb
535 pdb.set_trace()
536 elif reply == "q":
537 display = False
538 elif reply == "Q":
539 sys.exit(1)
540 elif reply == "s":
541 fileName = args.pop(0)
542 if not fileName:
543 print("Please provide a filename")
544 continue
545
546 print("Saving to %s" % fileName)
547 utils.saveSpatialCellSet(psfCellSet, fileName=fileName)
548 continue
549 break
550 else:
551 print("Unrecognised response: %s" % reply, file=sys.stderr)
552
553 if reply == "n":
554 break
555
556 # One last time, to take advantage of the last iteration
557 psf, eigenValues, nEigenComponents, fitChi2 = \
558 self._fitPsf(exposure, psfCellSet, actualKernelSize, nEigenComponents)
559
560 #
561 # Display code for debugging
562 #
563 if display and reply != "n":
564 disp = afwDisplay.Display(frame=0)
565 if displayExposure:
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)
573 if displayResiduals:
574 disp4 = afwDisplay.Display(frame=4)
575 utils.showPsfCandidates(exposure, psfCellSet, psf=psf, display=disp4,
576 normalize=normalizeResiduals,
577 showBadCandidates=showBadCandidates)
578
579 if displayPsfComponents:
580 disp6 = afwDisplay.Display(frame=6)
581 utils.showPsf(psf, eigenValues, display=disp6)
582
583 if displayPsfMosaic:
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)
591 #
592 # Generate some QA information
593 #
594 # Count PSF stars
595 #
596 numGoodStars = 0
597 numAvailStars = 0
598
599 avgX = 0.0
600 avgY = 0.0
601
602 for cell in psfCellSet.getCellList():
603 for cand in cell.begin(False): # don't ignore BAD stars
604 numAvailStars += 1
605
606 for cand in cell.begin(True): # do ignore BAD stars
607 src = cand.getSource()
608 if flagKey is not None:
609 src.set(flagKey, True)
610 avgX += src.getX()
611 avgY += src.getY()
612 numGoodStars += 1
613
614 avgX /= numGoodStars
615 avgY /= numGoodStars
616
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
623
624 psf = PcaPsf(psf.getKernel(), lsst.geom.Point2D(avgX, avgY))
625
626 return psf, psfCellSet
627
628
629def candidatesIter(psfCellSet, ignoreBad=True):
630 """Generator for Psf candidates.
631
632 This allows two 'for' loops to be reduced to one.
633
634 Parameters
635 ----------
636 psfCellSet : `lsst.afw.math.SpatialCellSet`
637 SpatialCellSet of PSF candidates.
638 ignoreBad : `bool`, optional
639 Ignore candidates flagged as BAD?
640
641 Yields
642 -------
644 A SpatialCell.
646 A PsfCandidate.
647 """
648 for cell in psfCellSet.getCellList():
649 for cand in cell.begin(ignoreBad):
650 yield (cell, cand)
651
652
653psfDeterminerRegistry.register("pca", PcaPsfDeterminerTask)
int min
int max
table::Key< int > type
Definition Detector.cc:163
An ellipse core with quadrupole moments as parameters.
Definition Quadrupole.h:47
A class to contain the data, WCS, and other information needed to describe an image of the sky.
Definition Exposure.h:72
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.
Definition PcaPsf.h:37
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.
Definition Runtime.h:76
daf::base::PropertyList * list
Definition fits.cc:928
bool strip
Definition fits.cc:926
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.