LSST Applications  21.0.0-147-g0e635eb1+1acddb5be5,22.0.0+052faf71bd,22.0.0+1ea9a8b2b2,22.0.0+6312710a6c,22.0.0+729191ecac,22.0.0+7589c3a021,22.0.0+9f079a9461,22.0.1-1-g7d6de66+b8044ec9de,22.0.1-1-g87000a6+536b1ee016,22.0.1-1-g8e32f31+6312710a6c,22.0.1-10-gd060f87+016f7cdc03,22.0.1-12-g9c3108e+df145f6f68,22.0.1-16-g314fa6d+c825727ab8,22.0.1-19-g93a5c75+d23f2fb6d8,22.0.1-19-gb93eaa13+aab3ef7709,22.0.1-2-g8ef0a89+b8044ec9de,22.0.1-2-g92698f7+9f079a9461,22.0.1-2-ga9b0f51+052faf71bd,22.0.1-2-gac51dbf+052faf71bd,22.0.1-2-gb66926d+6312710a6c,22.0.1-2-gcb770ba+09e3807989,22.0.1-20-g32debb5+b8044ec9de,22.0.1-23-gc2439a9a+fb0756638e,22.0.1-3-g496fd5d+09117f784f,22.0.1-3-g59f966b+1e6ba2c031,22.0.1-3-g849a1b8+f8b568069f,22.0.1-3-gaaec9c0+c5c846a8b1,22.0.1-32-g5ddfab5d3+60ce4897b0,22.0.1-4-g037fbe1+64e601228d,22.0.1-4-g8623105+b8044ec9de,22.0.1-5-g096abc9+d18c45d440,22.0.1-5-g15c806e+57f5c03693,22.0.1-7-gba73697+57f5c03693,master-g6e05de7fdc+c1283a92b8,master-g72cdda8301+729191ecac,w.2021.39
LSST Data Management Base Package
psfexPsfDeterminer.py
Go to the documentation of this file.
1 #
2 # LSST Data Management System
3 # Copyright 2008-2015 AURA/LSST.
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 <https://www.lsstcorp.org/LegalNotices/>.
21 #
22 import os
23 import numpy as np
24 
25 import lsst.daf.base as dafBase
26 import lsst.pex.config as pexConfig
27 import lsst.pex.exceptions as pexExcept
28 import lsst.geom as geom
29 import lsst.afw.geom.ellipses as afwEll
30 import lsst.afw.display as afwDisplay
31 import lsst.afw.image as afwImage
32 import lsst.afw.math as afwMath
33 import lsst.meas.algorithms as measAlg
34 import lsst.meas.algorithms.utils as maUtils
35 import lsst.meas.extensions.psfex as psfex
36 
37 
38 class PsfexPsfDeterminerConfig(measAlg.BasePsfDeterminerConfig):
39  spatialOrder = pexConfig.Field(
40  doc="specify spatial order for PSF kernel creation",
41  dtype=int,
42  default=2,
43  check=lambda x: x >= 1,
44  )
45  sizeCellX = pexConfig.Field(
46  doc="size of cell used to determine PSF (pixels, column direction)",
47  dtype=int,
48  default=256,
49  # minValue = 10,
50  check=lambda x: x >= 10,
51  )
52  sizeCellY = pexConfig.Field(
53  doc="size of cell used to determine PSF (pixels, row direction)",
54  dtype=int,
55  default=sizeCellX.default,
56  # minValue = 10,
57  check=lambda x: x >= 10,
58  )
59  samplingSize = pexConfig.Field(
60  doc="Resolution of the internal PSF model relative to the pixel size; "
61  "e.g. 0.5 is equal to 2x oversampling",
62  dtype=float,
63  default=0.5,
64  )
65  badMaskBits = pexConfig.ListField(
66  doc="List of mask bits which cause a source to be rejected as bad "
67  "N.b. INTRP is used specially in PsfCandidateSet; it means \"Contaminated by neighbour\"",
68  dtype=str,
69  default=["INTRP", "SAT"],
70  )
71  psfexBasis = pexConfig.ChoiceField(
72  doc="BASIS value given to psfex. PIXEL_AUTO will use the requested samplingSize only if "
73  "the FWHM < 3 pixels. Otherwise, it will use samplingSize=1. PIXEL will always use the "
74  "requested samplingSize",
75  dtype=str,
76  allowed={
77  "PIXEL": "Always use requested samplingSize",
78  "PIXEL_AUTO": "Only use requested samplingSize when FWHM < 3",
79  },
80  default='PIXEL_AUTO',
81  optional=False,
82  )
83  tolerance = pexConfig.Field(
84  doc="tolerance of spatial fitting",
85  dtype=float,
86  default=1e-2,
87  )
88  lam = pexConfig.Field(
89  doc="floor for variance is lam*data",
90  dtype=float,
91  default=0.05,
92  )
93  reducedChi2ForPsfCandidates = pexConfig.Field(
94  doc="for psf candidate evaluation",
95  dtype=float,
96  default=2.0,
97  )
98  spatialReject = pexConfig.Field(
99  doc="Rejection threshold (stdev) for candidates based on spatial fit",
100  dtype=float,
101  default=3.0,
102  )
103  recentroid = pexConfig.Field(
104  doc="Should PSFEX be permitted to recentroid PSF candidates?",
105  dtype=bool,
106  default=False,
107  )
108  kernelSize = pexConfig.Field(
109  doc=("Size of the postage stamp around each star that is extracted for fitting."
110  "Note: this reflects the oversampling setting of the psf, set by `samplingSize`;"
111  "e.g. `samplingSize=0.5` would require this value to be 2x what you expect."),
112  dtype=int,
113  default=81,
114  )
115 
116 
117 class PsfexPsfDeterminerTask(measAlg.BasePsfDeterminerTask):
118  ConfigClass = PsfexPsfDeterminerConfig
119 
120  def determinePsf(self, exposure, psfCandidateList, metadata=None, flagKey=None):
121  """Determine a PSFEX PSF model for an exposure given a list of PSF
122  candidates.
123 
124  Parameters
125  ----------
126  exposure: `lsst.afw.image.Exposure`
127  Exposure containing the PSF candidates.
128  psfCandidateList: iterable of `lsst.meas.algorithms.PsfCandidate`
129  Sequence of PSF candidates typically obtained by detecting sources
130  and then running them through a star selector.
131  metadata: metadata, optional
132  A home for interesting tidbits of information.
133  flagKey: `lsst.afw.table.Key`, optional
134  Schema key used to mark sources actually used in PSF determination.
135 
136  Returns
137  -------
138  psf: `lsst.meas.extensions.psfex.PsfexPsf`
139  The determined PSF.
140  """
141 
142  import lsstDebug
143  display = lsstDebug.Info(__name__).display
144  displayExposure = display and \
145  lsstDebug.Info(__name__).displayExposure # display the Exposure + spatialCells
146  displayPsfComponents = display and \
147  lsstDebug.Info(__name__).displayPsfComponents # show the basis functions
148  showBadCandidates = display and \
149  lsstDebug.Info(__name__).showBadCandidates # Include bad candidates (meaningless, methinks)
150  displayResiduals = display and \
151  lsstDebug.Info(__name__).displayResiduals # show residuals
152  displayPsfMosaic = display and \
153  lsstDebug.Info(__name__).displayPsfMosaic # show mosaic of reconstructed PSF(x,y)
154  normalizeResiduals = lsstDebug.Info(__name__).normalizeResiduals
155  afwDisplay.setDefaultMaskTransparency(75)
156  # Normalise residuals by object amplitude
157 
158  mi = exposure.getMaskedImage()
159 
160  nCand = len(psfCandidateList)
161  if nCand == 0:
162  raise RuntimeError("No PSF candidates supplied.")
163  #
164  # How big should our PSF models be?
165  #
166  if display: # only needed for debug plots
167  # construct and populate a spatial cell set
168  bbox = mi.getBBox(afwImage.PARENT)
169  psfCellSet = afwMath.SpatialCellSet(bbox, self.config.sizeCellX, self.config.sizeCellY)
170  else:
171  psfCellSet = None
172 
173  sizes = np.empty(nCand)
174  for i, psfCandidate in enumerate(psfCandidateList):
175  try:
176  if psfCellSet:
177  psfCellSet.insertCandidate(psfCandidate)
178  except Exception as e:
179  self.log.error("Skipping PSF candidate %d of %d: %s", i, len(psfCandidateList), e)
180  continue
181 
182  source = psfCandidate.getSource()
183  quad = afwEll.Quadrupole(source.getIxx(), source.getIyy(), source.getIxy())
184  rmsSize = quad.getTraceRadius()
185  sizes[i] = rmsSize
186 
187  if self.config.kernelSize >= 15:
188  self.log.warn("NOT scaling kernelSize by stellar quadrupole moment, but using absolute value")
189  actualKernelSize = self.config.kernelSize
190  else:
191  actualKernelSize = 2 * int(self.config.kernelSize * np.sqrt(np.median(sizes)) + 0.5) + 1
192  if actualKernelSize < self.config.kernelSizeMin:
193  actualKernelSize = self.config.kernelSizeMin
194  if actualKernelSize > self.config.kernelSizeMax:
195  actualKernelSize = self.config.kernelSizeMax
196  if display:
197  rms = np.median(sizes)
198  msg = "Median PSF RMS size=%.2f pixels (\"FWHM\"=%.2f)" % (rms, 2*np.sqrt(2*np.log(2))*rms)
199  self.log.debug(msg)
200 
201  # If we manually set the resolution then we need the size in pixel
202  # units
203  pixKernelSize = actualKernelSize
204  if self.config.samplingSize > 0:
205  pixKernelSize = int(actualKernelSize*self.config.samplingSize)
206  if pixKernelSize % 2 == 0:
207  pixKernelSize += 1
208  self.log.trace("Psfex Kernel size=%.2f, Image Kernel Size=%.2f", actualKernelSize, pixKernelSize)
209  psfCandidateList[0].setHeight(pixKernelSize)
210  psfCandidateList[0].setWidth(pixKernelSize)
211 
212  # -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- BEGIN PSFEX
213  #
214  # Insert the good candidates into the set
215  #
216  defaultsFile = os.path.join(os.environ["MEAS_EXTENSIONS_PSFEX_DIR"], "config", "default-lsst.psfex")
217  args_md = dafBase.PropertySet()
218  args_md.set("BASIS_TYPE", str(self.config.psfexBasis))
219  args_md.set("PSFVAR_DEGREES", str(self.config.spatialOrder))
220  args_md.set("PSF_SIZE", str(actualKernelSize))
221  args_md.set("PSF_SAMPLING", str(self.config.samplingSize))
222  prefs = psfex.Prefs(defaultsFile, args_md)
223  prefs.setCommandLine([])
224  prefs.addCatalog("psfexPsfDeterminer")
225 
226  prefs.use()
227  principalComponentExclusionFlag = bool(bool(psfex.Context.REMOVEHIDDEN)
228  if False else psfex.Context.KEEPHIDDEN)
229  context = psfex.Context(prefs.getContextName(), prefs.getContextGroup(),
230  prefs.getGroupDeg(), principalComponentExclusionFlag)
231  psfSet = psfex.Set(context)
232  psfSet.setVigSize(pixKernelSize, pixKernelSize)
233  psfSet.setFwhm(2*np.sqrt(2*np.log(2))*np.median(sizes))
234  psfSet.setRecentroid(self.config.recentroid)
235 
236  catindex, ext = 0, 0
237  backnoise2 = afwMath.makeStatistics(mi.getImage(), afwMath.VARIANCECLIP).getValue()
238  ccd = exposure.getDetector()
239  if ccd:
240  gain = np.mean(np.array([a.getGain() for a in ccd]))
241  else:
242  gain = 1.0
243  self.log.warn("Setting gain to %g" % (gain,))
244 
245  contextvalp = []
246  for i, key in enumerate(context.getName()):
247  if key[0] == ':':
248  try:
249  contextvalp.append(exposure.getMetadata().getScalar(key[1:]))
250  except KeyError as e:
251  raise RuntimeError("%s parameter not found in the header of %s" %
252  (key[1:], prefs.getContextName())) from e
253  else:
254  try:
255  contextvalp.append(np.array([psfCandidateList[_].getSource().get(key)
256  for _ in range(nCand)]))
257  except KeyError as e:
258  raise RuntimeError("%s parameter not found" % (key,)) from e
259  psfSet.setContextname(i, key)
260 
261  if display:
262  frame = 0
263  if displayExposure:
264  disp = afwDisplay.Display(frame=frame)
265  disp.mtv(exposure, title="psf determination")
266 
267  badBits = mi.getMask().getPlaneBitMask(self.config.badMaskBits)
268  fluxName = prefs.getPhotfluxRkey()
269  fluxFlagName = "base_" + fluxName + "_flag"
270 
271  xpos, ypos = [], []
272  for i, psfCandidate in enumerate(psfCandidateList):
273  source = psfCandidate.getSource()
274 
275  # skip sources with bad centroids
276  xc, yc = source.getX(), source.getY()
277  if not np.isfinite(xc) or not np.isfinite(yc):
278  continue
279  # skip flagged sources
280  if fluxFlagName in source.schema and source.get(fluxFlagName):
281  continue
282  # skip nonfinite and negative sources
283  flux = source.get(fluxName)
284  if flux < 0 or not np.isfinite(flux):
285  continue
286 
287  try:
288  pstamp = psfCandidate.getMaskedImage().clone()
289  except pexExcept.LengthError:
290  # Candidate is too close to the edge to get a stamp. Skip.
291  # TODO DM-27547: Replace with geometric condition
292  continue
293 
294  # From this point, we're configuring the "sample" (PSFEx's version
295  # of a PSF candidate).
296  # Having created the sample, we must proceed to configure it, and
297  # then fini (finalize), or it will be malformed.
298  try:
299  sample = psfSet.newSample()
300  sample.setCatindex(catindex)
301  sample.setExtindex(ext)
302  sample.setObjindex(i)
303 
304  imArray = pstamp.getImage().getArray()
305  imArray[np.where(np.bitwise_and(pstamp.getMask().getArray(), badBits))] = \
306  -2*psfex.BIG
307  sample.setVig(imArray)
308 
309  sample.setNorm(flux)
310  sample.setBacknoise2(backnoise2)
311  sample.setGain(gain)
312  sample.setX(xc)
313  sample.setY(yc)
314  sample.setFluxrad(sizes[i])
315 
316  for j in range(psfSet.getNcontext()):
317  sample.setContext(j, float(contextvalp[j][i]))
318  except Exception as e:
319  self.log.error("Exception when processing sample at (%f,%f): %s", xc, yc, e)
320  continue
321  else:
322  psfSet.finiSample(sample)
323 
324  xpos.append(xc) # for QA
325  ypos.append(yc)
326 
327  if displayExposure:
328  with disp.Buffering():
329  disp.dot("o", xc, yc, ctype=afwDisplay.CYAN, size=4)
330 
331  if psfSet.getNsample() == 0:
332  raise RuntimeError("No good PSF candidates to pass to PSFEx")
333 
334  # ---- Update min and max and then the scaling
335  for i in range(psfSet.getNcontext()):
336  cmin = contextvalp[i].min()
337  cmax = contextvalp[i].max()
338  psfSet.setContextScale(i, cmax - cmin)
339  psfSet.setContextOffset(i, (cmin + cmax)/2.0)
340 
341  # Don't waste memory!
342  psfSet.trimMemory()
343 
344  # =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- END PSFEX
345  #
346  # Do a PSFEX decomposition of those PSF candidates
347  #
348  fields = []
349  field = psfex.Field("Unknown")
350  field.addExt(exposure.getWcs(), exposure.getWidth(), exposure.getHeight(), psfSet.getNsample())
351  field.finalize()
352 
353  fields.append(field)
354 
355  sets = []
356  sets.append(psfSet)
357 
358  psfex.makeit(fields, sets)
359  psfs = field.getPsfs()
360 
361  # Flag which objects were actually used in psfex by
362  good_indices = []
363  for i in range(sets[0].getNsample()):
364  index = sets[0].getSample(i).getObjindex()
365  if index > -1:
366  good_indices.append(index)
367 
368  if flagKey is not None:
369  for i, psfCandidate in enumerate(psfCandidateList):
370  source = psfCandidate.getSource()
371  if i in good_indices:
372  source.set(flagKey, True)
373 
374  xpos = np.array(xpos)
375  ypos = np.array(ypos)
376  numGoodStars = len(good_indices)
377  avgX, avgY = np.mean(xpos), np.mean(ypos)
378 
379  psf = psfex.PsfexPsf(psfs[0], geom.Point2D(avgX, avgY))
380 
381  #
382  # Display code for debugging
383  #
384  if display:
385  assert psfCellSet is not None
386 
387  if displayExposure:
388  maUtils.showPsfSpatialCells(exposure, psfCellSet, showChi2=True,
389  symb="o", ctype=afwDisplay.YELLOW, ctypeBad=afwDisplay.RED,
390  size=8, display=disp)
391  if displayResiduals:
392  disp4 = afwDisplay.Display(frame=4)
393  maUtils.showPsfCandidates(exposure, psfCellSet, psf=psf, display=disp4,
394  normalize=normalizeResiduals,
395  showBadCandidates=showBadCandidates)
396  if displayPsfComponents:
397  disp6 = afwDisplay.Display(frame=6)
398  maUtils.showPsf(psf, display=disp6)
399  if displayPsfMosaic:
400  disp7 = afwDisplay.Display(frame=7)
401  maUtils.showPsfMosaic(exposure, psf, display=disp7, showFwhm=True)
402  disp.scale('linear', 0, 1)
403  #
404  # Generate some QA information
405  #
406  # Count PSF stars
407  #
408  if metadata is not None:
409  metadata.set("spatialFitChi2", np.nan)
410  metadata.set("numAvailStars", nCand)
411  metadata.set("numGoodStars", numGoodStars)
412  metadata.set("avgX", avgX)
413  metadata.set("avgY", avgY)
414 
415  return psf, psfCellSet
416 
417 
418 measAlg.psfDeterminerRegistry.register("psfex", PsfexPsfDeterminerTask)
int min
int max
A collection of SpatialCells covering an entire image.
Definition: SpatialCell.h:383
Class for storing generic metadata.
Definition: PropertySet.h:66
Represent a PSF as a linear combination of PSFEX (== Karhunen-Loeve) basis functions.
Definition: PsfexPsf.h:40
def determinePsf(self, exposure, psfCandidateList, metadata=None, flagKey=None)
Reports attempts to exceed implementation-defined length limits for some classes.
Definition: Runtime.h:76
Backwards-compatibility support for depersisting the old Calib (FluxMag0/FluxMag0Err) objects.
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:359
Fit spatial kernel using approximate fluxes for candidates, and solving a linear system of equations.