LSST Applications g063fba187b+fee0456c91,g0f08755f38+ea96e5a5a3,g1653933729+a8ce1bb630,g168dd56ebc+a8ce1bb630,g1a2382251a+90257ff92a,g20f6ffc8e0+ea96e5a5a3,g217e2c1bcf+937a289c59,g28da252d5a+daa7da44eb,g2bbee38e9b+253935c60e,g2bc492864f+253935c60e,g3156d2b45e+6e55a43351,g32e5bea42b+31359a2a7a,g347aa1857d+253935c60e,g35bb328faa+a8ce1bb630,g3a166c0a6a+253935c60e,g3b1af351f3+a8ce1bb630,g3e281a1b8c+c5dd892a6c,g414038480c+416496e02f,g41af890bb2+afe91b1188,g599934f4f4+0db33f7991,g7af13505b9+e36de7bce6,g80478fca09+da231ba887,g82479be7b0+a4516e59e3,g858d7b2824+ea96e5a5a3,g89c8672015+f4add4ffd5,g9125e01d80+a8ce1bb630,ga5288a1d22+bc6ab8dfbd,gb58c049af0+d64f4d3760,gc28159a63d+253935c60e,gcab2d0539d+3f2b72788c,gcf0d15dbbd+4ea9c45075,gda6a2b7d83+4ea9c45075,gdaeeff99f8+1711a396fd,ge79ae78c31+253935c60e,gef2f8181fd+3031e3cf99,gf0baf85859+c1f95f4921,gfa517265be+ea96e5a5a3,gfa999e8aa5+17cd334064,w.2024.50
LSST Data Management Base Package
Loading...
Searching...
No Matches
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 samplingSize = pexConfig.Field[float](
81 doc="Resolution of the internal PSF model relative to the pixel size; "
82 "e.g. 0.5 is equal to 2x oversampling. This affects only the size of "
83 "the PSF model stamp if piffPsfConfigYaml is set.",
84 default=1,
85 )
86 modelSize = pexConfig.Field[int](
87 doc="Internal model size for PIFF (typically odd, but not enforced). "
88 "Partially ignored if piffPsfConfigYaml is set.",
89 default=25,
90 )
91 outlierNSigma = pexConfig.Field[float](
92 doc="n sigma for chisq outlier rejection. "
93 "Ignored if piffPsfConfigYaml is set.",
94 default=4.0
95 )
96 outlierMaxRemove = pexConfig.Field[float](
97 doc="Max fraction of stars to remove as outliers each iteration. "
98 "Ignored if piffPsfConfigYaml is set.",
99 default=0.05
100 )
101 maxSNR = pexConfig.Field[float](
102 doc="Rescale the weight of bright stars such that their SNR is less "
103 "than this value.",
104 default=200.0
105 )
106 zeroWeightMaskBits = pexConfig.ListField[str](
107 doc="List of mask bits for which to set pixel weights to zero.",
108 default=['BAD', 'CR', 'INTRP', 'SAT', 'SUSPECT', 'NO_DATA']
109 )
110 minimumUnmaskedFraction = pexConfig.Field[float](
111 doc="Minimum fraction of unmasked pixels required to use star.",
112 default=0.5
113 )
114 interpolant = pexConfig.Field[str](
115 doc="GalSim interpolant name for Piff to use. "
116 "Options include 'Lanczos(N)', where N is an integer, along with "
117 "galsim.Cubic, galsim.Delta, galsim.Linear, galsim.Nearest, "
118 "galsim.Quintic, and galsim.SincInterpolant. Ignored if "
119 "piffPsfConfigYaml is set.",
120 check=_validateGalsimInterpolant,
121 default="Lanczos(11)",
122 )
123 debugStarData = pexConfig.Field[bool](
124 doc="Include star images used for fitting in PSF model object.",
125 default=False
126 )
127 useCoordinates = pexConfig.ChoiceField[str](
128 doc="Which spatial coordinates to regress against in PSF modeling.",
129 allowed=dict(
130 pixel='Regress against pixel coordinates',
131 field='Regress against field angles',
132 sky='Regress against RA/Dec'
133 ),
134 default='pixel'
135 )
136 piffLoggingLevel = pexConfig.RangeField[int](
137 doc="PIFF-specific logging level, from 0 (least chatty) to 3 (very verbose); "
138 "note that logs will come out with unrelated notations (e.g. WARNING does "
139 "not imply a warning.)",
140 default=0,
141 min=0,
142 max=3,
143 )
144 piffPsfConfigYaml = pexConfig.Field[str](
145 doc="Configuration file for PIFF in YAML format. This overrides many "
146 "other settings in this config and is not validated ahead of runtime.",
147 default=None,
148 optional=True,
149 )
150
151 def setDefaults(self):
152 super().setDefaults()
153 # stampSize should be at least 25 so that
154 # i) aperture flux with 12 pixel radius can be compared to PSF flux.
155 # ii) fake sources injected to match the 12 pixel aperture flux get
156 # measured correctly
157 # stampSize should also be at least sqrt(2)*modelSize/samplingSize so
158 # that rotated images have support in the model
159
160 self.stampSize = 25
161 # Resize the stamp to accommodate the model, but only if necessary.
162 if self.useCoordinates == "sky":
163 self.stampSize = max(25, 2*int(0.5*self.modelSize*np.sqrt(2)/self.samplingSize) + 1)
164
165 def validate(self):
166 super().validate()
167
168 if (stamp_size := self.stampSize) is not None:
169 model_size = self.modelSize
170 sampling_size = self.samplingSize
171 if self.useCoordinates == 'sky':
172 min_stamp_size = int(np.sqrt(2) * model_size / sampling_size)
173 else:
174 min_stamp_size = int(model_size / sampling_size)
175
176 if stamp_size < min_stamp_size:
177 msg = (f"PIFF model size of {model_size} is too large for stamp size {stamp_size}. "
178 f"Set stampSize >= {min_stamp_size}"
179 )
180 raise pexConfig.FieldValidationError(self.__class__.modelSize, self, msg)
181
182
183def getGoodPixels(maskedImage, zeroWeightMaskBits):
184 """Compute an index array indicating good pixels to use.
185
186 Parameters
187 ----------
188 maskedImage : `afw.image.MaskedImage`
189 PSF candidate postage stamp
190 zeroWeightMaskBits : `List[str]`
191 List of mask bits for which to set pixel weights to zero.
192
193 Returns
194 -------
195 good : `ndarray`
196 Index array indicating good pixels.
197 """
198 imArr = maskedImage.image.array
199 varArr = maskedImage.variance.array
200 bitmask = maskedImage.mask.getPlaneBitMask(zeroWeightMaskBits)
201 good = (
202 (varArr != 0)
203 & (np.isfinite(varArr))
204 & (np.isfinite(imArr))
205 & ((maskedImage.mask.array & bitmask) == 0)
206 )
207 return good
208
209
210def computeWeight(maskedImage, maxSNR, good):
211 """Derive a weight map without Poisson variance component due to signal.
212
213 Parameters
214 ----------
215 maskedImage : `afw.image.MaskedImage`
216 PSF candidate postage stamp
217 maxSNR : `float`
218 Maximum SNR applying variance floor.
219 good : `ndarray`
220 Index array indicating good pixels.
221
222 Returns
223 -------
224 weightArr : `ndarry`
225 Array to use for weight.
226
227 See Also
228 --------
229 `lsst.meas.algorithms.variance_plance.remove_signal_from_variance` :
230 Remove the Poisson contribution from sources in the variance plane of
231 an Exposure.
232 """
233 imArr = maskedImage.image.array
234 varArr = maskedImage.variance.array
235
236 # Fit a straight line to variance vs (sky-subtracted) signal.
237 # The evaluate that line at zero signal to get an estimate of the
238 # signal-free variance.
239 fit = np.polyfit(imArr[good], varArr[good], deg=1)
240 # fit is [1/gain, sky_var]
241 weightArr = np.zeros_like(imArr, dtype=float)
242 weightArr[good] = 1./fit[1]
243
244 applyMaxSNR(imArr, weightArr, good, maxSNR)
245 return weightArr
246
247
248def applyMaxSNR(imArr, weightArr, good, maxSNR):
249 """Rescale weight of bright stars to cap the computed SNR.
250
251 Parameters
252 ----------
253 imArr : `ndarray`
254 Signal (image) array of stamp.
255 weightArr : `ndarray`
256 Weight map array. May be rescaled in place.
257 good : `ndarray`
258 Index array of pixels to use when computing SNR.
259 maxSNR : `float`
260 Threshold for adjusting variance plane implementing maximum SNR.
261 """
262 # We define the SNR value following Piff. Here's the comment from that
263 # code base explaining the calculation.
264 #
265 # The S/N value that we use will be the weighted total flux where the
266 # weight function is the star's profile itself. This is the maximum S/N
267 # value that any flux measurement can possibly produce, which will be
268 # closer to an in-practice S/N than using all the pixels equally.
269 #
270 # F = Sum_i w_i I_i^2
271 # var(F) = Sum_i w_i^2 I_i^2 var(I_i)
272 # = Sum_i w_i I_i^2 <--- Assumes var(I_i) = 1/w_i
273 #
274 # S/N = F / sqrt(var(F))
275 #
276 # Note that if the image is pure noise, this will produce a "signal" of
277 #
278 # F_noise = Sum_i w_i 1/w_i = Npix
279 #
280 # So for a more accurate estimate of the S/N of the actual star itself, one
281 # should subtract off Npix from the measured F.
282 #
283 # The final formula then is:
284 #
285 # F = Sum_i w_i I_i^2
286 # S/N = (F-Npix) / sqrt(F)
287 F = np.sum(weightArr[good]*imArr[good]**2, dtype=float)
288 Npix = np.sum(good)
289 SNR = 0.0 if F < Npix else (F-Npix)/np.sqrt(F)
290 # rescale weight of bright stars. Essentially makes an error floor.
291 if SNR > maxSNR:
292 factor = (maxSNR / SNR)**2
293 weightArr[good] *= factor
294
295
296def _computeWeightAlternative(maskedImage, maxSNR):
297 """Alternative algorithm for creating weight map.
298
299 This version is equivalent to that used by Piff internally. The weight map
300 it produces tends to leave a residual when removing the Poisson component
301 due to the signal. We leave it here as a reference, but without intending
302 that it be used (or be maintained).
303 """
304 imArr = maskedImage.image.array
305 varArr = maskedImage.variance.array
306 good = (varArr != 0) & np.isfinite(varArr) & np.isfinite(imArr)
307
308 fit = np.polyfit(imArr[good], varArr[good], deg=1)
309 # fit is [1/gain, sky_var]
310 gain = 1./fit[0]
311 varArr[good] -= imArr[good] / gain
312 weightArr = np.zeros_like(imArr, dtype=float)
313 weightArr[good] = 1./varArr[good]
314
315 applyMaxSNR(imArr, weightArr, good, maxSNR)
316 return weightArr
317
318
320 """A measurePsfTask PSF estimator using Piff as the implementation.
321 """
322 ConfigClass = PiffPsfDeterminerConfig
323 _DefaultName = "psfDeterminer.Piff"
324
325 def __init__(self, config, schema=None, **kwds):
326 BasePsfDeterminerTask.__init__(self, config, schema=schema, **kwds)
327
328 piffLoggingLevels = {
329 0: logging.CRITICAL,
330 1: logging.WARNING,
331 2: logging.INFO,
332 3: logging.DEBUG,
333 }
334 self.piffLogger = lsst.utils.logging.getLogger(f"{self.log.name}.piff")
335 self.piffLogger.setLevel(piffLoggingLevels[self.config.piffLoggingLevel])
336
338 self, exposure, psfCandidateList, metadata=None, flagKey=None
339 ):
340 """Determine a Piff PSF model for an exposure given a list of PSF
341 candidates.
342
343 Parameters
344 ----------
345 exposure : `lsst.afw.image.Exposure`
346 Exposure containing the PSF candidates.
347 psfCandidateList : `list` of `lsst.meas.algorithms.PsfCandidate`
348 A sequence of PSF candidates typically obtained by detecting sources
349 and then running them through a star selector.
350 metadata : `lsst.daf.base import PropertyList` or `None`, optional
351 A home for interesting tidbits of information.
352 flagKey : `str` or `None`, optional
353 Schema key used to mark sources actually used in PSF determination.
354
355 Returns
356 -------
357 psf : `lsst.meas.extensions.piff.PiffPsf`
358 The measured PSF model.
359 psfCellSet : `None`
360 Unused by this PsfDeterminer.
361 """
362 psfCandidateList = self.downsampleCandidates(psfCandidateList)
363
364 if self.config.stampSize:
365 stampSize = self.config.stampSize
366 if stampSize > psfCandidateList[0].getWidth():
367 self.log.warning("stampSize is larger than the PSF candidate size. Using candidate size.")
368 stampSize = psfCandidateList[0].getWidth()
369 else: # TODO: Only the if block should stay after DM-36311
370 self.log.debug("stampSize not set. Using candidate size.")
371 stampSize = psfCandidateList[0].getWidth()
372
373 scale = exposure.getWcs().getPixelScale(exposure.getBBox().getCenter()).asArcseconds()
374
375 match self.config.useCoordinates:
376 case 'field':
377 detector = exposure.getDetector()
378 pix_to_field = detector.getTransform(PIXELS, FIELD_ANGLE)
379 gswcs = UVWcsWrapper(pix_to_field)
380 pointing = None
381
382 case 'sky':
383 gswcs = CelestialWcsWrapper(exposure.getWcs())
384 skyOrigin = exposure.getWcs().getSkyOrigin()
385 ra = skyOrigin.getLongitude().asDegrees()
386 dec = skyOrigin.getLatitude().asDegrees()
387 pointing = galsim.CelestialCoord(
388 ra*galsim.degrees,
389 dec*galsim.degrees
390 )
391
392 case 'pixel':
393 gswcs = galsim.PixelScale(scale)
394 pointing = None
395
396 stars = []
397 for candidate in psfCandidateList:
398 cmi = candidate.getMaskedImage(stampSize, stampSize)
399 good = getGoodPixels(cmi, self.config.zeroWeightMaskBits)
400 fracGood = np.sum(good)/good.size
401 if fracGood < self.config.minimumUnmaskedFraction:
402 continue
403 weight = computeWeight(cmi, self.config.maxSNR, good)
404
405 bbox = cmi.getBBox()
406 bds = galsim.BoundsI(
407 galsim.PositionI(*bbox.getMin()),
408 galsim.PositionI(*bbox.getMax())
409 )
410 gsImage = galsim.Image(bds, wcs=gswcs, dtype=float)
411 gsImage.array[:] = cmi.image.array
412 gsWeight = galsim.Image(bds, wcs=gswcs, dtype=float)
413 gsWeight.array[:] = weight
414
415 source = candidate.getSource()
416 image_pos = galsim.PositionD(source.getX(), source.getY())
417
418 data = piff.StarData(
419 gsImage,
420 image_pos,
421 weight=gsWeight,
422 pointing=pointing
423 )
424 stars.append(piff.Star(data, None))
425
426 if self.config.piffPsfConfigYaml is None:
427 piffConfig = {
428 'type': 'Simple',
429 'model': {
430 'type': 'PixelGrid',
431 'scale': scale * self.config.samplingSize,
432 'size': self.config.modelSize,
433 'interp': self.config.interpolant,
434 },
435 'interp': {
436 'type': 'BasisPolynomial',
437 'order': self.config.spatialOrder,
438 },
439 'outliers': {
440 'type': 'Chisq',
441 'nsigma': self.config.outlierNSigma,
442 'max_remove': self.config.outlierMaxRemove,
443 }
444 }
445 else:
446 piffConfig = yaml.safe_load(self.config.piffPsfConfigYaml)
447
448 piffResult = piff.PSF.process(piffConfig)
449 wcs = {0: gswcs}
450
451 piffResult.fit(stars, wcs, pointing, logger=self.piffLogger)
452 drawSize = 2*np.floor(0.5*stampSize/self.config.samplingSize) + 1
453
454 used_image_pos = [s.image_pos for s in piffResult.stars]
455 if flagKey:
456 for candidate in psfCandidateList:
457 source = candidate.getSource()
458 posd = galsim.PositionD(source.getX(), source.getY())
459 if posd in used_image_pos:
460 source.set(flagKey, True)
461
462 if metadata is not None:
463 metadata["spatialFitChi2"] = piffResult.chisq
464 metadata["numAvailStars"] = len(stars)
465 metadata["numGoodStars"] = len(piffResult.stars)
466 metadata["avgX"] = np.mean([p.x for p in piffResult.stars])
467 metadata["avgY"] = np.mean([p.y for p in piffResult.stars])
468
469 if not self.config.debugStarData:
470 for star in piffResult.stars:
471 # Remove large data objects from the stars
472 del star.fit.params
473 del star.fit.params_var
474 del star.fit.A
475 del star.fit.b
476 del star.data.image
477 del star.data.weight
478 del star.data.orig_weight
479
480 return PiffPsf(drawSize, drawSize, piffResult), None
481
482
483measAlg.psfDeterminerRegistry.register("piff", PiffPsfDeterminerTask)
int max
determinePsf(self, exposure, psfCandidateList, metadata=None, flagKey=None)
getGoodPixels(maskedImage, zeroWeightMaskBits)
applyMaxSNR(imArr, weightArr, good, maxSNR)