diff --git a/.gitignore b/.gitignore
index 983b472..1ce03dc 100644
--- a/.gitignore
+++ b/.gitignore
@@ -24,3 +24,6 @@ htmlcov/
 # Local virtual environments
 .venv/
 venv/
+
+# Data
+*.h5
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 388e439..b3adef2 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -12,7 +12,9 @@ find_package(nanobind CONFIG REQUIRED)
 nanobind_add_module(_core
     NB_STATIC
     src/bindings/module.cxx
+    src/bindings/segmentation.cxx
     src/bindings/utils.cxx
+    src/cpp/segmentation/mutex_watershed.cxx
     src/cpp/take_dict.cxx
 )
 
diff --git a/README.md b/README.md
index 23c7e9b..af49860 100644
--- a/README.md
+++ b/README.md
@@ -30,6 +30,17 @@ out = bic.utils.take_dict(relabeling, labels)
 # array([10, 30, 20, 10], dtype=uint64)
 ```
 
+```python
+affinities = np.ones((2, 4, 4), dtype=np.float32)
+offsets = [[0, 1], [1, 0]]
+
+segmentation = bic.segmentation.mutex_watershed(
+    affinities,
+    offsets,
+    number_of_attractive_channels=2,
+)
+```
+
 ## Scope
 
 The project is not a compatibility layer for `nifty`, `vigra`, or other large libraries. It keeps I/O and heavy dependencies out of the C++ core; callers should use existing Python packages for file formats and pass NumPy arrays into `bioimage-cpp`.
diff --git a/development/segmentation/_mutex_watershed_equivalence.py b/development/segmentation/_mutex_watershed_equivalence.py
new file mode 100644
index 0000000..ea51b3f
--- /dev/null
+++ b/development/segmentation/_mutex_watershed_equivalence.py
@@ -0,0 +1,227 @@
+from __future__ import annotations
+
+import argparse
+from pathlib import Path
+from statistics import median
+from time import perf_counter
+from typing import Callable
+
+import numpy as np
+
+
+PROJECT_ROOT = Path(__file__).resolve().parents[2]
+DEFAULT_DATA_PREFIX = PROJECT_ROOT / "examples" / "segmentation" / "isbi-data-"
+
+
+def load_problem(data_prefix: Path | str = DEFAULT_DATA_PREFIX):
+    from elf.segmentation.utils import load_mutex_watershed_problem
+
+    affinities, offsets = load_mutex_watershed_problem(prefix=str(data_prefix))
+    return np.ascontiguousarray(affinities), [tuple(offset) for offset in offsets]
+
+
+def prepare_2d_problem(
+    affinities: np.ndarray,
+    offsets: list[tuple[int, ...]],
+    z: int,
+    yx_shape: tuple[int, int],
+):
+    channels_2d = [i for i, offset in enumerate(offsets) if offset[0] == 0]
+    y, x = yx_shape
+    cropped = affinities[channels_2d, z, :y, :x]
+    offsets_2d = [offsets[i][1:] for i in channels_2d]
+    return np.ascontiguousarray(cropped), offsets_2d, 2
+
+
+def prepare_3d_problem(
+    affinities: np.ndarray,
+    offsets: list[tuple[int, ...]],
+    zyx_shape: tuple[int, int, int],
+):
+    z, y, x = zyx_shape
+    cropped = affinities[:, :z, :y, :x]
+    return np.ascontiguousarray(cropped), offsets, 3
+
+
+def run_bioimage_cpp(
+    affinities: np.ndarray,
+    offsets: list[tuple[int, ...]],
+    number_of_attractive_channels: int,
+) -> np.ndarray:
+    import bioimage_cpp as bic
+
+    affs = affinities.copy()
+    affs[:number_of_attractive_channels] *= -1
+    affs[:number_of_attractive_channels] += 1
+    return bic.segmentation.mutex_watershed(
+        affs,
+        offsets,
+        number_of_attractive_channels=number_of_attractive_channels,
+    )
+
+
+def run_affogato_reference(
+    affinities: np.ndarray,
+    offsets: list[tuple[int, ...]],
+    number_of_attractive_channels: int,
+) -> np.ndarray:
+    from elf.segmentation.mutex_watershed import mutex_watershed
+
+    spatial_ndim = len(offsets[0])
+    if number_of_attractive_channels != spatial_ndim:
+        raise ValueError(
+            "the elf mutex_watershed wrapper assumes one attractive channel per "
+            f"spatial axis, got ndim={spatial_ndim}, attractive channels="
+            f"{number_of_attractive_channels}"
+        )
+    return mutex_watershed(affinities.copy(), offsets, strides=[1] * spatial_ndim)
+
+
+def _load_validation_metrics():
+    try:
+        from elf.validation import rand_index, variation_of_information
+
+        return "elf.validation", rand_index, variation_of_information
+    except ImportError:
+        from elf.evaluation import rand_index, variation_of_information
+
+        return "elf.evaluation", rand_index, variation_of_information
+
+
+def compare_segmentations(
+    candidate: np.ndarray,
+    reference: np.ndarray,
+    *,
+    max_vi: float = 1.0e-10,
+    max_are: float = 1.0e-10,
+    min_rand_index: float = 1.0 - 1.0e-10,
+) -> dict[str, float | str | bool]:
+    source, rand_index, variation_of_information = _load_validation_metrics()
+
+    vi_split, vi_merge = variation_of_information(candidate, reference)
+    adapted_rand_error, ri = rand_index(candidate, reference)
+    exact_equal = bool(np.array_equal(candidate, reference))
+    equivalent = (
+        vi_split <= max_vi
+        and vi_merge <= max_vi
+        and adapted_rand_error <= max_are
+        and ri >= min_rand_index
+    )
+    metrics: dict[str, float | str | bool] = {
+        "validation_source": source,
+        "vi_split": float(vi_split),
+        "vi_merge": float(vi_merge),
+        "adapted_rand_error": float(adapted_rand_error),
+        "rand_index": float(ri),
+        "exact_label_equality": exact_equal,
+        "equivalent": equivalent,
+    }
+    if not equivalent:
+        raise AssertionError(
+            "mutex watershed results differ: "
+            f"VI split={vi_split:.6g}, VI merge={vi_merge:.6g}, "
+            f"adapted rand error={adapted_rand_error:.6g}, "
+            f"rand index={ri:.12g}, exact labels={exact_equal}"
+        )
+    return metrics
+
+
+def time_function(
+    run: Callable[[np.ndarray, list[tuple[int, ...]], int], np.ndarray],
+    affinities: np.ndarray,
+    offsets: list[tuple[int, ...]],
+    number_of_attractive_channels: int,
+    repeats: int,
+) -> tuple[list[float], np.ndarray]:
+    timings = []
+    result = None
+    for _ in range(repeats):
+        start = perf_counter()
+        result = run(affinities, offsets, number_of_attractive_channels)
+        timings.append(perf_counter() - start)
+    assert result is not None
+    return timings, result
+
+
+def print_report(
+    *,
+    ndim: int,
+    affinities: np.ndarray,
+    metrics: dict[str, float | str | bool],
+    bic_timings: list[float],
+    ref_timings: list[float],
+):
+    bic_median = median(bic_timings)
+    ref_median = median(ref_timings)
+    speedup = ref_median / bic_median if bic_median > 0 else float("inf")
+
+    print(f"Mutex watershed {ndim}D equivalence check")
+    print(f"affinities shape: {affinities.shape}, dtype: {affinities.dtype}")
+    print(f"validation metrics: {metrics['validation_source']}")
+    print(
+        "VI split/merge: "
+        f"{metrics['vi_split']:.6g} / {metrics['vi_merge']:.6g}"
+    )
+    print(
+        "adapted rand error / rand index: "
+        f"{metrics['adapted_rand_error']:.6g} / {metrics['rand_index']:.12g}"
+    )
+    print(f"exact label equality: {metrics['exact_label_equality']}")
+    print(f"bioimage-cpp median runtime: {bic_median:.6f} s")
+    print(f"affogato reference median runtime: {ref_median:.6f} s")
+    print(f"reference / bioimage-cpp runtime ratio: {speedup:.3f}x")
+
+
+def run_check(
+    *,
+    ndim: int,
+    repeats: int,
+    data_prefix: Path | str,
+    z: int,
+    yx_shape: tuple[int, int],
+    zyx_shape: tuple[int, int, int],
+):
+    affinities, offsets = load_problem(data_prefix)
+    if ndim == 2:
+        affs, used_offsets, attractive_channels = prepare_2d_problem(
+            affinities, offsets, z=z, yx_shape=yx_shape
+        )
+    elif ndim == 3:
+        affs, used_offsets, attractive_channels = prepare_3d_problem(
+            affinities, offsets, zyx_shape=zyx_shape
+        )
+    else:
+        raise ValueError(f"ndim must be 2 or 3, got {ndim}")
+
+    ref_timings, ref_seg = time_function(
+        run_affogato_reference, affs, used_offsets, attractive_channels, repeats
+    )
+    bic_timings, bic_seg = time_function(
+        run_bioimage_cpp, affs, used_offsets, attractive_channels, repeats
+    )
+    metrics = compare_segmentations(bic_seg, ref_seg)
+    print_report(
+        ndim=ndim,
+        affinities=affs,
+        metrics=metrics,
+        bic_timings=bic_timings,
+        ref_timings=ref_timings,
+    )
+
+
+def add_common_arguments(parser: argparse.ArgumentParser) -> None:
+    parser.add_argument(
+        "--data-prefix",
+        type=Path,
+        default=DEFAULT_DATA_PREFIX,
+        help=(
+            "Path prefix for the ISBI mutex watershed data. The loader expects "
+            "'test.h5' and 'train.h5' suffixes."
+        ),
+    )
+    parser.add_argument(
+        "--repeats",
+        type=int,
+        default=3,
+        help="Number of timed runs for each implementation.",
+    )
diff --git a/development/segmentation/check_mutex_watershed_2d.py b/development/segmentation/check_mutex_watershed_2d.py
new file mode 100644
index 0000000..8960639
--- /dev/null
+++ b/development/segmentation/check_mutex_watershed_2d.py
@@ -0,0 +1,40 @@
+from __future__ import annotations
+
+import argparse
+
+from _mutex_watershed_equivalence import add_common_arguments, run_check
+
+
+def main() -> None:
+    parser = argparse.ArgumentParser(
+        description="Compare bioimage-cpp and affogato mutex watershed in 2D."
+    )
+    add_common_arguments(parser)
+    parser.add_argument(
+        "--z",
+        type=int,
+        default=0,
+        help="Z slice from the example affinities used for the 2D check.",
+    )
+    parser.add_argument(
+        "--shape",
+        type=int,
+        nargs=2,
+        metavar=("Y", "X"),
+        default=(128, 128),
+        help="Spatial crop shape for the 2D check.",
+    )
+    args = parser.parse_args()
+
+    run_check(
+        ndim=2,
+        repeats=args.repeats,
+        data_prefix=args.data_prefix,
+        z=args.z,
+        yx_shape=tuple(args.shape),
+        zyx_shape=(0, 0, 0),
+    )
+
+
+if __name__ == "__main__":
+    main()
diff --git a/development/segmentation/check_mutex_watershed_3d.py b/development/segmentation/check_mutex_watershed_3d.py
new file mode 100644
index 0000000..45ae843
--- /dev/null
+++ b/development/segmentation/check_mutex_watershed_3d.py
@@ -0,0 +1,34 @@
+from __future__ import annotations
+
+import argparse
+
+from _mutex_watershed_equivalence import add_common_arguments, run_check
+
+
+def main() -> None:
+    parser = argparse.ArgumentParser(
+        description="Compare bioimage-cpp and affogato mutex watershed in 3D."
+    )
+    add_common_arguments(parser)
+    parser.add_argument(
+        "--shape",
+        type=int,
+        nargs=3,
+        metavar=("Z", "Y", "X"),
+        default=(6, 96, 96),
+        help="Spatial crop shape for the 3D check.",
+    )
+    args = parser.parse_args()
+
+    run_check(
+        ndim=3,
+        repeats=args.repeats,
+        data_prefix=args.data_prefix,
+        z=0,
+        yx_shape=(0, 0),
+        zyx_shape=tuple(args.shape),
+    )
+
+
+if __name__ == "__main__":
+    main()
diff --git a/examples/segmentation/mutex_watershed.py b/examples/segmentation/mutex_watershed.py
new file mode 100644
index 0000000..7477e48
--- /dev/null
+++ b/examples/segmentation/mutex_watershed.py
@@ -0,0 +1,63 @@
+import napari
+
+from elf.segmentation.utils import load_mutex_watershed_problem
+from elf.io import open_file
+
+
+def mws_bic(affinities, offsets):
+    import bioimage_cpp as bic
+    print("Start MWS ...")
+    affs = affinities.copy()
+    ndim = len(offsets[0])
+    affs[:ndim] *= -1
+    affs[:ndim] += 1
+    segmentation = bic.segmentation.mutex_watershed(
+        affs, offsets, number_of_attractive_channels=ndim
+    )
+    print("done ...")
+    return segmentation
+
+
+def mws_affogato(affinities, offsets):
+    from elf.segmentation.mutex_watershed import mutex_watershed as mws
+    print("Start MWS ...")
+    affs = affinities.copy()
+    segmentation = mws(affs, offsets, strides=[1, 1, 1])
+    print("done ...")
+    return segmentation
+
+
+def _filter_2d(affinities, offsets):
+    chans_2d = [i for i, off in enumerate(offsets) if off[0] == 0]
+    affinities = affinities[chans_2d][:, 0]
+    offsets = [off[1:] for i, off in enumerate(offsets) if i in chans_2d]
+    return affinities, offsets
+
+
+def main():
+    prefix = "isbi-data-"
+    data_path = f"{prefix}test.h5"
+    affinities, offsets = load_mutex_watershed_problem(prefix=prefix)
+
+    check_2d = True
+    if check_2d:
+        affinities, offsets = _filter_2d(affinities, offsets)
+
+    use_reference_impl = False
+    if use_reference_impl:
+        segmentation = mws_affogato(affinities, offsets)
+    else:
+        segmentation = mws_bic(affinities, offsets)
+
+    with open_file(data_path, "r") as f:
+        raw = f["raw"][0] if check_2d else f["raw"][:]
+
+    viewer = napari.Viewer()
+    viewer.add_image(raw, name="raw")
+    viewer.add_image(affinities, name="affinities")
+    viewer.add_labels(segmentation, name="mws-segmentation")
+    napari.run()
+
+
+if __name__ == "__main__":
+    main()
diff --git a/include/bioimage_cpp/segmentation/mutex_watershed.hxx b/include/bioimage_cpp/segmentation/mutex_watershed.hxx
new file mode 100644
index 0000000..b96eb71
--- /dev/null
+++ b/include/bioimage_cpp/segmentation/mutex_watershed.hxx
@@ -0,0 +1,236 @@
+#pragma once
+
+#include "bioimage_cpp/array_view.hxx"
+
+#include <algorithm>
+#include <cstddef>
+#include <cstdint>
+#include <numeric>
+#include <stdexcept>
+#include <string>
+#include <unordered_map>
+#include <unordered_set>
+#include <vector>
+
+namespace bioimage_cpp {
+
+class DisjointSets {
+public:
+    explicit DisjointSets(const std::size_t size) : parents_(size), ranks_(size, 0) {
+        std::iota(parents_.begin(), parents_.end(), std::uint64_t{0});
+    }
+
+    std::uint64_t find(const std::uint64_t node) {
+        if (parents_[node] != node) {
+            parents_[node] = find(parents_[node]);
+        }
+        return parents_[node];
+    }
+
+    std::uint64_t unite_roots(std::uint64_t first, std::uint64_t second) {
+        if (ranks_[first] < ranks_[second]) {
+            std::swap(first, second);
+        }
+        parents_[second] = first;
+        if (ranks_[first] == ranks_[second]) {
+            ++ranks_[first];
+        }
+        return first;
+    }
+
+private:
+    std::vector<std::uint64_t> parents_;
+    std::vector<std::uint64_t> ranks_;
+};
+
+using MutexStorage = std::vector<std::unordered_set<std::uint64_t>>;
+
+inline bool check_mutex(
+    const std::uint64_t first,
+    const std::uint64_t second,
+    const MutexStorage &mutexes
+) {
+    const auto &first_mutexes = mutexes[first];
+    const auto &second_mutexes = mutexes[second];
+    if (first_mutexes.size() < second_mutexes.size()) {
+        return first_mutexes.find(second) != first_mutexes.end();
+    }
+    return second_mutexes.find(first) != second_mutexes.end();
+}
+
+inline void insert_mutex(
+    const std::uint64_t first,
+    const std::uint64_t second,
+    MutexStorage &mutexes
+) {
+    mutexes[first].insert(second);
+    mutexes[second].insert(first);
+}
+
+inline void merge_mutexes(
+    const std::uint64_t root_from,
+    const std::uint64_t root_to,
+    MutexStorage &mutexes
+) {
+    auto &mutexes_from = mutexes[root_from];
+    auto &mutexes_to = mutexes[root_to];
+
+    for (const auto other_root : mutexes_from) {
+        auto &other_mutexes = mutexes[other_root];
+        other_mutexes.erase(root_from);
+        if (other_root != root_to) {
+            other_mutexes.insert(root_to);
+            mutexes_to.insert(other_root);
+        }
+    }
+    mutexes_to.erase(root_from);
+    mutexes_to.erase(root_to);
+    mutexes_from.clear();
+}
+
+inline std::vector<std::ptrdiff_t> c_order_strides(const std::vector<std::ptrdiff_t> &shape) {
+    std::vector<std::ptrdiff_t> strides(shape.size(), 1);
+    for (std::ptrdiff_t axis = static_cast<std::ptrdiff_t>(shape.size()) - 2; axis >= 0; --axis) {
+        strides[static_cast<std::size_t>(axis)] =
+            strides[static_cast<std::size_t>(axis + 1)] * shape[static_cast<std::size_t>(axis + 1)];
+    }
+    return strides;
+}
+
+inline bool is_valid_grid_edge(
+    const std::uint64_t node,
+    const std::vector<std::ptrdiff_t> &offset,
+    const std::vector<std::ptrdiff_t> &shape,
+    const std::vector<std::ptrdiff_t> &strides
+) {
+    for (std::size_t axis = 0; axis < shape.size(); ++axis) {
+        const auto coord = static_cast<std::ptrdiff_t>(node / static_cast<std::uint64_t>(strides[axis])) % shape[axis];
+        const auto neighbor = coord + offset[axis];
+        if (neighbor < 0 || neighbor >= shape[axis]) {
+            return false;
+        }
+    }
+    return true;
+}
+
+template <class T>
+void mutex_watershed_grid(
+    const ConstArrayView<T> &affinities,
+    const ConstArrayView<std::uint8_t> &valid_edges,
+    const std::vector<std::vector<std::ptrdiff_t>> &offsets,
+    const std::size_t number_of_attractive_channels,
+    const ArrayView<std::uint64_t> &out
+) {
+    if (affinities.ndim() != 3 && affinities.ndim() != 4) {
+        throw std::invalid_argument(
+            "affinities must have shape (channels, y, x) or (channels, z, y, x), got ndim=" +
+            std::to_string(affinities.ndim())
+        );
+    }
+    if (offsets.empty()) {
+        throw std::invalid_argument("offsets must not be empty");
+    }
+    if (valid_edges.shape != affinities.shape) {
+        throw std::invalid_argument("valid_edges shape must match affinities shape");
+    }
+
+    const auto number_of_channels = static_cast<std::size_t>(affinities.shape[0]);
+    const auto spatial_ndim = static_cast<std::size_t>(affinities.ndim() - 1);
+    if (offsets.size() != number_of_channels) {
+        throw std::invalid_argument(
+            "offsets length must match affinities channel count, got offsets length=" +
+            std::to_string(offsets.size()) + ", channels=" + std::to_string(number_of_channels)
+        );
+    }
+    if (number_of_attractive_channels > number_of_channels) {
+        throw std::invalid_argument("number_of_attractive_channels must be <= number of channels");
+    }
+    for (const auto &offset : offsets) {
+        if (offset.size() != spatial_ndim) {
+            throw std::invalid_argument(
+                "each offset must have length matching the spatial ndim, got spatial ndim=" +
+                std::to_string(spatial_ndim)
+            );
+        }
+    }
+
+    std::vector<std::ptrdiff_t> spatial_shape(
+        affinities.shape.begin() + 1,
+        affinities.shape.end()
+    );
+    if (out.shape != spatial_shape) {
+        throw std::invalid_argument("out shape must match affinities spatial shape");
+    }
+
+    const auto number_of_nodes = static_cast<std::uint64_t>(std::accumulate(
+        spatial_shape.begin(),
+        spatial_shape.end(),
+        std::ptrdiff_t{1},
+        [](const std::ptrdiff_t a, const std::ptrdiff_t b) { return a * b; }
+    ));
+    const auto spatial_strides = c_order_strides(spatial_shape);
+
+    std::vector<std::ptrdiff_t> offset_strides(number_of_channels, 0);
+    for (std::size_t channel = 0; channel < number_of_channels; ++channel) {
+        for (std::size_t axis = 0; axis < spatial_ndim; ++axis) {
+            offset_strides[channel] += offsets[channel][axis] * spatial_strides[axis];
+        }
+    }
+
+    const auto number_of_edges = number_of_nodes * number_of_channels;
+    std::vector<std::uint64_t> edge_order(static_cast<std::size_t>(number_of_edges));
+    std::iota(edge_order.begin(), edge_order.end(), std::uint64_t{0});
+    std::sort(edge_order.begin(), edge_order.end(), [&](const std::uint64_t first, const std::uint64_t second) {
+        const T first_weight = affinities.data[first];
+        const T second_weight = affinities.data[second];
+        if (first_weight == second_weight) {
+            return first < second;
+        }
+        return first_weight > second_weight;
+    });
+
+    DisjointSets sets(static_cast<std::size_t>(number_of_nodes));
+    MutexStorage mutexes(static_cast<std::size_t>(number_of_nodes));
+
+    for (const auto edge_id : edge_order) {
+        if (valid_edges.data[edge_id] == 0) {
+            continue;
+        }
+
+        const auto channel = static_cast<std::size_t>(edge_id / number_of_nodes);
+        const auto u = edge_id % number_of_nodes;
+        if (!is_valid_grid_edge(u, offsets[channel], spatial_shape, spatial_strides)) {
+            continue;
+        }
+
+        const auto v_signed = static_cast<std::int64_t>(u) + static_cast<std::int64_t>(offset_strides[channel]);
+        const auto v = static_cast<std::uint64_t>(v_signed);
+        std::uint64_t root_u = sets.find(u);
+        std::uint64_t root_v = sets.find(v);
+        if (root_u == root_v || check_mutex(root_u, root_v, mutexes)) {
+            continue;
+        }
+
+        const bool is_mutex_edge = channel >= number_of_attractive_channels;
+        if (is_mutex_edge) {
+            insert_mutex(root_u, root_v, mutexes);
+        } else {
+            const auto new_root = sets.unite_roots(root_u, root_v);
+            const auto old_root = (new_root == root_u) ? root_v : root_u;
+            merge_mutexes(old_root, new_root, mutexes);
+        }
+    }
+
+    std::unordered_map<std::uint64_t, std::uint64_t> label_map;
+    for (std::uint64_t node = 0; node < number_of_nodes; ++node) {
+        const auto root = sets.find(node);
+        auto found = label_map.find(root);
+        if (found == label_map.end()) {
+            const auto next_label = static_cast<std::uint64_t>(label_map.size() + 1);
+            found = label_map.emplace(root, next_label).first;
+        }
+        out.data[node] = found->second;
+    }
+}
+
+} // namespace bioimage_cpp
diff --git a/src/bindings/module.cxx b/src/bindings/module.cxx
index bc3a0be..ac1230e 100644
--- a/src/bindings/module.cxx
+++ b/src/bindings/module.cxx
@@ -1,3 +1,4 @@
+#include "segmentation.hxx"
 #include "utils.hxx"
 
 #include <nanobind/nanobind.h>
@@ -6,5 +7,6 @@ namespace nb = nanobind;
 
 NB_MODULE(_core, m) {
     m.doc() = "C++ extension module for bioimage_cpp.";
+    bioimage_cpp::bindings::bind_segmentation(m);
     bioimage_cpp::bindings::bind_utils(m);
 }
diff --git a/src/bindings/segmentation.cxx b/src/bindings/segmentation.cxx
new file mode 100644
index 0000000..d37bdc4
--- /dev/null
+++ b/src/bindings/segmentation.cxx
@@ -0,0 +1,113 @@
+#include "segmentation.hxx"
+
+#include "bioimage_cpp/array_view.hxx"
+#include "bioimage_cpp/segmentation/mutex_watershed.hxx"
+
+#include <nanobind/ndarray.h>
+#include <nanobind/stl/vector.h>
+
+#include <cstddef>
+#include <cstdint>
+#include <memory>
+#include <numeric>
+#include <vector>
+
+namespace nb = nanobind;
+
+namespace bioimage_cpp::bindings {
+namespace {
+
+template <class T>
+using AffinityArray = nb::ndarray<nb::numpy, const T, nb::c_contig>;
+
+using ValidEdgeArray = nb::ndarray<nb::numpy, const std::uint8_t, nb::c_contig>;
+using LabelArray = nb::ndarray<nb::numpy, std::uint64_t, nb::c_contig>;
+
+template <class T>
+LabelArray mutex_watershed_grid_t(
+    AffinityArray<T> affinities,
+    ValidEdgeArray valid_edges,
+    const std::vector<std::vector<std::ptrdiff_t>> &offsets,
+    const std::size_t number_of_attractive_channels
+) {
+    std::vector<std::ptrdiff_t> affinity_shape(affinities.ndim());
+    for (std::size_t axis = 0; axis < affinities.ndim(); ++axis) {
+        affinity_shape[axis] = static_cast<std::ptrdiff_t>(affinities.shape(axis));
+    }
+    std::vector<std::ptrdiff_t> valid_edges_shape(valid_edges.ndim());
+    for (std::size_t axis = 0; axis < valid_edges.ndim(); ++axis) {
+        valid_edges_shape[axis] = static_cast<std::ptrdiff_t>(valid_edges.shape(axis));
+    }
+
+    std::vector<std::size_t> label_shape;
+    label_shape.reserve(affinities.ndim() > 0 ? affinities.ndim() - 1 : 0);
+    std::vector<std::ptrdiff_t> label_view_shape;
+    label_view_shape.reserve(affinities.ndim() > 0 ? affinities.ndim() - 1 : 0);
+    for (std::size_t axis = 1; axis < affinities.ndim(); ++axis) {
+        label_shape.push_back(affinities.shape(axis));
+        label_view_shape.push_back(static_cast<std::ptrdiff_t>(affinities.shape(axis)));
+    }
+
+    const auto number_of_nodes = std::accumulate(
+        label_view_shape.begin(),
+        label_view_shape.end(),
+        std::ptrdiff_t{1},
+        [](const std::ptrdiff_t a, const std::ptrdiff_t b) { return a * b; }
+    );
+    auto *data = new std::uint64_t[static_cast<std::size_t>(number_of_nodes)]();
+    nb::capsule owner(data, [](void *p) noexcept { delete[] static_cast<std::uint64_t *>(p); });
+
+    ConstArrayView<T> affinities_view{
+        affinities.data(),
+        affinity_shape,
+        {},
+    };
+    ConstArrayView<std::uint8_t> valid_edges_view{
+        valid_edges.data(),
+        valid_edges_shape,
+        {},
+    };
+    ArrayView<std::uint64_t> out_view{
+        data,
+        label_view_shape,
+        {},
+    };
+
+    {
+        nb::gil_scoped_release release;
+        mutex_watershed_grid<T>(
+            affinities_view,
+            valid_edges_view,
+            offsets,
+            number_of_attractive_channels,
+            out_view
+        );
+    }
+
+    return LabelArray(data, label_shape.size(), label_shape.data(), owner);
+}
+
+} // namespace
+
+void bind_segmentation(nb::module_ &m) {
+    m.def(
+        "_mutex_watershed_grid_float32",
+        &mutex_watershed_grid_t<float>,
+        nb::arg("affinities"),
+        nb::arg("valid_edges"),
+        nb::arg("offsets"),
+        nb::arg("number_of_attractive_channels"),
+        "Run mutex watershed on a 2D or 3D image-derived grid graph with float32 affinities."
+    );
+    m.def(
+        "_mutex_watershed_grid_float64",
+        &mutex_watershed_grid_t<double>,
+        nb::arg("affinities"),
+        nb::arg("valid_edges"),
+        nb::arg("offsets"),
+        nb::arg("number_of_attractive_channels"),
+        "Run mutex watershed on a 2D or 3D image-derived grid graph with float64 affinities."
+    );
+}
+
+} // namespace bioimage_cpp::bindings
diff --git a/src/bindings/segmentation.hxx b/src/bindings/segmentation.hxx
new file mode 100644
index 0000000..8502e84
--- /dev/null
+++ b/src/bindings/segmentation.hxx
@@ -0,0 +1,9 @@
+#pragma once
+
+#include <nanobind/nanobind.h>
+
+namespace bioimage_cpp::bindings {
+
+void bind_segmentation(nanobind::module_ &m);
+
+} // namespace bioimage_cpp::bindings
diff --git a/src/bioimage_cpp/__init__.py b/src/bioimage_cpp/__init__.py
index 0b5ebdd..ba81546 100644
--- a/src/bioimage_cpp/__init__.py
+++ b/src/bioimage_cpp/__init__.py
@@ -1,6 +1,7 @@
 """Dependency-light bioimage analysis algorithms backed by C++."""
 
 from ._version import __version__
+from . import segmentation
 from . import utils
 
-__all__ = ["__version__", "utils"]
+__all__ = ["__version__", "segmentation", "utils"]
diff --git a/src/bioimage_cpp/segmentation/__init__.py b/src/bioimage_cpp/segmentation/__init__.py
new file mode 100644
index 0000000..cf4b303
--- /dev/null
+++ b/src/bioimage_cpp/segmentation/__init__.py
@@ -0,0 +1,5 @@
+"""Segmentation algorithms."""
+
+from .mutex_watershed import mutex_watershed
+
+__all__ = ["mutex_watershed"]
diff --git a/src/bioimage_cpp/segmentation/mutex_watershed.py b/src/bioimage_cpp/segmentation/mutex_watershed.py
new file mode 100644
index 0000000..d500dd0
--- /dev/null
+++ b/src/bioimage_cpp/segmentation/mutex_watershed.py
@@ -0,0 +1,227 @@
+"""Segmentation algorithms."""
+
+from __future__ import annotations
+
+from collections.abc import Sequence
+
+import numpy as np
+
+from .. import _core
+
+_MUTEX_WATERSHED_BY_DTYPE = {
+    np.dtype("float32"): _core._mutex_watershed_grid_float32,
+    np.dtype("float64"): _core._mutex_watershed_grid_float64,
+}
+
+
+def mutex_watershed(
+    affinities: np.ndarray,
+    offsets: Sequence[Sequence[int]],
+    number_of_attractive_channels: int,
+    *,
+    strides: Sequence[int] | None = None,
+    randomized_strides: bool = False,
+    mask: np.ndarray | None = None,
+) -> np.ndarray:
+    """Run mutex watershed on a 2D or 3D image-derived grid graph.
+
+    Parameters
+    ----------
+    affinities:
+        Array with shape ``(channels, y, x)`` for 2D data or
+        ``(channels, z, y, x)`` for 3D data. Supported dtypes are
+        ``float32`` and ``float64``. Non-contiguous arrays are copied.
+    offsets:
+        One offset per channel, in NumPy axis order. Each offset has length 2
+        for 2D affinities or length 3 for 3D affinities.
+    number_of_attractive_channels:
+        The first this many affinity channels are attractive merge edges. The
+        remaining channels are mutex edges.
+    strides:
+        Optional spatial sub-sampling strides for mutex edges. Attractive
+        channels are always kept. If given, it must have one positive integer
+        per spatial dimension.
+    randomized_strides:
+        If ``True``, sub-sample mutex edges randomly with probability
+        ``1 / np.prod(strides)`` instead of on a regular grid.
+    mask:
+        Optional boolean foreground mask with shape ``affinities.shape[1:]``.
+        Edges touching ``False`` pixels are ignored and ``False`` pixels are
+        labeled as background ``0`` in the output.
+
+    Returns
+    -------
+    np.ndarray
+        Consecutive 1-based ``uint64`` segmentation labels with shape
+        ``affinities.shape[1:]``.
+    """
+    array = np.asarray(affinities)
+    if array.ndim not in (3, 4):
+        raise ValueError(
+            "affinities must have shape (channels, y, x) or "
+            f"(channels, z, y, x), got ndim={array.ndim}"
+        )
+
+    dtype = array.dtype
+    try:
+        run = _MUTEX_WATERSHED_BY_DTYPE[dtype]
+    except KeyError as error:
+        supported = ", ".join(str(dtype) for dtype in _MUTEX_WATERSHED_BY_DTYPE)
+        raise TypeError(
+            f"affinities must have one of dtypes ({supported}), got dtype={dtype}"
+        ) from error
+
+    normalized_offsets = [tuple(int(value) for value in offset) for offset in offsets]
+    spatial_ndim = array.ndim - 1
+    if len(normalized_offsets) != array.shape[0]:
+        raise ValueError(
+            "offsets length must match affinities channel count, got "
+            f"offsets length={len(normalized_offsets)}, channels={array.shape[0]}"
+        )
+    if any(len(offset) != spatial_ndim for offset in normalized_offsets):
+        raise ValueError(
+            "each offset must have length matching the spatial ndim, got "
+            f"spatial ndim={spatial_ndim}"
+        )
+    if number_of_attractive_channels < 0:
+        raise ValueError("number_of_attractive_channels must be non-negative")
+    if number_of_attractive_channels > array.shape[0]:
+        raise ValueError("number_of_attractive_channels must be <= number of channels")
+
+    normalized_strides = _normalize_strides(strides, spatial_ndim, randomized_strides)
+    valid_edges = _compute_valid_edges(
+        array.shape,
+        normalized_offsets,
+        number_of_attractive_channels,
+        normalized_strides,
+        randomized_strides,
+        mask,
+    )
+
+    contiguous = np.ascontiguousarray(array)
+    labels = run(
+        contiguous,
+        valid_edges,
+        normalized_offsets,
+        number_of_attractive_channels,
+    )
+    if mask is not None:
+        labels[np.logical_not(np.asarray(mask))] = 0
+    return labels
+
+
+def _normalize_strides(
+    strides: Sequence[int] | None,
+    spatial_ndim: int,
+    randomized_strides: bool,
+) -> tuple[int, ...] | None:
+    if strides is None:
+        if randomized_strides:
+            raise ValueError("randomized_strides requires strides")
+        return None
+
+    normalized = tuple(int(stride) for stride in strides)
+    if len(normalized) != spatial_ndim:
+        raise ValueError(
+            "strides length must match the spatial ndim, got "
+            f"strides length={len(normalized)}, spatial ndim={spatial_ndim}"
+        )
+    if any(stride <= 0 for stride in normalized):
+        raise ValueError("strides must contain only positive integers")
+    return normalized
+
+
+def _valid_source_slices(
+    image_shape: tuple[int, ...],
+    offset: tuple[int, ...],
+) -> tuple[slice, ...] | None:
+    slices = []
+    for axis_size, step in zip(image_shape, offset, strict=True):
+        if step > 0:
+            if step >= axis_size:
+                return None
+            slices.append(slice(0, axis_size - step))
+        elif step < 0:
+            if -step >= axis_size:
+                return None
+            slices.append(slice(-step, axis_size))
+        else:
+            slices.append(slice(None))
+    return tuple(slices)
+
+
+def _neighbor_slices(
+    image_shape: tuple[int, ...],
+    offset: tuple[int, ...],
+) -> tuple[slice, ...] | None:
+    slices = []
+    for axis_size, step in zip(image_shape, offset, strict=True):
+        if step > 0:
+            if step >= axis_size:
+                return None
+            slices.append(slice(step, axis_size))
+        elif step < 0:
+            if -step >= axis_size:
+                return None
+            slices.append(slice(0, axis_size + step))
+        else:
+            slices.append(slice(None))
+    return tuple(slices)
+
+
+def _compute_valid_edges(
+    affinity_shape: tuple[int, ...],
+    offsets: Sequence[tuple[int, ...]],
+    number_of_attractive_channels: int,
+    strides: tuple[int, ...] | None,
+    randomized_strides: bool,
+    mask: np.ndarray | None,
+) -> np.ndarray:
+    image_shape = tuple(int(size) for size in affinity_shape[1:])
+    valid_edges = np.zeros(affinity_shape, dtype=bool)
+
+    for channel, offset in enumerate(offsets):
+        source_slices = _valid_source_slices(image_shape, offset)
+        if source_slices is not None:
+            valid_edges[(channel,) + source_slices] = True
+
+    if strides is not None:
+        stride_edges = np.zeros_like(valid_edges, dtype=bool)
+        stride_edges[:number_of_attractive_channels] = True
+        if randomized_strides:
+            stride_factor = 1.0 / np.prod(strides)
+            stride_edges[number_of_attractive_channels:] = (
+                np.random.random(
+                    valid_edges[number_of_attractive_channels:].shape
+                )
+                < stride_factor
+            )
+        else:
+            valid_slice = (slice(number_of_attractive_channels, None),) + tuple(
+                slice(None, None, stride) for stride in strides
+            )
+            stride_edges[valid_slice] = True
+        valid_edges &= stride_edges
+
+    if mask is not None:
+        mask_array = np.asarray(mask)
+        if mask_array.shape != image_shape:
+            raise ValueError(
+                "mask shape must match affinities spatial shape, got "
+                f"mask shape={mask_array.shape}, spatial shape={image_shape}"
+            )
+        if mask_array.dtype != np.dtype("bool"):
+            raise TypeError(f"mask must have dtype bool, got dtype={mask_array.dtype}")
+
+        invalid = np.logical_not(mask_array)
+        for channel, offset in enumerate(offsets):
+            source_slices = _valid_source_slices(image_shape, offset)
+            neighbor_slices = _neighbor_slices(image_shape, offset)
+            if source_slices is None or neighbor_slices is None:
+                continue
+            touches_invalid = invalid[source_slices] | invalid[neighbor_slices]
+            channel_valid = valid_edges[(channel,) + source_slices]
+            channel_valid[touches_invalid] = False
+            valid_edges[(channel,) + source_slices] = channel_valid
+
+    return np.ascontiguousarray(valid_edges, dtype=np.uint8)
diff --git a/src/cpp/segmentation/mutex_watershed.cxx b/src/cpp/segmentation/mutex_watershed.cxx
new file mode 100644
index 0000000..cfac854
--- /dev/null
+++ b/src/cpp/segmentation/mutex_watershed.cxx
@@ -0,0 +1,23 @@
+#include "bioimage_cpp/segmentation/mutex_watershed.hxx"
+
+#include <cstdint>
+
+namespace bioimage_cpp {
+
+template void mutex_watershed_grid<float>(
+    const ConstArrayView<float> &,
+    const ConstArrayView<std::uint8_t> &,
+    const std::vector<std::vector<std::ptrdiff_t>> &,
+    std::size_t,
+    const ArrayView<std::uint64_t> &
+);
+
+template void mutex_watershed_grid<double>(
+    const ConstArrayView<double> &,
+    const ConstArrayView<std::uint8_t> &,
+    const std::vector<std::vector<std::ptrdiff_t>> &,
+    std::size_t,
+    const ArrayView<std::uint64_t> &
+);
+
+} // namespace bioimage_cpp
diff --git a/tests/segmentation/test_mutex_watershed.py b/tests/segmentation/test_mutex_watershed.py
new file mode 100644
index 0000000..92640c8
--- /dev/null
+++ b/tests/segmentation/test_mutex_watershed.py
@@ -0,0 +1,179 @@
+import numpy as np
+import pytest
+
+import bioimage_cpp as bic
+
+
+def test_mutex_watershed_2d_attractive_edges_merge_all_pixels():
+    affinities = np.ones((2, 3, 4), dtype=np.float32)
+    offsets = [[0, 1], [1, 0]]
+
+    labels = bic.segmentation.mutex_watershed(
+        affinities, offsets, number_of_attractive_channels=2
+    )
+
+    assert labels.dtype == np.uint64
+    assert labels.shape == (3, 4)
+    np.testing.assert_array_equal(labels, np.ones((3, 4), dtype=np.uint64))
+
+
+def test_mutex_watershed_2d_mutex_edge_blocks_lower_attractive_edge():
+    affinities = np.array([[[0.5, 0.0]], [[0.9, 0.0]]], dtype=np.float64)
+    offsets = [[0, 1], [0, 1]]
+
+    labels = bic.segmentation.mutex_watershed(
+        affinities, offsets, number_of_attractive_channels=1
+    )
+
+    np.testing.assert_array_equal(labels, np.array([[1, 2]], dtype=np.uint64))
+
+
+def test_mutex_watershed_2d_higher_attractive_edge_wins_before_mutex_edge():
+    affinities = np.array([[[0.9, 0.0]], [[0.5, 0.0]]], dtype=np.float32)
+    offsets = [[0, 1], [0, 1]]
+
+    labels = bic.segmentation.mutex_watershed(
+        affinities, offsets, number_of_attractive_channels=1
+    )
+
+    np.testing.assert_array_equal(labels, np.array([[1, 1]], dtype=np.uint64))
+
+
+def test_mutex_watershed_2d_respects_grid_boundaries():
+    affinities = np.ones((1, 2, 3), dtype=np.float32)
+    offsets = [[0, 1]]
+
+    labels = bic.segmentation.mutex_watershed(
+        affinities, offsets, number_of_attractive_channels=1
+    )
+
+    np.testing.assert_array_equal(
+        labels, np.array([[1, 1, 1], [2, 2, 2]], dtype=np.uint64)
+    )
+
+
+def test_mutex_watershed_2d_accepts_negative_offsets():
+    affinities = np.ones((1, 2, 3), dtype=np.float32)
+    offsets = [[0, -1]]
+
+    labels = bic.segmentation.mutex_watershed(
+        affinities, offsets, number_of_attractive_channels=1
+    )
+
+    np.testing.assert_array_equal(
+        labels, np.array([[1, 1, 1], [2, 2, 2]], dtype=np.uint64)
+    )
+
+
+def test_mutex_watershed_3d_attractive_edges_merge_all_pixels():
+    affinities = np.ones((3, 2, 2, 2), dtype=np.float32)
+    offsets = [[1, 0, 0], [0, 1, 0], [0, 0, 1]]
+
+    labels = bic.segmentation.mutex_watershed(
+        affinities, offsets, number_of_attractive_channels=3
+    )
+
+    assert labels.shape == (2, 2, 2)
+    np.testing.assert_array_equal(labels, np.ones((2, 2, 2), dtype=np.uint64))
+
+
+def test_mutex_watershed_strides_subsample_mutex_edges():
+    affinities = np.zeros((2, 1, 5), dtype=np.float32)
+    affinities[0, 0, :] = 0.8
+    affinities[1, 0, :] = 0.9
+    offsets = [[0, 1], [0, 1]]
+
+    dense = bic.segmentation.mutex_watershed(
+        affinities, offsets, number_of_attractive_channels=1
+    )
+    strided = bic.segmentation.mutex_watershed(
+        affinities, offsets, number_of_attractive_channels=1, strides=[1, 2]
+    )
+
+    np.testing.assert_array_equal(dense, np.array([[1, 2, 3, 4, 5]], dtype=np.uint64))
+    np.testing.assert_array_equal(strided, np.array([[1, 2, 2, 3, 3]], dtype=np.uint64))
+
+
+def test_mutex_watershed_randomized_strides_use_numpy_random_state():
+    affinities = np.zeros((2, 1, 5), dtype=np.float32)
+    affinities[0, 0, :] = 0.8
+    affinities[1, 0, :] = 0.9
+    offsets = [[0, 1], [0, 1]]
+
+    np.random.seed(17)
+    first = bic.segmentation.mutex_watershed(
+        affinities,
+        offsets,
+        number_of_attractive_channels=1,
+        strides=[1, 2],
+        randomized_strides=True,
+    )
+    np.random.seed(17)
+    second = bic.segmentation.mutex_watershed(
+        affinities,
+        offsets,
+        number_of_attractive_channels=1,
+        strides=[1, 2],
+        randomized_strides=True,
+    )
+
+    np.testing.assert_array_equal(first, second)
+
+
+def test_mutex_watershed_mask_sets_background_zero_and_blocks_edges():
+    affinities = np.ones((1, 1, 5), dtype=np.float32)
+    offsets = [[0, 1]]
+    mask = np.array([[True, True, False, True, True]], dtype=bool)
+
+    labels = bic.segmentation.mutex_watershed(
+        affinities, offsets, number_of_attractive_channels=1, mask=mask
+    )
+
+    np.testing.assert_array_equal(labels, np.array([[1, 1, 0, 3, 3]], dtype=np.uint64))
+
+
+def test_mutex_watershed_rejects_wrong_affinity_dtype():
+    with pytest.raises(TypeError, match="affinities must have one of dtypes"):
+        bic.segmentation.mutex_watershed(
+            np.ones((1, 2, 2), dtype=np.uint8),
+            [[0, 1]],
+            number_of_attractive_channels=1,
+        )
+
+
+def test_mutex_watershed_rejects_wrong_offset_length():
+    with pytest.raises(ValueError, match="each offset must have length"):
+        bic.segmentation.mutex_watershed(
+            np.ones((1, 2, 2), dtype=np.float32),
+            [[0, 1, 0]],
+            number_of_attractive_channels=1,
+        )
+
+
+def test_mutex_watershed_rejects_wrong_number_of_offsets():
+    with pytest.raises(ValueError, match="offsets length must match"):
+        bic.segmentation.mutex_watershed(
+            np.ones((2, 2, 2), dtype=np.float32),
+            [[0, 1]],
+            number_of_attractive_channels=1,
+        )
+
+
+def test_mutex_watershed_rejects_invalid_strides():
+    with pytest.raises(ValueError, match="strides length must match"):
+        bic.segmentation.mutex_watershed(
+            np.ones((1, 2, 2), dtype=np.float32),
+            [[0, 1]],
+            number_of_attractive_channels=1,
+            strides=[1, 1, 1],
+        )
+
+
+def test_mutex_watershed_rejects_invalid_mask():
+    with pytest.raises(TypeError, match="mask must have dtype bool"):
+        bic.segmentation.mutex_watershed(
+            np.ones((1, 2, 2), dtype=np.float32),
+            [[0, 1]],
+            number_of_attractive_channels=1,
+            mask=np.ones((2, 2), dtype=np.uint8),
+        )