Add bilateral filtering for feature-preserving smoothing

[brendancol]

Author of Proposal: @brendan

Reason or problem

mean() blurs everything, edges included. Run it over a DEM and ridgelines soften, cliff faces round off, small channels vanish. There's currently no way to smooth noise while keeping sharp features.

Proposal

Add a bilateral() function for bilateral filtering on 2D rasters. The filter weights neighbors by spatial distance and value similarity, so pixels across an edge barely contribute. Flat areas get smoothed; edges stay sharp.

Design:

Two parameters beyond the input raster:

• sigma_spatial (float, default 1.0): Spatial Gaussian falloff. Larger = wider neighborhood.
• sigma_range (float, default 10.0): Range (value-similarity) Gaussian falloff. Larger = more smoothing across value jumps; smaller = sharper edge preservation.

Kernel radius is ceil(2 * sigma_spatial), so the window is 2*radius + 1 pixels wide.

Per-pixel algorithm:

1. Iterate over neighbors within kernel radius
2. spatial_weight = exp(-dist^2 / (2 * sigma_spatial^2))
3. range_weight = exp(-(value_diff)^2 / (2 * sigma_range^2))
4. weight = spatial_weight * range_weight
5. Output = weighted average of neighbor values

NaN pixels are skipped. If all neighbors are NaN, output is NaN.

Four backends via ArrayTypeFunctionMapping:

• NumPy: @ngjit Numba kernel over the 2D grid
• CuPy: @cuda.jit CUDA kernel, one thread per pixel
• Dask+NumPy: map_overlap wrapping the NumPy kernel with depth=radius
• Dask+CuPy: map_overlap wrapping the CuPy kernel with depth=radius

Supports boundary modes (nan, nearest, reflect, wrap), same as the rest of the focal API.

Usage:

from xrspatial import bilateral

# Light smoothing, strong edge preservation
smoothed = bilateral(dem, sigma_spatial=1.0, sigma_range=5.0)

# Heavier smoothing, moderate edge preservation
smoothed = bilateral(dem, sigma_spatial=2.0, sigma_range=20.0)

# Accessor style
smoothed = dem.xrs.bilateral(sigma_spatial=1.5, sigma_range=10.0)

Value: Right now it's either no smoothing or blur everything. Bilateral filtering is standard in image processing and GIS, and it's missing here.

Stakeholders and impacts

Useful for terrain analysis and raster preprocessing. Adds a new module (bilateral.py), doesn't touch existing functions.

Drawbacks

O(n * k^2) where k is kernel width. Large sigma_spatial will be slow, but that's the nature of the algorithm and consistent with other focal ops.

Alternatives

• Median filter: some edge preservation but different characteristics, harder to JIT efficiently.
• Anisotropic diffusion: iterative, more complex, better as a separate feature.
• Gaussian blur via convolve_2d(): already possible but doesn't preserve edges.

Unresolved questions

None. The bilateral filter is well-established and fits the existing focal op pattern.

Additional notes

Reference: Tomasi & Manduchi, "Bilateral Filtering for Gray and Color Images," ICCV 1998.