Coastal & Marine Spatial Analysis Pipelines

Comparing Median and Bilateral Filters on DEMs

Both filters are described as edge-preserving, both are one call, and they behave very differently on a seabed. A median filter is the right default for a surface with isolated bad cells and sharp steps; it is the wrong choice on a smooth slope, where it produces terracing that looks like real structure. A bilateral filter handles both but needs a noise measurement to be set correctly, and set wrongly it degrades to a Gaussian. This page picks between them on evidence, within the surface smoothing workflow.

Root cause: they preserve different things, and “edge” is ambiguous

A median replaces each cell with the median of its neighbourhood. That operation is exactly invariant to any monotonic transformation of the values, which gives it two properties that matter here: it passes a step unchanged, and it removes any feature smaller than half the window regardless of amplitude.

It also has a property that is rarely stated: on a linear ramp the median of a symmetric window is the centre value, so a pure slope passes unchanged — but on a discretised ramp, where the values are quantised or noisy, the median snaps to whichever quantisation level dominates the window, and the output is a staircase. On a bathymetric grid with 1 cm vertical resolution and a 0.3° slope, the terrace risers land every few metres and read as bedforms.

A bilateral filter weights each neighbour by two factors: distance, as a Gaussian does, and value difference, so a neighbour whose depth differs by much more than a range parameter contributes nothing. It preserves any edge larger than the range parameter, whatever its shape, and it smooths everything smaller — which is a ramp-friendly behaviour a median does not have.

import numpy as np

def terrace_metric(z: np.ndarray) -> float:
    """Fraction of cells whose value equals a neighbour's exactly.

    On a continuous surface this is near zero. A median filter drives it up,
    and the value is a direct measure of how terraced the output is.
    """
    same = np.zeros_like(z, dtype=bool)
    same[:, 1:] |= np.isclose(z[:, 1:], z[:, :-1], atol=1e-9)
    same[1:, :] |= np.isclose(z[1:, :], z[:-1, :], atol=1e-9)
    return float(same.mean())
One Profile, Three Filters, Three Different Survivors A source profile is drawn containing three features in sequence: a vertical step of about one metre, a long smooth ramp, and an isolated one-cell spike. Below it, three filtered versions are shown. The median-filtered profile keeps the step exactly, removes the spike completely, and converts the smooth ramp into a visible staircase of flat treads and abrupt risers. The bilateral-filtered profile keeps the step exactly, keeps the ramp smooth, and removes the spike because its amplitude exceeded the range sigma. The Gaussian-filtered profile rounds the step into an S-curve, keeps the ramp, and reduces but does not remove the spike, spreading it across the kernel width. source: step · ramp · spike median 7×7 step kept · spike gone · ramp terraced bilateral step kept · ramp kept · spike gone Gaussian step rounded · ramp kept · spike smeared

Step-by-step selection with production code

1. Measure the noise before choosing anything

Both filters need a scale, and the honest source for it is the surface itself over ground known to be flat:

import logging

import numpy as np

logger = logging.getLogger(__name__)


def noise_sigma_m(z: np.ndarray, flat_mask: np.ndarray) -> float:
    """Robust cell-to-cell noise scale over ground known to be flat."""
    if flat_mask.sum() < 500:
        raise ValueError("need at least 500 flat cells to estimate the noise")
    d = np.diff(z, axis=1)[flat_mask[:, 1:]]
    d = d[np.isfinite(d)]
    # MAD-based sigma: robust to the outliers this measurement is meant to survive.
    sigma = 1.4826 * float(np.median(np.abs(d - np.median(d)))) / np.sqrt(2.0)
    logger.info("noise sigma over %d flat cells: %.4f m", int(flat_mask.sum()), sigma)
    return sigma

The division by √2 turns a difference-of-two-cells spread into a per-cell one, which is the quantity both filters want. Skipping it inflates the estimate by 41% and produces a bilateral filter that smooths across genuine edges.

2. Run both on a patch that contains every feature type

The comparison is only meaningful over terrain that exercises the difference. A patch of flat mud tells you nothing:

def compare_on_patch(z: np.ndarray, sigma_m: float, *, size: int = 7) -> dict[str, float]:
    """Run both filters on the same patch and report what each did."""
    from scipy.ndimage import median_filter
    import cv2

    med = median_filter(z, size=size, mode="nearest")
    bil = cv2.bilateralFilter(z.astype("float32"), d=size,
                              sigmaColor=2.0 * sigma_m, sigmaSpace=size / 3.0)
    return {
        "median_terracing": terrace_metric(med),
        "bilateral_terracing": terrace_metric(bil),
        "median_noise_left": float(np.std(z - med)),
        "bilateral_noise_left": float(np.std(z - bil)),
        "median_max_change": float(np.nanmax(np.abs(z - med))),
        "bilateral_max_change": float(np.nanmax(np.abs(z - bil))),
    }

median_max_change is the number that catches a median filter eating a real feature. On a surface with a genuine 3 m pinnacle narrower than the window, the median’s maximum change is 3 m, and no amount of terracing statistics is as direct a statement as that.

3. Read the numbers against the terrain, not against each other

There is no universally better filter here, and the comparison output is meant to be read against what the product is for:

Terrain and product Prefer Because
Flat mud with dropouts, navigation surface median isolated bad cells vanish; nothing to terrace
Gently sloping shelf, geomorphology bilateral a median terraces the slope into false bedforms
Rock outcrops, habitat mapping bilateral preserves shoulders a median rounds when they are narrow
Dredge pocket with vertical walls median the walls are steps, which a median keeps exactly
Anything where noise σ is unknown median it needs no scale parameter and cannot be mis-tuned

That last row is the one most often decisive in practice. A bilateral filter with a wrong range sigma is a Gaussian, and a Gaussian on a rock outcrop is worse than a median on a slope.

The Bilateral Filter Has One Parameter and Three Regimes Three panels showing the same edge-and-noise profile filtered with different range sigma values relative to the measured noise sigma. At half the noise sigma, the filter's range weight rejects almost every neighbour, so the output is nearly identical to the noisy input and nothing has been smoothed. At twice the noise sigma, the flat regions are smoothed clean while the edge is preserved exactly, which is the intended behaviour. At twenty times the noise sigma, the range weight accepts every neighbour including those across the edge, so the filter has become a plain Gaussian and the edge is rounded into a ramp. A caption states that the middle regime is narrow and is found by measuring the noise, not by trying values. range σ = 0.5 × noise nothing smoothed range σ = 2 × noise flat smoothed, edge exact range σ = 20 × noise a Gaussian — the edge is a ramp now The useful regime is roughly 1.5 to 3 times the measured noise, and it is narrow. Outside it the filter does nothing or does the wrong thing, and both failures look plausible on a shaded relief. Which is why the parameter is derived from a noise measurement rather than found by trying values — a value tuned by eye on one survey does not transfer to the next, because the noise floor moves with sea state. Re-measure per survey; the measurement costs one pass over flat ground.

Verification and acceptance test 

import pytest


def test_the_chosen_filter_does_not_eat_a_known_feature(dem, pinnacle_mask, sigma_m):
    """A feature the product must keep has to survive the smoother."""
    out = smooth(dem, sigma_m)                     # whichever filter was chosen
    lost = float(np.nanmax(dem[pinnacle_mask]) - np.nanmax(out[pinnacle_mask]))
    assert lost < 0.10, f"the smoother removed {lost:.2f} m of a feature we must keep"


def test_terracing_did_not_increase_on_the_slope(dem, slope_mask, sigma_m):
    before = terrace_metric(dem[slope_mask])
    after = terrace_metric(smooth(dem, sigma_m)[slope_mask])
    assert after < before + 0.05, (
        f"terracing rose from {before:.3f} to {after:.3f} — a median is quantising the slope")

Both tests need a hand-marked mask over terrain the survey team cares about, and that mask is a five-minute job that turns filter selection from a matter of taste into a matter of record.

Cost, and where each filter belongs in the pipeline

The two filters differ by more than their edge behaviour, and on a survey-scale grid the practical differences decide as often as the theoretical ones.

Cost. A median is a rank operation and its naive cost grows with the square of the window; a 15×15 median over a 40,000 × 40,000 grid is measured in hours, not minutes. A bilateral filter is a weighted sum and, with a separable approximation, is close to linear in the window. For large windows the bilateral is faster, which inverts the usual assumption that the simpler filter is cheaper.

Determinism. Both are deterministic, but a median’s output depends on the tie-breaking rule for even-sized windows, and library defaults differ. Use odd windows and the question does not arise; use an even one and two libraries produce different surfaces from the same input.

Chunking. Both are neighbourhood operations and both need a halo of at least half the window, with the same seam consequences described in computing seafloor slope. The bilateral filter is the more dangerous of the two here, because its range weighting makes the output depend on the values in the halo as well as their positions — an under-sized halo changes the answer rather than merely truncating the kernel.

Placement. Both belong after artefact removal and before derivative computation. Smoothing before the spike removal stage spreads each spike over the kernel, turning a one-cell problem that a robust filter would have deleted into a broad bump that no filter will.

The Simpler Filter Is Not the Cheaper One A runtime curve against window size for a forty thousand by forty thousand cell grid. The median filter's curve grows quadratically, from a few minutes at a three by three window to several hours at fifteen by fifteen. The bilateral filter's curve, using a separable approximation, grows close to linearly and stays under an hour across the same range. The two cross at about a nine by nine window. An annotation notes that below the crossing the median is cheaper and above it the bilateral is, which inverts the common assumption, and that both need a halo of half the window when chunked. window size (cells) runtime median — grows with window² bilateral, separable — near linear they cross near 9×9 391521 Both need a halo of half the window when chunked — and the bilateral's answer, not just its edge, depends on it.

Edge cases and gotchas

  • cv2.bilateralFilter requires float32 and no NaNs. It silently treats NaN as a value and propagates it across the whole kernel. Fill with a normalised convolution first, or mask afterwards using the pre-filter validity — the approach in applying Gaussian filters.
  • A median window must be odd, and its cost grows fast. scipy.ndimage.median_filter is O(window²) per cell without an optimised kernel; a 15×15 median over a large grid is slow enough to change how the pipeline is scheduled.
  • Neither filter respects the grid’s uncertainty band. Both weight every valid cell equally regardless of how many soundings produced it. Where a per-cell count exists, weighting by it before smoothing is a small change with a large effect at the survey edges.

Up: Surface Smoothing Algorithms in Python