22__all__ = [
"PiffPsfDeterminerConfig",
"PiffPsfDeterminerTask"]
32from .piffPsf
import PiffPsf
35def _validateGalsimInterpolant(name: str) -> bool:
36 """A helper function to validate the GalSim interpolant at config time.
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
48 galsim.SincInterpolant
53 Whether the provided interpolant name
is valid.
56 for pattern
in (re.compile(
r"Lanczos\(\d+\)"), re.compile(
r"galsim.Lanczos\(\d+\)"),):
57 match = re.match(pattern, name)
60 return match.end() == len(name)
63 names = {f
"galsim.{interp}" for interp
in
64 (
"Cubic",
"Delta",
"Linear",
"Nearest",
"Quintic",
"SincInterpolant")
70 spatialOrder = pexConfig.Field[int](
71 doc=
"specify spatial order for PSF kernel creation",
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",
79 outlierNSigma = pexConfig.Field[float](
80 doc=
"n sigma for chisq outlier rejection",
83 outlierMaxRemove = pexConfig.Field[float](
84 doc=
"Max fraction of stars to remove as outliers each iteration",
87 maxSNR = pexConfig.Field[float](
88 doc=
"Rescale the weight of bright stars such that their SNR is less "
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']
96 minimumUnmaskedFraction = pexConfig.Field[float](
97 doc=
"Minimum fraction of unmasked pixels required to use star.",
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)",
108 debugStarData = pexConfig.Field[bool](
109 doc=
"Include star images used for fitting in PSF model object.",
123 """Compute an index array indicating good pixels to use.
128 PSF candidate postage stamp
129 zeroWeightMaskBits : `List[str]`
130 List of mask bits for which to set pixel weights to zero.
135 Index array indicating good pixels.
137 imArr = maskedImage.image.array
138 varArr = maskedImage.variance.array
139 bitmask = maskedImage.mask.getPlaneBitMask(zeroWeightMaskBits)
142 & (np.isfinite(varArr))
143 & (np.isfinite(imArr))
144 & ((maskedImage.mask.array & bitmask) == 0)
150 """Derive a weight map without Poisson variance component due to signal.
155 PSF candidate postage stamp
157 Maximum SNR applying variance floor.
159 Index array indicating good pixels.
164 Array to use for weight.
166 imArr = maskedImage.image.array
167 varArr = maskedImage.variance.array
172 fit = np.polyfit(imArr[good], varArr[good], deg=1)
174 weightArr = np.zeros_like(imArr, dtype=float)
175 weightArr[good] = 1./fit[1]
182 """Rescale weight of bright stars to cap the computed SNR.
187 Signal (image) array of stamp.
188 weightArr : `ndarray`
189 Weight map array. May be rescaled in place.
191 Index array of pixels to use when computing SNR.
193 Threshold
for adjusting variance plane implementing maximum SNR.
220 F = np.sum(weightArr[good]*imArr[good]**2, dtype=float)
222 SNR = 0.0
if F < Npix
else (F-Npix)/np.sqrt(F)
225 factor = (maxSNR / SNR)**2
226 weightArr[good] *= factor
229def _computeWeightAlternative(maskedImage, maxSNR):
230 """Alternative algorithm for creating weight map.
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).
237 imArr = maskedImage.image.array
238 varArr = maskedImage.variance.array
239 good = (varArr != 0) & np.isfinite(varArr) & np.isfinite(imArr)
241 fit = np.polyfit(imArr[good], varArr[good], deg=1)
244 varArr[good] -= imArr[good] / gain
245 weightArr = np.zeros_like(imArr, dtype=float)
246 weightArr[good] = 1./varArr[good]
253 """A measurePsfTask PSF estimator using Piff as the implementation.
255 ConfigClass = PiffPsfDeterminerConfig
256 _DefaultName = "psfDeterminer.Piff"
259 self, exposure, psfCandidateList, metadata=None, flagKey=None
261 """Determine a Piff PSF model for an exposure given a list of PSF
267 Exposure containing the PSF candidates.
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.
278 psf : `lsst.meas.extensions.piff.PiffPsf`
279 The measured PSF model.
281 Unused by this PsfDeterminer.
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()
289 self.log.debug(
"stampSize not set. Using candidate size.")
290 stampSize = psfCandidateList[0].getWidth()
295 for candidate
in psfCandidateList:
296 cmi = candidate.getMaskedImage(stampSize, stampSize)
298 fracGood = np.sum(good)/good.size
299 if fracGood < self.config.minimumUnmaskedFraction:
304 bds = galsim.BoundsI(
305 galsim.PositionI(*bbox.getMin()),
306 galsim.PositionI(*bbox.getMax())
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
313 source = candidate.getSource()
314 image_pos = galsim.PositionD(source.getX(), source.getY())
316 data = piff.StarData(
321 stars.append(piff.Star(data,
None))
327 'scale': self.config.samplingSize,
329 'interp': self.config.interpolant
332 'type':
'BasisPolynomial',
333 'order': self.config.spatialOrder
337 'nsigma': self.config.outlierNSigma,
338 'max_remove': self.config.outlierMaxRemove
342 piffResult = piff.PSF.process(piffConfig)
344 wcs = {0: galsim.PixelScale(1.0)}
347 piffResult.fit(stars, wcs, pointing, logger=self.log)
348 drawSize = 2*np.floor(0.5*stampSize/self.config.samplingSize) + 1
350 used_image_pos = [s.image_pos
for s
in piffResult.stars]
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)
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])
365 if not self.config.debugStarData:
366 for star
in piffResult.stars:
369 del star.fit.params_var
374 del star.data.orig_weight
376 return PiffPsf(drawSize, drawSize, piffResult),
None
380 def _validatePsfCandidates(psfCandidateList, stampSize):
381 """Raise if psfCandidates are smaller than the configured kernelSize.
386 Sequence of psf candidates to check.
388 Size of image model to use in PIFF.
393 Raised
if any psfCandidate has width
or height smaller than
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()}."
405measAlg.psfDeterminerRegistry.register(
"piff", PiffPsfDeterminerTask)
A class to contain the data, WCS, and other information needed to describe an image of the sky.
A class to manipulate images, masks, and variance as a single object.
Class stored in SpatialCells for spatial Psf fitting.
def _validatePsfCandidates(psfCandidateList, stampSize)
def determinePsf(self, exposure, psfCandidateList, metadata=None, flagKey=None)
def getGoodPixels(maskedImage, zeroWeightMaskBits)
def applyMaxSNR(imArr, weightArr, good, maxSNR)
def computeWeight(maskedImage, maxSNR, good)