Speed up rhef via sort-and-group inner loop

[GillySpace27]

Motivation

On large solar maps ~sunkit_image.radial.rhef spends most of its time bucketing pixels into bins, not equalizing ranks. The per-bin loop rebuilds a (map_r >= lo & map_r < hi) boolean mask on every iteration, so the work scales O(N × nbins) — at a 4096² image with 2048 bins that is ~33 G boolean operations before any ranking happens. The kernel cost shows up before the user sees any output and grows with how fine they want their bins, which is exactly the wrong scaling for an algorithm whose whole point is "rank pixels per radial bin."

This PR removes that bottleneck. As the original author of rhef in sunkit-image, this had been on my back burner for a while; the downstream PUNCH analysis I drive now hit the wall hard enough to force the fix.

Summary

Replace the per-bin mask loop with a single sort-and-group pass:

1. One searchsorted on the bin lower edges gives every pixel its bin index in O(N log nbins).
2. One stable argsort groups same-bin pixels into contiguous slices.
3. The per-bin ranking step is then identical to the old code, just indexed into the sorted-slice view instead of via a boolean mask.

Equivalence

Output is byte-identical to the original implementation across both supported ranking methods (scipy, numpy) and across application_radius, upsilon, and overlapping-bin configurations. A new test_rhef_matches_reference_loop equivalence suite re-implements the old per-bin loop inline as the oracle and asserts both np.array_equal(out[finite], ref[finite]) AND np.array_equal(np.isfinite(out), np.isfinite(ref)) — the latter catching the single-corner-pixel edge case where r == edges_hi[-1] lands in no bin under < hi semantics.

Performance

Measured on a synthetic Map with method="scipy", scipy 1.13, numpy 2.0, macOS arm64:

N	nbins	old	new	speedup
2048	256	1.84 s	0.90 s	2.0×
2048	720	3.75 s	0.92 s	4.1×
2048	2048	9.10 s	0.95 s	9.6×
4096	256	5.46 s	2.51 s	2.2×
4096	720	11.26 s	2.53 s	4.5×
4096	2048	27.53 s	2.61 s	10.5×

The kernel cost no longer scales with bin count; wins grow with both image size and bin density.

cProfile on the 4096²/720 case shows the remaining ~2.5 s is dominated by find_pixel_radii / blackout_pixels_above_radius's WCS pixel-to-world lookup (~2.1 s of the 2.5 s — called twice). The inner RHEF kernel itself is ~150 ms. Caching that WCS step across the two calls is the subject of a follow-up PR (#325).

Notes

• No changes to the public API: signature, parameter defaults, units, return type, and plot_settings["norm"] side-effect are all preserved.
• The towncrier fragment changelog/324.feature.rst is included.
• The PR also carries a small drive-by fix (separate commit) pinning skimage.morphology.white_tophat's mode='reflect' argument in sunkit_image/stara.py, which was failing on py314-devdeps against the scikit-image 2.0 deprecation. The fix is conditional on the kwarg's availability so it stays compatible with the project's minimum scikit-image (0.20). Happy to split it into its own PR if maintainers prefer.

cc @sunpy/sunkit-image-maintainers