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
piffPsfDeterminer.py
Go to the documentation of this file.
1# This file is part of meas_extensions_piff.
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__ = ["PiffPsfDeterminerConfig", "PiffPsfDeterminerTask"]
23
24import numpy as np
25import piff
26import galsim
27import re
28import logging
29import yaml
30
31import lsst.utils.logging
32from lsst.afw.cameraGeom import PIXELS, FIELD_ANGLE
33import lsst.pex.config as pexConfig
34import lsst.meas.algorithms as measAlg
35from lsst.meas.algorithms.psfDeterminer import BasePsfDeterminerTask
36from .piffPsf import PiffPsf
37from .wcs_wrapper import CelestialWcsWrapper, UVWcsWrapper
38
39
40def _validateGalsimInterpolant(name: str) -> bool:
41 """A helper function to validate the GalSim interpolant at config time.
42
43 Parameters
44 ----------
45 name : str
46 The name of the interpolant to use from GalSim. Valid options are:
47 galsim.Lanczos(N) or Lancsos(N), where N is a positive integer
48 galsim.Linear
49 galsim.Cubic
50 galsim.Quintic
51 galsim.Delta
52 galsim.Nearest
53 galsim.SincInterpolant
54
55 Returns
56 -------
57 is_valid : bool
58 Whether the provided interpolant name is valid.
59 """
60 # First, check if ``name`` is a valid Lanczos interpolant.
61 for pattern in (re.compile(r"Lanczos\‍(\d+\‍)"), re.compile(r"galsim.Lanczos\‍(\d+\‍)"),):
62 match = re.match(pattern, name) # Search from the start of the string.
63 if match is not None:
64 # Check that the pattern is also the end of the string.
65 return match.end() == len(name)
66
67 # If not, check if ``name`` is any other valid GalSim interpolant.
68 names = {f"galsim.{interp}" for interp in
69 ("Cubic", "Delta", "Linear", "Nearest", "Quintic", "SincInterpolant")
70 }
71 return name in names
72
73
74class PiffPsfDeterminerConfig(BasePsfDeterminerTask.ConfigClass):
75 spatialOrder = pexConfig.Field[int](
76 doc="Spatial order for PSF kernel creation. "
77 "Ignored if piffPsfConfigYaml is set.",
78 default=2,
79 )
80 piffBasisPolynomialSolver = pexConfig.ChoiceField[str](
81 doc="Solver to solve linear algebra for "
82 "the spatial interpolation of the PSF. "
83 "Ignore if piffPsfConfigYaml is not None.",
84 allowed=dict(
85 scipy="Default solver using scipy.",
86 cpp="C++ solver using eigen.",
87 jax="JAX solver.",
88 qr="QR decomposition using scipy/numpy.",
89 ),
90 default="scipy",
91 )
92 piffPixelGridFitCenter = pexConfig.Field[bool](
93 doc="PixelGrid can re-estimate center if this option is True. "
94 "Will use provided centroid if set to False. Default in Piff "
95 "is set to True. Depends on how input centroids can be trust. "
96 "Ignore if piffPsfConfigYaml is not None.",
97 default=True
98 )
99 samplingSize = pexConfig.Field[float](
100 doc="Resolution of the internal PSF model relative to the pixel size; "
101 "e.g. 0.5 is equal to 2x oversampling. This affects only the size of "
102 "the PSF model stamp if piffPsfConfigYaml is set.",
103 default=1,
104 )
105 modelSize = pexConfig.Field[int](
106 doc="Internal model size for PIFF (typically odd, but not enforced). "
107 "Partially ignored if piffPsfConfigYaml is set.",
108 default=25,
109 )
110 outlierNSigma = pexConfig.Field[float](
111 doc="n sigma for chisq outlier rejection. "
112 "Ignored if piffPsfConfigYaml is set.",
113 default=4.0
114 )
115 outlierMaxRemove = pexConfig.Field[float](
116 doc="Max fraction of stars to remove as outliers each iteration. "
117 "Ignored if piffPsfConfigYaml is set.",
118 default=0.05
119 )
120 maxSNR = pexConfig.Field[float](
121 doc="Rescale the weight of bright stars such that their SNR is less "
122 "than this value.",
123 default=200.0
124 )
125 zeroWeightMaskBits = pexConfig.ListField[str](
126 doc="List of mask bits for which to set pixel weights to zero.",
127 default=['BAD', 'CR', 'INTRP', 'SAT', 'SUSPECT', 'NO_DATA']
128 )
129 minimumUnmaskedFraction = pexConfig.Field[float](
130 doc="Minimum fraction of unmasked pixels required to use star.",
131 default=0.5
132 )
133 interpolant = pexConfig.Field[str](
134 doc="GalSim interpolant name for Piff to use. "
135 "Options include 'Lanczos(N)', where N is an integer, along with "
136 "galsim.Cubic, galsim.Delta, galsim.Linear, galsim.Nearest, "
137 "galsim.Quintic, and galsim.SincInterpolant. Ignored if "
138 "piffPsfConfigYaml is set.",
139 check=_validateGalsimInterpolant,
140 default="Lanczos(11)",
141 )
142 zerothOrderInterpNotEnoughStars = pexConfig.Field[bool](
143 doc="If True, use zeroth order interpolation if not enough stars are found.",
144 default=False
145 )
146 debugStarData = pexConfig.Field[bool](
147 doc="Include star images used for fitting in PSF model object.",
148 default=False
149 )
150 useCoordinates = pexConfig.ChoiceField[str](
151 doc="Which spatial coordinates to regress against in PSF modeling.",
152 allowed=dict(
153 pixel='Regress against pixel coordinates',
154 field='Regress against field angles',
155 sky='Regress against RA/Dec'
156 ),
157 default='pixel'
158 )
159 piffMaxIter = pexConfig.Field[int](
160 doc="Maximum iteration while doing outlier rejection."
161 "Might be overwrite if zerothOrderInterpNotEnoughStars is True and "
162 "ignore if piffPsfConfigYaml is not None.",
163 default=15
164 )
165 piffLoggingLevel = pexConfig.RangeField[int](
166 doc="PIFF-specific logging level, from 0 (least chatty) to 3 (very verbose); "
167 "note that logs will come out with unrelated notations (e.g. WARNING does "
168 "not imply a warning.)",
169 default=0,
170 min=0,
171 max=3,
172 )
173 piffPsfConfigYaml = pexConfig.Field[str](
174 doc="Configuration file for PIFF in YAML format. This overrides many "
175 "other settings in this config and is not validated ahead of runtime.",
176 default=None,
177 optional=True,
178 )
179
180 def setDefaults(self):
181 super().setDefaults()
182 # stampSize should be at least 25 so that
183 # i) aperture flux with 12 pixel radius can be compared to PSF flux.
184 # ii) fake sources injected to match the 12 pixel aperture flux get
185 # measured correctly
186 # stampSize should also be at least sqrt(2)*modelSize/samplingSize so
187 # that rotated images have support in the model
188
189 self.stampSize = 25
190 # Resize the stamp to accommodate the model, but only if necessary.
191 if self.useCoordinates == "sky":
192 self.stampSize = max(25, 2*int(0.5*self.modelSize*np.sqrt(2)/self.samplingSize) + 1)
193
194 def validate(self):
195 super().validate()
196
197 if (stamp_size := self.stampSize) is not None:
198 model_size = self.modelSize
199 sampling_size = self.samplingSize
200 if self.useCoordinates == 'sky':
201 min_stamp_size = int(np.sqrt(2) * model_size / sampling_size)
202 else:
203 min_stamp_size = int(model_size / sampling_size)
204
205 if stamp_size < min_stamp_size:
206 msg = (f"PIFF model size of {model_size} is too large for stamp size {stamp_size}. "
207 f"Set stampSize >= {min_stamp_size}"
208 )
209 raise pexConfig.FieldValidationError(self.__class__.modelSize, self, msg)
210
211
212def getGoodPixels(maskedImage, zeroWeightMaskBits):
213 """Compute an index array indicating good pixels to use.
214
215 Parameters
216 ----------
217 maskedImage : `afw.image.MaskedImage`
218 PSF candidate postage stamp
219 zeroWeightMaskBits : `List[str]`
220 List of mask bits for which to set pixel weights to zero.
221
222 Returns
223 -------
224 good : `ndarray`
225 Index array indicating good pixels.
226 """
227 imArr = maskedImage.image.array
228 varArr = maskedImage.variance.array
229 bitmask = maskedImage.mask.getPlaneBitMask(zeroWeightMaskBits)
230 good = (
231 (varArr != 0)
232 & (np.isfinite(varArr))
233 & (np.isfinite(imArr))
234 & ((maskedImage.mask.array & bitmask) == 0)
235 )
236 return good
237
238
239def computeWeight(maskedImage, maxSNR, good):
240 """Derive a weight map without Poisson variance component due to signal.
241
242 Parameters
243 ----------
244 maskedImage : `afw.image.MaskedImage`
245 PSF candidate postage stamp
246 maxSNR : `float`
247 Maximum SNR applying variance floor.
248 good : `ndarray`
249 Index array indicating good pixels.
250
251 Returns
252 -------
253 weightArr : `ndarry`
254 Array to use for weight.
255
256 See Also
257 --------
258 `lsst.meas.algorithms.variance_plance.remove_signal_from_variance` :
259 Remove the Poisson contribution from sources in the variance plane of
260 an Exposure.
261 """
262 imArr = maskedImage.image.array
263 varArr = maskedImage.variance.array
264
265 # Fit a straight line to variance vs (sky-subtracted) signal.
266 # The evaluate that line at zero signal to get an estimate of the
267 # signal-free variance.
268 fit = np.polyfit(imArr[good], varArr[good], deg=1)
269 # fit is [1/gain, sky_var]
270 weightArr = np.zeros_like(imArr, dtype=float)
271 weightArr[good] = 1./fit[1]
272
273 applyMaxSNR(imArr, weightArr, good, maxSNR)
274 return weightArr
275
276
277def applyMaxSNR(imArr, weightArr, good, maxSNR):
278 """Rescale weight of bright stars to cap the computed SNR.
279
280 Parameters
281 ----------
282 imArr : `ndarray`
283 Signal (image) array of stamp.
284 weightArr : `ndarray`
285 Weight map array. May be rescaled in place.
286 good : `ndarray`
287 Index array of pixels to use when computing SNR.
288 maxSNR : `float`
289 Threshold for adjusting variance plane implementing maximum SNR.
290 """
291 # We define the SNR value following Piff. Here's the comment from that
292 # code base explaining the calculation.
293 #
294 # The S/N value that we use will be the weighted total flux where the
295 # weight function is the star's profile itself. This is the maximum S/N
296 # value that any flux measurement can possibly produce, which will be
297 # closer to an in-practice S/N than using all the pixels equally.
298 #
299 # F = Sum_i w_i I_i^2
300 # var(F) = Sum_i w_i^2 I_i^2 var(I_i)
301 # = Sum_i w_i I_i^2 <--- Assumes var(I_i) = 1/w_i
302 #
303 # S/N = F / sqrt(var(F))
304 #
305 # Note that if the image is pure noise, this will produce a "signal" of
306 #
307 # F_noise = Sum_i w_i 1/w_i = Npix
308 #
309 # So for a more accurate estimate of the S/N of the actual star itself, one
310 # should subtract off Npix from the measured F.
311 #
312 # The final formula then is:
313 #
314 # F = Sum_i w_i I_i^2
315 # S/N = (F-Npix) / sqrt(F)
316 F = np.sum(weightArr[good]*imArr[good]**2, dtype=float)
317 Npix = np.sum(good)
318 SNR = 0.0 if F < Npix else (F-Npix)/np.sqrt(F)
319 # rescale weight of bright stars. Essentially makes an error floor.
320 if SNR > maxSNR:
321 factor = (maxSNR / SNR)**2
322 weightArr[good] *= factor
323
324
325def _computeWeightAlternative(maskedImage, maxSNR):
326 """Alternative algorithm for creating weight map.
327
328 This version is equivalent to that used by Piff internally. The weight map
329 it produces tends to leave a residual when removing the Poisson component
330 due to the signal. We leave it here as a reference, but without intending
331 that it be used (or be maintained).
332 """
333 imArr = maskedImage.image.array
334 varArr = maskedImage.variance.array
335 good = (varArr != 0) & np.isfinite(varArr) & np.isfinite(imArr)
336
337 fit = np.polyfit(imArr[good], varArr[good], deg=1)
338 # fit is [1/gain, sky_var]
339 gain = 1./fit[0]
340 varArr[good] -= imArr[good] / gain
341 weightArr = np.zeros_like(imArr, dtype=float)
342 weightArr[good] = 1./varArr[good]
343
344 applyMaxSNR(imArr, weightArr, good, maxSNR)
345 return weightArr
346
347
349 """A measurePsfTask PSF estimator using Piff as the implementation.
350 """
351 ConfigClass = PiffPsfDeterminerConfig
352 _DefaultName = "psfDeterminer.Piff"
353
354 def __init__(self, config, schema=None, **kwds):
355 BasePsfDeterminerTask.__init__(self, config, schema=schema, **kwds)
356
357 piffLoggingLevels = {
358 0: logging.CRITICAL,
359 1: logging.WARNING,
360 2: logging.INFO,
361 3: logging.DEBUG,
362 }
363 self.piffLogger = lsst.utils.logging.getLogger(f"{self.log.name}.piff")
364 self.piffLogger.setLevel(piffLoggingLevels[self.config.piffLoggingLevel])
365
367 self, exposure, psfCandidateList, metadata=None, flagKey=None
368 ):
369 """Determine a Piff PSF model for an exposure given a list of PSF
370 candidates.
371
372 Parameters
373 ----------
374 exposure : `lsst.afw.image.Exposure`
375 Exposure containing the PSF candidates.
376 psfCandidateList : `list` of `lsst.meas.algorithms.PsfCandidate`
377 A sequence of PSF candidates typically obtained by detecting sources
378 and then running them through a star selector.
379 metadata : `lsst.daf.base import PropertyList` or `None`, optional
380 A home for interesting tidbits of information.
381 flagKey : `str` or `None`, optional
382 Schema key used to mark sources actually used in PSF determination.
383
384 Returns
385 -------
386 psf : `lsst.meas.extensions.piff.PiffPsf`
387 The measured PSF model.
388 psfCellSet : `None`
389 Unused by this PsfDeterminer.
390 """
391 psfCandidateList = self.downsampleCandidates(psfCandidateList)
392
393 if self.config.stampSize:
394 stampSize = self.config.stampSize
395 if stampSize > psfCandidateList[0].getWidth():
396 self.log.warning("stampSize is larger than the PSF candidate size. Using candidate size.")
397 stampSize = psfCandidateList[0].getWidth()
398 else: # TODO: Only the if block should stay after DM-36311
399 self.log.debug("stampSize not set. Using candidate size.")
400 stampSize = psfCandidateList[0].getWidth()
401
402 scale = exposure.getWcs().getPixelScale(exposure.getBBox().getCenter()).asArcseconds()
403
404 match self.config.useCoordinates:
405 case 'field':
406 detector = exposure.getDetector()
407 pix_to_field = detector.getTransform(PIXELS, FIELD_ANGLE)
408 gswcs = UVWcsWrapper(pix_to_field)
409 pointing = None
410
411 case 'sky':
412 gswcs = CelestialWcsWrapper(exposure.getWcs())
413 skyOrigin = exposure.getWcs().getSkyOrigin()
414 ra = skyOrigin.getLongitude().asDegrees()
415 dec = skyOrigin.getLatitude().asDegrees()
416 pointing = galsim.CelestialCoord(
417 ra*galsim.degrees,
418 dec*galsim.degrees
419 )
420
421 case 'pixel':
422 gswcs = galsim.PixelScale(scale)
423 pointing = None
424
425 stars = []
426 for candidate in psfCandidateList:
427 cmi = candidate.getMaskedImage(stampSize, stampSize)
428 good = getGoodPixels(cmi, self.config.zeroWeightMaskBits)
429 fracGood = np.sum(good)/good.size
430 if fracGood < self.config.minimumUnmaskedFraction:
431 continue
432 weight = computeWeight(cmi, self.config.maxSNR, good)
433
434 bbox = cmi.getBBox()
435 bds = galsim.BoundsI(
436 galsim.PositionI(*bbox.getMin()),
437 galsim.PositionI(*bbox.getMax())
438 )
439 gsImage = galsim.Image(bds, wcs=gswcs, dtype=float)
440 gsImage.array[:] = cmi.image.array
441 gsWeight = galsim.Image(bds, wcs=gswcs, dtype=float)
442 gsWeight.array[:] = weight
443
444 source = candidate.getSource()
445 image_pos = galsim.PositionD(source.getX(), source.getY())
446
447 data = piff.StarData(
448 gsImage,
449 image_pos,
450 weight=gsWeight,
451 pointing=pointing
452 )
453 stars.append(piff.Star(data, None))
454
455 if self.config.piffPsfConfigYaml is None:
456 piffConfig = {
457 'type': 'Simple',
458 'model': {
459 'type': 'PixelGrid',
460 'scale': scale * self.config.samplingSize,
461 'size': self.config.modelSize,
462 'interp': self.config.interpolant,
463 'centered': self.config.piffPixelGridFitCenter,
464 },
465 'interp': {
466 'type': 'BasisPolynomial',
467 'order': self.config.spatialOrder,
468 'solver': self.config.piffBasisPolynomialSolver,
469 },
470 'outliers': {
471 'type': 'Chisq',
472 'nsigma': self.config.outlierNSigma,
473 'max_remove': self.config.outlierMaxRemove,
474 },
475 'max_iter': self.config.piffMaxIter
476 }
477 else:
478 piffConfig = yaml.safe_load(self.config.piffPsfConfigYaml)
479
480 def _get_threshold(nth_order):
481 # number of free parameter in the polynomial interpolation
482 freeParameters = ((nth_order + 1) * (nth_order + 2)) // 2
483 return freeParameters
484
485 if self.config.zerothOrderInterpNotEnoughStars:
486 if piffConfig['interp']['type'] in ['BasisPolynomial', 'Polynomial']:
487 threshold = _get_threshold(piffConfig['interp']['order'])
488 if len(stars) <= threshold:
489 self.log.warning(
490 f"Only {len(stars)} stars found, "
491 f"but {threshold} required. Using zeroth order interpolation."
492 )
493 piffConfig['interp']['order'] = 0
494 # No need to do any outlier rejection assume
495 # PSF to be average of few stars.
496 piffConfig['max_iter'] = 1
497
498 self._piffConfig = piffConfig
499 piffResult = piff.PSF.process(piffConfig)
500 wcs = {0: gswcs}
501
502 piffResult.fit(stars, wcs, pointing, logger=self.piffLogger)
503 drawSize = 2*np.floor(0.5*stampSize/self.config.samplingSize) + 1
504
505 used_image_pos = [s.image_pos for s in piffResult.stars]
506 if flagKey:
507 for candidate in psfCandidateList:
508 source = candidate.getSource()
509 posd = galsim.PositionD(source.getX(), source.getY())
510 if posd in used_image_pos:
511 source.set(flagKey, True)
512
513 if metadata is not None:
514 metadata["spatialFitChi2"] = piffResult.chisq
515 metadata["numAvailStars"] = len(stars)
516 metadata["numGoodStars"] = len(piffResult.stars)
517 metadata["avgX"] = np.mean([p.x for p in piffResult.stars])
518 metadata["avgY"] = np.mean([p.y for p in piffResult.stars])
519
520 if not self.config.debugStarData:
521 for star in piffResult.stars:
522 # Remove large data objects from the stars
523 del star.fit.params
524 del star.fit.params_var
525 del star.fit.A
526 del star.fit.b
527 del star.data.image
528 del star.data.weight
529 del star.data.orig_weight
530
531 return PiffPsf(drawSize, drawSize, piffResult), None
532
533
534measAlg.psfDeterminerRegistry.register("piff", PiffPsfDeterminerTask)
determinePsf(self, exposure, psfCandidateList, metadata=None, flagKey=None)
getGoodPixels(maskedImage, zeroWeightMaskBits)
applyMaxSNR(imArr, weightArr, good, maxSNR)