PSF fourth-order moments in the sky frame (spin-2 + rho4)

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,

src/shapepipe/modules/psfex_interp_package/psfex_interp.py

@@ -2,14 +2,15 @@
 
 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
 
 """
 
 import os
 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