Coverage for python / lsst / meas / extensions / shapeHSM / _hsm_moments.py: 25%
204 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-04-18 08:58 +0000
1# This file is part of meas_extensions_shapeHSM.
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/>.
22import logging
24import galsim
25import lsst.afw.geom as afwGeom
26import lsst.afw.image as afwImage
27import lsst.afw.math as afwMath
28import lsst.afw.table as afwTable
29import lsst.meas.base as measBase
30import lsst.pex.config as pexConfig
31import numpy as np
32from lsst.geom import Box2I, Point2D, Point2I
33from lsst.meas.base import colorExtractor
34from lsst.pex.exceptions import InvalidParameterError, NotFoundError
36__all__ = [
37 "HsmSourceMomentsConfig",
38 "HsmSourceMomentsPlugin",
39 "HsmSourceMomentsRoundConfig",
40 "HsmSourceMomentsRoundPlugin",
41 "HsmPsfMomentsConfig",
42 "HsmPsfMomentsPlugin",
43 "HsmPsfMomentsDebiasedConfig",
44 "HsmPsfMomentsDebiasedPlugin",
45]
48class HsmMomentsConfig(measBase.SingleFramePluginConfig):
49 """Base configuration for HSM adaptive moments measurement."""
51 roundMoments = pexConfig.Field[bool](doc="Use round weight function?", default=False)
52 addFlux = pexConfig.Field[bool](doc="Store measured flux?", default=False)
53 subtractCenter = pexConfig.Field[bool](doc="Subtract starting center from x/y outputs?", default=False)
56class HsmMomentsPlugin(measBase.SingleFramePlugin):
57 """Base plugin for HSM adaptive moments measurement."""
59 ConfigClass = HsmMomentsConfig
61 def __init__(self, config, name, schema, metadata, logName=None):
62 if logName is None:
63 logName = __name__
64 super().__init__(config, name, schema, metadata, logName=logName)
66 # Define flags for possible issues that might arise during measurement.
67 flagDefs = measBase.FlagDefinitionList()
68 self.FAILURE = flagDefs.addFailureFlag("General failure flag, set if anything went wrong")
69 self.NO_PIXELS = flagDefs.add("flag_no_pixels", "No pixels to measure")
70 self.NOT_CONTAINED = flagDefs.add(
71 "flag_not_contained", "Center not contained in footprint bounding box"
72 )
73 self.PARENT_SOURCE = flagDefs.add("flag_parent_source", "Parent source, ignored")
74 self.GALSIM = flagDefs.add("flag_galsim", "GalSim failure")
75 self.INVALID_PARAM = flagDefs.add("flag_invalid_param", "Invalid combination of moments")
76 self.EDGE = flagDefs.add("flag_edge", "Variance undefined outside image edge")
77 self.NO_PSF = flagDefs.add("flag_no_psf", "Exposure lacks PSF")
79 # Embed the flag definitions in the schema using a flag handler.
80 self.flagHandler = measBase.FlagHandler.addFields(schema, name, flagDefs)
82 # Utilize a safe centroid extractor that uses the detection footprint
83 # as a fallback if necessary.
84 self.centroidExtractor = measBase.SafeCentroidExtractor(schema, name)
85 self.log = logging.getLogger(self.logName)
87 @classmethod
88 def getExecutionOrder(cls):
89 return cls.SHAPE_ORDER
91 def _calculate(
92 self,
93 record: afwTable.SourceRecord,
94 *,
95 image: galsim.Image,
96 weight_image: galsim.Image,
97 centroid: Point2D,
98 sigma: float = 5.0,
99 precision: float = 1.0e-6,
100 ) -> None:
101 """
102 Calculate adaptive moments using GalSim's HSM and modify the record in
103 place.
105 Parameters
106 ----------
107 record : `~lsst.afw.table.SourceRecord`
108 Record to store measurements.
109 image : `~galsim.Image`
110 Image on which to perform measurements.
111 weight_image : `~galsim.Image`
112 The combined badpix/weight image for input to galsim HSM code.
113 centroid : `~lsst.geom.Point2D`
114 Centroid guess for HSM adaptive moments.
115 sigma : `float`, optional
116 Estimate of object's Gaussian sigma in pixels. Default is 5.0.
117 precision : `float`, optional
118 Precision for HSM adaptive moments. Default is 1.0e-6.
120 Raises
121 ------
122 MeasurementError
123 Raised for errors in measurement.
124 """
125 # Convert centroid to GalSim's PositionD type.
126 # Using _PositionD skips some sanity checks and argument parsing of the
127 # public PositionD class. Note that GalSim considers _PositionD as a
128 # public constructor despite the leading underscore.
129 guessCentroid = galsim._PositionD(centroid.x, centroid.y)
130 try:
131 # Attempt to compute HSM moments.
132 shape = galsim.hsm.FindAdaptiveMom(
133 image,
134 weight=weight_image,
135 badpix=None, # Already incorporated into `weight_image`.
136 guess_sig=sigma,
137 precision=precision,
138 guess_centroid=guessCentroid,
139 strict=True, # Raises GalSimHSMError if estimation fails.
140 check=False, # This speeds up the code!
141 round_moments=self.config.roundMoments,
142 hsmparams=None,
143 )
144 except galsim.GalSimHSMError as error:
145 raise measBase.MeasurementError(str(error), self.GALSIM.number)
147 # Retrieve computed moments sigma and centroid.
148 determinantRadius = shape.moments_sigma
149 centroidResult = shape.moments_centroid
151 # Subtract center if required by configuration.
152 if self.config.subtractCenter:
153 centroidResult.x -= centroid.getX()
154 centroidResult.y -= centroid.getY()
156 # Convert GalSim's `galsim.PositionD` to `lsst.geom.Point2D`.
157 centroidResult = Point2D(centroidResult.x, centroidResult.y)
159 # Populate the record with the centroid results.
160 record.set(self.centroidResultKey, centroidResult)
162 # Convert GalSim measurements to lsst measurements.
163 try:
164 # Create an ellipse for the shape.
165 observed_shape = shape.observed_shape
166 ellipse = afwGeom.ellipses.SeparableDistortionDeterminantRadius(
167 e1=observed_shape.e1,
168 e2=observed_shape.e2,
169 radius=determinantRadius,
170 normalize=True, # Fail if |e|>1.
171 )
172 # Get the quadrupole moments from the ellipse.
173 quad = afwGeom.ellipses.Quadrupole(ellipse)
174 except InvalidParameterError as error:
175 raise measBase.MeasurementError(error, self.INVALID_PARAM.number)
177 # Store the quadrupole moments in the record.
178 record.set(self.shapeKey, quad)
180 # Store the flux if required by configuration.
181 if self.config.addFlux:
182 record.set(self.fluxKey, shape.moments_amp)
184 def fail(self, record, error=None):
185 # Docstring inherited.
186 self.flagHandler.handleFailure(record)
187 if error:
188 centroid = self.centroidExtractor(record, self.flagHandler)
189 self.log.debug(
190 "Failed to measure shape for %d at (%f, %f): %s",
191 record.getId(),
192 centroid.getX(),
193 centroid.getY(),
194 error,
195 )
198class HsmSourceMomentsConfig(HsmMomentsConfig):
199 """Configuration for HSM adaptive moments measurement for sources."""
201 badMaskPlanes = pexConfig.ListField[str](
202 doc="Mask planes used to reject bad pixels.", default=["BAD", "SAT"]
203 )
206@measBase.register("ext_shapeHSM_HsmSourceMoments")
207class HsmSourceMomentsPlugin(HsmMomentsPlugin):
208 """Plugin for HSM adaptive moments measurement for sources."""
210 ConfigClass = HsmSourceMomentsConfig
212 def __init__(self, config, name, schema, metadata, logName=None):
213 super().__init__(config, name, schema, metadata, logName=logName)
214 self.centroidResultKey = afwTable.Point2DKey.addFields(
215 schema, name, "Centroid of the source via the HSM shape algorithm", "pixel"
216 )
217 self.shapeKey = afwTable.QuadrupoleKey.addFields(
218 schema,
219 name,
220 "Adaptive moments of the source via the HSM shape algorithm",
221 afwTable.CoordinateType.PIXEL,
222 )
223 if config.addFlux:
224 self.fluxKey = schema.addField(
225 schema.join(name, "Flux"), type=float, doc="Flux of the source via the HSM shape algorithm"
226 )
228 def measure(self, record, exposure):
229 """
230 Measure adaptive moments of sources given an exposure and set the
231 results in the record in place.
233 Parameters
234 ----------
235 record : `~lsst.afw.table.SourceRecord`
236 The record where measurement outputs will be stored.
237 exposure : `~lsst.afw.image.Exposure`
238 The exposure containing the source which needs measurement.
240 Raises
241 ------
242 MeasurementError
243 Raised for errors in measurement.
244 """
245 # Extract the centroid from the record.
246 center = self.centroidExtractor(record, self.flagHandler)
247 color = colorExtractor(record)
249 # Get the bounding box of the source's footprint.
250 bbox = record.getFootprint().getBBox()
252 # Check that the bounding box has non-zero area.
253 if bbox.getArea() == 0:
254 raise measBase.MeasurementError(self.NO_PIXELS.doc, self.NO_PIXELS.number)
256 # Ensure that the centroid is within the bounding box.
257 if not bbox.contains(Point2I(center)):
258 raise measBase.MeasurementError(self.NOT_CONTAINED.doc, self.NOT_CONTAINED.number)
260 # Get the trace radius of the PSF.
261 psfSigma = exposure.getPsf().computeShape(center, color).getTraceRadius()
263 # Turn bounding box corners into GalSim bounds.
264 xmin, xmax = bbox.getMinX(), bbox.getMaxX()
265 ymin, ymax = bbox.getMinY(), bbox.getMaxY()
266 bounds = galsim._BoundsI(xmin, xmax, ymin, ymax)
268 # Get the `lsst.meas.base` mask for bad pixels.
269 badpix = exposure.mask[bbox].array.copy()
270 bitValue = exposure.mask.getPlaneBitMask(self.config.badMaskPlanes)
271 badpix &= bitValue
273 # Extract the numpy array underlying the image within the bounding box
274 # of the source.
275 imageArray = exposure.image[bbox].array
277 # Create a GalSim image using the extracted array.
278 # NOTE: GalSim's HSM uses the FITS convention of 1,1 for the
279 # lower-left corner.
280 image = galsim._Image(imageArray, bounds, None)
282 # Convert the mask of bad pixels to a format suitable for galsim.
283 gd = badpix == 0
284 badpix[gd] = 1
285 badpix[~gd] = 0
287 weight_image = galsim._Image(badpix, bounds, None)
289 # Call the internal method to calculate adaptive moments using GalSim.
290 self._calculate(
291 record,
292 image=image,
293 weight_image=weight_image,
294 sigma=2.5 * psfSigma,
295 precision=1.0e-6,
296 centroid=center,
297 )
300class HsmSourceMomentsRoundConfig(HsmSourceMomentsConfig):
301 """Configuration for HSM adaptive moments measurement for sources using
302 round weight function.
303 """
305 def setDefaults(self):
306 super().setDefaults()
307 self.roundMoments = True
309 def validate(self):
310 if not self.roundMoments:
311 raise pexConfig.FieldValidationError(
312 self.__class__.roundMoments, self, "roundMoments should be set to `True`."
313 )
314 super().validate()
317@measBase.register("ext_shapeHSM_HsmSourceMomentsRound")
318class HsmSourceMomentsRoundPlugin(HsmSourceMomentsPlugin):
319 """Plugin for HSM adaptive moments measurement for sources using round
320 weight function.
321 """
323 ConfigClass = HsmSourceMomentsRoundConfig
326class HsmPsfMomentsConfig(HsmMomentsConfig):
327 """Configuration for HSM adaptive moments measurement for PSFs."""
329 useSourceCentroidOffset = pexConfig.Field[bool](
330 doc="If True, then draw the PSF to be measured in the coordinate "
331 "system of the original image (the PSF model origin - which is "
332 "commonly the PSF centroid - may end up near a pixel edge or corner). "
333 "If False, then draw the PSF to be measured in a shifted coordinate "
334 "system such that the PSF model origin lands precisely in the center "
335 "of the central pixel of the PSF image.",
336 default=False,
337 )
339 def setDefaults(self):
340 super().setDefaults()
341 self.subtractCenter = True
344@measBase.register("ext_shapeHSM_HsmPsfMoments")
345class HsmPsfMomentsPlugin(HsmMomentsPlugin):
346 """Plugin for HSM adaptive moments measurement for PSFs."""
348 ConfigClass = HsmPsfMomentsConfig
349 _debiased = False
351 def __init__(self, config, name, schema, metadata, logName=None):
352 super().__init__(config, name, schema, metadata, logName=logName)
353 docPrefix = "Debiased centroid" if self._debiased else "Centroid"
354 self.centroidResultKey = afwTable.Point2DKey.addFields(
355 schema, name, docPrefix + " of the PSF via the HSM shape algorithm", "pixel"
356 )
357 docPrefix = "Debiased adaptive" if self._debiased else "Adaptive"
358 self.shapeKey = afwTable.QuadrupoleKey.addFields(
359 schema,
360 name,
361 docPrefix + " moments of the PSF via the HSM shape algorithm",
362 afwTable.CoordinateType.PIXEL,
363 )
364 if config.addFlux:
365 self.fluxKey = schema.addField(
366 schema.join(name, "Flux"),
367 type=float,
368 doc="Flux of the PSF via the HSM shape algorithm",
369 )
371 def measure(self, record, exposure):
372 """
373 Measure adaptive moments of the PSF given an exposure and set the
374 results in the record in place.
376 Parameters
377 ----------
378 record : `~lsst.afw.table.SourceRecord`
379 The record where measurement outputs will be stored.
380 exposure : `~lsst.afw.image.Exposure`
381 The exposure containing the PSF which needs measurement.
383 Raises
384 ------
385 MeasurementError
386 Raised for errors in measurement.
387 """
388 # Extract the centroid from the record.
389 center = self.centroidExtractor(record, self.flagHandler)
390 color = colorExtractor(record)
392 # Retrieve the PSF from the exposure.
393 psf = exposure.getPsf()
395 # Check that the PSF is not None.
396 if not psf:
397 raise measBase.MeasurementError(self.NO_PSF.doc, self.NO_PSF.number)
399 # Get the bounding box of the PSF.
400 psfBBox = psf.computeImageBBox(center, color)
402 # Two methods for getting PSF image evaluated at the source centroid:
403 if self.config.useSourceCentroidOffset:
404 # 1. Using `computeImage()` returns an image in the same coordinate
405 # system as the pixelized image.
406 psfImage = psf.computeImage(center, color)
407 else:
408 psfImage = psf.computeKernelImage(center, color)
409 # 2. Using `computeKernelImage()` to return an image does not
410 # retain any information about the original bounding box of the
411 # PSF. We therefore reset the origin to be the same as the
412 # pixelized image.
413 psfImage.setXY0(psfBBox.getMin())
415 # Get the trace radius of the PSF.
416 psfSigma = psf.computeShape(center, color).getTraceRadius()
418 # Get the bounding box in the parent coordinate system.
419 bbox = psfImage.getBBox(afwImage.PARENT)
421 # Turn bounding box corners into GalSim bounds.
422 xmin, xmax = bbox.getMinX(), bbox.getMaxX()
423 ymin, ymax = bbox.getMinY(), bbox.getMaxY()
424 bounds = galsim._BoundsI(xmin, xmax, ymin, ymax)
426 # Adjust the psfImage for noise as needed, and retrieve the mask of bad
427 # pixels.
428 badpix = self._adjustNoise(psfImage, psfBBox, exposure, record)
430 # Extract the numpy array underlying the PSF image.
431 imageArray = psfImage.array
433 # Create a GalSim image using the PSF image array.
434 image = galsim._Image(imageArray, bounds, None)
436 if badpix is not None:
437 gd = badpix == 0
438 badpix[gd] = 1
439 badpix[~gd] = 0
441 weight_image = galsim._Image(badpix, bounds, None)
442 else:
443 arr = np.ones(imageArray.shape, dtype=np.int32)
444 weight_image = galsim._Image(arr, bounds, None)
446 # Decide on the centroid position based on configuration.
447 if self.config.useSourceCentroidOffset:
448 # If the source centroid offset should be used, use the source
449 # centroid.
450 centroid = center
451 else:
452 # Otherwise, use the center of the bounding box of psfImage.
453 centroid = Point2D(psfBBox.getMin() + psfBBox.getDimensions() // 2)
455 # Call the internal method to calculate adaptive moments using GalSim.
456 self._calculate(
457 record,
458 image=image,
459 weight_image=weight_image,
460 sigma=psfSigma,
461 centroid=centroid,
462 )
464 def _adjustNoise(self, *args) -> None:
465 """A noop in the base class, returning None for the bad pixel mask.
466 This method is designed to be overridden in subclasses."""
467 pass
470class HsmPsfMomentsDebiasedConfig(HsmPsfMomentsConfig):
471 """Configuration for debiased HSM adaptive moments measurement for PSFs."""
473 noiseSource = pexConfig.ChoiceField[str](
474 doc="Noise source. How to choose variance of the zero-mean Gaussian noise added to image.",
475 allowed={
476 "meta": "variance = the 'BGMEAN' metadata entry",
477 "variance": "variance = the image's variance plane",
478 },
479 default="variance",
480 )
481 seedOffset = pexConfig.Field[int](doc="Seed offset for random number generator.", default=0)
482 badMaskPlanes = pexConfig.ListField[str](
483 doc="Mask planes used to reject bad pixels.", default=["BAD", "SAT"]
484 )
486 def setDefaults(self):
487 super().setDefaults()
488 self.useSourceCentroidOffset = True
491@measBase.register("ext_shapeHSM_HsmPsfMomentsDebiased")
492class HsmPsfMomentsDebiasedPlugin(HsmPsfMomentsPlugin):
493 """Plugin for debiased HSM adaptive moments measurement for PSFs."""
495 ConfigClass = HsmPsfMomentsDebiasedConfig
496 _debiased = True
498 @classmethod
499 def getExecutionOrder(cls):
500 # Since the standard execution order increases in steps of 1, it's
501 # safer to keep the increase by hand to less than 1. The exact value
502 # does not matter.
503 return cls.FLUX_ORDER + 0.1
505 def _adjustNoise(
506 self,
507 psfImage: afwImage.Image,
508 psfBBox: Box2I,
509 exposure: afwImage.Exposure,
510 record: afwTable.SourceRecord,
511 ) -> np.ndarray:
512 """
513 Adjusts noise in the PSF image and updates the bad pixel mask based on
514 exposure data. This method modifies `psfImage` in place and returns a
515 newly created `badpix` mask.
517 Parameters
518 ----------
519 psfImage : `~lsst.afw.image.Image`
520 The PSF image to be adjusted. This image is modified in place.
521 psfBBox : `~lsst.geom.Box2I`
522 The bounding box of the PSF.
523 exposure : `~lsst.afw.image.Exposure`
524 The exposure object containing relevant metadata and mask
525 information.
526 record : `~lsst.afw.table.SourceRecord`
527 The source record where measurement outputs will be stored. May be
528 modified in place to set flags.
530 Returns
531 -------
532 badpix : `~np.ndarray`
533 Numpy image array (np.int32) representing bad pixels, where zero
534 indicates good pixels and any nonzero value denotes a bad pixel.
536 Raises
537 ------
538 MeasurementError
539 If there's an issue during the noise adjustment process.
540 FatalAlgorithmError
541 If BGMEAN is not present in the metadata when using the meta noise
542 source.
543 """
544 # Psf image crossing exposure edge is fine if we're getting the
545 # variance from metadata, but not okay if we're getting the
546 # variance from the variance plane. In both cases, set the EDGE
547 # flag, but only fail hard if using variance plane.
548 overlap = psfImage.getBBox()
549 overlap.clip(exposure.getBBox())
550 if overlap != psfImage.getBBox():
551 self.flagHandler.setValue(record, self.EDGE.number, True)
552 if self.config.noiseSource == "variance":
553 self.flagHandler.setValue(record, self.FAILURE.number, True)
554 raise measBase.MeasurementError(self.EDGE.doc, self.EDGE.number)
556 # Match PSF flux to source.
557 psfImage *= record.getPsfInstFlux()
559 # Add Gaussian noise to image in 4 steps:
560 # 1. Initialize the noise image and random number generator.
561 noise = afwImage.Image(psfImage.getBBox(), dtype=psfImage.dtype, initialValue=0.0)
562 seed = record.getId() + self.config.seedOffset
563 rand = afwMath.Random("MT19937", seed)
565 # 2. Generate Gaussian noise image.
566 afwMath.randomGaussianImage(noise, rand)
568 # 3. Determine the noise scaling based on the noise source.
569 if self.config.noiseSource == "meta":
570 # Retrieve BGMEAN from the exposure metadata.
571 try:
572 bgmean = exposure.getMetadata().getAsDouble("BGMEAN")
573 except NotFoundError as error:
574 raise measBase.FatalAlgorithmError(str(error))
575 # Scale the noise by the square root of the background mean.
576 noise *= np.sqrt(bgmean)
577 elif self.config.noiseSource == "variance":
578 # Get the variance image from the exposure and restrict to the
579 # PSF bounding box.
580 var = afwImage.Image(exposure.variance[psfImage.getBBox()], dtype=psfImage.dtype, deep=True)
581 # Scale the noise by the square root of the variance.
582 var.sqrt() # In-place square root.
583 noise *= var
585 # 4. Add the scaled noise to the PSF image.
586 psfImage += noise
588 # Masking is needed for debiased PSF moments.
589 badpix = afwImage.Mask(psfBBox)
590 # NOTE: We repeat the `overlap` calculation in the two lines below to
591 # align with the old C++ version and minimize potential discrepancies.
592 # There's zero chance this will be a time sink, and using the bbox from
593 # the image that's about to be cropped seems safer than using the bbox
594 # from a different image, even if they're nominally supposed to have
595 # the same bounds.
596 overlap = badpix.getBBox()
597 overlap.clip(exposure.getBBox())
598 badpix[overlap] = exposure.mask[overlap]
599 # Pull out the numpy view of the badpix mask image.
600 badpix = badpix.array
602 bitValue = exposure.mask.getPlaneBitMask(self.config.badMaskPlanes)
603 badpix &= bitValue
605 return badpix