Implement skeleton to graph conversion

Author: constantinpape

MIGRATION_GUIDE.md

@@ -2551,6 +2551,27 @@ skeletons = bic.skeleton.teasar_labels(
 # {original_label: (vertices, edges, radii), ...}
 ```
 
+Convert the topology of either result to the native undirected graph while
+keeping coordinates and radii in their NumPy arrays:
+
+```python
+import numpy as np
+
+graph = bic.skeleton.skeleton_to_graph(vertices, edges)
+
+# Node ids index vertices and radii directly.
+degrees = np.fromiter(
+    (len(graph.node_adjacency(node)) for node in range(graph.number_of_nodes)),
+    dtype=np.uint64,
+    count=graph.number_of_nodes,
+)
+endpoint_vertices = vertices[degrees <= 1]
+branch_vertices = vertices[degrees > 2]
+```
+
+The conversion also preserves disconnected forests and isolated vertices. The
+graph stores topology only; it does not copy or own ``vertices`` or ``radii``.
+
 Important differences and current scope:
 
 - `teasar(mask)` treats every nonzero value as one foreground class. It

examples/skeleton/teasar.py

@@ -5,17 +5,17 @@
 import napari
 import numpy as np
 
-from bioimage_cpp.skeleton import teasar
+from bioimage_cpp.skeleton import skeleton_to_graph, teasar
 
 
 def main():
     path = Path(__file__).with_name("00004_gt_mask.mrc")
 
-    full_crop = True
+    full_crop = False
     if full_crop:
         bb = np.s_[:]
     else:
-        bb = np.s_[100:200, 1024:1700, 1024:1700]
+        bb = np.s_[100:300, 1024:2048, 1024:2048]
 
     print("Load data ...")
     with mrcfile.mmap(path, mode="r", permissive=True) as mrc:
@@ -32,9 +32,17 @@ def main():
     vertices, edges, _ = teasar(mask, scale=3.0, number_of_threads=8)
     print("... done in:", time.time() - t0, "s")
 
-    # degrees = np.bincount(edges.ravel(), minlength=len(vertices))
-    # endpoints = degrees <= 1
-    # branch_points = degrees > 2
+    graph = skeleton_to_graph(vertices, edges)
+    degrees = np.fromiter(
+        (
+            len(graph.node_adjacency(node))
+            for node in range(graph.number_of_nodes)
+        ),
+        dtype=np.uint64,
+        count=graph.number_of_nodes,
+    )
+    endpoints = degrees <= 1
+    branch_points = degrees > 2
 
     viewer = napari.Viewer(ndisplay=3)
     viewer.add_labels(mask, name="mask", opacity=0.6)
@@ -47,26 +55,26 @@ def main():
         opacity=0.9,
         blending="translucent_no_depth",
     )
-    # viewer.add_points(
-    #     vertices[endpoints],
-    #     name="endpoints",
-    #     size=3,
-    #     face_color="#ffd166",
-    #     border_color="#1b1b1b",
-    #     border_width=0.15,
-    #     n_dimensional=True,
-    #     blending="translucent_no_depth",
-    # )
-    # viewer.add_points(
-    #     vertices[branch_points],
-    #     name="branch points",
-    #     size=5,
-    #     face_color="#ff4d6d",
-    #     border_color="#1b1b1b",
-    #     border_width=0.15,
-    #     n_dimensional=True,
-    #     blending="translucent_no_depth",
-    # )
+    viewer.add_points(
+        vertices[endpoints],
+        name="endpoints",
+        size=3,
+        face_color="#ffd166",
+        border_color="#1b1b1b",
+        border_width=0.15,
+        n_dimensional=True,
+        blending="translucent_no_depth",
+    )
+    viewer.add_points(
+        vertices[branch_points],
+        name="branch points",
+        size=5,
+        face_color="#ff4d6d",
+        border_color="#1b1b1b",
+        border_width=0.15,
+        n_dimensional=True,
+        blending="translucent_no_depth",
+    )
     napari.run()
 
 

src/bioimage_cpp/skeleton/__init__.py

@@ -9,6 +9,7 @@
 from .. import _core
 from .._validation import strict_index
 from ..distance._distance import _as_binary_input, _normalize_sampling, _normalize_threads
+from ._graph import skeleton_to_graph
 
 
 _TEASAR_LABELS_BY_DTYPE = {
@@ -108,6 +109,9 @@ def teasar(
         shape ``(V,)`` and contains physical distance-to-boundary values.
         With several components the tuple represents their skeleton forest in
         first-voxel C-order. An empty mask returns correctly typed empty arrays.
+        Use :func:`skeleton_to_graph` to convert the vertex and edge arrays to
+        an undirected graph; graph node ids continue to index ``vertices`` and
+        ``radii``.
 
     Raises
     ------
@@ -198,4 +202,4 @@ def teasar_labels(
     return run(labels_c, labels_array.dtype.type(background_value), *options)
 
 
-__all__ = ["teasar", "teasar_labels"]
+__all__ = ["skeleton_to_graph", "teasar", "teasar_labels"]

src/bioimage_cpp/skeleton/_graph.py

@@ -0,0 +1,57 @@
+"""Utilities for working with skeleton topology."""
+
+from __future__ import annotations
+
+import numpy as np
+
+from ..graph import UndirectedGraph
+
+
+def skeleton_to_graph(
+    vertices: np.ndarray,
+    edges: np.ndarray,
+) -> UndirectedGraph:
+    """Convert skeleton vertex and edge arrays to an undirected graph.
+
+    This accepts the ``vertices`` and ``edges`` arrays returned by
+    :func:`bioimage_cpp.skeleton.teasar`, or by an individual result from
+    :func:`bioimage_cpp.skeleton.teasar_labels`. Graph node ids correspond
+    directly to rows in ``vertices`` (and in the accompanying TEASAR
+    ``radii`` array). The returned graph stores topology only; coordinates and
+    radii remain in their original arrays.
+
+    Parameters
+    ----------
+    vertices:
+        Two-dimensional array with one row per skeleton vertex. Only the
+        number of rows is used to construct the graph.
+    edges:
+        Integer array with shape ``(E, 2)`` containing vertex ids. Edge
+        endpoints must be smaller than ``len(vertices)`` and self-edges are
+        not supported.
+
+    Returns
+    -------
+    bioimage_cpp.graph.UndirectedGraph
+        An undirected graph with one node per vertex and the topology given by
+        ``edges``. Empty skeletons and isolated vertices are preserved.
+
+    Raises
+    ------
+    ValueError
+        If ``vertices`` is not two-dimensional or the edge shape/topology is
+        invalid.
+    TypeError
+        If ``edges`` is not an integer array.
+    IndexError
+        If an edge endpoint is outside the vertex range.
+    """
+    vertex_array = np.asarray(vertices)
+    if vertex_array.ndim != 2:
+        raise ValueError(
+            f"vertices must be a 2D array, got ndim={vertex_array.ndim}"
+        )
+    return UndirectedGraph.from_edges(vertex_array.shape[0], edges)
+
+
+__all__ = ["skeleton_to_graph"]

tests/skeleton/test_skeleton_to_graph.py

@@ -0,0 +1,69 @@
+import numpy as np
+import pytest
+
+import bioimage_cpp as bic
+
+
+def _node_degrees(graph):
+    return np.fromiter(
+        (
+            len(graph.node_adjacency(node))
+            for node in range(graph.number_of_nodes)
+        ),
+        dtype=np.uint64,
+        count=graph.number_of_nodes,
+    )
+
+
+def test_skeleton_to_graph_preserves_teasar_topology_and_vertex_ids():
+    mask = np.zeros((9, 11, 10), dtype=np.uint8)
+    mask[4, 5, 1:5] = 1
+    for step in range(5):
+        mask[4, 5 - step, 4 + step] = 1
+        mask[4, 5 + step, 4 + step] = 1
+
+    vertices, edges, _ = bic.skeleton.teasar(mask)
+    graph = bic.skeleton.skeleton_to_graph(vertices, edges)
+
+    assert isinstance(graph, bic.graph.UndirectedGraph)
+    assert graph.number_of_nodes == len(vertices)
+    assert graph.number_of_edges == len(edges)
+    np.testing.assert_array_equal(graph.uv_ids(), edges)
+
+    degrees = _node_degrees(graph)
+    endpoints = {
+        tuple(vertex.astype(int)) for vertex in vertices[degrees <= 1]
+    }
+    assert endpoints == {(4, 5, 1), (4, 1, 8), (4, 9, 8)}
+    assert np.count_nonzero(degrees > 2) == 1
+
+
+def test_skeleton_to_graph_preserves_empty_and_isolated_vertices():
+    empty = bic.skeleton.skeleton_to_graph(
+        np.empty((0, 3), dtype=np.float64),
+        np.empty((0, 2), dtype=np.uint64),
+    )
+    assert empty.number_of_nodes == 0
+    assert empty.number_of_edges == 0
+
+    isolated = bic.skeleton.skeleton_to_graph(
+        np.array([[2.0, 3.0, 4.0]]),
+        np.empty((0, 2), dtype=np.uint64),
+    )
+    assert isolated.number_of_nodes == 1
+    assert isolated.number_of_edges == 0
+    np.testing.assert_array_equal(_node_degrees(isolated), [0])
+
+
+def test_skeleton_to_graph_rejects_invalid_inputs():
+    with pytest.raises(ValueError, match="vertices must be a 2D array"):
+        bic.skeleton.skeleton_to_graph(
+            np.array([0.0, 1.0, 2.0]),
+            np.empty((0, 2), dtype=np.uint64),
+        )
+
+    vertices = np.zeros((2, 3), dtype=np.float64)
+    with pytest.raises(ValueError, match="uvs must have shape"):
+        bic.skeleton.skeleton_to_graph(vertices, [0, 1])
+    with pytest.raises(IndexError, match="node id must be < number_of_nodes"):
+        bic.skeleton.skeleton_to_graph(vertices, [[0, 2]])