Coverage for tests / test_hsm.py: 12%
606 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-04-30 09:05 +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 itertools
23import os
24import unittest
26import galsim
27import lsst.afw.detection as afwDetection
28import lsst.afw.geom as afwGeom
29import lsst.afw.geom.ellipses as afwEll
30import lsst.afw.image as afwImage
31import lsst.afw.math as afwMath
32import lsst.afw.table as afwTable
33import lsst.geom as geom
34import lsst.meas.algorithms as algorithms
35import lsst.meas.base as base
36import lsst.meas.base.tests
37import lsst.meas.extensions.shapeHSM as shapeHSM
38import lsst.pex.config as pexConfig
39import lsst.utils.tests
40import numpy as np
41from lsst.daf.base import PropertySet
43SIZE_DECIMALS = 2 # Number of decimals for equality in sizes
44SHAPE_DECIMALS = 3 # Number of decimals for equality in shapes
46# The following values are pulled directly from GalSim's test_hsm.py:
47file_indices = [0, 2, 4, 6, 8]
48x_centroid = [35.888, 19.44, 8.74, 20.193, 57.94]
49y_centroid = [19.845, 25.047, 11.92, 38.93, 27.73]
50sky_var = [35.01188, 35.93418, 35.15456, 35.11146, 35.16454]
51correction_methods = ["KSB", "BJ", "LINEAR", "REGAUSS"]
52# Note: expected results give shear for KSB and distortion for others, but the results below have
53# converted KSB expected results to distortion for the sake of consistency
54e1_expected = np.array(
55 [
56 [0.467603106752, 0.381211727, 0.398856937, 0.401755571],
57 [0.28618443944, 0.199222784, 0.233883543, 0.234257525],
58 [0.271533794146, 0.158049396, 0.183517068, 0.184893412],
59 [-0.293754156071, -0.457024541, 0.123946584, -0.609233462],
60 [0.557720893779, 0.374143023, 0.714147448, 0.435404409],
61 ]
62)
63e2_expected = np.array(
64 [
65 [-0.867225166489, -0.734855778, -0.777027588, -0.774684891],
66 [-0.469354341577, -0.395520479, -0.502540961, -0.464466257],
67 [-0.519775291311, -0.471589061, -0.574750641, -0.529664935],
68 [0.345688365839, -0.342047099, 0.120603755, -0.44609129428863525],
69 [0.525728304099, 0.370691830, 0.702724807, 0.433999442],
70 ]
71)
72resolution_expected = np.array(
73 [
74 [0.796144249, 0.835624917, 0.835624917, 0.827796187],
75 [0.685023735, 0.699602704, 0.699602704, 0.659457638],
76 [0.634736458, 0.651040481, 0.651040481, 0.614663396],
77 [0.477027015, 0.477210752, 0.477210752, 0.423157447],
78 [0.595205998, 0.611824797, 0.611824797, 0.563582092],
79 ]
80)
81sigma_e_expected = np.array(
82 [
83 [0.016924826, 0.014637648, 0.014637648, 0.014465546],
84 [0.075769504, 0.073602324, 0.073602324, 0.064414520],
85 [0.110253112, 0.106222900, 0.106222900, 0.099357106],
86 [0.185276702, 0.184300955, 0.184300955, 0.173478300],
87 [0.073020065, 0.070270966, 0.070270966, 0.061856263],
88 ]
89)
90# End of GalSim's values
92# These values calculated using GalSim's HSM as part of GalSim
93galsim_e1 = np.array(
94 [
95 [0.399292618036, 0.381213068962, 0.398856908083, 0.401749581099],
96 [0.155929282308, 0.199228107929, 0.233882278204, 0.234371587634],
97 [0.150018423796, 0.158052951097, 0.183515056968, 0.184561833739],
98 [-2.6984937191, -0.457033962011, 0.123932465911, -0.60886412859],
99 [0.33959621191, 0.374140143394, 0.713756918907, 0.43560180068],
100 ]
101)
102galsim_e2 = np.array(
103 [
104 [-0.74053555727, -0.734855830669, -0.777024209499, -0.774700462818],
105 [-0.25573053956, -0.395517915487, -0.50251352787, -0.464388132095],
106 [-0.287168383598, -0.471584022045, -0.574719130993, -0.5296921134],
107 [3.1754450798, -0.342054128647, 0.120592080057, -0.446093201637],
108 [0.320115834475, 0.370669454336, 0.702303349972, 0.433968126774],
109 ]
110)
111galsim_resolution = np.array(
112 [
113 [0.79614430666, 0.835625052452, 0.835625052452, 0.827822327614],
114 [0.685023903847, 0.699601829052, 0.699601829052, 0.659438848495],
115 [0.634736537933, 0.651039719582, 0.651039719582, 0.614759743214],
116 [0.477026551962, 0.47721144557, 0.47721144557, 0.423227936029],
117 [0.595205545425, 0.611821532249, 0.611821532249, 0.563564240932],
118 ]
119)
120galsim_err = np.array(
121 [
122 [0.0169247947633, 0.0146376201883, 0.0146376201883, 0.0144661813974],
123 [0.0757696777582, 0.0736026018858, 0.0736026018858, 0.0644160583615],
124 [0.110252402723, 0.106222368777, 0.106222368777, 0.0993555411696],
125 [0.185278102756, 0.184301897883, 0.184301897883, 0.17346136272],
126 [0.0730196461082, 0.0702708885074, 0.0702708885074, 0.0618583671749],
127 ]
128)
130moments_expected = np.array(
131 [ # sigma, e1, e2
132 [2.24490427971, 0.336240686301, -0.627372910656],
133 [1.9031778574, 0.150566105384, -0.245272792302],
134 [1.77790760994, 0.112286123389, -0.286203939641],
135 [1.45464873314, -0.155597168978, -0.102008266223],
136 [1.63144648075, 0.22886961923, 0.228813588897],
137 ]
138)
139centroid_expected = np.array(
140 [ # x, y
141 [36.218247328, 20.5678722157],
142 [20.325744838, 25.4176650386],
143 [9.54257706283, 12.6134786199],
144 [20.6407850048, 39.5864802706],
145 [58.5008586442, 28.2850942049],
146 ]
147)
149round_moments_expected = np.array(
150 [ # sigma, e1, e2, flux, x, y
151 [2.40270376205, 0.197810277343, -0.372329413891, 3740.22436523, 36.4032272633, 20.4847916447],
152 [1.89714717865, 0.046496052295, -0.0987404286861, 776.709594727, 20.2893584046, 25.4230368047],
153 [1.77995181084, 0.0416346564889, -0.143147706985, 534.59197998, 9.51994111869, 12.6250775205],
154 [1.46549296379, -0.0831127092242, -0.0628845766187, 348.294403076, 20.6242279632, 39.5941625731],
155 [1.64031589031, 0.0867517963052, 0.0940798297524, 793.374450684, 58.4728765002, 28.2686937854],
156 ]
157)
160def makePluginAndCat(alg, name, control=None, metadata=False, centroid=None, psfflux=None, addFlux=False):
161 if control is None:
162 control = alg.ConfigClass()
163 if addFlux:
164 control.addFlux = True
165 schema = afwTable.SourceTable.makeMinimalSchema()
166 if centroid:
167 lsst.afw.table.Point2DKey.addFields(schema, centroid, "centroid", "pixel")
168 schema.getAliasMap().set("slot_Centroid", centroid)
169 if psfflux:
170 base.PsfFluxAlgorithm(base.PsfFluxControl(), psfflux, schema)
171 schema.getAliasMap().set("slot_PsfFlux", psfflux)
172 if metadata:
173 plugin = alg(control, name, schema, PropertySet())
174 else:
175 plugin = alg(control, name, schema)
176 cat = afwTable.SourceCatalog(schema)
177 if centroid:
178 cat.defineCentroid(centroid)
179 return plugin, cat
182class MomentsTestCase(unittest.TestCase):
183 """A test case for shape measurement"""
185 def setUp(self):
186 # load the known values
187 self.dataDir = os.path.join(os.getenv("MEAS_EXTENSIONS_SHAPEHSM_DIR"), "tests", "data")
188 self.bkgd = 1000.0 # standard for atlas image
189 self.offset = geom.Extent2I(1234, 1234)
190 self.xy0 = geom.Point2I(5678, 9876)
192 def tearDown(self):
193 del self.offset
194 del self.xy0
196 def runMeasurement(self, algorithmName, imageid, x, y, v, addFlux=False, maskAll=False):
197 """Run the measurement algorithm on an image"""
198 # load the test image
199 imgFile = os.path.join(self.dataDir, "image.%d.fits" % imageid)
200 img = afwImage.ImageF(imgFile)
201 img -= self.bkgd
202 nx, ny = img.getWidth(), img.getHeight()
203 msk = afwImage.Mask(geom.Extent2I(nx, ny), 0x0)
204 var = afwImage.ImageF(geom.Extent2I(nx, ny), v)
205 mimg = afwImage.MaskedImageF(img, msk, var)
206 msk.getArray()[:] = np.where(np.fabs(img.getArray()) < 1.0e-8, msk.getPlaneBitMask("BAD"), 0)
207 if maskAll:
208 msk.array[:] |= msk.getPlaneBitMask("BAD")
210 # Put it in a bigger image, in case it matters
211 big = afwImage.MaskedImageF(self.offset + mimg.getDimensions())
212 big.getImage().set(0)
213 big.getMask().set(0)
214 big.getVariance().set(v)
215 subBig = afwImage.MaskedImageF(big, geom.Box2I(big.getXY0() + self.offset, mimg.getDimensions()))
216 subBig.assign(mimg)
217 mimg = big
218 mimg.setXY0(self.xy0)
220 exposure = afwImage.makeExposure(mimg)
221 cdMatrix = np.array([1.0 / (2.53 * 3600.0), 0.0, 0.0, 1.0 / (2.53 * 3600.0)])
222 cdMatrix.shape = (2, 2)
223 exposure.setWcs(
224 afwGeom.makeSkyWcs(
225 crpix=geom.Point2D(1.0, 1.0), crval=geom.SpherePoint(0, 0, geom.degrees), cdMatrix=cdMatrix
226 )
227 )
229 # load the corresponding test psf
230 psfFile = os.path.join(self.dataDir, "psf.%d.fits" % imageid)
231 psfImg = afwImage.ImageD(psfFile)
232 psfImg -= self.bkgd
234 kernel = afwMath.FixedKernel(psfImg)
235 kernelPsf = algorithms.KernelPsf(kernel)
236 exposure.setPsf(kernelPsf)
238 # perform the shape measurement
239 msConfig = base.SingleFrameMeasurementConfig()
240 msConfig.plugins.names |= [algorithmName]
241 control = msConfig.plugins[algorithmName]
242 alg = base.SingleFramePlugin.registry[algorithmName].PluginClass
243 # NOTE: It is essential to remove the floating point part of the position for the
244 # Algorithm._apply. Otherwise, when the PSF is realised it will have been warped
245 # to account for the sub-pixel offset and we won't get *exactly* this PSF.
246 plugin, table = makePluginAndCat(
247 alg, algorithmName, control, centroid="centroid", metadata=True, addFlux=addFlux
248 )
249 center = geom.Point2D(int(x), int(y)) + geom.Extent2D(self.offset + geom.Extent2I(self.xy0))
250 source = table.makeRecord()
251 source.set("centroid_x", center.getX())
252 source.set("centroid_y", center.getY())
253 source.setFootprint(afwDetection.Footprint(afwGeom.SpanSet(exposure.getBBox(afwImage.PARENT))))
254 plugin.measure(source, exposure)
256 return source
258 def testValidateHsmSourceMomentsRoundConfig(self):
259 algorithmName = "ext_shapeHSM_HsmSourceMomentsRound"
260 msConfig = base.SingleFrameMeasurementConfig()
261 msConfig.plugins.names |= [algorithmName]
262 control = msConfig.plugins[algorithmName]
263 control.validate() # This should not raise any error.
264 with self.assertRaises(pexConfig.FieldValidationError):
265 control.roundMoments = False
266 control.validate()
268 def testHsmSourceMoments(self):
269 for i, imageid in enumerate(file_indices):
270 source = self.runMeasurement(
271 "ext_shapeHSM_HsmSourceMoments", imageid, x_centroid[i], y_centroid[i], sky_var[i]
272 )
273 x = source.get("ext_shapeHSM_HsmSourceMoments_x")
274 y = source.get("ext_shapeHSM_HsmSourceMoments_y")
275 xx = source.get("ext_shapeHSM_HsmSourceMoments_xx")
276 yy = source.get("ext_shapeHSM_HsmSourceMoments_yy")
277 xy = source.get("ext_shapeHSM_HsmSourceMoments_xy")
279 # Centroids from GalSim use the FITS lower-left corner of 1,1
280 offset = self.xy0 + self.offset
281 self.assertAlmostEqual(x - offset.getX(), centroid_expected[i][0] - 1, 3)
282 self.assertAlmostEqual(y - offset.getY(), centroid_expected[i][1] - 1, 3)
284 expected = afwEll.Quadrupole(
285 afwEll.SeparableDistortionDeterminantRadius(
286 moments_expected[i][1], moments_expected[i][2], moments_expected[i][0]
287 )
288 )
290 self.assertAlmostEqual(xx, expected.getIxx(), SHAPE_DECIMALS)
291 self.assertAlmostEqual(xy, expected.getIxy(), SHAPE_DECIMALS)
292 self.assertAlmostEqual(yy, expected.getIyy(), SHAPE_DECIMALS)
294 def testHsmSourceMomentsRound(self):
295 for i, imageid in enumerate(file_indices):
296 source = self.runMeasurement(
297 "ext_shapeHSM_HsmSourceMomentsRound",
298 imageid,
299 x_centroid[i],
300 y_centroid[i],
301 sky_var[i],
302 addFlux=True,
303 )
304 x = source.get("ext_shapeHSM_HsmSourceMomentsRound_x")
305 y = source.get("ext_shapeHSM_HsmSourceMomentsRound_y")
306 xx = source.get("ext_shapeHSM_HsmSourceMomentsRound_xx")
307 yy = source.get("ext_shapeHSM_HsmSourceMomentsRound_yy")
308 xy = source.get("ext_shapeHSM_HsmSourceMomentsRound_xy")
309 flux = source.get("ext_shapeHSM_HsmSourceMomentsRound_Flux")
311 # Centroids from GalSim use the FITS lower-left corner of 1,1
312 offset = self.xy0 + self.offset
313 self.assertAlmostEqual(x - offset.getX(), round_moments_expected[i][4] - 1, 3)
314 self.assertAlmostEqual(y - offset.getY(), round_moments_expected[i][5] - 1, 3)
316 expected = afwEll.Quadrupole(
317 afwEll.SeparableDistortionDeterminantRadius(
318 round_moments_expected[i][1], round_moments_expected[i][2], round_moments_expected[i][0]
319 )
320 )
321 self.assertAlmostEqual(xx, expected.getIxx(), SHAPE_DECIMALS)
322 self.assertAlmostEqual(xy, expected.getIxy(), SHAPE_DECIMALS)
323 self.assertAlmostEqual(yy, expected.getIyy(), SHAPE_DECIMALS)
325 self.assertAlmostEqual(flux, round_moments_expected[i][3], SHAPE_DECIMALS)
327 def testHsmSourceMomentsVsSdssShape(self):
328 # Initialize a config and activate the plugins.
329 sfmConfig = base.SingleFrameMeasurementConfig()
330 sfmConfig.plugins.names |= ["ext_shapeHSM_HsmSourceMoments", "base_SdssShape"]
332 # Create a minimal schema (columns).
333 schema = lsst.meas.base.tests.TestDataset.makeMinimalSchema()
335 # Instantiate the task.
336 sfmTask = base.SingleFrameMeasurementTask(config=sfmConfig, schema=schema)
338 # Create a simple, test dataset.
339 bbox = lsst.geom.Box2I(lsst.geom.Point2I(0, 0), lsst.geom.Extent2I(100, 100))
340 dataset = lsst.meas.base.tests.TestDataset(bbox)
342 # First source is a point.
343 dataset.addSource(100000.0, lsst.geom.Point2D(49.5, 49.5))
345 # Second source is a galaxy.
346 dataset.addSource(300000.0, lsst.geom.Point2D(76.3, 79.2), afwGeom.Quadrupole(2.0, 3.0, 0.5))
348 # Third source is also a galaxy.
349 dataset.addSource(250000.0, lsst.geom.Point2D(28.9, 41.35), afwGeom.Quadrupole(1.8, 3.5, 0.4))
351 # Get the exposure and catalog.
352 exposure, catalog = dataset.realize(10.0, sfmTask.schema, randomSeed=0)
354 # Run the measurement task to get the output catalog.
355 sfmTask.run(catalog, exposure)
356 cat = catalog.asAstropy()
358 # Get the moments from the catalog.
359 xSdss, ySdss = cat["base_SdssShape_x"], cat["base_SdssShape_y"]
360 xxSdss, xySdss, yySdss = cat["base_SdssShape_xx"], cat["base_SdssShape_xy"], cat["base_SdssShape_yy"]
361 xHsm, yHsm = cat["ext_shapeHSM_HsmSourceMoments_x"], cat["ext_shapeHSM_HsmSourceMoments_y"]
362 xxHsm, xyHsm, yyHsm = (
363 cat["ext_shapeHSM_HsmSourceMoments_xx"],
364 cat["ext_shapeHSM_HsmSourceMoments_xy"],
365 cat["ext_shapeHSM_HsmSourceMoments_yy"],
366 )
368 # Loop over the sources and check that the moments are the same.
369 for i in range(3):
370 self.assertAlmostEqual(xSdss[i], xHsm[i], 2)
371 self.assertAlmostEqual(ySdss[i], yHsm[i], 2)
372 self.assertAlmostEqual(xxSdss[i], xxHsm[i], SHAPE_DECIMALS)
373 self.assertAlmostEqual(xySdss[i], xyHsm[i], SHAPE_DECIMALS)
374 self.assertAlmostEqual(yySdss[i], yyHsm[i], SHAPE_DECIMALS)
376 def testHsmSourceMomentsAllMasked(self):
377 i = 0
378 imageid = file_indices[0]
379 with self.assertRaises(base.MeasurementError):
380 _ = self.runMeasurement(
381 "ext_shapeHSM_HsmSourceMoments",
382 imageid,
383 x_centroid[i],
384 y_centroid[i],
385 sky_var[i],
386 maskAll=True,
387 )
390class ShapeTestCase(unittest.TestCase):
391 """A test case for shape measurement"""
393 def setUp(self):
394 # load the known values
395 self.dataDir = os.path.join(os.getenv("MEAS_EXTENSIONS_SHAPEHSM_DIR"), "tests", "data")
396 self.bkgd = 1000.0 # standard for atlas image
397 self.offset = geom.Extent2I(1234, 1234)
398 self.xy0 = geom.Point2I(5678, 9876)
400 def tearDown(self):
401 del self.offset
402 del self.xy0
404 @staticmethod
405 def computeDirectShapeFromGalSim(record, exposure, config):
406 """
407 Retrieve the shape as estimated directly by GalSim for comparison
408 purposes.
410 Parameters
411 ----------
412 record : `~lsst.afw.table.SourceRecord`
413 The record containing the center and footprint of the source which
414 needs measurement.
415 exposure : `~lsst.afw.image.Exposure`
416 The exposure containing the source which needs measurement.
417 config : `~lsst.meas.extensions.shapeHSM._hsm_shape.\
418 HsmShapeConfig`
419 The configuration object containing parameters and settings for
420 this measurement. This needs to be a subclass in the format
421 HsmShape<Method>Config, where <Method> represents the name of the
422 correction method being utilized (e.g., Ksb, Regauss, etc.).
424 Returns
425 -------
426 shapeDirect : `~galsim.hsm.ShapeData`
427 An object containing the results of shape measurement.
428 """
430 # Get the center of the source as a Point2D.
431 center = geom.Point2D(record.get("centroid_x"), record.get("centroid_y"))
433 # Get the PSF image evaluated at the source centroid.
434 psfImage = exposure.getPsf().computeImage(center)
435 psfImage.setXY0(0, 0)
437 # Get the GalSim images to use in the EstimateShear call.
438 bbox = record.getFootprint().getBBox()
439 bounds = galsim.BoundsI(bbox.getMinX(), bbox.getMaxX(), bbox.getMinY(), bbox.getMaxY())
440 image = galsim.Image(exposure.image[bbox].array, bounds=bounds)
441 psfBBox = psfImage.getBBox(afwImage.PARENT)
442 psfBounds = galsim.BoundsI(psfBBox.getMinX(), psfBBox.getMaxX(), psfBBox.getMinY(), psfBBox.getMaxY())
443 psf = galsim.Image(psfImage.array, bounds=psfBounds)
445 # Get the mask of bad pixels.
446 subMask = exposure.mask[bbox]
447 badpix = subMask.array.copy() # Copy it since badpix gets modified.
448 bitValue = exposure.mask.getPlaneBitMask(config.badMaskPlanes)
449 badpix &= bitValue
450 badpix = galsim.Image(badpix, bounds=bounds)
452 # Estimate the sky variance.
453 sctrl = afwMath.StatisticsControl()
454 sctrl.setAndMask(bitValue)
455 variance = afwImage.Image(
456 exposure.variance[bbox],
457 dtype=exposure.variance.dtype,
458 deep=True,
459 )
460 stat = afwMath.makeStatistics(variance, subMask, afwMath.MEDIAN, sctrl)
461 skyvar = stat.getValue(afwMath.MEDIAN)
463 # Prepare various values for GalSim's EstimateShear.
464 recomputeFlux = "FIT"
465 precision = 1.0e-6
466 psfSigma = exposure.getPsf().computeShape(center).getTraceRadius()
467 guessCentroid = galsim.PositionD(center.x, center.y)
469 # Estimate the shape using GalSim's Python interface.
470 shapeDirect = galsim.hsm.EstimateShear(
471 gal_image=image,
472 PSF_image=psf,
473 weight=None,
474 badpix=badpix,
475 sky_var=skyvar,
476 shear_est=config.shearType.upper(),
477 recompute_flux=recomputeFlux.upper(),
478 guess_sig_gal=2.5 * psfSigma,
479 guess_sig_PSF=psfSigma,
480 precision=precision,
481 guess_centroid=guessCentroid,
482 hsmparams=None,
483 )
484 return shapeDirect
486 def runMeasurement(self, algorithmName, imageid, x, y, v):
487 """Run the measurement algorithm on an image"""
488 # load the test image
489 imgFile = os.path.join(self.dataDir, "image.%d.fits" % imageid)
490 img = afwImage.ImageF(imgFile)
491 img -= self.bkgd
492 nx, ny = img.getWidth(), img.getHeight()
493 msk = afwImage.Mask(geom.Extent2I(nx, ny), 0x0)
494 var = afwImage.ImageF(geom.Extent2I(nx, ny), v)
495 mimg = afwImage.MaskedImageF(img, msk, var)
496 msk.getArray()[:] = np.where(np.fabs(img.getArray()) < 1.0e-8, msk.getPlaneBitMask("BAD"), 0)
498 # Put it in a bigger image, in case it matters
499 big = afwImage.MaskedImageF(self.offset + mimg.getDimensions())
500 big.getImage().set(0)
501 big.getMask().set(0)
502 big.getVariance().set(v)
503 subBig = afwImage.MaskedImageF(big, geom.Box2I(big.getXY0() + self.offset, mimg.getDimensions()))
504 subBig.assign(mimg)
505 mimg = big
506 mimg.setXY0(self.xy0)
508 exposure = afwImage.makeExposure(mimg)
509 cdMatrix = np.array([1.0 / (2.53 * 3600.0), 0.0, 0.0, 1.0 / (2.53 * 3600.0)])
510 cdMatrix.shape = (2, 2)
511 exposure.setWcs(
512 afwGeom.makeSkyWcs(
513 crpix=geom.Point2D(1.0, 1.0), crval=geom.SpherePoint(0, 0, geom.degrees), cdMatrix=cdMatrix
514 )
515 )
517 # load the corresponding test psf
518 psfFile = os.path.join(self.dataDir, "psf.%d.fits" % imageid)
519 psfImg = afwImage.ImageD(psfFile)
520 psfImg -= self.bkgd
522 kernel = afwMath.FixedKernel(psfImg)
523 kernelPsf = algorithms.KernelPsf(kernel)
524 exposure.setPsf(kernelPsf)
526 # perform the shape measurement
527 msConfig = base.SingleFrameMeasurementConfig()
528 msConfig.plugins.names |= [algorithmName]
529 control = msConfig.plugins[algorithmName]
530 alg = base.SingleFramePlugin.registry[algorithmName].PluginClass
531 # NOTE: It is essential to remove the floating point part of the position for the
532 # Algorithm._apply. Otherwise, when the PSF is realised it will have been warped
533 # to account for the sub-pixel offset and we won't get *exactly* this PSF.
534 plugin, table = makePluginAndCat(alg, algorithmName, control, centroid="centroid", metadata=True)
535 center = geom.Point2D(int(x), int(y)) + geom.Extent2D(self.offset + geom.Extent2I(self.xy0))
536 source = table.makeRecord()
537 source.set("centroid_x", center.getX())
538 source.set("centroid_y", center.getY())
539 source.setFootprint(afwDetection.Footprint(afwGeom.SpanSet(exposure.getBBox(afwImage.PARENT))))
540 plugin.measure(source, exposure)
542 shapeDirect = self.computeDirectShapeFromGalSim(source, exposure, control)
544 return source, alg.measTypeSymbol, shapeDirect
546 def testHsmShape(self):
547 """Test that we can instantiate and play with a measureShape"""
549 nFail = 0
550 msg = ""
552 for (algNum, algName), (i, imageid) in itertools.product(
553 enumerate(correction_methods), enumerate(file_indices)
554 ):
555 algorithmName = "ext_shapeHSM_HsmShape" + algName[0:1].upper() + algName[1:].lower()
557 source, preEstimationMeasType, shapeDirect = self.runMeasurement(
558 algorithmName, imageid, x_centroid[i], y_centroid[i], sky_var[i]
559 )
561 postEstimationMeasType = shapeDirect.meas_type
563 # Check consistency with GalSim output
564 self.assertEqual(
565 preEstimationMeasType,
566 postEstimationMeasType,
567 "The plugin setup is incompatible with GalSim output.",
568 )
570 ##########################################
571 # see how we did
572 if algName in ("KSB"):
573 # Need to convert g1,g2 --> e1,e2 because GalSim has done that
574 # for the expected values ("for consistency")
575 g1 = source.get(algorithmName + "_g1")
576 g2 = source.get(algorithmName + "_g2")
577 scale = 2.0 / (1.0 + g1**2 + g2**2)
578 e1 = g1 * scale
579 e2 = g2 * scale
580 sigma = source.get(algorithmName + "_sigma")
581 # Ensure the values calculated are identical to those obtained
582 # from GalSim.
583 self.assertEqual(g1, shapeDirect.corrected_g1)
584 self.assertEqual(g2, shapeDirect.corrected_g2)
585 else:
586 e1 = source.get(algorithmName + "_e1")
587 e2 = source.get(algorithmName + "_e2")
588 sigma = 0.5 * source.get(algorithmName + "_sigma")
589 # Ensure the values calculated are identical to those obtained
590 # from GalSim.
591 self.assertEqual(e1, shapeDirect.corrected_e1)
592 self.assertEqual(e2, shapeDirect.corrected_e2)
594 resolution = source.get(algorithmName + "_resolution")
595 flags = source.get(algorithmName + "_flag")
597 # Check that the shape error and the resolution factor are the same
598 # as GalSim's.
599 self.assertEqual(sigma, shapeDirect.corrected_shape_err)
600 self.assertEqual(resolution, shapeDirect.resolution_factor)
602 tests = [
603 # label, known-value, measured, tolerance
604 ["e1", float(e1_expected[i][algNum]), e1, 0.5 * 10**-SHAPE_DECIMALS],
605 ["e2", float(e2_expected[i][algNum]), e2, 0.5 * 10**-SHAPE_DECIMALS],
606 ["resolution", float(resolution_expected[i][algNum]), resolution, 0.5 * 10**-SIZE_DECIMALS],
607 # sigma won't match exactly because
608 # we're using skyvar=mean(var) instead of measured value ... expected a difference
609 ["sigma", float(sigma_e_expected[i][algNum]), sigma, 0.07],
610 ["shapeStatus", 0, flags, 0],
611 ]
613 for test in tests:
614 label, know, hsm, limit = test
615 err = hsm - know
616 msgTmp = "%-12s %s %5s: %6.6f %6.6f (val-known) = %.3g\n" % (
617 algName,
618 imageid,
619 label,
620 know,
621 hsm,
622 err,
623 )
624 if not np.isfinite(err) or abs(err) > limit:
625 msg += msgTmp
626 nFail += 1
628 self.assertAlmostEqual(g1 if algName in ("KSB") else e1, galsim_e1[i][algNum], SHAPE_DECIMALS)
629 self.assertAlmostEqual(g2 if algName in ("KSB") else e2, galsim_e2[i][algNum], SHAPE_DECIMALS)
630 self.assertAlmostEqual(resolution, galsim_resolution[i][algNum], SIZE_DECIMALS)
631 self.assertAlmostEqual(sigma, galsim_err[i][algNum], delta=0.07)
633 self.assertEqual(nFail, 0, "\n" + msg)
635 @lsst.utils.tests.methodParametersProduct(
636 # Increasing the width beyond 4.5 leads to noticeable truncation of the
637 # PSF, i.e. a PSF that is too large for the box. While this truncated
638 # state leads to incorrect measurements, it is necessary for testing
639 # purposes to evaluate the behavior under these extreme conditions.
640 # Increasing the width beyond 41.3 and zeroPadding beyond 32 with
641 # everything else constant fails to converge for this particular test
642 # dataset.
643 width=(2.0, 3.0, 4.0, 10.0, 40.0),
644 zeroPadding=(2, 5, 10, 20, 30),
645 varyBBox=(True, False),
646 wrongBBox=(True, False),
647 algName=correction_methods,
648 )
649 def testHsmShapeWithVariousPsfsVsDirectGalsim(self, width, zeroPadding, varyBBox, wrongBBox, algName):
650 # Set the full algorithm name.
651 algorithmName = "ext_shapeHSM_HsmShape" + algName[0:1].upper() + algName[1:].lower()
653 # Initialize a config and activate the plugins.
654 sfmConfig = base.SingleFrameMeasurementConfig()
655 sfmConfig.plugins.names |= [algorithmName]
657 # Create a minimal schema (columns).
658 schema = lsst.meas.base.tests.TestDataset.makeMinimalSchema()
660 # Create a simple, test dataset.
661 bbox = lsst.geom.Box2I(lsst.geom.Point2I(0, 0), lsst.geom.Extent2I(60, 60))
662 dataset = lsst.meas.base.tests.TestDataset(bbox)
664 # Add a galaxy.
665 center = lsst.geom.Point2D(24.9, 32.5)
666 dataset.addSource(150000.0, center, afwGeom.Quadrupole(3.0, 4.0, 0.5))
668 # Get the exposure.
669 exposure, _ = dataset.realize(noise=10.0, schema=schema, randomSeed=1746)
671 # Create and set the PSF for the exposure.
672 psf = PyGaussianPsf(
673 35,
674 35,
675 width,
676 varyBBox=varyBBox,
677 wrongBBox=wrongBBox,
678 zeroPadding=zeroPadding,
679 )
680 exposure.getMaskedImage().set(1.0, 0, 1.0)
681 exposure.setPsf(psf)
683 # Conduct the measurement directly by GalSim.
684 alg = base.SingleFramePlugin.registry[algorithmName].PluginClass
685 plugin, table = makePluginAndCat(alg, algorithmName, centroid="centroid", metadata=True)
686 record = table.makeRecord()
687 record.set("centroid_x", center.x)
688 record.set("centroid_y", center.y)
689 record.setFootprint(afwDetection.Footprint(afwGeom.SpanSet(exposure.getBBox(afwImage.PARENT))))
690 shapeDirect = self.computeDirectShapeFromGalSim(record, exposure, plugin.config)
692 # Run the shapeHSM measurement task and update the record.
693 plugin.measure(record, exposure)
695 if algName in ("KSB"):
696 g1Direct, g2Direct = shapeDirect.corrected_g1, shapeDirect.corrected_g2
697 sigmaDirect = shapeDirect.corrected_shape_err
698 g1Hsm, g2Hsm = record[algorithmName + "_g1"], record[algorithmName + "_g2"]
699 sigmaHsm = record[algorithmName + "_sigma"]
700 # Check that the answers are "identical" between the two methods.
701 self.assertEqual(g1Direct, g1Hsm)
702 self.assertEqual(g2Direct, g2Hsm)
703 self.assertEqual(sigmaDirect, sigmaHsm)
704 else:
705 e1Direct, e2Direct = shapeDirect.corrected_e1, shapeDirect.corrected_e2
706 sigmaDirect = shapeDirect.corrected_shape_err
707 e1Hsm, e2Hsm = record[algorithmName + "_e1"], record[algorithmName + "_e2"]
708 # The factor of 0.5 is because shapeHSM returns
709 # `2 * corrected_shape_err` as sigma for e-type distortions.
710 sigmaHsm = 0.5 * record[algorithmName + "_sigma"]
711 # Check that the answers are "identical" between the two methods.
712 self.assertEqual(e1Direct, e1Hsm)
713 self.assertEqual(e2Direct, e2Hsm)
714 self.assertEqual(sigmaDirect, sigmaHsm)
716 resolutionDirect = shapeDirect.resolution_factor
717 flagsDirect = shapeDirect.correction_status
718 resolutionHSM = record[algorithmName + "_resolution"]
719 flagsHSM = record[algorithmName + "_flag"]
721 # Check that the resolution factor and the correction status are
722 # exactly the same as when using GalSim directly.
723 self.assertEqual(resolutionDirect, resolutionHSM)
724 self.assertEqual(flagsDirect, flagsHSM)
726 @lsst.utils.tests.methodParametersProduct(
727 registerAllShapePlugins=(True, False),
728 )
729 def testHsmShapeConfig(self, registerAllShapePlugins=False):
730 if registerAllShapePlugins:
731 msConfig = base.SingleFrameMeasurementConfig()
732 msConfig.plugins.names |= [
733 "ext_shapeHSM_HsmShape" + algName.capitalize() for algName in correction_methods
734 ]
736 for algName in correction_methods:
737 algorithmName = "ext_shapeHSM_HsmShape" + algName.capitalize()
738 if not registerAllShapePlugins:
739 msConfig = base.SingleFrameMeasurementConfig()
740 msConfig.plugins.names |= [algorithmName]
741 control = msConfig.plugins[algorithmName]
743 # Ensure that 'shearType' is not settable under any circumstances,
744 # i.e., expect AttributeError when trying to set it.
745 for shearType in correction_methods:
746 self.assertEqual(control.shearType, algName)
747 with self.assertRaises(AttributeError):
748 control.shearType = shearType
751class PyGaussianPsf(afwDetection.Psf):
752 # Like afwDetection.GaussianPsf, but handles computeImage exactly instead of
753 # via interpolation. This is a subminimal implementation. It works for the
754 # tests here but isn't fully functional as a Psf class.
756 def __init__(self, width, height, sigma, varyBBox=False, wrongBBox=False, zeroPadding=0):
757 afwDetection.Psf.__init__(self, isFixed=not varyBBox)
758 self.zeroPadding = int(zeroPadding) # To address DM-42294
759 self.dimensions = geom.Extent2I(width + self.zeroPadding, height + self.zeroPadding)
760 self.sigma = sigma
761 self.varyBBox = varyBBox # To address DM-29863
762 self.wrongBBox = wrongBBox # To address DM-30426
764 def zeroPad(self, img):
765 # The thickness of the zero padding to encase the image edges.
766 p = self.zeroPadding # p must be a positive integer
767 # Replace the outermost p pixels of the top, bottom, left, and right
768 # edges of the image array with zeros.
769 img.array[:p, :] = 0 # Top edge
770 img.array[-p:, :] = 0 # Bottom edge
771 img.array[:, :p] = 0 # Left edge
772 img.array[:, -p:] = 0 # Right edge
773 return img
775 def _doComputeKernelImage(self, position=None, color=None):
776 bbox = self.computeBBox(position, color)
777 img = afwImage.Image(bbox, dtype=np.float64)
778 x, y = np.ogrid[bbox.minY : bbox.maxY + 1, bbox.minX : bbox.maxX + 1]
779 rsqr = x**2 + y**2
780 img.array[:] = np.exp(-0.5 * rsqr / self.sigma**2)
781 if self.zeroPadding > 0:
782 img = self.zeroPad(img)
783 img.array /= np.sum(img.array)
784 return img
786 def _doComputeImage(self, position=None, color=None):
787 bbox = self.computeBBox(position, color)
788 if self.wrongBBox:
789 # For DM-30426:
790 # Purposely make computeImage.getBBox() and computeBBox()
791 # inconsistent. Old shapeHSM code attempted to infer the former
792 # from the latter, but was unreliable. New code infers the former
793 # directly, so this inconsistency no longer breaks things.
794 bbox.shift(geom.Extent2I(1, 2))
795 img = afwImage.Image(bbox, dtype=np.float64)
796 y, x = np.ogrid[float(bbox.minY) : bbox.maxY + 1, bbox.minX : bbox.maxX + 1]
797 x -= position.x - np.floor(position.x + 0.5)
798 y -= position.y - np.floor(position.y + 0.5)
799 rsqr = x**2 + y**2
800 img.array[:] = np.exp(-0.5 * rsqr / self.sigma**2)
801 if self.zeroPadding > 0:
802 img = self.zeroPad(img)
803 img.array /= np.sum(img.array)
804 img.setXY0(
805 geom.Point2I(img.getX0() + np.floor(position.x + 0.5), img.getY0() + np.floor(position.y + 0.5))
806 )
807 return img
809 def _doComputeBBox(self, position=None, color=None):
810 # Variable size bbox for addressing DM-29863.
811 dims = self.dimensions
812 if self.varyBBox:
813 if position.x > 20.0:
814 dims = dims + geom.Extent2I(2, 2)
815 return geom.Box2I(geom.Point2I(-dims / 2), dims)
817 def _doComputeShape(self, position=None, color=None):
818 return afwGeom.ellipses.Quadrupole(self.sigma**2, self.sigma**2, 0.0)
821class PsfMomentsTestCase(unittest.TestCase):
822 """A test case for PSF moments measurement"""
824 @staticmethod
825 def computeDirectPsfMomentsFromGalSim(psf, center, useSourceCentroidOffset=False):
826 """Directly from GalSim."""
827 psfBBox = psf.computeImageBBox(center)
828 psfSigma = psf.computeShape(center).getTraceRadius()
829 if useSourceCentroidOffset:
830 psfImage = psf.computeImage(center)
831 centroid = center
832 else:
833 psfImage = psf.computeKernelImage(center)
834 psfImage.setXY0(psfBBox.getMin())
835 centroid = geom.Point2D(psfBBox.getMin() + psfBBox.getDimensions() // 2)
836 bbox = psfImage.getBBox(afwImage.PARENT)
837 bounds = galsim.BoundsI(bbox.getMinX(), bbox.getMaxX(), bbox.getMinY(), bbox.getMaxY())
838 image = galsim.Image(psfImage.array, bounds=bounds)
839 guessCentroid = galsim.PositionD(centroid.x, centroid.y)
840 shape = galsim.hsm.FindAdaptiveMom(
841 image,
842 weight=None,
843 badpix=None,
844 guess_sig=psfSigma,
845 precision=1e-6,
846 guess_centroid=guessCentroid,
847 strict=True,
848 round_moments=False,
849 hsmparams=None,
850 )
851 ellipse = lsst.afw.geom.ellipses.SeparableDistortionDeterminantRadius(
852 e1=shape.observed_shape.e1,
853 e2=shape.observed_shape.e2,
854 radius=shape.moments_sigma,
855 normalize=True, # Fail if |e|>1.
856 )
857 quad = lsst.afw.geom.ellipses.Quadrupole(ellipse)
858 ixx = quad.getIxx()
859 iyy = quad.getIyy()
860 ixy = quad.getIxy()
861 return ixx, iyy, ixy
863 @lsst.utils.tests.methodParameters(
864 # Make Cartesian product of settings to feed to methodParameters
865 **dict(
866 list(
867 zip(
868 (
869 kwargs := dict(
870 # Increasing the width beyond 4.5 leads to noticeable
871 # truncation of the PSF, i.e. a PSF that is too large for the
872 # box. While this truncated state leads to incorrect
873 # measurements, it is necessary for testing purposes to
874 # evaluate the behavior under these extreme conditions.
875 width=(2.0, 3.0, 4.0, 10.0, 40.0, 100.0),
876 useSourceCentroidOffset=(True, False),
877 varyBBox=(True, False),
878 wrongBBox=(True, False),
879 center=(
880 (23.0, 34.0), # various offsets that might cause trouble
881 (23.5, 34.0),
882 (23.5, 34.5),
883 (23.15, 34.25),
884 (22.81, 34.01),
885 (22.81, 33.99),
886 (1.2, 1.3), # psfImage extends outside exposure; that's okay
887 (-100.0, -100.0),
888 (-100.5, -100.0),
889 (-100.5, -100.5),
890 ),
891 )
892 ).keys(),
893 zip(*itertools.product(*kwargs.values())),
894 )
895 )
896 )
897 )
898 def testHsmPsfMoments(self, width, useSourceCentroidOffset, varyBBox, wrongBBox, center):
899 psf = PyGaussianPsf(35, 35, width, varyBBox=varyBBox, wrongBBox=wrongBBox)
900 exposure = afwImage.ExposureF(45, 56)
901 exposure.getMaskedImage().set(1.0, 0, 1.0)
902 exposure.setPsf(psf)
904 # perform the moment measurement
905 algorithmName = "ext_shapeHSM_HsmPsfMoments"
906 msConfig = base.SingleFrameMeasurementConfig()
907 msConfig.algorithms.names = [algorithmName]
908 control = msConfig.plugins[algorithmName]
909 alg = base.SingleFramePlugin.registry[algorithmName].PluginClass
910 self.assertFalse(control.useSourceCentroidOffset)
911 control.useSourceCentroidOffset = useSourceCentroidOffset
912 plugin, cat = makePluginAndCat(
913 alg,
914 algorithmName,
915 centroid="centroid",
916 control=control,
917 metadata=True,
918 )
919 source = cat.addNew()
920 source.set("centroid_x", center[0])
921 source.set("centroid_y", center[1])
922 offset = geom.Point2I(*center)
923 tmpSpans = afwGeom.SpanSet.fromShape(int(width), offset=offset)
924 source.setFootprint(afwDetection.Footprint(tmpSpans))
925 plugin.measure(source, exposure)
926 x = source.get("ext_shapeHSM_HsmPsfMoments_x")
927 y = source.get("ext_shapeHSM_HsmPsfMoments_y")
928 xx = source.get("ext_shapeHSM_HsmPsfMoments_xx")
929 yy = source.get("ext_shapeHSM_HsmPsfMoments_yy")
930 xy = source.get("ext_shapeHSM_HsmPsfMoments_xy")
932 self.assertFalse(source.get("ext_shapeHSM_HsmPsfMoments_flag"))
933 self.assertFalse(source.get("ext_shapeHSM_HsmPsfMoments_flag_no_pixels"))
934 self.assertFalse(source.get("ext_shapeHSM_HsmPsfMoments_flag_not_contained"))
935 self.assertFalse(source.get("ext_shapeHSM_HsmPsfMoments_flag_parent_source"))
937 if width < 4.5:
938 # i.e., as long as the PSF is not truncated for our 35x35 box.
939 self.assertAlmostEqual(x, 0.0, 3)
940 self.assertAlmostEqual(y, 0.0, 3)
941 expected = afwEll.Quadrupole(afwEll.Axes(width, width, 0.0))
942 self.assertAlmostEqual(xx, expected.getIxx(), SHAPE_DECIMALS)
943 self.assertAlmostEqual(xy, expected.getIxy(), SHAPE_DECIMALS)
944 self.assertAlmostEqual(yy, expected.getIyy(), SHAPE_DECIMALS)
946 # Test schema documentation
947 for fieldName in cat.schema.extract("*HsmPsfMoments_[xy]"):
948 self.assertEqual(
949 cat.schema[fieldName].asField().getDoc(), "Centroid of the PSF via the HSM shape algorithm"
950 )
951 for fieldName in cat.schema.extract("*HsmPsfMoments_[xy][xy]*"):
952 self.assertEqual(
953 cat.schema[fieldName].asField().getDoc(),
954 "Adaptive moments of the PSF via the HSM shape algorithm",
955 )
957 # Test that the moments are identical to those obtained directly by
958 # GalSim. For `width` > 4.5 where the truncation becomes significant,
959 # the answer might not be 'correct' but should remain 'consistent'.
960 xxDirect, yyDirect, xyDirect = self.computeDirectPsfMomentsFromGalSim(
961 psf,
962 geom.Point2D(*center),
963 useSourceCentroidOffset=useSourceCentroidOffset,
964 )
965 self.assertEqual(xx, xxDirect)
966 self.assertEqual(yy, yyDirect)
967 self.assertEqual(xy, xyDirect)
969 @lsst.utils.tests.methodParameters(
970 # Make Cartesian product of settings to feed to methodParameters
971 **dict(
972 list(
973 zip(
974 (
975 kwargs := dict(
976 width=(2.0, 3.0, 4.0),
977 useSourceCentroidOffset=(True, False),
978 varyBBox=(True, False),
979 wrongBBox=(True, False),
980 center=(
981 (23.0, 34.0), # various offsets that might cause trouble
982 (23.5, 34.0),
983 (23.5, 34.5),
984 (23.15, 34.25),
985 (22.81, 34.01),
986 (22.81, 33.99),
987 ),
988 )
989 ).keys(),
990 zip(*itertools.product(*kwargs.values())),
991 )
992 )
993 )
994 )
995 def testHsmPsfMomentsDebiased(self, width, useSourceCentroidOffset, varyBBox, wrongBBox, center):
996 # As a note, it's really hard to actually unit test whether we've
997 # succesfully "debiased" these measurements. That would require a
998 # many-object comparison of moments with and without noise. So we just
999 # test similar to the biased moments above.
1000 var = 1.2
1001 # As we reduce the flux, our deviation from the expected value
1002 # increases, so decrease tolerance.
1003 for flux, decimals in [
1004 (1e6, 3),
1005 (1e4, 1),
1006 (1e3, 0),
1007 ]:
1008 psf = PyGaussianPsf(35, 35, width, varyBBox=varyBBox, wrongBBox=wrongBBox)
1009 exposure = afwImage.ExposureF(45, 56)
1010 exposure.getMaskedImage().set(1.0, 0, var)
1011 exposure.setPsf(psf)
1013 algorithmName = "ext_shapeHSM_HsmPsfMomentsDebiased"
1014 alg = base.SingleFramePlugin.registry[algorithmName].PluginClass
1016 # perform the shape measurement
1017 control = lsst.meas.extensions.shapeHSM.HsmPsfMomentsDebiasedConfig()
1018 self.assertTrue(control.useSourceCentroidOffset)
1019 self.assertEqual(control.noiseSource, "variance")
1020 control.useSourceCentroidOffset = useSourceCentroidOffset
1021 plugin, cat = makePluginAndCat(
1022 alg,
1023 algorithmName,
1024 centroid="centroid",
1025 psfflux="base_PsfFlux",
1026 control=control,
1027 metadata=True,
1028 )
1029 source = cat.addNew()
1030 source.set("centroid_x", center[0])
1031 source.set("centroid_y", center[1])
1032 offset = geom.Point2I(*center)
1033 source.set("base_PsfFlux_instFlux", flux)
1034 tmpSpans = afwGeom.SpanSet.fromShape(int(width), offset=offset)
1035 source.setFootprint(afwDetection.Footprint(tmpSpans))
1037 plugin.measure(source, exposure)
1038 x = source.get("ext_shapeHSM_HsmPsfMomentsDebiased_x")
1039 y = source.get("ext_shapeHSM_HsmPsfMomentsDebiased_y")
1040 xx = source.get("ext_shapeHSM_HsmPsfMomentsDebiased_xx")
1041 yy = source.get("ext_shapeHSM_HsmPsfMomentsDebiased_yy")
1042 xy = source.get("ext_shapeHSM_HsmPsfMomentsDebiased_xy")
1043 for flag in [
1044 "ext_shapeHSM_HsmPsfMomentsDebiased_flag",
1045 "ext_shapeHSM_HsmPsfMomentsDebiased_flag_no_pixels",
1046 "ext_shapeHSM_HsmPsfMomentsDebiased_flag_not_contained",
1047 "ext_shapeHSM_HsmPsfMomentsDebiased_flag_parent_source",
1048 "ext_shapeHSM_HsmPsfMomentsDebiased_flag_edge",
1049 ]:
1050 self.assertFalse(source.get(flag))
1052 expected = afwEll.Quadrupole(afwEll.Axes(width, width, 0.0))
1054 self.assertAlmostEqual(x, 0.0, decimals)
1055 self.assertAlmostEqual(y, 0.0, decimals)
1057 T = expected.getIxx() + expected.getIyy()
1058 self.assertAlmostEqual((xx - expected.getIxx()) / T, 0.0, decimals)
1059 self.assertAlmostEqual((xy - expected.getIxy()) / T, 0.0, decimals)
1060 self.assertAlmostEqual((yy - expected.getIyy()) / T, 0.0, decimals)
1062 # Repeat using noiseSource='meta'. Should get nearly the same
1063 # results if BGMEAN is set to `var` above.
1064 exposure2 = afwImage.ExposureF(45, 56)
1065 # set the variance plane to something else to ensure we're
1066 # ignoring it
1067 exposure2.getMaskedImage().set(1.0, 0, 2 * var + 1.1)
1068 exposure2.setPsf(psf)
1069 exposure2.getMetadata().set("BGMEAN", var)
1071 control2 = shapeHSM.HsmPsfMomentsDebiasedConfig()
1072 control2.noiseSource = "meta"
1073 control2.useSourceCentroidOffset = useSourceCentroidOffset
1074 plugin2, cat2 = makePluginAndCat(
1075 alg,
1076 algorithmName,
1077 centroid="centroid",
1078 psfflux="base_PsfFlux",
1079 control=control2,
1080 metadata=True,
1081 )
1082 source2 = cat2.addNew()
1083 source2.set("centroid_x", center[0])
1084 source2.set("centroid_y", center[1])
1085 offset2 = geom.Point2I(*center)
1086 source2.set("base_PsfFlux_instFlux", flux)
1087 tmpSpans2 = afwGeom.SpanSet.fromShape(int(width), offset=offset2)
1088 source2.setFootprint(afwDetection.Footprint(tmpSpans2))
1090 plugin2.measure(source2, exposure2)
1091 x2 = source2.get("ext_shapeHSM_HsmPsfMomentsDebiased_x")
1092 y2 = source2.get("ext_shapeHSM_HsmPsfMomentsDebiased_y")
1093 xx2 = source2.get("ext_shapeHSM_HsmPsfMomentsDebiased_xx")
1094 yy2 = source2.get("ext_shapeHSM_HsmPsfMomentsDebiased_yy")
1095 xy2 = source2.get("ext_shapeHSM_HsmPsfMomentsDebiased_xy")
1096 for flag in [
1097 "ext_shapeHSM_HsmPsfMomentsDebiased_flag",
1098 "ext_shapeHSM_HsmPsfMomentsDebiased_flag_no_pixels",
1099 "ext_shapeHSM_HsmPsfMomentsDebiased_flag_not_contained",
1100 "ext_shapeHSM_HsmPsfMomentsDebiased_flag_parent_source",
1101 "ext_shapeHSM_HsmPsfMomentsDebiased_flag_edge",
1102 ]:
1103 self.assertFalse(source.get(flag))
1105 # Would be identically equal, but variance input via "BGMEAN" is
1106 # consumed in c++ as a double, where variance from the variance
1107 # plane is a c++ float.
1108 self.assertAlmostEqual(x, x2, 8)
1109 self.assertAlmostEqual(y, y2, 8)
1110 self.assertAlmostEqual(xx, xx2, 5)
1111 self.assertAlmostEqual(xy, xy2, 5)
1112 self.assertAlmostEqual(yy, yy2, 5)
1114 # Test schema documentation
1115 for fieldName in cat.schema.extract("*HsmPsfMoments_[xy]"):
1116 self.assertEqual(
1117 cat.schema[fieldName].asField().getDoc(),
1118 "Debiased centroid of the PSF via the HSM shape algorithm",
1119 )
1120 for fieldName in cat.schema.extract("*HsmPsfMoments_[xy][xy]*"):
1121 self.assertEqual(
1122 cat.schema[fieldName].asField().getDoc(),
1123 "Debiased adaptive moments of the PSF via the HSM shape algorithm",
1124 )
1126 testHsmPsfMomentsDebiasedEdgeArgs = dict(
1127 width=(2.0, 3.0, 4.0), useSourceCentroidOffset=(True, False), center=((1.2, 1.3), (33.2, 50.1))
1128 )
1130 @lsst.utils.tests.methodParameters(
1131 # Make Cartesian product of settings to feed to methodParameters
1132 **dict(
1133 list(
1134 zip(
1135 (
1136 kwargs := dict(
1137 width=(2.0, 3.0, 4.0),
1138 useSourceCentroidOffset=(True, False),
1139 center=[(1.2, 1.3), (33.2, 50.1)],
1140 )
1141 ).keys(),
1142 zip(*itertools.product(*kwargs.values())),
1143 )
1144 )
1145 )
1146 )
1147 def testHsmPsfMomentsDebiasedEdge(self, width, useSourceCentroidOffset, center):
1148 # As we reduce the flux, our deviation from the expected value
1149 # increases, so decrease tolerance.
1150 var = 1.2
1151 for flux, decimals in [
1152 (1e6, 3),
1153 (1e4, 2),
1154 (1e3, 1),
1155 ]:
1156 psf = PyGaussianPsf(35, 35, width)
1157 exposure = afwImage.ExposureF(45, 56)
1158 exposure.getMaskedImage().set(1.0, 0, 2 * var + 1.1)
1159 exposure.setPsf(psf)
1161 algorithmName = "ext_shapeHSM_HsmPsfMomentsDebiased"
1162 alg = base.SingleFramePlugin.registry[algorithmName].PluginClass
1164 # perform the shape measurement
1165 control = shapeHSM.HsmPsfMomentsDebiasedConfig()
1166 control.useSourceCentroidOffset = useSourceCentroidOffset
1167 self.assertEqual(control.noiseSource, "variance")
1168 plugin, cat = makePluginAndCat(
1169 alg,
1170 algorithmName,
1171 centroid="centroid",
1172 psfflux="base_PsfFlux",
1173 control=control,
1174 metadata=True,
1175 )
1176 source = cat.addNew()
1177 source.set("centroid_x", center[0])
1178 source.set("centroid_y", center[1])
1179 offset = geom.Point2I(*center)
1180 source.set("base_PsfFlux_instFlux", flux)
1181 tmpSpans = afwGeom.SpanSet.fromShape(int(width), offset=offset)
1182 source.setFootprint(afwDetection.Footprint(tmpSpans))
1184 # Edge fails when setting noise from var plane
1185 with self.assertRaises(base.MeasurementError):
1186 plugin.measure(source, exposure)
1188 # Succeeds when noise is from meta
1189 exposure.getMetadata().set("BGMEAN", var)
1190 control.noiseSource = "meta"
1191 plugin, cat = makePluginAndCat(
1192 alg,
1193 algorithmName,
1194 centroid="centroid",
1195 psfflux="base_PsfFlux",
1196 control=control,
1197 metadata=True,
1198 )
1199 source = cat.addNew()
1200 source.set("centroid_x", center[0])
1201 source.set("centroid_y", center[1])
1202 offset = geom.Point2I(*center)
1203 source.set("base_PsfFlux_instFlux", flux)
1204 tmpSpans = afwGeom.SpanSet.fromShape(int(width), offset=offset)
1205 source.setFootprint(afwDetection.Footprint(tmpSpans))
1206 plugin.measure(source, exposure)
1208 x = source.get("ext_shapeHSM_HsmPsfMomentsDebiased_x")
1209 y = source.get("ext_shapeHSM_HsmPsfMomentsDebiased_y")
1210 xx = source.get("ext_shapeHSM_HsmPsfMomentsDebiased_xx")
1211 yy = source.get("ext_shapeHSM_HsmPsfMomentsDebiased_yy")
1212 xy = source.get("ext_shapeHSM_HsmPsfMomentsDebiased_xy")
1213 self.assertFalse(source.get("ext_shapeHSM_HsmPsfMomentsDebiased_flag"))
1214 self.assertFalse(source.get("ext_shapeHSM_HsmPsfMomentsDebiased_flag_no_pixels"))
1215 self.assertFalse(source.get("ext_shapeHSM_HsmPsfMomentsDebiased_flag_not_contained"))
1216 self.assertFalse(source.get("ext_shapeHSM_HsmPsfMomentsDebiased_flag_parent_source"))
1217 # but _does_ set EDGE flag in this case
1218 self.assertTrue(source.get("ext_shapeHSM_HsmPsfMomentsDebiased_flag_edge"))
1220 expected = afwEll.Quadrupole(afwEll.Axes(width, width, 0.0))
1222 self.assertAlmostEqual(x, 0.0, decimals)
1223 self.assertAlmostEqual(y, 0.0, decimals)
1225 T = expected.getIxx() + expected.getIyy()
1226 self.assertAlmostEqual((xx - expected.getIxx()) / T, 0.0, decimals)
1227 self.assertAlmostEqual((xy - expected.getIxy()) / T, 0.0, decimals)
1228 self.assertAlmostEqual((yy - expected.getIyy()) / T, 0.0, decimals)
1230 # But fails hard if meta doesn't contain BGMEAN
1231 exposure.getMetadata().remove("BGMEAN")
1232 plugin, cat = makePluginAndCat(
1233 alg,
1234 algorithmName,
1235 centroid="centroid",
1236 psfflux="base_PsfFlux",
1237 control=control,
1238 metadata=True,
1239 )
1240 source = cat.addNew()
1241 source.set("centroid_x", center[0])
1242 source.set("centroid_y", center[1])
1243 offset = geom.Point2I(*center)
1244 source.set("base_PsfFlux_instFlux", flux)
1245 tmpSpans = afwGeom.SpanSet.fromShape(int(width), offset=offset)
1246 source.setFootprint(afwDetection.Footprint(tmpSpans))
1247 with self.assertRaises(base.FatalAlgorithmError):
1248 plugin.measure(source, exposure)
1250 def testHsmPsfMomentsDebiasedBadNoiseSource(self):
1251 control = shapeHSM.HsmPsfMomentsDebiasedConfig()
1252 with self.assertRaises(pexConfig.FieldValidationError):
1253 control.noiseSource = "ACM"
1256class TestMemory(lsst.utils.tests.MemoryTestCase):
1257 pass
1260def setup_module(module):
1261 lsst.utils.tests.init()
1264if __name__ == "__main__": 1264 ↛ 1265line 1264 didn't jump to line 1265 because the condition on line 1264 was never true
1265 lsst.utils.tests.init()
1266 unittest.main()