From 662f5b42e23d15ffa149017c082b5d1e53797e7e Mon Sep 17 00:00:00 2001 From: GOESTERN-1107078 Date: Tue, 21 Jul 2026 17:19:20 +0200 Subject: [PATCH 1/3] instance segmentation pipeline - graph postprocessing, unit tests, example script --- examples/skeleton/instance_segmentation.py | 127 +++++ src/bioimage_cpp/skeleton/__init__.py | 3 +- src/bioimage_cpp/skeleton/postprocessing.py | 500 ++++++++++++++++++++ tests/skeleton/test_postprocessing.py | 196 ++++++++ 4 files changed, 825 insertions(+), 1 deletion(-) create mode 100644 examples/skeleton/instance_segmentation.py create mode 100644 src/bioimage_cpp/skeleton/postprocessing.py create mode 100644 tests/skeleton/test_postprocessing.py diff --git a/examples/skeleton/instance_segmentation.py b/examples/skeleton/instance_segmentation.py new file mode 100644 index 0000000..cf06d81 --- /dev/null +++ b/examples/skeleton/instance_segmentation.py @@ -0,0 +1,127 @@ +#!/usr/bin/env python +"""Split a skeletonized mask into filament instances with bioimage_cpp. + +Skeletonize a binary mask with TEASAR, clean the skeleton graph so each filament +becomes its own connected component, and save the labeled instances. +""" + +import argparse +from pathlib import Path + +import numpy as np + +from bioimage_cpp.graph import connected_components +from bioimage_cpp.skeleton import clean_graph, draw_instances, skeleton_to_graph, teasar + + +def parse_args(): + parser = argparse.ArgumentParser(description=__doc__.splitlines()[0]) + parser.add_argument("--mask_path", required=True, + help="Binary mask to skeletonize; an .h5/.hdf5 or MRC volume.") + parser.add_argument("--seg_key", default=None, + help="HDF5 dataset holding the mask; required for .h5/.hdf5 inputs.") + parser.add_argument("--output_dir", default=".", + help="Directory for the skeleton and instance .npz files.") + parser.add_argument("--tag", default="", + help="Suffix appended to the output filenames.") + parser.add_argument("--pixel_size", type=float, default=10.0, + help="Isotropic voxel spacing for teasar; graph coordinates use these units.") + parser.add_argument("--scale", type=float, default=0.0, + help="teasar invalidation-radius scale.") + parser.add_argument("--constant", type=float, default=70.0, + help="teasar invalidation-radius constant.") + parser.add_argument("--number_of_threads", type=int, default=8, + help="Threads for teasar.") + parser.add_argument("--direction_span", type=int, default=10, + help="Nodes walked along each arm to estimate its direction at a junction.") + parser.add_argument("--min_through_angle", type=float, default=170.0, + help="Min through-pair angle (deg) for a degree-4 crossing to split.") + parser.add_argument("--min_branch_angle", type=float, default=30.0, + help="Min branch angle (deg) for a degree-3 odd arm to be separated.") + parser.add_argument("--tick_length", type=float, default=50.0, + help="Prune dead-end branches shorter than this; 0 disables.") + parser.add_argument("--join_radius", type=float, default=50.0, + help="Join collinear endpoints across gaps up to this distance; 0 disables.") + parser.add_argument("--min_join_angle", type=float, default=175.0, + help="Min straightness (deg) for a join; 180 is collinear.") + parser.add_argument("--circle_size", type=float, default=40.0, + help="Instance tube diameter for the --view render, in --pixel_size units.") + parser.add_argument("--view", action="store_true", + help="Open the pipeline stages in napari.") + return parser.parse_args() + + +def load_mask(mask_path, seg_key): + path = Path(mask_path) + if path.suffix in (".h5", ".hdf5"): + import h5py + + if seg_key is None: + raise ValueError("--seg_key is required for HDF5 inputs") + with h5py.File(path, "r") as f: + if seg_key not in f: + raise KeyError(f"{seg_key} not in {path}") + return np.asarray(f[seg_key][:]) + + import mrcfile + + with mrcfile.mmap(path, mode="r", permissive=True) as mrc: + return np.asarray(mrc.data) + + +def view_stages(binary, stages, pixel_size, circle_size): + import napari + + radius = (circle_size / 2) / pixel_size + viewer = napari.Viewer(ndisplay=3) + viewer.add_labels(binary, name="mask", opacity=0.4) + for index, (name, vertices, edges, _) in enumerate(stages): + labels = connected_components(skeleton_to_graph(vertices, edges)) + volume = draw_instances(vertices / pixel_size, edges, labels, binary.shape, radius) + viewer.add_labels(volume, name=f"instances ({name})", visible=index == len(stages) - 1) + napari.run() + + +def main(): + args = parse_args() + mask = load_mask(args.mask_path, args.seg_key) + binary = (mask > 0).astype(np.uint8) + + spacing = (args.pixel_size,) * 3 + raw_vertices, raw_edges, raw_radii = teasar( + binary, spacing=spacing, scale=args.scale, constant=args.constant, + number_of_threads=args.number_of_threads, + ) + + stages = [] + vertices, edges, radii = clean_graph( + raw_vertices, raw_edges, radii=raw_radii, + direction_span=args.direction_span, min_through_angle=args.min_through_angle, + min_branch_angle=args.min_branch_angle, tick_length=args.tick_length, + join_radius=args.join_radius, min_join_angle=args.min_join_angle, + save_intermediates=stages, + ) + labels = connected_components(skeleton_to_graph(vertices, edges)) + + name = Path(args.mask_path).stem + if args.tag: + name = f"{name}_{args.tag}" + output_dir = Path(args.output_dir) + output_dir.mkdir(parents=True, exist_ok=True) + + skeleton = dict(vertices=raw_vertices, edges=raw_edges) + if raw_radii is not None: + skeleton["radii"] = raw_radii + instances = dict(vertices=vertices, edges=edges, labels=labels) + if radii is not None: + instances["radii"] = radii + np.savez_compressed(output_dir / f"{name}_skeleton.npz", **skeleton) + np.savez_compressed(output_dir / f"{name}_instances.npz", **instances) + print(f"{name}: {len(np.unique(labels))} instances -> {output_dir}") + + if args.view: + view_stages(binary, stages, args.pixel_size, args.circle_size) + + +if __name__ == "__main__": + main() diff --git a/src/bioimage_cpp/skeleton/__init__.py b/src/bioimage_cpp/skeleton/__init__.py index ef5aae9..dee457a 100644 --- a/src/bioimage_cpp/skeleton/__init__.py +++ b/src/bioimage_cpp/skeleton/__init__.py @@ -10,6 +10,7 @@ from .._validation import strict_index from ..distance._distance import _as_binary_input, _normalize_sampling, _normalize_threads from ._graph import skeleton_to_graph +from .postprocessing import clean_graph, draw_instances, join_close_components, remove_ticks _TEASAR_LABELS_BY_DTYPE = { @@ -202,4 +203,4 @@ def teasar_labels( return run(labels_c, labels_array.dtype.type(background_value), *options) -__all__ = ["skeleton_to_graph", "teasar", "teasar_labels"] +__all__ = ["clean_graph", "draw_instances", "join_close_components", "remove_ticks", "skeleton_to_graph", "teasar", "teasar_labels"] diff --git a/src/bioimage_cpp/skeleton/postprocessing.py b/src/bioimage_cpp/skeleton/postprocessing.py new file mode 100644 index 0000000..4fd9a8c --- /dev/null +++ b/src/bioimage_cpp/skeleton/postprocessing.py @@ -0,0 +1,500 @@ +"""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), removed. +""" + +from __future__ import annotations + +from collections import defaultdict + +import numpy as np + +from ..graph import connected_components +from ._graph import skeleton_to_graph + + +def _tangent(v, first_neighbor, graph, vertices, span): + prev, cur = v, int(first_neighbor) + for _ in range(span - 1): + nbrs = np.asarray(graph.node_adjacency(cur))[:, 0] + if len(nbrs) != 2: + break + nxt = int(nbrs[0]) if int(nbrs[0]) != prev else int(nbrs[1]) + prev, cur = cur, nxt + d = vertices[cur] - vertices[v] + n = np.linalg.norm(d) + return d / n if n > 0 else d + + +def _pair_angle(di, dj): + return np.degrees(np.arccos(np.clip(di @ dj, -1.0, 1.0))) + + +def split_degree3(v, graph, vertices, direction_span=1, min_branch_angle=30.0): + adj = np.asarray(graph.node_adjacency(int(v))) + if adj.shape[0] != 3: + return None + neighbors, edge_ids = adj[:, 0], adj[:, 1] + dirs = np.stack([_tangent(int(v), n, graph, vertices, direction_span) for n in neighbors]) + best, best_angle = None, -1.0 + for a, b in [(0, 1), (0, 2), (1, 2)]: + ang = _pair_angle(dirs[a], dirs[b]) + if ang > best_angle: + best_angle, best = ang, (a, b) + i, j = best + odd = ({0, 1, 2} - {i, j}).pop() + branch_angle = min(_pair_angle(dirs[odd], dirs[i]), _pair_angle(dirs[odd], dirs[j])) + if branch_angle < min_branch_angle: + return None + return [int(edge_ids[odd])] + + +def split_degree4(v, graph, vertices, direction_span=1, min_through_angle=160.0): + adj = np.asarray(graph.node_adjacency(int(v))) + if adj.shape[0] != 4: + return None + neighbors, edge_ids = adj[:, 0], adj[:, 1] + dirs = np.stack([_tangent(int(v), n, graph, vertices, direction_span) for n in neighbors]) + best, best_score, best_min = None, -np.inf, 0.0 + for (a, b), (c, d) in [((0, 1), (2, 3)), ((0, 2), (1, 3)), ((0, 3), (1, 2))]: + ang1, ang2 = _pair_angle(dirs[a], dirs[b]), _pair_angle(dirs[c], dirs[d]) + if ang1 + ang2 > best_score: + best_score, best, best_min = ang1 + ang2, ((a, b), (c, d)), min(ang1, ang2) + if best_min < min_through_angle: + return None + (_, pair_b) = best + return [int(edge_ids[k]) for k in pair_b] + + +def clean_graph( + vertices: np.ndarray, + edges: np.ndarray, + radii: np.ndarray | None = None, + *, + direction_span: int = 5, + min_through_angle: float = 160.0, + min_branch_angle: float = 30.0, + tick_length: float = 0.0, + join_radius: float = 0.0, + 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. + + 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``. + 3. Split each degree-4 crossing, separating its two through pairs when they + are collinear to within ``min_through_angle``. + 4. If ``join_radius > 0``, reconnect collinear endpoints across gaps less + than this distance via :func:`join_close_components`. + + Parameters + ---------- + vertices: + Float array with shape ``(V, D)`` of skeleton vertex coordinates. + edges: + Integer array with shape ``(E, 2)`` indexing ``vertices``. + radii: + Optional per-vertex radii, carried through the same remapping. + direction_span: + Number of nodes over which each arm's tangent is measured. + 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. + tick_length: + If > 0, prune dead-end branches shorter than this (physical) distance. + join_radius: + If > 0, reconnect collinear endpoints across gaps up to this (physical) + distance. + min_join_angle: + Minimum straightness (degrees) required for a join; 180 is collinear. + save_intermediates: + If a list is given, ``(name, vertices, edges, radii)`` snapshots are + appended after each step ("raw", "ticks", "split", "join"). + + Returns + ------- + vertices, edges, radii: + The cleaned arrays, reindexed with unused vertices dropped. ``radii`` is + ``None`` when no input radii were given. + """ + vertices = np.asarray(vertices, dtype=np.float64) + edges = np.asarray(edges, dtype=np.int64).copy() + + def _snapshot(name): + if save_intermediates is not None: + save_intermediates.append(( + name, vertices.copy(), edges.copy(), + None if radii is None else np.asarray(radii).copy(), + )) + + _snapshot("raw") + if tick_length > 0: + vertices, edges, radii = remove_ticks(vertices, edges, tick_length, radii=radii) + edges = edges.copy() + _snapshot("ticks") + graph = skeleton_to_graph(vertices, edges) + degrees = np.bincount(edges.reshape(-1), minlength=len(vertices)) + + splits, prune_edges = [], set() + for v in np.where(degrees == 3)[0]: + ids = split_degree3(v, graph, vertices, direction_span, min_branch_angle) + if ids: + prune_edges.update(ids) + for v in np.where(degrees == 4)[0]: + ids = split_degree4(v, graph, vertices, direction_span, min_through_angle) + if ids: + splits.append((int(v), ids)) + + extra_vertices, extra_radii = [], [] + next_id = len(vertices) + for node, edge_ids in splits: + dup = next_id + next_id += 1 + extra_vertices.append(vertices[node]) + if radii is not None: + extra_radii.append(radii[node]) + for e in edge_ids: + edges[e] = np.where(edges[e] == node, dup, edges[e]) + if extra_vertices: + vertices = np.concatenate([vertices, np.asarray(extra_vertices)], axis=0) + if radii is not None: + radii = np.concatenate([radii, np.asarray(extra_radii)]) + + if prune_edges: + keep = np.ones(len(edges), dtype=bool) + keep[list(prune_edges)] = False + edges = edges[keep] + + used = np.zeros(len(vertices), dtype=bool) + if len(edges): + used[edges.reshape(-1)] = True + remap = np.full(len(vertices), -1, dtype=np.int64) + remap[used] = np.arange(int(used.sum())) + vertices = vertices[used] + if radii is not None: + radii = radii[used] + edges = remap[edges] + _snapshot("split") + + if join_radius > 0: + vertices, edges, radii = join_close_components( + vertices, edges, join_radius, + min_join_angle=min_join_angle, direction_span=direction_span, radii=radii, + ) + _snapshot("join") + return vertices, edges, radii + + +def _adjacency(num_nodes, edges): + src = np.concatenate([edges[:, 0], edges[:, 1]]) + dst = np.concatenate([edges[:, 1], edges[:, 0]]) + eid = np.concatenate([np.arange(len(edges)), np.arange(len(edges))]) + order = np.argsort(src, kind="stable") + dst, eid = dst[order], eid[order] + degrees = np.bincount(src, minlength=num_nodes) + indptr = np.zeros(num_nodes + 1, dtype=np.int64) + np.cumsum(degrees, out=indptr[1:]) + return indptr, dst, eid, degrees + + +def remove_ticks(vertices, edges, tick_length, radii=None): + """Remove short dead-end branches ("ticks") from a skeleton graph. + + Ports kimimaro's `remove_ticks`. A distance graph is built over the critical + points (terminals, degree 1; branch points, degree >= 3), whose superedges + are the paths between them weighted by physical length. The shortest terminal + branch below ``tick_length`` is removed repeatedly; when a branch point drops + to degree 2 its two superedges are fused into one, so a real filament end is + re-measured rather than clipped. Standalone paths (both ends terminal) are + never removed. + + Parameters + ---------- + vertices: + Float array with shape ``(V, D)`` of skeleton vertex coordinates. + edges: + Integer array with shape ``(E, 2)`` indexing ``vertices``. + tick_length: + Maximum branch length (physical) that may be culled. + radii: + Optional per-vertex radii, carried through the same remapping. + + Returns + ------- + vertices, edges, radii: + The pruned arrays, reindexed with unused vertices dropped. ``radii`` is + ``None`` when no input radii were given. + """ + vertices = np.asarray(vertices, dtype=np.float64) + edges = np.asarray(edges, dtype=np.int64) + num_nodes = len(vertices) + 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() + 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 + if best is None: + break + a, b, length, path = supers.pop(best) + incident[a].discard(best) + incident[b].discard(best) + dropped.update(path) + 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() + + used = np.zeros(num_nodes, dtype=bool) + if len(edges): + used[edges.reshape(-1)] = True + remap = np.full(num_nodes, -1, dtype=np.int64) + remap[used] = np.arange(int(used.sum())) + vertices = vertices[used] + if radii is not None: + radii = radii[used] + edges = remap[edges] + return vertices, edges, radii + + +def _endpoint_tangent(endpoint, indptr, dst, degrees, vertices, span): + if degrees[endpoint] == 0: + return None + prev, cur = endpoint, int(dst[indptr[endpoint]]) + for _ in range(span - 1): + if degrees[cur] != 2: + break + s, e = indptr[cur], indptr[cur + 1] + nbrs = dst[s:e] + nxt = int(nbrs[0]) if int(nbrs[0]) != prev else int(nbrs[1]) + prev, cur = cur, nxt + direction = vertices[endpoint] - vertices[cur] # outward: interior -> tip + norm = np.linalg.norm(direction) + return direction / norm if norm > 0 else None + + +def join_close_components(vertices, edges, radius, *, min_join_angle=175.0, + direction_span=5, radii=None): + """Reconnect fragmented filaments by joining collinear endpoints across gaps. + + Endpoints (degree = 1) of different connected components are joined with a + new edge when they are within ``radius`` and the two fragments are nearly + collinear through the gap: the outward tangent at each endpoint must point + along the gap to within ``180 - min_join_angle`` degrees, so a straight + continuation reads ~180 (``min_join_angle``). Joins are made shortest-first + with a union-find so each pair of components is linked at most once and each + endpoint is used once. + + Parameters + ---------- + vertices, edges: + Skeleton graph; ``edges`` indexes ``vertices``. + radius: + Maximum gap (physical) across which endpoints may be joined. + min_join_angle: + Minimum straightness (degrees) of the joined path; 180 is perfectly + collinear. + direction_span: + Nodes over which each endpoint's tangent is measured. + radii: + Optional per-vertex radii, returned unchanged. + + Returns + ------- + vertices, edges, radii: + ``vertices`` and ``radii`` unchanged; ``edges`` has the join edges added. + """ + from scipy.spatial import cKDTree + + vertices = np.asarray(vertices, dtype=np.float64) + edges = np.asarray(edges, dtype=np.int64) + n = len(vertices) + if len(edges) == 0: + return vertices, edges.copy(), radii + + indptr, dst, _, degrees = _adjacency(n, edges) + labels = connected_components(skeleton_to_graph(vertices, edges)) + endpoints = np.where(degrees <= 1)[0] + if len(endpoints) < 2: + return vertices, edges.copy(), radii + + tol = np.deg2rad(180.0 - min_join_angle) + tree = cKDTree(vertices[endpoints]) + candidates = [] + for ia, ib in tree.query_pairs(radius): + a, b = int(endpoints[ia]), int(endpoints[ib]) + if labels[a] == labels[b]: + continue + gap = vertices[b] - vertices[a] + dist = float(np.linalg.norm(gap)) + if dist == 0.0: + continue + ta = _endpoint_tangent(a, indptr, dst, degrees, vertices, direction_span) + tb = _endpoint_tangent(b, indptr, dst, degrees, vertices, direction_span) + if ta is None or tb is None: + continue + gdir = gap / dist + bend_a = np.arccos(np.clip(ta @ gdir, -1.0, 1.0)) + bend_b = np.arccos(np.clip(tb @ -gdir, -1.0, 1.0)) + if bend_a <= tol and bend_b <= tol: + candidates.append((dist, a, b)) + + candidates.sort() + parent = {} + + def find(x): + parent.setdefault(x, x) + root = x + while parent[root] != root: + root = parent[root] + while parent[x] != root: + parent[x], x = root, parent[x] + return root + + used, new_edges = set(), [] + for _, a, b in candidates: + if a in used or b in used: + continue + ra, rb = find(int(labels[a])), find(int(labels[b])) + if ra == rb: + continue + parent[ra] = rb + new_edges.append([a, b]) + used.add(a) + used.add(b) + + if new_edges: + edges = np.concatenate([edges, np.asarray(new_edges, dtype=np.int64)], axis=0) + else: + edges = edges.copy() + return vertices, edges, radii + + +def _line_voxels(start, stop): + delta = stop - start + steps = int(np.abs(delta).max()) + 1 + t = np.linspace(0.0, 1.0, steps) + return np.rint(start[None, :] + t[:, None] * delta[None, :]).astype(np.int64) + + +def draw_instances(vertices, edges, labels, shape, radius=1): + """Rasterize a labeled skeleton graph into a dense instance volume. + + Each edge is drawn as a line between its two vertices and dilated by a ball + of the given radius. Every voxel on a component's tubes is set to that + component's label plus one, so background stays zero. + + Parameters + ---------- + vertices: + Float array with shape ``(V, 3)`` of vertex coordinates in voxel index + space matching ``shape`` (``(z, y, x)`` order). + edges: + Integer array with shape ``(E, 2)`` indexing ``vertices``. + labels: + Per-vertex integer labels, e.g. from + :func:`bioimage_cpp.graph.connected_components`. + shape: + Output volume shape ``(Z, Y, X)``. + radius: + Tube radius in voxels. + + Returns + ------- + numpy.ndarray + Integer volume of ``shape`` with background ``0`` and each tube voxel + set to ``labels[endpoint] + 1``. + """ + vertices = np.asarray(vertices) + edges = np.asarray(edges) + labels = np.asarray(labels) + shape = tuple(int(s) for s in shape) + + if len(edges) == 0: + return np.zeros(shape, dtype=np.uint16) + + r = int(radius) + zz, yy, xx = np.ogrid[-r:r + 1, -r:r + 1, -r:r + 1] + offsets = np.stack(np.where(zz ** 2 + yy ** 2 + xx ** 2 <= r * r), axis=1) - r + + vi = np.rint(vertices).astype(np.int64) + centers, center_labels = [], [] + for a, b in edges: + pts = _line_voxels(vi[a], vi[b]) + centers.append(pts) + center_labels.append(np.full(len(pts), int(labels[a]) + 1)) + centers = np.concatenate(centers) + center_labels = np.concatenate(center_labels) + + dtype = np.uint16 if int(center_labels.max()) < 2 ** 16 else np.uint32 + volume = np.zeros(shape, dtype=dtype) + shape_arr = np.asarray(shape) + for offset in offsets: + coords = centers + offset + inside = ((coords >= 0) & (coords < shape_arr)).all(axis=1) + c = coords[inside] + volume[c[:, 0], c[:, 1], c[:, 2]] = center_labels[inside] + return volume + + +__all__ = ["clean_graph", "draw_instances", "join_close_components", "remove_ticks"] diff --git a/tests/skeleton/test_postprocessing.py b/tests/skeleton/test_postprocessing.py new file mode 100644 index 0000000..42b96d2 --- /dev/null +++ b/tests/skeleton/test_postprocessing.py @@ -0,0 +1,196 @@ +import numpy as np +import pytest + +import bioimage_cpp as bic +from bioimage_cpp.skeleton.postprocessing import split_degree3, split_degree4 + + +def _graph(vertices, edges): + vertices = np.asarray(vertices, dtype=np.float64) + edges = np.asarray(edges, dtype=np.int64) + return vertices, edges, bic.skeleton.skeleton_to_graph(vertices, edges) + + +def _component_count(vertices, edges): + labels = bic.graph.connected_components(bic.skeleton.skeleton_to_graph(vertices, edges)) + return len(np.unique(labels)) + + +def _edge_rows(edges): + return {tuple(sorted(map(int, edge))) for edge in np.asarray(edges)} + + +def test_split_degree3_prunes_the_arm_that_diverges_from_the_through_pair(): + # Node 0 is a T: a straight through-pair along x (nodes 1, 2) plus an odd arm along y (node 3). + vertices, _, graph = _graph( + [[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, 1], [0, 2], [0, 3]], + ) + assert split_degree3(0, graph, vertices, direction_span=1, min_branch_angle=30.0) == [2] + + +def test_split_degree3_keeps_a_nearly_collinear_arm(): + # The third arm now runs almost parallel to the through pair, below min_branch_angle. + vertices, _, graph = _graph( + [[0.0, 0.0, 0.0], [0.0, 0.0, 1.0], [0.0, 0.0, -1.0], [0.0, 1.0, 10.0]], + [[0, 1], [0, 2], [0, 3]], + ) + assert split_degree3(0, graph, vertices, direction_span=1, min_branch_angle=30.0) is None + + +def test_split_degree4_detaches_one_collinear_through_pair(): + # A '+' crossing: edges 0,1 are the x through-pair, edges 2,3 the y through-pair. + vertices, _, graph = _graph( + [[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]], + [[0, 1], [0, 2], [0, 3], [0, 4]], + ) + result = split_degree4(0, graph, vertices, direction_span=1, min_through_angle=170.0) + assert set(result) in ({0, 1}, {2, 3}) + + +def test_split_degree4_keeps_a_crossing_without_a_straight_through_pair(): + # Four arms at 0, 60, 120, 180 degrees: no pair is collinear enough to split. + angle = np.deg2rad([0.0, 60.0, 120.0, 180.0]) + arms = np.stack([np.zeros_like(angle), np.sin(angle), np.cos(angle)], axis=1) + vertices, _, graph = _graph( + np.concatenate([[[0.0, 0.0, 0.0]], arms]), + [[0, 1], [0, 2], [0, 3], [0, 4]], + ) + assert split_degree4(0, graph, vertices, direction_span=1, min_through_angle=170.0) is None + + +@pytest.mark.parametrize( + "tick_length, n_vertices, n_edges", + [(1.5, 5, 4), (0.5, 6, 5)], +) +def test_remove_ticks_prunes_only_spurs_below_the_threshold(tick_length, n_vertices, n_edges): + # Backbone 0-1-2-3-4 along x with a length-1 spur (node 5) off the branch node 2. + vertices = np.array( + [[0.0, 0.0, 0.0], [0.0, 0.0, 1.0], [0.0, 0.0, 2.0], [0.0, 0.0, 3.0], + [0.0, 0.0, 4.0], [0.0, 1.0, 2.0]], + ) + edges = np.array([[0, 1], [1, 2], [2, 3], [3, 4], [2, 5]], dtype=np.int64) + + out_vertices, out_edges, _ = bic.skeleton.remove_ticks(vertices, edges, tick_length=tick_length) + + assert len(out_vertices) == n_vertices + assert len(out_edges) == n_edges + + +def test_remove_ticks_never_removes_a_standalone_path(): + vertices = np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 1.0], [0.0, 0.0, 2.0]]) + edges = np.array([[0, 1], [1, 2]], dtype=np.int64) + + out_vertices, out_edges, _ = bic.skeleton.remove_ticks(vertices, edges, tick_length=100.0) + + assert len(out_vertices) == 3 + assert len(out_edges) == 2 + + +def test_remove_ticks_on_empty_edges_returns_the_inputs(): + vertices = np.zeros((3, 3), dtype=np.float64) + edges = np.empty((0, 2), dtype=np.int64) + + out_vertices, out_edges, _ = bic.skeleton.remove_ticks(vertices, edges, tick_length=10.0) + + assert len(out_vertices) == 3 + assert len(out_edges) == 0 + + +def test_join_close_components_links_collinear_endpoints_within_radius(): + # 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]]) + edges = np.array([[0, 1], [2, 3]], dtype=np.int64) + + _, out_edges, _ = bic.skeleton.join_close_components(vertices, edges, radius=2.5) + + assert len(out_edges) == 3 + assert (1, 2) in _edge_rows(out_edges) + + +def test_join_close_components_ignores_endpoints_beyond_the_radius(): + 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]]) + edges = np.array([[0, 1], [2, 3]], dtype=np.int64) + + _, out_edges, _ = bic.skeleton.join_close_components(vertices, edges, radius=1.5) + + assert len(out_edges) == 2 + + +def test_join_close_components_ignores_non_collinear_endpoints(): + # Two fragments offset in y so the gap bends away from each endpoint tangent. + vertices = np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 1.0], [0.0, 3.0, 2.0], [0.0, 4.0, 2.0]]) + edges = np.array([[0, 1], [2, 3]], dtype=np.int64) + + _, out_edges, _ = bic.skeleton.join_close_components(vertices, edges, radius=3.5) + + assert len(out_edges) == 2 + + +def test_join_close_components_does_not_join_within_one_component(): + vertices = np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 1.0], [0.0, 0.0, 2.0]]) + edges = np.array([[0, 1], [1, 2]], dtype=np.int64) + + _, out_edges, _ = bic.skeleton.join_close_components(vertices, edges, radius=5.0) + + assert len(out_edges) == 2 + + +def test_clean_graph_splits_a_crossing_into_two_components(): + 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]], + ) + edges = np.array([[0, 1], [0, 2], [0, 3], [0, 4]], dtype=np.int64) + radii = np.arange(len(vertices), dtype=np.float64) + + assert _component_count(vertices, edges) == 1 + + out_vertices, out_edges, out_radii = bic.skeleton.clean_graph( + vertices, edges, radii=radii, direction_span=1, tick_length=0.0, join_radius=0.0, + ) + + assert _component_count(out_vertices, out_edges) == 2 + assert len(out_radii) == len(out_vertices) + + +def test_clean_graph_records_every_stage_snapshot(): + vertices = np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 1.0], [0.0, 0.0, 2.0]]) + edges = np.array([[0, 1], [1, 2]], dtype=np.int64) + stages = [] + + bic.skeleton.clean_graph(vertices, edges, save_intermediates=stages) + + assert [name for name, *_ in stages] == ["raw", "ticks", "split", "join"] + + +def test_draw_instances_rasterizes_the_edge_with_the_component_label(): + vertices = np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 3.0]]) + edges = np.array([[0, 1]], dtype=np.int64) + labels = np.array([2, 2]) + + volume = bic.skeleton.draw_instances(vertices, edges, labels, (1, 1, 6), radius=0) + + np.testing.assert_array_equal(volume[0, 0], [3, 3, 3, 3, 0, 0]) + assert set(np.unique(volume)) == {0, 3} + + +def test_draw_instances_dilates_by_the_radius(): + vertices = np.array([[2.0, 2.0, 0.0], [2.0, 2.0, 4.0]]) + edges = np.array([[0, 1]], dtype=np.int64) + labels = np.array([0, 0]) + + thin = bic.skeleton.draw_instances(vertices, edges, labels, (5, 5, 5), radius=0) + thick = bic.skeleton.draw_instances(vertices, edges, labels, (5, 5, 5), radius=1) + + assert int(thick.sum() > 0) and thick.astype(bool).sum() > thin.astype(bool).sum() + + +def test_draw_instances_on_empty_edges_returns_a_zero_volume(): + vertices = np.zeros((2, 3), dtype=np.float64) + edges = np.empty((0, 2), dtype=np.int64) + labels = np.zeros(2, dtype=np.int64) + + volume = bic.skeleton.draw_instances(vertices, edges, labels, (2, 3, 4)) + + assert volume.shape == (2, 3, 4) + assert volume.sum() == 0 From 503ff98675bfed8c4bebd11fd0ced86cd305cf37 Mon Sep 17 00:00:00 2001 From: GOESTERN-1107078 Date: Tue, 21 Jul 2026 18:09:26 +0200 Subject: [PATCH 2/3] skeleton postprocessing - edits to code style and documentation --- examples/skeleton/instance_segmentation.py | 54 +++++------- src/bioimage_cpp/skeleton/postprocessing.py | 94 +++++++++------------ tests/skeleton/test_postprocessing.py | 58 +++++-------- 3 files changed, 81 insertions(+), 125 deletions(-) diff --git a/examples/skeleton/instance_segmentation.py b/examples/skeleton/instance_segmentation.py index cf06d81..acd8d66 100644 --- a/examples/skeleton/instance_segmentation.py +++ b/examples/skeleton/instance_segmentation.py @@ -17,37 +17,35 @@ def parse_args(): parser = argparse.ArgumentParser(description=__doc__.splitlines()[0]) parser.add_argument("--mask_path", required=True, - help="Binary mask to skeletonize; an .h5/.hdf5 or MRC volume.") + help="Binary mask to skeletonize; an .h5 or MRC volume.") parser.add_argument("--seg_key", default=None, - help="HDF5 dataset holding the mask; required for .h5/.hdf5 inputs.") + help="HDF5 dataset holding the mask; required for .h5 inputs.") parser.add_argument("--output_dir", default=".", help="Directory for the skeleton and instance .npz files.") - parser.add_argument("--tag", default="", - help="Suffix appended to the output filenames.") parser.add_argument("--pixel_size", type=float, default=10.0, help="Isotropic voxel spacing for teasar; graph coordinates use these units.") - parser.add_argument("--scale", type=float, default=0.0, + parser.add_argument("--teasar_scale", type=float, default=0.0, help="teasar invalidation-radius scale.") - parser.add_argument("--constant", type=float, default=70.0, + parser.add_argument("--teasar_const", type=float, default=70.0, help="teasar invalidation-radius constant.") parser.add_argument("--number_of_threads", type=int, default=8, help="Threads for teasar.") parser.add_argument("--direction_span", type=int, default=10, help="Nodes walked along each arm to estimate its direction at a junction.") parser.add_argument("--min_through_angle", type=float, default=170.0, - help="Min through-pair angle (deg) for a degree-4 crossing to split.") + help="Min through-pair angle (degrees) for a degree-4 crossing to split.") parser.add_argument("--min_branch_angle", type=float, default=30.0, - help="Min branch angle (deg) for a degree-3 odd arm to be separated.") + help="Min branch angle (degrees) for a degree-3 odd arm to be separated.") parser.add_argument("--tick_length", type=float, default=50.0, - help="Prune dead-end branches shorter than this; 0 disables.") - parser.add_argument("--join_radius", type=float, default=50.0, + help="Prune dead-end branches shorter than this distance; 0 disables.") + parser.add_argument("--join_dist", type=float, default=50.0, help="Join collinear endpoints across gaps up to this distance; 0 disables.") parser.add_argument("--min_join_angle", type=float, default=175.0, - help="Min straightness (deg) for a join; 180 is collinear.") - parser.add_argument("--circle_size", type=float, default=40.0, + help="Min straightness (degrees) for a join; 180 is collinear.") + parser.add_argument("--circle_size", type=float, default=70.0, help="Instance tube diameter for the --view render, in --pixel_size units.") parser.add_argument("--view", action="store_true", - help="Open the pipeline stages in napari.") + help="Open the binary mask and instances in napari.") return parser.parse_args() @@ -69,19 +67,6 @@ def load_mask(mask_path, seg_key): return np.asarray(mrc.data) -def view_stages(binary, stages, pixel_size, circle_size): - import napari - - radius = (circle_size / 2) / pixel_size - viewer = napari.Viewer(ndisplay=3) - viewer.add_labels(binary, name="mask", opacity=0.4) - for index, (name, vertices, edges, _) in enumerate(stages): - labels = connected_components(skeleton_to_graph(vertices, edges)) - volume = draw_instances(vertices / pixel_size, edges, labels, binary.shape, radius) - viewer.add_labels(volume, name=f"instances ({name})", visible=index == len(stages) - 1) - napari.run() - - def main(): args = parse_args() mask = load_mask(args.mask_path, args.seg_key) @@ -89,23 +74,19 @@ def main(): spacing = (args.pixel_size,) * 3 raw_vertices, raw_edges, raw_radii = teasar( - binary, spacing=spacing, scale=args.scale, constant=args.constant, + binary, spacing=spacing, scale=args.teasar_scale, constant=args.teasar_const, number_of_threads=args.number_of_threads, ) - stages = [] vertices, edges, radii = clean_graph( raw_vertices, raw_edges, radii=raw_radii, direction_span=args.direction_span, min_through_angle=args.min_through_angle, min_branch_angle=args.min_branch_angle, tick_length=args.tick_length, - join_radius=args.join_radius, min_join_angle=args.min_join_angle, - save_intermediates=stages, + join_dist=args.join_dist, min_join_angle=args.min_join_angle, ) labels = connected_components(skeleton_to_graph(vertices, edges)) name = Path(args.mask_path).stem - if args.tag: - name = f"{name}_{args.tag}" output_dir = Path(args.output_dir) output_dir.mkdir(parents=True, exist_ok=True) @@ -120,7 +101,14 @@ def main(): print(f"{name}: {len(np.unique(labels))} instances -> {output_dir}") if args.view: - view_stages(binary, stages, args.pixel_size, args.circle_size) + import napari + + radius = (args.circle_size / 2) / args.pixel_size + volume = draw_instances(vertices / args.pixel_size, edges, labels, binary.shape, radius) + viewer = napari.Viewer(ndisplay=3) + viewer.add_labels(binary, name="mask", opacity=0.4) + viewer.add_labels(volume, name="instances") + napari.run() if __name__ == "__main__": diff --git a/src/bioimage_cpp/skeleton/postprocessing.py b/src/bioimage_cpp/skeleton/postprocessing.py index 4fd9a8c..a95d951 100644 --- a/src/bioimage_cpp/skeleton/postprocessing.py +++ b/src/bioimage_cpp/skeleton/postprocessing.py @@ -3,7 +3,7 @@ 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), removed. +the remaining arm(s) are either separated, or for short dead-ends (spurs), pruned. """ from __future__ import annotations @@ -13,6 +13,7 @@ import numpy as np from ..graph import connected_components +from ..utils import UnionFind from ._graph import skeleton_to_graph @@ -33,6 +34,17 @@ def _pair_angle(di, dj): return np.degrees(np.arccos(np.clip(di @ dj, -1.0, 1.0))) +def _compact(vertices, edges, radii): + used = np.zeros(len(vertices), dtype=bool) + if len(edges): + used[edges.reshape(-1)] = True + remap = np.full(len(vertices), -1, dtype=np.int64) + remap[used] = np.arange(int(used.sum())) + if radii is not None: + radii = radii[used] + return vertices[used], remap[edges], radii + + def split_degree3(v, graph, vertices, direction_span=1, min_branch_angle=30.0): adj = np.asarray(graph.node_adjacency(int(v))) if adj.shape[0] != 3: @@ -78,7 +90,7 @@ def clean_graph( min_through_angle: float = 160.0, min_branch_angle: float = 30.0, tick_length: float = 0.0, - join_radius: float = 0.0, + join_dist: float = 0.0, min_join_angle: float = 175.0, save_intermediates: list | None = None, ) -> tuple[np.ndarray, np.ndarray, np.ndarray | None]: @@ -92,8 +104,8 @@ def clean_graph( from the through pair by at least ``min_branch_angle``. 3. Split each degree-4 crossing, separating its two through pairs when they are collinear to within ``min_through_angle``. - 4. If ``join_radius > 0``, reconnect collinear endpoints across gaps less - than this distance via :func:`join_close_components`. + 4. If ``join_dist > 0``, join collinear endpoints across gaps up to + this distance via :func:`join_close_components`. Parameters ---------- @@ -112,8 +124,8 @@ def clean_graph( through pair for the odd arm to be separated. tick_length: If > 0, prune dead-end branches shorter than this (physical) distance. - join_radius: - If > 0, reconnect collinear endpoints across gaps up to this (physical) + join_dist: + If > 0, join collinear endpoints across gaps up to this (physical) distance. min_join_angle: Minimum straightness (degrees) required for a join; 180 is collinear. @@ -175,20 +187,12 @@ def _snapshot(name): keep[list(prune_edges)] = False edges = edges[keep] - used = np.zeros(len(vertices), dtype=bool) - if len(edges): - used[edges.reshape(-1)] = True - remap = np.full(len(vertices), -1, dtype=np.int64) - remap[used] = np.arange(int(used.sum())) - vertices = vertices[used] - if radii is not None: - radii = radii[used] - edges = remap[edges] + vertices, edges, radii = _compact(vertices, edges, radii) _snapshot("split") - if join_radius > 0: + if join_dist > 0: vertices, edges, radii = join_close_components( - vertices, edges, join_radius, + vertices, edges, join_dist, min_join_angle=min_join_angle, direction_span=direction_span, radii=radii, ) _snapshot("join") @@ -208,15 +212,15 @@ def _adjacency(num_nodes, edges): def remove_ticks(vertices, edges, tick_length, radii=None): - """Remove short dead-end branches ("ticks") from a skeleton graph. + """Prune short dead-end branches ("ticks") from a skeleton graph. Ports kimimaro's `remove_ticks`. A distance graph is built over the critical points (terminals, degree 1; branch points, degree >= 3), whose superedges are the paths between them weighted by physical length. The shortest terminal - branch below ``tick_length`` is removed repeatedly; when a branch point drops + branch below ``tick_length`` is pruned repeatedly; when a branch point drops to degree 2 its two superedges are fused into one, so a real filament end is re-measured rather than clipped. Standalone paths (both ends terminal) are - never removed. + never pruned. Parameters ---------- @@ -225,7 +229,7 @@ def remove_ticks(vertices, edges, tick_length, radii=None): edges: Integer array with shape ``(E, 2)`` indexing ``vertices``. tick_length: - Maximum branch length (physical) that may be culled. + Maximum branch length (physical) that may be pruned. radii: Optional per-vertex radii, carried through the same remapping. @@ -305,15 +309,7 @@ def remove_ticks(vertices, edges, tick_length, radii=None): else: edges = edges.copy() - used = np.zeros(num_nodes, dtype=bool) - if len(edges): - used[edges.reshape(-1)] = True - remap = np.full(num_nodes, -1, dtype=np.int64) - remap[used] = np.arange(int(used.sum())) - vertices = vertices[used] - if radii is not None: - radii = radii[used] - edges = remap[edges] + vertices, edges, radii = _compact(vertices, edges, radii) return vertices, edges, radii @@ -333,23 +329,23 @@ def _endpoint_tangent(endpoint, indptr, dst, degrees, vertices, span): return direction / norm if norm > 0 else None -def join_close_components(vertices, edges, radius, *, min_join_angle=175.0, +def join_close_components(vertices, edges, dist, *, min_join_angle=175.0, direction_span=5, radii=None): """Reconnect fragmented filaments by joining collinear endpoints across gaps. Endpoints (degree = 1) of different connected components are joined with a - new edge when they are within ``radius`` and the two fragments are nearly + new edge when they are within ``dist`` and the two fragments are nearly collinear through the gap: the outward tangent at each endpoint must point along the gap to within ``180 - min_join_angle`` degrees, so a straight continuation reads ~180 (``min_join_angle``). Joins are made shortest-first - with a union-find so each pair of components is linked at most once and each + with a union-find so each pair of components is joined at most once and each endpoint is used once. Parameters ---------- vertices, edges: Skeleton graph; ``edges`` indexes ``vertices``. - radius: + dist: Maximum gap (physical) across which endpoints may be joined. min_join_angle: Minimum straightness (degrees) of the joined path; 180 is perfectly @@ -381,7 +377,7 @@ def join_close_components(vertices, edges, radius, *, min_join_angle=175.0, tol = np.deg2rad(180.0 - min_join_angle) tree = cKDTree(vertices[endpoints]) candidates = [] - for ia, ib in tree.query_pairs(radius): + for ia, ib in tree.query_pairs(dist): a, b = int(endpoints[ia]), int(endpoints[ib]) if labels[a] == labels[b]: continue @@ -400,25 +396,15 @@ def join_close_components(vertices, edges, radius, *, min_join_angle=175.0, candidates.append((dist, a, b)) candidates.sort() - parent = {} - - def find(x): - parent.setdefault(x, x) - root = x - while parent[root] != root: - root = parent[root] - while parent[x] != root: - parent[x], x = root, parent[x] - return root + uf = UnionFind(int(labels.max()) + 1) used, new_edges = set(), [] for _, a, b in candidates: if a in used or b in used: continue - ra, rb = find(int(labels[a])), find(int(labels[b])) - if ra == rb: + if uf.find(int(labels[a])) == uf.find(int(labels[b])): continue - parent[ra] = rb + uf.merge(int(labels[a]), int(labels[b])) new_edges.append([a, b]) used.add(a) used.add(b) @@ -430,13 +416,6 @@ def find(x): return vertices, edges, radii -def _line_voxels(start, stop): - delta = stop - start - steps = int(np.abs(delta).max()) + 1 - t = np.linspace(0.0, 1.0, steps) - return np.rint(start[None, :] + t[:, None] * delta[None, :]).astype(np.int64) - - def draw_instances(vertices, edges, labels, shape, radius=1): """Rasterize a labeled skeleton graph into a dense instance volume. @@ -480,7 +459,10 @@ def draw_instances(vertices, edges, labels, shape, radius=1): vi = np.rint(vertices).astype(np.int64) centers, center_labels = [], [] for a, b in edges: - pts = _line_voxels(vi[a], vi[b]) + delta = vi[b] - vi[a] + steps = int(np.abs(delta).max()) + 1 + t = np.linspace(0.0, 1.0, steps) + pts = np.rint(vi[a][None, :] + t[:, None] * delta[None, :]).astype(np.int64) centers.append(pts) center_labels.append(np.full(len(pts), int(labels[a]) + 1)) centers = np.concatenate(centers) diff --git a/tests/skeleton/test_postprocessing.py b/tests/skeleton/test_postprocessing.py index 42b96d2..65f408d 100644 --- a/tests/skeleton/test_postprocessing.py +++ b/tests/skeleton/test_postprocessing.py @@ -16,11 +16,7 @@ def _component_count(vertices, edges): return len(np.unique(labels)) -def _edge_rows(edges): - return {tuple(sorted(map(int, edge))) for edge in np.asarray(edges)} - - -def test_split_degree3_prunes_the_arm_that_diverges_from_the_through_pair(): +def test_split_degree3_splits_odd_arm(): # Node 0 is a T: a straight through-pair along x (nodes 1, 2) plus an odd arm along y (node 3). vertices, _, graph = _graph( [[0.0, 0.0, 0.0], [0.0, 0.0, 1.0], [0.0, 0.0, -1.0], [0.0, 1.0, 0.0]], @@ -29,7 +25,7 @@ def test_split_degree3_prunes_the_arm_that_diverges_from_the_through_pair(): assert split_degree3(0, graph, vertices, direction_span=1, min_branch_angle=30.0) == [2] -def test_split_degree3_keeps_a_nearly_collinear_arm(): +def test_split_degree3_keeps_collinear_arm(): # The third arm now runs almost parallel to the through pair, below min_branch_angle. vertices, _, graph = _graph( [[0.0, 0.0, 0.0], [0.0, 0.0, 1.0], [0.0, 0.0, -1.0], [0.0, 1.0, 10.0]], @@ -38,7 +34,7 @@ def test_split_degree3_keeps_a_nearly_collinear_arm(): assert split_degree3(0, graph, vertices, direction_span=1, min_branch_angle=30.0) is None -def test_split_degree4_detaches_one_collinear_through_pair(): +def test_split_degree4_splits_through_pair(): # A '+' crossing: edges 0,1 are the x through-pair, edges 2,3 the y through-pair. vertices, _, graph = _graph( [[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]], @@ -48,7 +44,7 @@ def test_split_degree4_detaches_one_collinear_through_pair(): assert set(result) in ({0, 1}, {2, 3}) -def test_split_degree4_keeps_a_crossing_without_a_straight_through_pair(): +def test_split_degree4_keeps_non_collinear_crossing(): # Four arms at 0, 60, 120, 180 degrees: no pair is collinear enough to split. angle = np.deg2rad([0.0, 60.0, 120.0, 180.0]) arms = np.stack([np.zeros_like(angle), np.sin(angle), np.cos(angle)], axis=1) @@ -63,7 +59,7 @@ def test_split_degree4_keeps_a_crossing_without_a_straight_through_pair(): "tick_length, n_vertices, n_edges", [(1.5, 5, 4), (0.5, 6, 5)], ) -def test_remove_ticks_prunes_only_spurs_below_the_threshold(tick_length, n_vertices, n_edges): +def test_remove_ticks_prunes_spurs_below_threshold(tick_length, n_vertices, n_edges): # Backbone 0-1-2-3-4 along x with a length-1 spur (node 5) off the branch node 2. vertices = np.array( [[0.0, 0.0, 0.0], [0.0, 0.0, 1.0], [0.0, 0.0, 2.0], [0.0, 0.0, 3.0], @@ -77,7 +73,7 @@ def test_remove_ticks_prunes_only_spurs_below_the_threshold(tick_length, n_verti assert len(out_edges) == n_edges -def test_remove_ticks_never_removes_a_standalone_path(): +def test_remove_ticks_keeps_standalone_path(): vertices = np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 1.0], [0.0, 0.0, 2.0]]) edges = np.array([[0, 1], [1, 2]], dtype=np.int64) @@ -87,7 +83,7 @@ def test_remove_ticks_never_removes_a_standalone_path(): assert len(out_edges) == 2 -def test_remove_ticks_on_empty_edges_returns_the_inputs(): +def test_remove_ticks_empty_edges(): vertices = np.zeros((3, 3), dtype=np.float64) edges = np.empty((0, 2), dtype=np.int64) @@ -97,46 +93,46 @@ def test_remove_ticks_on_empty_edges_returns_the_inputs(): assert len(out_edges) == 0 -def test_join_close_components_links_collinear_endpoints_within_radius(): +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]]) edges = np.array([[0, 1], [2, 3]], dtype=np.int64) - _, out_edges, _ = bic.skeleton.join_close_components(vertices, edges, radius=2.5) + _, out_edges, _ = bic.skeleton.join_close_components(vertices, edges, dist=2.5) assert len(out_edges) == 3 - assert (1, 2) in _edge_rows(out_edges) + assert (1, 2) in {tuple(sorted(map(int, edge))) for edge in out_edges} -def test_join_close_components_ignores_endpoints_beyond_the_radius(): +def test_join_close_components_ignores_beyond_dist(): 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]]) edges = np.array([[0, 1], [2, 3]], dtype=np.int64) - _, out_edges, _ = bic.skeleton.join_close_components(vertices, edges, radius=1.5) + _, out_edges, _ = bic.skeleton.join_close_components(vertices, edges, dist=1.5) assert len(out_edges) == 2 -def test_join_close_components_ignores_non_collinear_endpoints(): +def test_join_close_components_ignores_non_collinear(): # Two fragments offset in y so the gap bends away from each endpoint tangent. vertices = np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 1.0], [0.0, 3.0, 2.0], [0.0, 4.0, 2.0]]) edges = np.array([[0, 1], [2, 3]], dtype=np.int64) - _, out_edges, _ = bic.skeleton.join_close_components(vertices, edges, radius=3.5) + _, out_edges, _ = bic.skeleton.join_close_components(vertices, edges, dist=3.5) assert len(out_edges) == 2 -def test_join_close_components_does_not_join_within_one_component(): +def test_join_close_components_skips_same_component(): vertices = np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 1.0], [0.0, 0.0, 2.0]]) edges = np.array([[0, 1], [1, 2]], dtype=np.int64) - _, out_edges, _ = bic.skeleton.join_close_components(vertices, edges, radius=5.0) + _, out_edges, _ = bic.skeleton.join_close_components(vertices, edges, dist=5.0) assert len(out_edges) == 2 -def test_clean_graph_splits_a_crossing_into_two_components(): +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]], ) @@ -146,24 +142,14 @@ def test_clean_graph_splits_a_crossing_into_two_components(): assert _component_count(vertices, edges) == 1 out_vertices, out_edges, out_radii = bic.skeleton.clean_graph( - vertices, edges, radii=radii, direction_span=1, tick_length=0.0, join_radius=0.0, + vertices, edges, radii=radii, direction_span=1, tick_length=0.0, join_dist=0.0, ) assert _component_count(out_vertices, out_edges) == 2 assert len(out_radii) == len(out_vertices) -def test_clean_graph_records_every_stage_snapshot(): - vertices = np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 1.0], [0.0, 0.0, 2.0]]) - edges = np.array([[0, 1], [1, 2]], dtype=np.int64) - stages = [] - - bic.skeleton.clean_graph(vertices, edges, save_intermediates=stages) - - assert [name for name, *_ in stages] == ["raw", "ticks", "split", "join"] - - -def test_draw_instances_rasterizes_the_edge_with_the_component_label(): +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) labels = np.array([2, 2]) @@ -174,7 +160,7 @@ def test_draw_instances_rasterizes_the_edge_with_the_component_label(): assert set(np.unique(volume)) == {0, 3} -def test_draw_instances_dilates_by_the_radius(): +def test_draw_instances_dilates_by_radius(): vertices = np.array([[2.0, 2.0, 0.0], [2.0, 2.0, 4.0]]) edges = np.array([[0, 1]], dtype=np.int64) labels = np.array([0, 0]) @@ -182,10 +168,10 @@ def test_draw_instances_dilates_by_the_radius(): thin = bic.skeleton.draw_instances(vertices, edges, labels, (5, 5, 5), radius=0) thick = bic.skeleton.draw_instances(vertices, edges, labels, (5, 5, 5), radius=1) - assert int(thick.sum() > 0) and thick.astype(bool).sum() > thin.astype(bool).sum() + assert thick.astype(bool).sum() > thin.astype(bool).sum() -def test_draw_instances_on_empty_edges_returns_a_zero_volume(): +def test_draw_instances_empty_edges(): vertices = np.zeros((2, 3), dtype=np.float64) edges = np.empty((0, 2), dtype=np.int64) labels = np.zeros(2, dtype=np.int64) From 54f32cdb236ef5fd5b71e79a81d576948c047360 Mon Sep 17 00:00:00 2001 From: GOESTERN-1107078 Date: Tue, 21 Jul 2026 18:49:33 +0200 Subject: [PATCH 3/3] fix init file merge conflict --- src/bioimage_cpp/skeleton/__init__.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/bioimage_cpp/skeleton/__init__.py b/src/bioimage_cpp/skeleton/__init__.py index a2632fe..f8de89f 100644 --- a/src/bioimage_cpp/skeleton/__init__.py +++ b/src/bioimage_cpp/skeleton/__init__.py @@ -203,8 +203,7 @@ def teasar_labels( return run(labels_c, labels_array.dtype.type(background_value), *options) -__all__ = ["clean_graph", "draw_instances", "join_close_components", "remove_ticks", "skeleton_to_graph", "teasar", "teasar_labels"] from . import distributed -__all__ = ["distributed", "skeleton_to_graph", "teasar", "teasar_labels"] +__all__ = ["clean_graph", "distributed", "draw_instances", "join_close_components", "remove_ticks", "skeleton_to_graph", "teasar", "teasar_labels"]