Loading [MathJax]/extensions/tex2jax.js
LSST Applications g04a91732dc+cc8870d3f5,g07dc498a13+5aa0b8792f,g0fba68d861+80045be308,g1409bbee79+5aa0b8792f,g1a7e361dbc+5aa0b8792f,g1fd858c14a+f64bc332a9,g208c678f98+1ae86710ed,g35bb328faa+fcb1d3bbc8,g4d2262a081+47ad8a29a8,g4d39ba7253+9633a327c1,g4e0f332c67+5d362be553,g53246c7159+fcb1d3bbc8,g60b5630c4e+9633a327c1,g668ecb457e+25d63fd678,g78460c75b0+2f9a1b4bcd,g786e29fd12+cf7ec2a62a,g7b71ed6315+fcb1d3bbc8,g8852436030+8b64ca622a,g89139ef638+5aa0b8792f,g89e1512fd8+04325574d3,g8d6b6b353c+9633a327c1,g9125e01d80+fcb1d3bbc8,g989de1cb63+5aa0b8792f,g9f33ca652e+b196626af7,ga9baa6287d+9633a327c1,gaaedd4e678+5aa0b8792f,gabe3b4be73+1e0a283bba,gb1101e3267+71e32094df,gb58c049af0+f03b321e39,gb90eeb9370+2807b1ad02,gcf25f946ba+8b64ca622a,gd315a588df+a39986a76f,gd6cbbdb0b4+c8606af20c,gd9a9a58781+fcb1d3bbc8,gde0f65d7ad+4e42d81ab7,ge278dab8ac+932305ba37,ge82c20c137+76d20ab76d,gfe73954cf8+a1301e4c20,w.2025.11
LSST Data Management Base Package
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Modules Pages
psfexPsfDeterminer.py
Go to the documentation of this file.
1# This file is part of meas_extensions_psfex.
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__ = (
23 "PsfexPsfDeterminerConfig",
24 "PsfexPsfDeterminerTask",
25 "PsfexNoStarsError",
26 "PsfexNoGoodStarsError",
27 "PsfexTooFewGoodStarsError",
28)
29
30import os
31import numpy as np
32
33import lsst.daf.base as dafBase
34import lsst.pex.config as pexConfig
35import lsst.pex.exceptions as pexExcept
36import lsst.geom as geom
37import lsst.afw.geom.ellipses as afwEll
38import lsst.afw.display as afwDisplay
39import lsst.afw.image as afwImage
40import lsst.afw.math as afwMath
41import lsst.meas.algorithms as measAlg
42import lsst.meas.algorithms.utils as maUtils
43import lsst.meas.extensions.psfex as psfex
44from lsst.pipe.base import AlgorithmError
45
46
47class PsfexNoStarsError(AlgorithmError):
48 """Raised if no stars are available for PSF determination."""
49
50 def __init__(self) -> None:
51 super().__init__("No psf candidates supplied.")
52
53
54class PsfexNoGoodStarsError(AlgorithmError):
55 """Raised if no "good" stars are available for PSF determination.
56
57 Parameters
58 ----------
59 num_available_stars : `int`
60 Number of available stars for PSF determination.
61 """
62
63 def __init__(self, num_available_stars) -> None:
64 self._num_available_stars = num_available_stars
65 super().__init__(f"No good psf candidates to pass to psfex out of {num_available_stars} available.")
66
67 @property
68 def metadata(self) -> dict:
69 return {
70 "num_available_stars": self._num_available_stars,
71 }
72
73
74class PsfexTooFewGoodStarsError(AlgorithmError):
75 """Raised if too few good stars are available for PSF determination.
76
77 Parameters
78 ----------
79 num_available_stars : `int`
80 Number of available stars for PSF determination.
81 num_good_stars : `int`
82 Number of good stars used for PSF determination.
83 poly_ndim_initial : `int`
84 Initial number of dependency parameters (dimensions) used in
85 polynomial fitting.
86 poly_ndim_final : `int`
87 Final number of dependency parameters (dimensions) set in the PSFEx
88 model after initializing the PSF structure.
89 """
90
92 self,
93 num_available_stars,
94 num_good_stars,
95 poly_ndim_initial,
96 poly_ndim_final,
97 ) -> None:
98 self._num_available_stars = num_available_stars
99 self._num_good_stars = num_good_stars
100 self._poly_ndim_initial = poly_ndim_initial
101 self._poly_ndim_final = poly_ndim_final
102 super().__init__(
103 f"Failed to determine psfex psf: too few good stars ({num_good_stars}) out of "
104 f"{num_available_stars} available. To accommodate this low count, the polynomial dimension was "
105 f"lowered from {poly_ndim_initial} to {poly_ndim_final}, which is not handled by the Science "
106 "Pipelines code."
107 )
108
109 @property
110 def metadata(self) -> dict:
111 return {
112 "num_available_stars": self._num_available_stars,
113 "num_good_stars": self._num_good_stars,
114 "poly_ndim_initial": self._poly_ndim_initial,
115 "poly_ndim_final": self._poly_ndim_final,
116 }
117
118
119class PsfexPsfDeterminerConfig(measAlg.BasePsfDeterminerConfig):
120 spatialOrder = pexConfig.Field[int](
121 doc="specify spatial order for PSF kernel creation",
122 default=2,
123 check=lambda x: x >= 1,
124 )
125 sizeCellX = pexConfig.Field[int](
126 doc="size of cell used to determine PSF (pixels, column direction)",
127 default=256,
128 # minValue = 10,
129 check=lambda x: x >= 10,
130 )
131 sizeCellY = pexConfig.Field[int](
132 doc="size of cell used to determine PSF (pixels, row direction)",
133 default=sizeCellX.default,
134 # minValue = 10,
135 check=lambda x: x >= 10,
136 )
137 samplingSize = pexConfig.Field[float](
138 doc="Resolution of the internal PSF model relative to the pixel size; "
139 "e.g. 0.5 is equal to 2x oversampling",
140 default=0.5,
141 )
142 badMaskBits = pexConfig.ListField[str](
143 doc="List of mask bits which cause a source to be rejected as bad "
144 "N.b. INTRP is used specially in PsfCandidateSet; it means \"Contaminated by neighbour\"",
145 default=["INTRP", "SAT"],
146 )
147 psfexBasis = pexConfig.ChoiceField[str](
148 doc="BASIS value given to psfex. PIXEL_AUTO will use the requested samplingSize only if "
149 "the FWHM < 3 pixels. Otherwise, it will use samplingSize=1. PIXEL will always use the "
150 "requested samplingSize",
151 allowed={
152 "PIXEL": "Always use requested samplingSize",
153 "PIXEL_AUTO": "Only use requested samplingSize when FWHM < 3",
154 },
155 default='PIXEL_AUTO',
156 optional=False,
157 )
158 tolerance = pexConfig.Field[float](
159 doc="tolerance of spatial fitting",
160 default=1e-2,
161 )
162 lam = pexConfig.Field[float](
163 doc="floor for variance is lam*data",
164 default=0.05,
165 )
166 reducedChi2ForPsfCandidates = pexConfig.Field[float](
167 doc="for psf candidate evaluation",
168 default=2.0,
169 )
170 spatialReject = pexConfig.Field[float](
171 doc="Rejection threshold (stdev) for candidates based on spatial fit",
172 default=3.0,
173 )
174 recentroid = pexConfig.Field[bool](
175 doc="Should PSFEX be permitted to recentroid PSF candidates?",
176 default=False,
177 )
178 photometricFluxField = pexConfig.Field[str](
179 doc="Flux field to use for photometric normalization. This overrides the "
180 "``PHOTFLUX_KEY`` field for psfex. The associated flux error is "
181 "derived by appending ``Err`` to this field.",
182 default="base_CircularApertureFlux_9_0_instFlux",
183 )
184
185 def setDefaults(self):
186 super().setDefaults()
187 self.stampSize = 41
188
189 def validate(self):
190 super().validate()
191 fluxFieldEnding = self.photometricFluxField.rpartition("_")[-1]
192 if "_" not in self.photometricFluxField or not (
193 fluxFieldEnding == "flux" or fluxFieldEnding.endswith("Flux")
194 ):
195 raise ValueError(
196 "Expected `photometricFluxField` to end with '_flux' "
197 f"or '_*Flux', but got '{self.photometricFluxField}'"
198 )
199
200
201class PsfexPsfDeterminerTask(measAlg.BasePsfDeterminerTask):
202 ConfigClass = PsfexPsfDeterminerConfig
203
204 def determinePsf(self, exposure, psfCandidateList, metadata=None, flagKey=None):
205 """Determine a PSFEX PSF model for an exposure given a list of PSF
206 candidates.
207
208 Parameters
209 ----------
210 exposure: `lsst.afw.image.Exposure`
211 Exposure containing the PSF candidates.
212 psfCandidateList: iterable of `lsst.meas.algorithms.PsfCandidate`
213 Sequence of PSF candidates typically obtained by detecting sources
214 and then running them through a star selector.
215 metadata: metadata, optional
216 A home for interesting tidbits of information.
217 flagKey: `lsst.afw.table.Key`, optional
218 Schema key used to mark sources actually used in PSF determination.
219
220 Returns
221 -------
222 psf: `lsst.meas.extensions.psfex.PsfexPsf`
223 The determined PSF.
224 """
225 psfCandidateList = self.downsampleCandidates(psfCandidateList)
226
227 import lsstDebug
228 display = lsstDebug.Info(__name__).display
229 displayExposure = display and \
230 lsstDebug.Info(__name__).displayExposure # display the Exposure + spatialCells
231 displayPsfComponents = display and \
232 lsstDebug.Info(__name__).displayPsfComponents # show the basis functions
233 showBadCandidates = display and \
234 lsstDebug.Info(__name__).showBadCandidates # Include bad candidates (meaningless, methinks)
235 displayResiduals = display and \
236 lsstDebug.Info(__name__).displayResiduals # show residuals
237 displayPsfMosaic = display and \
238 lsstDebug.Info(__name__).displayPsfMosaic # show mosaic of reconstructed PSF(x,y)
239 normalizeResiduals = lsstDebug.Info(__name__).normalizeResiduals
240 afwDisplay.setDefaultMaskTransparency(75)
241 # Normalise residuals by object amplitude
242
243 mi = exposure.getMaskedImage()
244
245 nCand = len(psfCandidateList)
246 if nCand == 0:
247 raise PsfexNoStarsError()
248 #
249 # How big should our PSF models be?
250 #
251 if display: # only needed for debug plots
252 # construct and populate a spatial cell set
253 bbox = mi.getBBox(afwImage.PARENT)
254 psfCellSet = afwMath.SpatialCellSet(bbox, self.config.sizeCellX, self.config.sizeCellY)
255 else:
256 psfCellSet = None
257
258 sizes = np.empty(nCand)
259 for i, psfCandidate in enumerate(psfCandidateList):
260 try:
261 if psfCellSet:
262 psfCellSet.insertCandidate(psfCandidate)
263 except Exception as e:
264 self.log.error("Skipping PSF candidate %d of %d: %s", i, len(psfCandidateList), e)
265 continue
266
267 source = psfCandidate.getSource()
268 quad = afwEll.Quadrupole(source.getIxx(), source.getIyy(), source.getIxy())
269 rmsSize = quad.getTraceRadius()
270 sizes[i] = rmsSize
271
272 pixKernelSize = self.config.stampSize
273 actualKernelSize = int(2*np.floor(0.5*pixKernelSize/self.config.samplingSize) + 1)
274
275 if display:
276 rms = np.median(sizes)
277 self.log.debug("Median PSF RMS size=%.2f pixels (\"FWHM\"=%.2f)",
278 rms, 2*np.sqrt(2*np.log(2))*rms)
279
280 self.log.trace("Psfex Kernel size=%.2f, Image Kernel Size=%.2f", actualKernelSize, pixKernelSize)
281
282 # -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- BEGIN PSFEX
283 #
284 # Insert the good candidates into the set
285 #
286 defaultsFile = os.path.join(os.environ["MEAS_EXTENSIONS_PSFEX_DIR"], "config", "default-lsst.psfex")
287 args_md = dafBase.PropertySet()
288 args_md.set("BASIS_TYPE", str(self.config.psfexBasis))
289 args_md.set("PSFVAR_DEGREES", str(self.config.spatialOrder))
290 args_md.set("PSF_SIZE", str(actualKernelSize))
291 args_md.set("PSF_SAMPLING", str(self.config.samplingSize))
292 args_md.set("PHOTFLUX_KEY", str(self.config.photometricFluxField))
293 args_md.set("PHOTFLUXERR_KEY", str(self.config.photometricFluxField) + "Err")
294 prefs = psfex.Prefs(defaultsFile, args_md)
295 prefs.setCommandLine([])
296 prefs.addCatalog("psfexPsfDeterminer")
297
298 prefs.use()
299 principalComponentExclusionFlag = bool(bool(psfex.Context.REMOVEHIDDEN)
300 if False else psfex.Context.KEEPHIDDEN)
301 context = psfex.Context(prefs.getContextName(), prefs.getContextGroup(),
302 prefs.getGroupDeg(), principalComponentExclusionFlag)
303 psfSet = psfex.Set(context)
304 psfSet.setVigSize(pixKernelSize, pixKernelSize)
305 psfSet.setFwhm(2*np.sqrt(2*np.log(2))*np.median(sizes))
306 psfSet.setRecentroid(self.config.recentroid)
307
308 catindex, ext = 0, 0
309 backnoise2 = afwMath.makeStatistics(mi.getImage(), afwMath.VARIANCECLIP).getValue()
310 ccd = exposure.getDetector()
311 if ccd:
312 gain = np.mean(np.array([a.getGain() for a in ccd]))
313 else:
314 gain = 1.0
315 self.log.warning("Setting gain to %g", gain)
316
317 contextvalp = []
318 for i, key in enumerate(context.getName()):
319 if key[0] == ':':
320 try:
321 contextvalp.append(exposure.getMetadata().getScalar(key[1:]))
322 except KeyError as e:
323 raise RuntimeError("%s parameter not found in the header of %s" %
324 (key[1:], prefs.getContextName())) from e
325 else:
326 try:
327 contextvalp.append(np.array([psfCandidateList[_].getSource().get(key)
328 for _ in range(nCand)]))
329 except KeyError as e:
330 raise RuntimeError("%s parameter not found" % (key,)) from e
331 psfSet.setContextname(i, key)
332
333 if display:
334 frame = 0
335 if displayExposure:
336 disp = afwDisplay.Display(frame=frame)
337 disp.mtv(exposure, title="psf determination")
338
339 badBits = mi.getMask().getPlaneBitMask(self.config.badMaskBits)
340 fluxName = prefs.getPhotfluxRkey()
341 fluxFlagName = fluxName.rpartition("_")[0] + "_flag"
342
343 xpos, ypos = [], []
344 for i, psfCandidate in enumerate(psfCandidateList):
345 source = psfCandidate.getSource()
346
347 # skip sources with bad centroids
348 xc, yc = source.getX(), source.getY()
349 if not np.isfinite(xc) or not np.isfinite(yc):
350 continue
351 # skip flagged sources
352 if source.get(fluxFlagName):
353 continue
354 # skip nonfinite and negative sources
355 flux = source.get(fluxName)
356 if flux < 0 or not np.isfinite(flux):
357 continue
358
359 try:
360 pstamp = psfCandidate.getMaskedImage(pixKernelSize, pixKernelSize).clone()
361 except pexExcept.LengthError:
362 self.log.warning("Could not get stamp image for psfCandidate: %s with kernel size: %s",
363 psfCandidate, pixKernelSize)
364 continue
365
366 # From this point, we're configuring the "sample" (PSFEx's version
367 # of a PSF candidate).
368 # Having created the sample, we must proceed to configure it, and
369 # then fini (finalize), or it will be malformed.
370 try:
371 sample = psfSet.newSample()
372 sample.setCatindex(catindex)
373 sample.setExtindex(ext)
374 sample.setObjindex(i)
375
376 imArray = pstamp.getImage().getArray()
377 imArray[np.where(np.bitwise_and(pstamp.getMask().getArray(), badBits))] = \
378 -2*psfex.BIG
379 sample.setVig(imArray)
380
381 sample.setNorm(flux)
382 sample.setBacknoise2(backnoise2)
383 sample.setGain(gain)
384 sample.setX(xc)
385 sample.setY(yc)
386 sample.setFluxrad(sizes[i])
387
388 for j in range(psfSet.getNcontext()):
389 sample.setContext(j, float(contextvalp[j][i]))
390 except Exception as e:
391 self.log.error("Exception when processing sample at (%f,%f): %s", xc, yc, e)
392 continue
393 else:
394 psfSet.finiSample(sample)
395
396 xpos.append(xc) # for QA
397 ypos.append(yc)
398
399 if displayExposure:
400 with disp.Buffering():
401 disp.dot("o", xc, yc, ctype=afwDisplay.CYAN, size=4)
402
403 if psfSet.getNsample() == 0:
404 raise PsfexNoGoodStarsError(num_available_stars=nCand)
405
406 # ---- Update min and max and then the scaling
407 for i in range(psfSet.getNcontext()):
408 cmin = contextvalp[i].min()
409 cmax = contextvalp[i].max()
410 psfSet.setContextScale(i, cmax - cmin)
411 psfSet.setContextOffset(i, (cmin + cmax)/2.0)
412
413 # Don't waste memory!
414 psfSet.trimMemory()
415
416 # =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- END PSFEX
417 #
418 # Do a PSFEX decomposition of those PSF candidates
419 #
420 fields = []
421 field = psfex.Field("Unknown")
422 field.addExt(exposure.getWcs(), exposure.getWidth(), exposure.getHeight(), psfSet.getNsample())
423 field.finalize()
424
425 fields.append(field)
426
427 sets = []
428 sets.append(psfSet)
429
430 psfex.makeit(fields, sets)
431 psfs = field.getPsfs()
432
433 # Flag which objects were actually used in psfex by
434 good_indices = []
435 for i in range(sets[0].getNsample()):
436 index = sets[0].getSample(i).getObjindex()
437 if index > -1:
438 good_indices.append(index)
439
440 if flagKey is not None:
441 for i, psfCandidate in enumerate(psfCandidateList):
442 source = psfCandidate.getSource()
443 if i in good_indices:
444 source.set(flagKey, True)
445
446 xpos = np.array(xpos)
447 ypos = np.array(ypos)
448 numGoodStars = len(good_indices)
449 avgX, avgY = np.mean(xpos), np.mean(ypos)
450
451 psf = psfex.PsfexPsf(psfs[0], geom.Point2D(avgX, avgY))
452
453 # If there are too few stars, the PSFEx psf model will reduce the order
454 # to 0, which the Science Pipelines code cannot handle (see
455 # https://github.com/lsst/meas_extensions_psfex/blob/f0d5218b5446faf5e39edc30e31d2e6f673ef294/src/PsfexPsf.cc#L118
456 # ).
457 if (ndim := psf.getNdim()) == 0:
459 num_available_stars=nCand,
460 num_good_stars=numGoodStars,
461 poly_ndim_initial=psfSet.getNcontext(),
462 poly_ndim_final=ndim,
463 )
464
465 #
466 # Display code for debugging
467 #
468 if display:
469 assert psfCellSet is not None
470
471 if displayExposure:
472 maUtils.showPsfSpatialCells(exposure, psfCellSet, showChi2=True,
473 symb="o", ctype=afwDisplay.YELLOW, ctypeBad=afwDisplay.RED,
474 size=8, display=disp)
475 if displayResiduals:
476 disp4 = afwDisplay.Display(frame=4)
477 maUtils.showPsfCandidates(exposure, psfCellSet, psf=psf, display=disp4,
478 normalize=normalizeResiduals,
479 showBadCandidates=showBadCandidates)
480 if displayPsfComponents:
481 disp6 = afwDisplay.Display(frame=6)
482 maUtils.showPsf(psf, display=disp6)
483 if displayPsfMosaic:
484 disp7 = afwDisplay.Display(frame=7)
485 maUtils.showPsfMosaic(exposure, psf, display=disp7, showFwhm=True)
486 disp.scale('linear', 0, 1)
487 #
488 # Generate some QA information
489 #
490 # Count PSF stars
491 #
492 if metadata is not None:
493 metadata["spatialFitChi2"] = np.nan
494 metadata["numAvailStars"] = nCand
495 metadata["numGoodStars"] = numGoodStars
496 metadata["avgX"] = avgX
497 metadata["avgY"] = avgY
498
499 return psf, psfCellSet
500
501
502measAlg.psfDeterminerRegistry.register("psfex", PsfexPsfDeterminerTask)
A collection of SpatialCells covering an entire image.
Class for storing generic metadata.
Definition PropertySet.h:67
Represent a PSF as a linear combination of PSFEX (== Karhunen-Loeve) basis functions.
Definition PsfexPsf.h:40
determinePsf(self, exposure, psfCandidateList, metadata=None, flagKey=None)
None __init__(self, num_available_stars, num_good_stars, poly_ndim_initial, poly_ndim_final)
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