PSF fourth-order moments in the sky frame (spin-2 + rho4) - #859
Open
cailmdaley wants to merge 10 commits into
Open
PSF fourth-order moments in the sky frame (spin-2 + rho4)#859cailmdaley wants to merge 10 commits into
cailmdaley wants to merge 10 commits into
Conversation
Attach the local (Jacobian) WCS at each object position to the PSF/star stamps and pass use_sky_coords=True to galsim FindAdaptiveMom, so PSFEx and MCCD HSM ellipticities and sizes come out already in world coordinates. This makes the downstream WCS-Jacobian shape rotation (convert_psf_pix2world) redundant. Threaded at every in-repo HSM PSF/star call site: - psfex_interp: multi-epoch (per-CCD WCS log), single-exposure validation and classic (WCS from the galcat FITS_LDAC image header via _galcat_wcs) - mccd_interpolation_script: multi-epoch (per-CCD WCS log) - shapepipe_auxiliary_mccd: interpolation (WCS from galcat LDAC header) local_wcs_list() builds the per-object local WCS: use_sky_coords evaluates the stamp WCS at true_center, so it must be the full WCS linearised at the object's image position, not at the stamp centre. SIGMA_*_HSM (hence HSM_T_* = 2 sigma^2) is now in sky units (arcsec) rather than pixels; this is the intended, physically correct size. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01Gao7Pi5QDBgnTuagV7HT2w
HSM ellipticities and sizes are measured directly in world coordinates upstream, so the WCS-Jacobian shape rotation (transform_shape, getDecomposition) is redundant. Remove it from both the PSFEx and MCCD branches; the script keeps its remaining job of collating positions (X/Y/ RA/DEC), assigning the MCCD focal-plane CCD id, and merging per-exposure catalogues into per-patch star catalogues. Drops the galsim dependency. Caveat documented in the module docstring: the MCCD PSF_MOM_LIST/STAR_MOM_LIST columns come from the external mccd fit-validation code, still measured in the pixel frame; they are now passed through unrotated. Update the CANFAR doc to describe the step as collation, not pix->world shape conversion. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01Gao7Pi5QDBgnTuagV7HT2w
Equivalence test for use_sky_coords=True: draw an elliptical Gaussian on a stamp with a nontrivial local WCS (rotation + shear + 0.187 arcsec/pix scale, positive determinant), and assert that measuring in the pixel frame then applying the removed convert_psf_pix2world.transform_shape rotation (copied in verbatim as reference) agrees with measuring directly in sky coordinates, to numerical precision, for e1/e2 and sigma. A second test pins the deliberate unit change: the sky-frame sigma is the pixel sigma times the WCS scale. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01Gao7Pi5QDBgnTuagV7HT2w
…ranch The MCCD PSF_MOM_LIST/STAR_MOM_LIST columns come from the external mccd fit-validation code (mccd.auxiliary_fun.mccd_validation), which still measures HSM adaptive moments in the pixel frame. Stripping the WCS-Jacobian rotation there (71d502f) would have handed MergeStarCatMCCD pixel-frame ellipticities. Restore transform_shape + the per-CCD getDecomposition rotation for the MCCD MOM_LIST branch only (and the galsim import it needs); the in-repo PSFEx and MCCD-interpolation paths stay unrotated, as they now measure in world coordinates upstream. Update the module docstring caveat: this branch keeps the rotation until mccd itself adopts use_sky_coords. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01Gao7Pi5QDBgnTuagV7HT2w
local_wcs_list wraps galsim.AstropyWCS, which follows the 1-indexed FITS pixel convention, but the multi-epoch paths fed it positions from all_world2pix(..., 0) (origin 0) — evaluating the local Jacobian one pixel off. Compute a separate origin=1 position array (gal_pos_wcs) for local_wcs_list; the origin-0 gal_pos stays for the PSF-model interpolation (interpsfex / interp_MCCD), which is its only other use. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01Gao7Pi5QDBgnTuagV7HT2w
…add patch-less v2.0 Re-scope the script to its actual remaining purpose: collating per-exposure PSF validation catalogues into star catalogues. The WCS-Jacobian shape rotation is gone (HSM moments are measured in sky coordinates upstream), so "convert" no longer describes it. The external-mccd PSF_MOM_LIST/STAR_MOM_LIST rotation branch stays until mccd adopts use_sky_coords. Add a -V/--version option mirroring get_ccds_with_psf. v1.x is unchanged (loop over -P patches, <base>/P<patch>/output, validation_psf_conv-<patch>-<idx>). v2.0 is patch-less: single <base>/output root, patch token dropped from the output filename (validation_psf_conv-<idx>). The downstream combine_runs.bash glob (validation_psf_conv-*) matches both. Path/filename logic extracted into pure collate_paths/output_filename helpers with unit tests; docs collate section updated with both invocations. The script command name follows the file (Dockerfile symlinks by basename). Co-Authored-By: Claude Fable 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01WqNFhq9uGpDpprYeuGeb8k
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 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01K98geWHTuT62whqjKocfCm
…ments 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 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01K98geWHTuT62whqjKocfCm
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 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01K98geWHTuT62whqjKocfCm
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 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01K98geWHTuT62whqjKocfCm
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Closes #697. Supersedes #698. Follow-up to #812; based on its branch so the diff shows only the fourth-moment work — retarget to
developafter #812 merges.What this does
Adds fourth-order moment columns to the PSF/star measurement chain:
HSM_M4_1/2_{PSF,STAR}(spin-2 combinationsM40 − M04,2(M13 + M31), PSFHOME-style) andHSM_RHO4_{PSF,STAR}(galsim's spin-0 kurtosis, previously computed and discarded).How
All computation happens in the world frame. The pixel grid maps to sky offsets through the local WCS Jacobian. The whitening matrix comes from the sky-frame second moments returned by
use_sky_coords=True. This makes the spin-2 columns frame-invariant; the pixel-frame computation in #698 was not (verified numerically — it rotates with the CCD orientation).Columns flow through all write paths: single-epoch, validation, multi-epoch, and
merge_starcatinto the merged star catalogue.Two latent #698 issues do not carry over: the residual imaginary part from
scipy.linalg.sqrtmis dropped, and themerge_starcatcopy-paste bug (M_4_PSF_2reading the_1column) is avoided.Tests
23 pass (13 new + masked-pixel and centroid coverage from review, plus the #812 suite — no regression). Analytic guards: round Gaussian → spin-2 ≈ 0, rho4 ≈ 2; frame invariance across four WCS orientations at rtol=1e-5 (actual agreement ~1e-8);
-1e30mask sentinels; off-center sources. A physics note encoded in the tests: whitening circularises any sheared circular profile, so a single elliptical Moffat also gives spin-2 ≈ 0 — the spin-2 signal traces radial ellipticity gradients. The non-trivial fixture is a two-component coaxial Gaussian.Known follow-up
sp_validation's residual-star plot reads the oldM_4_STAR_1-style names; it needs 4 one-line renames to theHSM_M4_*grammar (separate repo, separate PR).🤖 Generated with Claude Code
https://claude.ai/code/session_01K98geWHTuT62whqjKocfCm