fix: split traced filaments at sharp bends without a junction

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:

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",

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
 

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))

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)