diff --git a/development/utils/_mesh_smoothing_reference.py b/development/utils/_mesh_smoothing_reference.py
new file mode 100644
index 0000000..2ff5323
--- /dev/null
+++ b/development/utils/_mesh_smoothing_reference.py
@@ -0,0 +1,42 @@
+from typing import Tuple
+
+import nifty
+import numpy as np
+
+
+def smooth_mesh(
+    verts: np.ndarray, normals: np.ndarray, faces: np.ndarray, iterations: int
+) -> Tuple[np.ndarray, np.ndarray]:
+    """Smooth mesh surface via laplacian smoothing.
+
+    Args:
+        verts: The mesh vertices.
+        normals: The mesh normals.
+        faces: The mesh faces.
+        iterations: The number of smoothing iterations.
+
+    Returns:
+        The vertices after smoothing.
+        The normals after smoothing.
+    """
+    n_verts = len(verts)
+    g = nifty.graph.undirectedGraph(n_verts)
+
+    edges = np.concatenate([faces[:, :2], faces[:, 1:], faces[:, ::2]], axis=0)
+    g.insertEdges(edges)
+
+    current_verts = verts
+    current_normals = normals
+    new_verts = np.zeros_like(verts, dtype=verts.dtype)
+    new_normals = np.zeros_like(normals, dtype=normals.dtype)
+
+    # Implement this directly in nifty for speed up?
+    for it in range(iterations):
+        for vert in range(n_verts):
+            nbrs = np.array([vert] + [nbr[0] for nbr in g.nodeAdjacency(vert)], dtype="int")
+            new_verts[vert] = np.mean(current_verts[nbrs], axis=0)
+            new_normals[vert] = np.mean(current_normals[nbrs], axis=0)
+        current_verts = new_verts
+        current_normals = new_normals
+
+    return new_verts, new_normals
diff --git a/development/utils/benchmark_mesh_smoothing.py b/development/utils/benchmark_mesh_smoothing.py
new file mode 100644
index 0000000..3dd0d6a
--- /dev/null
+++ b/development/utils/benchmark_mesh_smoothing.py
@@ -0,0 +1,177 @@
+"""Benchmark Laplacian mesh smoothing: bioimage_cpp vs the nifty-based reference.
+
+The mesh is built by triangulating ``n`` random points sampled on the unit
+sphere via ``scipy.spatial.ConvexHull`` — a closed manifold of ``2n - 4``
+triangles is a fair workload for the algorithm. The Python reference is
+imported from ``_mesh_smoothing_reference.py`` and uses ``nifty`` for graph
+adjacency; if either library is missing, the corresponding column is skipped.
+
+Run::
+
+    python development/utils/benchmark_mesh_smoothing.py --n-vertices 5000 --iterations 5 --repeats 3
+    python development/utils/benchmark_mesh_smoothing.py --n-vertices 20000 --threads 1,2,4,0
+
+Notes
+-----
+The reference does in-place Gauss-Seidel smoothing after the first iteration
+(an aliasing accident in the Python source); bioimage_cpp does textbook Jacobi
+smoothing. Correctness against the reference is only checked for
+``iterations=1``, which is where the two implementations agree exactly.
+"""
+
+from __future__ import annotations
+
+import argparse
+import csv
+import os
+import sys
+from dataclasses import dataclass
+from statistics import median
+from time import perf_counter
+
+import numpy as np
+
+import bioimage_cpp as bic
+
+
+def parse_args() -> argparse.Namespace:
+    parser = argparse.ArgumentParser(description=__doc__.splitlines()[0])
+    parser.add_argument("--n-vertices", type=int, default=5000)
+    parser.add_argument("--iterations", type=int, default=5)
+    parser.add_argument("--repeats", type=int, default=3)
+    parser.add_argument("--warmup", type=int, default=1)
+    parser.add_argument("--seed", type=int, default=2026)
+    parser.add_argument(
+        "--threads",
+        type=str,
+        default="1,0",
+        help="Comma-separated list of n_threads values to benchmark (0=auto).",
+    )
+    parser.add_argument("--csv", type=str, default=None, help="Optional CSV output path.")
+    return parser.parse_args()
+
+
+@dataclass(frozen=True)
+class Mesh:
+    verts: np.ndarray
+    normals: np.ndarray
+    faces: np.ndarray
+
+
+def build_sphere_mesh(n_vertices: int, seed: int) -> Mesh:
+    from scipy.spatial import ConvexHull
+
+    rng = np.random.default_rng(seed)
+    raw = rng.standard_normal((n_vertices, 3))
+    points = (raw / np.linalg.norm(raw, axis=1, keepdims=True)).astype(np.float64)
+    hull = ConvexHull(points)
+    return Mesh(verts=points, normals=points.copy(), faces=hull.simplices.astype(np.int64))
+
+
+def time_call(fn, repeats: int, warmup: int) -> list[float]:
+    for _ in range(warmup):
+        fn()
+    times = []
+    for _ in range(repeats):
+        t0 = perf_counter()
+        fn()
+        times.append(perf_counter() - t0)
+    return times
+
+
+def import_reference():
+    here = os.path.dirname(os.path.abspath(__file__))
+    if here not in sys.path:
+        sys.path.insert(0, here)
+    try:
+        import nifty  # noqa: F401
+    except ImportError:
+        return None
+    try:
+        from _mesh_smoothing_reference import smooth_mesh as ref
+    except ImportError:
+        return None
+    return ref
+
+
+def main() -> int:
+    args = parse_args()
+    thread_values = [int(t) for t in args.threads.split(",") if t.strip()]
+
+    try:
+        mesh = build_sphere_mesh(args.n_vertices, args.seed)
+    except ImportError:
+        print("scipy is required to build the benchmark mesh.")
+        return 1
+
+    n_faces = mesh.faces.shape[0]
+    print(
+        f"Mesh: {args.n_vertices} vertices, {n_faces} faces "
+        f"({args.iterations} smoothing iterations, {args.repeats} repeats)"
+    )
+
+    rows: list[dict] = []
+
+    # bioimage_cpp variants (one per --threads value).
+    for n_threads in thread_values:
+        label = f"bioimage_cpp[n_threads={n_threads}]"
+
+        def run(verts=mesh.verts, normals=mesh.normals, faces=mesh.faces, n=n_threads):
+            bic.utils.smooth_mesh(verts, normals, faces, iterations=args.iterations, n_threads=n)
+
+        times = time_call(run, args.repeats, args.warmup)
+        rows.append(
+            {"impl": label, "median_s": median(times), "best_s": min(times)}
+        )
+
+    # Reference (single-threaded by definition).
+    reference = import_reference()
+    if reference is not None:
+        def run_ref():
+            reference(mesh.verts, mesh.normals, mesh.faces, args.iterations)
+
+        times = time_call(run_ref, args.repeats, args.warmup)
+        rows.append(
+            {"impl": "reference[nifty]", "median_s": median(times), "best_s": min(times)}
+        )
+    else:
+        print("(reference skipped: nifty not installed)")
+
+    # Correctness check against the reference at iterations=1 (only point of
+    # exact agreement; see module docstring).
+    if reference is not None:
+        ref_v, ref_n = reference(mesh.verts, mesh.normals, mesh.faces, 1)
+        ours_v, ours_n = bic.utils.smooth_mesh(mesh.verts, mesh.normals, mesh.faces, iterations=1)
+        v_max_diff = float(np.max(np.abs(ours_v - ref_v)))
+        n_max_diff = float(np.max(np.abs(ours_n - ref_n)))
+        print(f"Correctness vs reference (iterations=1): "
+              f"max|Δverts|={v_max_diff:.2e}, max|Δnormals|={n_max_diff:.2e}")
+
+    # Print result table.
+    print()
+    print(f"{'impl':<32} {'median (s)':>12} {'best (s)':>12} {'speedup':>10}")
+    baseline = None
+    for row in rows:
+        if row["impl"] == "reference[nifty]":
+            baseline = row["median_s"]
+            break
+    for row in rows:
+        speedup = baseline / row["median_s"] if baseline else float("nan")
+        print(
+            f"{row['impl']:<32} {row['median_s']:>12.4f} {row['best_s']:>12.4f} "
+            f"{speedup:>10.2f}x"
+        )
+
+    if args.csv:
+        with open(args.csv, "w", newline="") as f:
+            writer = csv.DictWriter(f, fieldnames=["impl", "median_s", "best_s"])
+            writer.writeheader()
+            for row in rows:
+                writer.writerow(row)
+        print(f"\nWrote {args.csv}")
+
+    return 0
+
+
+if __name__ == "__main__":
+    raise SystemExit(main())
diff --git a/include/bioimage_cpp/mesh_smoothing.hxx b/include/bioimage_cpp/mesh_smoothing.hxx
new file mode 100644
index 0000000..9637e5a
--- /dev/null
+++ b/include/bioimage_cpp/mesh_smoothing.hxx
@@ -0,0 +1,187 @@
+#pragma once
+
+#include "bioimage_cpp/array_view.hxx"
+#include "bioimage_cpp/detail/threading.hxx"
+
+#include <algorithm>
+#include <cstddef>
+#include <cstdint>
+#include <numeric>
+#include <stdexcept>
+#include <string>
+#include <utility>
+#include <vector>
+
+namespace bioimage_cpp {
+
+namespace detail::mesh_smoothing {
+
+template <class I>
+struct Adjacency {
+    std::vector<std::size_t> offsets;
+    std::vector<I> neighbours;
+};
+
+template <class I>
+Adjacency<I> build_adjacency(const ConstArrayView<I> &faces, std::ptrdiff_t n_verts) {
+    const std::ptrdiff_t n_faces = faces.shape[0];
+
+    std::vector<std::pair<I, I>> edges;
+    edges.reserve(static_cast<std::size_t>(n_faces) * 3);
+
+    const auto signed_n_verts = static_cast<std::int64_t>(n_verts);
+    auto check_index = [&](const I value) {
+        const auto signed_value = static_cast<std::int64_t>(value);
+        if (signed_value < 0 || signed_value >= signed_n_verts) {
+            throw std::invalid_argument(
+                "faces contains index " + std::to_string(signed_value)
+                + " out of range [0, " + std::to_string(n_verts) + ")"
+            );
+        }
+    };
+
+    for (std::ptrdiff_t f = 0; f < n_faces; ++f) {
+        const I a = faces.data[f * 3 + 0];
+        const I b = faces.data[f * 3 + 1];
+        const I c = faces.data[f * 3 + 2];
+        check_index(a);
+        check_index(b);
+        check_index(c);
+        edges.emplace_back(std::min(a, b), std::max(a, b));
+        edges.emplace_back(std::min(b, c), std::max(b, c));
+        edges.emplace_back(std::min(a, c), std::max(a, c));
+    }
+
+    std::sort(edges.begin(), edges.end());
+    edges.erase(std::unique(edges.begin(), edges.end()), edges.end());
+
+    std::vector<std::size_t> offsets(static_cast<std::size_t>(n_verts) + 1, 0);
+    for (const auto &edge : edges) {
+        ++offsets[static_cast<std::size_t>(edge.first) + 1];
+        ++offsets[static_cast<std::size_t>(edge.second) + 1];
+    }
+    for (std::size_t i = 1; i < offsets.size(); ++i) {
+        offsets[i] += offsets[i - 1];
+    }
+
+    std::vector<I> neighbours(offsets.back());
+    std::vector<std::size_t> insert_pos(offsets.begin(), offsets.end() - 1);
+    for (const auto &edge : edges) {
+        neighbours[insert_pos[static_cast<std::size_t>(edge.first)]++] = edge.second;
+        neighbours[insert_pos[static_cast<std::size_t>(edge.second)]++] = edge.first;
+    }
+
+    return Adjacency<I>{std::move(offsets), std::move(neighbours)};
+}
+
+} // namespace detail::mesh_smoothing
+
+// Laplacian smoothing of a triangular mesh: each vertex (and corresponding
+// normal) is replaced by the mean of itself and its 1-ring neighbours,
+// repeated for `iterations` passes. `verts` and `normals` are (n_verts, dim);
+// `faces` is (n_faces, 3) with values in [0, n_verts). Adjacency is built once
+// and reused across iterations. For `iterations == 0`, inputs are copied to
+// outputs unchanged.
+template <class V, class I>
+void smooth_mesh(
+    const ConstArrayView<V> &verts,
+    const ConstArrayView<V> &normals,
+    const ConstArrayView<I> &faces,
+    std::size_t iterations,
+    int n_threads,
+    const ArrayView<V> &out_verts,
+    const ArrayView<V> &out_normals
+) {
+    if (verts.shape.size() != 2) {
+        throw std::invalid_argument(
+            "verts must have ndim=2, got ndim=" + std::to_string(verts.shape.size())
+        );
+    }
+    if (normals.shape != verts.shape) {
+        throw std::invalid_argument("normals shape must match verts shape");
+    }
+    if (faces.shape.size() != 2 || faces.shape[1] != 3) {
+        throw std::invalid_argument("faces must have shape (n_faces, 3)");
+    }
+    if (out_verts.shape != verts.shape) {
+        throw std::invalid_argument("out_verts shape must match verts shape");
+    }
+    if (out_normals.shape != verts.shape) {
+        throw std::invalid_argument("out_normals shape must match verts shape");
+    }
+
+    const std::ptrdiff_t n_verts = verts.shape[0];
+    const std::ptrdiff_t dim = verts.shape[1];
+    const std::size_t n_total = static_cast<std::size_t>(n_verts) * static_cast<std::size_t>(dim);
+
+    if (iterations == 0) {
+        std::copy(verts.data, verts.data + n_total, out_verts.data);
+        std::copy(normals.data, normals.data + n_total, out_normals.data);
+        return;
+    }
+
+    const auto adjacency = detail::mesh_smoothing::build_adjacency<I>(faces, n_verts);
+
+    std::vector<V> scratch_verts(n_total);
+    std::vector<V> scratch_normals(n_total);
+
+    // Choose initial write target so the final write lands in out_*.
+    // iter 0 reads from `verts`/`normals` (input); subsequent iters alternate
+    // between out_* and scratch. The last write must be to out_*, so when
+    // iterations is even, iter 0 writes scratch (and iter 1 writes out).
+    V *buf_a_verts = out_verts.data;
+    V *buf_a_normals = out_normals.data;
+    V *buf_b_verts = scratch_verts.data();
+    V *buf_b_normals = scratch_normals.data();
+    if (iterations % 2 == 0) {
+        std::swap(buf_a_verts, buf_b_verts);
+        std::swap(buf_a_normals, buf_b_normals);
+    }
+
+    const auto &offsets = adjacency.offsets;
+    const auto &neighbours = adjacency.neighbours;
+    const std::size_t threads = detail::normalize_thread_count(
+        n_threads < 0 ? 0 : static_cast<std::size_t>(n_threads),
+        static_cast<std::size_t>(n_verts)
+    );
+
+    auto smooth_pass = [&](const V *src_verts, const V *src_normals, V *dst_verts, V *dst_normals) {
+        detail::parallel_for_chunks(
+            threads,
+            static_cast<std::size_t>(n_verts),
+            [&](std::size_t, std::size_t begin, std::size_t end) {
+                for (std::size_t v = begin; v < end; ++v) {
+                    const std::size_t beg = offsets[v];
+                    const std::size_t fin = offsets[v + 1];
+                    const auto count = static_cast<V>(fin - beg + 1);
+                    const std::size_t row = v * static_cast<std::size_t>(dim);
+                    for (std::ptrdiff_t d = 0; d < dim; ++d) {
+                        V sum_v = src_verts[row + static_cast<std::size_t>(d)];
+                        V sum_n = src_normals[row + static_cast<std::size_t>(d)];
+                        for (std::size_t k = beg; k < fin; ++k) {
+                            const std::size_t nbr_row =
+                                static_cast<std::size_t>(neighbours[k]) * static_cast<std::size_t>(dim);
+                            sum_v += src_verts[nbr_row + static_cast<std::size_t>(d)];
+                            sum_n += src_normals[nbr_row + static_cast<std::size_t>(d)];
+                        }
+                        dst_verts[row + static_cast<std::size_t>(d)] = sum_v / count;
+                        dst_normals[row + static_cast<std::size_t>(d)] = sum_n / count;
+                    }
+                }
+            }
+        );
+    };
+
+    // First iteration: read from inputs, write to buf_a.
+    smooth_pass(verts.data, normals.data, buf_a_verts, buf_a_normals);
+
+    for (std::size_t it = 1; it < iterations; ++it) {
+        const V *src_verts = (it % 2 == 1) ? buf_a_verts : buf_b_verts;
+        const V *src_normals = (it % 2 == 1) ? buf_a_normals : buf_b_normals;
+        V *dst_verts = (it % 2 == 1) ? buf_b_verts : buf_a_verts;
+        V *dst_normals = (it % 2 == 1) ? buf_b_normals : buf_a_normals;
+        smooth_pass(src_verts, src_normals, dst_verts, dst_normals);
+    }
+}
+
+} // namespace bioimage_cpp
diff --git a/src/bindings/utils.cxx b/src/bindings/utils.cxx
index b51de32..c279c07 100644
--- a/src/bindings/utils.cxx
+++ b/src/bindings/utils.cxx
@@ -1,9 +1,11 @@
 #include "utils.hxx"
 
 #include "bioimage_cpp/array_view.hxx"
+#include "bioimage_cpp/mesh_smoothing.hxx"
 #include "bioimage_cpp/take_dict.hxx"
 
 #include <nanobind/ndarray.h>
+#include <nanobind/stl/pair.h>
 #include <nanobind/stl/string.h>
 
 #include <cstddef>
@@ -13,6 +15,7 @@
 #include <stdexcept>
 #include <string>
 #include <unordered_map>
+#include <utility>
 #include <vector>
 
 namespace nb = nanobind;
@@ -77,6 +80,77 @@ OutputArray<T> take_dict_t(const nb::dict &relabeling, InputArray<T> to_relabel)
     return OutputArray<T>(data, ndarray_shape.size(), ndarray_shape.data(), owner);
 }
 
+using FacesArray = nb::ndarray<nb::numpy, const std::int64_t, nb::c_contig>;
+
+template <class V>
+std::pair<OutputArray<V>, OutputArray<V>> smooth_mesh_t(
+    InputArray<V> verts,
+    InputArray<V> normals,
+    FacesArray faces,
+    std::size_t iterations,
+    int n_threads
+) {
+    if (verts.ndim() != 2) {
+        throw std::invalid_argument(
+            "verts must have ndim=2, got ndim=" + std::to_string(verts.ndim())
+        );
+    }
+    if (normals.ndim() != 2 || normals.shape(0) != verts.shape(0)
+        || normals.shape(1) != verts.shape(1)) {
+        throw std::invalid_argument("normals must have the same shape as verts");
+    }
+    if (faces.ndim() != 2 || faces.shape(1) != 3) {
+        throw std::invalid_argument("faces must have shape (n_faces, 3)");
+    }
+
+    const std::size_t n_verts_size = verts.shape(0);
+    const std::size_t dim_size = verts.shape(1);
+    const std::size_t n_faces_size = faces.shape(0);
+    const std::size_t n_total = n_verts_size * dim_size;
+
+    std::vector<std::size_t> verts_ndarray_shape{n_verts_size, dim_size};
+    std::vector<std::ptrdiff_t> verts_view_shape{
+        static_cast<std::ptrdiff_t>(n_verts_size),
+        static_cast<std::ptrdiff_t>(dim_size),
+    };
+    std::vector<std::ptrdiff_t> faces_view_shape{
+        static_cast<std::ptrdiff_t>(n_faces_size),
+        std::ptrdiff_t{3},
+    };
+
+    auto *verts_data = new V[n_total]();
+    nb::capsule verts_owner(verts_data, [](void *p) noexcept { delete[] static_cast<V *>(p); });
+    auto *normals_data = new V[n_total]();
+    nb::capsule normals_owner(normals_data, [](void *p) noexcept { delete[] static_cast<V *>(p); });
+
+    ConstArrayView<V> verts_view{verts.data(), verts_view_shape, {}};
+    ConstArrayView<V> normals_view{normals.data(), verts_view_shape, {}};
+    ConstArrayView<std::int64_t> faces_view{faces.data(), faces_view_shape, {}};
+    ArrayView<V> out_verts_view{verts_data, verts_view_shape, {}};
+    ArrayView<V> out_normals_view{normals_data, verts_view_shape, {}};
+
+    {
+        nb::gil_scoped_release release;
+        smooth_mesh<V, std::int64_t>(
+            verts_view,
+            normals_view,
+            faces_view,
+            iterations,
+            n_threads,
+            out_verts_view,
+            out_normals_view
+        );
+    }
+
+    OutputArray<V> out_verts(
+        verts_data, verts_ndarray_shape.size(), verts_ndarray_shape.data(), verts_owner
+    );
+    OutputArray<V> out_normals(
+        normals_data, verts_ndarray_shape.size(), verts_ndarray_shape.data(), normals_owner
+    );
+    return {std::move(out_verts), std::move(out_normals)};
+}
+
 } // namespace
 
 void bind_utils(nb::module_ &m) {
@@ -108,6 +182,26 @@ void bind_utils(nb::module_ &m) {
         nb::arg("to_relabel"),
         "Map a contiguous int64 array through an integer dictionary."
     );
+    m.def(
+        "_smooth_mesh_float32",
+        &smooth_mesh_t<float>,
+        nb::arg("verts"),
+        nb::arg("normals"),
+        nb::arg("faces"),
+        nb::arg("iterations"),
+        nb::arg("n_threads"),
+        "Laplacian smoothing of a triangular mesh with float32 vertices and normals."
+    );
+    m.def(
+        "_smooth_mesh_float64",
+        &smooth_mesh_t<double>,
+        nb::arg("verts"),
+        nb::arg("normals"),
+        nb::arg("faces"),
+        nb::arg("iterations"),
+        nb::arg("n_threads"),
+        "Laplacian smoothing of a triangular mesh with float64 vertices and normals."
+    );
 }
 
 } // namespace bioimage_cpp::bindings
diff --git a/src/bioimage_cpp/utils.py b/src/bioimage_cpp/utils.py
index cd61d02..cb43e57 100644
--- a/src/bioimage_cpp/utils.py
+++ b/src/bioimage_cpp/utils.py
@@ -19,6 +19,11 @@
     np.dtype("int64"): _core._take_dict_int64,
 }
 
+_SMOOTH_MESH_BY_DTYPE = {
+    np.dtype("float32"): _core._smooth_mesh_float32,
+    np.dtype("float64"): _core._smooth_mesh_float64,
+}
+
 _COUNT_TABLE_DTYPE = np.dtype([("label", np.uint64), ("count", np.uint64)])
 _OVERLAP_TABLE_DTYPE = np.dtype(
     [("label_a", np.uint64), ("label_b", np.uint64), ("count", np.uint64)]
@@ -302,6 +307,111 @@ def take_dict(relabeling: Mapping[Any, Any], to_relabel: np.ndarray) -> np.ndarr
         raise
 
 
+def smooth_mesh(
+    verts: np.ndarray,
+    normals: np.ndarray,
+    faces: np.ndarray,
+    iterations: int,
+    *,
+    n_threads: int = 0,
+) -> tuple[np.ndarray, np.ndarray]:
+    """Laplacian smoothing of a triangular mesh.
+
+    Each vertex (and the corresponding normal) is iteratively replaced by the
+    mean of itself and its 1-ring neighbours in the mesh. Updates are applied
+    in Jacobi (double-buffered) style: every vertex in an iteration is updated
+    from the same input state, independent of vertex order. The adjacency is
+    built once from ``faces``; each interior edge is treated as a single
+    undirected edge regardless of how many faces share it.
+
+    Parameters
+    ----------
+    verts:
+        2D ``float32`` or ``float64`` array of shape ``(n_verts, dim)``.
+    normals:
+        Array with the same shape and dtype as ``verts``.
+    faces:
+        2D integer array of shape ``(n_faces, 3)`` with values in
+        ``[0, n_verts)``. Non-``int64`` inputs are converted internally.
+    iterations:
+        Number of smoothing passes. ``0`` returns copies of the inputs.
+    n_threads:
+        Number of threads for the inner smoothing loop. ``0`` (the default)
+        picks the hardware concurrency; ``1`` forces serial execution.
+
+    Returns
+    -------
+    (np.ndarray, np.ndarray)
+        Smoothed vertices and normals with the same shape and dtype as the
+        inputs.
+    """
+    verts_array = np.asarray(verts)
+    normals_array = np.asarray(normals)
+    faces_array = np.asarray(faces)
+
+    if verts_array.ndim != 2:
+        raise ValueError(
+            f"verts must have ndim=2, got ndim={verts_array.ndim}"
+        )
+    if normals_array.shape != verts_array.shape:
+        raise ValueError(
+            "normals must have the same shape as verts, got "
+            f"verts shape={verts_array.shape}, normals shape={normals_array.shape}"
+        )
+    if verts_array.dtype != normals_array.dtype:
+        raise TypeError(
+            "verts and normals must have the same dtype, got "
+            f"verts dtype={verts_array.dtype}, normals dtype={normals_array.dtype}"
+        )
+
+    try:
+        smooth = _SMOOTH_MESH_BY_DTYPE[verts_array.dtype]
+    except KeyError as error:
+        supported = ", ".join(str(dtype) for dtype in _SMOOTH_MESH_BY_DTYPE)
+        raise TypeError(
+            f"verts must have one of dtypes ({supported}), got dtype={verts_array.dtype}"
+        ) from error
+
+    if faces_array.ndim != 2 or faces_array.shape[1] != 3:
+        raise ValueError(
+            f"faces must have shape (n_faces, 3), got shape={faces_array.shape}"
+        )
+    if not np.issubdtype(faces_array.dtype, np.integer):
+        raise TypeError(
+            f"faces must have an integer dtype, got dtype={faces_array.dtype}"
+        )
+
+    n_verts = verts_array.shape[0]
+    if faces_array.size > 0:
+        face_min = int(faces_array.min())
+        face_max = int(faces_array.max())
+        if face_min < 0 or face_max >= n_verts:
+            raise ValueError(
+                "faces must contain indices in [0, n_verts), got values in "
+                f"[{face_min}, {face_max}] with n_verts={n_verts}"
+            )
+
+    iterations_int = int(iterations)
+    if iterations_int < 0:
+        raise ValueError(f"iterations must be non-negative, got {iterations_int}")
+
+    n_threads_int = int(n_threads)
+    if n_threads_int < 0:
+        raise ValueError(f"n_threads must be non-negative, got {n_threads_int}")
+
+    verts_contiguous = np.ascontiguousarray(verts_array)
+    normals_contiguous = np.ascontiguousarray(normals_array)
+    faces_contiguous = np.ascontiguousarray(faces_array, dtype=np.int64)
+
+    return smooth(
+        verts_contiguous,
+        normals_contiguous,
+        faces_contiguous,
+        iterations_int,
+        n_threads_int,
+    )
+
+
 def _normalize_labels(labels: np.ndarray, name: str) -> np.ndarray:
     array = np.asarray(labels)
     if not np.issubdtype(array.dtype, np.integer):
@@ -364,5 +474,6 @@ def _validate_normalize_by(normalize_by: str) -> None:
     "SegmentationOverlap",
     "UnionFind",
     "segmentation_overlap",
+    "smooth_mesh",
     "take_dict",
 ]
diff --git a/tests/test_utils_mesh_smoothing.py b/tests/test_utils_mesh_smoothing.py
new file mode 100644
index 0000000..7f3b46e
--- /dev/null
+++ b/tests/test_utils_mesh_smoothing.py
@@ -0,0 +1,247 @@
+import os
+import sys
+
+import numpy as np
+import pytest
+
+import bioimage_cpp as bic
+
+
+def _tetrahedron(dtype):
+    verts = np.array(
+        [[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]],
+        dtype=dtype,
+    )
+    normals = verts.copy() + dtype(0.5)
+    faces = np.array(
+        [[0, 1, 2], [0, 1, 3], [0, 2, 3], [1, 2, 3]], dtype=np.int64
+    )
+    return verts, normals, faces
+
+
+def _octahedron():
+    verts = np.array(
+        [
+            [1.0, 0.0, 0.0],
+            [-1.0, 0.0, 0.0],
+            [0.0, 1.0, 0.0],
+            [0.0, -1.0, 0.0],
+            [0.0, 0.0, 1.0],
+            [0.0, 0.0, -1.0],
+        ],
+        dtype=np.float64,
+    )
+    faces = np.array(
+        [
+            [0, 2, 4],
+            [2, 1, 4],
+            [1, 3, 4],
+            [3, 0, 4],
+            [0, 5, 2],
+            [2, 5, 1],
+            [1, 5, 3],
+            [3, 5, 0],
+        ],
+        dtype=np.int64,
+    )
+    return verts, faces
+
+
+@pytest.mark.parametrize("dtype", [np.float32, np.float64])
+def test_tetrahedron_one_iteration_collapses_to_centroid(dtype):
+    verts, normals, faces = _tetrahedron(dtype)
+    out_verts, out_normals = bic.utils.smooth_mesh(verts, normals, faces, iterations=1)
+
+    assert out_verts.dtype == verts.dtype
+    assert out_normals.dtype == normals.dtype
+    assert out_verts.shape == verts.shape
+    assert out_normals.shape == normals.shape
+
+    centroid_verts = verts.mean(axis=0, dtype=dtype)
+    centroid_normals = normals.mean(axis=0, dtype=dtype)
+    rtol = 1e-5 if dtype == np.float32 else 1e-12
+    np.testing.assert_allclose(out_verts, np.broadcast_to(centroid_verts, verts.shape), rtol=rtol)
+    np.testing.assert_allclose(out_normals, np.broadcast_to(centroid_normals, normals.shape), rtol=rtol)
+
+
+def test_iterations_zero_returns_copy_of_inputs():
+    verts, normals, faces = _tetrahedron(np.float64)
+    out_verts, out_normals = bic.utils.smooth_mesh(verts, normals, faces, iterations=0)
+
+    np.testing.assert_array_equal(out_verts, verts)
+    np.testing.assert_array_equal(out_normals, normals)
+    # Outputs must be independent copies (mutation does not affect inputs).
+    out_verts[0, 0] = 42.0
+    assert verts[0, 0] != 42.0
+
+
+def test_octahedron_converges_to_centroid():
+    verts, faces = _octahedron()
+    normals = verts.copy()
+    out_verts, _ = bic.utils.smooth_mesh(verts, normals, faces, iterations=100)
+    centroid = verts.mean(axis=0)
+    # Every vertex should be near the centroid after enough iterations.
+    np.testing.assert_allclose(out_verts, np.broadcast_to(centroid, verts.shape), atol=1e-6)
+
+
+@pytest.mark.parametrize("dim", [2, 4])
+def test_non_3d_feature_vectors(dim):
+    rng = np.random.default_rng(42)
+    n_verts = 8
+    verts = rng.standard_normal((n_verts, dim)).astype(np.float64)
+    normals = rng.standard_normal((n_verts, dim)).astype(np.float64)
+    faces = np.array(
+        [[0, 1, 2], [1, 2, 3], [4, 5, 6], [5, 6, 7], [0, 4, 7], [2, 3, 5]],
+        dtype=np.int64,
+    )
+    out_verts, out_normals = bic.utils.smooth_mesh(verts, normals, faces, iterations=3)
+    assert out_verts.shape == (n_verts, dim)
+    assert out_normals.shape == (n_verts, dim)
+
+
+def test_parity_with_iterations_alternation():
+    # Run iterations one-at-a-time and compare against a single multi-iteration call.
+    verts, normals, faces = _tetrahedron(np.float64)
+    cur_v, cur_n = verts, normals
+    for _ in range(5):
+        cur_v, cur_n = bic.utils.smooth_mesh(cur_v, cur_n, faces, iterations=1)
+    multi_v, multi_n = bic.utils.smooth_mesh(verts, normals, faces, iterations=5)
+    np.testing.assert_allclose(multi_v, cur_v, rtol=1e-12)
+    np.testing.assert_allclose(multi_n, cur_n, rtol=1e-12)
+
+
+def test_threading_matches_serial():
+    rng = np.random.default_rng(7)
+    n_verts = 256
+    verts = rng.standard_normal((n_verts, 3)).astype(np.float64)
+    normals = rng.standard_normal((n_verts, 3)).astype(np.float64)
+    # Triangle strip across the vertices.
+    faces = np.stack(
+        [np.arange(0, n_verts - 2), np.arange(1, n_verts - 1), np.arange(2, n_verts)],
+        axis=1,
+    ).astype(np.int64)
+
+    serial = bic.utils.smooth_mesh(verts, normals, faces, iterations=8, n_threads=1)
+    parallel = bic.utils.smooth_mesh(verts, normals, faces, iterations=8, n_threads=4)
+    auto = bic.utils.smooth_mesh(verts, normals, faces, iterations=8, n_threads=0)
+    np.testing.assert_array_equal(serial[0], parallel[0])
+    np.testing.assert_array_equal(serial[1], parallel[1])
+    np.testing.assert_array_equal(serial[0], auto[0])
+    np.testing.assert_array_equal(serial[1], auto[1])
+
+
+def test_non_contiguous_inputs_are_accepted():
+    verts, normals, faces = _tetrahedron(np.float64)
+    # Make non-contiguous views by replicating and slicing.
+    big_verts = np.repeat(verts, 2, axis=0)
+    big_normals = np.repeat(normals, 2, axis=0)
+    verts_view = big_verts[::2]
+    normals_view = big_normals[::2]
+    assert not verts_view.flags.c_contiguous
+
+    out_verts, out_normals = bic.utils.smooth_mesh(verts_view, normals_view, faces, iterations=1)
+    expected_v, expected_n = bic.utils.smooth_mesh(verts, normals, faces, iterations=1)
+    np.testing.assert_array_equal(out_verts, expected_v)
+    np.testing.assert_array_equal(out_normals, expected_n)
+
+
+def test_faces_dtype_is_normalised_to_int64():
+    verts, normals, faces = _tetrahedron(np.float64)
+    out_int64, _ = bic.utils.smooth_mesh(verts, normals, faces.astype(np.int64), iterations=1)
+    out_int32, _ = bic.utils.smooth_mesh(verts, normals, faces.astype(np.int32), iterations=1)
+    out_uint32, _ = bic.utils.smooth_mesh(verts, normals, faces.astype(np.uint32), iterations=1)
+    np.testing.assert_array_equal(out_int64, out_int32)
+    np.testing.assert_array_equal(out_int64, out_uint32)
+
+
+def test_rejects_mismatched_shapes():
+    verts = np.zeros((4, 3), dtype=np.float64)
+    normals_bad = np.zeros((4, 2), dtype=np.float64)
+    faces = np.array([[0, 1, 2]], dtype=np.int64)
+    with pytest.raises(ValueError, match="same shape as verts"):
+        bic.utils.smooth_mesh(verts, normals_bad, faces, iterations=1)
+
+
+def test_rejects_1d_verts():
+    verts = np.zeros(4, dtype=np.float64)
+    normals = verts.copy()
+    faces = np.array([[0, 1, 2]], dtype=np.int64)
+    with pytest.raises(ValueError, match="verts must have ndim=2"):
+        bic.utils.smooth_mesh(verts, normals, faces, iterations=1)
+
+
+def test_rejects_bad_faces_shape():
+    verts, normals, _ = _tetrahedron(np.float64)
+    bad_faces = np.array([[0, 1, 2, 3]], dtype=np.int64)
+    with pytest.raises(ValueError, match=r"faces must have shape"):
+        bic.utils.smooth_mesh(verts, normals, bad_faces, iterations=1)
+
+
+def test_rejects_out_of_range_faces():
+    verts, normals, _ = _tetrahedron(np.float64)
+    bad_faces = np.array([[0, 1, 10]], dtype=np.int64)
+    with pytest.raises(ValueError, match=r"\[0, n_verts\)"):
+        bic.utils.smooth_mesh(verts, normals, bad_faces, iterations=1)
+
+
+def test_rejects_negative_iterations():
+    verts, normals, faces = _tetrahedron(np.float64)
+    with pytest.raises(ValueError, match="iterations must be non-negative"):
+        bic.utils.smooth_mesh(verts, normals, faces, iterations=-1)
+
+
+def test_rejects_negative_n_threads():
+    verts, normals, faces = _tetrahedron(np.float64)
+    with pytest.raises(ValueError, match="n_threads must be non-negative"):
+        bic.utils.smooth_mesh(verts, normals, faces, iterations=1, n_threads=-1)
+
+
+def test_rejects_unsupported_verts_dtype():
+    verts = np.zeros((4, 3), dtype=np.int32)
+    normals = verts.copy()
+    faces = np.array([[0, 1, 2]], dtype=np.int64)
+    with pytest.raises(TypeError, match="verts must have one of dtypes"):
+        bic.utils.smooth_mesh(verts, normals, faces, iterations=1)
+
+
+def test_rejects_mismatched_verts_normals_dtype():
+    verts = np.zeros((4, 3), dtype=np.float32)
+    normals = np.zeros((4, 3), dtype=np.float64)
+    faces = np.array([[0, 1, 2]], dtype=np.int64)
+    with pytest.raises(TypeError, match="same dtype"):
+        bic.utils.smooth_mesh(verts, normals, faces, iterations=1)
+
+
+def test_parity_with_python_reference():
+    """Compare against the nifty-based Python reference. Skipped if nifty missing.
+
+    Only one iteration is compared. The reference has an aliasing quirk
+    (``current_verts = new_verts`` makes them refer to the same buffer) that
+    turns iterations 1+ into in-place Gauss-Seidel smoothing, whereas this
+    implementation does textbook Jacobi smoothing (independent read/write
+    buffers). The two agree exactly at ``iterations=1``.
+    """
+    pytest.importorskip("nifty")
+
+    repo_root = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))
+    dev_path = os.path.join(repo_root, "development", "utils")
+    sys.path.insert(0, dev_path)
+    try:
+        from _mesh_smoothing_reference import smooth_mesh as smooth_mesh_reference
+    finally:
+        sys.path.remove(dev_path)
+
+    rng = np.random.default_rng(2026)
+    scipy_spatial = pytest.importorskip("scipy.spatial")
+    n_points = 80
+    raw = rng.standard_normal((n_points, 3))
+    points = raw / np.linalg.norm(raw, axis=1, keepdims=True)
+    hull = scipy_spatial.ConvexHull(points)
+    verts = points.astype(np.float64)
+    faces = hull.simplices.astype(np.int64)
+    normals = verts.copy()
+
+    out_v, out_n = bic.utils.smooth_mesh(verts, normals, faces, iterations=1)
+    ref_v, ref_n = smooth_mesh_reference(verts, normals, faces, 1)
+    np.testing.assert_allclose(out_v, ref_v, rtol=1e-10, atol=1e-12)
+    np.testing.assert_allclose(out_n, ref_n, rtol=1e-10, atol=1e-12)