skeletor.post

 29def radii(s, mesh=None, method='knn', aggregate='mean', validate=False, **kwargs):
 30    """Extract radii for given skeleton table.
 31
 32    Important
 33    ---------
 34    This function really only produces useful radii if the skeleton is centered
 35    inside the mesh. `by_wavefront` does that by default whereas all other
 36    skeletonization methods don't. Your best bet to get centered skeletons is
 37    to contract the mesh first (`sk.pre.contract`).
 38
 39    Parameters
 40    ----------
 41    s :         skeletor.Skeleton
 42                Skeleton to clean up.
 43    mesh :      trimesh.Trimesh, optional
 44                Original mesh (e.g. before contraction). If not provided will
 45                use the mesh associated with ``s``.
 46    method :    "knn" | "ray"
 47                Whether and how to add radius information to each node::
 48
 49                    - "knn" uses k-nearest-neighbors to get radii: fast but
 50                      potential for being very wrong
 51                    - "ray" uses ray-casting to get radii: slower but sometimes
 52                      less wrong
 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

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

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

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

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

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

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