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