Skip to content

Fix regridding.fill() next to NaN elements, at the edges of the array, and on non-square arrays. - #71

Merged
roytsmart merged 3 commits into
mainfrom
bugfix/fill-boundaries
Oct 8, 2026
Merged

roytsmart merged 3 commits into
mainfrom
bugfix/fill-boundaries

Conversation

@roytsmart

@roytsmart roytsmart commented Oct 7, 2026 •

Copy link
Copy Markdown
Collaborator

The Gauss-Seidel method of regridding.fill() had five problems, found while filling the dust spots in IRIS slit-jaw images.

  1. NaN neighbors spread into the fill. A missing element next to a NaN element that is not in where (for example, outside an image's field of view) became NaN. The NaN then spread through the rest of the missing region. On 31 IRIS frames, 2,057 filled pixels came out NaN. Infinite elements and NaN values in guess spread the same way.
  2. The edges were periodic. The neighbors were indexed modulo the size of the array, so a missing element on one edge was filled using values from the opposite edge.
  3. Non-square arrays were weighted unequally. The grid spacing came from mapping each axis onto [-1, 1]. On a (num_y, num_x) array this weighted the two axes in the ratio ((num_x - 1) / (num_y - 1))**2: 1.12 for a 548 × 518 image and 2.4 for 548 × 845. Square arrays were unaffected.
  4. Slices of a broadcast a shared memory. a was copied before it was broadcast against where, so when a had fewer dimensions than where, every slice was filled in the same memory. The slices then overwrote each other from parallel threads, giving results off by up to 0.76.
  5. A where with fewer dimensions than a raised a FutureWarning. It was a writeable view from np.broadcast_arrays, and numba inspects that flag.

Each missing element is now replaced by the mean of its valid nearest neighbors.

  • A neighbor is valid if it is finite. Neighbors outside the array and non-finite elements not in where are left out, so the gradient of the filled elements is zero across those boundaries (a Neumann boundary).
  • Non-finite elements not in where are left as they are.
  • A missing element whose guess is not finite is valid only once its neighbors give it a value. As a result, guess=np.nan leaves elements cut off from every valid element as NaN. The default guess, the median of the known elements, leaves them at the median, and the docstring now says so.
  • The elements are assumed to be equally spaced along both axes.
  • a is copied after it is broadcast, and where is broadcast with np.broadcast_to, which returns a read-only view.
  • Validity is tracked in a separate boolean array, because the kernel is compiled with fastmath, which assumes there are no NaNs or infinities.

The number of axes is now checked before any work is done, and the two error messages are fixed: one was missing an f-prefix, the other a space.

Tests

  • test_fill_gauss_seidel_nan_neighbors: a harmonic function with zero gradient across a NaN region is recovered exactly, and the NaN region stays NaN.
  • test_fill_gauss_seidel_infinite_neighbors: the same with ±inf, which is kept as it is.
  • test_fill_gauss_seidel_edges: a strip on the edge of the array relaxes to the adjacent valid row or column, not toward the far edge. It is checked along both axes.
  • test_fill_gauss_seidel_isotropic: on a 5 × 9 array, one missing element becomes the plain mean of its four neighbors.
  • test_fill_gauss_seidel_no_valid_neighbors: an element with no valid neighbors is left at guess.
  • test_fill_gauss_seidel_broadcast: with a 2D a and a 3D where, and the other way round, each slice equals a fill of that slice alone, with warnings as errors.
  • test_fill_gauss_seidel_guess_nan: a guess array with one NaN, and guess=np.nan, both still recover a harmonic function.
  • test_fill_gauss_seidel_cut_off: elements enclosed by unknown elements are left at the median, or as NaN with guess=np.nan.
  • test_fill_gauss_seidel_axis, test_fill_method: the two errors.

All but test_fill_gauss_seidel_axis fail on main. The existing tests pass unchanged, including the exact recovery of a harmonic function inside a block.

Checked locally on the IRIS frames: none of the 36,945 filled pixels come out NaN, and the NaN outside the field is kept. 1000 iterations over 31 frames of 548 × 518 take 0.8 s.

This changes the filled values of elements within the relaxation length of an edge, and of elements in non-square arrays.

Regions of missing elements that touch an edge, or touch unknown elements, now converge more slowly than interior regions. At the default 100 iterations, a 15-column strip along the edge of a 51 × 51 array is still 0.33 from converged. The docstring now says these regions need more iterations. #72 tracks over-relaxation or a stopping tolerance, plus two other improvements from the review.

🤖 Generated with Claude Code

https://claude.ai/code/session_011VGPHLGJi8acbuLYSfBHTk

…ay, and on non-square arrays.

The Gauss-Seidel method of `regridding.fill()` now replaces each missing
element with the mean of its valid nearest neighbors.

- NaN elements which are not in `where` are no longer used as neighbors,
  so they no longer spread NaN into the filled elements, and are left as NaN.
- The edges of the array are no longer periodic.
  Neighbors outside the array are left out instead.
- The elements are assumed to be equally spaced along both axes.
  Before, the spacing came from mapping each axis onto [-1, 1],
  which weighted the two axes of a non-square array differently.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_011VGPHLGJi8acbuLYSfBHTk
@codecov

codecov Bot commented Oct 7, 2026 •

Copy link
Copy Markdown

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 99.72%. Comparing base (91ef979) to head (8c6020c).

Additional details and impacted files
@@           Coverage Diff           @@
##             main      #71   +/-   ##
=======================================
  Coverage   99.71%   99.72%           
=======================================
  Files          65       65           
  Lines        4251     4346   +95     
=======================================
+ Hits         4239     4334   +95     
  Misses         12       12           
Flag Coverage Δ
unittests 99.72% <100.00%> (+<0.01%) ⬆️

Flags with carried forward coverage won't be shown. Click here to find out more.

☔ View full report in Codecov by Harness.
📢 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.

roytsmart added a commit to sun-data/interface-region-imaging-spectrograph that referenced this pull request Oct 8, 2026
Revert this commit once a release of `regridding` includes #71.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_011VGPHLGJi8acbuLYSfBHTk
roytsmart added a commit to sun-data/interface-region-imaging-spectrograph that referenced this pull request Oct 8, 2026
Revert this commit once a release of `regridding` includes #71.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_011VGPHLGJi8acbuLYSfBHTk
… broadcast `where` without warnings.

- Neighbors are valid only if they are finite, so infinite elements that are
  not in `where` are left out like NaNs, and are left as they are.
- A missing element whose guess is not finite becomes valid only once its
  neighbors give it a value, so a NaN in `guess` no longer spreads through
  the fill, and `guess=np.nan` leaves elements cut off from every valid
  element as NaN. The default guess, the median, ignores infinite elements.
- `where` is broadcast with `np.broadcast_to`, a read-only view, so a mask
  with fewer dimensions than `a` no longer raises a FutureWarning when numba
  inspects it.
- The number of axes is checked before any work is done, and both error
  messages are fixed.
- The docstring says which elements keep the guess, and that regions that
  touch the edges need more iterations.
- New tests cover the boundary along the first axis, broadcasting `a` and
  `where` against each other, infinite neighbors, NaN guesses, elements that
  are cut off, and the two errors.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_011VGPHLGJi8acbuLYSfBHTk
The cut-off elements average each other, which rounds the median in about
a quarter of the random arrays.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_011VGPHLGJi8acbuLYSfBHTk
@roytsmart
roytsmart merged commit e3acd29 into main Oct 8, 2026
19 checks passed
roytsmart added a commit to sun-data/interface-region-imaging-spectrograph that referenced this pull request Oct 8, 2026
This reverts commit c0f52a3, which
installed `regridding` from the branch of #71 until it was released.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_011VGPHLGJi8acbuLYSfBHTk
roytsmart added a commit to sun-data/interface-region-imaging-spectrograph that referenced this pull request Oct 8, 2026
…tion.remove_dust()`. (#58)

* Find the dust on the slit-jaw CCD, and remove it with `SlitJawObservation.remove_dust()`.

`SlitJawObservation.dust` marks the pixels of each frame darkened by dust.
The dust comes from the bad-pixel map which the IRIS team measures from the
flat field of each channel and publishes in the SolarSoft database, and it
is moved into each level 2 frame the way `IRIS_DUSTBUSTER` moves it:
grown to cover the dark rim of each particle, rotated and scaled about the
center of the slit, placed relative to the slit of each frame, and moved by
the shift of up to 7 pixels which makes the dusty pixels darkest.

`SlitJawObservation.remove_dust()` replaces each dusty pixel with the median
of the same pixel in the two frames before and after it which are not dusty,
scaled by exposure time, as `IRIS_DUSTBUSTER` does, and fills the pixels no
neighboring frame sees with `regridding.fill()`.

The dust is found by default, and `dust=False` skips the download of the
bad-pixel maps.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_011VGPHLGJi8acbuLYSfBHTk

* Temporarily install `regridding` from sun-data/regridding#71.

Revert this commit once a release of `regridding` includes #71.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_011VGPHLGJi8acbuLYSfBHTk

* Pass `dust` through `iris.sji.open` when it loads the whole observation.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_011VGPHLGJi8acbuLYSfBHTk

* Require `regridding>=3.5.1`, which includes sun-data/regridding#71.

This reverts commit c0f52a3, which
installed `regridding` from the branch of #71 until it was released.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_011VGPHLGJi8acbuLYSfBHTk

---------

Co-authored-by: Claude Opus 5.5 <noreply@anthropic.com>
@roytsmart
roytsmart deleted the bugfix/fill-boundaries branch October 8, 2026 22:50
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

1 participant