rasterize: supercover boundary burn for all_touched=True (#2169)#2196
Merged
Conversation
…ue (#2169) The previous boundary burn drew rings with a Bresenham line, which only steps one pixel per dominant-axis step. Shallow and diagonal edges then miss the off-axis cells the segment crosses, leaving gaps along polygon outlines under all_touched=True. Switch to an Amanatides & Woo grid traversal that visits every cell a ring segment passes through. Float pixel-space endpoints are clipped to the raster with Liang-Barsky, then the walker advances one cell at a time along whichever axis hits its next boundary first. Applied to all four backends (numpy, cupy, dask+numpy, dask+cupy) so they agree pixel-for-pixel with rasterio.features.rasterize(all_touched=True).
brendancol
commented
May 20, 2026
Contributor
Author
brendancol
left a comment
There was a problem hiding this comment.
PR Review: rasterize: supercover boundary burn for all_touched=True (#2169)
Blockers (must fix before merge)
None.
Suggestions (should fix, not blocking)
None.
Nits (optional improvements)
xrspatial/rasterize.py:1467-- the GPU kernel usesINF = 1.0e308instead ofmath.inf. numba.cuda acceptsmath.infand reads clearer about intent. Not blocking; the magic constant still saturates thet_max_x < t_max_ycomparisons correctly.xrspatial/rasterize.py:780--np.array([-dx, dx, -dy, dy])builds a 4xN temporary on every call. For very large polygon counts,np.stack(..., axis=0)or computing each plane inline would shave a copy. Probably not measurable for typical polygon counts; mentioning it only because the line is repeated from the existing line-clipping path.- The new test file
test_rasterize_all_touched_supercover_2169.pyruns all four backends per case, but the cupy/dask+cupy classes silently skip on CI without a GPU. That matches existing test conventions in this repo so it is fine; just noting that GPU parity is covered locally only until CI grows a GPU runner.
What looks good
- Algorithm: Amanatides & Woo grid traversal is implemented correctly.
t_maxinitialization handles both positive and negativedx/dy, thet_deltaformula uses1/|d|, and the tie-break case advances both axes together so the supercover sweep does not over-burn diagonal-neighbour cells. - Bounds: world-space Liang-Barsky clipping happens before the pixel-space walk, so the kernel never iterates over far-away cells.
- The internal bounds check (
0 <= cy < height and 0 <= cx < width) is a sensible guard for floating-point endpoints that floor towidthorheight. max_steps = abs(end_cx - cx) + abs(end_cy - cy) + 2is a tight upper bound for the cell-walk length and guards against pathological floating-point loops.- All four backends (numpy, cupy, dask+numpy, dask+cupy) call the same float-coordinate extraction and the same supercover walker, so backend parity is structural rather than coincidental.
- The non-all-touched path (scanline fill), point burning, and line burning are unchanged. The Bresenham
_burn_lines_cpuand_burn_lines_gpukernels are still used forLineString/MultiLineStringgeometries. - Docstring update on
rasterize(... all_touched=True ...)accurately describes the new behaviour and references rasterio for parity expectations. - Test coverage: 4 polygon shapes (diagonal strip, thin triangle, edges-on-pixel-centers, concave subpixel) x 4 backends, plus a regression-specific "diagonal was undercounted before" case and a dask-matches-numpy parity test. All 18 pass locally. Broader rasterize suite (298 tests) still green.
Checklist
- Algorithm matches reference (Amanatides & Woo 1987, "A Fast Voxel Traversal Algorithm for Ray Tracing")
- All implemented backends produce consistent results
- NaN handling unchanged from existing path
- Edge cases covered by tests (subpixel vertices, edges on pixel centres, concave shapes)
- Dask chunk boundaries handled (segments clipped per-tile, parity test asserts dask == eager numpy)
- No premature materialization
- No benchmark needed (bug fix, no new public API)
- README feature matrix unchanged (rasterize was already listed)
- Docstring updated
The kernel had a local INF = 1.0e308 sentinel for the t_max / t_delta no-step case. numba.cuda lowers math.inf to the proper IEEE-754 infinity bit pattern, so swap the magic constant for math.inf and drop the local.
brendancol
commented
May 20, 2026
Contributor
Author
brendancol
left a comment
There was a problem hiding this comment.
Follow-up review
Dispositions for the previous review:
- Nit 1 (
INF = 1.0e308->math.inf): fixed in 5795f96. Tests still pass on all four backends locally. - Nit 2 (
np.array([-dx, dx, -dy, dy])4xN allocation): deferred. The same pattern lives in the older_extract_lines_vectorizedLiang-Barsky code; rewriting only the new copy at line 780 would leave the codebase inconsistent. Better handled as a follow-up sweep across both call sites. - Nit 3 (CI lacks a GPU runner for cupy/dask+cupy classes): dismissed. The test gating uses the same pattern as the rest of the repo, and CI GPU coverage is out of scope for this PR.
No new findings on the latest commit.
2 tasks
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
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.
Summary
all_touched=True, polygon boundaries were drawn with a Bresenham line, which only steps one cell per dominant-axis step. Shallow and diagonal edges then skipped the off-axis cells the segment crossed, soxrspatial.rasterize(...)returned fewer pixels thanrasterio.features.rasterize(..., all_touched=True)on the same input.all_touched=False), line burning, and point burning are unchanged.Test plan
xrspatial/tests/test_rasterize_all_touched_supercover_2169.pypins parity againstrasterio.features.rasterize(all_touched=True)for four polygons (diagonal strip, thin triangle, edges on pixel centres, concave subpixel) on all four backends. 18 tests pass locally (4 backends x 4 polygons + 2 regression cases).xrspatial/tests/test_rasterize_accuracy.pyrasterio reference tests still pass.Closes #2169