Preserving Seafloor Features While Smoothing #

Smoothing sits at the end of the surface smoothing workflow, and it is where good bathymetry quietly goes wrong. An isotropic Gaussian or mean filter that successfully erases the fine ripple of acquisition noise also erases the very things a hydrographic surface exists to represent: the sharp shoulders of a dredged channel, the crest of a sand ridge, the vertical wall of a wreck, the rim of a scour pit. A linear low-pass filter cannot tell noise from signal — it attenuates every high-frequency component, and real seafloor edges are high-frequency. The result passes a visual smell test but fails hydrographic fidelity: feature amplitudes shrink, positions drift, and least-depth over a hazard grows deeper than the truth. The fix is to stop using filters that treat all gradients equally and switch to edge-preserving smoothing — bilateral filtering and Perona-Malik anisotropic diffusion — tuned so that noise is averaged away while real breaks are protected.

Why Linear Low-Pass Filters Flatten Features #

A Gaussian filter convolves the surface with a fixed kernel: each output cell is a weighted average of its neighbours, and the weights depend only on distance, never on the values. Across a sharp ridge the kernel happily averages the crest with the flanks, so the crest loses height and the toe gains it — the feature is smeared symmetrically outward. The wider the kernel needed to suppress noise, the more feature amplitude it costs, and there is no kernel that removes noise without touching edges because both live in the same frequency band. The reproduction below builds a clean ridge, adds noise, and measures how much crest a Gaussian steals:

import numpy as np
from scipy.ndimage import gaussian_filter

x = np.linspace(-50, 50, 400)
ridge = 6.0 * np.exp(-(x ** 2) / (2 * 4.0 ** 2))   # a 6 m, ~4 m-wide ridge
rng = np.random.default_rng(3)
noisy = ridge + rng.normal(0.0, 0.4, size=x.shape)  # 0.4 m acquisition noise

smoothed = gaussian_filter(noisy, sigma=6.0)         # enough to kill the noise
lost = ridge.max() - smoothed.max()
print(round(float(lost), 2))  # metres of crest amplitude destroyed

The Gaussian removes the noise but shaves a large fraction off the crest — an unacceptable trade for a navigation surface where the shoalest point is safety-critical. Edge-preserving filters break the tie by making the averaging value-aware: they diffuse freely across gentle, noisy terrain and refuse to diffuse across strong gradients. This is the same discipline that makes the NaN-safe Gaussian workflow respect voids — here we extend value-awareness to feature edges as well.

Step-by-Step Fix #

1. Apply a bilateral filter #

A bilateral filter weights each neighbour by two Gaussians: one on spatial distance and one on depth difference. Neighbours that differ sharply in depth — across a channel wall — get near-zero weight, so the wall is preserved while the flat, noisy bed is smoothed:

import logging
import numpy as np
from skimage.restoration import denoise_bilateral

logger = logging.getLogger(__name__)


def bilateral_smooth(dem: np.ndarray,
                     spatial_sigma: float,
                     depth_sigma_m: float) -> np.ndarray:
    """Edge-preserving bilateral smoothing of a bathymetric grid.

    Args:
        dem: 2-D depth array (metres); must be finite (mask voids first).
        spatial_sigma: spatial Gaussian sigma in pixels (noise footprint).
        depth_sigma_m: depth-range sigma; set near the noise std so real
            steps larger than this survive.

    Returns:
        Smoothed float32 array, same shape as ``dem``.

    Raises:
        ValueError: if the array contains non-finite values or bad sigmas.
    """
    if not np.all(np.isfinite(dem)):
        raise ValueError("mask NaN/nodata to a finite value before bilateral")
    if spatial_sigma <= 0 or depth_sigma_m <= 0:
        raise ValueError("both sigmas must be > 0")

    out = denoise_bilateral(dem.astype(np.float64),
                            sigma_color=depth_sigma_m,
                            sigma_spatial=spatial_sigma,
                            channel_axis=None)
    logger.info("bilateral: spatial_sigma=%.1f px, depth_sigma=%.2f m",
                spatial_sigma, depth_sigma_m)
    return out.astype(np.float32)

The depth_sigma_m parameter is the feature/noise threshold: set it a little above the acquisition-noise standard deviation and below the smallest real relief you must keep, so a $0.4\,\mathrm{m}$ noise floor is averaged while a $6\,\mathrm{m}$ ridge is untouched.

2. Apply Perona-Malik anisotropic diffusion #

For fine control, anisotropic diffusion iteratively evolves the surface under $\partial_t I = \nabla \cdot (c(\lVert\nabla I\rVert)\,\nabla I)$, where the conductance $c$ falls toward zero at strong gradients, so diffusion halts at edges and continues on smooth ground. It preserves edges more aggressively than the bilateral filter and lets you dial the number of iterations to the noise scale:

import logging
import numpy as np

logger = logging.getLogger(__name__)


def anisotropic_diffusion(dem: np.ndarray,
                          n_iter: int,
                          kappa: float,
                          gamma: float = 0.2) -> np.ndarray:
    """Perona-Malik anisotropic diffusion for feature-preserving smoothing.

    Args:
        dem: 2-D finite depth array (metres).
        n_iter: diffusion iterations; more smoothing, more staircase risk.
        kappa: gradient threshold (metres) separating noise from edges.
        gamma: step size; keep <= 0.25 for 2-D numerical stability.

    Returns:
        Smoothed float32 array.

    Raises:
        ValueError: on non-finite input, non-positive kappa, or gamma > 0.25.
    """
    if not np.all(np.isfinite(dem)):
        raise ValueError("mask voids before diffusing")
    if kappa <= 0:
        raise ValueError("kappa must be > 0")
    if gamma > 0.25:
        raise ValueError("gamma > 0.25 is numerically unstable in 2-D")

    img = dem.astype(np.float64).copy()
    for i in range(n_iter):
        dn = np.diff(img, axis=0, prepend=img[:1, :])
        ds = -np.diff(img, axis=0, append=img[-1:, :])
        de = np.diff(img, axis=1, prepend=img[:, :1])
        dw = -np.diff(img, axis=1, append=img[:, -1:])
        # Perona-Malik exponential conductance: stops flux at strong gradients.
        cn, cs = np.exp(-(dn / kappa) ** 2), np.exp(-(ds / kappa) ** 2)
        ce, cw = np.exp(-(de / kappa) ** 2), np.exp(-(dw / kappa) ** 2)
        img += gamma * (cn * dn + cs * ds + ce * de + cw * dw)
    logger.info("anisotropic diffusion: %d iters, kappa=%.2f m", n_iter, kappa)
    return img.astype(np.float32)

Here kappa plays the role the depth sigma played above: gradients below kappa are treated as noise and smoothed, gradients above it are treated as features and preserved. Set it from the measured noise level, not by eye.

Verification and Acceptance Test #

The acceptance test must prove both halves of the trade: features survived and noise fell. Measure feature amplitude against a known reference (a benchmark tile or the raw least-depth over a charted feature) and background variance over a flat control patch, then assert both:

import logging
import numpy as np

logger = logging.getLogger(__name__)


def verify_feature_retention(reference_feature: np.ndarray,
                             smoothed_feature: np.ndarray,
                             flat_before: np.ndarray,
                             flat_after: np.ndarray,
                             amp_tol_frac: float = 0.1,
                             min_noise_drop: float = 0.5) -> dict:
    """Assert feature amplitude retained and noise variance reduced."""
    ref_amp = float(reference_feature.max() - reference_feature.min())
    got_amp = float(smoothed_feature.max() - smoothed_feature.min())
    retained = got_amp / ref_amp if ref_amp else 0.0
    if retained < 1.0 - amp_tol_frac:
        raise AssertionError(
            f"feature amplitude fell to {retained:.0%} of reference")

    var_before = float(np.var(flat_before))
    var_after = float(np.var(flat_after))
    drop = 1.0 - (var_after / var_before if var_before else 1.0)
    if drop < min_noise_drop:
        raise AssertionError(f"noise variance only dropped {drop:.0%}")

    logger.info("retained %.0f%% amplitude, noise -%.0f%%",
                100 * retained, 100 * drop)
    return {"amplitude_retained": retained, "noise_variance_drop": drop}

A filter that fails either assertion is mistuned: too little smoothing leaves the noise, too much eats the feature. Run this over each protected feature class before publishing.

Edge Cases and Gotchas #

• Over-iteration causes staircasing. Perona-Malik diffusion, pushed too far, converts smooth slopes into flat terraces separated by artificial steps (the well-known staircase effect), which a hydrographer will read as false micro-scarps. Cap n_iter, prefer a slightly larger kappa over more iterations, and inspect a slope map for banding before shipping.
• Nodata edges bias the diffusion. Both filters treat a sentinel or a NaN-filled void as a real value and let it bleed into the feature. Mask voids to a neutral value only for the duration of the filter, restore them afterward, and keep gaps intact per handling data gaps in sparse bathymetry interpolation so smoothing never invents depth across a hole.
• S-44 fidelity is the real acceptance bar. Visual smoothness is not the objective; the IHO S-44 feature-detection and total-vertical-uncertainty requirements are. Difference the smoothed surface against the source soundings over every charted feature and confirm the least depth has not deepened beyond the allowable uncertainty before the grid enters a navigation product.

Related #

• Applying Gaussian Filters to Marine DEMs — the NaN-safe isotropic baseline this page improves on for feature-rich terrain.
• Handling Data Gaps in Sparse Bathymetry Interpolation — keeping voids intact so smoothing never diffuses across them.
• Surface Smoothing Algorithms in Python — parent workflow: filter selection, kernel sizing, and chunked application.

Up: Surface Smoothing Algorithms in Python