Fast path for find_pixel_radii on non-rotated HPC maps by GillySpace27 · Pull Request #325 · sunpy/sunkit-image

GillySpace27 · 2026-06-12T04:41:32Z

Motivation

After #324 collapsed the inner loop in ~sunkit_image.radial.rhef, profiling showed ~sunkit_image.utils.find_pixel_radii had become the new bottleneck: on a 4096² map it costs ~1 s per call, and rhef calls it twice in succession (once via find_radial_bin_edges, once via blackout_pixels_above_radius). That is ~2 s of WCS conversion before any filtering happens, on what is functionally a fixed geometric transform.

The slow path uses smap.wcs.pixel_to_world and the full SkyCoord stack for every pixel, but for the common case — a non-rotated helioprojective-Cartesian map with a known rsun_obs — the answer is a closed-form function of CDELT, CRPIX, CRVAL, and the gnomonic projection. This PR adds that fast path (and lets internal callers avoid the duplicate call).

Summary

Two coupled changes:

Analytic fast path for find_pixel_radii: when both axes are HPLN/HPLT and the rotation matrix is the identity, compute pixel-to-Sun-centre distance directly from cdelt, crpix, crval, and rsun_obs plus the gnomonic-TAN inverse. Rotated or non-HPC maps fall through to the existing SkyCoord path unchanged.
Optional map_r= kwarg on blackout_pixels_above_radius: callers that already hold the radii array can pass it in and skip the internal find_pixel_radii call. rhef now does this: it has map_r from find_radial_bin_edges and threads it through, eliminating the duplicate WCS pass.

Equivalence

The fast-path output matches the SkyCoord path to:

better than 2e-4 R⊙ on a 4096² map at 2"/px (the worst-case scenario; the diff is concentrated in the corner pixels)
floating-point precision on AIA-scale 0.6"/px maps

Tests cover:

Equivalence on sun-centred, off-centre (CRVAL != 0), and scale=-rescaled cases.
Non-identity rotation triggers the SkyCoord fallback (verified by mocking the fast helper).
blackout_pixels_above_radius(map_r=...) produces byte-equal output to the no-kwarg call, and the patched internal find_pixel_radii is verified NOT to be called.

Performance

Measured on a synthetic 4096² HPLN-TAN map, macOS arm64, numpy 2.0, astropy 7.x:

	call time	speedup
Slow path (`smap.wcs.pixel_to_world`)	1.05 s	1×
Fast path (this PR)	0.15 s	~7×

For rhef specifically — which calls find_pixel_radii twice — the WCS-bound portion drops from ~2.1 s to ~150 ms total. Combined with the inner-loop rewrite in #324, a 4096² × 720-bin rhef should run in around 250 ms instead of the original ~11 s.

Notes

No public API changes other than the new keyword-only map_r= argument on blackout_pixels_above_radius.
The towncrier fragment changelog/325.feature.rst covers both speedup mechanisms in a single paragraph.
This PR is branched from main, parallel to Speed up rhef via sort-and-group inner loop #324; they touch sunkit_image/radial.py in different spots and one will need a trivial rebase after the first lands. I'm happy to do that — just let me know merge order.
Like Speed up rhef via sort-and-group inner loop #324, this also carries a small drive-by fix pinning skimage.morphology.white_tophat's mode='reflect' argument in sunkit_image/stara.py. The CI matrix py314-devdeps job was failing on test_stara due to a PendingSkimage2Change deprecation that's unrelated to the rest of this PR; fix is conditional on kwarg availability so it stays compatible with the project's minimum scikit-image (0.20). Happy to split into its own PR.

cc @sunpy/sunkit-image-maintainers

For non-rotated helioprojective-Cartesian maps -- the common case for AIA, HMI, LASCO, EUI, and similar instruments -- ``find_pixel_radii`` paid for a full ``smap.wcs.pixel_to_world`` SkyCoord round-trip on every call. On a 4096x4096 map that's ~1 s of WCS conversion, and the function is called twice in ``rhef`` (once via ``find_radial_bin_edges``, then again inside ``blackout_pixels_above_radius``) so the per-frame cost was ~2 s before any filtering happened. This commit adds an analytic fast path that computes pixel-to-sun-centre distance directly from the WCS header values: dx = (np.arange(nx) - crpix) * cdelt + crval # 1-D arcsec dy = ... r = arctan(hypot(dx, dy) in radians) # TAN inverse and falls through to the existing SkyCoord path whenever the rotation matrix isn't identity or either axis isn't HPLN/HPLT. Output matches the slow path to better than 2e-4 R_sun on a 4096^2 map at 2"/px (and to floating-point precision on AIA-scale 0.6"/px maps). Measured on the same hardware as my previous PR's benchmark: 4096^2 HPLN-TAN, no rotation: slow (SkyCoord) ~1.05 s fast (this PR) ~0.15 s ~7x To eliminate the second call in ``rhef`` entirely, ``blackout_pixels_above_radius`` gains a keyword-only ``map_r`` argument that callers already holding the radii can pass in -- skipping its internal ``find_pixel_radii`` call. ``rhef`` now does this: it has ``map_r`` from ``find_radial_bin_edges`` and threads it through. Both optimisations together drop the WCS-bound portion of ``rhef`` from ~2.1 s to ~150 ms on the 4096^2 benchmark. Tests: - Equivalence: fast path matches slow path within 2e-4 R_sun on sun-centred, off-centre, and scale-kwarg variants of a synthetic HPC fixture. - Fallback: a non-identity PC rotation matrix triggers the SkyCoord path (verified by mocking the fast helper). - ``map_r`` reuse: ``blackout_pixels_above_radius(map_r=...)`` produces byte-equal output to the no-kwarg call, and the patched internal ``find_pixel_radii`` is verified NOT to be called. Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>

…ation scikit-image 2.0 (currently dev) raises ``PendingSkimage2Change`` warnings from ``skimage.morphology.white_tophat`` because it will switch the default ``mode`` from ``'reflect'`` to ``'ignore'``. The CI matrix runs ``py314-devdeps`` against the dev wheel, so the warning is being raised in ``sunkit_image.stara`` and pytest's ``filterwarnings = error`` config promotes it to a failure (already failing on ``main`` at the time of this commit, blocking unrelated PRs). Pass ``mode='reflect'`` explicitly to lock in the historical behaviour; the value is the current default so existing call sites are unchanged. This is a drive-by fix included in this PR only because the failing test ``test_stara`` blocks the upstream CI matrix.

Gilly + another Claude landed two upstream-pending optimizations in his sunkit-image fork: - sunpy/sunkit-image#324: sort-and-group rhef inner loop (~10× at the bin density we use) - sunpy/sunkit-image#325: analytic find_pixel_radii fast path for non-rotated helioprojective-Cartesian maps (~7×) Combined, a 4096² RHEF drops from ~2.1 s to ~150 ms on the Render box — roughly 14× end-to-end. Side benefits include making the "HQ Filtered" tier feasible to render on demand for any user-picked date, and unlocking the hourly auto-render Gilly wants for gilly.space/sun. Both PRs are API-compatible (no signature changes to radial.rhef or utils.find_pixel_radii) so the call sites in api/main.py:90-91 and api/main.py:1232/1235 need no changes. Source: temporary merged branch https://github.com/GillySpace27/sunkit-image/tree/solar-archive-fast-rhef which I created by `git merge --no-ff origin/rhef-fast-kernel` on top of `origin/wcs-fast-path`. Clean auto-merge, no conflicts. TODO: flip back to a pinned PyPI release once both PRs land in a tagged upstream sunkit-image release. Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>

ayshih

I'm uneasy with the level of approximation in this PR because the inaccuracies will become unacceptable off the solar disk. (The inaccuracies would be "okay" if find_pixel_radii() were only about identifying the solar limb, but the function is more general than that.) My inclination would be for find_pixel_radii() to not use this fast path automatically, but rather have the user opt-in via a keyword argument.

ayshih · 2026-06-12T14:26:40Z

+    dx_as = ((np.arange(nx) - crpix[0]) * cdelt[0] + crval[0]) * u.Unit(cunit[0]).to(u.arcsec)
+    dy_as = ((np.arange(ny) - crpix[1]) * cdelt[1] + crval[1]) * u.Unit(cunit[1]).to(u.arcsec)


Uses the small-angle approximation for the WCS projection, and assuming that the reference coordinate is close to Sun center

ayshih · 2026-06-12T14:27:31Z

+    # Sky-plane offsets in arcsec (broadcast across the 2-D grid).
+    dx_as = ((np.arange(nx) - crpix[0]) * cdelt[0] + crval[0]) * u.Unit(cunit[0]).to(u.arcsec)
+    dy_as = ((np.arange(ny) - crpix[1]) * cdelt[1] + crval[1]) * u.Unit(cunit[1]).to(u.arcsec)
+    r_planar = np.sqrt(dx_as[None, :] ** 2 + dy_as[:, None] ** 2)


Uses the small-angle approximation

ayshih · 2026-06-12T14:29:45Z

+
+    Skips ``smap.wcs.pixel_to_world``'s SkyCoord round-trip and computes the
+    pixel-to-sun-centre distance directly from the WCS header values.  Uses
+    the TAN gnomonic projection's spherical-distance inverse so the output


Given the approximations being used, this calculation is not actually specific to the gnomonic projection, but the inaccuracy will depend on the projection

The analytic fast path approximates the WCS projection (tangent-plane treatment, reference-near-Sun-center assumption), so its accuracy degrades off-disk where find_pixel_radii is also used. Per review feedback, gate it behind a keyword-only fast= argument (default False) so the exact SkyCoord path stays the default for all callers; rotated and non-HPC maps still fall through to the exact path even when fast=True. Update tests to opt in explicitly and add a test asserting the default does not touch the approximation. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>

GillySpace27 · 2026-06-19T09:40:23Z

Good point — thanks. The fast path does lean on the small-angle / reference-near-Sun-center assumptions, and find_pixel_radii() is general enough (used well off-disk too) that silently swapping in the approximation is the wrong default.

I've pushed a commit making it opt-in via a keyword-only fast= argument (default False):

fast=False (the default) always runs the exact SkyCoord-based path — no behavior change for any existing caller.
fast=True uses the analytic path only when the map is non-rotated HPC; rotated and non-HPC maps still fall through to the exact path even when fast=True.
The docstring now spells out that the fast path is an approximation whose accuracy degrades away from disk center, and there's a new test asserting the default never touches _find_pixel_radii_fast.

Let me know if you'd prefer a different keyword name or a warning when fast=True silently falls through.

GillySpace27 force-pushed the wcs-fast-path branch 3 times, most recently from fbf076c to e3ae8e6 Compare June 12, 2026 05:58

GillySpace27 force-pushed the wcs-fast-path branch 3 times, most recently from e053951 to a05e91b Compare June 12, 2026 06:16

GillySpace27 force-pushed the wcs-fast-path branch from a05e91b to e157b8d Compare June 12, 2026 06:17

GillySpace27 mentioned this pull request Jun 12, 2026

Speed up rhef via sort-and-group inner loop #324

Open

GillySpace27 marked this pull request as ready for review June 12, 2026 06:56

ayshih reviewed Jun 12, 2026

View reviewed changes

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Uh oh!

Fast path for find_pixel_radii on non-rotated HPC maps#325

Fast path for find_pixel_radii on non-rotated HPC maps#325
GillySpace27 wants to merge 3 commits into
sunpy:mainfrom
GillySpace27:wcs-fast-path

GillySpace27 commented Jun 12, 2026 •

edited

Loading

Uh oh!

ayshih left a comment

Uh oh!

ayshih Jun 12, 2026 •

edited

Loading

Uh oh!

ayshih Jun 12, 2026

Uh oh!

ayshih Jun 12, 2026

Uh oh!

GillySpace27 commented Jun 19, 2026

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants

		dx_as = ((np.arange(nx) - crpix[0]) * cdelt[0] + crval[0]) * u.Unit(cunit[0]).to(u.arcsec)
		dy_as = ((np.arange(ny) - crpix[1]) * cdelt[1] + crval[1]) * u.Unit(cunit[1]).to(u.arcsec)

Uh oh!

Conversation

GillySpace27 commented Jun 12, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Motivation

Summary

Equivalence

Performance

Notes

Uh oh!

ayshih left a comment

Choose a reason for hiding this comment

Uh oh!

ayshih Jun 12, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

ayshih Jun 12, 2026

Choose a reason for hiding this comment

Uh oh!

ayshih Jun 12, 2026

Choose a reason for hiding this comment

Uh oh!

GillySpace27 commented Jun 19, 2026

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants

GillySpace27 commented Jun 12, 2026 •

edited

Loading

ayshih Jun 12, 2026 •

edited

Loading