diff --git a/src/shapepipe/modules/merge_starcat_package/merge_starcat.py b/src/shapepipe/modules/merge_starcat_package/merge_starcat.py index 7bc0cb76b..58fdd3767 100644 --- a/src/shapepipe/modules/merge_starcat_package/merge_starcat.py +++ b/src/shapepipe/modules/merge_starcat_package/merge_starcat.py @@ -559,7 +559,9 @@ def process(self): """ x, y, ra, dec = [], [], [], [] g1_psf, g2_psf, size_psf = [], [], [] + m4_1_psf, m4_2_psf, rho4_psf = [], [], [] g1, g2, size = [], [], [] + m4_1_star, m4_2_star, rho4_star = [], [], [] flag_psf, flag_star = [], [] mag, snr, psfex_acc = [], [], [] ccd_nb = [] @@ -588,9 +590,15 @@ def process(self): g1_psf += list(data_j["HSM_G1_PSF"]) g2_psf += list(data_j["HSM_G2_PSF"]) size_psf += list(data_j["HSM_T_PSF"]) + m4_1_psf += list(data_j["HSM_M4_1_PSF"]) + m4_2_psf += list(data_j["HSM_M4_2_PSF"]) + rho4_psf += list(data_j["HSM_RHO4_PSF"]) g1 += list(data_j["HSM_G1_STAR"]) g2 += list(data_j["HSM_G2_STAR"]) size += list(data_j["HSM_T_STAR"]) + m4_1_star += list(data_j["HSM_M4_1_STAR"]) + m4_2_star += list(data_j["HSM_M4_2_STAR"]) + rho4_star += list(data_j["HSM_RHO4_STAR"]) # flags flag_psf += list(data_j["HSM_FLAG_PSF"]) @@ -636,9 +644,15 @@ def process(self): "HSM_G1_PSF": g1_psf, "HSM_G2_PSF": g2_psf, "HSM_T_PSF": size_psf, + "HSM_M4_1_PSF": m4_1_psf, + "HSM_M4_2_PSF": m4_2_psf, + "HSM_RHO4_PSF": rho4_psf, "HSM_G1_STAR": g1, "HSM_G2_STAR": g2, "HSM_T_STAR": size, + "HSM_M4_1_STAR": m4_1_star, + "HSM_M4_2_STAR": m4_2_star, + "HSM_RHO4_STAR": rho4_star, "HSM_FLAG_PSF": flag_psf, "HSM_FLAG_STAR": flag_star, "MAG": mag, diff --git a/src/shapepipe/modules/psfex_interp_package/psfex_interp.py b/src/shapepipe/modules/psfex_interp_package/psfex_interp.py index b6da6623f..044ed0092 100644 --- a/src/shapepipe/modules/psfex_interp_package/psfex_interp.py +++ b/src/shapepipe/modules/psfex_interp_package/psfex_interp.py @@ -2,7 +2,7 @@ This module computes the PSFs from a PSFEx model at several galaxy positions. -:Authors: Morgan Schmitz, Axel Guinot, Martin Kilbinger +:Authors: Morgan Schmitz, Axel Guinot, Martin Kilbinger, Sacha Guerrini """ @@ -10,6 +10,7 @@ import re import numpy as np +import scipy.linalg as alg from astropy.io import fits from cs_util import size as cs_size from sqlitedict import SqliteDict @@ -62,6 +63,117 @@ def local_wcs_list(wcs, positions): ] +def _fourth_moments(image, moms, wcs=None): + r"""Fourth-Order Moments. + + Compute the spin-2 fourth-moment combinations and galsim's spin-0 + ``moments_rho4`` for a single object, given its HSM adaptive-moment result. + + The spin-2 combinations follow the PSFHOME convention: whiten the pixel grid + with the object's own second-moment matrix ``M`` (so a matched elliptical + Gaussian becomes a unit circle), Gaussian-weight it, form the centred + ``p + q = 4`` moments ``M_pq`` over the whitened coordinates ``(u, v)``, and + return :: + + fourth_moment_1 = M_40 - M_04 + fourth_moment_2 = 2 * (M_13 + M_31) + + **Frame handling.** Whitening with the *symmetric* matrix square root + ``sqrt(inv(M))`` leaves the whitened axes aligned with the coordinate frame + in which ``M`` and the grid are built (it rescales along the object's + principal axes but does not rotate the frame). The spin-2 combinations are + therefore measured relative to *those* axes. To make them a property of the + sky and not of the CCD orientation, everything must live in one frame — the + world frame: + + * ``moms`` comes from ``FindAdaptiveMom(use_sky_coords=True)``, so + ``observed_shape``, ``moments_sigma`` (arcsec) and ``moments_centroid`` + are already in world ``(u, v)`` coordinates. ``M`` built from them is the + world-frame second-moment matrix. + * The pixel grid is mapped to world offsets about the centroid with the local + WCS Jacobian ``J``: ``world = J @ (pixel - true_center)``. galsim reports + ``moments_centroid`` relative to the stamp ``true_center``, so subtracting + it gives ``J @ (pixel - pixel_centroid)`` — the world offset from the + centroid, in the same arcsec units as ``M``. + + The result is invariant under a re-orientation of the pixel frame (a + different WCS rotation viewing the same sky object) — the science guarantee + that motivates the sky-coordinate measurement. When ``wcs`` is ``None`` the + computation is done consistently in the pixel frame instead (the axes are + then the CCD axes and the result is *not* frame-invariant); this path exists + only for the legacy pixel-frame branch. + + Parameters + ---------- + image : numpy.ndarray + Object postage stamp (PSF or star vignet), masked pixels already zeroed. + moms : galsim.hsm.ShapeData + Adaptive-moment result for ``image``. Measured with + ``use_sky_coords=True`` when ``wcs`` is provided. + wcs : galsim.JacobianWCS, optional + Local WCS at the object position, matching the one attached to the image + for the ``use_sky_coords`` measurement. ``None`` selects the pixel frame. + + Returns + ------- + tuple of float + ``(fourth_moment_1, fourth_moment_2, moments_rho4)``. On an HSM failure + (``moms.error_message`` set) the spin-2 terms are filled with ``0.0`` and + ``moments_rho4`` is passed through (``-1`` on failure); the accompanying + ``FLAG`` column is the source of truth for validity. + + """ + if moms.error_message: + return 0.0, 0.0, moms.moments_rho4 + + ny, nx = image.shape + # 1-indexed pixel-centre coordinates (galsim Image origin is (1, 1)). + y_grid, x_grid = np.mgrid[:ny, :nx] + 1.0 + + # Second-moment matrix M with det(M) = sigma**4, built in the measurement + # frame from the (distortion) ellipticity and size. In sky mode these are + # world-frame quantities, so M is the world-frame second-moment matrix. + e1 = moms.observed_shape.e1 + e2 = moms.observed_shape.e2 + sigma4 = moms.moments_sigma**4 + c = (1 + e1) / (1 - e1) + M = np.zeros((2, 2)) + M[1, 1] = np.sqrt(sigma4 / (c - 0.25 * e2**2 * (1 + c) ** 2)) + M[0, 0] = c * M[1, 1] + M[0, 1] = M[1, 0] = 0.5 * e2 * (M[0, 0] + M[1, 1]) + + if wcs is not None: + # Map the pixel grid to world offsets about the world centroid. + x0 = 0.5 * (nx + 1) + y0 = 0.5 * (ny + 1) + pix = np.array([x_grid - x0, y_grid - y0]) + jac = wcs.jacobian().getMatrix() # [[dudx, dudy], [dvdx, dvdy]] + world = np.einsum("ij,jqp->iqp", jac, pix) + u = world[0] - moms.moments_centroid.x + v = world[1] - moms.moments_centroid.y + else: + # Pixel frame: offsets straight from the pixel-frame centroid. + u = x_grid - moms.moments_centroid.x + v = y_grid - moms.moments_centroid.y + + # Whiten with the symmetric square root of inv(M); M is SPD so the principal + # square root is real (drop the ~1e-16 imaginary residue sqrtm returns). + sqrt_inv_M = np.real(alg.sqrtm(np.linalg.inv(M))) + std_pos = np.einsum("ij,jqp->iqp", sqrt_inv_M, np.array([u, v])) + std_x, std_y = std_pos[0], std_pos[1] + + weight = np.exp(-0.5 * (std_x**2 + std_y**2)) + image_weight = weight * image + normalization = np.sum(image_weight) + + def _m(p, q): + return np.sum(image_weight * std_x**p * std_y**q) / normalization + + fourth_moment_1 = _m(4, 0) - _m(0, 4) + fourth_moment_2 = 2 * (_m(1, 3) + _m(3, 1)) + return fourth_moment_1, fourth_moment_2, moms.moments_rho4 + + class PSFExInterpolator(object): """The PSFEx Interpolator Class. @@ -355,6 +467,14 @@ def _get_psfshapes(self, wcs_list=None): returned are already in sky coordinates and need no post-hoc rotation. When ``None``, moments are measured in the pixel frame. + Notes + ----- + In addition to the second-moment shape, the spin-2 fourth-moment + combinations and galsim's spin-0 ``moments_rho4`` are stored (columns 4, + 5, 6 of ``psf_shapes``). With a ``wcs_list`` these are computed in the + world frame so they are invariant to the CCD orientation; see + :func:`_fourth_moments`. + """ if import_fail: raise ImportError("Galsim is required to get shapes information") @@ -371,6 +491,7 @@ def _get_psfshapes(self, wcs_list=None): hsm.FindAdaptiveMom(Image(psf), strict=False) for psf in self.interp_PSFs ] + wcs_list = [None] * len(self.interp_PSFs) self.psf_shapes = np.array( [ @@ -379,8 +500,9 @@ def _get_psfshapes(self, wcs_list=None): moms.observed_shape.g2, moms.moments_sigma, int(bool(moms.error_message)), + *_fourth_moments(psf, moms, wcs), ] - for moms in psf_moms + for psf, moms, wcs in zip(self.interp_PSFs, psf_moms, wcs_list) ] ) @@ -403,6 +525,9 @@ def _write_output(self): "HSM_G2_PSF": self.psf_shapes[:, 1], "HSM_T_PSF": cs_size.sigma_to_T(self.psf_shapes[:, 2]), "HSM_FLAG_PSF": self.psf_shapes[:, 3].astype(int), + "HSM_M4_1_PSF": self.psf_shapes[:, 4], + "HSM_M4_2_PSF": self.psf_shapes[:, 5], + "HSM_RHO4_PSF": self.psf_shapes[:, 6], } else: data = {"VIGNET": self.interp_PSFs} @@ -484,6 +609,13 @@ def _get_starshapes(self, star_vign, wcs_list=None): provided, moments are measured in world coordinates (``use_sky_coords=True``). + Notes + ----- + As in :meth:`_get_psfshapes`, the spin-2 fourth-moment combinations and + galsim's ``moments_rho4`` are stored (columns 4, 5, 6 of ``star_shapes``), + world-frame when a ``wcs_list`` is given. Masked pixels are zeroed before + the fourth-moment sum so they do not contribute. + """ if import_fail: raise ImportError("Galsim is required to get shapes information") @@ -508,6 +640,7 @@ def _get_starshapes(self, star_vign, wcs_list=None): ) for star, mask in zip(star_vign, masks) ] + wcs_list = [None] * len(star_vign) self.star_shapes = np.array( [ @@ -516,8 +649,11 @@ def _get_starshapes(self, star_vign, wcs_list=None): moms.observed_shape.g2, moms.moments_sigma, int(bool(moms.error_message)), + *_fourth_moments(np.where(mask == 1, 0.0, star), moms, wcs), ] - for moms in star_moms + for star, mask, moms, wcs in zip( + star_vign, masks, star_moms, wcs_list + ) ] ) @@ -598,10 +734,16 @@ def _write_output_validation(self, star_dict, psfex_cat_dict): "HSM_G2_PSF": self.psf_shapes[:, 1], "HSM_T_PSF": cs_size.sigma_to_T(self.psf_shapes[:, 2]), "HSM_FLAG_PSF": self.psf_shapes[:, 3].astype(int), + "HSM_M4_1_PSF": self.psf_shapes[:, 4], + "HSM_M4_2_PSF": self.psf_shapes[:, 5], + "HSM_RHO4_PSF": self.psf_shapes[:, 6], "HSM_G1_STAR": self.star_shapes[:, 0], "HSM_G2_STAR": self.star_shapes[:, 1], "HSM_T_STAR": cs_size.sigma_to_T(self.star_shapes[:, 2]), "HSM_FLAG_STAR": self.star_shapes[:, 3].astype(int), + "HSM_M4_1_STAR": self.star_shapes[:, 4], + "HSM_M4_2_STAR": self.star_shapes[:, 5], + "HSM_RHO4_STAR": self.star_shapes[:, 6], } data = {**data, **star_dict} @@ -874,6 +1016,15 @@ def _interpolate_me(self): shape_dict["HSM_FLAG_PSF"] = final_list[j][2][ where_res[0] ][3] + shape_dict["HSM_M4_1_PSF"] = final_list[j][2][ + where_res[0] + ][4] + shape_dict["HSM_M4_2_PSF"] = final_list[j][2][ + where_res[0] + ][5] + shape_dict["HSM_RHO4_PSF"] = final_list[j][2][ + where_res[0] + ][6] output_dict[id_tmp][final_list[j][3][where_res[0]]][ "SHAPES" ] = shape_dict diff --git a/tests/module/test_psf_fourth_moments.py b/tests/module/test_psf_fourth_moments.py new file mode 100644 index 000000000..4782ca5f4 --- /dev/null +++ b/tests/module/test_psf_fourth_moments.py @@ -0,0 +1,255 @@ +"""UNIT TESTS FOR PSF/STAR FOURTH-ORDER MOMENTS. + +Cover the spin-2 fourth-moment combinations (``HSM_M4_1_*`` / ``HSM_M4_2_*``) +and galsim's spin-0 ``moments_rho4`` (``HSM_RHO4_*``) added to the PSFEx-interp +sky-coordinate shape measurement (:mod:`shapepipe.modules.psfex_interp_package`). + +The measurement whitens the object by its own second-moment matrix, so the +predictions are analytic for elliptically-symmetric profiles: + +* a round or elliptical **Gaussian** matches its own adaptive weight, so the + whitened profile is an isotropic Gaussian: the spin-2 fourth moments vanish + and ``rho4`` equals the Gaussian value 2. +* to get *non-zero* spin-2 fourth moments the profile must not be a single + sheared circular profile (whitening circularises any such profile). We use a + sum of two coaxial Gaussians of different ellipticity, whose isophote shape + changes with radius. + +The science guarantee is **frame invariance**: because the whitening and the +grid are built in the world frame (see ``_fourth_moments``), one sky object +rendered onto two differently oriented pixel grids gives the same fourth +moments. That is the key test below. + +The WCS construction mirrors ``test_hsm_sky_coords.py``: a CD-matrix-like local +Jacobian with a realistic scale, rotation and small shear. +""" + +import galsim +import numpy as np +import numpy.testing as npt +import pytest + +from shapepipe.modules.psfex_interp_package.psfex_interp import ( + PSFExInterpolator, + _fourth_moments, +) + +# Column layout of psf_shapes / star_shapes. +M4_1, M4_2, RHO4 = 4, 5, 6 + +_STAMP = 101 + + +def make_wcs(theta_deg, scale=0.187, g1=0.06, g2=-0.04): + """A CD-matrix-like local WCS: scale, rotation ``theta_deg``, small shear.""" + th = np.deg2rad(theta_deg) + c, s = np.cos(th), np.sin(th) + mat = galsim.Shear(g1=g1, g2=g2).getMatrix() * scale @ np.array( + [[c, -s], [s, c]] + ) + return galsim.JacobianWCS(mat[0, 0], mat[0, 1], mat[1, 0], mat[1, 1]) + + +def render(profile, wcs): + """Render a world-frame profile onto the pixel grid defined by ``wcs``.""" + return profile.drawImage( + nx=_STAMP, ny=_STAMP, wcs=wcs, method="no_pixel" + ).array + + +def psf_row(profile, wcs): + """Run the real ``_get_psfshapes`` on one stamp; return its shape row.""" + interp = object.__new__(PSFExInterpolator) + interp.interp_PSFs = [render(profile, wcs)] + interp._get_psfshapes([wcs]) + return interp.psf_shapes[0] + + +def star_row(profile, wcs, mask=None): + """Run the real ``_get_starshapes`` on one star vignet; return its row.""" + stamp = render(profile, wcs) + if mask is not None: + stamp = np.where(mask, -1e30, stamp) + interp = object.__new__(PSFExInterpolator) + interp._get_starshapes(np.array([stamp]), [wcs]) + return interp.star_shapes[0] + + +# Sum of two coaxial Gaussians of different ellipticity: not elliptically +# symmetric, so its whitened spin-2 fourth moments are non-zero. +def composite(beta_deg=0.0): + prof = galsim.Gaussian(sigma=0.4).shear(e1=0.5) + galsim.Gaussian( + sigma=1.0 + ).shear(e1=0.1) + if beta_deg: + prof = prof.rotate(beta_deg * galsim.degrees) + return prof + + +# --------------------------------------------------------------------------- +# Analytic checks on Gaussians. +# --------------------------------------------------------------------------- + + +def test_round_gaussian_moments_vanish_rho4_is_gaussian(): + """Round Gaussian: spin-2 fourth moments ~ 0, rho4 ~ 2 (Gaussian value).""" + row = psf_row(galsim.Gaussian(sigma=0.6), make_wcs(0.0)) + npt.assert_allclose(row[M4_1], 0.0, atol=1e-4) + npt.assert_allclose(row[M4_2], 0.0, atol=1e-4) + npt.assert_allclose(row[RHO4], 2.0, rtol=1e-3) + + +@pytest.mark.parametrize( + "e1, e2", + [(0.3, 0.0), (0.0, 0.25), (0.3, 0.15), (-0.2, -0.3)], +) +def test_elliptical_gaussian_whitening_vanishes(e1, e2): + """Elliptical Gaussian: whitening by its own 2nd moment circularises it, so + the spin-2 fourth moments vanish regardless of (e1, e2).""" + row = psf_row(galsim.Gaussian(sigma=0.6).shear(e1=e1, e2=e2), make_wcs(0.0)) + npt.assert_allclose(row[M4_1], 0.0, atol=1e-4) + npt.assert_allclose(row[M4_2], 0.0, atol=1e-4) + npt.assert_allclose(row[RHO4], 2.0, rtol=1e-3) + + +def test_composite_has_nonzero_spin2(): + """Coaxial two-Gaussian composite: non-Gaussian radial shape leaves a real + spin-2 fourth moment. Axis-aligned, so only M4_1 is excited (M4_2 ~ 0).""" + row = psf_row(composite(), make_wcs(0.0)) + assert abs(row[M4_1]) > 1e-2 + npt.assert_allclose(row[M4_2], 0.0, atol=1e-4) + + +# --------------------------------------------------------------------------- +# Frame invariance -- the science guarantee. +# --------------------------------------------------------------------------- + + +@pytest.mark.parametrize("theta", [28.0, 63.0, -45.0, 90.0]) +def test_psf_fourth_moments_frame_invariant(theta): + """One sky object rendered under two WCS orientations gives the same + world-frame fourth moments (PSF path).""" + prof = composite(beta_deg=30.0) # rotated so both M4_1 and M4_2 are excited + ref = psf_row(prof, make_wcs(0.0)) + rot = psf_row(prof, make_wcs(theta)) + + assert abs(ref[M4_1]) > 1e-2 and abs(ref[M4_2]) > 1e-2 # non-trivial + npt.assert_allclose(rot[M4_1], ref[M4_1], rtol=1e-5, atol=1e-6) + npt.assert_allclose(rot[M4_2], ref[M4_2], rtol=1e-5, atol=1e-6) + npt.assert_allclose(rot[RHO4], ref[RHO4], rtol=1e-5) + + +@pytest.mark.parametrize("theta", [28.0, -45.0]) +def test_star_fourth_moments_frame_invariant(theta): + """Same guarantee on the star path (which also handles a bad-pixel mask).""" + prof = composite(beta_deg=30.0) + ref = star_row(prof, make_wcs(0.0)) + rot = star_row(prof, make_wcs(theta)) + + npt.assert_allclose(rot[M4_1], ref[M4_1], rtol=1e-5, atol=1e-6) + npt.assert_allclose(rot[M4_2], ref[M4_2], rtol=1e-5, atol=1e-6) + npt.assert_allclose(rot[RHO4], ref[RHO4], rtol=1e-5) + + +@pytest.mark.parametrize("theta", [0.0, 40.0]) +def test_fourth_moments_invariant_under_world_shift(theta): + """Off-center object: the sky-frame centroid subtraction makes the fourth + moments invariant under a world-frame translation of the source. + + Every other fixture is drawn centered, so ``moments_centroid`` is ~0 and the + recentering in ``_fourth_moments`` is a no-op. Real star vignets are only + approximately centered, so this case shifts the source in the *world* frame + (arcsec) and asserts the fourth moments are unchanged -- exercising the + centroid transform that is otherwise dead in the suite. + """ + prof = composite(beta_deg=30.0) + wcs = make_wcs(theta) + ref = psf_row(prof, wcs) + shifted = psf_row(prof.shift(0.9, -1.3), wcs) # world-frame shift, arcsec + + assert abs(ref[M4_1]) > 1e-2 and abs(ref[M4_2]) > 1e-2 # non-trivial + npt.assert_allclose(shifted[M4_1], ref[M4_1], rtol=1e-4, atol=1e-6) + npt.assert_allclose(shifted[M4_2], ref[M4_2], rtol=1e-4, atol=1e-6) + npt.assert_allclose(shifted[RHO4], ref[RHO4], rtol=1e-4) + + +def test_psf_and_star_paths_agree(): + """PSF and star measurement of the same clean stamp agree.""" + prof = composite(beta_deg=30.0) + wcs = make_wcs(17.0) + p = psf_row(prof, wcs) + s = star_row(prof, wcs) + npt.assert_allclose(s[M4_1], p[M4_1], rtol=1e-4, atol=1e-6) + npt.assert_allclose(s[M4_2], p[M4_2], rtol=1e-4, atol=1e-6) + npt.assert_allclose(s[RHO4], p[RHO4], rtol=1e-4) + + +# --------------------------------------------------------------------------- +# Bad-pixel masking on the star path. +# --------------------------------------------------------------------------- + + +def test_masked_pixels_are_zeroed_before_fourth_moment_sum(): + """Star path: pixels carrying the ``-1e30`` bad-pixel sentinel are zeroed + before the fourth-moment sum, so a masked vignet reproduces the clean + result -- and the zeroing is load-bearing. + + Without the zeroing the sum runs over the raw ``-1e30`` sentinels and the + fourth moments become garbage (the sum over an un-zeroed stamp is order + ``-1e30``). We mask low-flux outskirt pixels: dropping them leaves the clean + measurement essentially unchanged, while feeding the same sentinels through + un-zeroed blows the result up by many orders of magnitude. + """ + prof = composite(beta_deg=30.0) + wcs = make_wcs(17.0) + clean = star_row(prof, wcs) + + # A block of outskirt pixels, offset off-axis so it does not cancel in M4_1. + mask = np.zeros((_STAMP, _STAMP), dtype=bool) + c = _STAMP // 2 + mask[c + 18 : c + 22, c + 12 : c + 15] = True + + masked = star_row(prof, wcs, mask=mask) + + # Zeroing => the masked measurement reproduces the clean one. + npt.assert_allclose(masked[M4_1], clean[M4_1], rtol=1e-5, atol=1e-6) + npt.assert_allclose(masked[M4_2], clean[M4_2], rtol=1e-5, atol=1e-6) + npt.assert_allclose(masked[RHO4], clean[RHO4], rtol=1e-5) + + # Contrast: the same sentinel stamp with the zeroing removed is garbage. + sentinel = np.where(mask, -1e30, render(prof, wcs)) + moms = galsim.hsm.FindAdaptiveMom( + galsim.Image(sentinel, wcs=wcs), + badpix=galsim.Image(mask.astype(float)), + strict=False, + use_sky_coords=True, + ) + raw = _fourth_moments(sentinel, moms, wcs) # no zeroing applied + assert abs(raw[0] - clean[M4_1]) > 1.0 + + +# --------------------------------------------------------------------------- +# Failure handling. +# --------------------------------------------------------------------------- + + +def test_hsm_failure_fills_sentinels(): + """A flat stamp fails HSM: FLAG is set, spin-2 filled with 0, rho4 = -1.""" + interp = object.__new__(PSFExInterpolator) + interp.interp_PSFs = [np.zeros((21, 21))] + interp._get_psfshapes([make_wcs(0.0)]) + row = interp.psf_shapes[0] + assert int(row[3]) == 1 # FLAG + assert row[M4_1] == 0.0 + assert row[M4_2] == 0.0 + assert row[RHO4] == -1.0 + + +def test_fourth_moments_failure_guard(): + """``_fourth_moments`` short-circuits on a ShapeData carrying an error.""" + failed = galsim.hsm.ShapeData(error_message="boom", moments_rho4=-1.0) + assert _fourth_moments(np.zeros((5, 5)), failed, make_wcs(0.0)) == ( + 0.0, + 0.0, + -1.0, + )