Refactor to avoid duplicate functionality

Author: constantinpape

include/bioimage_cpp/detail/edge_hash.hxx

@@ -0,0 +1,27 @@
+#pragma once
+
+#include <cstddef>
+#include <cstdint>
+#include <utility>
+
+namespace bioimage_cpp::detail {
+
+using NodeId = std::uint64_t;
+using Edge = std::pair<NodeId, NodeId>;
+
+inline Edge edge_key(NodeId u, NodeId v) {
+    if (v < u) {
+        std::swap(u, v);
+    }
+    return {u, v};
+}
+
+struct EdgeHash {
+    std::size_t operator()(const Edge &edge) const {
+        const auto first = static_cast<std::size_t>(edge.first);
+        const auto second = static_cast<std::size_t>(edge.second);
+        return first ^ (second + 0x9e3779b97f4a7c15ULL + (first << 6U) + (first >> 2U));
+    }
+};
+
+} // namespace bioimage_cpp::detail

include/bioimage_cpp/detail/grid.hxx

@@ -0,0 +1,59 @@
+#pragma once
+
+#include <cstddef>
+#include <cstdint>
+#include <vector>
+
+namespace bioimage_cpp::detail {
+
+// C-order strides for a row-major array of the given shape, in units of array
+// elements (not bytes). The innermost (last) axis has stride 1.
+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;
+}
+
+// Translate a flat node index by a per-axis offset on a row-major grid.
+//
+// Returns true when the offset keeps the result inside the grid, in which case
+// `target_out` is set to the neighbor's flat index. Returns false otherwise and
+// leaves `target_out` unchanged. `strides` must match `shape` and is typically
+// `c_order_strides(shape)`.
+inline bool valid_offset_target(
+    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,
+    std::uint64_t &target_out
+) {
+    std::int64_t target_signed = static_cast<std::int64_t>(node);
+    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;
+        }
+        target_signed += static_cast<std::int64_t>(offset[axis] * strides[axis]);
+    }
+    target_out = static_cast<std::uint64_t>(target_signed);
+    return true;
+}
+
+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
+) {
+    std::uint64_t unused = 0;
+    return valid_offset_target(node, offset, shape, strides, unused);
+}
+
+} // namespace bioimage_cpp::detail

include/bioimage_cpp/detail/relabel.hxx

@@ -0,0 +1,25 @@
+#pragma once
+
+#include <cstddef>
+#include <cstdint>
+#include <unordered_map>
+#include <vector>
+
+namespace bioimage_cpp::detail {
+
+// Map an arbitrary labeling to the dense range [0, k) preserving the first
+// occurrence order of each distinct input label.
+inline std::vector<std::uint64_t> dense_relabel(const std::vector<std::uint64_t> &labels) {
+    std::unordered_map<std::uint64_t, std::uint64_t> relabeling;
+    std::vector<std::uint64_t> result(labels.size());
+    for (std::size_t index = 0; index < labels.size(); ++index) {
+        auto found = relabeling.find(labels[index]);
+        if (found == relabeling.end()) {
+            found = relabeling.emplace(labels[index], static_cast<std::uint64_t>(relabeling.size())).first;
+        }
+        result[index] = found->second;
+    }
+    return result;
+}
+
+} // namespace bioimage_cpp::detail

include/bioimage_cpp/detail/threading.hxx

@@ -0,0 +1,60 @@
+#pragma once
+
+#include <algorithm>
+#include <cstddef>
+#include <thread>
+#include <utility>
+#include <vector>
+
+namespace bioimage_cpp::detail {
+
+inline std::size_t normalize_thread_count(
+    const std::size_t requested,
+    const std::size_t number_of_work_items
+) {
+    if (number_of_work_items == 0) {
+        return 1;
+    }
+    std::size_t n_threads = requested;
+    if (n_threads == 0) {
+        n_threads = std::thread::hardware_concurrency();
+        if (n_threads == 0) {
+            n_threads = 1;
+        }
+    }
+    return std::max<std::size_t>(1, std::min(n_threads, number_of_work_items));
+}
+
+// Split [0, number_of_work_items) into n_threads contiguous chunks and invoke
+// `chunk(thread_id, begin, end)` on each chunk. Thread 0 runs on the calling
+// thread; threads 1..n_threads-1 run in parallel std::threads. The caller is
+// responsible for picking n_threads via normalize_thread_count and for the
+// thread safety of `chunk`.
+template <class Chunk>
+void parallel_for_chunks(
+    const std::size_t n_threads,
+    const std::size_t number_of_work_items,
+    Chunk &&chunk
+) {
+    const auto bounds = [&](const std::size_t thread_id) {
+        const auto begin = thread_id * number_of_work_items / n_threads;
+        const auto end = (thread_id + 1) * number_of_work_items / n_threads;
+        return std::pair<std::size_t, std::size_t>{begin, end};
+    };
+
+    std::vector<std::thread> threads;
+    threads.reserve(n_threads > 0 ? n_threads - 1 : 0);
+    for (std::size_t thread_id = 1; thread_id < n_threads; ++thread_id) {
+        const auto [begin, end] = bounds(thread_id);
+        threads.emplace_back([thread_id, begin, end, &chunk]() {
+            chunk(thread_id, begin, end);
+        });
+    }
+    const auto [begin, end] = bounds(0);
+    chunk(0, begin, end);
+    for (auto &thread : threads) {
+        thread.join();
+    }
+}
+
+} // namespace bioimage_cpp::detail

include/bioimage_cpp/detail/union_find.hxx

@@ -0,0 +1,78 @@
+#pragma once
+
+#include <cstddef>
+#include <cstdint>
+#include <numeric>
+#include <utility>
+#include <vector>
+
+namespace bioimage_cpp::detail {
+
+// Disjoint-set / union-find with path compression and union-by-rank.
+//
+// Three merge variants are provided so both rank-driven and caller-driven
+// merge directions are expressible:
+//   * merge(u, v):        find both, then union by rank.
+//   * merge_to(stable, removed): find both, then force `stable` to be the
+//                                parent regardless of rank. The caller picks
+//                                which root survives.
+//   * unite_roots(a, b):  preconditions: a and b are roots and a != b. Union
+//                         by rank without an additional find().
+class UnionFind {
+public:
+    explicit UnionFind(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) {
+        auto current = static_cast<std::size_t>(node);
+        while (parents_[current] != current) {
+            parents_[current] = parents_[parents_[current]];
+            current = static_cast<std::size_t>(parents_[current]);
+        }
+        return static_cast<std::uint64_t>(current);
+    }
+
+    std::uint64_t merge(std::uint64_t first, std::uint64_t second) {
+        first = find(first);
+        second = find(second);
+        if (first == second) {
+            return first;
+        }
+        return unite_roots(first, second);
+    }
+
+    std::uint64_t merge_to(std::uint64_t stable, std::uint64_t removed) {
+        stable = find(stable);
+        removed = find(removed);
+        if (stable == removed) {
+            return stable;
+        }
+        parents_[static_cast<std::size_t>(removed)] = stable;
+        if (ranks_[static_cast<std::size_t>(stable)] <= ranks_[static_cast<std::size_t>(removed)]) {
+            ranks_[static_cast<std::size_t>(stable)] = ranks_[static_cast<std::size_t>(removed)] + 1;
+        }
+        return stable;
+    }
+
+    std::uint64_t unite_roots(std::uint64_t first, std::uint64_t second) {
+        if (ranks_[static_cast<std::size_t>(first)] < ranks_[static_cast<std::size_t>(second)]) {
+            std::swap(first, second);
+        }
+        parents_[static_cast<std::size_t>(second)] = first;
+        if (ranks_[static_cast<std::size_t>(first)] == ranks_[static_cast<std::size_t>(second)]) {
+            ++ranks_[static_cast<std::size_t>(first)];
+        }
+        return first;
+    }
+
+    [[nodiscard]] std::size_t size() const {
+        return parents_.size();
+    }
+
+private:
+    std::vector<std::uint64_t> parents_;
+    std::vector<std::uint64_t> ranks_;
+};
+
+} // namespace bioimage_cpp::detail

include/bioimage_cpp/graph/connected_components.hxx

@@ -1,63 +1,17 @@
 #pragma once
 
+#include "bioimage_cpp/detail/union_find.hxx"
 #include "bioimage_cpp/graph/undirected_graph.hxx"
 
 #include <cstddef>
 #include <cstdint>
-#include <numeric>
 #include <unordered_map>
-#include <utility>
 #include <vector>
 
 namespace bioimage_cpp::graph {
 
-class UnionFind {
-public:
-    explicit UnionFind(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_[static_cast<std::size_t>(node)] != node) {
-            parents_[static_cast<std::size_t>(node)] = find(parents_[static_cast<std::size_t>(node)]);
-        }
-        return parents_[static_cast<std::size_t>(node)];
-    }
-
-    void merge(std::uint64_t first, std::uint64_t second) {
-        first = find(first);
-        second = find(second);
-        if (first == second) {
-            return;
-        }
-        if (ranks_[static_cast<std::size_t>(first)] < ranks_[static_cast<std::size_t>(second)]) {
-            std::swap(first, second);
-        }
-        parents_[static_cast<std::size_t>(second)] = first;
-        if (ranks_[static_cast<std::size_t>(first)] == ranks_[static_cast<std::size_t>(second)]) {
-            ++ranks_[static_cast<std::size_t>(first)];
-        }
-    }
-
-    void merge_to(std::uint64_t stable, std::uint64_t removed) {
-        stable = find(stable);
-        removed = find(removed);
-        if (stable == removed) {
-            return;
-        }
-        parents_[static_cast<std::size_t>(removed)] = stable;
-        if (ranks_[static_cast<std::size_t>(stable)] <= ranks_[static_cast<std::size_t>(removed)]) {
-            ranks_[static_cast<std::size_t>(stable)] = ranks_[static_cast<std::size_t>(removed)] + 1;
-        }
-    }
-
-private:
-    std::vector<std::uint64_t> parents_;
-    std::vector<std::uint64_t> ranks_;
-};
-
 inline std::vector<std::uint64_t> dense_labels_from_union_find(
-    UnionFind &sets,
+    detail::UnionFind &sets,
     const std::uint64_t number_of_nodes
 ) {
     std::unordered_map<std::uint64_t, std::uint64_t> relabeling;
@@ -77,7 +31,7 @@ inline std::vector<std::uint64_t> connected_components(
     const UndirectedGraph &graph,
     const std::uint8_t *edge_mask = nullptr
 ) {
-    UnionFind sets(static_cast<std::size_t>(graph.number_of_nodes()));
+    detail::UnionFind sets(static_cast<std::size_t>(graph.number_of_nodes()));
     for (std::uint64_t edge = 0; edge < graph.number_of_edges(); ++edge) {
         if (edge_mask != nullptr && edge_mask[static_cast<std::size_t>(edge)] == 0) {
             continue;