diff --git a/docs/source/pipeline_canfar.md b/docs/source/pipeline_canfar.md index 7e8bb29c..96add575 100644 --- a/docs/source/pipeline_canfar.md +++ b/docs/source/pipeline_canfar.md @@ -333,21 +333,43 @@ shapepipe_run -c $SP_CONFIG/config_Ms_$psf.ini shapepipe_run -c $SP_CONFIG/config_Pl_$psf.ini ``` -#### Convert star catalogue to wCS +#### Collate star catalogues -Convert all input validation PSF files and create directories per patch `P?`. -Create files `validation_psf_conv--.fits` (for the v1.4 setup only one file): +Collate all input validation PSF files into star catalogues, gathering positions +(X/Y/RA/DEC) and the MCCD CCD id. + +Note: HSM shapes are no longer rotated into world coordinates at this step. The +PSF/star ellipticities and sizes are now measured directly in sky coordinates +during PSF interpolation (galsim `FindAdaptiveMom(use_sky_coords=True)`), so +`collate_star_cat.py` only collates and passes the shapes through. + +> **v2.0 is patch-less.** Runs up to `v1.6` are organised in sky patches +> `P1`..`P`, each patch a run directory. `v2.0` removes the patch concept: a +> single run root, outputs directly under it. `v2.0` is the default; select an +> older layout with `-V` (e.g. `-V v1.6`). For `v2.0` the patch loop and the +> `-P` option no longer apply, and the patch token drops from the output +> filename (`validation_psf_conv-.fits` instead of +> `validation_psf_conv--.fits`). ```bash cd /path/to/version -mkdir stat_car +mkdir star_cat cd star_cat ``` -For each patch run +For `v2.0` (the default), run once against the patch-less run root, producing +files `validation_psf_conv-.fits`: + +```bash +collate_star_cat.py -i .. -v +``` + +For `v1.x`, pass the version explicitly and run once per patch, creating a +directory per patch `P?` and producing files +`validation_psf_conv--.fits` (for the v1.4 setup only one file): ```bash -convert_psf_pix2world.py -i .. -P $patchnum -v +collate_star_cat.py -i .. -V v1.6 -P $patchnum -v ``` Combine previously created files as links within one ShapePipe run directory (for the v1.4 setup only one link). diff --git a/scripts/python/convert_psf_pix2world.py b/scripts/python/collate_star_cat.py similarity index 79% rename from scripts/python/convert_psf_pix2world.py rename to scripts/python/collate_star_cat.py index dea99656..f7ab25b4 100755 --- a/scripts/python/convert_psf_pix2world.py +++ b/scripts/python/collate_star_cat.py @@ -1,5 +1,32 @@ #! /usr/bin/env python3 +"""COLLATE STAR CATALOGUES. + +Collate the per-exposure PSF validation catalogues into star catalogues: gather +positions (X/Y/RA/DEC), assign the MCCD focal-plane CCD id, and merge into +``validation_psf_conv`` FITS files. + +Catalogue layout depends on the major version (``-V``). Runs up to ``v1.6`` are +organised in sky patches (``P1``..``P``): input runs live under +``/P/output`` and outputs are named +``validation_psf_conv--.fits``. ``v2.0`` removes the patch concept: +input runs live under a single ``/output`` root and outputs drop the patch +token (``validation_psf_conv-.fits``). + +HSM ellipticities and sizes are no longer rotated here. PSFEx and the in-repo +MCCD interpolation now measure adaptive moments directly in world coordinates +(galsim ``FindAdaptiveMom(use_sky_coords=True)``), so the WCS-Jacobian shape +rotation this script used to perform is redundant and has been removed. + +Caveat: the MCCD ``PSF_MOM_LIST``/``STAR_MOM_LIST`` columns are produced by the +external ``mccd`` fit-validation code (``mccd.auxiliary_fun.mccd_validation``), +which still measures HSM moments in the pixel frame. Those shapes are therefore +still rotated into world coordinates here, via the WCS-Jacobian rotation, until +``mccd`` itself adopts ``use_sky_coords``. This is the one branch that keeps the +rotation; the in-repo PSFEx and MCCD-interpolation paths measure adaptive +moments directly in world coordinates upstream and pass them straight through. +""" + import sys import os import re @@ -16,6 +43,58 @@ from cs_util import logging +def collate_paths(input_base_dir, output_base_dir, patch): + """Collate Paths. + + Return the ``(input run dir, output dir)`` for a patch. ``patch`` is None + for the patch-less v2.0 layout, which drops the ``P`` token; v1.x + passes the patch number. + + Parameters + ---------- + input_base_dir : str + input base directory + output_base_dir : str + output base directory + patch : str or None + patch number, or None for the patch-less v2.0 layout + + Returns + ------- + tuple + input run directory and output directory + + """ + if patch is None: + return f"{input_base_dir}/output/", output_base_dir + return f"{input_base_dir}/P{patch}/output/", f"{output_base_dir}/P{patch}" + + +def output_filename(file_pattern, patch, idx): + """Output Filename. + + Build the collated catalogue filename. ``patch`` is None for the patch-less + v2.0 layout, which drops the patch token. + + Parameters + ---------- + file_pattern : str + input file pattern (e.g. ``validation_psf``) + patch : str or None + patch number, or None for the patch-less v2.0 layout + idx : int + exposure run index + + Returns + ------- + str + output catalogue file name + + """ + patch_token = "" if patch is None else f"{patch}-" + return f"{file_pattern}_conv-{patch_token}{idx}.fits" + + def transform_shape(mom_list, jac): """Transform Shape. @@ -378,6 +457,7 @@ def params_default(self): self._params = { "input_base_dir": ".", "output_base_dir": ".", + "version_cat": "v2.0", "mode": "merge", "patches": "", "psf": "psfex", @@ -386,6 +466,7 @@ def params_default(self): self._short_options = { "input_base_dir": "-i", + "version_cat": "-V", "mode": "-m", "psf": "-p", "patches": "-P", @@ -395,15 +476,19 @@ def params_default(self): self._help_strings = { "input_base_dir": ( - "input base dir, runs are expected in" - + " /P/output;" - " default is {}" + "input base dir; for v1.x runs are expected in" + + " /P/output, for v2.0 (patch-less) in" + + " /output; default is {}" + ), + "version_cat": ( + "catalogue major version, allowed are v1.3, v1.4, v1.5, v1.6," + + " v2.0; v2.0 is patch-less; default is {}" ), "mode": ( "run mode, allowed are 'merge', 'test'; default is" + " '{}'" ), "psf": "PSF model, allowed are 'psfex' and 'mccd'; default is {}", - "patches": "(list of) input patches", + "patches": "(list of) input patches; ignored for v2.0", } # Output column names with types @@ -454,7 +539,22 @@ def run(self): Main processing function. """ - if self._params["mode"] == "test": + # Guard against a mistyped version silently falling through to the + # v1.x patch loop (e.g. ``-V v2`` or ``-V 2.0``). + allowed_versions = ("v1.3", "v1.4", "v1.5", "v1.6", "v2.0") + if self._params["version_cat"] not in allowed_versions: + raise ValueError( + f"Invalid version {self._params['version_cat']}; allowed are" + + f" {', '.join(allowed_versions)}" + ) + + # v2.0 removes the patch concept: a single patch-less run root. For + # v1.x, iterate over the requested sky patches as before. ``patch`` is + # None in the patch-less case, which drops the patch token from the + # input path and the output filename. + if self._params["version_cat"] == "v2.0": + patch_nums = [None] + elif self._params["mode"] == "test": patch_nums = ["3", "4"] else: patch_nums = cs_args.my_string_split(self._params["patches"]) @@ -464,13 +564,16 @@ def run(self): # Loop over patches for patch in patch_nums: - output_dir = f"{self._params['output_base_dir']}/P{patch}" - if not os.path.isdir(output_dir): - os.makedirs(output_dir, exist_ok=False) + patch_dir, output_dir = collate_paths( + self._params["input_base_dir"], + self._params["output_base_dir"], + patch, + ) + print("Running patch-less (v2.0)" if patch is None else f"Running patch: {patch}") - print("Running patch:", patch) + if not os.path.isdir(output_dir): + os.makedirs(output_dir, exist_ok=True) - patch_dir = f"{self._params['input_base_dir']}/P{patch}/output/" subdirs = f"{patch_dir}/{self._params['sub_dir_pattern']}*" exp_run_dirs = glob.glob(subdirs) n_exp_runs = len(exp_run_dirs) @@ -516,8 +619,10 @@ def transform_exposures(self, output_dir, patch, idx, exp_run_dir): """ output_path = ( - f"{output_dir}/{self._params['file_pattern_psfint']}_conv" - + f"-{patch}-{idx}.fits" + f"{output_dir}/" + + output_filename( + self._params["file_pattern_psfint"], patch, idx + ) ) if os.path.exists(output_path): print(f"Skipping transform_exposures, file {output_path} exists") @@ -555,7 +660,6 @@ def transform_exposures(self, output_dir, patch, idx, exp_run_dir): if self._params["psf"] == "psfex": psf_file_hdus = fits.open(psf_file_path, memmap=False) psf_file = psf_file_hdus[2].data - header_file = psf_file_hdus[1].data psf_file_hdus.close() mod = "RA" else: @@ -564,60 +668,10 @@ def transform_exposures(self, output_dir, patch, idx, exp_run_dir): except Exception: continue - new_e1_psf = np.zeros_like(psf_file[mod]) - new_e2_psf = np.zeros_like(psf_file[mod]) - new_sig_psf = np.zeros_like(psf_file[mod]) - new_e1_star = np.zeros_like(psf_file[mod]) - new_e2_star = np.zeros_like(psf_file[mod]) - new_sig_star = np.zeros_like(psf_file[mod]) - if self._params["psf"] == "psfex": - header = fits.Header.fromstring( - "\n".join(header_file[0][0]), sep="\n" - ) - wcs = galsim.AstropyWCS(header=header) - - k = 0 - for ind, obj in enumerate(psf_file): - try: - # jac = wcs.jacobian(world_pos=galsim.CelestialCoord( - # ra=obj["RA"]*galsim.degrees, - # dec=obj["DEC"]*galsim.degrees - # )) - jac = wcs.jacobian( - image_pos=galsim.PositionD( - obj["X"], - obj["Y"], - ) - ) - except Exception: - continue - g1_psf_tmp, g2_psf_tmp, sig_psf_tmp = transform_shape( - [ - obj["E1_PSF_HSM"], - obj["E2_PSF_HSM"], - obj["SIGMA_PSF_HSM"], - ], - jac, - ) - - new_e1_psf[ind] = g1_psf_tmp - new_e2_psf[ind] = g2_psf_tmp - new_sig_psf[ind] = sig_psf_tmp - - g1_star_tmp, g2_star_tmp, sig_star_tmp = transform_shape( - [ - obj["E1_STAR_HSM"], - obj["E2_STAR_HSM"], - obj["SIGMA_STAR_HSM"], - ], - jac, - ) - new_e1_star[ind] = g1_star_tmp - new_e2_star[ind] = g2_star_tmp - new_sig_star[ind] = sig_star_tmp - k += 1 - + # HSM ellipticities and sizes are measured directly in world + # coordinates upstream (FindAdaptiveMom use_sky_coords=True), so + # they are passed straight through; only positions are collated. exp_cat = np.array( list( map( @@ -628,13 +682,13 @@ def transform_exposures(self, output_dir, patch, idx, exp_run_dir): psf_file["Y"], psf_file["RA"], psf_file["DEC"], - new_e1_psf, - new_e2_psf, - new_sig_psf, + psf_file["E1_PSF_HSM"], + psf_file["E2_PSF_HSM"], + psf_file["SIGMA_PSF_HSM"], psf_file["FLAG_PSF_HSM"], - new_e1_star, - new_e2_star, - new_sig_star, + psf_file["E1_STAR_HSM"], + psf_file["E2_STAR_HSM"], + psf_file["SIGMA_STAR_HSM"], psf_file["FLAG_STAR_HSM"], np.ones_like(psf_file["RA"], dtype=int) * ccd_id, @@ -661,8 +715,23 @@ def transform_exposures(self, output_dir, patch, idx, exp_run_dir): ] ) + # Local-to-CCD position: subtract each CCD's focal-plane shift. new_x = np.zeros_like(psf_file[mod]) new_y = np.zeros_like(psf_file[mod]) + + # The MCCD PSF_MOM_LIST/STAR_MOM_LIST columns come from the + # external mccd fit-validation code, which still measures HSM + # moments in the pixel frame; rotate them into world coordinates + # via the per-CCD WCS Jacobian. This rotation stays until mccd + # itself adopts use_sky_coords (see the module docstring). The + # in-repo PSFEx / MCCD-interpolation paths are already in world + # coordinates and are passed through unrotated. + new_e1_psf = np.zeros_like(psf_file[mod]) + new_e2_psf = np.zeros_like(psf_file[mod]) + new_sig_psf = np.zeros_like(psf_file[mod]) + new_e1_star = np.zeros_like(psf_file[mod]) + new_e2_star = np.zeros_like(psf_file[mod]) + new_sig_star = np.zeros_like(psf_file[mod]) new_flag_psf = np.zeros_like(psf_file[mod]) new_flag_star = np.zeros_like(psf_file[mod]) for ccd_id in range(40): diff --git a/src/shapepipe/modules/mccd_package/mccd_interpolation_script.py b/src/shapepipe/modules/mccd_package/mccd_interpolation_script.py index bfd9f341..a01ca095 100644 --- a/src/shapepipe/modules/mccd_package/mccd_interpolation_script.py +++ b/src/shapepipe/modules/mccd_package/mccd_interpolation_script.py @@ -14,6 +14,7 @@ from sqlitedict import SqliteDict from shapepipe.pipeline import file_io +from shapepipe.modules.psfex_interp_package.psfex_interp import local_wcs_list try: import galsim.hsm as hsm @@ -204,18 +205,34 @@ def _get_galaxy_positions(self): raise KeyError(pos_param_err) galcat.close() - def _get_psfshapes(self): + def _get_psfshapes(self, wcs_list=None): """Get PSF Shapes. Compute shapes of PSF at galaxy positions using HSM. + Parameters + ---------- + wcs_list : list of galsim.BaseWCS, optional + Local WCS at each object position. When provided, HSM adaptive + moments are measured in world coordinates (``use_sky_coords=True``), + so the returned shapes and sizes are already in sky coordinates. + """ if import_fail: raise ImportError("Galsim is required to get shapes information") - psf_moms = [ - hsm.FindAdaptiveMom(Image(psf), strict=False) for psf in self.interp_PSFs - ] + if wcs_list is not None: + psf_moms = [ + hsm.FindAdaptiveMom( + Image(psf, wcs=wcs), strict=False, use_sky_coords=True + ) + for psf, wcs in zip(self.interp_PSFs, wcs_list) + ] + else: + psf_moms = [ + hsm.FindAdaptiveMom(Image(psf), strict=False) + for psf in self.interp_PSFs + ] self.psf_shapes = np.array( [ @@ -253,7 +270,7 @@ def _write_output(self): data = {"VIGNET": self.interp_PSFs} output.save_as_fits(data, sex_cat_path=self._galcat_path) - def _get_starshapes(self, star_vign): + def _get_starshapes(self, star_vign, wcs_list=None): """Get Star Shapes. Compute shapes of stars at stars positions using HSM. @@ -262,6 +279,10 @@ def _get_starshapes(self, star_vign): ---------- star_vign : numpy.ndarray Array containing the star's vignets + wcs_list : list of galsim.BaseWCS, optional + Local WCS at each star position; see :meth:`_get_psfshapes`. When + provided, moments are measured in world coordinates + (``use_sky_coords=True``). """ if import_fail: @@ -270,10 +291,23 @@ def _get_starshapes(self, star_vign): masks = np.zeros_like(star_vign) masks[np.where(star_vign == -1e30)] = 1 - star_moms = [ - hsm.FindAdaptiveMom(Image(star), badpix=Image(mask), strict=False) - for star, mask in zip(star_vign, masks) - ] + if wcs_list is not None: + star_moms = [ + hsm.FindAdaptiveMom( + Image(star, wcs=wcs), + badpix=Image(mask), + strict=False, + use_sky_coords=True, + ) + for star, mask, wcs in zip(star_vign, masks, wcs_list) + ] + else: + star_moms = [ + hsm.FindAdaptiveMom( + Image(star), badpix=Image(mask), strict=False + ) + for star, mask in zip(star_vign, masks) + ] self.star_shapes = np.array( [ @@ -451,6 +485,15 @@ def _interpolate_me(self): 0, ) ).T + # 1-indexed (FITS) positions for local_wcs_list, which wraps + # galsim.AstropyWCS; gal_pos (origin 0) stays for MCCD interp. + gal_pos_wcs = np.array( + self._f_wcs_file[exp_name][ccd]["WCS"].all_world2pix( + self.gal_pos[:, 0][ind_obj], + self.gal_pos[:, 1][ind_obj], + 1, + ) + ).T self.interp_PSFs = interp_MCCD(mccd_model_path, gal_pos, ccd) @@ -497,7 +540,11 @@ def _interpolate_me(self): array_id = np.concatenate((array_id, np.copy(obj_id))) if self._compute_shape: - self._get_psfshapes() + self._get_psfshapes( + local_wcs_list( + self._f_wcs_file[exp_name][ccd]["WCS"], gal_pos_wcs + ) + ) if array_shape is None: array_shape = np.copy(self.psf_shapes) else: diff --git a/src/shapepipe/modules/mccd_package/shapepipe_auxiliary_mccd.py b/src/shapepipe/modules/mccd_package/shapepipe_auxiliary_mccd.py index 943433e8..f0aea4f6 100644 --- a/src/shapepipe/modules/mccd_package/shapepipe_auxiliary_mccd.py +++ b/src/shapepipe/modules/mccd_package/shapepipe_auxiliary_mccd.py @@ -16,6 +16,7 @@ from astropy.io import fits from cs_util import size as cs_size +from shapepipe.modules.psfex_interp_package.psfex_interp import local_wcs_list from shapepipe.pipeline import file_io NOT_ENOUGH_STARS = "Not enough stars to train the model." @@ -403,6 +404,10 @@ def mccd_interpolation_pipeline( x_pos = galcat[2].data[pos_params[0]] y_pos = galcat[2].data[pos_params[1]] interp_pos = np.array([x_pos, y_pos]).T + # CCD image WCS, stored as a card list in the FITS_LDAC LDAC_IMHEAD HDU + wcs_header = fits.Header.fromstring( + "\n".join(galcat[1].data[0][0]), sep="\n" + ) # Close catalog galcat.close() @@ -411,9 +416,18 @@ def mccd_interpolation_pipeline( if interp_PSFs is not None: if get_shapes: + # Local WCS at each object position lets HSM measure adaptive + # moments directly in world coordinates (use_sky_coords=True). + wcs_list = local_wcs_list( + galsim.AstropyWCS(header=wcs_header), interp_pos + ) PSF_moms = [ - galsim.hsm.FindAdaptiveMom(galsim.Image(psf), strict=False) - for psf in interp_PSFs + galsim.hsm.FindAdaptiveMom( + galsim.Image(psf, wcs=wcs), + strict=False, + use_sky_coords=True, + ) + for psf, wcs in zip(interp_PSFs, wcs_list) ] PSF_shapes = np.array( diff --git a/src/shapepipe/modules/merge_starcat_package/merge_starcat.py b/src/shapepipe/modules/merge_starcat_package/merge_starcat.py index 7bc0cb76..58fdd376 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 c14aeb13..044ed009 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 @@ -17,6 +18,7 @@ from shapepipe.pipeline import file_io try: + import galsim import galsim.hsm as hsm from galsim import Image except ImportError: @@ -30,6 +32,148 @@ FILE_NOT_FOUND = "File_not_found" +def local_wcs_list(wcs, positions): + """Local WCS List. + + Build a list of local (Jacobian) WCS from an astropy/galsim WCS, one per + object position. Attaching the local WCS to an object-centred stamp lets + ``FindAdaptiveMom(use_sky_coords=True)`` measure moments directly in world + coordinates: at the stamp's ``true_center`` galsim evaluates this local WCS, + which must therefore be the linearisation of the full WCS *at the object's + position in the original image*, not at the stamp centre. + + Parameters + ---------- + wcs : galsim.BaseWCS or astropy.wcs.WCS + Full (celestial) WCS for the CCD. + positions : numpy.ndarray + Object positions ``(N, 2)`` in image (pixel) coordinates. + + Returns + ------- + list of galsim.JacobianWCS + One local WCS per object position. + + """ + if not isinstance(wcs, galsim.BaseWCS): + wcs = galsim.AstropyWCS(wcs=wcs) + return [ + wcs.local(image_pos=galsim.PositionD(float(x), float(y))) + for x, y in 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. @@ -145,7 +289,9 @@ def process(self): self._w_log.info(f"Psf model file {self._dotpsf_path} not found.") else: if self._compute_shape: - self._get_psfshapes() + self._get_psfshapes( + local_wcs_list(self._galcat_wcs(), self.gal_pos) + ) self._write_output() @@ -307,18 +453,45 @@ def _interpolate(self): self.gal_pos, ) - def _get_psfshapes(self): + def _get_psfshapes(self, wcs_list=None): """Get PSF Shapes. Compute shapes of PSF at galaxy positions using HSM. + Parameters + ---------- + wcs_list : list of galsim.BaseWCS, optional + Local WCS at each object position (see :func:`local_wcs_list`). When + provided, HSM adaptive moments are measured directly in world + coordinates (``use_sky_coords=True``), so the ``g1``/``g2``/``sigma`` + 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") - psf_moms = [ - hsm.FindAdaptiveMom(Image(psf), strict=False) for psf in self.interp_PSFs - ] + if wcs_list is not None: + psf_moms = [ + hsm.FindAdaptiveMom( + Image(psf, wcs=wcs), strict=False, use_sky_coords=True + ) + for psf, wcs in zip(self.interp_PSFs, wcs_list) + ] + else: + psf_moms = [ + hsm.FindAdaptiveMom(Image(psf), strict=False) + for psf in self.interp_PSFs + ] + wcs_list = [None] * len(self.interp_PSFs) self.psf_shapes = np.array( [ @@ -327,8 +500,9 @@ def _get_psfshapes(self): 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) ] ) @@ -351,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} @@ -407,21 +584,37 @@ def process_validation(self, psfex_cat_path): star_dict["SNR"] = np.copy(star_cat.get_data()["SNR_WIN"]) star_cat.close() - self._get_psfshapes() - self._get_starshapes(star_vign) + wcs_list = local_wcs_list( + self._galcat_wcs(), + np.column_stack((star_dict["X"], star_dict["Y"])), + ) + + self._get_psfshapes(wcs_list) + self._get_starshapes(star_vign, wcs_list) psfex_cat_dict = self._get_psfexcatdict(psfex_cat_path) self._write_output_validation(star_dict, psfex_cat_dict) - def _get_starshapes(self, star_vign): + def _get_starshapes(self, star_vign, wcs_list=None): """Get Star Shapes. Compute shapes of stars at stars positions using HSM. Parameters ---------- - numpy.ndarray + star_vign : numpy.ndarray Array containing the star's vignets. + wcs_list : list of galsim.BaseWCS, optional + Local WCS at each star position; see :meth:`_get_psfshapes`. When + 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: @@ -430,10 +623,24 @@ def _get_starshapes(self, star_vign): masks = np.zeros_like(star_vign) masks[np.where(star_vign == -1e30)] = 1 - star_moms = [ - hsm.FindAdaptiveMom(Image(star), badpix=Image(mask), strict=False) - for star, mask in zip(star_vign, masks) - ] + if wcs_list is not None: + star_moms = [ + hsm.FindAdaptiveMom( + Image(star, wcs=wcs), + badpix=Image(mask), + strict=False, + use_sky_coords=True, + ) + for star, mask, wcs in zip(star_vign, masks, wcs_list) + ] + else: + star_moms = [ + hsm.FindAdaptiveMom( + Image(star), badpix=Image(mask), strict=False + ) + for star, mask in zip(star_vign, masks) + ] + wcs_list = [None] * len(star_vign) self.star_shapes = np.array( [ @@ -442,11 +649,34 @@ def _get_starshapes(self, star_vign): 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 + ) ] ) + def _galcat_wcs(self): + """Galaxy Catalogue WCS. + + Read the CCD image WCS from the SExtractor FITS_LDAC galaxy catalogue. + The original image header is stored as a card list in the ``LDAC_IMHEAD`` + HDU (HDU 1); this is the same header the (now retired) pixel-to-world + conversion reconstructed to rotate shapes. + + Returns + ------- + galsim.AstropyWCS + WCS of the CCD the stars/PSFs live on. + + """ + with fits.open(self._galcat_path, memmap=False) as hdu_list: + header = fits.Header.fromstring( + "\n".join(hdu_list[1].data[0][0]), sep="\n" + ) + return galsim.AstropyWCS(header=header) + def _get_psfexcatdict(self, psfex_cat_path): """Get PSFEx Catalogue Dictionary. @@ -504,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} @@ -656,6 +892,15 @@ def _interpolate_me(self): 0, ) ).T + # 1-indexed (FITS) positions for local_wcs_list, which wraps + # galsim.AstropyWCS; gal_pos (origin 0) stays for PSFEx interp. + gal_pos_wcs = np.array( + self._f_wcs_file[exp_name][ccd]["WCS"].all_world2pix( + self.gal_pos[:, 0][ind_obj], + self.gal_pos[:, 1][ind_obj], + 1, + ) + ).T self.interp_PSFs = self.interpsfex( dot_psf_path, @@ -700,7 +945,11 @@ def _interpolate_me(self): array_id = np.concatenate((array_id, np.copy(obj_id))) if self._compute_shape: - self._get_psfshapes() + self._get_psfshapes( + local_wcs_list( + self._f_wcs_file[exp_name][ccd]["WCS"], gal_pos_wcs + ) + ) if array_shape is None: array_shape = np.copy(self.psf_shapes) else: @@ -767,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_collate_star_cat.py b/tests/module/test_collate_star_cat.py new file mode 100644 index 00000000..20d52ce2 --- /dev/null +++ b/tests/module/test_collate_star_cat.py @@ -0,0 +1,70 @@ +"""UNIT TESTS FOR STAR-CATALOGUE COLLATION PATHS. + +Pin the patch vs patch-less (v2.0) path and filename convention of +``scripts/python/collate_star_cat.py``. Runs up to v1.6 carry a ``P`` +token in both the input run directory and the output filename; v2.0 is +patch-less (``patch is None``) and drops that token, reading from a single +``/output`` root and writing ``validation_psf_conv-.fits`` — the +name still matched by the downstream ``validation_psf_conv-*`` glob. +""" + +import importlib.util +from pathlib import Path + +import pytest + +# The collation script lives under scripts/python (not an importable package), +# so load it by path. +_SCRIPT = ( + Path(__file__).resolve().parents[2] + / "scripts" + / "python" + / "collate_star_cat.py" +) +_spec = importlib.util.spec_from_file_location("collate_star_cat", _SCRIPT) +collate_star_cat = importlib.util.module_from_spec(_spec) +_spec.loader.exec_module(collate_star_cat) + + +@pytest.mark.parametrize( + "patch, exp_input, exp_output", + [ + ("3", "in/P3/output/", "out/P3"), + (None, "in/output/", "out"), + ], +) +def test_collate_paths(patch, exp_input, exp_output): + """v1.x carries the P token; v2.0 (patch None) drops it.""" + assert collate_star_cat.collate_paths("in", "out", patch) == ( + exp_input, + exp_output, + ) + + +@pytest.mark.parametrize( + "patch, expected", + [ + ("3", "validation_psf_conv-3-0.fits"), + (None, "validation_psf_conv-0.fits"), + ], +) +def test_output_filename(patch, expected): + """The patch token is present for v1.x and absent for v2.0.""" + assert collate_star_cat.output_filename("validation_psf", patch, 0) == expected + + +def test_output_filename_matches_downstream_glob(): + """Both layouts stay under the downstream ``validation_psf_conv-*`` glob.""" + for patch in ("1", None): + assert collate_star_cat.output_filename( + "validation_psf", patch, 5 + ).startswith("validation_psf_conv-") + + +@pytest.mark.parametrize("bad", ["v2", "2.0", "v1.7", ""]) +def test_invalid_version_raises(bad): + """A mistyped -V is rejected rather than falling through to v1.x.""" + obj = collate_star_cat.Convert() + obj._params["version_cat"] = bad + with pytest.raises(ValueError): + obj.run() diff --git a/tests/module/test_hsm_sky_coords.py b/tests/module/test_hsm_sky_coords.py new file mode 100644 index 00000000..3d19c0af --- /dev/null +++ b/tests/module/test_hsm_sky_coords.py @@ -0,0 +1,128 @@ +"""UNIT TESTS FOR HSM SKY-COORDINATE MOMENTS. + +Pin the equivalence that motivates measuring PSF/star HSM adaptive moments +directly in world coordinates (``FindAdaptiveMom(use_sky_coords=True)``) instead +of measuring in the pixel frame and rotating afterwards with a WCS Jacobian. + +The reference ``_transform_shape`` below is a verbatim copy of the shape-rotation +math that ``scripts/python/collate_star_cat.py`` used to apply (``transform_shape``, +removed in this change). The test draws an elliptical Gaussian on a stamp with a +nontrivial local WCS (rotation + shear + scale, matching a real CD matrix), then +checks that: + + pixel-frame moments → _transform_shape(local WCS Jacobian) + +agrees with the sky-coordinate moments to numerical precision. + +The local WCS here has positive determinant (no parity flip): production CFHT +runs decomposed with ``flip=False`` (the old ``transform_shape`` flip branch, +which is not a correct galsim reflection, never fired), so the faithful +equivalence to pin is the no-flip one. +""" + +import galsim +import numpy as np +import numpy.testing as npt +import pytest + + +def _transform_shape(mom_list, jac): + """Reference: the old pixel-frame → world shape rotation (verbatim). + + Copied from the removed ``collate_star_cat.transform_shape``. ``mom_list`` + is ``[g1, g2, sigma]`` measured in the pixel frame; ``jac`` is the local WCS. + """ + scale, shear, theta, flip = jac.getDecomposition() + + sig_tmp = mom_list[2] * scale + shape = galsim.Shear(g1=mom_list[0], g2=mom_list[1]) + if flip: + shape = galsim.Shear(g1=-shape.g1, g2=shape.g2) + shape = galsim.Shear(g=shape.g, beta=shape.beta + theta) + shape = shear + shape + + return shape.g1, shape.g2, sig_tmp + + +# A CD-matrix-like local WCS: 0.187 arcsec/pix scale, ~28 deg rotation and a +# small shear, positive determinant (no parity flip). +_SCALE = 0.187 +_THETA = 28.0 * np.pi / 180.0 +_c, _s = np.cos(_THETA), np.sin(_THETA) +_SHEAR = galsim.Shear(g1=0.06, g2=-0.04) +_A = _SHEAR.getMatrix() * _SCALE @ np.array([[_c, -_s], [_s, _c]]) +LOCAL_WCS = galsim.JacobianWCS(_A[0, 0], _A[0, 1], _A[1, 0], _A[1, 1]) + + +@pytest.mark.parametrize( + "e1, e2, sigma_pix", + [ + (0.3, 0.1, 3.0), + (-0.2, 0.25, 2.5), + (0.15, -0.35, 4.0), + (0.0, 0.0, 3.5), + ], +) +def test_sky_coords_moments_match_jacobian_rotation(e1, e2, sigma_pix): + """use_sky_coords=True reproduces pixel-frame moments + Jacobian rotation. + + The local WCS has positive determinant, so ``getDecomposition`` returns + ``flip=False`` and the reference path is a pure scale+rotation+shear — the + transformation ``use_sky_coords`` performs internally. + """ + # Sanity: the reference decomposition sees no parity flip. + assert LOCAL_WCS.getDecomposition()[3] is False + + # An elliptical Gaussian PSF, drawn on a unit-pixel-scale stamp; sigma is in + # pixels because the pixel scale is 1. + stamp = ( + galsim.Gaussian(sigma=sigma_pix) + .shear(e1=e1, e2=e2) + .drawImage(nx=51, ny=51, scale=1.0, method="no_pixel") + ) + + # Pixel-frame measurement, then the old Jacobian rotation. + moms_pix = galsim.hsm.FindAdaptiveMom( + galsim.Image(stamp.array), strict=False + ) + g1_ref, g2_ref, sigma_ref = _transform_shape( + [ + moms_pix.observed_shape.g1, + moms_pix.observed_shape.g2, + moms_pix.moments_sigma, + ], + LOCAL_WCS, + ) + + # Sky-coordinate measurement: same stamp, local WCS attached. + moms_sky = galsim.hsm.FindAdaptiveMom( + galsim.Image(stamp.array, wcs=LOCAL_WCS), + strict=False, + use_sky_coords=True, + ) + + npt.assert_allclose(moms_sky.observed_shape.g1, g1_ref, atol=1e-9, rtol=1e-6) + npt.assert_allclose(moms_sky.observed_shape.g2, g2_ref, atol=1e-9, rtol=1e-6) + npt.assert_allclose(moms_sky.moments_sigma, sigma_ref, rtol=1e-6) + + +def test_sky_coords_sigma_is_scaled_to_world_units(): + """The sky-frame size is the pixel size times the WCS scale (arcsec/pix). + + Pins the deliberate unit change: SIGMA_*_HSM leaves the measurement in world + (arcsec) units rather than pixels once use_sky_coords is on. + """ + stamp = ( + galsim.Gaussian(sigma=3.0) + .drawImage(nx=51, ny=51, scale=1.0, method="no_pixel") + ) + sigma_pix = galsim.hsm.FindAdaptiveMom( + galsim.Image(stamp.array), strict=False + ).moments_sigma + sigma_sky = galsim.hsm.FindAdaptiveMom( + galsim.Image(stamp.array, wcs=LOCAL_WCS), + strict=False, + use_sky_coords=True, + ).moments_sigma + + npt.assert_allclose(sigma_sky, sigma_pix * _SCALE, rtol=1e-6) diff --git a/tests/module/test_psf_fourth_moments.py b/tests/module/test_psf_fourth_moments.py new file mode 100644 index 00000000..4782ca5f --- /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, + )