From 425778ca3d76ca63290446495b0778d9241525c0 Mon Sep 17 00:00:00 2001 From: Constantin Pape Date: Sat, 25 Jul 2026 20:57:26 +0200 Subject: [PATCH 1/2] Implement contraction graph and use it to simplify remove ticks --- CMakeLists.txt | 1 + MIGRATION_GUIDE.md | 79 +++ .../bioimage_cpp/graph/contraction_graph.hxx | 662 ++++++++++++++++++ .../graph/contraction_options.hxx | 10 + .../graph/detail/contraction_topology.hxx | 526 ++++++++++++++ src/bindings/graph.hxx | 1 + src/bindings/graph_contraction.cxx | 332 +++++++++ src/bindings/module.cxx | 1 + src/bioimage_cpp/graph/__init__.py | 261 ++++++- src/bioimage_cpp/skeleton/postprocessing.py | 146 ++-- tests/graph/test_contraction_graph.py | 298 ++++++++ tests/skeleton/test_postprocessing.py | 94 +++ 12 files changed, 2345 insertions(+), 66 deletions(-) create mode 100644 include/bioimage_cpp/graph/contraction_graph.hxx create mode 100644 include/bioimage_cpp/graph/contraction_options.hxx create mode 100644 include/bioimage_cpp/graph/detail/contraction_topology.hxx create mode 100644 src/bindings/graph_contraction.cxx create mode 100644 tests/graph/test_contraction_graph.py diff --git a/CMakeLists.txt b/CMakeLists.txt index 8ca97ec..ec87538 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -18,6 +18,7 @@ nanobind_add_module(_core src/bindings/filters.cxx src/bindings/flow.cxx src/bindings/graph.cxx + src/bindings/graph_contraction.cxx src/bindings/ground_truth.cxx src/bindings/label_multiset.cxx src/bindings/mesh.cxx diff --git a/MIGRATION_GUIDE.md b/MIGRATION_GUIDE.md index adf74b7..4bd3007 100644 --- a/MIGRATION_GUIDE.md +++ b/MIGRATION_GUIDE.md @@ -460,6 +460,69 @@ Common method/property mapping: | `extractSubgraphFromNodes` | `extract_subgraph_from_nodes` | | `edgesFromNodeList` | `edges_from_node_list` | +### Mutable Graph Contractions + +`ContractionGraph` supports custom graph postprocessing with stable ids and +named value maps: + +```python +import numpy as np +import bioimage_cpp as bic + +graph = bic.graph.UndirectedGraph.from_edges( + 3, + [[0, 1], [1, 2], [0, 2]], +) +work = bic.graph.ContractionGraph(graph) +work.add_node_values( + "position", + np.array([[0.0], [2.0], [10.0]], dtype=np.float32), + reduction="mean", +) +work.add_edge_values( + "weight", + np.array([100.0, 2.0, 3.0], dtype=np.float64), + reduction="sum", +) + +work.contract_edge(0, keep_node=0) +result = work.materialize() + +contracted = result.graph +positions = result.node_values["position"] +weights = result.edge_values["weight"] +``` + +Register all value maps before the first mutation. Maps can use `float32` or +`float64` independently. The supported reductions are `sum`, `mean`, `min`, +and `max`. + +Node values reduce when an edge contraction merges two nodes. Edge values +reduce when parallel edges merge. The contracted edge becomes internal and has +no output value. `suppress_node` removes a degree-2 node and reduces its two +incident edge values into the replacement edge. + +The input must be a simple `UndirectedGraph`. Use +`parallel_edges="keep"` to keep parallel edges that mutations create: + +```python +work = bic.graph.ContractionGraph(graph, parallel_edges="keep") +``` + +The default `parallel_edges="merge"` folds these edges immediately. In keep +mode, `find_edges(u, v)` returns all matching active edge ids. + +`materialize()` always returns a simple-graph snapshot. In keep mode, it folds +parallel edges in the snapshot with each edge map's reduction. It does not +change the mutable graph. Rows in `node_values` and `edge_values` match the ids +in the materialized graph. `node_mapping` and `edge_mapping` map input ids to +materialized ids. Multiple input edges can map to one output edge. The mappings +use `-1` when no output item represents an input item. + +The mutable object is not an `UndirectedGraph`. Active ids stay stable and can +have gaps after mutations. Materialization creates the dense +`UndirectedGraph`. + ### Region Adjacency Graphs Nifty: @@ -2606,6 +2669,22 @@ Important differences and current scope: heap. Compact IDs avoid full-volume shortest-path fields; the public Dijkstra functions remain dense and exact FP64. +### Skeleton graph postprocessing + +`remove_ticks(vertices, edges, tick_length, radii=None)` removes short terminal +branches from a skeleton graph. It compresses degree-2 paths, sums their +physical lengths, and repeatedly removes the shortest terminal branch below +the strict threshold. + +The mutable topology and length reduction run in the C++ contraction backend. +The Python wrapper selects branches and compacts the output arrays. Distinct +compressed paths between the same branch nodes remain separate, so cyclic and +theta-shaped skeletons keep their degree semantics. + +Standalone paths and components that contain only a cycle remain unchanged. +The result keeps the existing array contract: vertices are `float64`, edges are +`int64`, and optional radii preserve their dtype. + ### Blockwise skeletonization and exact stitching `bioimage_cpp.skeleton.distributed` provides the dependency-light primitives diff --git a/include/bioimage_cpp/graph/contraction_graph.hxx b/include/bioimage_cpp/graph/contraction_graph.hxx new file mode 100644 index 0000000..0440c35 --- /dev/null +++ b/include/bioimage_cpp/graph/contraction_graph.hxx @@ -0,0 +1,662 @@ +#pragma once + +#include "bioimage_cpp/graph/contraction_options.hxx" +#include "bioimage_cpp/graph/detail/contraction_topology.hxx" +#include "bioimage_cpp/graph/undirected_graph.hxx" +#include "bioimage_cpp/util/union_find.hxx" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace bioimage_cpp::graph { + +enum class Reduction { + Sum, + Mean, + Minimum, + Maximum, +}; + +template +struct MaterializedValues { + std::string name; + std::vector values; + std::vector shape; +}; + +using MaterializedValueMap = std::variant< + MaterializedValues, + MaterializedValues +>; + +struct MaterializedContraction { + UndirectedGraph graph; + std::vector node_values; + std::vector edge_values; + std::vector node_mapping; + std::vector edge_mapping; +}; + +template +class ReducedValueMap { +public: + ReducedValueMap( + std::string name, + const std::span values, + std::vector component_shape, + const std::size_t item_count, + const Reduction reduction + ) + : name_(std::move(name)), + component_shape_(std::move(component_shape)), + component_size_(shape_size(component_shape_)), + values_(values.begin(), values.end()), + counts_(reduction == Reduction::Mean ? item_count : 0, 1), + mean_sums_(reduction == Reduction::Mean ? values.size() : 0), + reduction_(reduction) { + if (values_.size() != checked_product(item_count, component_size_)) { + throw std::invalid_argument( + "values size does not match the item count and component shape" + ); + } + if (reduction_ == Reduction::Mean) { + std::transform( + values_.begin(), + values_.end(), + mean_sums_.begin(), + [](const T value) { + return static_cast(value); + } + ); + } + } + + [[nodiscard]] const std::string &name() const { + return name_; + } + + [[nodiscard]] const std::vector &component_shape() const { + return component_shape_; + } + + [[nodiscard]] std::size_t component_size() const { + return component_size_; + } + + [[nodiscard]] Reduction reduction() const { + return reduction_; + } + + void reduce(const std::uint64_t kept, const std::uint64_t removed) { + const auto kept_index = static_cast(kept); + const auto removed_index = static_cast(removed); + if (kept_index == removed_index) { + return; + } + const auto kept_offset = checked_product(kept_index, component_size_); + const auto removed_offset = checked_product(removed_index, component_size_); + + if (reduction_ == Reduction::Mean) { + const auto kept_count = counts_[kept_index]; + const auto removed_count = counts_[removed_index]; + const auto total = kept_count + removed_count; + const auto total_weight = static_cast(total); + for (std::size_t component = 0; component < component_size_; ++component) { + auto &sum = mean_sums_[kept_offset + component]; + sum += mean_sums_[removed_offset + component]; + values_[kept_offset + component] = + static_cast(sum / total_weight); + } + counts_[kept_index] = total; + return; + } + + for (std::size_t component = 0; component < component_size_; ++component) { + auto &a = values_[kept_offset + component]; + const auto b = values_[removed_offset + component]; + if (std::isnan(a) || std::isnan(b)) { + a = std::numeric_limits::quiet_NaN(); + } else if (reduction_ == Reduction::Sum) { + a += b; + } else if (reduction_ == Reduction::Minimum) { + a = std::min(a, b); + } else { + a = std::max(a, b); + } + } + } + + [[nodiscard]] MaterializedValues gather( + const std::span item_ids + ) const { + std::vector shape; + shape.reserve(component_shape_.size() + 1); + shape.push_back(item_ids.size()); + shape.insert(shape.end(), component_shape_.begin(), component_shape_.end()); + + std::vector output; + output.reserve(checked_product(item_ids.size(), component_size_)); + for (const auto item : item_ids) { + const auto offset = + checked_product(static_cast(item), component_size_); + output.insert( + output.end(), + values_.begin() + static_cast(offset), + values_.begin() + static_cast(offset + component_size_) + ); + } + return {name_, std::move(output), std::move(shape)}; + } + + [[nodiscard]] MaterializedValues value(const std::uint64_t item) const { + std::vector shape = component_shape_; + const auto offset = + checked_product(static_cast(item), component_size_); + std::vector output( + values_.begin() + static_cast(offset), + values_.begin() + static_cast(offset + component_size_) + ); + return {name_, std::move(output), std::move(shape)}; + } + +private: + static std::size_t checked_product( + const std::size_t first, + const std::size_t second + ) { + if (second != 0 && + first > std::numeric_limits::max() / second) { + throw std::overflow_error("value-map shape overflows size_t"); + } + return first * second; + } + + static std::size_t shape_size(const std::vector &shape) { + std::size_t size = 1; + for (const auto extent : shape) { + size = checked_product(size, extent); + } + return size; + } + + std::string name_; + std::vector component_shape_; + std::size_t component_size_; + std::vector values_; + std::vector counts_; + std::vector mean_sums_; + Reduction reduction_; +}; + +using ReducedValueMapVariant = std::variant< + ReducedValueMap, + ReducedValueMap +>; + +// Stateful graph contraction with named floating-point value maps. The input +// graph is simple. Parallel edges created by mutations can merge immediately or +// remain separate until materialization. +class ContractionGraph { +public: + using NodeId = std::uint64_t; + using EdgeId = std::uint64_t; + using Edge = UndirectedGraph::Edge; + + explicit ContractionGraph( + const UndirectedGraph &graph, + const ParallelEdgePolicy parallel_edge_policy = ParallelEdgePolicy::Merge + ) + : topology_(graph, parallel_edge_policy), + node_sets_(static_cast(graph.number_of_nodes())), + edge_parent_(static_cast(graph.number_of_edges())), + edge_state_( + static_cast(graph.number_of_edges()), + EdgeState::Active + ), + parallel_edge_policy_(parallel_edge_policy) { + initialize_edge_parents(); + } + + [[nodiscard]] std::size_t number_of_nodes() const { + return topology_.number_of_nodes(); + } + + [[nodiscard]] std::size_t number_of_edges() const { + return topology_.number_of_edges(); + } + + [[nodiscard]] std::vector nodes() const { + return topology_.active_nodes(); + } + + [[nodiscard]] std::vector edges() const { + return topology_.active_edges(); + } + + [[nodiscard]] Edge uv(const EdgeId edge) const { + return topology_.uv(edge); + } + + [[nodiscard]] std::int64_t find_edge(const NodeId u, const NodeId v) const { + return topology_.find_edge(u, v); + } + + [[nodiscard]] std::vector find_edges( + const NodeId u, + const NodeId v + ) const { + return topology_.find_edges(u, v); + } + + [[nodiscard]] ParallelEdgePolicy parallel_edge_policy() const { + return parallel_edge_policy_; + } + + [[nodiscard]] const std::vector &node_adjacency( + const NodeId node + ) const { + return topology_.node_adjacency(node); + } + + [[nodiscard]] std::size_t degree(const NodeId node) const { + return topology_.degree(node); + } + + [[nodiscard]] bool node_active(const NodeId node) const { + return topology_.node_active(node); + } + + [[nodiscard]] bool edge_active(const EdgeId edge) const { + return topology_.edge_active(edge); + } + + NodeId representative(const NodeId original_node) { + if (original_node >= node_sets_.size()) { + throw std::out_of_range("original_node is out of range"); + } + return node_sets_.find(original_node); + } + + template + void add_node_values( + std::string name, + const std::span values, + std::vector component_shape, + const Reduction reduction + ) { + require_registration_allowed(name, node_values_, "node"); + node_values_.emplace_back( + ReducedValueMap( + std::move(name), + values, + std::move(component_shape), + topology_.node_capacity(), + reduction + ) + ); + } + + template + void add_edge_values( + std::string name, + const std::span values, + std::vector component_shape, + const Reduction reduction + ) { + require_registration_allowed(name, edge_values_, "edge"); + edge_values_.emplace_back( + ReducedValueMap( + std::move(name), + values, + std::move(component_shape), + topology_.edge_capacity(), + reduction + ) + ); + } + + NodeId contract_edge( + const EdgeId edge, + const std::optional keep_node = std::nullopt + ) { + mutation_started_ = true; + const auto kept = topology_.contract_edge(edge, keep_node, change_); + apply_node_reduction(change_.kept_node, change_.removed_node); + node_sets_.merge_to(change_.kept_node, change_.removed_node); + apply_removed_edge_states(EdgeState::Internal); + apply_edge_folds(); + return kept; + } + + void erase_edge(const EdgeId edge) { + mutation_started_ = true; + topology_.erase_edge(edge, change_); + set_edge_state(edge, EdgeState::Deleted); + } + + EdgeId suppress_node(const NodeId node) { + mutation_started_ = true; + const auto replacement = topology_.suppress_node(node, change_); + apply_edge_folds(); + return replacement; + } + + [[nodiscard]] const ReducedValueMapVariant &node_value_map( + const std::string &name + ) const { + return find_value_map(node_values_, name, "node"); + } + + [[nodiscard]] const ReducedValueMapVariant &edge_value_map( + const std::string &name + ) const { + return find_value_map(edge_values_, name, "edge"); + } + + [[nodiscard]] MaterializedValueMap node_value( + const std::string &name, + const NodeId node + ) const { + if (!topology_.node_active(node)) { + throw std::invalid_argument("node is inactive"); + } + return std::visit([node](const auto &map) { + return MaterializedValueMap(map.value(node)); + }, node_value_map(name)); + } + + [[nodiscard]] MaterializedValueMap edge_value( + const std::string &name, + const EdgeId edge + ) const { + if (!topology_.edge_active(edge)) { + throw std::invalid_argument("edge is inactive"); + } + return std::visit([edge](const auto &map) { + return MaterializedValueMap(map.value(edge)); + }, edge_value_map(name)); + } + + [[nodiscard]] MaterializedValueMap active_node_values( + const std::string &name + ) const { + const auto ids = topology_.active_nodes(); + return std::visit([&ids](const auto &map) { + return MaterializedValueMap(map.gather(ids)); + }, node_value_map(name)); + } + + [[nodiscard]] MaterializedValueMap active_edge_values( + const std::string &name + ) const { + const auto ids = topology_.active_edges(); + return std::visit([&ids](const auto &map) { + return MaterializedValueMap(map.gather(ids)); + }, edge_value_map(name)); + } + + // Build an independent simple-graph snapshot. Parallel active edges fold + // only in the snapshot. Multiple input items can map to one output id. + MaterializedContraction materialize() { + const auto active_nodes = topology_.active_nodes(); + std::vector dense_node(topology_.node_capacity(), -1); + for (std::size_t dense = 0; dense < active_nodes.size(); ++dense) { + dense_node[static_cast(active_nodes[dense])] = + static_cast(dense); + } + + struct MaterializedEdge { + Edge uv; + EdgeId stable_edge; + }; + std::vector materialized_edges; + materialized_edges.reserve(topology_.number_of_edges()); + for (const auto edge : topology_.active_edges()) { + const auto uv = topology_.uv(edge); + const auto dense_u = + static_cast(dense_node[static_cast(uv.first)]); + const auto dense_v = + static_cast(dense_node[static_cast(uv.second)]); + if (dense_u == dense_v) { + throw std::runtime_error( + "cannot materialize a graph with active self edges" + ); + } + materialized_edges.push_back( + {{std::min(dense_u, dense_v), std::max(dense_u, dense_v)}, edge} + ); + } + std::sort( + materialized_edges.begin(), + materialized_edges.end(), + [](const auto &a, const auto &b) { + return a.uv != b.uv + ? a.uv < b.uv + : a.stable_edge < b.stable_edge; + } + ); + + std::vector output_edges; + std::vector active_edges; + std::vector> snapshot_folds; + output_edges.reserve(materialized_edges.size()); + active_edges.reserve(materialized_edges.size()); + snapshot_folds.reserve(materialized_edges.size()); + std::vector dense_edge(topology_.edge_capacity(), -1); + for (const auto &edge : materialized_edges) { + if (output_edges.empty() || output_edges.back() != edge.uv) { + output_edges.push_back(edge.uv); + active_edges.push_back(edge.stable_edge); + } else { + snapshot_folds.emplace_back(active_edges.back(), edge.stable_edge); + } + dense_edge[static_cast(edge.stable_edge)] = + static_cast(output_edges.size() - 1); + } + + auto graph = UndirectedGraph::from_sorted_unique_edges( + static_cast(active_nodes.size()), + std::move(output_edges), + true + ); + + std::vector node_mapping(topology_.node_capacity(), -1); + for (NodeId original = 0; original < topology_.node_capacity(); ++original) { + const auto root = node_sets_.find(original); + if (topology_.node_active(root)) { + node_mapping[static_cast(original)] = + dense_node[static_cast(root)]; + } + } + + std::vector edge_mapping(topology_.edge_capacity(), -1); + for (EdgeId original = 0; original < topology_.edge_capacity(); ++original) { + const auto root = edge_root(original); + if (edge_state_[static_cast(root)] == EdgeState::Active && + topology_.edge_active(root)) { + edge_mapping[static_cast(original)] = + dense_edge[static_cast(root)]; + } + } + + std::vector node_values; + node_values.reserve(node_values_.size()); + for (const auto &map : node_values_) { + node_values.push_back(std::visit( + [&active_nodes](const auto &typed) { + return MaterializedValueMap(typed.gather(active_nodes)); + }, + map + )); + } + + std::vector edge_values; + edge_values.reserve(edge_values_.size()); + for (const auto &map : edge_values_) { + edge_values.push_back(std::visit( + [&active_edges, &snapshot_folds](const auto &typed) { + auto snapshot = typed; + for (const auto &[kept, removed] : snapshot_folds) { + snapshot.reduce(kept, removed); + } + return MaterializedValueMap(snapshot.gather(active_edges)); + }, + map + )); + } + + return { + std::move(graph), + std::move(node_values), + std::move(edge_values), + std::move(node_mapping), + std::move(edge_mapping), + }; + } + +private: + enum class EdgeState : unsigned char { + Active, + Internal, + Deleted, + }; + + template + static decltype(auto) visit_value( + const ReducedValueMapVariant &map, + Function &&function + ) { + return std::visit(std::forward(function), map); + } + + static const std::string &map_name(const ReducedValueMapVariant &map) { + return visit_value(map, [](const auto &typed) -> const std::string & { + return typed.name(); + }); + } + + static const ReducedValueMapVariant &find_value_map( + const std::vector &maps, + const std::string &name, + const char *kind + ) { + for (const auto &map : maps) { + if (map_name(map) == name) { + return map; + } + } + throw std::invalid_argument( + std::string("unknown ") + kind + " value map: " + name + ); + } + + void require_registration_allowed( + const std::string &name, + const std::vector &maps, + const char *kind + ) const { + if (mutation_started_) { + throw std::invalid_argument( + "value maps must be registered before the first graph mutation" + ); + } + if (name.empty()) { + throw std::invalid_argument("value map name must not be empty"); + } + for (const auto &map : maps) { + if (map_name(map) == name) { + throw std::invalid_argument( + std::string(kind) + " value map already exists: " + name + ); + } + } + } + + void initialize_edge_parents() { + for (EdgeId edge = 0; edge < edge_parent_.size(); ++edge) { + edge_parent_[static_cast(edge)] = edge; + } + } + + EdgeId edge_root(EdgeId edge) { + auto current = static_cast(edge); + while (edge_parent_[current] != current) { + edge_parent_[current] = + edge_parent_[static_cast(edge_parent_[current])]; + current = static_cast(edge_parent_[current]); + } + return static_cast(current); + } + + void set_edge_state(const EdgeId edge, const EdgeState state) { + const auto root = edge_root(edge); + edge_state_[static_cast(root)] = state; + } + + void apply_removed_edge_states(const EdgeState state) { + for (const auto edge : change_.deactivated_edges) { + const auto folded = std::any_of( + change_.folded_edges.begin(), + change_.folded_edges.end(), + [edge](const auto &fold) { + return fold.removed_edge == edge; + } + ); + if (!folded) { + set_edge_state(edge, state); + } + } + } + + void apply_node_reduction(const NodeId kept, const NodeId removed) { + for (auto &map : node_values_) { + std::visit( + [kept, removed](auto &typed) { + typed.reduce(kept, removed); + }, + map + ); + } + } + + void apply_edge_folds() { + for (const auto &fold : change_.folded_edges) { + for (auto &map : edge_values_) { + std::visit( + [&fold](auto &typed) { + typed.reduce(fold.kept_edge, fold.removed_edge); + }, + map + ); + } + const auto kept_root = edge_root(fold.kept_edge); + const auto removed_root = edge_root(fold.removed_edge); + if (kept_root != removed_root) { + edge_parent_[static_cast(removed_root)] = kept_root; + } + } + } + + detail::ContractionTopology topology_; + util::UnionFind node_sets_; + std::vector edge_parent_; + std::vector edge_state_; + std::vector node_values_; + std::vector edge_values_; + detail::ContractionChange change_; + ParallelEdgePolicy parallel_edge_policy_ = ParallelEdgePolicy::Merge; + bool mutation_started_ = false; +}; + +} // namespace bioimage_cpp::graph diff --git a/include/bioimage_cpp/graph/contraction_options.hxx b/include/bioimage_cpp/graph/contraction_options.hxx new file mode 100644 index 0000000..516a394 --- /dev/null +++ b/include/bioimage_cpp/graph/contraction_options.hxx @@ -0,0 +1,10 @@ +#pragma once + +namespace bioimage_cpp::graph { + +enum class ParallelEdgePolicy { + Merge, + Keep, +}; + +} // namespace bioimage_cpp::graph diff --git a/include/bioimage_cpp/graph/detail/contraction_topology.hxx b/include/bioimage_cpp/graph/detail/contraction_topology.hxx new file mode 100644 index 0000000..2047a95 --- /dev/null +++ b/include/bioimage_cpp/graph/detail/contraction_topology.hxx @@ -0,0 +1,526 @@ +#pragma once + +#include "bioimage_cpp/graph/contraction_options.hxx" +#include "bioimage_cpp/graph/undirected_graph.hxx" + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace bioimage_cpp::graph::detail { + +inline constexpr std::uint64_t invalid_contraction_id = + std::numeric_limits::max(); + +struct ContractionAdjacency { + std::uint64_t node; + std::uint64_t edge; +}; + +struct ContractionEdge { + std::uint64_t u = 0; + std::uint64_t v = 0; + bool active = false; +}; + +struct EdgeRekey { + std::uint64_t edge; + std::uint64_t old_u; + std::uint64_t old_v; + std::uint64_t new_u; + std::uint64_t new_v; +}; + +struct EdgeFold { + std::uint64_t kept_edge; + std::uint64_t removed_edge; +}; + +// Callers reuse this record across mutations to avoid one allocation per edge +// contraction. Solver adapters can use the same record to update heaps and +// objective-specific state after the topology changes. +struct ContractionChange { + std::uint64_t kept_node = invalid_contraction_id; + std::uint64_t removed_node = invalid_contraction_id; + std::uint64_t removed_edge = invalid_contraction_id; + std::uint64_t replacement_edge = invalid_contraction_id; + std::vector deactivated_edges; + std::vector rekeyed_edges; + std::vector folded_edges; + + void clear() { + kept_node = invalid_contraction_id; + removed_node = invalid_contraction_id; + removed_edge = invalid_contraction_id; + replacement_edge = invalid_contraction_id; + deactivated_edges.clear(); + rekeyed_edges.clear(); + folded_edges.clear(); + } +}; + +// Mutable topology with stable sparse ids. Parallel-edge merging is optional. +// The class owns no values, heaps, or component labels. Callers update that +// state from the reusable ContractionChange record after each mutation. +class ContractionTopology { +public: + using NodeId = std::uint64_t; + using EdgeId = std::uint64_t; + using Edge = UndirectedGraph::Edge; + + ContractionTopology() = default; + + explicit ContractionTopology( + const UndirectedGraph &graph, + const ParallelEdgePolicy parallel_policy = ParallelEdgePolicy::Merge + ) + : ContractionTopology( + graph.number_of_nodes(), + graph.uv_ids(), + parallel_policy + ) { + } + + ContractionTopology( + const NodeId number_of_nodes, + const std::vector &edges, + const ParallelEdgePolicy parallel_policy + ) + : adjacency_(static_cast(number_of_nodes)), + nodes_active_(static_cast(number_of_nodes), true), + active_node_count_(static_cast(number_of_nodes)), + parallel_policy_(parallel_policy), + scratch_edge_(static_cast(number_of_nodes), invalid_contraction_id) { + edges_.reserve(edges.size()); + for (const auto &edge : edges) { + add_initial_edge(edge.first, edge.second); + } + active_edge_count_ = edges_.size(); + } + + [[nodiscard]] NodeId node_id_upper_bound() const { + return adjacency_.empty() ? 0 : static_cast(adjacency_.size() - 1); + } + + [[nodiscard]] EdgeId edge_id_upper_bound() const { + return edges_.empty() ? 0 : static_cast(edges_.size() - 1); + } + + [[nodiscard]] std::size_t node_capacity() const { + return adjacency_.size(); + } + + [[nodiscard]] std::size_t edge_capacity() const { + return edges_.size(); + } + + [[nodiscard]] std::size_t number_of_nodes() const { + return active_node_count_; + } + + [[nodiscard]] std::size_t number_of_edges() const { + return active_edge_count_; + } + + [[nodiscard]] bool node_active(const NodeId node) const { + validate_node_id(node); + return nodes_active_[static_cast(node)]; + } + + [[nodiscard]] bool edge_active(const EdgeId edge) const { + validate_edge_id(edge); + return edges_[static_cast(edge)].active; + } + + [[nodiscard]] Edge uv(const EdgeId edge) const { + validate_active_edge(edge); + const auto &entry = edges_[static_cast(edge)]; + return {entry.u, entry.v}; + } + + [[nodiscard]] std::size_t degree(const NodeId node) const { + validate_active_node(node); + return adjacency_[static_cast(node)].size(); + } + + [[nodiscard]] const std::vector &node_adjacency( + const NodeId node + ) const { + validate_active_node(node); + return adjacency_[static_cast(node)]; + } + + [[nodiscard]] std::vector active_nodes() const { + std::vector result; + result.reserve(active_node_count_); + for (NodeId node = 0; node < adjacency_.size(); ++node) { + if (nodes_active_[static_cast(node)]) { + result.push_back(node); + } + } + return result; + } + + [[nodiscard]] std::vector active_edges() const { + std::vector result; + result.reserve(active_edge_count_); + for (EdgeId edge = 0; edge < edges_.size(); ++edge) { + if (edges_[static_cast(edge)].active) { + result.push_back(edge); + } + } + return result; + } + + [[nodiscard]] std::int64_t find_edge(const NodeId u, const NodeId v) const { + validate_active_node(u); + validate_active_node(v); + auto result = invalid_contraction_id; + for (const auto &adjacency : adjacency_[static_cast(u)]) { + if (adjacency.node == v) { + result = std::min(result, adjacency.edge); + } + } + return result == invalid_contraction_id + ? -1 + : static_cast(result); + } + + [[nodiscard]] std::vector find_edges( + const NodeId u, + const NodeId v + ) const { + validate_active_node(u); + validate_active_node(v); + std::vector result; + for (const auto &adjacency : adjacency_[static_cast(u)]) { + if (adjacency.node == v) { + result.push_back(adjacency.edge); + } + } + std::sort(result.begin(), result.end()); + return result; + } + + NodeId contract_edge( + const EdgeId edge, + const std::optional keep_node, + ContractionChange &change + ) { + change.clear(); + validate_active_edge(edge); + const auto contracted = edges_[static_cast(edge)]; + if (contracted.u == contracted.v) { + throw std::invalid_argument("cannot contract a self edge"); + } + + auto kept = contracted.u; + auto removed = contracted.v; + if (keep_node.has_value()) { + if (*keep_node != contracted.u && *keep_node != contracted.v) { + throw std::invalid_argument( + "keep_node must be an endpoint of the contracted edge" + ); + } + kept = *keep_node; + removed = kept == contracted.u ? contracted.v : contracted.u; + } else { + const auto degree_u = adjacency_[static_cast(contracted.u)].size(); + const auto degree_v = adjacency_[static_cast(contracted.v)].size(); + if (degree_v > degree_u || + (degree_v == degree_u && contracted.v < contracted.u)) { + kept = contracted.v; + removed = contracted.u; + } + } + + change.kept_node = kept; + change.removed_node = removed; + change.removed_edge = edge; + + deactivate_edge(edge); + change.deactivated_edges.push_back(edge); + + if (parallel_policy_ == ParallelEdgePolicy::Merge) { + stamp_neighbors(kept); + } + + const auto removed_adjacency = + adjacency_[static_cast(removed)]; + for (const auto &adjacency : removed_adjacency) { + const auto current_edge = adjacency.edge; + if (!edges_[static_cast(current_edge)].active) { + continue; + } + const auto neighbor = adjacency.node; + if (neighbor == kept) { + deactivate_edge(current_edge); + change.deactivated_edges.push_back(current_edge); + continue; + } + + const auto existing = + parallel_policy_ == ParallelEdgePolicy::Merge + ? scratch_edge_[static_cast(neighbor)] + : invalid_contraction_id; + if (existing != invalid_contraction_id) { + deactivate_edge(current_edge); + change.deactivated_edges.push_back(current_edge); + change.folded_edges.push_back({existing, current_edge}); + continue; + } + + const auto old_edge = edges_[static_cast(current_edge)]; + rekey_edge_endpoint(current_edge, removed, kept); + adjacency_[static_cast(kept)].push_back( + {neighbor, current_edge} + ); + if (parallel_policy_ == ParallelEdgePolicy::Merge) { + scratch_edge_[static_cast(neighbor)] = current_edge; + } + const auto &new_edge = edges_[static_cast(current_edge)]; + change.rekeyed_edges.push_back( + { + current_edge, + old_edge.u, + old_edge.v, + new_edge.u, + new_edge.v, + } + ); + } + + if (parallel_policy_ == ParallelEdgePolicy::Merge) { + clear_stamps(kept); + } + adjacency_[static_cast(removed)].clear(); + nodes_active_[static_cast(removed)] = false; + --active_node_count_; + return kept; + } + + void erase_edge(const EdgeId edge, ContractionChange &change) { + change.clear(); + validate_active_edge(edge); + change.removed_edge = edge; + deactivate_edge(edge); + change.deactivated_edges.push_back(edge); + } + + EdgeId suppress_node(const NodeId node, ContractionChange &change) { + change.clear(); + validate_active_node(node); + const auto incident = adjacency_[static_cast(node)]; + if (incident.size() != 2) { + throw std::invalid_argument( + "suppress_node requires degree 2, got degree=" + + std::to_string(incident.size()) + ); + } + if (incident[0].edge == incident[1].edge) { + throw std::invalid_argument("cannot suppress a node incident to a self edge"); + } + + const auto edge_a = incident[0].edge; + const auto edge_b = incident[1].edge; + const auto neighbor_a = incident[0].node; + const auto neighbor_b = incident[1].node; + if (neighbor_a == neighbor_b) { + throw std::invalid_argument( + "suppress_node would create a self edge" + ); + } + + change.removed_node = node; + deactivate_edge(edge_a); + deactivate_edge(edge_b); + change.deactivated_edges.push_back(edge_a); + change.deactivated_edges.push_back(edge_b); + + EdgeId replacement = invalid_contraction_id; + if (parallel_policy_ == ParallelEdgePolicy::Merge && + neighbor_a != neighbor_b) { + const auto found = find_edge(neighbor_a, neighbor_b); + if (found >= 0) { + replacement = static_cast(found); + change.folded_edges.push_back({replacement, edge_a}); + change.folded_edges.push_back({replacement, edge_b}); + } + } + + if (replacement == invalid_contraction_id) { + replacement = std::min(edge_a, edge_b); + const auto folded = replacement == edge_a ? edge_b : edge_a; + auto &replacement_entry = edges_[static_cast(replacement)]; + const auto old_u = replacement_entry.u; + const auto old_v = replacement_entry.v; + replacement_entry.u = neighbor_a; + replacement_entry.v = neighbor_b; + replacement_entry.active = true; + ++active_edge_count_; + add_adjacency(neighbor_a, neighbor_b, replacement); + if (old_u != node && old_v != node) { + throw std::runtime_error( + "suppressed edge is not incident to the suppressed node" + ); + } + change.rekeyed_edges.push_back( + { + replacement, + old_u, + old_v, + neighbor_a, + neighbor_b, + } + ); + change.folded_edges.push_back({replacement, folded}); + } + + adjacency_[static_cast(node)].clear(); + nodes_active_[static_cast(node)] = false; + --active_node_count_; + change.replacement_edge = replacement; + return replacement; + } + +private: + void validate_node_id(const NodeId node) const { + if (node >= adjacency_.size()) { + throw std::out_of_range( + "node id " + std::to_string(node) + " is out of range" + ); + } + } + + void validate_edge_id(const EdgeId edge) const { + if (edge >= edges_.size()) { + throw std::out_of_range( + "edge id " + std::to_string(edge) + " is out of range" + ); + } + } + + void validate_active_node(const NodeId node) const { + validate_node_id(node); + if (!nodes_active_[static_cast(node)]) { + throw std::invalid_argument( + "node " + std::to_string(node) + " is inactive" + ); + } + } + + void validate_active_edge(const EdgeId edge) const { + validate_edge_id(edge); + if (!edges_[static_cast(edge)].active) { + throw std::invalid_argument( + "edge " + std::to_string(edge) + " is inactive" + ); + } + } + + void add_initial_edge(const NodeId u, const NodeId v) { + if (u >= adjacency_.size() || v >= adjacency_.size()) { + throw std::out_of_range("edge endpoint is out of range"); + } + if (u == v) { + throw std::invalid_argument("self edges are not supported"); + } + if (parallel_policy_ == ParallelEdgePolicy::Merge) { + for (const auto &adjacency : adjacency_[static_cast(u)]) { + if (adjacency.node == v) { + throw std::invalid_argument("parallel edges are not supported"); + } + } + } + const auto edge = static_cast(edges_.size()); + edges_.push_back({u, v, true}); + add_adjacency(u, v, edge); + } + + void add_adjacency(const NodeId u, const NodeId v, const EdgeId edge) { + adjacency_[static_cast(u)].push_back({v, edge}); + adjacency_[static_cast(v)].push_back({u, edge}); + } + + static void erase_adjacency_edge( + std::vector &adjacency, + const EdgeId edge + ) { + for (std::size_t index = 0; index < adjacency.size(); ++index) { + if (adjacency[index].edge == edge) { + adjacency[index] = adjacency.back(); + adjacency.pop_back(); + return; + } + } + } + + void deactivate_edge(const EdgeId edge) { + auto &entry = edges_[static_cast(edge)]; + if (!entry.active) { + return; + } + erase_adjacency_edge(adjacency_[static_cast(entry.u)], edge); + erase_adjacency_edge(adjacency_[static_cast(entry.v)], edge); + entry.active = false; + --active_edge_count_; + } + + void rekey_edge_endpoint( + const EdgeId edge, + const NodeId old_node, + const NodeId new_node + ) { + auto &entry = edges_[static_cast(edge)]; + if (entry.u == old_node) { + entry.u = new_node; + } else if (entry.v == old_node) { + entry.v = new_node; + } else { + throw std::runtime_error("edge is not incident to the rekeyed node"); + } + const auto other = entry.u == new_node ? entry.v : entry.u; + auto &other_adjacency = adjacency_[static_cast(other)]; + for (auto &adjacency : other_adjacency) { + if (adjacency.edge == edge) { + adjacency.node = new_node; + return; + } + } + throw std::runtime_error("edge adjacency is inconsistent"); + } + + void stamp_neighbors(const NodeId node) { + for (const auto &adjacency : adjacency_[static_cast(node)]) { + if (adjacency.node != node) { + scratch_edge_[static_cast(adjacency.node)] = + adjacency.edge; + } + } + } + + void clear_stamps(const NodeId node) { + for (const auto &adjacency : adjacency_[static_cast(node)]) { + if (adjacency.node != node) { + scratch_edge_[static_cast(adjacency.node)] = + invalid_contraction_id; + } + } + } + + std::vector> adjacency_; + std::vector edges_; + std::vector nodes_active_; + std::size_t active_node_count_ = 0; + std::size_t active_edge_count_ = 0; + ParallelEdgePolicy parallel_policy_ = ParallelEdgePolicy::Merge; + std::vector scratch_edge_; +}; + +} // namespace bioimage_cpp::graph::detail diff --git a/src/bindings/graph.hxx b/src/bindings/graph.hxx index 8badc5e..dd4c48f 100644 --- a/src/bindings/graph.hxx +++ b/src/bindings/graph.hxx @@ -5,5 +5,6 @@ namespace bioimage_cpp::bindings { void bind_graph(nanobind::module_ &m); +void bind_graph_contraction(nanobind::module_ &m); } // namespace bioimage_cpp::bindings diff --git a/src/bindings/graph_contraction.cxx b/src/bindings/graph_contraction.cxx new file mode 100644 index 0000000..9bceea4 --- /dev/null +++ b/src/bindings/graph_contraction.cxx @@ -0,0 +1,332 @@ +#include "graph.hxx" +#include "ndarray.hxx" + +#include "bioimage_cpp/graph/contraction_graph.hxx" + +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace nb = nanobind; + +namespace bioimage_cpp::bindings { +namespace { + +using Graph = graph::UndirectedGraph; +using ContractionGraph = graph::ContractionGraph; +using UInt64Array = nb::ndarray; + +template +using ConstFloatingArray = nb::ndarray; + +graph::ParallelEdgePolicy parallel_edge_policy_from_code(const int code) { + switch (code) { + case 0: + return graph::ParallelEdgePolicy::Merge; + case 1: + return graph::ParallelEdgePolicy::Keep; + default: + throw std::invalid_argument( + "unknown parallel edge policy code: " + std::to_string(code) + ); + } +} + +graph::Reduction reduction_from_code(const int code) { + switch (code) { + case 0: + return graph::Reduction::Sum; + case 1: + return graph::Reduction::Mean; + case 2: + return graph::Reduction::Minimum; + case 3: + return graph::Reduction::Maximum; + default: + throw std::invalid_argument("unknown reduction code: " + std::to_string(code)); + } +} + +template +void add_values( + ContractionGraph &contraction, + const std::string &name, + const ConstFloatingArray values, + const int reduction_code, + const bool node_values +) { + const auto expected = node_values + ? contraction.nodes().size() + : contraction.edges().size(); + if (values.ndim() == 0) { + throw std::invalid_argument("values must have at least one dimension"); + } + if (values.shape(0) != expected) { + throw std::invalid_argument( + "values.shape[0] must match the original number of " + + std::string(node_values ? "nodes" : "edges") + + ", got " + std::to_string(values.shape(0)) + + " for " + std::to_string(expected) + ); + } + + std::vector component_shape; + component_shape.reserve(values.ndim() - 1); + std::size_t size = 1; + for (std::size_t axis = 0; axis < values.ndim(); ++axis) { + size *= values.shape(axis); + if (axis > 0) { + component_shape.push_back(values.shape(axis)); + } + } + const auto input = std::span(values.data(), size); + if (node_values) { + contraction.add_node_values( + name, + input, + std::move(component_shape), + reduction_from_code(reduction_code) + ); + } else { + contraction.add_edge_values( + name, + input, + std::move(component_shape), + reduction_from_code(reduction_code) + ); + } +} + +UInt64Array ids_to_array(const std::vector &ids) { + return detail::copy_vector_to_array(ids); +} + +UInt64Array adjacency_to_array( + const std::vector &adjacency +) { + auto output = detail::make_array_for_overwrite( + {adjacency.size(), std::size_t{2}} + ); + for (std::size_t index = 0; index < adjacency.size(); ++index) { + output.data()[2 * index] = adjacency[index].node; + output.data()[2 * index + 1] = adjacency[index].edge; + } + return output; +} + +template +nb::object values_to_array(const graph::MaterializedValues &values) { + return nb::cast(detail::copy_vector_to_array(values.values, values.shape)); +} + +nb::object values_to_array(const graph::MaterializedValueMap &values) { + return std::visit( + [](const auto &typed) { + return values_to_array(typed); + }, + values + ); +} + +nb::tuple active_node_values( + const ContractionGraph &contraction, + const std::string &name +) { + return nb::make_tuple( + ids_to_array(contraction.nodes()), + values_to_array(contraction.active_node_values(name)) + ); +} + +nb::tuple active_edge_values( + const ContractionGraph &contraction, + const std::string &name +) { + return nb::make_tuple( + ids_to_array(contraction.edges()), + values_to_array(contraction.active_edge_values(name)) + ); +} + +UInt64Array graph_uv_ids(const Graph &graph) { + auto output = detail::make_array_for_overwrite( + {static_cast(graph.number_of_edges()), std::size_t{2}} + ); + for (std::uint64_t edge = 0; edge < graph.number_of_edges(); ++edge) { + const auto uv = graph.uv(edge); + output.data()[2 * edge] = uv.first; + output.data()[2 * edge + 1] = uv.second; + } + return output; +} + +nb::dict materialized_values_to_dict( + const std::vector &values +) { + nb::dict output; + for (const auto &value : values) { + std::visit( + [&output](const auto &typed) { + output[nb::str(typed.name.c_str())] = values_to_array(typed); + }, + value + ); + } + return output; +} + +nb::tuple materialize(ContractionGraph &contraction) { + auto result = contraction.materialize(); + const auto number_of_nodes = result.graph.number_of_nodes(); + auto uvs = graph_uv_ids(result.graph); + auto node_values = materialized_values_to_dict(result.node_values); + auto edge_values = materialized_values_to_dict(result.edge_values); + auto node_mapping = detail::copy_vector_to_array(result.node_mapping); + auto edge_mapping = detail::copy_vector_to_array(result.edge_mapping); + return nb::make_tuple( + number_of_nodes, + std::move(uvs), + std::move(node_values), + std::move(edge_values), + std::move(node_mapping), + std::move(edge_mapping) + ); +} + +} // namespace + +void bind_graph_contraction(nb::module_ &m) { + nb::class_(m, "_ContractionGraph") + .def( + "__init__", + [](ContractionGraph *self, const Graph &graph, const int parallel_edges) { + new (self) ContractionGraph( + graph, + parallel_edge_policy_from_code(parallel_edges) + ); + }, + nb::arg("graph"), + nb::arg("parallel_edges") = 0 + ) + .def_prop_ro("number_of_nodes", &ContractionGraph::number_of_nodes) + .def_prop_ro("number_of_edges", &ContractionGraph::number_of_edges) + .def("nodes", &ContractionGraph::nodes) + .def("edges", &ContractionGraph::edges) + .def("uv", &ContractionGraph::uv, nb::arg("edge")) + .def("find_edge", &ContractionGraph::find_edge, nb::arg("u"), nb::arg("v")) + .def( + "find_edges", + [](const ContractionGraph &self, const std::uint64_t u, + const std::uint64_t v) { + return ids_to_array(self.find_edges(u, v)); + }, + nb::arg("u"), + nb::arg("v") + ) + .def( + "node_adjacency", + [](const ContractionGraph &self, const std::uint64_t node) { + return adjacency_to_array(self.node_adjacency(node)); + }, + nb::arg("node") + ) + .def("degree", &ContractionGraph::degree, nb::arg("node")) + .def("is_node_active", &ContractionGraph::node_active, nb::arg("node")) + .def("is_edge_active", &ContractionGraph::edge_active, nb::arg("edge")) + .def( + "representative", + &ContractionGraph::representative, + nb::arg("original_node") + ) + .def( + "_add_node_values_float32", + [](ContractionGraph &self, const std::string &name, + const ConstFloatingArray values, const int reduction) { + add_values(self, name, values, reduction, true); + }, + nb::arg("name"), + nb::arg("values"), + nb::arg("reduction") + ) + .def( + "_add_node_values_float64", + [](ContractionGraph &self, const std::string &name, + const ConstFloatingArray values, const int reduction) { + add_values(self, name, values, reduction, true); + }, + nb::arg("name"), + nb::arg("values"), + nb::arg("reduction") + ) + .def( + "_add_edge_values_float32", + [](ContractionGraph &self, const std::string &name, + const ConstFloatingArray values, const int reduction) { + add_values(self, name, values, reduction, false); + }, + nb::arg("name"), + nb::arg("values"), + nb::arg("reduction") + ) + .def( + "_add_edge_values_float64", + [](ContractionGraph &self, const std::string &name, + const ConstFloatingArray values, const int reduction) { + add_values(self, name, values, reduction, false); + }, + nb::arg("name"), + nb::arg("values"), + nb::arg("reduction") + ) + .def( + "contract_edge", + &ContractionGraph::contract_edge, + nb::arg("edge"), + nb::arg("keep_node") = std::nullopt + ) + .def("erase_edge", &ContractionGraph::erase_edge, nb::arg("edge")) + .def("suppress_node", &ContractionGraph::suppress_node, nb::arg("node")) + .def( + "node_value", + [](const ContractionGraph &self, const std::string &name, + const std::uint64_t node) { + return values_to_array(self.node_value(name, node)); + }, + nb::arg("name"), + nb::arg("node") + ) + .def( + "edge_value", + [](const ContractionGraph &self, const std::string &name, + const std::uint64_t edge) { + return values_to_array(self.edge_value(name, edge)); + }, + nb::arg("name"), + nb::arg("edge") + ) + .def( + "active_node_values", + &active_node_values, + nb::arg("name") + ) + .def( + "active_edge_values", + &active_edge_values, + nb::arg("name") + ) + .def("materialize", &materialize); +} + +} // namespace bioimage_cpp::bindings diff --git a/src/bindings/module.cxx b/src/bindings/module.cxx index 1c8b56e..8a95fbe 100644 --- a/src/bindings/module.cxx +++ b/src/bindings/module.cxx @@ -25,6 +25,7 @@ NB_MODULE(_core, m) { bioimage_cpp::bindings::bind_filters(m); bioimage_cpp::bindings::bind_flow(m); bioimage_cpp::bindings::bind_graph(m); + bioimage_cpp::bindings::bind_graph_contraction(m); bioimage_cpp::bindings::bind_ground_truth(m); bioimage_cpp::bindings::bind_label_multiset(m); bioimage_cpp::bindings::bind_mesh(m); diff --git a/src/bioimage_cpp/graph/__init__.py b/src/bioimage_cpp/graph/__init__.py index e6fc841..674e2aa 100644 --- a/src/bioimage_cpp/graph/__init__.py +++ b/src/bioimage_cpp/graph/__init__.py @@ -2,8 +2,9 @@ Top-level surface: -- Graph structures: :class:`UndirectedGraph`, :class:`GridGraph2D`, - :class:`GridGraph3D`, :class:`RegionAdjacencyGraph`. +- Graph structures: :class:`UndirectedGraph`, :class:`ContractionGraph`, + :class:`GridGraph2D`, :class:`GridGraph3D`, + :class:`RegionAdjacencyGraph`. - Constructors: :func:`undirected_graph`, :func:`grid_graph`, :func:`region_adjacency_graph`. - Algorithms: :func:`connected_components`, :func:`breadth_first_search`, @@ -40,6 +41,8 @@ from __future__ import annotations +from dataclasses import dataclass + import numpy as np from .. import _core @@ -54,6 +57,7 @@ _as_uv_array, _normalize_labels, _normalize_number_of_threads, + _resolve_weight_dtype, ) @@ -120,6 +124,257 @@ def deserialize(cls, serialization): return cls.from_edges(number_of_nodes, uvs) +@dataclass(frozen=True) +class ContractionResult: + """Materialized graph and values from a :class:`ContractionGraph`. + + Rows in each node or edge value array match the corresponding ids in + ``graph``. ``node_mapping`` and ``edge_mapping`` map ids from the input + graph to these dense ids. Multiple input ids can map to one output id. An + entry is ``-1`` when no materialized item represents the input item. + """ + + graph: UndirectedGraph + node_values: dict[str, np.ndarray] + edge_values: dict[str, np.ndarray] + node_mapping: np.ndarray + edge_mapping: np.ndarray + + +_CONTRACTION_REDUCTIONS = { + "sum": 0, + "mean": 1, + "min": 2, + "max": 3, +} + +_CONTRACTION_PARALLEL_EDGES = { + "merge": 0, + "keep": 1, +} + + +def _contraction_values(values, name: str, expected_items: int) -> np.ndarray: + array = _resolve_weight_dtype(values, name) + if array.ndim == 0: + raise ValueError(f"{name} must have at least one dimension") + if array.shape[0] != expected_items: + raise ValueError( + f"{name}.shape[0] must be {expected_items}, got {array.shape[0]}" + ) + return np.ascontiguousarray(array) + + +def _contraction_reduction(reduction: str) -> int: + try: + return _CONTRACTION_REDUCTIONS[reduction] + except KeyError as error: + supported = ", ".join(_CONTRACTION_REDUCTIONS) + raise ValueError( + f"reduction must be one of ({supported}), got {reduction!r}" + ) from error + + +def _contraction_parallel_edges(parallel_edges: str) -> int: + try: + return _CONTRACTION_PARALLEL_EDGES[parallel_edges] + except KeyError as error: + supported = ", ".join(_CONTRACTION_PARALLEL_EDGES) + raise ValueError( + f"parallel_edges must be one of ({supported}), got {parallel_edges!r}" + ) from error + + +class ContractionGraph: + """Mutable graph for edge contractions and value reduction. + + The graph copies the input topology. Node and edge ids stay stable while + they are active. Mutations can leave gaps in the id range. The input graph + must be simple. + + Register value maps before the first mutation. A map must have shape + ``(number_of_items, ...)`` and a floating dtype. ``float32`` and + ``float64`` are preserved. The reduction is component-wise: + + - ``"sum"`` adds represented values. + - ``"mean"`` computes the mean over represented input items. + - ``"min"`` and ``"max"`` select the extrema. + + Node maps reduce when an edge contraction merges its endpoints. Edge maps + reduce when two surviving edges are folded into one. The contracted edge + has no surviving target, so its value is discarded. + + ``parallel_edges="merge"`` folds parallel edges during each mutation. + ``parallel_edges="keep"`` keeps them separate in the mutable graph. + :meth:`materialize` always creates a simple graph and folds remaining + parallel edges only in the returned snapshot. + + This class is not thread-safe. + """ + + def __init__( + self, + graph: UndirectedGraph, + *, + parallel_edges: str = "merge", + ): + parallel_edges_code = _contraction_parallel_edges(parallel_edges) + self._core = _core._ContractionGraph(graph, parallel_edges_code) + self._original_number_of_nodes = int(graph.number_of_nodes) + self._original_number_of_edges = int(graph.number_of_edges) + self._parallel_edges = parallel_edges + self._mutation_started = False + + @property + def parallel_edges(self) -> str: + """Policy for parallel edges created by mutations.""" + return self._parallel_edges + + @property + def number_of_nodes(self) -> int: + """Number of active nodes.""" + return int(self._core.number_of_nodes) + + @property + def number_of_edges(self) -> int: + """Number of active edges.""" + return int(self._core.number_of_edges) + + def nodes(self) -> np.ndarray: + """Return active stable node ids in ascending order.""" + return np.asarray(self._core.nodes(), dtype=np.uint64) + + def edges(self) -> np.ndarray: + """Return active stable edge ids in ascending order.""" + return np.asarray(self._core.edges(), dtype=np.uint64) + + def uv(self, edge: int) -> tuple[int, int]: + u, v = self._core.uv(int(edge)) + return int(u), int(v) + + def find_edge(self, u: int, v: int) -> int: + """Return the smallest matching active edge id, or ``-1``.""" + return int(self._core.find_edge(int(u), int(v))) + + def find_edges(self, u: int, v: int) -> np.ndarray: + """Return all matching active edge ids in ascending order.""" + return np.asarray(self._core.find_edges(int(u), int(v)), dtype=np.uint64) + + def node_adjacency(self, node: int) -> np.ndarray: + """Return ``(neighbor, edge)`` rows for an active node.""" + return self._core.node_adjacency(int(node)) + + def degree(self, node: int) -> int: + return int(self._core.degree(int(node))) + + def is_node_active(self, node: int) -> bool: + return bool(self._core.is_node_active(int(node))) + + def is_edge_active(self, edge: int) -> bool: + return bool(self._core.is_edge_active(int(edge))) + + def representative(self, original_node: int) -> int: + """Return the current representative of an input node.""" + return int(self._core.representative(int(original_node))) + + def add_node_values(self, name: str, values, *, reduction: str) -> None: + """Register a node value map before the first mutation.""" + if self._mutation_started: + raise ValueError( + "value maps must be registered before the first graph mutation" + ) + array = _contraction_values( + values, + name, + self._original_number_of_nodes, + ) + code = _contraction_reduction(reduction) + getattr(self._core, f"_add_node_values_{array.dtype.name}")( + str(name), array, code + ) + + def add_edge_values(self, name: str, values, *, reduction: str) -> None: + """Register an edge value map before the first mutation.""" + if self._mutation_started: + raise ValueError( + "value maps must be registered before the first graph mutation" + ) + array = _contraction_values( + values, + name, + self._original_number_of_edges, + ) + code = _contraction_reduction(reduction) + getattr(self._core, f"_add_edge_values_{array.dtype.name}")( + str(name), array, code + ) + + def contract_edge(self, edge: int, *, keep_node: int | None = None) -> int: + """Contract an active edge and return the retained node id.""" + keep = None if keep_node is None else int(keep_node) + result = int(self._core.contract_edge(int(edge), keep)) + self._mutation_started = True + return result + + def erase_edge(self, edge: int) -> None: + """Erase an active edge without merging its endpoints.""" + self._core.erase_edge(int(edge)) + self._mutation_started = True + + def suppress_node(self, node: int) -> int: + """Suppress a degree-2 node and return the replacement edge id. + + The operation removes the node and its two incident edges. It connects + the two neighbors and reduces both edge value rows into the replacement + edge. The suppressed node value has no target and is discarded. + """ + result = int(self._core.suppress_node(int(node))) + self._mutation_started = True + return result + + def node_value(self, name: str, node: int) -> np.ndarray: + """Return a copy of one active node value.""" + return self._core.node_value(str(name), int(node)) + + def edge_value(self, name: str, edge: int) -> np.ndarray: + """Return a copy of one active edge value.""" + return self._core.edge_value(str(name), int(edge)) + + def active_node_values(self, name: str) -> tuple[np.ndarray, np.ndarray]: + """Return active node ids and their value rows.""" + return self._core.active_node_values(str(name)) + + def active_edge_values(self, name: str) -> tuple[np.ndarray, np.ndarray]: + """Return active edge ids and their value rows.""" + return self._core.active_edge_values(str(name)) + + def materialize(self) -> ContractionResult: + """Create a dense graph and aligned copies of all registered values. + + The output graph is simple. The snapshot folds active parallel edges + with each edge map's registered reduction. + + Materialization does not change this contraction graph. The result + remains valid after later mutations. + """ + ( + number_of_nodes, + uvs, + node_values, + edge_values, + node_mapping, + edge_mapping, + ) = self._core.materialize() + graph = UndirectedGraph.from_unique_edges(int(number_of_nodes), uvs) + return ContractionResult( + graph=graph, + node_values=dict(node_values), + edge_values=dict(edge_values), + node_mapping=node_mapping, + edge_mapping=edge_mapping, + ) + + def _normalize_projection_offsets(offsets, ndim: int) -> list[list[int]]: return [list(offset) for offset in strict_offsets(offsets, ndim)] @@ -678,6 +933,8 @@ def project_node_labels_to_pixels( __all__ = [ + "ContractionGraph", + "ContractionResult", "GridGraph2D", "GridGraph3D", "RagCoordinates", diff --git a/src/bioimage_cpp/skeleton/postprocessing.py b/src/bioimage_cpp/skeleton/postprocessing.py index d70f802..d5b0042 100644 --- a/src/bioimage_cpp/skeleton/postprocessing.py +++ b/src/bioimage_cpp/skeleton/postprocessing.py @@ -8,11 +8,9 @@ from __future__ import annotations -from collections import defaultdict - import numpy as np -from ..graph import connected_components +from ..graph import ContractionGraph, UndirectedGraph, connected_components from ..utils import UnionFind from ._graph import skeleton_to_graph @@ -211,6 +209,30 @@ def _adjacency(num_nodes, edges): return indptr, dst, eid, degrees +def _nodes_in_noncycle_components(indptr, dst, degrees): + eligible = np.zeros(len(degrees), dtype=bool) + visited = np.zeros(len(degrees), dtype=bool) + for start in range(len(degrees)): + if visited[start]: + continue + component = [] + stack = [start] + visited[start] = True + has_critical_node = False + while stack: + node = stack.pop() + component.append(node) + has_critical_node |= degrees[node] != 2 + for index in range(indptr[node], indptr[node + 1]): + neighbor = int(dst[index]) + if not visited[neighbor]: + visited[neighbor] = True + stack.append(neighbor) + if has_critical_node: + eligible[component] = True + return eligible + + def remove_ticks(vertices, edges, tick_length, radii=None): """Prune short dead-end branches ("ticks") from a skeleton graph. @@ -227,7 +249,8 @@ def remove_ticks(vertices, edges, tick_length, radii=None): vertices: Float array with shape ``(V, D)`` of skeleton vertex coordinates. edges: - Integer array with shape ``(E, 2)`` indexing ``vertices``. + Integer array with shape ``(E, 2)`` indexing ``vertices``. Self-edges + and duplicate undirected edges are not supported. tick_length: Maximum branch length (physical) that may be pruned. radii: @@ -240,75 +263,70 @@ def remove_ticks(vertices, edges, tick_length, radii=None): ``None`` when no input radii were given. """ vertices = np.asarray(vertices, dtype=np.float64) - edges = np.asarray(edges, dtype=np.int64) + if vertices.ndim != 2: + raise ValueError(f"vertices must be a 2D array, got ndim={vertices.ndim}") + raw_edges = np.asarray(edges) + if not np.issubdtype(raw_edges.dtype, np.integer): + raise TypeError(f"edges must have an integer dtype, got dtype={raw_edges.dtype}") + if raw_edges.ndim != 2 or raw_edges.shape[1] != 2: + raise ValueError(f"edges must have shape (E, 2), got shape={raw_edges.shape}") + if np.issubdtype(raw_edges.dtype, np.signedinteger) and np.any(raw_edges < 0): + raise ValueError("edges must not contain negative node ids") + edges_u64 = np.ascontiguousarray(raw_edges, dtype=np.uint64) num_nodes = len(vertices) + if edges_u64.size and int(edges_u64.max()) >= num_nodes: + raise IndexError( + f"edge endpoint must be smaller than len(vertices)={num_nodes}" + ) + if radii is not None and len(radii) != num_nodes: + raise ValueError( + f"radii length must be {num_nodes}, got length={len(radii)}" + ) + edges = edges_u64.astype(np.int64, copy=False) if len(edges) == 0: return vertices, edges.copy(), radii - indptr, dst, eid, degrees = _adjacency(num_nodes, edges) - - # Distance supergraph: sid -> [end_a, end_b, length, edge_ids]. - supers = {} - incident = defaultdict(set) - edge_used = np.zeros(len(edges), dtype=bool) - sid = 0 - for node in map(int, np.where(degrees != 2)[0]): - for k in range(indptr[node], indptr[node + 1]): - if edge_used[eid[k]]: - continue - prev, cur = node, int(dst[k]) - path = [int(eid[k])] - length = float(np.linalg.norm(vertices[cur] - vertices[node])) - while degrees[cur] == 2: - s, e = indptr[cur], indptr[cur + 1] - nbrs, eids = dst[s:e], eid[s:e] - pick = 0 if int(nbrs[0]) != prev else 1 - path.append(int(eids[pick])) - nxt = int(nbrs[pick]) - length += float(np.linalg.norm(vertices[nxt] - vertices[cur])) - prev, cur = cur, nxt - for pe in path: - edge_used[pe] = True - supers[sid] = [node, cur, length, path] - incident[node].add(sid) - incident[cur].add(sid) - sid += 1 - - dropped = set() + graph = UndirectedGraph.from_edges(num_nodes, edges_u64) + if graph.number_of_edges != len(edges): + raise ValueError("edges must not contain duplicate undirected edges") + + indptr, dst, _, degrees = _adjacency(num_nodes, edges) + lengths = np.linalg.norm( + vertices[edges[:, 1]] - vertices[edges[:, 0]], + axis=1, + ) + work = ContractionGraph(graph, parallel_edges="keep") + work.add_edge_values("length", lengths, reduction="sum") + + eligible = _nodes_in_noncycle_components(indptr, dst, degrees) + for node in np.where((degrees == 2) & eligible)[0]: + node = int(node) + if work.is_node_active(node) and work.degree(node) == 2: + work.suppress_node(node) + + threshold = float(tick_length) while True: - best, best_len = None, tick_length - for s, (a, b, length, _) in supers.items(): - terminal_a, terminal_b = len(incident[a]) == 1, len(incident[b]) == 1 - if (terminal_a ^ terminal_b) and length < best_len: - best, best_len = s, length + edge_ids, active_lengths = work.active_edge_values("length") + best = None + for edge, length in zip(edge_ids, active_lengths, strict=True): + edge = int(edge) + length = float(length) + a, b = work.uv(edge) + terminal_a = work.degree(a) == 1 + terminal_b = work.degree(b) == 1 + if terminal_a ^ terminal_b and length < threshold: + candidate = (length, edge, a, b) + if best is None or candidate[:2] < best[:2]: + best = candidate if best is None: break - a, b, length, path = supers.pop(best) - incident[a].discard(best) - incident[b].discard(best) - dropped.update(path) + _, edge, a, b = best + work.erase_edge(edge) for node in (a, b): - if len(incident[node]) == 2: - s1, s2 = incident[node] - a1, b1, l1, p1 = supers.pop(s1) - a2, b2, l2, p2 = supers.pop(s2) - far1 = b1 if a1 == node else a1 - far2 = b2 if a2 == node else a2 - for x in (far1, far2, node): - incident[x].discard(s1) - incident[x].discard(s2) - supers[sid] = [far1, far2, l1 + l2, p1 + p2] - incident[far1].add(sid) - incident[far2].add(sid) - sid += 1 - - if dropped: - keep = np.ones(len(edges), dtype=bool) - keep[list(dropped)] = False - edges = edges[keep] - else: - edges = edges.copy() + if work.is_node_active(node) and work.degree(node) == 2: + work.suppress_node(node) + edges = edges[work.materialize().edge_mapping >= 0] vertices, edges, radii = _compact(vertices, edges, radii) return vertices, edges, radii diff --git a/tests/graph/test_contraction_graph.py b/tests/graph/test_contraction_graph.py new file mode 100644 index 0000000..54dd73a --- /dev/null +++ b/tests/graph/test_contraction_graph.py @@ -0,0 +1,298 @@ +import numpy as np +import pytest + +import bioimage_cpp as bic + + +def _triangle(): + return bic.graph.UndirectedGraph.from_edges( + 3, + np.array([[0, 1], [1, 2], [0, 2]], dtype=np.uint64), + ) + + +def _parallel_paths(): + return bic.graph.UndirectedGraph.from_edges( + 4, + np.array([[0, 2], [1, 2], [0, 3], [1, 3]], dtype=np.uint64), + ) + + +def test_contract_edge_materializes_aligned_values_and_mappings(): + graph = _triangle() + work = bic.graph.ContractionGraph(graph) + work.add_node_values( + "position", + np.array([[0.0, 2.0], [2.0, 4.0], [10.0, 12.0]], dtype=np.float32), + reduction="mean", + ) + work.add_edge_values( + "weight", + np.array([100.0, 2.0, 3.0], dtype=np.float64), + reduction="sum", + ) + + assert work.contract_edge(0, keep_node=0) == 0 + result = work.materialize() + + assert isinstance(result.graph, bic.graph.UndirectedGraph) + np.testing.assert_array_equal(result.graph.uv_ids(), [[0, 1]]) + np.testing.assert_allclose( + result.node_values["position"], + [[1.0, 3.0], [10.0, 12.0]], + ) + assert result.node_values["position"].dtype == np.float32 + np.testing.assert_array_equal(result.edge_values["weight"], [5.0]) + assert result.edge_values["weight"].dtype == np.float64 + np.testing.assert_array_equal(result.node_mapping, [0, 0, 1]) + np.testing.assert_array_equal(result.edge_mapping, [-1, 0, 0]) + + +@pytest.mark.parametrize("dtype", [np.float32, np.float64]) +@pytest.mark.parametrize( + "reduction, expected", + [ + ("sum", 8.0), + ("mean", 4.0), + ("min", 3.0), + ("max", 5.0), + ], +) +def test_edge_reductions_preserve_dtype(dtype, reduction, expected): + work = bic.graph.ContractionGraph(_triangle()) + work.add_edge_values( + "value", + np.array([100.0, 3.0, 5.0], dtype=dtype), + reduction=reduction, + ) + + work.contract_edge(0, keep_node=0) + values = work.materialize().edge_values["value"] + + assert values.dtype == np.dtype(dtype) + np.testing.assert_allclose(values, [expected]) + + +@pytest.mark.parametrize( + "reduction, expected", + [ + ("sum", 6.0), + ("mean", 3.0), + ("min", 2.0), + ("max", 4.0), + ], +) +def test_node_reductions(reduction, expected): + graph = bic.graph.UndirectedGraph.from_edges(2, [[0, 1]]) + work = bic.graph.ContractionGraph(graph) + work.add_node_values( + "value", + np.array([2.0, 4.0], dtype=np.float32), + reduction=reduction, + ) + + work.contract_edge(0) + + np.testing.assert_allclose(work.materialize().node_values["value"], [expected]) + + +def test_mean_tracks_the_number_of_represented_edges(): + graph = bic.graph.UndirectedGraph.from_edges( + 4, + [[0, 2], [1, 2], [0, 1], [0, 3], [2, 3]], + ) + work = bic.graph.ContractionGraph(graph) + work.add_edge_values( + "value", + np.array([1.0, 3.0, 100.0, 10.0, 100.0], dtype=np.float64), + reduction="mean", + ) + + work.contract_edge(2, keep_node=0) + work.contract_edge(4, keep_node=2) + + np.testing.assert_allclose( + work.materialize().edge_values["value"], + [(1.0 + 3.0 + 10.0) / 3.0], + ) + + +def test_suppress_node_reduces_edges_and_drops_node_value(): + graph = bic.graph.UndirectedGraph.from_edges(3, [[0, 1], [1, 2]]) + work = bic.graph.ContractionGraph(graph) + work.add_node_values( + "node", + np.array([1.0, 2.0, 3.0], dtype=np.float64), + reduction="sum", + ) + work.add_edge_values( + "length", + np.array([2.0, 3.0], dtype=np.float32), + reduction="sum", + ) + + assert work.suppress_node(1) == 0 + result = work.materialize() + + np.testing.assert_array_equal(result.graph.uv_ids(), [[0, 1]]) + np.testing.assert_array_equal(result.node_values["node"], [1.0, 3.0]) + np.testing.assert_array_equal(result.edge_values["length"], [5.0]) + np.testing.assert_array_equal(result.node_mapping, [0, -1, 1]) + np.testing.assert_array_equal(result.edge_mapping, [0, 0]) + + +def test_materialized_snapshot_is_independent(): + graph = bic.graph.UndirectedGraph.from_edges(3, [[0, 1], [1, 2]]) + work = bic.graph.ContractionGraph(graph) + work.add_edge_values( + "length", + np.array([2.0, 3.0], dtype=np.float64), + reduction="sum", + ) + replacement = work.suppress_node(1) + first = work.materialize() + + work.erase_edge(replacement) + second = work.materialize() + + assert first.graph.number_of_edges == 1 + np.testing.assert_array_equal(first.edge_values["length"], [5.0]) + assert second.graph.number_of_edges == 0 + assert second.edge_values["length"].shape == (0,) + np.testing.assert_array_equal(second.edge_mapping, [-1, -1]) + + +def test_active_value_queries_align_with_stable_ids(): + work = bic.graph.ContractionGraph(_triangle()) + work.add_edge_values( + "value", + np.array([10.0, 20.0, 30.0], dtype=np.float32), + reduction="sum", + ) + work.contract_edge(0, keep_node=0) + + edge_ids, values = work.active_edge_values("value") + + np.testing.assert_array_equal(edge_ids, [2]) + np.testing.assert_array_equal(values, [50.0]) + assert work.edge_value("value", 2).shape == () + assert float(work.edge_value("value", 2)) == 50.0 + + +def test_input_graph_is_copied(): + graph = bic.graph.UndirectedGraph.from_edges(3, [[0, 1]]) + work = bic.graph.ContractionGraph(graph) + + graph.insert_edge(1, 2) + + assert work.number_of_edges == 1 + assert work.find_edge(1, 2) == -1 + + +@pytest.mark.parametrize("dtype", [np.float32, np.float64]) +@pytest.mark.parametrize( + "reduction, expected_active, expected_materialized", + [ + ("sum", [3.0, 7.0], 10.0), + ("mean", [1.5, 3.5], 2.5), + ("min", [1.0, 3.0], 1.0), + ("max", [2.0, 4.0], 4.0), + ], +) +def test_keep_parallel_edges_fold_only_in_materialized_snapshot( + dtype, + reduction, + expected_active, + expected_materialized, +): + work = bic.graph.ContractionGraph( + _parallel_paths(), + parallel_edges="keep", + ) + work.add_edge_values( + "value", + np.array([1.0, 2.0, 3.0, 4.0], dtype=dtype), + reduction=reduction, + ) + + work.suppress_node(2) + work.suppress_node(3) + + assert work.parallel_edges == "keep" + assert work.find_edge(0, 1) == 0 + np.testing.assert_array_equal(work.find_edges(0, 1), [0, 2]) + edge_ids, active_values = work.active_edge_values("value") + np.testing.assert_array_equal(edge_ids, [0, 2]) + np.testing.assert_allclose(active_values, expected_active) + + result = work.materialize() + + np.testing.assert_array_equal(result.graph.uv_ids(), [[0, 1]]) + np.testing.assert_allclose( + result.edge_values["value"], + [expected_materialized], + ) + assert result.edge_values["value"].dtype == np.dtype(dtype) + np.testing.assert_array_equal(result.node_mapping, [0, 1, -1, -1]) + np.testing.assert_array_equal(result.edge_mapping, [0, 0, 0, 0]) + + assert work.number_of_edges == 2 + np.testing.assert_array_equal(work.edges(), [0, 2]) + np.testing.assert_allclose( + work.active_edge_values("value")[1], + expected_active, + ) + + +def test_keep_parallel_edges_rejects_self_edge_from_suppression(): + work = bic.graph.ContractionGraph(_triangle(), parallel_edges="keep") + work.suppress_node(0) + + with pytest.raises(ValueError, match="self edge"): + work.suppress_node(1) + + +def test_parallel_edge_policy_validation_and_public_surface(): + with pytest.raises(ValueError, match="parallel_edges"): + bic.graph.ContractionGraph(_triangle(), parallel_edges="invalid") + + assert not hasattr(bic.graph.ContractionGraph, "_from_edges") + assert not hasattr( + bic.graph.ContractionGraph(_triangle()), + "_deleted_original_edges", + ) + + +def test_nan_values_propagate(): + work = bic.graph.ContractionGraph(_triangle()) + work.add_edge_values( + "value", + np.array([0.0, np.nan, 1.0], dtype=np.float64), + reduction="min", + ) + work.contract_edge(0, keep_node=0) + + assert np.isnan(work.materialize().edge_values["value"][0]) + + +def test_invalid_value_registration_and_inactive_ids(): + work = bic.graph.ContractionGraph(_triangle()) + + with pytest.raises(TypeError, match="floating dtype"): + work.add_edge_values("value", np.ones(3, dtype=np.int64), reduction="sum") + with pytest.raises(ValueError, match=r"shape\[0\]"): + work.add_edge_values("value", np.ones(2, dtype=np.float32), reduction="sum") + with pytest.raises(ValueError, match="reduction"): + work.add_edge_values("value", np.ones(3, dtype=np.float32), reduction="median") + + work.add_edge_values("value", np.ones(3, dtype=np.float32), reduction="sum") + with pytest.raises(ValueError, match="already exists"): + work.add_edge_values("value", np.ones(3, dtype=np.float32), reduction="sum") + + work.contract_edge(0, keep_node=0) + with pytest.raises(ValueError, match="before the first graph mutation"): + work.add_node_values("late", np.ones(3, dtype=np.float32), reduction="sum") + with pytest.raises(ValueError, match="inactive"): + work.edge_value("value", 0) + with pytest.raises(ValueError, match="degree 2"): + work.suppress_node(2) diff --git a/tests/skeleton/test_postprocessing.py b/tests/skeleton/test_postprocessing.py index f245220..6cb9a41 100644 --- a/tests/skeleton/test_postprocessing.py +++ b/tests/skeleton/test_postprocessing.py @@ -93,6 +93,100 @@ def test_remove_ticks_empty_edges(): assert len(out_edges) == 0 +def test_remove_ticks_rejects_duplicate_edges(): + vertices = np.array([[0.0, 0.0], [1.0, 0.0]]) + edges = np.array([[0, 1], [1, 0]], dtype=np.int64) + + with pytest.raises(ValueError, match="duplicate undirected edges"): + bic.skeleton.remove_ticks(vertices, edges, tick_length=1.0) + + +def test_remove_ticks_prunes_spurs_at_two_branch_points(): + vertices = np.array( + [ + [0.0, 0.0], + [1.0, 0.0], + [2.0, 0.0], + [3.0, 0.0], + [4.0, 0.0], + [1.0, 1.0], + [3.0, 1.0], + ] + ) + edges = np.array( + [[0, 1], [1, 2], [2, 3], [3, 4], [1, 5], [3, 6]], + dtype=np.int64, + ) + + out_vertices, out_edges, _ = bic.skeleton.remove_ticks( + vertices, + edges, + tick_length=1.5, + ) + + assert len(out_vertices) == 5 + assert len(out_edges) == 4 + + +def test_remove_ticks_keeps_pure_cycle(): + vertices = np.array( + [[0.0, 0.0], [0.0, 1.0], [1.0, 1.0], [1.0, 0.0]] + ) + edges = np.array([[0, 1], [1, 2], [2, 3], [3, 0]], dtype=np.int64) + + out_vertices, out_edges, _ = bic.skeleton.remove_ticks( + vertices, + edges, + tick_length=100.0, + ) + + np.testing.assert_array_equal(out_vertices, vertices) + np.testing.assert_array_equal(out_edges, edges) + + +def test_remove_ticks_preserves_parallel_paths_while_pruning_spur(): + vertices = np.array( + [ + [0.0, 0.0], + [4.0, 0.0], + [2.0, 1.0], + [2.0, 0.0], + [2.0, -1.0], + [0.0, 0.5], + ] + ) + edges = np.array( + [[0, 2], [2, 1], [0, 3], [3, 1], [0, 4], [4, 1], [0, 5]], + dtype=np.int64, + ) + + out_vertices, out_edges, _ = bic.skeleton.remove_ticks( + vertices, + edges, + tick_length=1.0, + ) + + assert len(out_vertices) == 5 + assert len(out_edges) == 6 + + +def test_remove_ticks_preserves_radii_dtype(): + vertices = np.array([[0.0, 0.0], [1.0, 0.0], [1.0, 0.5], [2.0, 0.0]]) + edges = np.array([[0, 1], [1, 2], [1, 3]], dtype=np.int64) + radii = np.arange(4, dtype=np.float32) + + out_vertices, out_edges, out_radii = bic.skeleton.remove_ticks( + vertices, + edges, + tick_length=1.0, + radii=radii, + ) + + assert len(out_vertices) == 3 + assert len(out_edges) == 2 + assert out_radii.dtype == np.float32 + + def test_join_close_components_links_collinear_endpoints(): # Two collinear fragments along x with a gap of 2 between node 1 and node 2. vertices = np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 1.0], [0.0, 0.0, 3.0], [0.0, 0.0, 4.0]]) From f59d8db871608b6c3ff719fc5f88f0d5497523e6 Mon Sep 17 00:00:00 2001 From: Constantin Pape Date: Sat, 25 Jul 2026 21:42:17 +0200 Subject: [PATCH 2/2] Update skeleton post-processing functionality --- MIGRATION_GUIDE.md | 6 ++ .../bioimage_cpp/graph/contraction_graph.hxx | 4 + .../graph/detail/contraction_topology.hxx | 11 +++ .../bioimage_cpp/graph/undirected_graph.hxx | 10 +++ src/bindings/graph.cxx | 9 +++ src/bindings/graph_contraction.cxx | 5 ++ src/bioimage_cpp/graph/__init__.py | 4 + src/bioimage_cpp/skeleton/postprocessing.py | 36 +++++---- tests/graph/test_contraction_graph.py | 3 + tests/graph/test_undirected_graph.py | 13 ++++ tests/skeleton/test_postprocessing.py | 78 +++++++++++++++++++ 11 files changed, 165 insertions(+), 14 deletions(-) diff --git a/MIGRATION_GUIDE.md b/MIGRATION_GUIDE.md index 4bd3007..3199ff3 100644 --- a/MIGRATION_GUIDE.md +++ b/MIGRATION_GUIDE.md @@ -460,6 +460,8 @@ Common method/property mapping: | `extractSubgraphFromNodes` | `extract_subgraph_from_nodes` | | `edgesFromNodeList` | `edges_from_node_list` | +`node_degrees()` returns the degree of every node as a `uint64` array. + ### Mutable Graph Contractions `ContractionGraph` supports custom graph postprocessing with stable ids and @@ -502,6 +504,10 @@ reduce when parallel edges merge. The contracted edge becomes internal and has no output value. `suppress_node` removes a degree-2 node and reduces its two incident edge values into the replacement edge. +Use `can_suppress_node(node)` before suppression when the mutable graph can +contain parallel edges. The method also rejects a suppression that would create +a self-edge. + The input must be a simple `UndirectedGraph`. Use `parallel_edges="keep"` to keep parallel edges that mutations create: diff --git a/include/bioimage_cpp/graph/contraction_graph.hxx b/include/bioimage_cpp/graph/contraction_graph.hxx index 0440c35..d740b8d 100644 --- a/include/bioimage_cpp/graph/contraction_graph.hxx +++ b/include/bioimage_cpp/graph/contraction_graph.hxx @@ -273,6 +273,10 @@ public: return topology_.degree(node); } + [[nodiscard]] bool can_suppress_node(const NodeId node) const { + return topology_.can_suppress_node(node); + } + [[nodiscard]] bool node_active(const NodeId node) const { return topology_.node_active(node); } diff --git a/include/bioimage_cpp/graph/detail/contraction_topology.hxx b/include/bioimage_cpp/graph/detail/contraction_topology.hxx index 2047a95..69acc09 100644 --- a/include/bioimage_cpp/graph/detail/contraction_topology.hxx +++ b/include/bioimage_cpp/graph/detail/contraction_topology.hxx @@ -149,6 +149,17 @@ public: return adjacency_[static_cast(node)].size(); } + [[nodiscard]] bool can_suppress_node(const NodeId node) const { + validate_node_id(node); + if (!nodes_active_[static_cast(node)]) { + return false; + } + const auto &incident = adjacency_[static_cast(node)]; + return incident.size() == 2 && + incident[0].edge != incident[1].edge && + incident[0].node != incident[1].node; + } + [[nodiscard]] const std::vector &node_adjacency( const NodeId node ) const { diff --git a/include/bioimage_cpp/graph/undirected_graph.hxx b/include/bioimage_cpp/graph/undirected_graph.hxx index 715cccd..40fe000 100644 --- a/include/bioimage_cpp/graph/undirected_graph.hxx +++ b/include/bioimage_cpp/graph/undirected_graph.hxx @@ -166,6 +166,16 @@ public: ); } + [[nodiscard]] std::vector node_degrees() const { + ensure_adjacency_built(); + std::vector result(static_cast(number_of_nodes_)); + for (NodeId node = 0; node < number_of_nodes_; ++node) { + const auto index = static_cast(node); + result[index] = adjacency_offsets_[index + 1] - adjacency_offsets_[index]; + } + return result; + } + virtual EdgeId insert_edge(const NodeId u, const NodeId v) { validate_node(u); validate_node(v); diff --git a/src/bindings/graph.cxx b/src/bindings/graph.cxx index 76f2764..3be5510 100644 --- a/src/bindings/graph.cxx +++ b/src/bindings/graph.cxx @@ -546,6 +546,10 @@ UInt64Array graph_node_adjacency(const Graph &graph, const std::uint64_t node) { return result; } +UInt64Array graph_node_degrees(const Graph &graph) { + return detail::copy_vector_to_array(graph.node_degrees()); +} + UInt64Array graph_serialize(const Graph &graph) { auto result = make_uint64_array({static_cast(graph.serialization_size())}); auto *data = result.data(); @@ -1895,6 +1899,11 @@ void bind_graph(nb::module_ &m) { .def("insert_edges", &graph_insert_edges, nb::arg("uvs")) .def("find_edges", &graph_find_edges, nb::arg("uvs")) .def("node_adjacency", &graph_node_adjacency, nb::arg("node")) + .def( + "node_degrees", + &graph_node_degrees, + "Return the degree of each node as a uint64 array." + ) .def_prop_ro("serialization_size", &Graph::serialization_size) .def("serialize", &graph_serialize) .def( diff --git a/src/bindings/graph_contraction.cxx b/src/bindings/graph_contraction.cxx index 9bceea4..16480e6 100644 --- a/src/bindings/graph_contraction.cxx +++ b/src/bindings/graph_contraction.cxx @@ -243,6 +243,11 @@ void bind_graph_contraction(nb::module_ &m) { nb::arg("node") ) .def("degree", &ContractionGraph::degree, nb::arg("node")) + .def( + "can_suppress_node", + &ContractionGraph::can_suppress_node, + nb::arg("node") + ) .def("is_node_active", &ContractionGraph::node_active, nb::arg("node")) .def("is_edge_active", &ContractionGraph::edge_active, nb::arg("edge")) .def( diff --git a/src/bioimage_cpp/graph/__init__.py b/src/bioimage_cpp/graph/__init__.py index 674e2aa..b070d7f 100644 --- a/src/bioimage_cpp/graph/__init__.py +++ b/src/bioimage_cpp/graph/__init__.py @@ -267,6 +267,10 @@ def node_adjacency(self, node: int) -> np.ndarray: def degree(self, node: int) -> int: return int(self._core.degree(int(node))) + def can_suppress_node(self, node: int) -> bool: + """Return whether :meth:`suppress_node` accepts this node.""" + return bool(self._core.can_suppress_node(int(node))) + def is_node_active(self, node: int) -> bool: return bool(self._core.is_node_active(int(node))) diff --git a/src/bioimage_cpp/skeleton/postprocessing.py b/src/bioimage_cpp/skeleton/postprocessing.py index d5b0042..fe35ee7 100644 --- a/src/bioimage_cpp/skeleton/postprocessing.py +++ b/src/bioimage_cpp/skeleton/postprocessing.py @@ -1,9 +1,7 @@ """Post-processing for skeleton graphs. -Resolve junction nodes so each filament becomes its own connected component. -Degree-3 and degree-4 nodes are split or pruned based on the angles between their -incident edges: the straightest pair is kept as the through-going filament and -the remaining arm(s) are either separated, or for short dead-ends (spurs), pruned. +Cut degree-3 branches and split degree-4 crossings based on the angles between +their incident edges. Short dead-end branches can also be pruned by length. """ from __future__ import annotations @@ -92,14 +90,17 @@ def clean_filament_graph( min_join_angle: float = 175.0, save_intermediates: list | None = None, ) -> tuple[np.ndarray, np.ndarray, np.ndarray | None]: - """Split skeleton junctions so each filament is its own component. + """Clean a skeleton graph by cutting branches and splitting crossings. Postprocessing steps, in order: 1. If ``tick_length > 0``, prune dead-end branches shorter than it via :func:`remove_ticks`. - 2. Split each degree-3 junction, separating the odd arm when it diverges - from the through pair by at least ``min_branch_angle``. + 2. At each degree-3 junction, identify the straightest pair of arms. If the + remaining arm diverges by at least ``min_branch_angle``, remove its edge + adjacent to the junction. This cuts the arm from the through-going + filament and removes its first segment. A one-edge arm is removed + completely. 3. Split each degree-4 crossing, separating its two through pairs when they are collinear to within ``min_through_angle``. 4. If ``join_dist > 0``, join collinear endpoints across gaps up to @@ -118,8 +119,9 @@ def clean_filament_graph( min_through_angle: Minimum through-pair angle (degrees) for a degree-4 crossing to split. min_branch_angle: - Minimum angle (degrees) between a degree-3 node's odd arm and its - through pair for the odd arm to be separated. + Minimum angle in degrees between a degree-3 node's odd arm and its + through pair. When the angle meets this threshold, the edge adjacent + to the junction on the odd arm is removed. tick_length: If > 0, prune dead-end branches shorter than this (physical) distance. join_dist: @@ -153,7 +155,7 @@ def _snapshot(name): edges = edges.copy() _snapshot("ticks") graph = skeleton_to_graph(vertices, edges) - degrees = np.bincount(edges.reshape(-1), minlength=len(vertices)) + degrees = graph.node_degrees() splits, prune_edges = [], set() for v in np.where(degrees == 3)[0]: @@ -301,7 +303,7 @@ def remove_ticks(vertices, edges, tick_length, radii=None): eligible = _nodes_in_noncycle_components(indptr, dst, degrees) for node in np.where((degrees == 2) & eligible)[0]: node = int(node) - if work.is_node_active(node) and work.degree(node) == 2: + if work.can_suppress_node(node): work.suppress_node(node) threshold = float(tick_length) @@ -323,7 +325,7 @@ def remove_ticks(vertices, edges, tick_length, radii=None): _, edge, a, b = best work.erase_edge(edge) for node in (a, b): - if work.is_node_active(node) and work.degree(node) == 2: + if work.can_suppress_node(node): work.suppress_node(node) edges = edges[work.materialize().edge_mapping >= 0] @@ -378,7 +380,12 @@ def join_close_components(vertices, edges, dist, *, min_join_angle=175.0, vertices, edges, radii: ``vertices`` and ``radii`` unchanged; ``edges`` has the join edges added. """ - from scipy.spatial import cKDTree + try: + from scipy.spatial import cKDTree + except ImportError as error: + raise ImportError( + "join_close_components requires scipy; install scipy to use this function" + ) from error vertices = np.asarray(vertices, dtype=np.float64) edges = np.asarray(edges, dtype=np.int64) @@ -450,7 +457,8 @@ def draw_instances(vertices, edges, labels, shape, radius=1): Integer array with shape ``(E, 2)`` indexing ``vertices``. labels: Per-vertex integer labels, e.g. from - :func:`bioimage_cpp.graph.connected_components`. + :func:`bioimage_cpp.graph.connected_components`. Both endpoints of + each edge should have the same label. shape: Output volume shape ``(Z, Y, X)``. radius: diff --git a/tests/graph/test_contraction_graph.py b/tests/graph/test_contraction_graph.py index 54dd73a..e02eb5a 100644 --- a/tests/graph/test_contraction_graph.py +++ b/tests/graph/test_contraction_graph.py @@ -246,8 +246,11 @@ def test_keep_parallel_edges_fold_only_in_materialized_snapshot( def test_keep_parallel_edges_rejects_self_edge_from_suppression(): work = bic.graph.ContractionGraph(_triangle(), parallel_edges="keep") + assert work.can_suppress_node(0) work.suppress_node(0) + assert not work.can_suppress_node(0) + assert not work.can_suppress_node(1) with pytest.raises(ValueError, match="self edge"): work.suppress_node(1) diff --git a/tests/graph/test_undirected_graph.py b/tests/graph/test_undirected_graph.py index 51286ee..512e89d 100644 --- a/tests/graph/test_undirected_graph.py +++ b/tests/graph/test_undirected_graph.py @@ -55,6 +55,19 @@ def test_undirected_graph_node_adjacency(): ) +def test_undirected_graph_node_degrees_updates_after_insertion(): + graph = bic.graph.UndirectedGraph.from_edges(5, [[0, 1], [1, 2], [1, 3]]) + + degrees = graph.node_degrees() + + assert degrees.dtype == np.uint64 + np.testing.assert_array_equal(degrees, [1, 3, 1, 1, 0]) + + graph.insert_edge(3, 4) + + np.testing.assert_array_equal(graph.node_degrees(), [1, 3, 1, 2, 1]) + + def test_undirected_graph_serialize_and_deserialize(): graph = bic.graph.UndirectedGraph.from_edges(3, [[0, 2], [0, 1]]) diff --git a/tests/skeleton/test_postprocessing.py b/tests/skeleton/test_postprocessing.py index 6cb9a41..7f2f970 100644 --- a/tests/skeleton/test_postprocessing.py +++ b/tests/skeleton/test_postprocessing.py @@ -1,3 +1,5 @@ +import builtins + import numpy as np import pytest @@ -170,6 +172,31 @@ def test_remove_ticks_preserves_parallel_paths_while_pruning_spur(): assert len(out_edges) == 6 +def test_remove_ticks_preserves_two_parallel_paths_while_pruning_spur(): + vertices = np.array( + [ + [0.0, 0.0], + [4.0, 0.0], + [2.0, 1.0], + [2.0, -1.0], + [0.0, 0.5], + ] + ) + edges = np.array( + [[0, 2], [2, 1], [0, 3], [3, 1], [0, 4]], + dtype=np.int64, + ) + + out_vertices, out_edges, _ = bic.skeleton.remove_ticks( + vertices, + edges, + tick_length=1.0, + ) + + assert len(out_vertices) == 4 + assert len(out_edges) == 4 + + def test_remove_ticks_preserves_radii_dtype(): vertices = np.array([[0.0, 0.0], [1.0, 0.0], [1.0, 0.5], [2.0, 0.0]]) edges = np.array([[0, 1], [1, 2], [1, 3]], dtype=np.int64) @@ -226,6 +253,24 @@ def test_join_close_components_skips_same_component(): assert len(out_edges) == 2 +def test_join_close_components_reports_missing_scipy(monkeypatch): + original_import = builtins.__import__ + + def import_without_scipy(name, *args, **kwargs): + if name == "scipy.spatial": + raise ModuleNotFoundError("No module named 'scipy'") + return original_import(name, *args, **kwargs) + + monkeypatch.setattr(builtins, "__import__", import_without_scipy) + + with pytest.raises(ImportError, match="requires scipy"): + bic.skeleton.join_close_components( + np.array([[0.0, 0.0], [1.0, 0.0]]), + np.array([[0, 1]], dtype=np.int64), + dist=1.0, + ) + + def test_clean_graph_splits_crossing(): vertices = np.array( [[0.0, 0.0, 0.0], [0.0, 0.0, 1.0], [0.0, 0.0, -1.0], [0.0, 1.0, 0.0], [0.0, -1.0, 0.0]], @@ -243,6 +288,39 @@ def test_clean_graph_splits_crossing(): assert len(out_radii) == len(out_vertices) +def test_clean_graph_cuts_degree3_junction_edge(): + vertices = np.array( + [[0.0, 0.0], [1.0, 0.0], [-1.0, 0.0], [0.0, 1.0], [0.0, 2.0]] + ) + edges = np.array([[0, 1], [0, 2], [0, 3], [3, 4]], dtype=np.int64) + + out_vertices, out_edges, _ = bic.skeleton.clean_filament_graph( + vertices, + edges, + direction_span=1, + ) + + np.testing.assert_array_equal(out_vertices, vertices) + np.testing.assert_array_equal(out_edges, [[0, 1], [0, 2], [3, 4]]) + assert _component_count(out_vertices, out_edges) == 2 + + +def test_clean_graph_removes_one_edge_degree3_arm(): + vertices = np.array( + [[0.0, 0.0], [1.0, 0.0], [-1.0, 0.0], [0.0, 1.0]] + ) + edges = np.array([[0, 1], [0, 2], [0, 3]], dtype=np.int64) + + out_vertices, out_edges, _ = bic.skeleton.clean_filament_graph( + vertices, + edges, + direction_span=1, + ) + + np.testing.assert_array_equal(out_vertices, vertices[:3]) + np.testing.assert_array_equal(out_edges, [[0, 1], [0, 2]]) + + def test_draw_instances_labels_edge(): vertices = np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 3.0]]) edges = np.array([[0, 1]], dtype=np.int64)