LSST Applications 27.0.0,g0265f82a02+469cd937ee,g02d81e74bb+21ad69e7e1,g1470d8bcf6+cbe83ee85a,g2079a07aa2+e67c6346a6,g212a7c68fe+04a9158687,g2305ad1205+94392ce272,g295015adf3+81dd352a9d,g2bbee38e9b+469cd937ee,g337abbeb29+469cd937ee,g3939d97d7f+72a9f7b576,g487adcacf7+71499e7cba,g50ff169b8f+5929b3527e,g52b1c1532d+a6fc98d2e7,g591dd9f2cf+df404f777f,g5a732f18d5+be83d3ecdb,g64a986408d+21ad69e7e1,g858d7b2824+21ad69e7e1,g8a8a8dda67+a6fc98d2e7,g99cad8db69+f62e5b0af5,g9ddcbc5298+d4bad12328,ga1e77700b3+9c366c4306,ga8c6da7877+71e4819109,gb0e22166c9+25ba2f69a1,gb6a65358fc+469cd937ee,gbb8dafda3b+69d3c0e320,gc07e1c2157+a98bf949bb,gc120e1dc64+615ec43309,gc28159a63d+469cd937ee,gcf0d15dbbd+72a9f7b576,gdaeeff99f8+a38ce5ea23,ge6526c86ff+3a7c1ac5f1,ge79ae78c31+469cd937ee,gee10cc3b42+a6fc98d2e7,gf1cff7945b+21ad69e7e1,gfbcc870c63+9a11dc8c8f
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
29from lsst.afw.cameraGeom import PIXELS, FIELD_ANGLE
30import lsst.pex.config as pexConfig
31import lsst.meas.algorithms as measAlg
32from lsst.meas.algorithms.psfDeterminer import BasePsfDeterminerTask
33from .piffPsf import PiffPsf
34from .wcs_wrapper import CelestialWcsWrapper, UVWcsWrapper
35
36
37def _validateGalsimInterpolant(name: str) -> bool:
38 """A helper function to validate the GalSim interpolant at config time.
39
40 Parameters
41 ----------
42 name : str
43 The name of the interpolant to use from GalSim. Valid options are:
44 galsim.Lanczos(N) or Lancsos(N), where N is a positive integer
45 galsim.Linear
46 galsim.Cubic
47 galsim.Quintic
48 galsim.Delta
49 galsim.Nearest
50 galsim.SincInterpolant
51
52 Returns
53 -------
54 is_valid : bool
55 Whether the provided interpolant name is valid.
56 """
57 # First, check if ``name`` is a valid Lanczos interpolant.
58 for pattern in (re.compile(r"Lanczos\‍(\d+\‍)"), re.compile(r"galsim.Lanczos\‍(\d+\‍)"),):
59 match = re.match(pattern, name) # Search from the start of the string.
60 if match is not None:
61 # Check that the pattern is also the end of the string.
62 return match.end() == len(name)
63
64 # If not, check if ``name`` is any other valid GalSim interpolant.
65 names = {f"galsim.{interp}" for interp in
66 ("Cubic", "Delta", "Linear", "Nearest", "Quintic", "SincInterpolant")
67 }
68 return name in names
69
70
71class PiffPsfDeterminerConfig(BasePsfDeterminerTask.ConfigClass):
72 spatialOrder = pexConfig.Field[int](
73 doc="specify spatial order for PSF kernel creation",
74 default=2,
75 )
76 samplingSize = pexConfig.Field[float](
77 doc="Resolution of the internal PSF model relative to the pixel size; "
78 "e.g. 0.5 is equal to 2x oversampling",
79 default=1,
80 )
81 outlierNSigma = pexConfig.Field[float](
82 doc="n sigma for chisq outlier rejection",
83 default=4.0
84 )
85 outlierMaxRemove = pexConfig.Field[float](
86 doc="Max fraction of stars to remove as outliers each iteration",
87 default=0.05
88 )
89 maxSNR = pexConfig.Field[float](
90 doc="Rescale the weight of bright stars such that their SNR is less "
91 "than this value.",
92 default=200.0
93 )
94 zeroWeightMaskBits = pexConfig.ListField[str](
95 doc="List of mask bits for which to set pixel weights to zero.",
96 default=['BAD', 'CR', 'INTRP', 'SAT', 'SUSPECT', 'NO_DATA']
97 )
98 minimumUnmaskedFraction = pexConfig.Field[float](
99 doc="Minimum fraction of unmasked pixels required to use star.",
100 default=0.5
101 )
102 interpolant = pexConfig.Field[str](
103 doc="GalSim interpolant name for Piff to use. "
104 "Options include 'Lanczos(N)', where N is an integer, along with "
105 "galsim.Cubic, galsim.Delta, galsim.Linear, galsim.Nearest, "
106 "galsim.Quintic, and galsim.SincInterpolant.",
107 check=_validateGalsimInterpolant,
108 default="Lanczos(11)",
109 )
110 debugStarData = pexConfig.Field[bool](
111 doc="Include star images used for fitting in PSF model object.",
112 default=False
113 )
114 useCoordinates = pexConfig.ChoiceField[str](
115 doc="Which spatial coordinates to regress against in PSF modeling.",
116 allowed=dict(
117 pixel='Regress against pixel coordinates',
118 field='Regress against field angles',
119 sky='Regress against RA/Dec'
120 ),
121 default='pixel'
122 )
123
124 def setDefaults(self):
125 super().setDefaults()
126 # stampSize should be at least 25 so that
127 # i) aperture flux with 12 pixel radius can be compared to PSF flux.
128 # ii) fake sources injected to match the 12 pixel aperture flux get
129 # measured correctly
130 self.stampSize = 25
131
132
133def getGoodPixels(maskedImage, zeroWeightMaskBits):
134 """Compute an index array indicating good pixels to use.
135
136 Parameters
137 ----------
138 maskedImage : `afw.image.MaskedImage`
139 PSF candidate postage stamp
140 zeroWeightMaskBits : `List[str]`
141 List of mask bits for which to set pixel weights to zero.
142
143 Returns
144 -------
145 good : `ndarray`
146 Index array indicating good pixels.
147 """
148 imArr = maskedImage.image.array
149 varArr = maskedImage.variance.array
150 bitmask = maskedImage.mask.getPlaneBitMask(zeroWeightMaskBits)
151 good = (
152 (varArr != 0)
153 & (np.isfinite(varArr))
154 & (np.isfinite(imArr))
155 & ((maskedImage.mask.array & bitmask) == 0)
156 )
157 return good
158
159
160def computeWeight(maskedImage, maxSNR, good):
161 """Derive a weight map without Poisson variance component due to signal.
162
163 Parameters
164 ----------
165 maskedImage : `afw.image.MaskedImage`
166 PSF candidate postage stamp
167 maxSNR : `float`
168 Maximum SNR applying variance floor.
169 good : `ndarray`
170 Index array indicating good pixels.
171
172 Returns
173 -------
174 weightArr : `ndarry`
175 Array to use for weight.
176
177 See Also
178 --------
179 `lsst.meas.algorithms.variance_plance.remove_signal_from_variance` :
180 Remove the Poisson contribution from sources in the variance plane of
181 an Exposure.
182 """
183 imArr = maskedImage.image.array
184 varArr = maskedImage.variance.array
185
186 # Fit a straight line to variance vs (sky-subtracted) signal.
187 # The evaluate that line at zero signal to get an estimate of the
188 # signal-free variance.
189 fit = np.polyfit(imArr[good], varArr[good], deg=1)
190 # fit is [1/gain, sky_var]
191 weightArr = np.zeros_like(imArr, dtype=float)
192 weightArr[good] = 1./fit[1]
193
194 applyMaxSNR(imArr, weightArr, good, maxSNR)
195 return weightArr
196
197
198def applyMaxSNR(imArr, weightArr, good, maxSNR):
199 """Rescale weight of bright stars to cap the computed SNR.
200
201 Parameters
202 ----------
203 imArr : `ndarray`
204 Signal (image) array of stamp.
205 weightArr : `ndarray`
206 Weight map array. May be rescaled in place.
207 good : `ndarray`
208 Index array of pixels to use when computing SNR.
209 maxSNR : `float`
210 Threshold for adjusting variance plane implementing maximum SNR.
211 """
212 # We define the SNR value following Piff. Here's the comment from that
213 # code base explaining the calculation.
214 #
215 # The S/N value that we use will be the weighted total flux where the
216 # weight function is the star's profile itself. This is the maximum S/N
217 # value that any flux measurement can possibly produce, which will be
218 # closer to an in-practice S/N than using all the pixels equally.
219 #
220 # F = Sum_i w_i I_i^2
221 # var(F) = Sum_i w_i^2 I_i^2 var(I_i)
222 # = Sum_i w_i I_i^2 <--- Assumes var(I_i) = 1/w_i
223 #
224 # S/N = F / sqrt(var(F))
225 #
226 # Note that if the image is pure noise, this will produce a "signal" of
227 #
228 # F_noise = Sum_i w_i 1/w_i = Npix
229 #
230 # So for a more accurate estimate of the S/N of the actual star itself, one
231 # should subtract off Npix from the measured F.
232 #
233 # The final formula then is:
234 #
235 # F = Sum_i w_i I_i^2
236 # S/N = (F-Npix) / sqrt(F)
237 F = np.sum(weightArr[good]*imArr[good]**2, dtype=float)
238 Npix = np.sum(good)
239 SNR = 0.0 if F < Npix else (F-Npix)/np.sqrt(F)
240 # rescale weight of bright stars. Essentially makes an error floor.
241 if SNR > maxSNR:
242 factor = (maxSNR / SNR)**2
243 weightArr[good] *= factor
244
245
246def _computeWeightAlternative(maskedImage, maxSNR):
247 """Alternative algorithm for creating weight map.
248
249 This version is equivalent to that used by Piff internally. The weight map
250 it produces tends to leave a residual when removing the Poisson component
251 due to the signal. We leave it here as a reference, but without intending
252 that it be used (or be maintained).
253 """
254 imArr = maskedImage.image.array
255 varArr = maskedImage.variance.array
256 good = (varArr != 0) & np.isfinite(varArr) & np.isfinite(imArr)
257
258 fit = np.polyfit(imArr[good], varArr[good], deg=1)
259 # fit is [1/gain, sky_var]
260 gain = 1./fit[0]
261 varArr[good] -= imArr[good] / gain
262 weightArr = np.zeros_like(imArr, dtype=float)
263 weightArr[good] = 1./varArr[good]
264
265 applyMaxSNR(imArr, weightArr, good, maxSNR)
266 return weightArr
267
268
270 """A measurePsfTask PSF estimator using Piff as the implementation.
271 """
272 ConfigClass = PiffPsfDeterminerConfig
273 _DefaultName = "psfDeterminer.Piff"
274
276 self, exposure, psfCandidateList, metadata=None, flagKey=None
277 ):
278 """Determine a Piff PSF model for an exposure given a list of PSF
279 candidates.
280
281 Parameters
282 ----------
283 exposure : `lsst.afw.image.Exposure`
284 Exposure containing the PSF candidates.
285 psfCandidateList : `list` of `lsst.meas.algorithms.PsfCandidate`
286 A sequence of PSF candidates typically obtained by detecting sources
287 and then running them through a star selector.
288 metadata : `lsst.daf.base import PropertyList` or `None`, optional
289 A home for interesting tidbits of information.
290 flagKey : `str` or `None`, optional
291 Schema key used to mark sources actually used in PSF determination.
292
293 Returns
294 -------
295 psf : `lsst.meas.extensions.piff.PiffPsf`
296 The measured PSF model.
297 psfCellSet : `None`
298 Unused by this PsfDeterminer.
299 """
300 psfCandidateList = self.downsampleCandidates(psfCandidateList)
301
302 if self.config.stampSize:
303 stampSize = self.config.stampSize
304 if stampSize > psfCandidateList[0].getWidth():
305 self.log.warning("stampSize is larger than the PSF candidate size. Using candidate size.")
306 stampSize = psfCandidateList[0].getWidth()
307 else: # TODO: Only the if block should stay after DM-36311
308 self.log.debug("stampSize not set. Using candidate size.")
309 stampSize = psfCandidateList[0].getWidth()
310
311 scale = exposure.getWcs().getPixelScale().asArcseconds()
312 match self.config.useCoordinates:
313 case 'field':
314 detector = exposure.getDetector()
315 pix_to_field = detector.getTransform(PIXELS, FIELD_ANGLE)
316 gswcs = UVWcsWrapper(pix_to_field)
317 pointing = None
318 case 'sky':
319 gswcs = CelestialWcsWrapper(exposure.getWcs())
320 skyOrigin = exposure.getWcs().getSkyOrigin()
321 ra = skyOrigin.getLongitude().asDegrees()
322 dec = skyOrigin.getLatitude().asDegrees()
323 pointing = galsim.CelestialCoord(
324 ra*galsim.degrees,
325 dec*galsim.degrees
326 )
327 case 'pixel':
328 gswcs = galsim.PixelScale(scale)
329 pointing = None
330
331 stars = []
332 for candidate in psfCandidateList:
333 cmi = candidate.getMaskedImage(stampSize, stampSize)
334 good = getGoodPixels(cmi, self.config.zeroWeightMaskBits)
335 fracGood = np.sum(good)/good.size
336 if fracGood < self.config.minimumUnmaskedFraction:
337 continue
338 weight = computeWeight(cmi, self.config.maxSNR, good)
339
340 bbox = cmi.getBBox()
341 bds = galsim.BoundsI(
342 galsim.PositionI(*bbox.getMin()),
343 galsim.PositionI(*bbox.getMax())
344 )
345 gsImage = galsim.Image(bds, wcs=gswcs, dtype=float)
346 gsImage.array[:] = cmi.image.array
347 gsWeight = galsim.Image(bds, wcs=gswcs, dtype=float)
348 gsWeight.array[:] = weight
349
350 source = candidate.getSource()
351 image_pos = galsim.PositionD(source.getX(), source.getY())
352
353 data = piff.StarData(
354 gsImage,
355 image_pos,
356 weight=gsWeight,
357 pointing=pointing
358 )
359 stars.append(piff.Star(data, None))
360
361 piffConfig = {
362 'type': "Simple",
363 'model': {
364 'type': 'PixelGrid',
365 'scale': scale * self.config.samplingSize,
366 'size': stampSize,
367 'interp': self.config.interpolant
368 },
369 'interp': {
370 'type': 'BasisPolynomial',
371 'order': self.config.spatialOrder
372 },
373 'outliers': {
374 'type': 'Chisq',
375 'nsigma': self.config.outlierNSigma,
376 'max_remove': self.config.outlierMaxRemove
377 }
378 }
379
380 piffResult = piff.PSF.process(piffConfig)
381 wcs = {0: gswcs}
382
383 piffResult.fit(stars, wcs, pointing, logger=self.log)
384 drawSize = 2*np.floor(0.5*stampSize/self.config.samplingSize) + 1
385
386 used_image_pos = [s.image_pos for s in piffResult.stars]
387 if flagKey:
388 for candidate in psfCandidateList:
389 source = candidate.getSource()
390 posd = galsim.PositionD(source.getX(), source.getY())
391 if posd in used_image_pos:
392 source.set(flagKey, True)
393
394 if metadata is not None:
395 metadata["spatialFitChi2"] = piffResult.chisq
396 metadata["numAvailStars"] = len(stars)
397 metadata["numGoodStars"] = len(piffResult.stars)
398 metadata["avgX"] = np.mean([p.x for p in piffResult.stars])
399 metadata["avgY"] = np.mean([p.y for p in piffResult.stars])
400
401 if not self.config.debugStarData:
402 for star in piffResult.stars:
403 # Remove large data objects from the stars
404 del star.fit.params
405 del star.fit.params_var
406 del star.fit.A
407 del star.fit.b
408 del star.data.image
409 del star.data.weight
410 del star.data.orig_weight
411
412 return PiffPsf(drawSize, drawSize, piffResult), None
413
414
415measAlg.psfDeterminerRegistry.register("piff", PiffPsfDeterminerTask)
determinePsf(self, exposure, psfCandidateList, metadata=None, flagKey=None)
getGoodPixels(maskedImage, zeroWeightMaskBits)
applyMaxSNR(imArr, weightArr, good, maxSNR)