diff --git a/MIGRATION_GUIDE.md b/MIGRATION_GUIDE.md
index 99f5dad..7c4277f 100644
--- a/MIGRATION_GUIDE.md
+++ b/MIGRATION_GUIDE.md
@@ -899,6 +899,18 @@ Important differences:
   `grid_affinity_features_with_lifted(...)`, which returns local graph weights
   plus explicit long-range `uv_ids` and weights for lifted multicut or mutex
   watershed style workflows.
+- Nifty's `projectEdgeIdsToPixels` and `projectEdgeIdsToPixelsWithOffsets` map
+  to `graph.project_edge_ids_to_pixels()` and
+  `graph.project_edge_ids_to_pixels_with_offsets(offsets, *, strides=None, mask=None)`
+  on `GridGraph2D` / `GridGraph3D`. The basic form returns an `int64` array of
+  shape `(ndim, *graph.shape)` with each grid edge id written at its pivot
+  pixel and `-1` elsewhere. The offsets form returns
+  `(array, n_valid)`: an `int64` array of shape `(len(offsets), *graph.shape)`
+  whose non-`-1` entries are a sequential counter over the in-bounds (and
+  filter-accepted) targets, plus the total count. `strides` keeps only coords
+  aligned along every axis; `mask` keeps only coords where a boolean array
+  of shape `(len(offsets), *graph.shape)` is true. Like in nifty, `strides`
+  and `mask` are mutually exclusive — passing both raises `ValueError`.
 - The three grid feature functions preserve `float32` and `float64` input
   dtype end-to-end (no internal copy to `float64`); other dtypes are
   promoted to `float64`. Output weight arrays match the input dtype.
diff --git a/include/bioimage_cpp/graph/grid_edge_projection.hxx b/include/bioimage_cpp/graph/grid_edge_projection.hxx
new file mode 100644
index 0000000..19770cd
--- /dev/null
+++ b/include/bioimage_cpp/graph/grid_edge_projection.hxx
@@ -0,0 +1,276 @@
+#pragma once
+
+#include "bioimage_cpp/array_view.hxx"
+#include "bioimage_cpp/graph/grid_graph.hxx"
+
+#include <array>
+#include <cstddef>
+#include <cstdint>
+#include <stdexcept>
+#include <vector>
+
+namespace bioimage_cpp::graph {
+
+namespace detail_grid_edge_projection {
+
+// Walk a sub-shape in C-order, computing a flat output offset incrementally
+// using `strides` (the strides of the *full* containing array, not of the
+// sub-shape itself). Calls `callback(flat)` at every leaf. With `Axis` known
+// at compile time the compiler unrolls the entire nest for D = 2 / D = 3.
+template <std::size_t Axis, std::size_t D, class Callback>
+void enumerate_in_c_order(
+    const std::array<std::uint64_t, D> &subshape,
+    const std::array<std::ptrdiff_t, D> &strides,
+    std::ptrdiff_t base,
+    Callback &&callback
+) {
+    if constexpr (Axis == D) {
+        callback(base);
+    } else {
+        for (std::uint64_t i = 0; i < subshape[Axis]; ++i) {
+            enumerate_in_c_order<Axis + 1, D>(subshape, strides, base, callback);
+            base += strides[Axis];
+        }
+    }
+}
+
+// Walk every coordinate in C-order over `shape`, tracking whether
+// `coord + offset` is in bounds, and call `callback(in_bounds, coord)` at the
+// leaf. `coord` is filled in-place and re-used across calls so the sampler
+// can read it without copying.
+template <std::size_t Axis, std::size_t D, class Callback>
+void enumerate_with_offset(
+    const std::array<std::uint64_t, D> &shape,
+    const std::array<std::ptrdiff_t, D> &offset,
+    std::array<std::ptrdiff_t, D> &coord,
+    bool in_bounds,
+    Callback &&callback
+) {
+    if constexpr (Axis == D) {
+        callback(in_bounds, coord);
+    } else {
+        for (std::uint64_t c = 0; c < shape[Axis]; ++c) {
+            coord[Axis] = static_cast<std::ptrdiff_t>(c);
+            const std::ptrdiff_t target = coord[Axis] + offset[Axis];
+            const bool axis_ok =
+                target >= 0 &&
+                target < static_cast<std::ptrdiff_t>(shape[Axis]);
+            enumerate_with_offset<Axis + 1, D>(
+                shape, offset, coord, in_bounds && axis_ok, callback
+            );
+        }
+    }
+}
+
+template <std::size_t D>
+std::array<std::ptrdiff_t, D> full_c_strides(
+    const std::array<std::uint64_t, D> &shape
+) {
+    std::array<std::ptrdiff_t, D> strides{};
+    std::ptrdiff_t s = 1;
+    for (std::size_t i = D; i-- > 0; ) {
+        strides[i] = s;
+        s *= static_cast<std::ptrdiff_t>(shape[i]);
+    }
+    return strides;
+}
+
+template <std::size_t D>
+std::ptrdiff_t shape_product(const std::array<std::uint64_t, D> &shape) {
+    std::ptrdiff_t p = 1;
+    for (std::size_t i = 0; i < D; ++i) {
+        p *= static_cast<std::ptrdiff_t>(shape[i]);
+    }
+    return p;
+}
+
+template <std::size_t D>
+void require_projection_output_shape(
+    const std::array<std::uint64_t, D> &graph_shape,
+    const std::vector<std::ptrdiff_t> &out_shape,
+    const std::ptrdiff_t expected_leading,
+    const char *argument_name
+) {
+    if (out_shape.size() != D + 1) {
+        throw std::invalid_argument(
+            std::string(argument_name) + " must have ndim == graph.ndim + 1"
+        );
+    }
+    if (out_shape[0] != expected_leading) {
+        throw std::invalid_argument(
+            std::string(argument_name) + " leading dimension does not match expected length"
+        );
+    }
+    for (std::size_t d = 0; d < D; ++d) {
+        if (out_shape[d + 1] != static_cast<std::ptrdiff_t>(graph_shape[d])) {
+            throw std::invalid_argument(
+                std::string(argument_name) + " spatial shape must match graph shape"
+            );
+        }
+    }
+}
+
+} // namespace detail_grid_edge_projection
+
+// Writes graph edge ids into `output` of shape (D, *graph.shape) and dtype
+// int64. The caller pre-fills `output` with -1; this function then sets
+// output[axis, c_0, ..., c_{D-1}] = edge_id for every edge whose spanning
+// axis is `axis` and whose pivot (smaller-endpoint) coordinate is
+// (c_0, ..., c_{D-1}). Slots where c_axis == shape[axis] - 1 are not visited
+// and keep their pre-filled -1.
+//
+// Exploits the fact that the GridGraph stores edges axis-major + C-order
+// within each axis: the k-th edge of axis `d` corresponds exactly to the
+// k-th pivot coordinate of the sub-shape (s_0, ..., s_d - 1, ..., s_{D-1})
+// walked in C-order. So we walk that sub-shape with the *full* output
+// strides and emit consecutive edge ids — no per-edge coordinate division.
+template <std::size_t D>
+void project_edge_ids_to_pixels(
+    const GridGraph<D> &graph,
+    const ArrayView<std::int64_t> &output
+) {
+    namespace dgep = detail_grid_edge_projection;
+    const auto &shape = graph.shape();
+    dgep::require_projection_output_shape<D>(
+        shape, output.shape, static_cast<std::ptrdiff_t>(D), "output"
+    );
+
+    const auto strides = dgep::full_c_strides<D>(shape);
+    const std::ptrdiff_t plane = dgep::shape_product<D>(shape);
+
+    auto *data = output.data;
+    std::int64_t edge_id = 0;
+    for (std::size_t axis = 0; axis < D; ++axis) {
+        auto pivot_shape = shape;
+        if (pivot_shape[axis] == 0) {
+            continue;
+        }
+        --pivot_shape[axis];
+        const std::ptrdiff_t base =
+            static_cast<std::ptrdiff_t>(axis) * plane;
+        dgep::enumerate_in_c_order<0, D>(
+            pivot_shape, strides, base,
+            [data, &edge_id](std::ptrdiff_t flat) {
+                data[flat] = edge_id++;
+            }
+        );
+    }
+}
+
+// Enumerates "lifted" edges defined by a set of per-channel grid offsets.
+// Walks (offset_idx, *coord) in C-order over (offsets.size(), *graph.shape).
+// For each coord whose target `coord + offsets[offset_idx]` stays in bounds
+// AND the sampler accepts, writes a sequential counter starting at 0;
+// rejected slots get -1. Returns the total number of valid entries written.
+//
+// The counter is NOT a graph edge id — it indexes into the implicit array
+// of lifted edges derived from an affinity volume.
+template <std::size_t D, class Sampler>
+std::uint64_t project_edge_ids_to_pixels_with_offsets_impl(
+    const GridGraph<D> &graph,
+    const std::vector<std::array<std::ptrdiff_t, D>> &offsets,
+    Sampler &&sampler,
+    const ArrayView<std::int64_t> &output
+) {
+    namespace dgep = detail_grid_edge_projection;
+    const auto &shape = graph.shape();
+    dgep::require_projection_output_shape<D>(
+        shape, output.shape, static_cast<std::ptrdiff_t>(offsets.size()), "output"
+    );
+
+    auto *out = output.data;
+    std::int64_t edge_id = 0;
+    std::size_t out_idx = 0;
+    std::array<std::ptrdiff_t, D> coord{};
+    for (std::size_t off_idx = 0; off_idx < offsets.size(); ++off_idx) {
+        const auto &off = offsets[off_idx];
+        dgep::enumerate_with_offset<0, D>(
+            shape, off, coord, /*in_bounds=*/true,
+            [&](bool in_bounds, const std::array<std::ptrdiff_t, D> &c) {
+                if (in_bounds && sampler(off_idx, c)) {
+                    out[out_idx++] = edge_id++;
+                } else {
+                    out[out_idx++] = -1;
+                }
+            }
+        );
+    }
+    return static_cast<std::uint64_t>(edge_id);
+}
+
+// Convenience overload: no filter.
+template <std::size_t D>
+std::uint64_t project_edge_ids_to_pixels_with_offsets(
+    const GridGraph<D> &graph,
+    const std::vector<std::array<std::ptrdiff_t, D>> &offsets,
+    const ArrayView<std::int64_t> &output
+) {
+    return project_edge_ids_to_pixels_with_offsets_impl<D>(
+        graph, offsets,
+        [](std::size_t, const std::array<std::ptrdiff_t, D> &) { return true; },
+        output
+    );
+}
+
+// Convenience overload: only coords aligned with `strides` along every axis.
+template <std::size_t D>
+std::uint64_t project_edge_ids_to_pixels_with_offsets(
+    const GridGraph<D> &graph,
+    const std::vector<std::array<std::ptrdiff_t, D>> &offsets,
+    const std::array<std::ptrdiff_t, D> &strides,
+    const ArrayView<std::int64_t> &output
+) {
+    for (std::size_t d = 0; d < D; ++d) {
+        if (strides[d] <= 0) {
+            throw std::invalid_argument("strides must be positive");
+        }
+    }
+    return project_edge_ids_to_pixels_with_offsets_impl<D>(
+        graph, offsets,
+        [&strides](std::size_t, const std::array<std::ptrdiff_t, D> &c) {
+            for (std::size_t d = 0; d < D; ++d) {
+                if (c[d] % strides[d] != 0) {
+                    return false;
+                }
+            }
+            return true;
+        },
+        output
+    );
+}
+
+// Convenience overload: only coords where the (n_offsets, *graph.shape)
+// `mask` is non-zero.
+template <std::size_t D>
+std::uint64_t project_edge_ids_to_pixels_with_offsets(
+    const GridGraph<D> &graph,
+    const std::vector<std::array<std::ptrdiff_t, D>> &offsets,
+    const ConstArrayView<std::uint8_t> &mask,
+    const ArrayView<std::int64_t> &output
+) {
+    namespace dgep = detail_grid_edge_projection;
+    const auto &shape = graph.shape();
+    dgep::require_projection_output_shape<D>(
+        shape, mask.shape, static_cast<std::ptrdiff_t>(offsets.size()), "mask"
+    );
+
+    const auto mask_spatial_strides = dgep::full_c_strides<D>(shape);
+    const std::ptrdiff_t mask_plane = dgep::shape_product<D>(shape);
+    const auto *mdata = mask.data;
+    return project_edge_ids_to_pixels_with_offsets_impl<D>(
+        graph, offsets,
+        [mdata, mask_spatial_strides, mask_plane](
+            std::size_t off_idx, const std::array<std::ptrdiff_t, D> &c
+        ) {
+            std::ptrdiff_t flat =
+                static_cast<std::ptrdiff_t>(off_idx) * mask_plane;
+            for (std::size_t d = 0; d < D; ++d) {
+                flat += c[d] * mask_spatial_strides[d];
+            }
+            return mdata[flat] != 0;
+        },
+        output
+    );
+}
+
+} // namespace bioimage_cpp::graph
diff --git a/src/bindings/graph.cxx b/src/bindings/graph.cxx
index 5cd0f99..7109879 100644
--- a/src/bindings/graph.cxx
+++ b/src/bindings/graph.cxx
@@ -5,6 +5,7 @@
 #include "bioimage_cpp/graph/connected_components.hxx"
 #include "bioimage_cpp/graph/edge_weighted_watershed.hxx"
 #include "bioimage_cpp/graph/feature_accumulation.hxx"
+#include "bioimage_cpp/graph/grid_edge_projection.hxx"
 #include "bioimage_cpp/graph/grid_features.hxx"
 #include "bioimage_cpp/graph/grid_graph.hxx"
 #include "bioimage_cpp/graph/label_accumulation.hxx"
@@ -279,6 +280,120 @@ std::int64_t grid_offset_target(
     return static_cast<std::int64_t>(target);
 }
 
+template <std::size_t D>
+Int64Array grid_project_edge_ids_to_pixels(const graph::GridGraph<D> &graph) {
+    const auto &shape = graph.shape();
+    std::vector<std::size_t> out_shape(D + 1);
+    out_shape[0] = D;
+    for (std::size_t d = 0; d < D; ++d) {
+        out_shape[d + 1] = shape[d];
+    }
+    auto result = make_int64_array(out_shape);
+
+    std::size_t total = 1;
+    for (const auto axis_size : out_shape) {
+        total *= axis_size;
+    }
+    auto *data = result.data();
+    std::fill(data, data + total, static_cast<std::int64_t>(-1));
+
+    std::vector<std::ptrdiff_t> view_shape(out_shape.size());
+    for (std::size_t i = 0; i < out_shape.size(); ++i) {
+        view_shape[i] = static_cast<std::ptrdiff_t>(out_shape[i]);
+    }
+    ArrayView<std::int64_t> view{data, view_shape, {}};
+
+    {
+        nb::gil_scoped_release release;
+        graph::project_edge_ids_to_pixels<D>(graph, view);
+    }
+    return result;
+}
+
+template <std::size_t D>
+nb::tuple grid_project_edge_ids_to_pixels_with_offsets(
+    const graph::GridGraph<D> &graph,
+    const std::vector<std::vector<std::ptrdiff_t>> &offsets,
+    const std::optional<std::vector<std::ptrdiff_t>> &strides,
+    const std::optional<ConstUInt8Array> &mask
+) {
+    if (strides.has_value() && mask.has_value()) {
+        throw std::invalid_argument("strides and mask cannot be given together");
+    }
+
+    std::vector<std::array<std::ptrdiff_t, D>> off_arr(offsets.size());
+    for (std::size_t i = 0; i < offsets.size(); ++i) {
+        if (offsets[i].size() != D) {
+            throw std::invalid_argument(
+                "each offset must have length matching graph ndim"
+            );
+        }
+        for (std::size_t d = 0; d < D; ++d) {
+            off_arr[i][d] = offsets[i][d];
+        }
+    }
+
+    const auto &shape = graph.shape();
+    std::vector<std::size_t> out_shape(D + 1);
+    out_shape[0] = offsets.size();
+    for (std::size_t d = 0; d < D; ++d) {
+        out_shape[d + 1] = shape[d];
+    }
+    auto result = make_int64_array(out_shape);
+
+    std::vector<std::ptrdiff_t> view_shape(out_shape.size());
+    for (std::size_t i = 0; i < out_shape.size(); ++i) {
+        view_shape[i] = static_cast<std::ptrdiff_t>(out_shape[i]);
+    }
+    ArrayView<std::int64_t> view{result.data(), view_shape, {}};
+
+    std::uint64_t n_valid = 0;
+    if (strides.has_value()) {
+        if (strides->size() != D) {
+            throw std::invalid_argument(
+                "strides must have length matching graph ndim"
+            );
+        }
+        std::array<std::ptrdiff_t, D> stride_arr{};
+        for (std::size_t d = 0; d < D; ++d) {
+            stride_arr[d] = (*strides)[d];
+        }
+        nb::gil_scoped_release release;
+        n_valid = graph::project_edge_ids_to_pixels_with_offsets<D>(
+            graph, off_arr, stride_arr, view
+        );
+    } else if (mask.has_value()) {
+        const auto &m = *mask;
+        if (m.ndim() != D + 1 || m.shape(0) != offsets.size()) {
+            throw std::invalid_argument(
+                "mask shape must be (n_offsets, *graph.shape)"
+            );
+        }
+        for (std::size_t d = 0; d < D; ++d) {
+            if (m.shape(d + 1) != shape[d]) {
+                throw std::invalid_argument(
+                    "mask shape must be (n_offsets, *graph.shape)"
+                );
+            }
+        }
+        std::vector<std::ptrdiff_t> mask_shape(m.ndim());
+        for (std::size_t i = 0; i < m.ndim(); ++i) {
+            mask_shape[i] = static_cast<std::ptrdiff_t>(m.shape(i));
+        }
+        ConstArrayView<std::uint8_t> mview{m.data(), mask_shape, {}};
+        nb::gil_scoped_release release;
+        n_valid = graph::project_edge_ids_to_pixels_with_offsets<D>(
+            graph, off_arr, mview, view
+        );
+    } else {
+        nb::gil_scoped_release release;
+        n_valid = graph::project_edge_ids_to_pixels_with_offsets<D>(
+            graph, off_arr, view
+        );
+    }
+    return nb::make_tuple(result, n_valid);
+}
+
 template <class T, std::size_t D>
 FloatingArray<T> grid_boundary_features_t(
     const graph::GridGraph<D> &graph,
@@ -1654,6 +1769,14 @@ void bind_graph(nb::module_ &m) {
         .def("edge_axis", &GridGraph2D::edge_axis, nb::arg("edge"))
         .def("edge_coordinates", &grid_edge_coordinates<2>, nb::arg("edge"))
         .def("offset_target", &grid_offset_target<2>, nb::arg("node"), nb::arg("offset"))
+        .def("project_edge_ids_to_pixels", &grid_project_edge_ids_to_pixels<2>)
+        .def(
+            "project_edge_ids_to_pixels_with_offsets",
+            &grid_project_edge_ids_to_pixels_with_offsets<2>,
+            nb::arg("offsets"),
+            nb::arg("strides") = std::nullopt,
+            nb::arg("mask") = std::nullopt
+        )
         .def_prop_ro("numberOfDimensions", &GridGraph2D::ndim)
         .def("nodeId", &grid_node_id<2>, nb::arg("coordinate"))
         .def("edgeAxis", &GridGraph2D::edge_axis, nb::arg("edge"))
@@ -1670,6 +1793,14 @@ void bind_graph(nb::module_ &m) {
         .def("edge_axis", &GridGraph3D::edge_axis, nb::arg("edge"))
         .def("edge_coordinates", &grid_edge_coordinates<3>, nb::arg("edge"))
         .def("offset_target", &grid_offset_target<3>, nb::arg("node"), nb::arg("offset"))
+        .def("project_edge_ids_to_pixels", &grid_project_edge_ids_to_pixels<3>)
+        .def(
+            "project_edge_ids_to_pixels_with_offsets",
+            &grid_project_edge_ids_to_pixels_with_offsets<3>,
+            nb::arg("offsets"),
+            nb::arg("strides") = std::nullopt,
+            nb::arg("mask") = std::nullopt
+        )
         .def_prop_ro("numberOfDimensions", &GridGraph3D::ndim)
         .def("nodeId", &grid_node_id<3>, nb::arg("coordinate"))
         .def("edgeAxis", &GridGraph3D::edge_axis, nb::arg("edge"))
diff --git a/src/bioimage_cpp/graph/__init__.py b/src/bioimage_cpp/graph/__init__.py
index 65e4f33..cae232e 100644
--- a/src/bioimage_cpp/graph/__init__.py
+++ b/src/bioimage_cpp/graph/__init__.py
@@ -106,6 +106,39 @@ def deserialize(cls, serialization):
         return cls.from_edges(number_of_nodes, uvs)
 
 
+def _normalize_projection_offsets(offsets, ndim: int) -> list[list[int]]:
+    normalized: list[list[int]] = []
+    for index, offset in enumerate(offsets):
+        values = [int(v) for v in offset]
+        if len(values) != ndim:
+            raise ValueError(
+                f"offsets[{index}] must have length {ndim}, got length={len(values)}"
+            )
+        normalized.append(values)
+    return normalized
+
+
+def _normalize_projection_strides(strides, ndim: int) -> list[int]:
+    values = [int(v) for v in strides]
+    if len(values) != ndim:
+        raise ValueError(
+            f"strides must have length {ndim}, got length={len(values)}"
+        )
+    if any(v <= 0 for v in values):
+        raise ValueError("strides must be positive")
+    return values
+
+
+def _normalize_projection_mask(mask, n_offsets: int, shape: tuple[int, ...]) -> np.ndarray:
+    array = np.asarray(mask)
+    expected = (n_offsets, *shape)
+    if array.shape != expected:
+        raise ValueError(
+            f"mask shape must be {expected}, got shape={array.shape}"
+        )
+    return np.ascontiguousarray(array.astype(np.uint8, copy=False))
+
+
 class GridGraph2D(_core.GridGraph2D):
     """Regular 2D nearest-neighbor grid graph.
 
@@ -128,6 +161,54 @@ def offset_target(self, node: int, offset):
     def offsetTarget(self, node: int, offset):
         return self.offset_target(node, offset)
 
+    def project_edge_ids_to_pixels(self) -> np.ndarray:
+        """Project edge ids onto a per-pixel array of shape ``(2, *shape)``.
+
+        Returns an ``int64`` array initialised to ``-1``. For every grid edge
+        with spanning axis ``d`` and pivot (smaller-endpoint) coordinate
+        ``c``, sets ``out[d, *c] = edge_id``. Slots where the pivot lies at
+        the last index of axis ``d`` remain ``-1``.
+        """
+        return super().project_edge_ids_to_pixels()
+
+    def project_edge_ids_to_pixels_with_offsets(
+        self,
+        offsets,
+        *,
+        strides=None,
+        mask=None,
+    ) -> tuple[np.ndarray, int]:
+        """Enumerate lifted edges defined by per-channel ``offsets``.
+
+        Walks ``(offset_idx, *coord)`` in C-order over
+        ``(len(offsets), *shape)``. For every coord whose target
+        ``coord + offsets[offset_idx]`` is in bounds (and survives the
+        optional ``strides`` or ``mask`` filter), writes a sequential
+        counter starting at 0; rejected slots get ``-1``.
+
+        ``strides`` and ``mask`` are mutually exclusive. ``strides`` keeps
+        only coords where every ``coord[d] % strides[d] == 0``; ``mask``
+        keeps only coords where ``mask[offset_idx, *coord]`` is true.
+
+        Returns ``(array, n_valid)`` — the ``int64`` array and the total
+        number of valid entries written. The counter is **not** a graph
+        edge id; it indexes into the implicit array of lifted edges.
+        """
+        if strides is not None and mask is not None:
+            raise ValueError("strides and mask cannot be given together")
+        offsets = _normalize_projection_offsets(offsets, 2)
+        normalized_strides = (
+            None if strides is None else _normalize_projection_strides(strides, 2)
+        )
+        normalized_mask = (
+            None
+            if mask is None
+            else _normalize_projection_mask(mask, len(offsets), tuple(self.shape))
+        )
+        return super().project_edge_ids_to_pixels_with_offsets(
+            offsets, normalized_strides, normalized_mask
+        )
+
 
 class GridGraph3D(_core.GridGraph3D):
     """Regular 3D nearest-neighbor grid graph.
@@ -151,6 +232,42 @@ def offset_target(self, node: int, offset):
     def offsetTarget(self, node: int, offset):
         return self.offset_target(node, offset)
 
+    def project_edge_ids_to_pixels(self) -> np.ndarray:
+        """Project edge ids onto a per-pixel array of shape ``(3, *shape)``.
+
+        See :meth:`GridGraph2D.project_edge_ids_to_pixels` for the
+        semantics; the only difference is the leading dimension matches the
+        graph rank.
+        """
+        return super().project_edge_ids_to_pixels()
+
+    def project_edge_ids_to_pixels_with_offsets(
+        self,
+        offsets,
+        *,
+        strides=None,
+        mask=None,
+    ) -> tuple[np.ndarray, int]:
+        """Enumerate lifted edges defined by per-channel ``offsets``.
+
+        See :meth:`GridGraph2D.project_edge_ids_to_pixels_with_offsets` for
+        the semantics; the only difference is the input shape.
+        """
+        if strides is not None and mask is not None:
+            raise ValueError("strides and mask cannot be given together")
+        offsets = _normalize_projection_offsets(offsets, 3)
+        normalized_strides = (
+            None if strides is None else _normalize_projection_strides(strides, 3)
+        )
+        normalized_mask = (
+            None
+            if mask is None
+            else _normalize_projection_mask(mask, len(offsets), tuple(self.shape))
+        )
+        return super().project_edge_ids_to_pixels_with_offsets(
+            offsets, normalized_strides, normalized_mask
+        )
+
 
 RegionAdjacencyGraph = _core.RegionAdjacencyGraph
 
diff --git a/tests/graph/test_grid_graph.py b/tests/graph/test_grid_graph.py
index 8f25994..54e19ad 100644
--- a/tests/graph/test_grid_graph.py
+++ b/tests/graph/test_grid_graph.py
@@ -172,6 +172,139 @@ def test_grid_graph_freeze_is_idempotent_and_safe_to_call_repeatedly():
     assert _adjacency_pairs(graph.node_adjacency(4)) == {(1, 1), (7, 4), (3, 8), (5, 9)}
 
 
+def _bounds_ok(coord, offset, shape):
+    return all(0 <= coord[d] + offset[d] < shape[d] for d in range(len(shape)))
+
+
+def test_project_edge_ids_to_pixels_2d_writes_each_edge_at_its_pivot():
+    graph = bic.graph.GridGraph2D((2, 3))
+    out = graph.project_edge_ids_to_pixels()
+
+    assert out.shape == (2, 2, 3)
+    assert out.dtype == np.int64
+    # Exactly graph.number_of_edges valid entries; the rest stay -1.
+    assert int((out != -1).sum()) == graph.number_of_edges
+    # Every edge id appears exactly once, at the slot
+    # (edge_axis, *edge_coordinates(edge)).
+    for edge in range(graph.number_of_edges):
+        axis = graph.edge_axis(edge)
+        pivot = tuple(int(c) for c in graph.edge_coordinates(edge))
+        assert out[(axis,) + pivot] == edge
+    # Slots on the "missing" row/column for each spanning axis stay -1.
+    np.testing.assert_array_equal(out[0, 1, :], -1)
+    np.testing.assert_array_equal(out[1, :, 2], -1)
+
+
+def test_project_edge_ids_to_pixels_3d_round_trip_through_edge_coordinates():
+    graph = bic.graph.GridGraph3D((2, 2, 2))
+    out = graph.project_edge_ids_to_pixels()
+
+    assert out.shape == (3, 2, 2, 2)
+    assert int((out != -1).sum()) == graph.number_of_edges
+    for edge in range(graph.number_of_edges):
+        axis = graph.edge_axis(edge)
+        pivot = tuple(int(c) for c in graph.edge_coordinates(edge))
+        assert out[(axis,) + pivot] == edge
+
+
+def test_project_edge_ids_to_pixels_with_offsets_2d_numbers_in_bounds_targets():
+    graph = bic.graph.GridGraph2D((2, 3))
+    offsets = [(1, 0), (0, 1), (2, 0), (0, -1)]
+    out, n_valid = graph.project_edge_ids_to_pixels_with_offsets(offsets)
+
+    assert out.shape == (4, 2, 3)
+    assert out.dtype == np.int64
+
+    # Walk (off_idx, *coord) in C-order, predict counter, compare to output.
+    expected = np.full((4, 2, 3), -1, dtype=np.int64)
+    counter = 0
+    for off_idx, offset in enumerate(offsets):
+        for r in range(2):
+            for c in range(3):
+                if _bounds_ok((r, c), offset, (2, 3)):
+                    expected[off_idx, r, c] = counter
+                    counter += 1
+    np.testing.assert_array_equal(out, expected)
+    assert n_valid == counter
+    # Offset (2, 0) on a height-2 grid is never in bounds.
+    np.testing.assert_array_equal(out[2], -1)
+
+
+def test_project_edge_ids_to_pixels_with_offsets_3d_long_range_offset_partially_out_of_bounds():
+    graph = bic.graph.GridGraph3D((3, 2, 2))
+    offsets = [(2, 0, 0), (0, 1, 0)]
+    out, n_valid = graph.project_edge_ids_to_pixels_with_offsets(offsets)
+
+    assert out.shape == (2, 3, 2, 2)
+    # Offset (2, 0, 0): valid only when z + 2 < 3, i.e. z == 0 — that's 4 entries.
+    # Offset (0, 1, 0): valid when y + 1 < 2, i.e. y == 0 — that's 3 * 1 * 2 = 6 entries.
+    assert n_valid == 4 + 6
+    # Sequential numbering: entries in C-order are 0..n_valid-1, rest are -1.
+    valid = out[out != -1]
+    np.testing.assert_array_equal(valid, np.arange(n_valid, dtype=np.int64))
+
+
+def test_project_edge_ids_to_pixels_with_offsets_strides_filter_keeps_aligned_coords():
+    graph = bic.graph.GridGraph2D((4, 4))
+    offsets = [(1, 0), (0, 1)]
+    out, n_valid = graph.project_edge_ids_to_pixels_with_offsets(
+        offsets, strides=(2, 2)
+    )
+
+    # Strides (2, 2) keep only (r, c) with r%2 == c%2 == 0: (0,0), (0,2),
+    # (2,0), (2,2). All four are in bounds for both offsets — 8 entries.
+    assert n_valid == 8
+    valid = out[out != -1]
+    np.testing.assert_array_equal(valid, np.arange(8, dtype=np.int64))
+    # Non-aligned coords are -1.
+    for r in range(4):
+        for c in range(4):
+            if r % 2 != 0 or c % 2 != 0:
+                assert out[0, r, c] == -1
+                assert out[1, r, c] == -1
+
+
+def test_project_edge_ids_to_pixels_with_offsets_mask_filter_keeps_true_positions():
+    graph = bic.graph.GridGraph2D((2, 3))
+    offsets = [(1, 0), (0, 1)]
+    mask = np.zeros((2, 2, 3), dtype=bool)
+    mask[0, 0, 1] = True  # (off=0, coord=(0,1)) — target (1,1) in bounds
+    mask[1, 1, 0] = True  # (off=1, coord=(1,0)) — target (1,1) in bounds
+    mask[1, 0, 2] = True  # (off=1, coord=(0,2)) — target (0,3) OUT OF BOUNDS
+    out, n_valid = graph.project_edge_ids_to_pixels_with_offsets(
+        offsets, mask=mask
+    )
+
+    # The out-of-bounds entry is suppressed even though the mask says true.
+    assert n_valid == 2
+    expected = np.full((2, 2, 3), -1, dtype=np.int64)
+    expected[0, 0, 1] = 0
+    expected[1, 1, 0] = 1
+    np.testing.assert_array_equal(out, expected)
+
+
+def test_project_edge_ids_to_pixels_with_offsets_rejects_strides_and_mask_together():
+    graph = bic.graph.GridGraph2D((2, 3))
+    mask = np.zeros((1, 2, 3), dtype=bool)
+    with pytest.raises(ValueError, match="strides and mask"):
+        graph.project_edge_ids_to_pixels_with_offsets(
+            [(1, 0)], strides=(1, 1), mask=mask
+        )
+
+
+def test_project_edge_ids_to_pixels_with_offsets_validates_shapes():
+    graph = bic.graph.GridGraph2D((2, 3))
+    with pytest.raises(ValueError, match=r"offsets\[0\] must have length 2"):
+        graph.project_edge_ids_to_pixels_with_offsets([(1, 0, 0)])
+    with pytest.raises(ValueError, match="strides must have length 2"):
+        graph.project_edge_ids_to_pixels_with_offsets([(1, 0)], strides=(1,))
+    with pytest.raises(ValueError, match="strides must be positive"):
+        graph.project_edge_ids_to_pixels_with_offsets([(1, 0)], strides=(0, 1))
+    with pytest.raises(ValueError, match="mask shape must be"):
+        bad_mask = np.zeros((1, 2, 4), dtype=bool)
+        graph.project_edge_ids_to_pixels_with_offsets([(1, 0)], mask=bad_mask)
+
+
 def test_grid_graph_rejects_invalid_shapes_and_coordinates():
     with pytest.raises(ValueError, match="length 2"):
         bic.graph.GridGraph2D((2, 3, 4))