Add a polygon-clipping kernel for 2D conservative weights

[roytsmart]

Summary

A new kernel for 2D first-order conservative weights, based on polygon clipping instead of the line sweep, plus the plumbing to use it.

weights_conservative_2d sweeps the grid lines of both grids and accumulates boundary integrals. That is fully general, but it forces a sequential walk along each line and a pass over the output grid whose cost does not shrink as the input grid shrinks.

When the output grid is a uniform, axis-aligned lattice, the overlap of an input cell with an output cell can be computed directly by clipping, and the output grid is never swept. Each input cell is clipped only against the output cells its bounding box touches, so the work per cell is bounded and every cell is independent of every other.

weights(method="conservative") now selects between the two automatically: clipping when every output grid qualifies, the sweep otherwise. The gate inspects every orthogonal element, so a call uses one kernel throughout.

Performance

Through weights(), scattering a rotated and distorted scene onto a lattice, both kernels warm:

case	sweep	clipping	speedup
201 -> 512x512	0.319 s	0.041 s	7.7x
401 -> 1024x1024	1.155 s	0.148 s	7.8x
401 -> 2048x2080 (an ESIS sensor)	3.454 s	0.324 s	10.7x

Row sums are exactly 1.0 for both. Applying both operators to the same scene agrees to about 1e-5 of the peak; that difference is the perturbation, which the clipping path no longer needs (see below). At the largest size the sweep emits 6 negative weights and clipping emits none — negative weights are impossible for overlap volumes, so this is a correctness difference as well.

The crossover is well below any realistic size: an earlier probe found the sweep ahead on a 40x40 output grid, but by 512x512 clipping is already 7.7x faster and the margin grows with the output grid, since the advantage comes from never sweeping it.

End to end, through optika

Measured on the ESIS-I model: optika.systems.LinearSystem.weights() with a 401x401 scene onto a 2048x1040 sensor, 12 orthogonal elements (4 channels x 3 spectral cells).

stage	sweep	clipping
kernel (x12)	6.81 s	0.76 s
_coalesce (argsort + reduceat)	~10.2 s	~1.3 s
_normalize_input_output_coordinates (incl. perturb)	0.95 s	0.00 s
optika's distortion / vignetting / effective area	~0.55 s	~0.55 s
total	17.9 s	2.05 s

9x end to end, but note where it comes from: more of it is in _coalesce than in the kernel. _coalesce merges duplicate (input, output) pairs with an np.argsort, so its cost scales with the number of triples emitted, and the sweep emits about four times as many (7 994 348 against 1 998 583 on the synthetic case above) because it accumulates a signed contribution from each grid separately.

So roughly 6 s of the saving is the kernel, roughly 9 s is the cheaper coalesce that follows from emitting fewer fragments, and about 1 s is skipping the perturbation. The sweep cannot recover that last part on this geometry, since the sensor lattice is exactly integer-aligned and therefore degenerate (see open question 3).

The clipping kernel emits no duplicates

Each input cell visits each candidate output cell exactly once, so every (input, output) pair is emitted once and _weights_to_arrays has nothing to merge. Measured on the same ESIS case:

	build	triples
sweep, coalesced (today's default)	18.05 s	22 531 238
sweep, raw	8.20 s	90 124 984
clipping, coalesced	2.05 s	22 531 162
clipping, raw	1.16 s	22 531 162

The sweep's raw output is exactly 4.0x its coalesced output; clipping's is 1.0x. So the coalesce costs clipping 0.89 s and removes nothing, and skipping it for this kernel would take the end-to-end ratio against today's default from 9x to about 15x.

weights() now takes a coalesce keyword (default True, so existing behaviour is unchanged) to make this a choice rather than a fixed cost. It pays for itself when the weights are reused (the Level-4 MART inversion caches and re-applies them), but for single-use weights it is dead loss, since regrid_from_weights sums duplicates during the scatter-add anyway. Like-for-like with coalesce=False, the ratio here is 8.20 s against 1.16 s, or 7.1x.

coalesce=False returns the fragments as built. The result is equivalent, since regrid_from_weights sums duplicates during the scatter-add either way, and there are tests covering that both forms preserve each input cell's total weight and regrid a scene to the same answer.

Incidentally the two kernels agree on the coalesced triple count to 76 parts in 22.5 million, which is a reassuring independent check at full scale.

Two further caveats on how far this generalises. The sweep is faster than clipping on a small output grid (see the table under open question 3); the advantage here comes from never sweeping the output grid, so it grows with the number of output cells relative to input cells. And optika.systems.LinearSystem.weights(), which is the caller that motivated this, takes about 20 s on a comparable problem, but the kernel is only about 1 s of that per orthogonal element: it calls the kernel once per channel and wavelength and also evaluates its distortion, vignetting, and effective-area models. So the end-to-end speedup for that caller will be well short of 13x.

Two correctness differences

Both are covered by tests, and both are cases where the sweep gives a wrong answer that a conservation check does not catch, since row sums stay at 1.0.

1. Grids wound in the opposite sense. Reversing one axis of the vertex array does not move any cell, it only renumbers them and flips the sign of their areas, so the weights must simply permute. Against that ground truth:

clip  vs expected : max|diff| = 0.0e+00
sweep vs expected : max|diff| = 3.959e+00

2. Degenerate grids. Resampling a grid onto itself must give the identity matrix. Without perturbation:

clip  vs identity : max|diff| = 1.776e-15
sweep vs identity : max|diff| = 4.0

Clipping handles coincident vertices and collinear edges exactly, so it does not need perturb at all. That matters beyond correctness: the perturbation is a source of run-to-run variation for anything downstream that seeds the global RNG, and this kernel removes the need for it rather than inheriting it.

Implementation notes

• Signed areas throughout, which is what makes case 1 work. A cell whose edges cross each other (a "bowtie") is still unsupported, since Sutherland-Hodgman assumes a simple polygon. This is documented in the docstring.
• Each input cell writes into its own slice of the output arrays, given by a prefix sum of per-cell upper bounds. The result therefore does not depend on the thread schedule, and there is no typed.List accumulation.
• Each cell is shifted onto its own candidate block before clipping, so the coordinates are of order one rather than of order the grid width and the shoelace differences keep their significant figures. This trick is borrowed from a GPU prototype Jacob Parker wrote for the ESIS Level-4 pipeline, which uses the same clipping approach in torch.
• grid_is_uniform_rectilinear() gates the kernel and raises if the output grid does not qualify.
• The kernel choice and the perturbation are decided together. Perturbing would destroy the exact lattice the clipping kernel needs, and clipping resolves degeneracies on its own, so the coordinates are normalized without perturbation, the gate runs, and they are only normalized again with it when falling back to the sweep.
• Pure numba, no new dependencies.

regrid() is left alone deliberately: it builds weights and applies them once, so it is the clearest case for defaulting to coalesce=False, but that is a behaviour change for every existing caller and belongs with question 1 below.

Open questions for review

1. Auto-detection is now implemented; the consequence is that results move. On the clipping path the perturbation is no longer applied, so weights differ by roughly 1e-5 of the peak from what the same call returns today. Anything committed downstream that was derived through method="conservative" onto a lattice will shift by that much. Flag if you would rather this were opt-in behind a keyword for a release first.
2. The bowtie case. If self-intersecting cells need supporting, decomposing each quad into two signed triangles handles them, at roughly 1.5-2x the clipping cost.
3. The sweep's output blows up on grid-aligned geometry. Chasing a head-to-head timing turned this up. Resampling a 4x4 distorted scene onto a 40x40 lattice, with perturb off:

output lattice	sweep	clipping
40x40, integer-aligned	520 s, 1 133 237 136 triples	1.0 s, 580 triples
40x40, shifted by 0.5	0.01 s, 2 336 triples	0.02 s, 584 triples
8x8, integer-aligned	2 517 s, then MemoryError	0.004 s, 68 triples

Shifting the lattice by half a cell breaks the degeneracy and changes nothing else, and it takes the sweep from 1.1 billion triples to 2 336. The 8x8 case is the striking one: a problem whose correct answer is 68 triples ran for 42 minutes and then exhausted memory. So the sweep terminates rather than looping forever, but on exactly-aligned grids it emits roughly a million times more weight triples than the problem has, which is what produced the MemoryError I hit at larger sizes.

The public weights() enables perturb by default in 2D, so callers do not normally meet this, but perturb=False on aligned grids looks unsafe. Happy to split this out into its own issue.

Note also from that table that the sweep is faster than clipping on a small output grid. The advantage here comes from never sweeping the output grid, so it grows with the number of output cells relative to input cells and only becomes decisive on something detector-sized.

Checks

• 350 passed, 2 xfailed
• _clipping.py, _weights.py, and _weights_conservative.py all at 100% patch coverage; the sweep kernel remains at 100% too, so the new routing did not cost it coverage
• black --check, ruff check, and pyright 1.1.411 all clean

🤖 Generated with Claude Code

[codecov]

Codecov Report

❌ Patch coverage is 99.76526% with 1 line in your changes missing coverage. Please review.
✅ Project coverage is 99.51%. Comparing base (aa5d2fe) to head (7febf51).

Files with missing lines	Patch %	Lines
regridding/_regrid/_tests/test_regrid.py	95.00%	1 Missing ⚠️

Additional details and impacted files

@@           Coverage Diff            @@
##             main      #59    +/-   ##
========================================
  Coverage   99.51%   99.51%            
========================================
  Files          44       46     +2     
  Lines        2049     2472   +423     
========================================
+ Hits         2039     2460   +421     
- Misses         10       12     +2     

Flag	Coverage Δ
unittests	99.51% <99.76%> (+<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.