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)