Implement RLE and bump version (#52)

Author: constantinpape

MIGRATION_GUIDE.md

@@ -770,6 +770,31 @@ Notes:
 - Every value in the input must be present in the mapping.
 - Non-contiguous inputs are copied before entering C++.
 
+### Run-Length Encoding
+
+`nifty.tools.computeRLE` computes the COCO-style binary run-length encoding of an
+array. Use `bic.utils.compute_rle`:
+
+```python
+import nifty.tools as nt
+
+mask = np.array([[0, 0, 1], [1, 1, 0], [0, 1, 1]], dtype=np.uint8)
+
+rle = nt.computeRLE(mask)            # nifty: Python list [2, 3, 2, 2]
+rle = bic.utils.compute_rle(mask)    # bioimage-cpp: np.int64 array([2, 3, 2, 2])
+```
+
+Notes:
+
+- The array is flattened in **C-order** and interpreted as binary
+  (zero vs. nonzero). Counts always start with a run of zeros and then
+  alternate; a leading `0` is emitted when the first element is nonzero.
+- Supported input dtypes are `bool`, `uint8`, `uint16`, `uint32`, `uint64`,
+  `int32`, and `int64`.
+- **Behavioral difference:** nifty returns a Python `list`; `bioimage-cpp`
+  returns a 1-D NumPy `int64` array.
+- Non-contiguous inputs are copied before entering C++.
+
 ### Edge-Weighted Watershed
 
 Nifty:

include/bioimage_cpp/run_length.hxx

@@ -0,0 +1,48 @@
+#pragma once
+
+#include "bioimage_cpp/array_view.hxx"
+
+#include <cstdint>
+#include <numeric>
+#include <vector>
+
+namespace bioimage_cpp {
+
+// Compute the COCO-style binary run-length encoding of `mask`.
+//
+// The array is treated as a flat C-order buffer and interpreted as binary
+// (zero vs. nonzero). The returned run lengths always start with a run of
+// zeros and then alternate (zeros, ones, zeros, ...). If the first element is
+// nonzero a leading 0 is emitted. An empty input yields an empty result.
+template <class T>
+std::vector<std::int64_t> compute_rle(const ConstArrayView<T> &mask) {
+    const auto n = std::accumulate(
+        mask.shape.begin(),
+        mask.shape.end(),
+        std::ptrdiff_t{1},
+        [](const std::ptrdiff_t a, const std::ptrdiff_t b) { return a * b; }
+    );
+
+    std::vector<std::int64_t> counts;
+    if (n <= 0) {
+        return counts;
+    }
+
+    bool current_value = false;  // runs start with zeros
+    std::int64_t run_length = 0;
+    for (std::ptrdiff_t i = 0; i < n; ++i) {
+        const bool value = mask.data[i] != T{0};
+        if (value == current_value) {
+            ++run_length;
+        } else {
+            counts.push_back(run_length);
+            current_value = value;
+            run_length = 1;
+        }
+    }
+    counts.push_back(run_length);
+
+    return counts;
+}
+
+} // namespace bioimage_cpp

src/bindings/utils.cxx

@@ -2,12 +2,14 @@
 
 #include "bioimage_cpp/array_view.hxx"
 #include "bioimage_cpp/mesh_smoothing.hxx"
+#include "bioimage_cpp/run_length.hxx"
 #include "bioimage_cpp/take_dict.hxx"
 
 #include <nanobind/ndarray.h>
 #include <nanobind/stl/pair.h>
 #include <nanobind/stl/string.h>
 
+#include <algorithm>
 #include <cstddef>
 #include <cstdint>
 #include <memory>
@@ -151,6 +153,35 @@ std::pair<OutputArray<V>, OutputArray<V>> smooth_mesh_t(
     return {std::move(out_verts), std::move(out_normals)};
 }
 
+template <class T>
+OutputArray<std::int64_t> compute_rle_t(InputArray<T> mask) {
+    std::vector<std::ptrdiff_t> view_shape(mask.ndim());
+    for (std::size_t axis = 0; axis < mask.ndim(); ++axis) {
+        view_shape[axis] = static_cast<std::ptrdiff_t>(mask.shape(axis));
+    }
+
+    ConstArrayView<T> input{
+        mask.data(),
+        view_shape,
+        {},
+    };
+
+    std::vector<std::int64_t> counts;
+    {
+        nb::gil_scoped_release release;
+        counts = compute_rle<T>(input);
+    }
+
+    auto *data = new std::int64_t[counts.size()];
+    std::copy(counts.begin(), counts.end(), data);
+    nb::capsule owner(data, [](void *p) noexcept {
+        delete[] static_cast<std::int64_t *>(p);
+    });
+
+    std::size_t shape[1] = {counts.size()};
+    return OutputArray<std::int64_t>(data, 1, shape, owner);
+}
+
 } // namespace
 
 void bind_utils(nb::module_ &m) {
@@ -202,6 +233,20 @@ void bind_utils(nb::module_ &m) {
         nb::arg("n_threads"),
         "Laplacian smoothing of a triangular mesh with float64 vertices and normals."
     );
+    m.def("_compute_rle_bool", &compute_rle_t<bool>, nb::arg("mask"),
+          "COCO-style binary run-length encoding of a contiguous bool array.");
+    m.def("_compute_rle_uint8", &compute_rle_t<std::uint8_t>, nb::arg("mask"),
+          "COCO-style binary run-length encoding of a contiguous uint8 array.");
+    m.def("_compute_rle_uint16", &compute_rle_t<std::uint16_t>, nb::arg("mask"),
+          "COCO-style binary run-length encoding of a contiguous uint16 array.");
+    m.def("_compute_rle_uint32", &compute_rle_t<std::uint32_t>, nb::arg("mask"),
+          "COCO-style binary run-length encoding of a contiguous uint32 array.");
+    m.def("_compute_rle_uint64", &compute_rle_t<std::uint64_t>, nb::arg("mask"),
+          "COCO-style binary run-length encoding of a contiguous uint64 array.");
+    m.def("_compute_rle_int32", &compute_rle_t<std::int32_t>, nb::arg("mask"),
+          "COCO-style binary run-length encoding of a contiguous int32 array.");
+    m.def("_compute_rle_int64", &compute_rle_t<std::int64_t>, nb::arg("mask"),
+          "COCO-style binary run-length encoding of a contiguous int64 array.");
 }
 
 } // namespace bioimage_cpp::bindings

src/bioimage_cpp/_version.py

@@ -1 +1 @@
-__version__ = "0.2.0"
+__version__ = "0.3.0"

src/bioimage_cpp/utils.py

@@ -24,6 +24,16 @@
     np.dtype("float64"): _core._smooth_mesh_float64,
 }
 
+_COMPUTE_RLE_BY_DTYPE = {
+    np.dtype("bool"): _core._compute_rle_bool,
+    np.dtype("uint8"): _core._compute_rle_uint8,
+    np.dtype("uint16"): _core._compute_rle_uint16,
+    np.dtype("uint32"): _core._compute_rle_uint32,
+    np.dtype("uint64"): _core._compute_rle_uint64,
+    np.dtype("int32"): _core._compute_rle_int32,
+    np.dtype("int64"): _core._compute_rle_int64,
+}
+
 _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)]
@@ -307,6 +317,47 @@ def take_dict(relabeling: Mapping[Any, Any], to_relabel: np.ndarray) -> np.ndarr
         raise
 
 
+def compute_rle(mask: np.ndarray) -> np.ndarray:
+    """Compute the COCO-style binary run-length encoding of an array.
+
+    The array is flattened in C-order and interpreted as binary (zero vs.
+    nonzero). The returned run lengths always start with a run of zeros and then
+    alternate (zeros, ones, zeros, ...); if the first element is nonzero a
+    leading ``0`` is emitted. An empty array yields an empty result.
+
+    This is the equivalent of ``nifty.tools.computeRLE``. Unlike nifty, which
+    returns a Python list, this returns a NumPy ``int64`` array.
+
+    Parameters
+    ----------
+    mask:
+        NumPy array with dtype ``bool``, ``uint8``, ``uint16``, ``uint32``,
+        ``uint64``, ``int32``, or ``int64``. Non-contiguous inputs are copied to
+        a contiguous array before calling the C++ kernel.
+
+    Returns
+    -------
+    np.ndarray
+        1-D ``int64`` array of run lengths.
+
+    Raises
+    ------
+    TypeError
+        If ``mask`` has an unsupported dtype.
+    """
+    array = np.asarray(mask)
+    dtype = array.dtype
+    try:
+        run = _COMPUTE_RLE_BY_DTYPE[dtype]
+    except KeyError as error:
+        supported = ", ".join(str(dtype) for dtype in _COMPUTE_RLE_BY_DTYPE)
+        raise TypeError(
+            f"mask must have one of dtypes ({supported}), got dtype={dtype}"
+        ) from error
+
+    return run(np.ascontiguousarray(array))
+
+
 def smooth_mesh(
     verts: np.ndarray,
     normals: np.ndarray,
@@ -473,6 +524,7 @@ def _validate_normalize_by(normalize_by: str) -> None:
     "BlockWithHalo",
     "SegmentationOverlap",
     "UnionFind",
+    "compute_rle",
     "segmentation_overlap",
     "smooth_mesh",
     "take_dict",

tests/test_utils_compute_rle.py

@@ -0,0 +1,101 @@
+import numpy as np
+import pytest
+
+import bioimage_cpp as bic
+
+SUPPORTED_DTYPES = [
+    np.bool_,
+    np.uint8,
+    np.uint16,
+    np.uint32,
+    np.uint64,
+    np.int32,
+    np.int64,
+]
+
+
+def _rle(values, dtype=np.uint8):
+    return bic.utils.compute_rle(np.asarray(values, dtype=dtype))
+
+
+@pytest.mark.parametrize(
+    "values, expected",
+    [
+        ([0, 0, 1, 1, 1, 0, 0, 1, 1], [2, 3, 2, 2]),
+        ([1, 1, 0, 1], [0, 2, 1, 1]),
+        ([1, 1, 1, 1, 1], [0, 5]),
+        ([0, 0, 0, 0], [4]),
+        ([1], [0, 1]),
+        ([0], [1]),
+    ],
+)
+def test_known_cases(values, expected):
+    result = _rle(values)
+    np.testing.assert_array_equal(result, np.asarray(expected, dtype=np.int64))
+
+
+def test_empty_input():
+    result = _rle([])
+    assert result.dtype == np.int64
+    assert result.shape == (0,)
+
+
+def test_output_is_int64_1d():
+    result = _rle([0, 1, 1])
+    assert result.dtype == np.int64
+    assert result.ndim == 1
+
+
+@pytest.mark.parametrize("dtype", SUPPORTED_DTYPES)
+def test_dtype_variants(dtype):
+    mask = np.array([0, 1, 1, 0, 0, 1], dtype=dtype)
+    result = bic.utils.compute_rle(mask)
+    np.testing.assert_array_equal(result, np.array([1, 2, 2, 1], dtype=np.int64))
+
+
+def _reference_rle(flat):
+    # COCO-style: counts of alternating runs starting with zeros.
+    counts = []
+    current = 0
+    run = 0
+    for value in flat:
+        binary = int(value != 0)
+        if binary == current:
+            run += 1
+        else:
+            counts.append(run)
+            current = binary
+            run = 1
+    counts.append(run)
+    return np.asarray(counts, dtype=np.int64)
+
+
+@pytest.mark.parametrize("shape", [(6, 8), (3, 4, 5)])
+def test_c_order_flatten(shape):
+    rng = np.random.default_rng(0)
+    mask = (rng.integers(0, 2, size=shape)).astype(np.uint8)
+    result = bic.utils.compute_rle(mask)
+    expected = _reference_rle(mask.ravel(order="C"))
+    np.testing.assert_array_equal(result, expected)
+
+
+def test_non_contiguous_input():
+    rng = np.random.default_rng(1)
+    base = (rng.integers(0, 2, size=(5, 6))).astype(np.uint8)
+    view = base.T  # non-contiguous (Fortran-ordered view)
+    assert not view.flags["C_CONTIGUOUS"]
+    result = bic.utils.compute_rle(view)
+    expected = _reference_rle(np.ascontiguousarray(view).ravel(order="C"))
+    np.testing.assert_array_equal(result, expected)
+
+
+def test_nonzero_values_are_binary():
+    # Values other than 0/1 are treated as "set".
+    mask = np.array([0, 7, 3, 0], dtype=np.int32)
+    result = bic.utils.compute_rle(mask)
+    np.testing.assert_array_equal(result, np.array([1, 2, 1], dtype=np.int64))
+
+
+def test_invalid_dtype_raises():
+    with pytest.raises(TypeError):
+        bic.utils.compute_rle(np.array([0.0, 1.0], dtype=np.float32))