Skip to content

Improve accuracy for high-resolution images in reproject_exact#589

Merged
astrofrog merged 8 commits intoastropy:mainfrom
astrofrog:fix-precision-reproject-exact
Mar 24, 2026
Merged

Improve accuracy for high-resolution images in reproject_exact#589
astrofrog merged 8 commits intoastropy:mainfrom
astrofrog:fix-precision-reproject-exact

Conversation

@astrofrog
Copy link
Copy Markdown
Member

As described in #199, the reproject_exact function suffers from precision issues when reprojecting images with sub-arcsecond pixel scales, which is particularly problematic for instruments like JWST NIRCam (0.03" pixels). There are two reasons for this:

  • Girard's theorem computes spherical polygon areas as sum_of_angles - (N-2)*π, which loses precision catastrophically when the polygon is tiny and the sum is nearly equal to (N-2)*π.
  • The spherical polygon intersection algorithm in the relies on a TOLERANCE constant (~9e-4 arcsec) that causes incorrect intersection results when pixel features are smaller than this threshold.

This PR fixes both issues by adding a planar approximation path in computeOverlap. When all vertices of both the input and output pixels are within ~3.4 arcmin of their centroid, the code projects them onto a local tangent plane and performs the entire overlap computation in 2D — using Sutherland-Hodgman polygon clipping and the shoelace formula instead of the spherical intersection and Girard's theorem. This avoids all sources of precision loss, since coordinate values in the tangent plane are proportional to the polygon size rather than being tiny perturbations on a unit sphere. The same planar approximation is also applied as a fallback within the Girard function itself, to handle the edge case of tiny sliver overlaps between large pixels. With these changes, flux conservation is accurate down to ~1e-5 arcsec, compared to the previous limit of ~0.07 arcsec.

Here's an updated plot of accuracy vs resolution based on this PR:

ratios

For comparison, the original plot looked like:

Use planar approximation (shoelace formula) instead of Girard's theorem
for small spherical polygons, avoiding catastrophic cancellation in the
area = sum_of_angles - (N-2)*pi computation. This fixes flux conservation
for pixel scales down to ~0.001 arcsec (previously failed below ~0.07
arcsec), which is important for instruments like JWST NIRCam.
…xels

For small pixels, bypass both the spherical polygon intersection
(ComputeIntersection) and Girard's theorem by doing the full overlap
computation in a local tangent plane using Sutherland-Hodgman 2D
polygon clipping and the shoelace formula. This eliminates all
precision issues from the TOLERANCE constant and spherical geometry
for sub-arcsecond pixels, extending accurate flux conservation down
to ~1e-6 arcsec (previously failed below ~0.07 arcsec).
@astrofrog
Copy link
Copy Markdown
Member Author

@acikota - your example from Dropbox now returns:

pixel area ratio between template and input image:  226.47954059928594

reprojecting...
sum initial:  100.0
sum reprojected:  0.4415410594776611
ratio initial/reprojected sum: 226.47950366903376

which seems ok?

@codecov
Copy link
Copy Markdown

codecov bot commented Mar 21, 2026

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 91.06%. Comparing base (bc02fca) to head (64ffcea).
⚠️ Report is 11 commits behind head on main.

Additional details and impacted files
@@           Coverage Diff           @@
##             main     #589   +/-   ##
=======================================
  Coverage   91.06%   91.06%           
=======================================
  Files          30       30           
  Lines        1858     1858           
=======================================
  Hits         1692     1692           
  Misses        166      166           

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:
  • ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

@astrofrog
Copy link
Copy Markdown
Member Author

So the interesting thing is that the planar approximation is actually better all the way up to 20" pixels. The comparison plot above looks like this if I zoom in to the y axis before the changes:

ratios

and like this after the changes:

ratios_new

So this means that actually for typical images the plane approximation is actually better as the numerical precision improvements are more important than the algorithmic differences.

@astrofrog
Copy link
Copy Markdown
Member Author

I've corrected the calculation of the comparison at large resolutions, and the curves are now - before:

ratios

after:

ratios

@astrofrog astrofrog added the bug label Mar 22, 2026
@astrofrog astrofrog marked this pull request as ready for review March 22, 2026 23:03
@astrofrog astrofrog merged commit 400e42d into astropy:main Mar 24, 2026
25 checks passed
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

Projects

None yet

Development

Successfully merging this pull request may close these issues.

1 participant