LSST Applications g070148d5b3+33e5256705,g0d53e28543+25c8b88941,g0da5cf3356+2dd1178308,g1081da9e2a+62d12e78cb,g17e5ecfddb+7e422d6136,g1c76d35bf8+ede3a706f7,g295839609d+225697d880,g2e2c1a68ba+cc1f6f037e,g2ffcdf413f+853cd4dcde,g38293774b4+62d12e78cb,g3b44f30a73+d953f1ac34,g48ccf36440+885b902d19,g4b2f1765b6+7dedbde6d2,g5320a0a9f6+0c5d6105b6,g56b687f8c9+ede3a706f7,g5c4744a4d9+ef6ac23297,g5ffd174ac0+0c5d6105b6,g6075d09f38+66af417445,g667d525e37+2ced63db88,g670421136f+2ced63db88,g71f27ac40c+2ced63db88,g774830318a+463cbe8d1f,g7876bc68e5+1d137996f1,g7985c39107+62d12e78cb,g7fdac2220c+0fd8241c05,g96f01af41f+368e6903a7,g9ca82378b8+2ced63db88,g9d27549199+ef6ac23297,gabe93b2c52+e3573e3735,gb065e2a02a+3dfbe639da,gbc3249ced9+0c5d6105b6,gbec6a3398f+0c5d6105b6,gc9534b9d65+35b9f25267,gd01420fc67+0c5d6105b6,geee7ff78d7+a14128c129,gf63283c776+ede3a706f7,gfed783d017+0c5d6105b6,w.2022.47
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
28
29import lsst.pex.config as pexConfig
30import lsst.meas.algorithms as measAlg
31from lsst.meas.algorithms.psfDeterminer import BasePsfDeterminerTask
32from .piffPsf import PiffPsf
33
34
35def _validateGalsimInterpolant(name: str) -> bool:
36 """A helper function to validate the GalSim interpolant at config time.
37
38 Parameters
39 ----------
40 name : str
41 The name of the interpolant to use from GalSim. Valid options are:
42 galsim.Lanczos(N) or Lancsos(N), where N is a positive integer
43 galsim.Linear
44 galsim.Cubic
45 galsim.Quintic
46 galsim.Delta
47 galsim.Nearest
48 galsim.SincInterpolant
49
50 Returns
51 -------
52 is_valid : bool
53 Whether the provided interpolant name is valid.
54 """
55 # First, check if ``name`` is a valid Lanczos interpolant.
56 for pattern in (re.compile(r"Lanczos\‍(\d+\‍)"), re.compile(r"galsim.Lanczos\‍(\d+\‍)"),):
57 match = re.match(pattern, name) # Search from the start of the string.
58 if match is not None:
59 # Check that the pattern is also the end of the string.
60 return match.end() == len(name)
61
62 # If not, check if ``name`` is any other valid GalSim interpolant.
63 names = {f"galsim.{interp}" for interp in
64 ("Cubic", "Delta", "Linear", "Nearest", "Quintic", "SincInterpolant")
65 }
66 return name in names
67
68
69class PiffPsfDeterminerConfig(BasePsfDeterminerTask.ConfigClass):
70 spatialOrder = pexConfig.Field[int](
71 doc="specify spatial order for PSF kernel creation",
72 default=2,
73 )
74 samplingSize = pexConfig.Field[float](
75 doc="Resolution of the internal PSF model relative to the pixel size; "
76 "e.g. 0.5 is equal to 2x oversampling",
77 default=1,
78 )
79 outlierNSigma = pexConfig.Field[float](
80 doc="n sigma for chisq outlier rejection",
81 default=4.0
82 )
83 outlierMaxRemove = pexConfig.Field[float](
84 doc="Max fraction of stars to remove as outliers each iteration",
85 default=0.05
86 )
87 maxSNR = pexConfig.Field[float](
88 doc="Rescale the weight of bright stars such that their SNR is less "
89 "than this value.",
90 default=200.0
91 )
92 zeroWeightMaskBits = pexConfig.ListField[str](
93 doc="List of mask bits for which to set pixel weights to zero.",
94 default=['BAD', 'CR', 'INTRP', 'SAT', 'SUSPECT', 'NO_DATA']
95 )
96 minimumUnmaskedFraction = pexConfig.Field[float](
97 doc="Minimum fraction of unmasked pixels required to use star.",
98 default=0.5
99 )
100 interpolant = pexConfig.Field[str](
101 doc="GalSim interpolant name for Piff to use. "
102 "Options include 'Lanczos(N)', where N is an integer, along with "
103 "galsim.Cubic, galsim.Delta, galsim.Linear, galsim.Nearest, "
104 "galsim.Quintic, and galsim.SincInterpolant.",
105 check=_validateGalsimInterpolant,
106 default="Lanczos(11)",
107 )
108 debugStarData = pexConfig.Field[bool](
109 doc="Include star images used for fitting in PSF model object.",
110 default=False
111 )
112
113 def setDefaults(self):
114 super().setDefaults()
115 # stampSize should be at least 25 so that
116 # i) aperture flux with 12 pixel radius can be compared to PSF flux.
117 # ii) fake sources injected to match the 12 pixel aperture flux get
118 # measured correctly
119 self.stampSize = 25
120
121
122def getGoodPixels(maskedImage, zeroWeightMaskBits):
123 """Compute an index array indicating good pixels to use.
124
125 Parameters
126 ----------
127 maskedImage : `afw.image.MaskedImage`
128 PSF candidate postage stamp
129 zeroWeightMaskBits : `List[str]`
130 List of mask bits for which to set pixel weights to zero.
131
132 Returns
133 -------
134 good : `ndarray`
135 Index array indicating good pixels.
136 """
137 imArr = maskedImage.image.array
138 varArr = maskedImage.variance.array
139 bitmask = maskedImage.mask.getPlaneBitMask(zeroWeightMaskBits)
140 good = (
141 (varArr != 0)
142 & (np.isfinite(varArr))
143 & (np.isfinite(imArr))
144 & ((maskedImage.mask.array & bitmask) == 0)
145 )
146 return good
147
148
149def computeWeight(maskedImage, maxSNR, good):
150 """Derive a weight map without Poisson variance component due to signal.
151
152 Parameters
153 ----------
154 maskedImage : `afw.image.MaskedImage`
155 PSF candidate postage stamp
156 maxSNR : `float`
157 Maximum SNR applying variance floor.
158 good : `ndarray`
159 Index array indicating good pixels.
160
161 Returns
162 -------
163 weightArr : `ndarry`
164 Array to use for weight.
165 """
166 imArr = maskedImage.image.array
167 varArr = maskedImage.variance.array
168
169 # Fit a straight line to variance vs (sky-subtracted) signal.
170 # The evaluate that line at zero signal to get an estimate of the
171 # signal-free variance.
172 fit = np.polyfit(imArr[good], varArr[good], deg=1)
173 # fit is [1/gain, sky_var]
174 weightArr = np.zeros_like(imArr, dtype=float)
175 weightArr[good] = 1./fit[1]
176
177 applyMaxSNR(imArr, weightArr, good, maxSNR)
178 return weightArr
179
180
181def applyMaxSNR(imArr, weightArr, good, maxSNR):
182 """Rescale weight of bright stars to cap the computed SNR.
183
184 Parameters
185 ----------
186 imArr : `ndarray`
187 Signal (image) array of stamp.
188 weightArr : `ndarray`
189 Weight map array. May be rescaled in place.
190 good : `ndarray`
191 Index array of pixels to use when computing SNR.
192 maxSNR : `float`
193 Threshold for adjusting variance plane implementing maximum SNR.
194 """
195 # We define the SNR value following Piff. Here's the comment from that
196 # code base explaining the calculation.
197 #
198 # The S/N value that we use will be the weighted total flux where the
199 # weight function is the star's profile itself. This is the maximum S/N
200 # value that any flux measurement can possibly produce, which will be
201 # closer to an in-practice S/N than using all the pixels equally.
202 #
203 # F = Sum_i w_i I_i^2
204 # var(F) = Sum_i w_i^2 I_i^2 var(I_i)
205 # = Sum_i w_i I_i^2 <--- Assumes var(I_i) = 1/w_i
206 #
207 # S/N = F / sqrt(var(F))
208 #
209 # Note that if the image is pure noise, this will produce a "signal" of
210 #
211 # F_noise = Sum_i w_i 1/w_i = Npix
212 #
213 # So for a more accurate estimate of the S/N of the actual star itself, one
214 # should subtract off Npix from the measured F.
215 #
216 # The final formula then is:
217 #
218 # F = Sum_i w_i I_i^2
219 # S/N = (F-Npix) / sqrt(F)
220 F = np.sum(weightArr[good]*imArr[good]**2, dtype=float)
221 Npix = np.sum(good)
222 SNR = 0.0 if F < Npix else (F-Npix)/np.sqrt(F)
223 # rescale weight of bright stars. Essentially makes an error floor.
224 if SNR > maxSNR:
225 factor = (maxSNR / SNR)**2
226 weightArr[good] *= factor
227
228
229def _computeWeightAlternative(maskedImage, maxSNR):
230 """Alternative algorithm for creating weight map.
231
232 This version is equivalent to that used by Piff internally. The weight map
233 it produces tends to leave a residual when removing the Poisson component
234 due to the signal. We leave it here as a reference, but without intending
235 that it be used (or be maintained).
236 """
237 imArr = maskedImage.image.array
238 varArr = maskedImage.variance.array
239 good = (varArr != 0) & np.isfinite(varArr) & np.isfinite(imArr)
240
241 fit = np.polyfit(imArr[good], varArr[good], deg=1)
242 # fit is [1/gain, sky_var]
243 gain = 1./fit[0]
244 varArr[good] -= imArr[good] / gain
245 weightArr = np.zeros_like(imArr, dtype=float)
246 weightArr[good] = 1./varArr[good]
247
248 applyMaxSNR(imArr, weightArr, good, maxSNR)
249 return weightArr
250
251
253 """A measurePsfTask PSF estimator using Piff as the implementation.
254 """
255 ConfigClass = PiffPsfDeterminerConfig
256 _DefaultName = "psfDeterminer.Piff"
257
259 self, exposure, psfCandidateList, metadata=None, flagKey=None
260 ):
261 """Determine a Piff PSF model for an exposure given a list of PSF
262 candidates.
263
264 Parameters
265 ----------
266 exposure : `lsst.afw.image.Exposure`
267 Exposure containing the PSF candidates.
268 psfCandidateList : `list` of `lsst.meas.algorithms.PsfCandidate`
269 A sequence of PSF candidates typically obtained by detecting sources
270 and then running them through a star selector.
271 metadata : `lsst.daf.base import PropertyList` or `None`, optional
272 A home for interesting tidbits of information.
273 flagKey : `str` or `None`, optional
274 Schema key used to mark sources actually used in PSF determination.
275
276 Returns
277 -------
278 psf : `lsst.meas.extensions.piff.PiffPsf`
279 The measured PSF model.
280 psfCellSet : `None`
281 Unused by this PsfDeterminer.
282 """
283 if self.config.stampSize:
284 stampSize = self.config.stampSize
285 if stampSize > psfCandidateList[0].getWidth():
286 self.log.warning("stampSize is larger than the PSF candidate size. Using candidate size.")
287 stampSize = psfCandidateList[0].getWidth()
288 else: # TODO: Only the if block should stay after DM-36311
289 self.log.debug("stampSize not set. Using candidate size.")
290 stampSize = psfCandidateList[0].getWidth()
291
292 self._validatePsfCandidates(psfCandidateList, stampSize)
293
294 stars = []
295 for candidate in psfCandidateList:
296 cmi = candidate.getMaskedImage(stampSize, stampSize)
297 good = getGoodPixels(cmi, self.config.zeroWeightMaskBits)
298 fracGood = np.sum(good)/good.size
299 if fracGood < self.config.minimumUnmaskedFraction:
300 continue
301 weight = computeWeight(cmi, self.config.maxSNR, good)
302
303 bbox = cmi.getBBox()
304 bds = galsim.BoundsI(
305 galsim.PositionI(*bbox.getMin()),
306 galsim.PositionI(*bbox.getMax())
307 )
308 gsImage = galsim.Image(bds, scale=1.0, dtype=float)
309 gsImage.array[:] = cmi.image.array
310 gsWeight = galsim.Image(bds, scale=1.0, dtype=float)
311 gsWeight.array[:] = weight
312
313 source = candidate.getSource()
314 image_pos = galsim.PositionD(source.getX(), source.getY())
315
316 data = piff.StarData(
317 gsImage,
318 image_pos,
319 weight=gsWeight
320 )
321 stars.append(piff.Star(data, None))
322
323 piffConfig = {
324 'type': "Simple",
325 'model': {
326 'type': 'PixelGrid',
327 'scale': self.config.samplingSize,
328 'size': stampSize,
329 'interp': self.config.interpolant
330 },
331 'interp': {
332 'type': 'BasisPolynomial',
333 'order': self.config.spatialOrder
334 },
335 'outliers': {
336 'type': 'Chisq',
337 'nsigma': self.config.outlierNSigma,
338 'max_remove': self.config.outlierMaxRemove
339 }
340 }
341
342 piffResult = piff.PSF.process(piffConfig)
343 # Run on a single CCD, and in image coords rather than sky coords.
344 wcs = {0: galsim.PixelScale(1.0)}
345 pointing = None
346
347 piffResult.fit(stars, wcs, pointing, logger=self.log)
348 drawSize = 2*np.floor(0.5*stampSize/self.config.samplingSize) + 1
349
350 used_image_pos = [s.image_pos for s in piffResult.stars]
351 if flagKey:
352 for candidate in psfCandidateList:
353 source = candidate.getSource()
354 posd = galsim.PositionD(source.getX(), source.getY())
355 if posd in used_image_pos:
356 source.set(flagKey, True)
357
358 if metadata is not None:
359 metadata["spatialFitChi2"] = piffResult.chisq
360 metadata["numAvailStars"] = len(stars)
361 metadata["numGoodStars"] = len(piffResult.stars)
362 metadata["avgX"] = np.mean([p.x for p in piffResult.stars])
363 metadata["avgY"] = np.mean([p.y for p in piffResult.stars])
364
365 if not self.config.debugStarData:
366 for star in piffResult.stars:
367 # Remove large data objects from the stars
368 del star.fit.params
369 del star.fit.params_var
370 del star.fit.A
371 del star.fit.b
372 del star.data.image
373 del star.data.weight
374 del star.data.orig_weight
375
376 return PiffPsf(drawSize, drawSize, piffResult), None
377
378 # TODO: DM-36311: This method can be removed.
379 @staticmethod
380 def _validatePsfCandidates(psfCandidateList, stampSize):
381 """Raise if psfCandidates are smaller than the configured kernelSize.
382
383 Parameters
384 ----------
385 psfCandidateList : `list` of `lsst.meas.algorithms.PsfCandidate`
386 Sequence of psf candidates to check.
387 stampSize : `int`
388 Size of image model to use in PIFF.
389
390 Raises
391 ------
392 RuntimeError
393 Raised if any psfCandidate has width or height smaller than
394 ``stampSize``.
395 """
396 # All candidates will necessarily have the same dimensions.
397 candidate = psfCandidateList[0]
398 if (candidate.getHeight() < stampSize
399 or candidate.getWidth() < stampSize):
400 raise RuntimeError(f"PSF candidates must be at least {stampSize=} pixels per side; "
401 f"found {candidate.getWidth()}x{candidate.getHeight()}."
402 )
403
404
405measAlg.psfDeterminerRegistry.register("piff", PiffPsfDeterminerTask)
A class to contain the data, WCS, and other information needed to describe an image of the sky.
Definition: Exposure.h:72
A class to manipulate images, masks, and variance as a single object.
Definition: MaskedImage.h:74
Class stored in SpatialCells for spatial Psf fitting.
Definition: PsfCandidate.h:55
def determinePsf(self, exposure, psfCandidateList, metadata=None, flagKey=None)
def getGoodPixels(maskedImage, zeroWeightMaskBits)
def applyMaxSNR(imArr, weightArr, good, maxSNR)
def computeWeight(maskedImage, maxSNR, good)