Pin edge_detection against scipy.ndimage with golden-value tests

[brendancol]

Reason or Problem

test_edge_detection.py checks properties (constant Sobel response on a ramp, zero Laplacian on a constant field) and cross-backend agreement, but never compares against an independent reference. A regression in the kernels or in convolve_2d that changed the numbers uniformly would keep all backends agreeing with each other and could slip past the property checks.

I ran all five functions (sobel_x, sobel_y, prewitt_x, prewitt_y, laplacian) against scipy.ndimage 1.16.1 on three synthetic inputs: a 64x64 Gaussian hill, a tilted plane with known analytic gradient, and a seeded random-walk surface. Results:

• interior pixels match ndi.sobel / ndi.prewitt / ndi.laplace to ~1e-13, and match ndi.correlate with the same 3x3 stencils bit for bit
• the tilted plane matches the analytic response exactly (Sobel = 8x the per-cell gradient, Prewitt = 6x, Laplacian = 0)
• boundary='nearest', 'reflect', and 'wrap' match the scipy modes of the same names on every pixel, edges included
• cupy, dask+numpy, and dask+cupy outputs are bit-identical to numpy

The parity holds today, but nothing in the test suite would catch it breaking.

Proposal

Add a golden-value test class to xrspatial/tests/test_edge_detection.py: a small 5x6 integer-valued input with hardcoded expected arrays generated from scipy.ndimage 1.16.1 (default mode='reflect'), asserted against the boundary='reflect' full output and the default boundary='nan' interior. Test-only change, no source edits.

Comparison script used for the validation

"""Reference validation: xrspatial.edge_detection vs scipy.ndimage.

Compares sobel_x/sobel_y/prewitt_x/prewitt_y/laplacian against
scipy.ndimage.sobel/prewitt/laplace (and ndi.correlate with the same
stencils) on three synthetic inputs. Interior pixels compared first,
edge rows/cols separately per boundary mode. Also checks the tilted
plane against the analytic response (Sobel=8*g, Prewitt=6*g, Lap=0)
and cross-backend parity (cupy, dask+numpy, dask+cupy vs numpy).
"""
import sys

import numpy as np
import scipy
import scipy.ndimage as ndi

import xrspatial
from xrspatial.edge_detection import (
    LAPLACIAN_KERNEL, PREWITT_X, PREWITT_Y, SOBEL_X, SOBEL_Y,
    laplacian, prewitt_x, prewitt_y, sobel_x, sobel_y,
)
from xrspatial.tests.general_checks import create_test_raster

print("xrspatial:", xrspatial.__file__)
print("scipy:", scipy.__version__)

RTOL = 1e-6
ATOL = 1e-9

# ---------------------------------------------------------------- inputs
rng = np.random.RandomState(20260718)

# 1) 64x64 Gaussian hill, 30 m cellsize (coords realistic; operator is
#    index-space so cellsize must not affect the result)
n = 64
ii, jj = np.mgrid[0:n, 0:n]
hill = 100.0 * np.exp(-(((ii - 32) ** 2 + (jj - 32) ** 2) / (2 * 10.0 ** 2)))

# 2) tilted plane, per-cell gradients gx (along columns), gy (along rows)
gx, gy = 2.5, -1.7
plane = gx * jj.astype(np.float64) + gy * ii.astype(np.float64)

# 3) rough surface: seeded 2D random walk (cumsum both axes)
rough = np.cumsum(np.cumsum(rng.randn(n, n), axis=0), axis=1)

INPUTS = {"hill": hill, "plane": plane, "rough": rough}

FUNCS = {
    "sobel_x": (sobel_x, SOBEL_X, lambda d, m: ndi.sobel(d, axis=1, mode=m)),
    "sobel_y": (sobel_y, SOBEL_Y, lambda d, m: ndi.sobel(d, axis=0, mode=m)),
    "prewitt_x": (prewitt_x, PREWITT_X, lambda d, m: ndi.prewitt(d, axis=1, mode=m)),
    "prewitt_y": (prewitt_y, PREWITT_Y, lambda d, m: ndi.prewitt(d, axis=0, mode=m)),
    "laplacian": (laplacian, LAPLACIAN_KERNEL, lambda d, m: ndi.laplace(d, mode=m)),
}

# xrspatial boundary -> candidate scipy modes (checked empirically on edges)
BOUNDARY_MODES = {
    "nearest": ["nearest"],
    "reflect": ["reflect", "mirror"],
    "wrap": ["wrap"],
}

failures = []


def check(label, a, b, rtol=RTOL, atol=ATOL):
    a = np.asarray(a, dtype=np.float64)
    b = np.asarray(b, dtype=np.float64)
    ok = np.allclose(a, b, rtol=rtol, atol=atol, equal_nan=True)
    maxdiff = float(np.nanmax(np.abs(a - b))) if a.size else 0.0
    status = "OK  " if ok else "FAIL"
    print(f"  {status} {label}: max|diff|={maxdiff:.3e}")
    if not ok:
        failures.append(label)
    return ok


# ------------------------------------------------- interior vs scipy
print("\n=== interior pixels (boundary='nan', vs scipy mode='mirror') ===")
for iname, data in INPUTS.items():
    agg = create_test_raster(data, backend="numpy", attrs={"res": (30, 30)})
    print(f"[{iname}]")
    for fname, (func, kern, sfunc) in FUNCS.items():
        r = func(agg).data
        ref = sfunc(data, "mirror")          # scipy named operator
        ref2 = ndi.correlate(data, kern, mode="mirror")  # direct stencil
        check(f"{iname}/{fname} interior vs ndi.{fname.split('_')[0]}",
              r[1:-1, 1:-1], ref[1:-1, 1:-1])
        check(f"{iname}/{fname} interior vs ndi.correlate(kernel)",
              r[1:-1, 1:-1], ref2[1:-1, 1:-1])

# ------------------------------------------------- analytic ground truth
print("\n=== analytic ground truth: tilted plane ===")
agg = create_test_raster(plane, backend="numpy", attrs={"res": (30, 30)})
expected = {
    "sobel_x": 8 * gx, "sobel_y": 8 * gy,
    "prewitt_x": 6 * gx, "prewitt_y": 6 * gy,
    "laplacian": 0.0,
}
for fname, (func, _, _) in FUNCS.items():
    r = func(agg).data[1:-1, 1:-1]
    check(f"plane/{fname} analytic={expected[fname]}",
          r, np.full_like(r, expected[fname]), atol=1e-9)

# ------------------------------------------------- edges per boundary mode
print("\n=== full array incl. edges, per boundary mode ===")
for iname, data in INPUTS.items():
    agg = create_test_raster(data, backend="numpy", attrs={"res": (30, 30)})
    for bmode, scipy_modes in BOUNDARY_MODES.items():
        for fname, (func, kern, _) in FUNCS.items():
            r = func(agg, boundary=bmode).data
            matched = None
            diffs = {}
            for smode in scipy_modes:
                ref = ndi.correlate(data, kern, mode=smode)
                d = float(np.max(np.abs(r - ref)))
                diffs[smode] = d
                if np.allclose(r, ref, rtol=RTOL, atol=ATOL):
                    matched = smode
                    break
            if matched:
                print(f"  OK   {iname}/{fname} boundary={bmode} == scipy mode={matched}")
            else:
                print(f"  FAIL {iname}/{fname} boundary={bmode}: diffs={diffs}")
                failures.append(f"{iname}/{fname}/boundary={bmode}")

# ------------------------------------------------- coords must not matter
print("\n=== cellsize/coords independence (30m vs 0.5m coords) ===")
agg30 = create_test_raster(hill, backend="numpy", attrs={"res": (30, 30)})
agg05 = create_test_raster(hill, backend="numpy", attrs={"res": (0.5, 0.5)})
for fname, (func, _, _) in FUNCS.items():
    check(f"hill/{fname} res=30 vs res=0.5",
          func(agg30).data, func(agg05).data)

# ------------------------------------------------- backend parity
print("\n=== backend parity vs numpy (hill + rough, boundary nan & reflect) ===")
BACKENDS = ["dask+numpy", "cupy", "dask+cupy"]
for iname in ["hill", "rough"]:
    data = INPUTS[iname]
    np_agg = create_test_raster(data, backend="numpy", attrs={"res": (30, 30)})
    for backend in BACKENDS:
        b_agg = create_test_raster(data, backend=backend,
                                   attrs={"res": (30, 30)}, chunks=(25, 25))
        for fname, (func, _, _) in FUNCS.items():
            for bmode in ["nan", "reflect"]:
                r_np = func(np_agg, boundary=bmode).data
                r_b = func(b_agg, boundary=bmode).data
                if hasattr(r_b, "compute"):
                    r_b = r_b.compute()
                if "cupy" in str(type(r_b)).lower() or hasattr(r_b, "get"):
                    r_b = r_b.get()
                check(f"{iname}/{fname}/{backend}/boundary={bmode}", r_np, r_b)

print("\n==============================")
if failures:
    print(f"{len(failures)} FAILURES:")
    for f in failures:
        print(" -", f)
    sys.exit(1)
print("ALL COMPARISONS PASSED")