skeletor.post

 28def radii(s, mesh=None, method='knn', aggregate='mean', validate=False, **kwargs):
 29    """Extract radii for given skeleton table.
 30
 31    Important
 32    ---------
 33    This function really only produces useful radii if the skeleton is centered
 34    inside the mesh. `by_wavefront` does that by default whereas all other
 35    skeletonization methods don't. Your best bet to get centered skeletons is
 36    to contract the mesh first (`sk.pre.contract`).
 37
 38    Parameters
 39    ----------
 40    s :         skeletor.Skeleton
 41                Skeleton to clean up.
 42    mesh :      trimesh.Trimesh, optional
 43                Original mesh (e.g. before contraction). If not provided will
 44                use the mesh associated with ``s``.
 45    method :    "knn" | "ray"
 46                Whether and how to add radius information to each node::
 47
 48                    - "knn" uses k-nearest-neighbors to get radii: fast but
 49                      potential for being very wrong
 50                    - "ray" uses ray-casting to get radii: slower but sometimes
 51                      less wrong
 52                      (requires `ncollpyde`)
 53
 54    aggregate : "mean" | "median" | "max" | "min" | "percentile75"
 55                Function used to aggregate radii over sample (i.e. across
 56                k nearest-neighbors or ray intersections)
 57    validate :  bool
 58                If True, will try to fix potential issues with the mesh
 59                (e.g. infinite values, duplicate vertices, degenerate faces)
 60                before skeletonization. Note that this might make changes to
 61                your mesh inplace!
 62    **kwargs
 63                Keyword arguments are passed to the respective method:
 64
 65                For method "knn"::
 66
 67                    n :             int (default 5)
 68                                    Radius will be the mean over n nearest-neighbors.
 69
 70                For method "ray"::
 71
 72                    n_rays :        int (default 20)
 73                                    Number of rays to cast for each node.
 74                    projection :    "sphere" (default) | "tangents"
 75                                    Whether to cast rays in a sphere around each node or in a
 76                                    circle orthogonally to the node's tangent vector.
 77                    fallback :      "knn" (default) | None | number
 78                                    If a point is outside or right on the surface of the mesh
 79                                    the raycasting will return nonesense results. We can either
 80                                    ignore those cases (``None``), assign a arbitrary number or
 81                                    we can fall back to radii from k-nearest-neighbors (``knn``).
 82
 83    Returns
 84    -------
 85    None
 86                    But attaches `radius` to the skeleton's SWC table. Existing
 87                    values are replaced!
 88
 89    """
 90    if isinstance(mesh, type(None)):
 91        mesh = s.mesh
 92
 93    mesh = make_trimesh(mesh, validate=True)
 94
 95    if method == 'knn':
 96        radius = get_radius_knn(s.swc[['x', 'y', 'z']].values,
 97                                aggregate=aggregate,
 98                                mesh=mesh, **kwargs)
 99    elif method == 'ray':
100        radius = get_radius_ray(s.swc,
101                                mesh=mesh,
102                                aggregate=aggregate,
103                                **kwargs)
104    else:
105        raise ValueError(f'Unknown method "{method}"')
106
107    s.swc['radius'] = radius
108
109    return

29def clean_up(s, mesh=None, validate=False, inplace=False, **kwargs):
30    """Clean up the skeleton.
31
32    Requires `ncollpyde` (`pip install ncollpyde`).
33
34    This function bundles a bunch of procedures to clean up the skeleton:
35
36      1. Remove twigs that are running parallel to their parent branch
37      2. Move nodes outside the mesh back inside (or at least snap to surface)
38
39    Note that this is not a magic bullet and some of this will not work (well)
40    if the original mesh was degenerate (e.g. internal faces or not watertight)
41    to begin with.
42
43    Parameters
44    ----------
45    s :         skeletor.Skeleton
46                Skeleton to clean up.
47    mesh :      trimesh.Trimesh, optional
48                Original mesh (e.g. before contraction). If not provided will
49                use the mesh associated with ``s``.
50    validate :  bool
51                If True, will try to fix potential issues with the mesh
52                (e.g. infinite values, duplicate vertices, degenerate faces)
53                before cleaning up. Note that this might change your mesh
54                inplace!
55    inplace :   bool
56                If False will make and return a copy of the skeleton. If True,
57                will modify the `s` inplace.
58
59    **kwargs
60                Keyword arguments are passed to the bundled function.
61
62                For `skeletor.postprocessing.drop_parallel_twigs`::
63
64                 theta :     float (default 0.01)
65                             For each twig we generate the dotproduct between the tangent
66                             vectors of it and its parents. If these line up perfectly the
67                             dotproduct will equal 1. ``theta`` determines how much that
68                             value can differ from 1 for us to still prune the twig: higher
69                             theta = more pruning.
70
71    Returns
72    -------
73    s_clean :   skeletor.Skeleton
74                Hopefully improved skeleton.
75
76    """
77    if isinstance(mesh, type(None)):
78        mesh = s.mesh
79
80    mesh = make_trimesh(mesh, validate=validate)
81
82    if not inplace:
83        s = s.copy()
84
85    # Drop parallel twigs
86    _ = drop_parallel_twigs(s, theta=kwargs.get("theta", 0.01), inplace=True)
87
88    # Recenter vertices
89    _ = recenter_vertices(s, mesh, inplace=True)
90
91    return s

846def smooth(
847    s, window: int = 3, to_smooth: list = ["x", "y", "z"], inplace: bool = False
848):
849    """Smooth skeleton using rolling windows.
850
851    Parameters
852    ----------
853    s :             skeletor.Skeleton
854                    Skeleton to be processed.
855    window :        int, optional
856                    Size (N observations) of the rolling window in number of
857                    nodes.
858    to_smooth :     list
859                    Columns of the node table to smooth. Should work with any
860                    numeric column (e.g. 'radius').
861    inplace :       bool
862                    If False will make and return a copy of the skeleton. If
863                    True, will modify the `s` inplace.
864
865    Returns
866    -------
867    s :             skeletor.Skeleton
868                    Skeleton with smoothed node table.
869
870    """
871    if not inplace:
872        s = s.copy()
873
874    # Prepare nodes (add parent_dist for later, set index)
875    nodes = s.swc.set_index("node_id", inplace=False).copy()
876
877    to_smooth = np.array(to_smooth)
878    miss = to_smooth[~np.isin(to_smooth, nodes.columns)]
879    if len(miss):
880        raise ValueError(f"Column(s) not found in node table: {miss}")
881
882    # Go over each segment and smooth
883    for seg in s.get_segments():
884        # Get this segment's parent distances and get cumsum
885        this_co = nodes.loc[seg, to_smooth]
886
887        interp = this_co.rolling(window, min_periods=1).mean()
888
889        nodes.loc[seg, to_smooth] = interp.values
890
891    # Reassign nodes
892    s.swc = nodes.reset_index(drop=False, inplace=False)
893
894    return s

897def despike(s, sigma=5, max_spike_length=1, inplace=False, reverse=False):
898    r"""Remove spikes in skeleton.
899
900    For each node A, the euclidean distance to its next successor (parent)
901    B and that node's successor C (i.e A->B->C) is computed. If
902    :math:`\\frac{dist(A,B)}{dist(A,C)}>sigma`, node B is considered a spike
903    and realigned between A and C.
904
905    Parameters
906    ----------
907    x :                 skeletor.Skeleton
908                        Skeleton to be processed.
909    sigma :             float | int, optional
910                        Threshold for spike detection. Smaller sigma = more
911                        aggressive spike detection.
912    max_spike_length :  int, optional
913                        Determines how long (# of nodes) a spike can be.
914    inplace :           bool, optional
915                        If False, a copy of the neuron is returned.
916    reverse :           bool, optional
917                        If True, will **also** walk the segments from proximal
918                        to distal. Use this to catch spikes on e.g. terminal
919                        nodes.
920
921    Returns
922    -------
923    s                   skeletor.Skeleton
924                        Despiked neuron.
925
926    """
927    if not inplace:
928        s = s.copy()
929
930    # Index nodes table by node ID
931    this_nodes = s.swc.set_index("node_id", inplace=False)
932
933    segments = s.get_segments()
934    segs_to_walk = segments
935
936    if reverse:
937        segs_to_walk += segs_to_walk[::-1]
938
939    # For each spike length do -> do this in reverse to correct the long
940    # spikes first
941    for l in list(range(1, max_spike_length + 1))[::-1]:
942        # Go over all segments
943        for seg in segs_to_walk:
944            # Get nodes A, B and C of this segment
945            this_A = this_nodes.loc[seg[: -l - 1]]
946            this_B = this_nodes.loc[seg[l:-1]]
947            this_C = this_nodes.loc[seg[l + 1 :]]
948
949            # Get coordinates
950            A = this_A[["x", "y", "z"]].values
951            B = this_B[["x", "y", "z"]].values
952            C = this_C[["x", "y", "z"]].values
953
954            # Calculate euclidian distances A->B and A->C
955            dist_AB = np.linalg.norm(A - B, axis=1)
956            dist_AC = np.linalg.norm(A - C, axis=1)
957
958            # Get the spikes
959            spikes_ix = np.where(
960                np.divide(dist_AB, dist_AC, where=dist_AC != 0) > sigma
961            )[0]
962            spikes = this_B.iloc[spikes_ix]
963
964            if not spikes.empty:
965                # Interpolate new position(s) between A and C
966                new_positions = A[spikes_ix] + (C[spikes_ix] - A[spikes_ix]) / 2
967
968                this_nodes.loc[spikes.index, ["x", "y", "z"]] = new_positions
969
970    # Reassign node table
971    s.swc = this_nodes.reset_index(drop=False, inplace=False)
972
973    return s

 94def remove_bristles(s, mesh=None, los_only=False, inplace=False):
 95    """Remove "bristles" that sometimes occurr along the backbone.
 96
 97    Works by finding terminal twigs that consist of only a single node.
 98
 99    Parameters
100    ----------
101    s :         skeletor.Skeleton
102                Skeleton to clean up.
103    mesh :      trimesh.Trimesh, optional
104                Original mesh (e.g. before contraction). If not provided will
105                use the mesh associated with ``s``.
106    los_only :  bool
107                If True, will only remove bristles that are in line of sight of
108                their parent (requires `ncollpyde`). If False, will remove all
109                single-node bristles.
110    inplace :   bool
111                If False will make and return a copy of the skeleton. If True,
112                will modify the `s` inplace.
113
114    Returns
115    -------
116    s :         skeletor.Skeleton
117                Skeleton with single-node twigs removed.
118
119    """
120    if isinstance(mesh, type(None)):
121        mesh = s.mesh
122
123    # Make a copy of the skeleton
124    if not inplace:
125        s = s.copy()
126
127    # Find branch points
128    pcount = s.swc[s.swc.parent_id >= 0].groupby("parent_id").size()
129    bp = pcount[pcount > 1].index
130
131    # Find terminal twigs
132    twigs = s.swc[~s.swc.node_id.isin(s.swc.parent_id)]
133    twigs = twigs[twigs.parent_id.isin(bp)]
134
135    if twigs.empty:
136        return s
137
138    if los_only:
139        coll = make_volume(mesh)
140
141        # Remove twigs that aren't inside the volume
142        twigs = twigs[coll.contains(twigs[["x", "y", "z"]].values)]
143
144        # Generate rays between all pairs and their parents
145        sources = twigs[["x", "y", "z"]].values
146        targets = (
147            s.swc.set_index("node_id").loc[twigs.parent_id, ["x", "y", "z"]].values
148        )
149
150        # Get intersections: `ix` points to index of line segment; `loc` is the
151        #  x/y/z coordinate of the intersection and `is_backface` is True if
152        # intersection happened at the inside of a mesh
153        ix, loc, is_backface = coll.intersections(sources, targets)
154
155        # Find pairs of twigs with no intersection - i.e. with line of sight
156        los = ~np.isin(np.arange(sources.shape[0]), ix)
157
158        # To remove: have line of sight
159        to_remove = twigs[los]
160    else:
161        to_remove = twigs
162
163    s.swc = s.swc[~s.swc.node_id.isin(to_remove.node_id)].copy()
164
165    # Update the mesh map
166    mesh_map = getattr(s, "mesh_map", None)
167    if not isinstance(mesh_map, type(None)):
168        for t in to_remove.itertuples():
169            mesh_map[mesh_map == t.node_id] = t.parent_id
170
171    # Reindex nodes
172    s.reindex(inplace=True)
173
174    return s

177def recenter_vertices(s, mesh=None, inplace=False):
178    """Move nodes that ended up outside the mesh back inside.
179
180    Requires `ncollpyde` (`pip install ncollpyde`).
181
182    Nodes can end up outside the original mesh e.g. if the mesh contraction
183    didn't do a good job (most likely because of internal/degenerate faces that
184    messed up the normals). This function rectifies this by snapping those nodes
185    nodes back to the closest vertex and then tries to move them into the
186    mesh's center. That second step is not guaranteed to work but at least you
187    won't have any more nodes outside the mesh.
188
189    Please note that if connected (!) nodes end up on the same position (i.e
190    because they snapped to the same vertex), we will collapse them.
191
192    Parameters
193    ----------
194    s :         skeletor.Skeleton
195    mesh :      trimesh.Trimesh
196                Original mesh.
197    inplace :   bool
198                If False will make and return a copy of the skeleton. If True,
199                will modify the `s` inplace.
200
201    Returns
202    -------
203    s :         skeletor.Skeleton
204                Skeleton with vertices recentered.
205
206    """
207    if isinstance(mesh, type(None)):
208        mesh = s.mesh
209
210    # Copy skeleton
211    if not inplace:
212        s = s.copy()
213
214    # Find nodes that are outside the mesh
215    coll = make_volume(mesh)
216    outside = ~coll.contains(s.vertices)
217
218    # Skip if all inside
219    if not any(outside):
220        return s
221
222    # For each outside find the closest vertex
223    tree = scipy.spatial.cKDTree(mesh.vertices)
224
225    # Find nodes that are right on top of original vertices
226    dist, ix = tree.query(s.vertices[outside])
227
228    # We don't want to just snap them back to the closest vertex but try to find
229    # the center. For this we will:
230    # 1. Move each vertex inside the mesh by just a bit
231    # 2. Cast a ray along the vertices' normals and find the opposite sides of the mesh
232    # 3. Calculate the distance
233
234    # Get the closest vertex...
235    closest_vertex = mesh.vertices[ix]
236    # .. and offset the vertex positions by just a bit so they "should" be
237    # inside the mesh. In reality that doesn't always happen if the mesh is not
238    # watertight
239    vnormals = mesh.vertex_normals[ix]
240    sources = closest_vertex - vnormals
241
242    # Prepare rays to cast
243    targets = sources - vnormals * 1e4
244
245    # Cast rays
246    hit_ix, loc, is_backface = coll.intersections(sources, targets)
247
248    # center-point if ray hits, otherwise keep closest_vertex
249    final_pos = closest_vertex.copy()
250
251    if len(loc) != 0:
252        # Get half-vector
253        halfvec = np.zeros(sources.shape)
254        halfvec[hit_ix] = (loc - closest_vertex[hit_ix]) / 2
255
256        # Offset vertices
257        candidate = closest_vertex + halfvec
258
259        # Keep only those that are properly inside the mesh
260        now_inside = coll.contains(candidate)
261        final_pos[now_inside] = candidate[now_inside]
262
263    # try harder to get strictly inside
264    still_outside = ~coll.contains(final_pos)
265    if still_outside.any():
266        edge_len = get_edges_unique(mesh, lengths=True)[1]
267        push = (
268            edge_len.mean() / 100
269            if edge_len.size
270            else 1e-4
271        )  # for high-res meshes
272        push = max(push, 1e-6)
273
274        candidate_in = closest_vertex - vnormals * push
275        candidate_out = closest_vertex + vnormals * push
276
277        in_ok = coll.contains(candidate_in)
278        out_ok = coll.contains(candidate_out)
279
280        use_in = still_outside & in_ok
281        use_out = still_outside & ~in_ok & out_ok
282
283        final_pos[use_in] = candidate_in[use_in]
284        final_pos[use_out] = candidate_out[use_out]
285
286    # Keep only those that are properly inside the mesh and fall back to the
287    # closest vertex if that's not the case
288    now_inside = coll.contains(final_pos)
289    final_pos[~now_inside] = closest_vertex[~now_inside]
290
291    # Replace coordinates
292    s.swc.loc[outside, "x"] = final_pos[:, 0]
293    s.swc.loc[outside, "y"] = final_pos[:, 1]
294    s.swc.loc[outside, "z"] = final_pos[:, 2]
295
296    # At this point we may have nodes that snapped to the same vertex and
297    # therefore end up at the same position. We will collapse those nodes
298    # - but only if they are actually connected!
299    # First find duplicate locations
300    u, i, c = np.unique(s.vertices, return_counts=True, return_inverse=True, axis=0)
301
302    # If any coordinates have counter higher 1
303    if c.max() > 1:
304        rewire = {}
305        # Find out which unique coordinates are duplicated
306        dupl = np.where(c > 1)[0]
307
308        # Go over each duplicated coordinate
309        for ix in dupl:
310            # Find the nodes on this duplicate coordinate
311            node_ix = np.where(i == ix)[0]
312
313            # Get their edges
314            edges = s.edges[np.all(np.isin(s.edges, node_ix), axis=1)]
315
316            # We will work on the graph to collapse nodes sequentially A->B->C
317            G = nx.DiGraph()
318            G.add_edges_from(edges)
319            for cc in nx.connected_components(G.to_undirected()):
320                # Root is the node without any outdegree in this subgraph
321                root = [n for n in cc if G.out_degree[n] == 0][0]
322                # We don't want to collapse into `root` because it's not actually
323                # among the nodes with the same coordinates but rather the "last"
324                # nodes parent
325                clps_into = next(G.predecessors(root))
326                # Keep track of how we need to rewire
327                rewire.update({c: clps_into for c in cc if c not in {root, clps_into}})
328
329        # Only mess with the skeleton if there were nodes to be merged
330        if rewire:
331            # Rewire
332            s.swc["parent_id"] = s.swc.parent_id.map(lambda x: rewire.get(x, x))
333
334            # Drop nodes that were collapsed
335            s.swc = s.swc.loc[~s.swc.node_id.isin(rewire)]
336
337            # Update mesh map
338            if not isinstance(s.mesh_map, type(None)):
339                s.mesh_map = [rewire.get(x, x) for x in s.mesh_map]
340
341            # Reindex to make vertex IDs continous again
342            s.reindex(inplace=True)
343
344            # This prevents future SettingsWithCopy Warnings:
345            if not inplace:
346                s.swc = s.swc.copy()
347
348    return s

351def fix_outside_edges(
352    s, mesh=None, inplace=False, max_iter=8, smooth_iters=1, eps=1e-6
353):
354    """Fix edges that cross outside the mesh boundary.
355
356    Requires `ncollpyde` (`pip install ncollpyde`).
357
358    This function detects skeleton edges that intersect the mesh boundary and
359    fixes them by iteratively splitting crossing edges (inserting new vertices
360    along the edge) and then recentering any vertices that end up outside the
361    mesh using `recenter_vertices()`.
362
363    Notes
364    -----
365    This will also modify original vertices positions (via `skeletor.post.recenter_vertices()`).
366    Splitting edges inserts new skeleton nodes that are not represented in "skel_map". Currently,
367    we invalidate any existing "mesh_map" and "skel_map" (setting it to `None`). In the
368    future, we may add functionality to update the mapping instead.
369
370    Parameters
371    ----------
372    s :         skeletor.Skeleton
373    mesh :      trimesh.Trimesh
374                Original mesh. If mesh is None, will use the mesh associated with input
375                skeleton (`s.mesh`).
376    inplace :   bool
377                If False will make and return a copy of the skeleton. If True,
378                will modify the `s` inplace.
379    max_iter :  int
380                Max split iterations for boundary-crossing edges.
381    smooth_iters : int
382                Number of smoothing iterations for degree-2 chain nodes.
383    eps :       float or {'auto'}
384                Ignore intersections within eps of either endpoint.
385                If "auto", uses mesh mean unique edge length * 1e-4.
386
387    Returns
388    -------
389    s :         skeletor.Skeleton
390
391    """
392    if isinstance(mesh, type(None)):
393        mesh = s.mesh
394
395    if mesh is None:
396        raise ValueError(
397            "Mesh is required for fixing outside edges. Please provide a mesh or ensure `s.mesh` is set."
398        )
399
400    if not inplace:
401        s = s.copy()
402
403    if s.swc is None or s.swc.empty:
404        return s
405
406    # Determine eps (scale-aware)
407    if isinstance(eps, str):
408        if eps.lower() != "auto":
409            raise ValueError("Invalid value for `eps`. Must be a number or 'auto'.")
410        try:
411            mean_length = float(np.nanmean(get_edges_unique(mesh, lengths=True)[1]))
412        except Exception:
413            mean_length = np.nan
414        eps = (
415            mean_length * 1e-4
416            if (np.isfinite(mean_length) and mean_length > 0)
417            else 1e-6
418        )
419    else:
420        eps = float(eps)
421
422    max_iter = int(max_iter)
423    if max_iter < 0:
424        raise ValueError("`max_iter` must be >= 0")
425
426    smooth_iters = int(smooth_iters)
427    if smooth_iters < 0:
428        raise ValueError("`smooth_iters` must be >= 0")
429
430    coll = make_volume(mesh)
431
432    # 1. Recenter any nodes outside the mesh
433    if (~coll.contains(s.vertices)).any():
434        recenter_vertices(s, mesh=mesh, inplace=True)
435
436    has_radius = "radius" in s.swc.columns
437
438    # 2. Iteratively split crossing edges
439    for _ in range(max_iter):
440        swc = s.swc
441        edge_rows = np.where(swc.parent_id.values >= 0)[0]
442        if edge_rows.size == 0:
443            break
444
445        sources = swc.loc[edge_rows, ["x", "y", "z"]].values
446        parent_ids = swc.loc[edge_rows, "parent_id"].values
447        targets = swc.set_index("node_id").loc[parent_ids, ["x", "y", "z"]].values
448
449        ix, loc, _ = coll.intersections(sources, targets)
450
451        crossing = np.zeros(edge_rows.shape[0], dtype=bool)
452
453        if len(ix):
454            d_src = np.linalg.norm(loc - sources[ix], axis=1)
455            d_tgt = np.linalg.norm(loc - targets[ix], axis=1)
456            real_crossings = (d_src > eps) & (d_tgt > eps)
457            if real_crossings.any():
458                crossing[np.unique(ix[real_crossings])] = True
459
460        if not crossing.any():
461            break
462
463        to_split = edge_rows[crossing]
464        next_node_id = int(swc.node_id.max()) + 1
465
466        new_rows = []
467        nodes = swc.set_index("node_id")
468
469        # Cache arrays for cheap positional access inside the loop
470        parent_id_arr = swc["parent_id"].to_numpy(copy=True)
471        parent_col = swc.columns.get_loc("parent_id")
472
473        xyz_arr = swc[["x", "y", "z"]].to_numpy(copy=False)
474        if not np.issubdtype(xyz_arr.dtype, np.number):
475            xyz_arr = xyz_arr.astype(float)
476
477        for edge_row in to_split:
478            # edge_row is a positional row index (from np.where)
479            parent_id = int(parent_id_arr[edge_row])
480            if parent_id < 0:
481                continue
482
483            # Midpoint (no projection; recenter will handle)
484            child_co = xyz_arr[edge_row].astype(float, copy=False)
485            parent_co = nodes.loc[parent_id, ["x", "y", "z"]].values.astype(
486                float, copy=False
487            )
488            midpoint = (child_co + parent_co) / 2.0
489
490            # Rewire child -> new node (positional write; avoids boolean mask on node_id)
491            swc.iat[edge_row, parent_col] = next_node_id
492            parent_id_arr[edge_row] = next_node_id  # keep cached view consistent
493
494            # Create new node row
495            row = {col: np.nan for col in swc.columns}
496            row["node_id"] = next_node_id
497            row["parent_id"] = parent_id
498            row["x"], row["y"], row["z"] = midpoint
499
500            if has_radius:
501                child_r = pd.to_numeric(
502                    pd.Series([swc.at[edge_row, "radius"]]), errors="coerce"
503                ).iloc[0]
504                parent_r = pd.to_numeric(
505                    pd.Series([nodes.loc[parent_id, "radius"]]), errors="coerce"
506                ).iloc[0]
507                row["radius"] = np.nanmean(np.array([child_r, parent_r], dtype=float))
508
509            new_rows.append(row)
510            next_node_id += 1
511
512        if not new_rows:
513            break
514
515        swc = pd.concat(
516            [swc, pd.DataFrame(new_rows, columns=swc.columns)], ignore_index=True
517        )
518        s.swc = swc
519
520        coords = s.swc[["x", "y", "z"]].values
521        if (~coll.contains(coords)).any():
522            recenter_vertices(s, mesh=mesh, inplace=True)
523
524    # 3. Smoothing (degree-2 chain nodes), then recenter
525    for _ in range(smooth_iters):
526        swc = s.swc
527
528        child_counts = swc[swc.parent_id >= 0].groupby("parent_id").size()
529        is_chain = (swc.parent_id >= 0) & (
530            swc.node_id.map(child_counts).fillna(0).astype(int) == 1
531        )
532        chain_nodes = swc.loc[is_chain, "node_id"].values.astype(int)
533
534        if chain_nodes.size == 0:
535            break
536
537        only_child = (
538            swc[swc.parent_id >= 0].groupby("parent_id").node_id.first().to_dict()
539        )
540        nodes = swc.set_index("node_id")
541
542        parent_ids = nodes.loc[chain_nodes, "parent_id"].values.astype(int)
543        child_ids = np.array([only_child[n] for n in chain_nodes], dtype=int)
544
545        parent_co = nodes.loc[parent_ids, ["x", "y", "z"]].values.astype(float)
546        child_co = nodes.loc[child_ids, ["x", "y", "z"]].values.astype(float)
547        smoothed = (parent_co + child_co) / 2.0
548
549        swc.loc[is_chain, ["x", "y", "z"]] = smoothed
550        s.swc = swc
551
552        coords = s.swc[["x", "y", "z"]].values
553        if (~coll.contains(coords)).any():
554            recenter_vertices(s, mesh=mesh, inplace=True)
555
556    swc = s.swc
557    edge_rows = np.where(swc.parent_id.values >= 0)[0]
558    remaining_crossings = 0
559    # Detect crossing edges again for double-checking
560    if edge_rows.size:
561        sources = swc.loc[edge_rows, ["x", "y", "z"]].values
562        parent_ids = swc.loc[edge_rows, "parent_id"].values
563        nodes = swc.set_index("node_id")
564        try:
565            targets = nodes.loc[parent_ids, ["x", "y", "z"]].values
566            ix, loc, _ = coll.intersections(sources, targets)
567            if len(ix):
568                d_src = np.linalg.norm(loc - sources[ix], axis=1)
569                d_tgt = np.linalg.norm(loc - targets[ix], axis=1)
570                real = (d_src > eps) & (d_tgt > eps)
571                if np.any(real):
572                    crossing = np.zeros(edge_rows.shape[0], dtype=bool)
573                    crossing[np.unique(ix[real])] = True
574                    remaining_crossings = int(crossing.sum())
575        except KeyError:
576            remaining_crossings = 0
577
578    if remaining_crossings > 0:
579        warnings.warn(
580            f"{remaining_crossings} crossing edges remain after {max_iter} "
581            "fix iteration(s); returning best-effort result. Consider increasing "
582            "`max_iter`, adjusting `eps`, and/or running `post.clean_up` / "
583            "`post.remove_bristles` first. Also check mesh quality (e.g. non-watertight "
584            "or degenerate faces).",
585            RuntimeWarning,
586        )
587
588    # Invalidate mesh_map
589    s.mesh_map = None
590
591    return s