From 552204de1dde81ad8a82983e362c3cab108cae58 Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Tue, 21 Jul 2026 14:37:33 +0200 Subject: [PATCH 1/4] feat(psfex-interp): add spin-2 fourth moments and rho4 in sky frame Extend the sky-coordinate HSM measurement (PR #812) to also store fourth-order moments for PSF and star: - spin-2 combinations HSM_M4_1_* = M40 - M04 and HSM_M4_2_* = 2*(M13 + M31), ported from the PSFHOME/moment4 math, and - galsim's spin-0 HSM_RHO4_* (moments_rho4). Frame handling is the key change over the moment4 port. Whitening with the symmetric sqrt(inv(M)) leaves the whitened axes aligned with the frame in which M and the pixel grid are built, so the spin-2 combinations are defined relative to those axes. Building them in the pixel frame ties them to the CCD orientation. Instead _fourth_moments builds everything in the world frame: M from the use_sky_coords world ellipticity/size, and the grid via the local WCS Jacobian mapping pixel offsets to world offsets about the world centroid (moments_centroid is relative to the stamp true_center). The result is invariant to a re-orientation of the pixel frame viewing the same sky object. HSM failures fill the spin-2 terms with 0 and pass rho4 through (-1); the FLAG column stays the source of truth. Masked star pixels are zeroed before the moment sum. Columns are added to the single-epoch, validation and multi-epoch output paths. Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_01K98geWHTuT62whqjKocfCm --- .../psfex_interp_package/psfex_interp.py | 157 +++++++++++++++++- 1 file changed, 154 insertions(+), 3 deletions(-) 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 From 351b985dc14f2496f388cd53569a7dcd38767eed Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Tue, 21 Jul 2026 14:37:33 +0200 Subject: [PATCH 2/4] test(psfex-interp): analytic and frame-invariance tests for fourth moments Drive the real _get_psfshapes / _get_starshapes on galsim stamps rendered through a nontrivial local WCS (scale + rotation + shear): - round and elliptical Gaussians: whitened spin-2 fourth moments ~ 0 and rho4 ~ 2 (whitening circularises any single sheared circular profile, so an elliptical Moffat would also vanish -- a two-Gaussian composite is used to get a genuine non-zero spin-2 signal), - frame invariance: one sky object rendered under WCS orientations of 0/28/63/-45/90 deg gives the same world-frame fourth moments, on both the PSF and star paths, - PSF/star agreement on a clean stamp, and HSM-failure sentinel handling. Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_01K98geWHTuT62whqjKocfCm --- tests/module/test_psf_fourth_moments.py | 189 ++++++++++++++++++++++++ 1 file changed, 189 insertions(+) create mode 100644 tests/module/test_psf_fourth_moments.py diff --git a/tests/module/test_psf_fourth_moments.py b/tests/module/test_psf_fourth_moments.py new file mode 100644 index 000000000..0cd2a82bd --- /dev/null +++ b/tests/module/test_psf_fourth_moments.py @@ -0,0 +1,189 @@ +"""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-3, atol=1e-5) + npt.assert_allclose(rot[M4_2], ref[M4_2], rtol=1e-3, atol=1e-5) + npt.assert_allclose(rot[RHO4], ref[RHO4], rtol=1e-4) + + +@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-3, atol=1e-5) + npt.assert_allclose(rot[M4_2], ref[M4_2], rtol=1e-3, atol=1e-5) + npt.assert_allclose(rot[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) + + +# --------------------------------------------------------------------------- +# 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, + ) From e5b7f5d25bba7a7d823ba883e160e5e8afe10fd4 Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Tue, 21 Jul 2026 14:53:38 +0200 Subject: [PATCH 3/4] fix(merge-starcat): carry HSM fourth moments through PSFEx merge MergeStarCatPSFEX read and wrote only HSM_G1/G2/T/FLAG for PSF and STAR, so the spin-2 fourth moments and rho4 written per-tile by psfex_interp were dropped at the merge step. The merged full_starcat is the catalogue PSF validation consumes, so the feature was inert past the per-tile stage. Read and write HSM_M4_1/HSM_M4_2/HSM_RHO4 for both PSF and STAR, matching the per-tile column grammar. Avoids the moment4 reference's copy-paste bug (m4_2_psf now reads the _2 column) and adds RHO4, which the reference omitted. Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_01K98geWHTuT62whqjKocfCm --- .../modules/merge_starcat_package/merge_starcat.py | 14 ++++++++++++++ 1 file changed, 14 insertions(+) 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, From b07f0cd43a15c914333f7f139259004e75f1a787 Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Tue, 21 Jul 2026 14:53:52 +0200 Subject: [PATCH 4/4] test(psfex-interp): cover bad-pixel masking and centroid transform Two paths the summary claimed but the suite left unexercised: - Bad-pixel masking: every star fixture had an all-zero mask, so zeroing masked pixels before the fourth-moment sum was a no-op in the suite. Add a test that plants the -1e30 sentinel in outskirt pixels, asserts the masked vignet reproduces the clean result, and asserts the same sentinels fed through un-zeroed blow the fourth moments up by ~5 orders (mutation: removing the zeroing makes M4_1 go -0.03 -> -4278). - Centroid transform: every fixture was drawn centered, so moments_centroid ~ 0 and the sky-frame recentering was dead. Add a world-frame-shift case; the fourth moments stay invariant only because the transformed centroid is subtracted (mutation: dropping the subtraction fails both angles). Also tighten the frame-invariance tolerances from rtol=1e-3/atol=1e-5 to rtol=1e-5/atol=1e-6 (measured cross-orientation agreement is ~1e-8), so a subtle partial-frame error at the 1e-4 level now bites. Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_01K98geWHTuT62whqjKocfCm --- tests/module/test_psf_fourth_moments.py | 78 +++++++++++++++++++++++-- 1 file changed, 72 insertions(+), 6 deletions(-) diff --git a/tests/module/test_psf_fourth_moments.py b/tests/module/test_psf_fourth_moments.py index 0cd2a82bd..4782ca5f4 100644 --- a/tests/module/test_psf_fourth_moments.py +++ b/tests/module/test_psf_fourth_moments.py @@ -134,9 +134,9 @@ def test_psf_fourth_moments_frame_invariant(theta): 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-3, atol=1e-5) - npt.assert_allclose(rot[M4_2], ref[M4_2], rtol=1e-3, atol=1e-5) - npt.assert_allclose(rot[RHO4], ref[RHO4], rtol=1e-4) + 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]) @@ -146,9 +146,31 @@ def test_star_fourth_moments_frame_invariant(theta): 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-3, atol=1e-5) - npt.assert_allclose(rot[M4_2], ref[M4_2], rtol=1e-3, atol=1e-5) - npt.assert_allclose(rot[RHO4], ref[RHO4], rtol=1e-4) + 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(): @@ -162,6 +184,50 @@ def test_psf_and_star_paths_agree(): 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. # ---------------------------------------------------------------------------