Coverage for python/lsst/meas/extensions/piff/piffPsfDeterminer.py: 84%
245 statements
« prev ^ index » next coverage.py v7.16.2, created at 2026-11-10 01:38 -0800
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/>.
22__all__ = ["PiffPsfDeterminerConfig", "PiffPsfDeterminerTask"]
24import numpy as np
25import piff
26import galsim
27import re
28import logging
29import yaml
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 lsst.pipe.base import AlgorithmError
37from .piffPsf import PiffPsf
38from .wcs_wrapper import CelestialWcsWrapper, UVWcsWrapper
41def _validateGalsimInterpolant(name: str) -> bool:
42 """A helper function to validate the GalSim interpolant at config time.
44 Parameters
45 ----------
46 name : str
47 The name of the interpolant to use from GalSim. Valid options are:
48 galsim.Lanczos(N) or Lancsos(N), where N is a positive integer
49 galsim.Linear
50 galsim.Cubic
51 galsim.Quintic
52 galsim.Delta
53 galsim.Nearest
54 galsim.SincInterpolant
56 Returns
57 -------
58 is_valid : bool
59 Whether the provided interpolant name is valid.
60 """
61 # First, check if ``name`` is a valid Lanczos interpolant.
62 for pattern in (re.compile(r"Lanczos\(\d+\)"), re.compile(r"galsim.Lanczos\(\d+\)"),):
63 match = re.match(pattern, name) # Search from the start of the string.
64 if match is not None:
65 # Check that the pattern is also the end of the string.
66 return match.end() == len(name)
68 # If not, check if ``name`` is any other valid GalSim interpolant.
69 names = {f"galsim.{interp}" for interp in
70 ("Cubic", "Delta", "Linear", "Nearest", "Quintic", "SincInterpolant")
71 }
72 return name in names
75class PiffTooFewGoodStarsError(AlgorithmError):
76 """Raised if too few good stars are available for PSF determination.
78 Parameters
79 ----------
80 num_good_stars : `int`
81 Number of good stars used for PSF determination.
82 poly_ndim : `int`
83 Number of dependency parameters (dimensions) used in
84 polynomial fitting.
85 minimum_dof : `int`
86 Minimum number of degree of freedom to do the fit.
87 """
89 def __init__(
90 self,
91 num_good_stars,
92 minimum_dof,
93 poly_ndim,
94 ):
95 self._num_good_stars = num_good_stars
96 self._poly_ndim = poly_ndim
97 self._minimum_dof = minimum_dof
98 super().__init__(
99 f"Failed to determine piff psf: too few good stars ({num_good_stars}) and minimum dof to fit "
100 f"a {poly_ndim} order polynomial is {minimum_dof}."
101 )
103 @property
104 def metadata(self):
105 return {
106 "num_good_stars": self._num_good_stars,
107 "poly_ndim": self._poly_ndim,
108 "minimum_dof": self._minimum_dof,
109 }
112class PiffPsfDeterminerConfig(BasePsfDeterminerTask.ConfigClass):
113 spatialOrderPerBand = pexConfig.DictField(
114 doc="Per-band spatial order for PSF kernel creation. "
115 "Ignored if piffPsfConfigYaml is set.",
116 keytype=str,
117 itemtype=int,
118 default={},
119 )
120 spatialOrder = pexConfig.Field[int](
121 doc="Spatial order for PSF kernel creation. "
122 "Ignored if piffPsfConfigYaml is set or if the current "
123 "band is in the spatialOrderPerBand dict config.",
124 default=2,
125 )
126 piffBasisPolynomialSolver = pexConfig.ChoiceField[str](
127 doc="Solver to solve linear algebra for "
128 "the spatial interpolation of the PSF. "
129 "Ignore if piffPsfConfigYaml is not None.",
130 allowed=dict(
131 scipy="Default solver using scipy.",
132 cpp="C++ solver using eigen.",
133 jax="JAX solver.",
134 qr="QR decomposition using scipy/numpy.",
135 ),
136 default="scipy",
137 )
138 piffPixelGridFitCenter = pexConfig.Field[bool](
139 doc="PixelGrid can re-estimate center if this option is True. "
140 "Will use provided centroid if set to False. Default in Piff "
141 "is set to True. Depends on how input centroids can be trust. "
142 "Ignore if piffPsfConfigYaml is not None.",
143 default=True
144 )
145 samplingSize = pexConfig.Field[float](
146 doc="Resolution of the internal PSF model relative to the pixel size; "
147 "e.g. 0.5 is equal to 2x oversampling. This affects only the size of "
148 "the PSF model stamp if piffPsfConfigYaml is set.",
149 default=1,
150 )
151 modelSize = pexConfig.Field[int](
152 doc="Internal model size for PIFF (typically odd, but not enforced). "
153 "Partially ignored if piffPsfConfigYaml is set.",
154 default=25,
155 )
156 outlierNSigma = pexConfig.Field[float](
157 doc="n sigma for chisq outlier rejection. "
158 "Ignored if piffPsfConfigYaml is set.",
159 default=4.0
160 )
161 outlierMaxRemove = pexConfig.Field[float](
162 doc="Max fraction of stars to remove as outliers each iteration. "
163 "Ignored if piffPsfConfigYaml is set.",
164 default=0.05
165 )
166 maxSNR = pexConfig.Field[float](
167 doc="Rescale the weight of bright stars such that their SNR is less "
168 "than this value.",
169 default=200.0
170 )
171 zeroWeightMaskBits = pexConfig.ListField[str](
172 doc="List of mask bits for which to set pixel weights to zero.",
173 default=['BAD', 'CR', 'INTRP', 'SAT', 'SUSPECT', 'NO_DATA']
174 )
175 minimumUnmaskedFraction = pexConfig.Field[float](
176 doc="Minimum fraction of unmasked pixels required to use star.",
177 default=0.5
178 )
179 useColor = pexConfig.Field[bool](
180 doc="Use color information to correct for amtospheric turbulences and "
181 "differential chromatic refraction."
182 "Ignored if piffPsfConfigYaml is set.",
183 default=False,
184 )
185 color = pexConfig.DictField(
186 doc="The bands to use for calculating color."
187 "Ignored if piffPsfConfigYaml is set.",
188 default={},
189 keytype=str,
190 itemtype=str,
191 )
192 colorOrder = pexConfig.Field[int](
193 doc="Color order for PSF kernel creation. Only pure color terms are "
194 "fitted; cross terms between color and spatial coordinates are "
195 "not included. Ignored if piffPsfConfigYaml is set.",
196 default=1,
197 )
198 interpolant = pexConfig.Field[str](
199 doc="GalSim interpolant name for Piff to use. "
200 "Options include 'Lanczos(N)', where N is an integer, along with "
201 "galsim.Cubic, galsim.Delta, galsim.Linear, galsim.Nearest, "
202 "galsim.Quintic, and galsim.SincInterpolant. Ignored if "
203 "piffPsfConfigYaml is set.",
204 check=_validateGalsimInterpolant,
205 default="Lanczos(11)",
206 )
207 zerothOrderInterpNotEnoughStars = pexConfig.Field[bool](
208 doc="If True, use zeroth order interpolation if not enough stars are found.",
209 default=False
210 )
211 debugStarData = pexConfig.Field[bool](
212 doc="Include star images used for fitting in PSF model object.",
213 default=False
214 )
215 useCoordinates = pexConfig.ChoiceField[str](
216 doc="Which spatial coordinates to regress against in PSF modeling.",
217 allowed=dict(
218 pixel='Regress against pixel coordinates',
219 field='Regress against field angles',
220 sky='Regress against RA/Dec'
221 ),
222 default='pixel'
223 )
224 piffMaxIter = pexConfig.Field[int](
225 doc="Maximum iteration while doing outlier rejection."
226 "Might be overwrite if zerothOrderInterpNotEnoughStars is True and "
227 "ignore if piffPsfConfigYaml is not None.",
228 default=15
229 )
230 piffLoggingLevel = pexConfig.RangeField[int](
231 doc="PIFF-specific logging level, from 0 (least chatty) to 3 (very verbose); "
232 "note that logs will come out with unrelated notations (e.g. WARNING does "
233 "not imply a warning.)",
234 default=0,
235 min=0,
236 max=3,
237 )
238 piffPsfConfigYaml = pexConfig.Field[str](
239 doc="Configuration file for PIFF in YAML format. This overrides many "
240 "other settings in this config and is not validated ahead of runtime.",
241 default=None,
242 optional=True,
243 )
245 def setDefaults(self):
246 super().setDefaults()
247 # stampSize should be at least 25 so that
248 # i) aperture flux with 12 pixel radius can be compared to PSF flux.
249 # ii) fake sources injected to match the 12 pixel aperture flux get
250 # measured correctly
251 # stampSize should also be at least sqrt(2)*modelSize/samplingSize so
252 # that rotated images have support in the model
254 self.stampSize = 27
255 # Resize the stamp to accommodate the model, but only if necessary.
256 if self.useCoordinates == "sky": 256 ↛ 257line 256 didn't jump to line 257 because the condition on line 256 was never true
257 self.stampSize = max(25, 2*int(0.5*self.modelSize*np.sqrt(2)/self.samplingSize) + 1)
259 def validate(self):
260 super().validate()
262 if (stamp_size := self.stampSize) is not None:
263 model_size = self.modelSize
264 sampling_size = self.samplingSize
265 if self.useCoordinates == 'sky':
266 min_stamp_size = int(np.sqrt(2) * model_size / sampling_size)
267 else:
268 min_stamp_size = int(model_size / sampling_size)
270 if stamp_size < min_stamp_size:
271 msg = (f"PIFF model size of {model_size} is too large for stamp size {stamp_size}. "
272 f"Set stampSize >= {min_stamp_size}"
273 )
274 raise pexConfig.FieldValidationError(self.__class__.modelSize, self, msg)
277def getGoodPixels(maskedImage, zeroWeightMaskBits):
278 """Compute an index array indicating good pixels to use.
280 Parameters
281 ----------
282 maskedImage : `afw.image.MaskedImage`
283 PSF candidate postage stamp
284 zeroWeightMaskBits : `List[str]`
285 List of mask bits for which to set pixel weights to zero.
287 Returns
288 -------
289 good : `ndarray`
290 Index array indicating good pixels.
291 """
292 imArr = maskedImage.image.array
293 varArr = maskedImage.variance.array
294 bitmask = maskedImage.mask.getPlaneBitMask(zeroWeightMaskBits)
295 good = (
296 (varArr != 0)
297 & (np.isfinite(varArr))
298 & (np.isfinite(imArr))
299 & ((maskedImage.mask.array & bitmask) == 0)
300 )
301 return good
304def computeWeight(maskedImage, maxSNR, good):
305 """Derive a weight map without Poisson variance component due to signal.
307 Parameters
308 ----------
309 maskedImage : `afw.image.MaskedImage`
310 PSF candidate postage stamp
311 maxSNR : `float`
312 Maximum SNR applying variance floor.
313 good : `ndarray`
314 Index array indicating good pixels.
316 Returns
317 -------
318 weightArr : `ndarry`
319 Array to use for weight.
321 See Also
322 --------
323 `lsst.meas.algorithms.variance_plance.remove_signal_from_variance` :
324 Remove the Poisson contribution from sources in the variance plane of
325 an Exposure.
326 """
327 imArr = maskedImage.image.array
328 varArr = maskedImage.variance.array
330 # Fit a straight line to variance vs (sky-subtracted) signal.
331 # The evaluate that line at zero signal to get an estimate of the
332 # signal-free variance.
333 fit = np.polyfit(imArr[good], varArr[good], deg=1)
334 # fit is [1/gain, sky_var]
335 weightArr = np.zeros_like(imArr, dtype=float)
336 weightArr[good] = 1./fit[1]
338 applyMaxSNR(imArr, weightArr, good, maxSNR)
339 return weightArr
342def applyMaxSNR(imArr, weightArr, good, maxSNR):
343 """Rescale weight of bright stars to cap the computed SNR.
345 Parameters
346 ----------
347 imArr : `ndarray`
348 Signal (image) array of stamp.
349 weightArr : `ndarray`
350 Weight map array. May be rescaled in place.
351 good : `ndarray`
352 Index array of pixels to use when computing SNR.
353 maxSNR : `float`
354 Threshold for adjusting variance plane implementing maximum SNR.
355 """
356 # We define the SNR value following Piff. Here's the comment from that
357 # code base explaining the calculation.
358 #
359 # The S/N value that we use will be the weighted total flux where the
360 # weight function is the star's profile itself. This is the maximum S/N
361 # value that any flux measurement can possibly produce, which will be
362 # closer to an in-practice S/N than using all the pixels equally.
363 #
364 # F = Sum_i w_i I_i^2
365 # var(F) = Sum_i w_i^2 I_i^2 var(I_i)
366 # = Sum_i w_i I_i^2 <--- Assumes var(I_i) = 1/w_i
367 #
368 # S/N = F / sqrt(var(F))
369 #
370 # Note that if the image is pure noise, this will produce a "signal" of
371 #
372 # F_noise = Sum_i w_i 1/w_i = Npix
373 #
374 # So for a more accurate estimate of the S/N of the actual star itself, one
375 # should subtract off Npix from the measured F.
376 #
377 # The final formula then is:
378 #
379 # F = Sum_i w_i I_i^2
380 # S/N = (F-Npix) / sqrt(F)
381 F = np.sum(weightArr[good]*imArr[good]**2, dtype=float)
382 Npix = np.sum(good)
383 SNR = 0.0 if F < Npix else (F-Npix)/np.sqrt(F)
384 # rescale weight of bright stars. Essentially makes an error floor.
385 if SNR > maxSNR: 385 ↛ exitline 385 didn't return from function 'applyMaxSNR' because the condition on line 385 was always true
386 factor = (maxSNR / SNR)**2
387 weightArr[good] *= factor
390def _computeWeightAlternative(maskedImage, maxSNR):
391 """Alternative algorithm for creating weight map.
393 This version is equivalent to that used by Piff internally. The weight map
394 it produces tends to leave a residual when removing the Poisson component
395 due to the signal. We leave it here as a reference, but without intending
396 that it be used (or be maintained).
397 """
398 imArr = maskedImage.image.array
399 varArr = maskedImage.variance.array
400 good = (varArr != 0) & np.isfinite(varArr) & np.isfinite(imArr)
402 fit = np.polyfit(imArr[good], varArr[good], deg=1)
403 # fit is [1/gain, sky_var]
404 gain = 1./fit[0]
405 varArr[good] -= imArr[good] / gain
406 weightArr = np.zeros_like(imArr, dtype=float)
407 weightArr[good] = 1./varArr[good]
409 applyMaxSNR(imArr, weightArr, good, maxSNR)
410 return weightArr
413class PiffPsfDeterminerTask(BasePsfDeterminerTask):
414 """A measurePsfTask PSF estimator using Piff as the implementation.
415 """
416 ConfigClass = PiffPsfDeterminerConfig
417 _DefaultName = "psfDeterminer.Piff"
419 def __init__(self, config, schema=None, **kwds):
420 BasePsfDeterminerTask.__init__(self, config, schema=schema, **kwds)
422 piffLoggingLevels = {
423 0: logging.CRITICAL,
424 1: logging.WARNING,
425 2: logging.INFO,
426 3: logging.DEBUG,
427 }
428 self.piffLogger = lsst.utils.logging.getLogger(f"{self.log.name}.piff")
429 self.piffLogger.setLevel(piffLoggingLevels[self.config.piffLoggingLevel])
431 def determinePsf(
432 self, exposure, psfCandidateList, metadata=None, flagKey=None,
433 ):
434 """Determine a Piff PSF model for an exposure given a list of PSF
435 candidates.
437 Parameters
438 ----------
439 exposure : `lsst.afw.image.Exposure`
440 Exposure containing the PSF candidates.
441 psfCandidateList : `list` of `lsst.meas.algorithms.PsfCandidate`
442 A sequence of PSF candidates typically obtained by detecting sources
443 and then running them through a star selector.
444 metadata : `lsst.daf.base import PropertyList` or `None`, optional
445 A home for interesting tidbits of information.
446 flagKey : `str` or `None`, optional
447 Schema key used to mark sources actually used in PSF determination.
449 Returns
450 -------
451 psf : `lsst.meas.extensions.piff.PiffPsf`
452 The measured PSF model.
453 psfCellSet : `None`
454 Unused by this PsfDeterminer.
455 """
456 psfCandidateList = self.downsampleCandidates(psfCandidateList)
458 if self.config.stampSize:
459 stampSize = self.config.stampSize
460 if stampSize > psfCandidateList[0].getWidth(): 460 ↛ 461line 460 didn't jump to line 461 because the condition on line 460 was never true
461 self.log.warning("stampSize is larger than the PSF candidate size. Using candidate size.")
462 stampSize = psfCandidateList[0].getWidth()
463 else: # TODO: Only the if block should stay after DM-36311
464 self.log.debug("stampSize not set. Using candidate size.")
465 stampSize = psfCandidateList[0].getWidth()
467 scale = exposure.getWcs().getPixelScale(exposure.getBBox().getCenter()).asArcseconds()
469 match self.config.useCoordinates:
470 case 'field': 470 ↛ 471line 470 didn't jump to line 471 because the pattern on line 470 never matched
471 detector = exposure.getDetector()
472 pix_to_field = detector.getTransform(PIXELS, FIELD_ANGLE)
473 gswcs = UVWcsWrapper(pix_to_field)
474 pointing = None
475 keys = ['u', 'v']
477 case 'sky':
478 gswcs = CelestialWcsWrapper(exposure.getWcs())
479 skyOrigin = exposure.getWcs().getSkyOrigin()
480 ra = skyOrigin.getLongitude().asDegrees()
481 dec = skyOrigin.getLatitude().asDegrees()
482 pointing = galsim.CelestialCoord(
483 ra*galsim.degrees,
484 dec*galsim.degrees
485 )
486 keys = ['u', 'v']
488 case 'pixel': 488 ↛ 493line 488 didn't jump to line 493 because the pattern on line 488 always matched
489 gswcs = galsim.PixelScale(scale)
490 pointing = None
491 keys = ['x', 'y']
493 stars = []
494 for candidate in psfCandidateList:
495 cmi = candidate.getMaskedImage(stampSize, stampSize)
496 good = getGoodPixels(cmi, self.config.zeroWeightMaskBits)
497 fracGood = np.sum(good)/good.size
498 if fracGood < self.config.minimumUnmaskedFraction: 498 ↛ 499line 498 didn't jump to line 499 because the condition on line 498 was never true
499 continue
500 weight = computeWeight(cmi, self.config.maxSNR, good)
502 bbox = cmi.getBBox()
503 bds = galsim.BoundsI(
504 galsim.PositionI(*bbox.getMin()),
505 galsim.PositionI(*bbox.getMax())
506 )
507 gsImage = galsim.Image(bds, wcs=gswcs, dtype=float)
508 gsImage.array[:] = cmi.image.array
509 gsWeight = galsim.Image(bds, wcs=gswcs, dtype=float)
510 gsWeight.array[:] = weight
512 source = candidate.getSource()
513 starId = source.getId()
514 image_pos = galsim.PositionD(source.getX(), source.getY())
516 data = piff.StarData(
517 gsImage,
518 image_pos,
519 weight=gsWeight,
520 pointing=pointing
521 )
522 star = piff.Star(data, None)
523 star.data.properties['starId'] = starId
524 star.data.properties['colorValue'] = candidate.getPsfColorValue()
525 star.data.properties['colorType'] = candidate.getPsfColorType()
526 stars.append(star)
528 # The following is mostly accommodating unittests that don't have
529 # the filter attribute set on the mock exposure.
530 if hasattr(exposure.filter, "bandLabel"): 530 ↛ 531line 530 didn't jump to line 531 because the condition on line 530 was never true
531 band = exposure.filter.bandLabel
532 else:
533 band = None
534 spatialOrder = self.config.spatialOrderPerBand.get(band, self.config.spatialOrder)
535 orders = [spatialOrder] * len(keys)
537 if self.config.useColor:
538 colors = [s.data.properties['colorValue'] for s in stars
539 if np.isfinite(s.data.properties['colorValue'])]
540 colorTypes = [s.data.properties['colorType'] for s in stars
541 if np.isfinite(s.data.properties['colorValue'])]
542 if len(colors) == 0: 542 ↛ 543line 542 didn't jump to line 543 because the condition on line 542 was never true
543 self.log.warning("No color informations for PSF candidates, Set PSF colors to 0s.")
544 meanColors = 0.
545 colorTypes = [stars[0].data.properties['colorType']]
546 else:
547 meanColors = np.mean(colors)
548 colorType = list(set(colorTypes))
549 if len(colorType) > 1: 549 ↛ 550line 549 didn't jump to line 550 because the condition on line 549 was never true
550 raise ValueError(f"More than one colorType was provided:{colorType}")
551 colorType = colorType[0]
552 for s in stars:
553 if not np.isfinite(s.data.properties['colorValue']): 553 ↛ 554line 553 didn't jump to line 554 because the condition on line 553 was never true
554 s.data.properties['colorValue'] = meanColors
555 s.data.properties['colorType'] = colorType
556 keys.append('colorValue')
557 orders.append(self.config.colorOrder)
559 if self.config.piffPsfConfigYaml is None:
560 piffConfig = {
561 'type': 'Simple',
562 'model': {
563 'type': 'PixelGrid',
564 'scale': scale * self.config.samplingSize,
565 'size': self.config.modelSize,
566 'interp': self.config.interpolant,
567 'centered': self.config.piffPixelGridFitCenter,
568 },
569 'interp': {
570 'type': 'BasisPolynomial',
571 'order': orders,
572 'keys': keys,
573 'solver': self.config.piffBasisPolynomialSolver,
574 },
575 'outliers': {
576 'type': 'Chisq',
577 'nsigma': self.config.outlierNSigma,
578 'max_remove': self.config.outlierMaxRemove,
579 },
580 'max_iter': self.config.piffMaxIter
581 }
582 if self.config.useColor:
583 # Cross terms between the two spatial coordinates are kept up
584 # to spatialOrder. Cross terms between color and the spatial
585 # coordinates are set to zero, i.e. only pure color terms up
586 # to colorOrder are fitted, so the basis is e.g.
587 # [1, c, y, y^2, x, x*y, x^2] for spatialOrder=2, colorOrder=1.
588 #
589 # In Piff's BasisPolynomial, a pair of keys omitted from the
590 # max_order dict gets no cross terms. The color/spatial pairs
591 # must be omitted rather than listed with an explicit 0:
592 # Piff <= 1.7 indexes dict pairs within the pair tuple instead
593 # of within `keys`, so an explicit ('x', 'colorValue'): 0
594 # entry would also remove the x*y term.
595 spatialKeys = (keys[0], keys[1])
596 piffConfig['interp']['max_order'] = {spatialKeys: spatialOrder}
597 else:
598 piffConfig = yaml.safe_load(self.config.piffPsfConfigYaml)
600 def _get_threshold(nth_order):
601 if isinstance(nth_order, list):
602 # right now, nth_order[0] and nth_order[1] are the same.
603 freeParameters = ((nth_order[0] + 1) * (nth_order[0] + 2)) // 2
604 if len(nth_order) == 3: # when color correction
605 freeParameters += nth_order[2]
606 else:
607 # number of free parameter in the polynomial interpolation
608 freeParameters = ((nth_order + 1) * (nth_order + 2)) // 2
609 return freeParameters
611 if piffConfig['interp']['type'] in ['BasisPolynomial', 'Polynomial']: 611 ↛ 629line 611 didn't jump to line 629 because the condition on line 611 was always true
612 threshold = _get_threshold(piffConfig['interp']['order'])
613 if len(stars) < threshold:
614 if self.config.zerothOrderInterpNotEnoughStars:
615 self.log.warning(
616 "Only %d stars found, "
617 "but %d required. Using zeroth order interpolation."%((len(stars), threshold))
618 )
619 piffConfig['interp']['order'] = [0] * len(keys)
620 # No need to do any outlier rejection assume
621 # PSF to be average of few stars.
622 piffConfig['max_iter'] = 1
623 else:
624 raise PiffTooFewGoodStarsError(
625 num_good_stars=len(stars),
626 minimum_dof=threshold,
627 poly_ndim=piffConfig['interp']['order'],
628 )
629 self._piffConfig = piffConfig
630 piffResult = piff.PSF.process(piffConfig)
631 wcs = {0: gswcs}
633 piffResult.fit(stars, wcs, pointing, logger=self.piffLogger)
635 nUsedStars = len([s for s in piffResult.stars if not s.is_flagged and not s.is_reserve])
637 if piffConfig['interp']['type'] in ['BasisPolynomial', 'Polynomial']: 637 ↛ 658line 637 didn't jump to line 658 because the condition on line 637 was always true
638 threshold = _get_threshold(piffConfig['interp']['order'])
639 if nUsedStars < threshold: 639 ↛ 640line 639 didn't jump to line 640 because the condition on line 639 was never true
640 if self.config.zerothOrderInterpNotEnoughStars:
641 self.log.warning(
642 "Only %d after outlier rejection, "
643 "but %d required. Using zeroth order interpolation."%((nUsedStars, threshold))
644 )
645 piffConfig['interp']['order'] = [0] * len(keys)
646 # No need to do any outlier rejection assume
647 # PSF to be average of few stars.
648 piffConfig['max_iter'] = 1
649 piffResult.fit(stars, wcs, pointing, logger=self.piffLogger)
650 nUsedStars = len(stars)
651 else:
652 raise PiffTooFewGoodStarsError(
653 num_good_stars=nUsedStars,
654 minimum_dof=threshold,
655 poly_ndim=piffConfig['interp']['order'],
656 )
658 drawSize = 2*np.floor(0.5*stampSize/self.config.samplingSize) + 1
660 used_image_starId = {s.data.properties['starId'] for s in piffResult.stars
661 if not s.is_flagged and not s.is_reserve}
663 assert len(used_image_starId) == nUsedStars, "Star IDs are not unique"
665 if flagKey: 665 ↛ 672line 665 didn't jump to line 672 because the condition on line 665 was always true
666 for candidate in psfCandidateList:
667 source = candidate.getSource()
668 starId = source.getId()
669 if starId in used_image_starId: 669 ↛ 666line 669 didn't jump to line 666 because the condition on line 669 was always true
670 source.set(flagKey, True)
672 if metadata is not None: 672 ↛ 681line 672 didn't jump to line 681 because the condition on line 672 was always true
673 metadata["spatialFitChi2"] = piffResult.chisq
674 metadata["numAvailStars"] = len(stars)
675 metadata["numGoodStars"] = nUsedStars
676 metadata["avgX"] = np.mean([s.x for s in piffResult.stars
677 if not s.is_flagged and not s.is_reserve])
678 metadata["avgY"] = np.mean([s.y for s in piffResult.stars
679 if not s.is_flagged and not s.is_reserve])
681 if not self.config.debugStarData:
682 for star in piffResult.stars:
683 # Remove large data objects from the stars
684 del star.fit.params
685 del star.fit.params_var
686 del star.fit.A
687 del star.fit.b
688 del star.data.image
689 del star.data.weight
690 del star.data.orig_weight
692 return PiffPsf(drawSize, drawSize, piffResult), None
695measAlg.psfDeterminerRegistry.register("piff", PiffPsfDeterminerTask)