LSSTApplications  10.0-2-g4f67435,11.0.rc2+1,11.0.rc2+12,11.0.rc2+3,11.0.rc2+4,11.0.rc2+5,11.0.rc2+6,11.0.rc2+7,11.0.rc2+8
LSSTDataManagementBasePackage
pcaPsfDeterminer.py
Go to the documentation of this file.
1 #
2 # LSST Data Management System
3 # Copyright 2008, 2009, 2010 LSST Corporation.
4 #
5 # This product includes software developed by the
6 # LSST Project (http://www.lsst.org/).
7 #
8 # This program is free software: you can redistribute it and/or modify
9 # it under the terms of the GNU General Public License as published by
10 # the Free Software Foundation, either version 3 of the License, or
11 # (at your option) any later version.
12 #
13 # This program is distributed in the hope that it will be useful,
14 # but WITHOUT ANY WARRANTY; without even the implied warranty of
15 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 # GNU General Public License for more details.
17 #
18 # You should have received a copy of the LSST License Statement and
19 # the GNU General Public License along with this program. If not,
20 # see <http://www.lsstcorp.org/LegalNotices/>.
21 #
22 import math
23 import sys
24 
25 import numpy
26 
27 import lsst.pex.config as pexConfig
28 import lsst.pex.exceptions as pexExceptions
29 import lsst.pex.logging as pexLog
30 import lsst.afw.geom as afwGeom
31 import lsst.afw.geom.ellipses as afwEll
32 import lsst.afw.display.ds9 as ds9
33 import lsst.afw.math as afwMath
34 from . import algorithmsLib
35 from . import utils as maUtils
36 
37 class PcaPsfDeterminerConfig(pexConfig.Config):
38  nonLinearSpatialFit = pexConfig.Field(
39  doc = "Use non-linear fitter for spatial variation of Kernel",
40  dtype = bool,
41  default = False,
42  )
43  nEigenComponents = pexConfig.Field(
44  doc = "number of eigen components for PSF kernel creation",
45  dtype = int,
46  default = 4,
47  )
48  spatialOrder = pexConfig.Field(
49  doc = "specify spatial order for PSF kernel creation",
50  dtype = int,
51  default = 2,
52  )
53  sizeCellX = pexConfig.Field(
54  doc = "size of cell used to determine PSF (pixels, column direction)",
55  dtype = int,
56  default = 256,
57 # minValue = 10,
58  check = lambda x: x >= 10,
59  )
60  sizeCellY = pexConfig.Field(
61  doc = "size of cell used to determine PSF (pixels, row direction)",
62  dtype = int,
63  default = sizeCellX.default,
64 # minValue = 10,
65  check = lambda x: x >= 10,
66  )
67  nStarPerCell = pexConfig.Field(
68  doc = "number of stars per psf cell for PSF kernel creation",
69  dtype = int,
70  default = 3,
71  )
72  kernelSize = pexConfig.Field(
73  doc = "radius of the kernel to create, relative to the square root of the stellar quadrupole moments",
74  dtype = float,
75  default = 10.0,
76  )
77  kernelSizeMin = pexConfig.Field(
78  doc = "Minimum radius of the kernel",
79  dtype = int,
80  default = 25,
81  )
82  kernelSizeMax = pexConfig.Field(
83  doc = "Maximum radius of the kernel",
84  dtype = int,
85  default = 45,
86  )
87  borderWidth = pexConfig.Field(
88  doc = "Number of pixels to ignore around the edge of PSF candidate postage stamps",
89  dtype = int,
90  default = 0,
91  )
92  nStarPerCellSpatialFit = pexConfig.Field(
93  doc = "number of stars per psf Cell for spatial fitting",
94  dtype = int,
95  default = 5,
96  )
97  constantWeight = pexConfig.Field(
98  doc = "Should each PSF candidate be given the same weight, independent of magnitude?",
99  dtype = bool,
100  default = True,
101  )
102  nIterForPsf = pexConfig.Field(
103  doc = "number of iterations of PSF candidate star list",
104  dtype = int,
105  default = 3,
106  )
107  tolerance = pexConfig.Field(
108  doc = "tolerance of spatial fitting",
109  dtype = float,
110  default = 1e-2,
111  )
112  lam = pexConfig.Field(
113  doc = "floor for variance is lam*data",
114  dtype = float,
115  default = 0.05,
116  )
117  reducedChi2ForPsfCandidates = pexConfig.Field(
118  doc = "for psf candidate evaluation",
119  dtype = float,
120  default = 2.0,
121  )
122  spatialReject = pexConfig.Field(
123  doc = "Rejection threshold (stdev) for candidates based on spatial fit",
124  dtype = float,
125  default = 3.0,
126  )
127  pixelThreshold = pexConfig.Field(
128  doc = "Threshold (stdev) for rejecting extraneous pixels around candidate; applied if positive",
129  dtype = float,
130  default = 0.0,
131  )
132  doRejectBlends = pexConfig.Field(
133  doc = "Reject candidates that are blended?",
134  dtype = bool,
135  default = False,
136  )
137  doMaskBlends = pexConfig.Field(
138  doc = "Mask blends in image?",
139  dtype = bool,
140  default = True,
141  )
142 
143 class PcaPsfDeterminer(object):
144  """!
145  A measurePsfTask psf estimator
146  """
147  ConfigClass = PcaPsfDeterminerConfig
148 
149  def __init__(self, config):
150  """!Construct a PCA PSF Fitter
151 
152  \param[in] config instance of PcaPsfDeterminerConfig
153  """
154  self.config = config
155  # N.b. name of component is meas.algorithms.psfDeterminer so you can turn on psf debugging
156  # independent of which determiner is active
157  self.debugLog = pexLog.Debug("meas.algorithms.psfDeterminer")
158  self.warnLog = pexLog.Log(pexLog.getDefaultLog(), "meas.algorithms.psfDeterminer")
159 
160  def _fitPsf(self, exposure, psfCellSet, kernelSize, nEigenComponents):
161  algorithmsLib.PsfCandidateF.setPixelThreshold(self.config.pixelThreshold)
162  algorithmsLib.PsfCandidateF.setMaskBlends(self.config.doMaskBlends)
163  #
164  # Loop trying to use nEigenComponents, but allowing smaller numbers if necessary
165  #
166  for nEigen in range(nEigenComponents, 0, -1):
167  # Determine KL components
168  try:
169  kernel, eigenValues = algorithmsLib.createKernelFromPsfCandidates(
170  psfCellSet, exposure.getDimensions(), exposure.getXY0(), nEigen,
171  self.config.spatialOrder, kernelSize, self.config.nStarPerCell,
172  bool(self.config.constantWeight))
173 
174  break # OK, we can get nEigen components
175  except pexExceptions.LengthError as e:
176  if nEigen == 1: # can't go any lower
177  raise
178 
179  self.warnLog.log(pexLog.Log.WARN, "%s: reducing number of eigen components" % e.what())
180  #
181  # We got our eigen decomposition so let's use it
182  #
183  # Express eigenValues in units of reduced chi^2 per star
184  size = kernelSize + 2*self.config.borderWidth
185  nu = size*size - 1 # number of degrees of freedom/star for chi^2
186  eigenValues = [l/float(algorithmsLib.countPsfCandidates(psfCellSet, self.config.nStarPerCell)*nu)
187  for l in eigenValues]
188 
189  # Fit spatial model
190  status, chi2 = algorithmsLib.fitSpatialKernelFromPsfCandidates(
191  kernel, psfCellSet, bool(self.config.nonLinearSpatialFit),
192  self.config.nStarPerCellSpatialFit, self.config.tolerance, self.config.lam)
193 
194  psf = algorithmsLib.PcaPsf(kernel)
195 
196  return psf, eigenValues, nEigen, chi2
197 
198 
199  def determinePsf(self, exposure, psfCandidateList, metadata=None, flagKey=None):
200  """!Determine a PCA PSF model for an exposure given a list of PSF candidates
201 
202  \param[in] exposure exposure containing the psf candidates (lsst.afw.image.Exposure)
203  \param[in] psfCandidateList a sequence of PSF candidates (each an lsst.meas.algorithms.PsfCandidate);
204  typically obtained by detecting sources and then running them through a star selector
205  \param[in,out] metadata a home for interesting tidbits of information
206  \param[in] flagKey schema key used to mark sources actually used in PSF determination
207 
208  \return a list of
209  - psf: the measured PSF, an lsst.meas.algorithms.PcaPsf
210  - cellSet: an lsst.afw.math.SpatialCellSet containing the PSF candidates
211  """
212  import lsstDebug
213  display = lsstDebug.Info(__name__).display
214  displayExposure = lsstDebug.Info(__name__).displayExposure # display the Exposure + spatialCells
215  displayPsfCandidates = lsstDebug.Info(__name__).displayPsfCandidates # show the viable candidates
216  displayIterations = lsstDebug.Info(__name__).displayIterations # display on each PSF iteration
217  displayPsfComponents = lsstDebug.Info(__name__).displayPsfComponents # show the PCA components
218  displayResiduals = lsstDebug.Info(__name__).displayResiduals # show residuals
219  displayPsfMosaic = lsstDebug.Info(__name__).displayPsfMosaic # show mosaic of reconstructed PSF(x,y)
220  matchKernelAmplitudes = lsstDebug.Info(__name__).matchKernelAmplitudes # match Kernel amplitudes for spatial plots
221  keepMatplotlibPlots = lsstDebug.Info(__name__).keepMatplotlibPlots # Keep matplotlib alive post mortem
222  displayPsfSpatialModel = lsstDebug.Info(__name__).displayPsfSpatialModel # Plot spatial model?
223  showBadCandidates = lsstDebug.Info(__name__).showBadCandidates # Include bad candidates
224  normalizeResiduals = lsstDebug.Info(__name__).normalizeResiduals # Normalise residuals by object amplitude
225  pause = lsstDebug.Info(__name__).pause # Prompt user after each iteration?
226 
227  if display > 1:
228  pause = True
229 
230  mi = exposure.getMaskedImage()
231 
232  if len(psfCandidateList) == 0:
233  raise RuntimeError("No PSF candidates supplied.")
234 
235  # construct and populate a spatial cell set
236  bbox = mi.getBBox()
237  psfCellSet = afwMath.SpatialCellSet(bbox, self.config.sizeCellX, self.config.sizeCellY)
238  sizes = []
239  for i, psfCandidate in enumerate(psfCandidateList):
240  try:
241  psfCellSet.insertCandidate(psfCandidate)
242  except Exception, e:
243  self.debugLog.debug(2, "Skipping PSF candidate %d of %d: %s" % (i, len(psfCandidateList), e))
244  continue
245  source = psfCandidate.getSource()
246 
247  quad = afwEll.Quadrupole(source.getIxx(), source.getIyy(), source.getIxy())
248  axes = afwEll.Axes(quad)
249  sizes.append(axes.getA())
250  if len(sizes) == 0:
251  raise RuntimeError("No usable PSF candidates supplied")
252  nEigenComponents = self.config.nEigenComponents # initial version
253 
254  if self.config.kernelSize >= 15:
255  self.debugLog.debug(1, \
256  "WARNING: NOT scaling kernelSize by stellar quadrupole moment " +
257  "because config.kernelSize=%s >= 15; using config.kernelSize as as the width, instead" \
258  % (self.config.kernelSize,)
259  )
260  actualKernelSize = int(self.config.kernelSize)
261  else:
262  medSize = numpy.median(sizes)
263  actualKernelSize = 2 * int(self.config.kernelSize * math.sqrt(medSize) + 0.5) + 1
264  if actualKernelSize < self.config.kernelSizeMin:
265  actualKernelSize = self.config.kernelSizeMin
266  if actualKernelSize > self.config.kernelSizeMax:
267  actualKernelSize = self.config.kernelSizeMax
268 
269  if display:
270  print "Median size=%s" % (medSize,)
271  self.debugLog.debug(3, "Kernel size=%s" % (actualKernelSize,))
272 
273  # Set size of image returned around candidate
274  psfCandidateList[0].setHeight(actualKernelSize)
275  psfCandidateList[0].setWidth(actualKernelSize)
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  frame = 0
293  if displayExposure:
294  ds9.mtv(exposure, frame=frame, title="psf determination")
295  maUtils.showPsfSpatialCells(exposure, psfCellSet, self.config.nStarPerCell,
296  symb="o", ctype=ds9.CYAN, ctypeUnused=ds9.YELLOW,
297  size=4, frame=frame)
298 
299  #
300  # Do a PCA decomposition of those PSF candidates
301  #
302  reply = "y" # used in interactive mode
303  for iter in range(self.config.nIterForPsf):
304  if display and displayPsfCandidates: # Show a mosaic of usable PSF candidates
305  #
306  import lsst.afw.display.utils as displayUtils
307 
308  stamps = []
309  for cell in psfCellSet.getCellList():
310  for cand in cell.begin(not showBadCandidates): # maybe include bad candidates
311  cand = algorithmsLib.cast_PsfCandidateF(cand)
312 
313  try:
314  im = cand.getMaskedImage()
315 
316  chi2 = cand.getChi2()
317  if chi2 > 1e100:
318  chi2 = numpy.nan
319 
320  stamps.append((im, "%d%s" %
321  (maUtils.splitId(cand.getSource().getId(), True)["objId"], chi2),
322  cand.getStatus()))
323  except Exception, e:
324  continue
325 
326  mos = displayUtils.Mosaic()
327  for im, label, status in stamps:
328  im = type(im)(im, True)
329  try:
330  im /= afwMath.makeStatistics(im, afwMath.MAX).getValue()
331  except NotImplementedError:
332  pass
333 
334  mos.append(im, label,
335  ds9.GREEN if status == afwMath.SpatialCellCandidate.GOOD else
336  ds9.YELLOW if status == afwMath.SpatialCellCandidate.UNKNOWN else ds9.RED)
337 
338 
339  if mos.images:
340  mos.makeMosaic(frame=8, title="Psf Candidates")
341 
342  # Re-fit until we don't have any candidates with naughty chi^2 values influencing the fit
343  cleanChi2 = False # Any naughty (negative/NAN) chi^2 values?
344  while not cleanChi2:
345  cleanChi2 = True
346  #
347  # First, estimate the PSF
348  #
349  psf, eigenValues, nEigenComponents, fitChi2 = \
350  self._fitPsf(exposure, psfCellSet, actualKernelSize, nEigenComponents)
351  #
352  # In clipping, allow all candidates to be innocent until proven guilty on this iteration.
353  # Throw out any prima facie guilty candidates (naughty chi^2 values)
354  #
355  for cell in psfCellSet.getCellList():
356  awfulCandidates = []
357  for cand in cell.begin(False): # include bad candidates
358  cand = algorithmsLib.cast_PsfCandidateF(cand)
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.debugLog.debug(2, "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  cand = algorithmsLib.cast_PsfCandidateF(cand)
380  rchi2 = cand.getChi2() # reduced chi^2 when fitting PSF to candidate
381  assert rchi2 > 0
382  if rchi2 > self.config.reducedChi2ForPsfCandidates:
383  badCandidates.append(cand)
384 
385  badCandidates.sort(key=lambda x: x.getChi2(), reverse=True)
386  numBad = int(len(badCandidates) * (iter + 1) / self.config.nIterForPsf + 0.5)
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 = afwMath.cast_LinearCombinationKernel(psf.getKernel())
409  for cell in psfCellSet.getCellList():
410  for cand in cell.begin(False):
411  cand = algorithmsLib.cast_PsfCandidateF(cand)
412  candCenter = afwGeom.PointD(cand.getXCenter(), cand.getYCenter())
413  try:
414  im = cand.getMaskedImage(kernel.getWidth(), kernel.getHeight())
415  except Exception, e:
416  continue
417 
418  fit = algorithmsLib.fitKernelParamsToImage(noSpatialKernel, im, candCenter)
419  params = fit[0]
420  kernels = fit[1]
421  amp = 0.0
422  for p, k in zip(params, kernels):
423  amp += p * afwMath.cast_FixedKernel(k).getSum()
424 
425  predict = [kernel.getSpatialFunction(k)(candCenter.getX(), candCenter.getY()) for
426  k in range(kernel.getNKernelParameters())]
427 
428  #print cand.getSource().getId(), [a / amp for a in params], predict
429 
430  residuals.append([a / amp - p for a, p in zip(params, predict)])
431  candidates.append(cand)
432 
433  residuals = numpy.array(residuals)
434 
435  for k in range(kernel.getNKernelParameters()):
436  if False:
437  # Straight standard deviation
438  mean = residuals[:,k].mean()
439  rms = residuals[:,k].std()
440  elif False:
441  # Using interquartile range
442  sr = numpy.sort(residuals[:,k])
443  mean = sr[int(0.5*len(sr))] if len(sr) % 2 else \
444  0.5 * (sr[int(0.5*len(sr))] + sr[int(0.5*len(sr))+1])
445  rms = 0.74 * (sr[int(0.75*len(sr))] - sr[int(0.25*len(sr))])
446  else:
447  stats = afwMath.makeStatistics(residuals[:,k], afwMath.MEANCLIP | afwMath.STDEVCLIP)
448  mean = stats.getValue(afwMath.MEANCLIP)
449  rms = stats.getValue(afwMath.STDEVCLIP)
450 
451  rms = max(1.0e-4, rms) # Don't trust RMS below this due to numerical issues
452 
453  if display:
454  print "Mean for component %d is %f" % (k, mean)
455  print "RMS for component %d is %f" % (k, rms)
456  badCandidates = list()
457  for i, cand in enumerate(candidates):
458  if numpy.fabs(residuals[i,k] - mean) > self.config.spatialReject * rms:
459  badCandidates.append(i)
460 
461  badCandidates.sort(key=lambda x: numpy.fabs(residuals[x,k] - mean), reverse=True)
462 
463  numBad = int(len(badCandidates) * (iter + 1) / self.config.nIterForPsf + 0.5)
464 
465  for i, c in zip(range(min(len(badCandidates), numBad)), badCandidates):
466  cand = candidates[c]
467  if display:
468  print "Spatial clipping %d (%f,%f) based on %d: %f vs %f" % \
469  (cand.getSource().getId(), cand.getXCenter(), cand.getYCenter(), k,
470  residuals[badCandidates[i],k], self.config.spatialReject * rms)
471  cand.setStatus(afwMath.SpatialCellCandidate.BAD)
472 
473  #
474  # Display results
475  #
476  if display and displayIterations:
477  if displayExposure:
478  if iter > 0:
479  ds9.erase(frame=frame)
480  maUtils.showPsfSpatialCells(exposure, psfCellSet, self.config.nStarPerCell, showChi2=True,
481  symb="o", size=8, frame=frame,
482  ctype=ds9.YELLOW, ctypeBad=ds9.RED, ctypeUnused=ds9.MAGENTA)
483  if self.config.nStarPerCellSpatialFit != self.config.nStarPerCell:
484  maUtils.showPsfSpatialCells(exposure, psfCellSet, self.config.nStarPerCellSpatialFit,
485  symb="o", size=10, frame=frame,
486  ctype=ds9.YELLOW, ctypeBad=ds9.RED)
487  if displayResiduals:
488  while True:
489  try:
490  maUtils.showPsfCandidates(exposure, psfCellSet, psf=psf, frame=4,
491  normalize=normalizeResiduals,
492  showBadCandidates=showBadCandidates)
493  maUtils.showPsfCandidates(exposure, psfCellSet, psf=psf, frame=5,
494  normalize=normalizeResiduals,
495  showBadCandidates=showBadCandidates,
496  variance=True)
497  except:
498  if not showBadCandidates:
499  showBadCandidates = True
500  continue
501  break
502 
503  if displayPsfComponents:
504  maUtils.showPsf(psf, eigenValues, frame=6)
505  if displayPsfMosaic:
506  maUtils.showPsfMosaic(exposure, psf, frame=7, showFwhm=True)
507  ds9.scale('linear', 0, 1, frame=7)
508  if displayPsfSpatialModel:
509  maUtils.plotPsfSpatialModel(exposure, psf, psfCellSet, showBadCandidates=True,
510  matchKernelAmplitudes=matchKernelAmplitudes,
511  keepPlots=keepMatplotlibPlots)
512 
513  if pause:
514  while True:
515  try:
516  reply = raw_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] s[ave fileName] y[es]"
531  continue
532  elif reply == "p":
533  import pdb; pdb.set_trace()
534  elif reply == "q":
535  display = False
536  elif reply == "Q":
537  sys.exit(1)
538  elif reply == "s":
539  fileName = args.pop(0)
540  if not fileName:
541  print "Please provide a filename"
542  continue
543 
544  print "Saving to %s" % fileName
545  maUtils.saveSpatialCellSet(psfCellSet, fileName=fileName)
546  continue
547  break
548  else:
549  print >> sys.stderr, "Unrecognised response: %s" % reply
550 
551  if reply == "n":
552  break
553 
554  # One last time, to take advantage of the last iteration
555  psf, eigenValues, nEigenComponents, fitChi2 = \
556  self._fitPsf(exposure, psfCellSet, actualKernelSize, nEigenComponents)
557 
558  #
559  # Display code for debugging
560  #
561  if display and reply != "n":
562  if displayExposure:
563  maUtils.showPsfSpatialCells(exposure, psfCellSet, self.config.nStarPerCell, showChi2=True,
564  symb="o", ctype=ds9.YELLOW, ctypeBad=ds9.RED, size=8, frame=frame)
565  if self.config.nStarPerCellSpatialFit != self.config.nStarPerCell:
566  maUtils.showPsfSpatialCells(exposure, psfCellSet, self.config.nStarPerCellSpatialFit,
567  symb="o", ctype=ds9.YELLOW, ctypeBad=ds9.RED,
568  size=10, frame=frame)
569  if displayResiduals:
570  maUtils.showPsfCandidates(exposure, psfCellSet, psf=psf, frame=4,
571  normalize=normalizeResiduals,
572  showBadCandidates=showBadCandidates)
573 
574  if displayPsfComponents:
575  maUtils.showPsf(psf, eigenValues, frame=6)
576 
577  if displayPsfMosaic:
578  maUtils.showPsfMosaic(exposure, psf, frame=7, showFwhm=True)
579  ds9.scale("linear", 0, 1, frame=7)
580  if displayPsfSpatialModel:
581  maUtils.plotPsfSpatialModel(exposure, psf, psfCellSet, showBadCandidates=True,
582  matchKernelAmplitudes=matchKernelAmplitudes,
583  keepPlots=keepMatplotlibPlots)
584  #
585  # Generate some QA information
586  #
587  # Count PSF stars
588  #
589  numGoodStars = 0
590  numAvailStars = 0
591 
592  avgX = 0.0
593  avgY = 0.0
594 
595  for cell in psfCellSet.getCellList():
596  for cand in cell.begin(False): # don't ignore BAD stars
597  numAvailStars += 1
598 
599  for cand in cell.begin(True): # do ignore BAD stars
600  cand = algorithmsLib.cast_PsfCandidateF(cand)
601  src = cand.getSource()
602  if flagKey is not None:
603  src.set(flagKey, True)
604  avgX += src.getX()
605  avgY += src.getY()
606  numGoodStars += 1
607 
608  avgX /= numGoodStars
609  avgY /= numGoodStars
610 
611  if metadata != None:
612  metadata.set("spatialFitChi2", fitChi2)
613  metadata.set("numGoodStars", numGoodStars)
614  metadata.set("numAvailStars", numAvailStars)
615  metadata.set("avgX", avgX)
616  metadata.set("avgY", avgY)
617 
618  psf = algorithmsLib.PcaPsf(psf.getKernel(), afwGeom.Point2D(avgX, avgY))
619 
620  return psf, psfCellSet
621 
622 
623 def candidatesIter(psfCellSet, ignoreBad=True):
624  """!Generator for Psf candidates
625 
626  This allows two 'for' loops to be reduced to one.
627 
628  \param psfCellSet SpatialCellSet of PSF candidates
629  \param ignoreBad Ignore candidates flagged as BAD?
630  \return SpatialCell, PsfCandidate
631  """
632  for cell in psfCellSet.getCellList():
633  for cand in cell.begin(ignoreBad):
634  yield (cell, algorithmsLib.cast_PsfCandidateF(cand))
STL namespace.
a place to record messages and descriptions of the state of processing.
Definition: Log.h:154
boost::enable_if< typename ExpressionTraits< Scalar >::IsScalar, Scalar >::type sum(Scalar const &scalar)
Definition: operators.h:1250
A collection of SpatialCells covering an entire image.
Definition: SpatialCell.h:378
def determinePsf
Determine a PCA PSF model for an exposure given a list of PSF candidates.
Statistics makeStatistics(afwImage::Mask< afwImage::MaskPixel > const &msk, int const flags, StatisticsControl const &sctrl)
Specialization to handle Masks.
Definition: Statistics.cc:1082
def candidatesIter
Generator for Psf candidates.