Fix yet another off-by-one bug - #1013
Open
astronomyk wants to merge 6 commits into
Open
Conversation
overlay_image computed the array origin at which to place small_im as ceil(coords) - small_im.shape // 2. That is only correct when coords sits on the pixel lattice small_im's shape implies: half-integer for even shapes, integer for odd ones. For anything in between it is off by up to a whole pixel, and the direction of the error flips with the parity of the shape. Chunked FOVs hit this whenever the detector footprint is not a whole number of pixels across, as MICADO's is not (189.32 mm / 0.015 mm = 12621.33 px in y). The last chunk of each axis is then a fraction of a chunk wide, so it has a different size and a pixel lattice offset from its neighbours', and it got displaced the opposite way from them. The seam between the two picked up a duplicated or a dropped row/column and the image content shifted by a pixel across it. Snap the origin instead, which is parity-independent. Half-integer origins - a genuine sub-pixel shift, which sub_pixel=True would have to handle properly - round towards the lower index, matching what the old expression did for even shapes.
extract_area_from_imagehdu slices the source array on the source's own pixel grid, growing the slice outwards to whole pixels with floor/ceil so it covers the FOV footprint. It then built the accompanying WCS with create_wcs_from_points([xy0s, xy1s]) - the world corners from *before* that growing, rounded independently. The two disagreed in both size and phase. NAXISn could be a pixel smaller than data.shape, which fits.ImageHDU silently corrects on construction, leaving CRPIXn describing an array of the wrong size; and the WCS put the cutout on the FOV's pixel lattice rather than the source's, which for chunked FOVs is offset from it by a fraction of a pixel. Projecting such a cutout mis-registered it by up to a whole pixel, by differing amounts for neighbouring FOVs, so the checkerboard phase flipped across the chunk seams and a row/column was duplicated or dropped there. The FIXME this replaces described the size half of the problem. Build the cutout WCS by shifting the input WCS's CRPIX by the slice origin instead. The cutout's pixels are input pixels, so this is exact: the cutout keeps the input's pixel lattice and NAXISn always agrees with data.shape. The added integration test reassembles a checkerboard from chunked FOVs with chunk_size=8, over four detector footprints including two that are not a whole number of pixels across. All four fail without this commit and the one before it.
This is the last of the three off-by-one causes and the most invasive of
them: it changes CRVALn for any region that is not a whole number of pixels
across. It is separable from the two commits before it - drop it and the
dropped row and the phase flips stay fixed, but MICADO keeps one duplicated
row at the last chunk boundary in y.
create_wcs_from_points rounds the requested extent to a whole NAXISn and
then centres the grid on the region, folding that region's own rounding
remainder into CRVAL. Each region therefore gets its own sub-pixel lattice.
Chunked FOVs are all carved from one origin in whole-pixel steps - the image
plane, the full-size chunks and the smaller last chunk all measure from the
same detector edge - so they ought to share a lattice, and instead ended up
offset from one another by fractions of a pixel. The projection then has to
snap them to whole pixels, and the last chunk could land a whole pixel off
its neighbours.
MICADO's detector plane is 189.32 mm / 0.015 mm = 12621.33 px tall, so it
hits this on every run: with chunk_size 2048 the last chunk in y is 333.33
px and its lattice is offset from the six full chunks'.
When NAXISn truncates the region, anchor the lower pixel edge at points.min
and let the leftover fraction fall off the far edge; regions carved from the
same origin then coincide. When NAXISn covers the region there is nothing to
truncate, so keep centring - that is what a region smaller than one pixel
needs, where NAXISn was clamped up to 1.
Verified on the full MICADO array (12621 x 12888 image plane, 4 mas
checkerboard, chunk_size 2048), checking that the board's row and column
phase alternates strictly across the whole image plane:
origin/main rows: 10407 bad, cols: 10408 bad,
plus a dropped row at y = 12288
first two commits rows: 166 bad, cols: 0 bad
with this commit rows: 0 bad, cols: 0 bad
astronomyk
added
the
bugfix
astronomyk
marked this pull request as ready for review
September 4, 2026 19:03
Codecov Report✅ All modified and coverable lines are covered by tests. Additional details and impacted files@@ Coverage Diff @@
## main #1013 +/- ##
==========================================
+ Coverage 76.35% 76.52% +0.17%
==========================================
Files 69 69
Lines 9086 9071 -15
==========================================
+ Hits 6938 6942 +4
+ Misses 2148 2129 -19 ☔ View full report in Codecov by Harness. 🚀 New features to boost your workflow:
|
FOVManager._get_splits used np.arange(v_min, v_max, step) for the chunk
boundaries. That returns ceil((v_max - v_min) / step) elements, which is one
too many whenever the division rounds just above an integer - as it does for
a span that is an exact multiple of step:
(1.2000000000000002 - -1.2000000000000002) / 0.8 == 3.0000000000000004
Whether the surplus element lands exactly on v_max or one ulp below it
depends on how the platform rounds v_min + i * step. It lands on v_max on
x86, where FovVolumeList.split then discards it (v_max <= value). On arm64
macOS it lands below v_max, which is a perfectly legal split point, so split
accepts it and carves off a volume ~1e-15 px wide. create_wcs_from_points
rounds that up to a whole pixel, so it becomes a spurious one-pixel FOV,
off the pixel lattice, adding flux at the edge of the array.
A detector footprint that is an exact multiple of chunk_size therefore got
4 extra degenerate FOVs on macOS and none on Linux or Windows. Emulating
the arm64 rounding locally, a 24 px footprint at chunk_size 8 gives:
np.arange, x86 rounding 9 volumes, none degenerate
np.arange, arm64 rounding 16 volumes, 4 of width 0
chunk_edges 9 volumes, none degenerate
Compute the edges from the chunk count instead and return only interior
points, so no edge can coincide with either bound on any platform. step
still need not divide the span or be an integer, so this is not a range().
The two test helpers added in the previous commits duplicated the arange
expression and so inherited the same fragility - that is what failed CI on
macos-latest, on exactly the exact-multiple parametrisations. They now call
chunk_edges rather than reimplementing it.
Behaviour on x86 is unchanged: the full MICADO image plane (12621 x 12888)
comes out bit-identical to the previous commit.
teutoburg
approved these changes
Sep 5, 2026
teutoburg
left a comment
teutoburg
left a comment
Contributor
There was a problem hiding this comment.
Good to have more tests about this. Did you run the checkerboard on this branch an visually check if the original artifacts are now gone @astronomyk?
Comment on lines
+238
to
+251
| _pxs = float(pixel_scale.value if isinstance(pixel_scale, u.Quantity) | ||
| else pixel_scale) | ||
| _pts = points.value if isinstance(points, u.Quantity) else np.asarray(points) | ||
| _naxis = np.asarray(naxis, dtype=float) | ||
| _extent = np.atleast_1d(np.asarray(extent, dtype=float)) | ||
|
|
||
| crval = np.where( | ||
| _naxis < _extent, | ||
| _pts.min(axis=0) + (_naxis / 2) * _pxs, # anchored | ||
| (_pts.min(axis=0) + _pts.max(axis=0) + | ||
| np.asarray(offset, dtype=float)) / 2, # centred | ||
| ) | ||
| if isinstance(points, u.Quantity): | ||
| crval = crval * points.unit |
Contributor
There was a problem hiding this comment.
This is not pretty. I tried to fix it, which led me to discover much deeper unit inconsistencies with calls to this function. It works, but everything is super fragile. I decided that rather than changing things here, I'd prefer to merge this imperfect-but-limited fix, and then do the rest separately, which will also make these lines readable again.
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
2 participants
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.
@teutoburg found that the off-by-one errors were still present when FOVs that were smaller than a full detector were projected onto the image plane. This was found by creating a alternating pixel checkerboard grid and putting it in a source object. Then running it through optical train with FOV chunks smaller than image_plane. The pattern shown below was the result.
Claude's solution is to anchor FOV chunks to their origin rather than their centres, so that we avoid the nuances of int/floor/ceil etc when working out where to tplace the chunks on the image plane
Extra testing outside the unit tests in this PR gave the results:
origin/main rows: 10407 bad, cols: 10408 bad, plus a dropped row at y = 12288
first two commits rows: 166 bad, cols: 0 bad
with this commit rows: 0 bad, cols: 0 bad