Make cubic resample prefilter explicit so chunk seams stay sub-eps#1478
Make cubic resample prefilter explicit so chunk seams stay sub-eps#1478brendancol merged 3 commits intomainfrom
Conversation
…1464) scipy.ndimage.map_coordinates with order >= 2 silently runs spline_filter on a 12-pixel edge-padded copy of the input. Calling it implicitly per block (and twice per block in the NaN-aware case, once for filled and once for weights) couples the IIR boundary transient to the chunk geometry; the eager and chunked paths could disagree by ~1e-6 near seams. Run the prepad+spline_filter step ourselves, shift the sample coordinates by the pad, and pass prefilter=False to map_coordinates. Mirror the change in the cupy path. Bump cubic overlap depth from 10 to 16 so the residual transient (0.268^16 ~ 7e-10) rounds to zero in the float32 output. Add a TestCubicPrefilterParity suite that asserts atol=1e-10 between eager numpy and dask+numpy paths on polynomial, random, and NaN inputs.
|
@copilot resolve the merge conflicts in this pull request |
… prefilter + updated threshold Co-authored-by: brendancol <433221+brendancol@users.noreply.github.com>
Merge conflicts resolved in commit 5d358a4. The three conflicts in |
|
@copilot resolve the merge conflicts in this pull request |
…eep explicit prefilter Co-authored-by: brendancol <433221+brendancol@users.noreply.github.com>
Merge conflicts resolved again in commit 2d29663. The new commit |
Closes #1464.
Summary
scipy.ndimage.map_coordinates(order >= 2, prefilter=True)does three things behind the scenes formode='nearest':spline_filteron the padded array.The whole step happens inside each
map_coordinatescall. In the NaN-aware path it runs twice per block (once forfilled, once forweights). In the dask path it runs once per chunk on the overlapped block, with the prefilter's IIR boundary transient anchored to the chunk's overlapped edge rather than the global edge. Withdepth=10overlap the transient leaked into chunk-interior samples at ~1e-6 near chunk seams.This PR runs the prepad +
spline_filterstep ourselves, then callsmap_coordinates(prefilter=False). Boundary handling is now the same between numpy and dask paths, and the per-block prefilter onfilledandweightsruns once each instead of twice.The IIR transient still depends on overlap depth. The previous
_INTERP_DEPTH['cubic'] = 10left ~2e-6 residual, which is what surfaced this issue. Bumped to 16, which leaves ~7e-10 residual, below float32 epsilon, so the output rounds bit-identical between eager and chunked paths.Same change in the cupy path with
cupyx.scipy.ndimage.spline_filter.Test plan
TestCubicPrefilterParitysuite (6 tests) assertsatol=1e-10between eager numpy and dask+numpy for cubic resample of: