From 37b55e02de7fb08495b4977dee0b4d84ec40fedc Mon Sep 17 00:00:00 2001 From: uermel Date: Tue, 6 Oct 2026 04:15:10 -0700 Subject: [PATCH] fix: split traced filaments at sharp bends without a junction Two filaments that touch end to end skeletonize as one bent branch with no junction, so seg2fil traced them as one hairpin filament. A chain from a semantic segmentation is now split wherever it turns by more than --max-bend over the direction window, the rule junctions use. Closed chains and instances are kept whole, and the run summary counts the splits. Each piece is stored with the default Catmull-Rom curve like any other filament. On 10426 actin, kinks over 45 degrees drop from 4.9% to 1.5% of filaments. None of the 102 traced 10521 microtubules changes. Co-Authored-By: Claude Opus 5.5 --- src/copick_utils/cli/seg2fil.py | 7 +- src/copick_utils/cli/util.py | 2 +- .../converters/filaments_from_segmentation.py | 4 +- src/copick_utils/process/filament_tracing.py | 67 +++++++++++++++++-- tests/test_filaments.py | 14 ++++ 5 files changed, 85 insertions(+), 9 deletions(-) diff --git a/src/copick_utils/cli/seg2fil.py b/src/copick_utils/cli/seg2fil.py index c6ca9f9..a495164 100644 --- a/src/copick_utils/cli/seg2fil.py +++ b/src/copick_utils/cli/seg2fil.py @@ -103,7 +103,8 @@ def seg2fil( Holes up to `--fill-lumen` are filled first, so a tube labelled by its wall alone traces as one filament. Side branches shorter than `--prune-length` are pruned, and at each junction the branches that continue most nearly straight (within `--max-bend`) are joined, so two crossing - filaments stay two filaments. Free ends are extended to the edge of the segmentation, and + filaments stay two filaments; a path that bends more sharply than `--max-bend` (two filaments + touching end to end) is split there. Free ends are extended to the edge of the segmentation, and filaments shorter than `--min-length`, shorter than `--min-aspect` label diameters (blobs), or thinner than `--min-radius` (slivers of segmentation noise) are rejected. Options left unset are derived from the segmentation's own thickness or the object's radius. @@ -217,6 +218,7 @@ def seg2fil( short = sum(result.get("rejected_length", 0) for result in results.values() if result) blobs = sum(result.get("rejected_aspect", 0) for result in results.values() if result) thin = sum(result.get("rejected_radius", 0) for result in results.values() if result) + bends = sum(result.get("split_at_bends", 0) for result in results.values() if result) all_errors = [] for result in results.values(): @@ -225,7 +227,8 @@ def seg2fil( logger.info(f"Completed: {successful}/{len(results)} runs processed successfully") logger.info( - f"Total filaments: {filaments} ({length / 10:.0f} nm); rejected {short} short, {blobs} blob-like, {thin} thin", + f"Total filaments: {filaments} ({length / 10:.0f} nm); rejected {short} short, {blobs} blob-like, {thin} thin; " + f"split {bends} at sharp bends", ) if all_errors: diff --git a/src/copick_utils/cli/util.py b/src/copick_utils/cli/util.py index 05e78eb..68f784e 100644 --- a/src/copick_utils/cli/util.py +++ b/src/copick_utils/cli/util.py @@ -812,7 +812,7 @@ def add_filament_tracing_options_decorator(_func: click.Command) -> click.Comman type=float, default=45.0, help="Largest deviation from a straight line, in degrees, for a filament to continue through a " - "junction.", + "junction or around a bend; a sharper bend (two filaments touching end to end) splits it.", ), optgroup.option( "--smoothing", diff --git a/src/copick_utils/converters/filaments_from_segmentation.py b/src/copick_utils/converters/filaments_from_segmentation.py index 240cb45..6cf5d96 100644 --- a/src/copick_utils/converters/filaments_from_segmentation.py +++ b/src/copick_utils/converters/filaments_from_segmentation.py @@ -94,7 +94,7 @@ def filaments_from_segmentation( (the filament's tube radius), or no radius filter if the object has none. prune_length: Prune skeleton spurs shorter than this (Å). None: one label diameter. junction_merge: Merge junctions joined by a bridge up to this long (Å). None: one label diameter. - max_bend: Largest deviation from a straight line (degrees) to continue through a junction. + max_bend: Largest deviation from a straight line (degrees) to continue through a junction or around a bend. smoothing: RMS deviation of the fitted spline from the skeleton (Å). None: half a voxel. extend_ends: Extend free ends to the edge of the label. fill_lumen: Fill holes up to this radius (Å) in each axis-aligned slice before tracing, the empty lumen of @@ -199,7 +199,7 @@ def filaments_from_segmentation( logger.info( f"Traced {len(centrelines)} filaments ({report.total_length / 10:.0f} nm) in {run.name}; " f"rejected {report.rejected_length} short, {report.rejected_aspect} blob-like and " - f"{report.rejected_radius} thin", + f"{report.rejected_radius} thin; split {report.split_at_bends} at sharp bends", ) return filaments, stats diff --git a/src/copick_utils/process/filament_tracing.py b/src/copick_utils/process/filament_tracing.py index 95dbb1a..9cb3c9d 100644 --- a/src/copick_utils/process/filament_tracing.py +++ b/src/copick_utils/process/filament_tracing.py @@ -21,9 +21,11 @@ 3. At every junction, branches are paired by straightest continuation: the pair whose directions deviate least from a straight line is joined first, as long as the deviation is at most ``max_bend`` degrees. A crossing therefore gives two filaments, and a Y gives one filament through the junction and one that ends there. -4. Paired branches are chained into filaments. At a free end, the last label radius of skeleton is replaced by the - continuation of the filament's direction up to the edge of the label (thinning shortens a rod's skeleton by about - one radius at each end and often bends its last voxels). +4. Paired branches are chained into filaments. A chain that turns by more than ``max_bend`` degrees between its + directions over ``direction_window`` before and after a point is split there: two filaments that touch end to end + (a V) skeletonize as one bent branch with no junction, which the pairing never sees. At a free end, the last label + radius of skeleton is replaced by the continuation of the filament's direction up to the edge of the label + (thinning shortens a rod's skeleton by about one radius at each end and often bends its last voxels). 5. Filaments shorter than ``min_length``, shorter than ``min_aspect`` label diameters (blobs), or thinner than ``min_radius`` (slivers of segmentation noise, much thinner than the real filaments) are rejected. 6. Each filament is fitted with a smoothing B-spline (scipy ``splprep``) whose RMS deviation from the skeleton is @@ -73,7 +75,12 @@ class TraceParameters: tube radius (40 Å there). prune_length: Prune skeleton spurs shorter than this. None: one label diameter. junction_merge: Merge junctions joined by a bridge up to this long. None: one label diameter. - max_bend: Largest deviation from a straight line, in degrees, for two branches to continue through a junction. + max_bend: Largest deviation from a straight line, in degrees, for two branches to continue through a junction, + and for a chain to continue around a bend over ``direction_window`` (a sharper one is split). Basis: + measured that way, a filament turns by more than 45 degrees only on a radius below about 2.5 label + diameters. On dataset 10521 splitting changes none of the 102 microtubules traced on 58 tomograms; on + dataset 10426 (easymode actin) it splits paths where two filaments touch, turning by 45-159 degrees + (median 55, 51 splits on three tomograms). smoothing: RMS deviation of the fitted spline from the skeleton. None: half a voxel. extend_ends: Extend free ends along their direction to the edge of the label. direction_window: Length over which a branch's direction at a junction is measured. None: two label @@ -139,6 +146,7 @@ class TraceReport: rejected_aspect: int = 0 rejected_radius: int = 0 junctions_resolved: int = 0 + split_at_bends: int = 0 total_length: float = 0.0 dropped_branch_length: float = 0.0 @@ -351,6 +359,48 @@ def follow(start_branch: int, enter_at_start: bool) -> Chain: return chains, len(link) // 2 +def split_at_bends(chain: Chain, max_bend: float, window: float, trim: float = 0.0) -> List[Chain]: + """Split a chain wherever it turns by more than ``max_bend`` degrees, the turn at a point being the angle between + the chain's directions over ``window`` voxels of path before and after it. + + Two filaments that touch end to end (a V, or one ending against another at an acute angle) skeletonize as one + bent branch with no junction, so pairing at junctions never sees the turn. The new ends are not free: they lie + against the other filament, so they are not extended. Like a branch at a junction, each piece loses ``trim`` + voxels of path at the split, where thinning bends it towards the other. Closed chains, and chains shorter than + two windows, are kept whole. + + Returns: + The pieces, in order along the chain (the chain itself if it has no such bend). + """ + coords = chain.coords.astype(float) + if chain.closed or len(coords) < 3 or window <= 0: + return [chain] + reach = np.concatenate([[0.0], np.cumsum(np.linalg.norm(np.diff(coords, axis=0), axis=1))]) + inner = np.nonzero((reach >= window) & (reach <= reach[-1] - window))[0] + if inner.size == 0: + return [chain] + before = coords[inner] - coords[np.searchsorted(reach, reach[inner] - window)] + after = coords[np.minimum(np.searchsorted(reach, reach[inner] + window), len(coords) - 1)] - coords[inner] + norms = np.linalg.norm(before, axis=1) * np.linalg.norm(after, axis=1) + cosine = np.where(norms > 1e-9, np.sum(before * after, axis=1) / np.maximum(norms, 1e-9), 1.0) + bend = np.degrees(np.arccos(np.clip(cosine, -1.0, 1.0))) + k = int(np.argmax(bend)) + if bend[k] <= max_bend: + return [chain] + at = int(inner[k]) + head = Chain( + coords=_trim(chain.coords[: at + 1], trim, at_start=False) if trim > 0 else chain.coords[: at + 1], + free_ends=(chain.free_ends[0], False), + junctions=chain.junctions, + ) + tail = Chain( + coords=_trim(chain.coords[at:], trim, at_start=True) if trim > 0 else chain.coords[at:], + free_ends=(False, chain.free_ends[1]), + junctions=chain.junctions, + ) + return split_at_bends(head, max_bend, window, trim) + split_at_bends(tail, max_bend, window, trim) + + # --------------------------------------------------------------------------------------------------------------------- # Fitting # --------------------------------------------------------------------------------------------------------------------- @@ -527,6 +577,15 @@ def voxels(value: Optional[float], default: float) -> float: junction_trim=radius, ) report.junctions_resolved += resolved + if not single: + # An instance is one filament by definition; in a semantic label, a sharp bend is two filaments touching. + window = voxels(params.direction_window, 2.0 * diameter) + pieces = [] + for chain in chains: + parts = split_at_bends(chain, params.max_bend, window, trim=radius) + report.split_at_bends += len(parts) - 1 + pieces.extend(parts) + chains = pieces if single and len(chains) > 1: lengths = [np.sum(np.linalg.norm(np.diff(c.coords, axis=0), axis=1)) for c in chains] keep = int(np.argmax(lengths)) diff --git a/tests/test_filaments.py b/tests/test_filaments.py index 75de5ac..780f9cf 100644 --- a/tests/test_filaments.py +++ b/tests/test_filaments.py @@ -132,6 +132,20 @@ def test_y_junction_continues_the_straightest_pair(): assert report.junctions_resolved == 1 +def test_two_filaments_touching_end_to_end_are_two_filaments(): + # A V: two filaments that meet at their ends (actin on dataset 10426) skeletonize as one bent branch with no + # junction, so pairing at junctions never sees the turn. + volume = tube(np.zeros(SHAPE, np.uint8), (20, 30, 10), (20, 45, 60), 3) + tube(volume, (20, 45, 60), (20, 60, 10), 3) # 33 degrees between the arms + centrelines, _, report = trace_centrelines(volume, VS, TraceParameters(min_length=150)) + assert len(centrelines) == 2 and report.split_at_bends == 1 + for c in centrelines: + points, tangents, _ = evaluate_spline(c.tck, step=VS / 2) + overall = (points[-1] - points[0]) / np.linalg.norm(points[-1] - points[0]) + assert np.degrees(np.arccos(np.clip(tangents @ overall, -1, 1))).max() < 30.0 + assert c.length > 400 + + def test_ring_is_one_closed_filament(): volume = np.zeros(SHAPE, np.uint8) zz, yy, xx = np.indices(SHAPE)