skeletor

What is skeletor?

Unlike its namesake, this skeletor does not (yet) seek to conquer Eternia but to turn meshes into skeletons.

Before we get started some terminology:

  • a mesh is something that consists of vertices and faces
  • a skeleton is a (hierarchical) tree-like structure consisting of vertices (also called nodes) and edges that connect them

Skeletons are useful for a range of reasons. For example:

  1. Typically smaller (less vertices) than the mesh
  2. Have an implicit sense of topology (e.g. "this node is distal to that node")

Extracting skeletons from meshes (or other types of data such as voxels) is non-trivial and there are a great many research papers exploring various different approaches (see Google scholar).

skeletor implements some algorithms that I found useful in my work with neurons. In my experience there is unfortuntely no magic bullet when it comes to skeletonization and chances are you will have to fiddle around a bit to get decent results.

Installation

From PyPI:

pip3 install skeletor

For the bleeding-edge version from Github:

pip3 install git+https://github.com/navis-org/skeletor@master

Getting started

A skeletonization pipeline typically consists of:

  1. Some pre-processing of the mesh (e.g. fixing some potential errors like degenerate faces, unreferenced vertices, etc.)
  2. The skeletonization itself
  3. Some post-processing of the skeleton (e.g. adding radius information, smoothing, etc.)

Here is a complete list of available functions:

function description
example data
skeletor.example_mesh() load an example mesh
pre-processing
skeletor.pre.fix_mesh() fix some common errors found in meshes
skeletor.pre.remesh() re-generate mesh (uses Blender 3D)
skeletor.pre.simplify() reduce mesh complexity (uses Blender 3D)
skeletor.pre.contract() contract mesh to facilitate skeletonization [1]
skeletonization
skeletor.skeletonize.by_wavefront() very fast, works well for tubular meshes (like neurons)
skeletor.skeletonize.by_vertex_clusters() very fast but needs mesh to be contracted (see above)
skeletor.skeletonize.by_edge_collapse() presented in [1] but never got this to work well
skeletor.skeletonize.by_teasar() very fast and robust, works on mesh surface
skeletor.skeletonize.by_tangent_ball() very fast, best on smooth meshes
postprocessing
skeletor.post.clean_up() fix some potential errors in the skeleton
skeletor.post.radii() add radius information using various method
skeletor.post.smooth() smooth the skeleton
skeletor.post.remove_bristles() remove single-node bristles from the skeleton
skeletor.post.despike() smooth out jumps in the skeleton

See docstrings of the respective functions for details.

A pipeline might look like this:

  1. skeletor.pre.fix_mesh() to fix the mesh
  2. skeletor.pre.simplify() to simplify the mesh
  3. skeletor.pre.contract() to contract the mesh [1]
  4. skeletor.skeletonize.by_vertex_clusters() to generate a skeleton
  5. skeletor.post.clean_up() to clean up some potential issues with the skeleton
  6. skeletor.post.smooth() to smooth the skeleton
  7. skeletor.post.radii() to extract radii either by k-nearest neighbours or ray-casting

In my experience there is no one-size-fits-all. You will have to play around to find the right approach and parameters to get nice skeletons for your meshes. If you need help just open an issue.

Also check out the Gotchas below!

Examples

First load the example mesh (a fruit fly neuron):

>>> import skeletor as sk
>>> mesh = sk.example_mesh()
>>> mesh
<trimesh.Trimesh(vertices.shape=(6582, 3), faces.shape=(13772, 3))>

Next see if there is stuff to fix in the mesh (degenerate faces, duplicate vertices, etc.):

>>> fixed = sk.pre.fix_mesh(mesh, remove_disconnected=5, inplace=False)
>>> fixed
<trimesh.Trimesh(vertices.shape=(6213, 3), faces.shape=(12805, 3))>

Now for tubular meshes like this neuron, the "wave front" skeletonization method performs really well: it works by casting waves across the mesh and collapsing the resulting rings into a skeleton (kinda like when you throw a stone in a pond and track the expanding ripples).

>>> skel = sk.skeletonize.by_wavefront(fixed, waves=1, step_size=1)
>>> skel
<Skeleton(vertices=(1258, 3), edges=(1194, 2), method=wavefront)>

All skeletonization methods return a Skeleton object. These are just convenient objects to bundle the various outputs of the skeletonization.

>>> # x/y/z location of skeleton vertices (nodes)
>>> skel.vertices
array([[16744, 36720, 26407],
       ...,
       [22076, 23217, 24472]])
>>> # child -> parent edges
>>> skel.edges
array([[  64,   31],
       ...,
       [1257, 1252]])
>>> # Mapping for mesh to skeleton vertex indices
>>> skel.mesh_map
array([ 157,  158, 1062, ...,  525,  474,  547])
>>> # SWC table
>>> skel.swc.head()
   node_id  parent_id             x             y             z    radius
0        0         -1  16744.005859  36720.058594  26407.902344  0.000000
1        1         -1   5602.751953  22266.756510  15799.991211  7.542587
2        2         -1  16442.666667  14999.978516  10887.916016  5.333333

SWC is a commonly used format for saving skeletons. Skeleton objects have a method for quickly saving a correctly formatted SWC file:

>>> skel.save_swc('~/Documents/my_skeleton.swc')

If you installed pyglet (see above) you can also use trimesh's plotting capabilities to inspect the results:

>>> skel.show(mesh=True)

skeletor_example

That looks pretty good already but let's run some pro-forma postprocessing.

>>> sk.post.clean_up(skel, inplace=True)
<Skeleton(vertices=(1071, 3), edges=(1070, 2))>

So that would be a full pipeline mesh to skeleton. Don't expect your own meshes to produce such nice results off the bat though. Chances are you will need to play around to find the right recipe. If you don't know where to start, I suggest you try out mesh contraction + vertex clustering first:

>>> import skeletor as sk
>>> # Load the example mesh that ships with skeletor
>>> mesh = sk.example_mesh()
>>> # Alternatively use trimesh to load/construct your own mesh:
>>> # import trimesh as tm
>>> # mesh = tm.Trimesh(vertices, faces)
>>> # mesh = tm.load_mesh('some/mesh.obj')
>>> # Run some general clean-up (see docstring for details)
>>> fixed = sk.pre.fix_mesh(mesh, remove_disconnected=5, inplace=False)
>>> # Contract mesh to 10% (0.1) of original volume
>>> cont = sk.pre.contract(fixed, epsilon=0.1)
>>> # Skeletonize
>>> skel = sk.skeletonize.by_vertex_clusters(cont, sampling_dist=100)
>>> # Replace contracted mesh with original for postprocessing and plotting
>>> skel.mesh = fixed
>>> # Add radii (vertex cluster method does not do that automatically)
>>> sk.post.radii(skel, method='knn')
>>> skel.show(mesh=True)

Gotchas

  • while this is a general purpose library, my personal focus is on neurons and this has certainly influenced things like default parameter values and certain post-processing steps
  • meshes need to be triangular (we are using trimesh)
  • use sk.pre.simplify if your mesh is very complex (half a million vertices is where things start getting sluggish)
  • a good mesh contraction is often half the battle but it can be tricky to get to work
  • if the mesh consists of multiple disconnected pieces the skeleton will likewise be fragmented (i.e. will have multiple roots)
  • it's often a good idea to fix issues with the skeleton in postprocessing rather than trying to get the skeletonization to be perfect

Benchmarks

skeletor_benchmark

Benchmarks were run on a 2018 MacBook Pro (2.2 GHz Core i7, 32Gb memory) with optional fastremap dependency installed. Note some of these functions (e.g. contraction and TEASAR/vertex cluster skeletonization) vary a lot in speed based on parameterization.

What about algorithm X?

skeletor contains some algorithms that I found easy enough to implement and useful for my work with neurons. If you have some interesting paper/approach that could make a nice addition to skeletor, please get in touch on Github. Pull requests are always welcome!

References

[1] Au OK, Tai CL, Chu HK, Cohen-Or D, Lee TY. Skeleton extraction by mesh contraction. ACM Transactions on Graphics (TOG). 2008 Aug 1;27(3):44.

The abstract and the paper can be found here. Also see this YouTube video.

Some of the code in skeletor was modified from the Py_BL_MeshSkeletonization addon created by #0K Srinivasan Ramachandran and published under GPL3.

Top-level functions and classes

At top-level we only expose example_mesh() and the Skeleton class (which you probably won't ever need to touch manually). Everything else is neatly tucked away into submodules (see side-bar or above table).

  1#    This script is part of skeletor (http://www.github.com/navis-org/skeletor).
  2#    Copyright (C) 2018 Philipp Schlegel
  3#
  4#    This program is free software: you can redistribute it and/or modify
  5#    it under the terms of the GNU General Public License as published by
  6#    the Free Software Foundation, either version 3 of the License, or
  7#    (at your option) any later version.
  8#
  9#    This program is distributed in the hope that it will be useful,
 10#    but WITHOUT ANY WARRANTY; without even the implied warranty of
 11#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 12#    GNU General Public License for more details.
 13#
 14#    You should have received a copy of the GNU General Public License
 15#    along with this program.
 16
 17r"""
 18# What is skeletor?
 19
 20Unlike its [namesake](https://en.wikipedia.org/wiki/Skeletor), this `skeletor`
 21does not (yet) seek to conquer Eternia but to turn meshes into skeletons.
 22
 23Before we get started some terminology:
 24
 25- a *mesh* is something that consists of *vertices* and *faces*
 26- a *skeleton* is a (hierarchical) tree-like structure consisting of *vertices*
 27  (also called *nodes*) and *edges* that connect them
 28
 29Skeletons are useful for a range of reasons. For example:
 30
 311. Typically smaller (less vertices) than the mesh
 322. Have an implicit sense of topology (e.g. "this node is distal to that node")
 33
 34Extracting skeletons from meshes (or other types of data such as voxels) is
 35non-trivial and there are a great many research papers exploring various
 36different approaches (see [Google scholar](https://scholar.google.com/scholar
 37?hl=en&as_sdt=0%2C5&q=skeleton+extraction&btnG=)).
 38
 39`skeletor` implements some algorithms that I found useful in my work with
 40neurons. In my experience there is unfortuntely no magic bullet when it
 41comes to skeletonization and chances are you will have to fiddle around a bit
 42to get decent results.
 43
 44# Installation
 45
 46From PyPI:
 47```bash
 48pip3 install skeletor
 49```
 50
 51For the bleeding-edge version from Github:
 52```bash
 53pip3 install git+https://github.com/navis-org/skeletor@master
 54```
 55
 56# Getting started
 57
 58A skeletonization pipeline typically consists of:
 59
 601. Some pre-processing of the mesh (e.g. fixing some potential errors like
 61   degenerate faces, unreferenced vertices, etc.)
 622. The skeletonization itself
 633. Some post-processing of the skeleton (e.g. adding radius information, smoothing, etc.)
 64
 65------
 66
 67Here is a complete list of available functions:
 68
 69| function                                    | description                                                 |
 70| ------------------------------------------- | ----------------------------------------------------------- |
 71| **example data**                            |                                                             |
 72| `skeletor.example_mesh()`                   | load an example mesh                                        |
 73| **pre-processing**                          |                                                             |
 74| `skeletor.pre.fix_mesh()`                   | fix some common errors found in meshes                      |
 75| `skeletor.pre.remesh()`                     | re-generate mesh (uses Blender 3D)                          |
 76| `skeletor.pre.simplify()`                   | reduce mesh complexity (uses Blender 3D)                    |
 77| `skeletor.pre.contract()`                   | contract mesh to facilitate skeletonization [1]             |
 78| **skeletonization**                         |                                                             |
 79| `skeletor.skeletonize.by_wavefront()`       | very fast, works well for tubular meshes (like neurons)     |
 80| `skeletor.skeletonize.by_vertex_clusters()` | very fast but needs mesh to be contracted (see above)       |
 81| `skeletor.skeletonize.by_edge_collapse()`   | presented in [1] but never got this to work well            |
 82| `skeletor.skeletonize.by_teasar()`          | very fast and robust, works on mesh surface                 |
 83| `skeletor.skeletonize.by_tangent_ball()`    | very fast, best on smooth meshes                            |
 84| **postprocessing**                          |                                                             |
 85| `skeletor.post.clean_up()`                  | fix some potential errors in the skeleton                   |
 86| `skeletor.post.radii()`                     | add radius information using various method                 |
 87| `skeletor.post.smooth()`                    | smooth the skeleton                                         |
 88| `skeletor.post.remove_bristles()`           | remove single-node bristles from the skeleton               |
 89| `skeletor.post.despike()`                   | smooth out jumps in the skeleton                            |
 90
 91------
 92
 93See docstrings of the respective functions for details.
 94
 95A pipeline might look like this:
 96
 97 1. `skeletor.pre.fix_mesh()` to fix the mesh
 98 2. `skeletor.pre.simplify()` to simplify the mesh
 99 3. `skeletor.pre.contract()` to contract the mesh [1]
100 4. `skeletor.skeletonize.by_vertex_clusters()` to generate a skeleton
101 5. `skeletor.post.clean_up()` to clean up some potential issues with the skeleton
102 6. `skeletor.post.smooth()` to smooth the skeleton
103 7. `skeletor.post.radii()` to extract radii either by k-nearest neighbours or ray-casting
104
105In my experience there is no one-size-fits-all. You will have to play around to
106find the right approach and parameters to get nice skeletons for your meshes.
107If you need help just open an [issue](https://github.com/navis-org/skeletor/issues).
108
109Also check out the Gotchas below!
110
111# Examples
112
113First load the example mesh (a fruit fly neuron):
114
115```Python
116>>> import skeletor as sk
117>>> mesh = sk.example_mesh()
118>>> mesh
119<trimesh.Trimesh(vertices.shape=(6582, 3), faces.shape=(13772, 3))>
120```
121
122Next see if there is stuff to fix in the mesh (degenerate faces, duplicate
123vertices, etc.):
124
125```Python
126>>> fixed = sk.pre.fix_mesh(mesh, remove_disconnected=5, inplace=False)
127>>> fixed
128<trimesh.Trimesh(vertices.shape=(6213, 3), faces.shape=(12805, 3))>
129```
130
131Now for tubular meshes like this neuron, the "wave front" skeletonization method
132performs really well: it works by casting waves across the mesh and collapsing
133the resulting rings into a skeleton (kinda like when you throw a stone in a
134pond and track the expanding ripples).
135
136```Python
137>>> skel = sk.skeletonize.by_wavefront(fixed, waves=1, step_size=1)
138>>> skel
139<Skeleton(vertices=(1258, 3), edges=(1194, 2), method=wavefront)>
140```
141
142All skeletonization methods return a `Skeleton` object. These are just
143convenient objects to bundle the various outputs of the skeletonization.
144
145```Python
146>>> # x/y/z location of skeleton vertices (nodes)
147>>> skel.vertices
148array([[16744, 36720, 26407],
149       ...,
150       [22076, 23217, 24472]])
151>>> # child -> parent edges
152>>> skel.edges
153array([[  64,   31],
154       ...,
155       [1257, 1252]])
156>>> # Mapping for mesh to skeleton vertex indices
157>>> skel.mesh_map
158array([ 157,  158, 1062, ...,  525,  474,  547])
159>>> # SWC table
160>>> skel.swc.head()
161   node_id  parent_id             x             y             z    radius
1620        0         -1  16744.005859  36720.058594  26407.902344  0.000000
1631        1         -1   5602.751953  22266.756510  15799.991211  7.542587
1642        2         -1  16442.666667  14999.978516  10887.916016  5.333333
165```
166
167SWC is a commonly used format for saving skeletons. `Skeleton` objects
168have a method for quickly saving a correctly formatted SWC file:
169
170```Python
171>>> skel.save_swc('~/Documents/my_skeleton.swc')
172```
173
174If you installed `pyglet` (see above) you can also use `trimesh`'s plotting
175capabilities to inspect the results:
176
177```Python
178>>> skel.show(mesh=True)
179```
180
181<img src="https://github.com/navis-org/skeletor/raw/master/_static/example1.png" alt="skeletor_example" width="100%"/>
182
183That looks pretty good already but let's run some pro-forma postprocessing.
184
185```Python
186>>> sk.post.clean_up(skel, inplace=True)
187<Skeleton(vertices=(1071, 3), edges=(1070, 2))>
188```
189
190So that would be a full pipeline mesh to skeleton. Don't expect your own meshes
191to produce such nice results off the bat though. Chances are you will need to
192play around to find the right recipe. If you don't know where to start, I suggest
193you try out mesh contraction + vertex clustering first:
194
195```Python
196>>> import skeletor as sk
197>>> # Load the example mesh that ships with skeletor
198>>> mesh = sk.example_mesh()
199>>> # Alternatively use trimesh to load/construct your own mesh:
200>>> # import trimesh as tm
201>>> # mesh = tm.Trimesh(vertices, faces)
202>>> # mesh = tm.load_mesh('some/mesh.obj')
203>>> # Run some general clean-up (see docstring for details)
204>>> fixed = sk.pre.fix_mesh(mesh, remove_disconnected=5, inplace=False)
205>>> # Contract mesh to 10% (0.1) of original volume
206>>> cont = sk.pre.contract(fixed, epsilon=0.1)
207>>> # Skeletonize
208>>> skel = sk.skeletonize.by_vertex_clusters(cont, sampling_dist=100)
209>>> # Replace contracted mesh with original for postprocessing and plotting
210>>> skel.mesh = fixed
211>>> # Add radii (vertex cluster method does not do that automatically)
212>>> sk.post.radii(skel, method='knn')
213>>> skel.show(mesh=True)
214```
215
216# Gotchas
217
218- while this is a general purpose library, my personal focus is on neurons and
219  this has certainly influenced things like default parameter values and certain
220  post-processing steps
221- meshes need to be triangular (we are using `trimesh`)
222- use `sk.pre.simplify` if your mesh is very complex (half a million vertices is
223  where things start getting sluggish)
224- a good mesh contraction is often half the battle but it can be tricky to get
225  to work
226- if the mesh consists of multiple disconnected pieces the skeleton will
227  likewise be fragmented (i.e. will have multiple roots)
228- it's often a good idea to fix issues with the skeleton in postprocessing rather
229  than trying to get the skeletonization to be perfect
230
231# Benchmarks
232
233<img src="https://github.com/navis-org/skeletor/raw/master/benchmarks/benchmark_2.png" alt="skeletor_benchmark" width="100%"/>
234
235[Benchmarks](https://github.com/navis-org/skeletor/blob/master/benchmarks/skeletor_benchmark.ipynb)
236were run on a 2018 MacBook Pro (2.2 GHz Core i7, 32Gb memory) with optional
237`fastremap` dependency installed. Note some of these functions (e.g.
238contraction and TEASAR/vertex cluster skeletonization) vary a lot in
239speed based on parameterization.
240
241# What about algorithm `X`?
242
243`skeletor` contains some algorithms that I found easy enough to implement
244and useful for my work with neurons. If you have some interesting paper/approach
245that could make a nice addition to `skeletor`, please get in touch on Github.
246Pull requests are always welcome!
247
248# References
249
250`[1] Au OK, Tai CL, Chu HK, Cohen-Or D, Lee TY. Skeleton extraction by mesh contraction. ACM Transactions on Graphics (TOG). 2008 Aug 1;27(3):44.`
251
252The abstract and the paper can be found [here](http://visgraph.cse.ust.hk/projects/skeleton/).
253Also see [this](https://www.youtube.com/watch?v=-H7n59YQCRM&feature=youtu.be) YouTube video.
254
255Some of the code in skeletor was modified from the
256[Py_BL_MeshSkeletonization](https://github.com/aalavandhaann/Py_BL_MeshSkeletonization)
257addon created by #0K Srinivasan Ramachandran and published under GPL3.
258
259# Top-level functions and classes
260At top-level we only expose `example_mesh()` and the `Skeleton` class (which
261you probably won't ever need to touch manually). Everything else is neatly
262tucked away into submodules (see side-bar or above table).
263
264"""
265
266__version__ = "1.3.0"
267__version_vector__ = (1, 3, 0)
268
269from . import skeletonize
270from . import pre
271from . import post
272
273from .skeletonize.base import Skeleton
274from .data import example_mesh
275
276__docformat__ = "numpy"
277
278__all__ = ['Skeleton', 'example_mesh', 'pre', 'post', 'skeletonize']
class Skeleton: 
 33class Skeleton:
 34    """Class representing a skeleton.
 35
 36    Typically returned as results from a skeletonization.
 37
 38    Attributes
 39    ----------
 40    swc :       pd.DataFrame, optional
 41                SWC table.
 42    vertices :  (N, 3) array
 43                Vertex (node) positions.
 44    edges :     (M, 2) array
 45                Indices of connected vertex pairs.
 46    radii :     (N, ) array, optional
 47                Radii for each vertex (node) in the skeleton.
 48    mesh :      trimesh, optional
 49                The original mesh.
 50    mesh_map :  array, optional
 51                Same length as ``mesh``. Maps mesh vertices to vertices (nodes)
 52                in the skeleton.
 53    skel_map :  array of arrays, optional
 54                Inverse of `mesh_map`: maps skeleton vertices (nodes) to mesh
 55                vertices.
 56    method :    str, optional
 57                Which method was used to generate the skeleton.
 58
 59    """
 60
 61    def __init__(self, swc, mesh=None, mesh_map=None, method=None):
 62        self.swc = swc
 63        self.mesh = mesh
 64        self.mesh_map = mesh_map
 65        self.method = method
 66
 67    def __str__(self):
 68        """Summary."""
 69        return self.__repr__()
 70
 71    def __repr__(self):
 72        """Return quick summary of the skeleton's geometry."""
 73        elements = []
 74        if hasattr(self, 'vertices'):
 75            elements.append(f'vertices={self.vertices.shape}')
 76        if hasattr(self, 'edges'):
 77            elements.append(f'edges={self.edges.shape}')
 78        if hasattr(self, 'method'):
 79            elements.append(f'method={self.method}')
 80        return f'<Skeleton({", ".join(elements)})>'
 81
 82    @property
 83    def edges(self):
 84        """Return skeleton edges."""
 85        return self.swc.loc[self.swc.parent_id >= 0,
 86                            ['node_id', 'parent_id']].values
 87
 88    @property
 89    def vertices(self):
 90        """Return skeleton vertices (nodes)."""
 91        return self.swc[['x', 'y', 'z']].values
 92
 93    @property
 94    def radius(self):
 95        """Return radii."""
 96        if 'radius' not in self.swc.columns:
 97            raise ValueError('No radius info found. Run `skeletor.post.radii()`'
 98                             ' to get them.')
 99        return self.swc['radius'].values,
100
101    @property
102    def skeleton(self):
103        """Skeleton as trimesh Path3D."""
104        if not hasattr(self, '_skeleton'):
105            lines = [tm.path.entities.Line(e) for e in self.edges]
106
107            self._skeleton = tm.path.Path3D(entities=lines,
108                                            vertices=self.vertices,
109                                            process=False)
110        return self._skeleton
111
112    @property
113    def skel_map(self):
114        """Skeleton vertex (nodes) to mesh vertices. Based on `mesh_map`."""
115        if isinstance(self.mesh_map, type(None)):
116            return None
117        return pd.DataFrame(self.mesh_map
118                            ).reset_index(drop=False
119                                          ).groupby(0)['index'].apply(np.array
120                                                                      ).values
121
122    @property
123    def leafs(self):
124        """Leaf nodes (includes root)."""
125        swc = self.swc
126        leafs = swc[~swc.node_id.isin(swc.parent_id.values) | (swc.parent_id < 0)]
127        return leafs.copy()
128
129    def reindex(self, inplace=False):
130        """Clean up skeleton."""
131        x = self
132        if not inplace:
133            x = x.copy()
134
135        # Re-index to make node IDs continous again
136        x.swc, new_ids = reindex_swc(x.swc)
137
138        # Update mesh map
139        if not isinstance(x.mesh_map, type(None)):
140            x.mesh_map = np.array([new_ids.get(i, i) for i in x.mesh_map])
141
142        if not inplace:
143            return x
144
145    def copy(self):
146        """Return copy of the skeleton."""
147        return Skeleton(swc=self.swc.copy() if not isinstance(self.swc, type(None)) else None,
148                        mesh=self.mesh.copy() if not isinstance(self.mesh, type(None)) else None,
149                        mesh_map=self.mesh_map.copy() if not isinstance(self.mesh_map, type(None)) else None)
150
151    def get_graph(self):
152        """Generate networkX representation of the skeletons.
153
154        Distance between nodes will be used as edge weights.
155
156        Returns
157        -------
158        networkx.DiGraph
159
160        """
161        not_root = self.swc.parent_id >= 0
162        nodes = self.swc.loc[not_root]
163        parents = self.swc.set_index('node_id').loc[self.swc.loc[not_root, 'parent_id'].values]
164
165        dists = nodes[['x', 'y', 'z']].values - parents[['x', 'y', 'z']].values
166        dists = np.sqrt((dists ** 2).sum(axis=1))
167
168        G = nx.DiGraph()
169        G.add_nodes_from(self.swc.node_id.values)
170        G.add_weighted_edges_from(zip(nodes.node_id.values, nodes.parent_id.values, dists))
171
172        return G
173
174    def get_segments(self,
175                     weight = 'weight',
176                     return_lengths = False):
177        """Generate a list of linear segments while maximizing segment lengths.
178
179        Parameters
180        ----------
181        weight :    'weight' | None, optional
182                    If ``"weight"`` use physical, geodesic length to determine
183                    segment length. If ``None`` use number of nodes (faster).
184        return_lengths : bool
185                    If True, also return lengths of segments according to ``weight``.
186
187        Returns
188        -------
189        segments :  list
190                    Segments as list of lists containing node IDs. List is
191                    sorted by segment lengths.
192        lengths :   list
193                    Length for each segment according to ``weight``. Only provided
194                    if `return_lengths` is True.
195
196        """
197        assert weight in ('weight', None), f'Unable to use weight "{weight}"'
198
199        # Get graph representation
200        G = self.get_graph()
201
202        # Get distances to root
203        dists = {}
204        for root in self.swc[self.swc.parent_id < 0].node_id.values:
205            dists.update(nx.shortest_path_length(G, target=root, weight=weight))
206
207        # Sort leaf nodes
208        endNodeIDs = self.leafs[self.leafs.parent_id >= 0].node_id.values
209        endNodeIDs = sorted(endNodeIDs, key=lambda x: dists.get(x, 0), reverse=True)
210
211        seen: set = set()
212        sequences = []
213        for nodeID in endNodeIDs:
214            sequence = [nodeID]
215            parents = list(G.successors(nodeID))
216            while True:
217                if not parents:
218                    break
219                parentID = parents[0]
220                sequence.append(parentID)
221                if parentID in seen:
222                    break
223                seen.add(parentID)
224                parents = list(G.successors(parentID))
225
226            if len(sequence) > 1:
227                sequences.append(sequence)
228
229        # Sort sequences by length
230        lengths = [dists[s[0]] - dists[s[-1]] for s in sequences]
231        sequences = [x for _, x in sorted(zip(lengths, sequences), reverse=True)]
232
233        if return_lengths:
234            return sequences, sorted(lengths, reverse=True)
235        else:
236            return sequences
237
238    def save_swc(self, filepath):
239        """Save skeleton in SWC format.
240
241        Parameters
242        ----------
243        filepath :      path-like
244                        Filepath to save SWC to.
245
246        """
247        header = dedent(f"""\
248        # SWC format file
249        # based on specifications at http://www.neuronland.org/NLMorphologyConverter/MorphologyFormats/SWC/Spec.html
250        # Created on {datetime.date.today()} using skeletor (https://github.com/navis-org/skeletor)
251        # PointNo Label X Y Z Radius Parent
252        # Labels:
253        # 0 = undefined, 1 = soma, 5 = fork point, 6 = end point
254        """)
255
256        # Make copy of SWC table
257        swc = self.swc.copy()
258
259        # Set all labels to undefined
260        swc['label'] = 0
261        swc.loc[~swc.node_id.isin(swc.parent_id.values), 'label'] = 6
262        n_childs = swc.groupby('parent_id').size()
263        bp = n_childs[n_childs > 1].index.values
264        swc.loc[swc.node_id.isin(bp), 'label'] = 5
265
266        # Add radius if missing
267        if 'radius' not in swc.columns:
268            swc['radius'] = 0
269
270        # Get things in order
271        swc = swc[['node_id', 'label', 'x', 'y', 'z', 'radius', 'parent_id']]
272
273        # Adjust column titles
274        swc.columns = ['PointNo', 'Label', 'X', 'Y', 'Z', 'Radius', 'Parent']
275
276        with open(filepath, 'w') as file:
277            # Write header
278            file.write(header)
279
280            # Write data
281            writer = csv.writer(file, delimiter=' ')
282            writer.writerows(swc.astype(str).values)
283
284    def scene(self, mesh=False, **kwargs):
285        """Return a Scene object containing the skeleton.
286
287        Returns
288        -------
289        scene :     trimesh.scene.scene.Scene
290                    Contains the skeleton and optionally the mesh.
291
292        """
293        if mesh:
294            if isinstance(self.mesh, type(None)):
295                raise ValueError('Skeleton has no mesh.')
296
297            self.mesh.visual.face_colors = [100, 100, 100, 100]
298
299            # Note the copy(): without it the transform in show() changes
300            # the original meshes
301            sc = tm.Scene([self.mesh.copy(), self.skeleton.copy()], **kwargs)
302        else:
303            sc = tm.Scene(self.skeleton.copy(), **kwargs)
304
305        return sc
306
307    def show(self, mesh=False, **kwargs):
308        """Render the skeleton in an opengl window. Requires pyglet.
309
310        Parameters
311        ----------
312        mesh :      bool
313                    If True, will render transparent mesh on top of the
314                    skeleton.
315
316        Returns
317        --------
318        scene :     trimesh.scene.Scene
319                    Scene with skeleton in it.
320
321        """
322        scene = self.scene(mesh=mesh)
323
324        # I encountered some issues if object space is big and the easiest
325        # way to work around this is to apply a transform such that the
326        # coordinates have -5 to +5 bounds
327        fac = 5 / np.fabs(self.skeleton.bounds).max()
328        scene.apply_transform(np.diag([fac, fac, fac, 1]))
329
330        return scene.show(**kwargs)

Class representing a skeleton.

Typically returned as results from a skeletonization.

Attributes
  • swc (pd.DataFrame, optional): SWC table.
  • vertices ((N, 3) array): Vertex (node) positions.
  • edges ((M, 2) array): Indices of connected vertex pairs.
  • radii ((N, ) array, optional): Radii for each vertex (node) in the skeleton.
  • mesh (trimesh, optional): The original mesh.
  • mesh_map (array, optional): Same length as mesh. Maps mesh vertices to vertices (nodes) in the skeleton.
  • skel_map (array of arrays, optional): Inverse of mesh_map: maps skeleton vertices (nodes) to mesh vertices.
  • method (str, optional): Which method was used to generate the skeleton.
Skeleton(swc, mesh=None, mesh_map=None, method=None) 
61    def __init__(self, swc, mesh=None, mesh_map=None, method=None):
62        self.swc = swc
63        self.mesh = mesh
64        self.mesh_map = mesh_map
65        self.method = method
swc
mesh
mesh_map
method
edges 
82    @property
83    def edges(self):
84        """Return skeleton edges."""
85        return self.swc.loc[self.swc.parent_id >= 0,
86                            ['node_id', 'parent_id']].values

Return skeleton edges.

vertices 
88    @property
89    def vertices(self):
90        """Return skeleton vertices (nodes)."""
91        return self.swc[['x', 'y', 'z']].values

Return skeleton vertices (nodes).

radius 
93    @property
94    def radius(self):
95        """Return radii."""
96        if 'radius' not in self.swc.columns:
97            raise ValueError('No radius info found. Run `skeletor.post.radii()`'
98                             ' to get them.')
99        return self.swc['radius'].values,

Return radii.

skeleton 
101    @property
102    def skeleton(self):
103        """Skeleton as trimesh Path3D."""
104        if not hasattr(self, '_skeleton'):
105            lines = [tm.path.entities.Line(e) for e in self.edges]
106
107            self._skeleton = tm.path.Path3D(entities=lines,
108                                            vertices=self.vertices,
109                                            process=False)
110        return self._skeleton

Skeleton as trimesh Path3D.

skel_map 
112    @property
113    def skel_map(self):
114        """Skeleton vertex (nodes) to mesh vertices. Based on `mesh_map`."""
115        if isinstance(self.mesh_map, type(None)):
116            return None
117        return pd.DataFrame(self.mesh_map
118                            ).reset_index(drop=False
119                                          ).groupby(0)['index'].apply(np.array
120                                                                      ).values

Skeleton vertex (nodes) to mesh vertices. Based on mesh_map.

leafs 
122    @property
123    def leafs(self):
124        """Leaf nodes (includes root)."""
125        swc = self.swc
126        leafs = swc[~swc.node_id.isin(swc.parent_id.values) | (swc.parent_id < 0)]
127        return leafs.copy()

Leaf nodes (includes root).

def reindex(self, inplace=False): 
129    def reindex(self, inplace=False):
130        """Clean up skeleton."""
131        x = self
132        if not inplace:
133            x = x.copy()
134
135        # Re-index to make node IDs continous again
136        x.swc, new_ids = reindex_swc(x.swc)
137
138        # Update mesh map
139        if not isinstance(x.mesh_map, type(None)):
140            x.mesh_map = np.array([new_ids.get(i, i) for i in x.mesh_map])
141
142        if not inplace:
143            return x

Clean up skeleton.

def copy(self): 
145    def copy(self):
146        """Return copy of the skeleton."""
147        return Skeleton(swc=self.swc.copy() if not isinstance(self.swc, type(None)) else None,
148                        mesh=self.mesh.copy() if not isinstance(self.mesh, type(None)) else None,
149                        mesh_map=self.mesh_map.copy() if not isinstance(self.mesh_map, type(None)) else None)

Return copy of the skeleton.

def get_graph(self): 
151    def get_graph(self):
152        """Generate networkX representation of the skeletons.
153
154        Distance between nodes will be used as edge weights.
155
156        Returns
157        -------
158        networkx.DiGraph
159
160        """
161        not_root = self.swc.parent_id >= 0
162        nodes = self.swc.loc[not_root]
163        parents = self.swc.set_index('node_id').loc[self.swc.loc[not_root, 'parent_id'].values]
164
165        dists = nodes[['x', 'y', 'z']].values - parents[['x', 'y', 'z']].values
166        dists = np.sqrt((dists ** 2).sum(axis=1))
167
168        G = nx.DiGraph()
169        G.add_nodes_from(self.swc.node_id.values)
170        G.add_weighted_edges_from(zip(nodes.node_id.values, nodes.parent_id.values, dists))
171
172        return G

Generate networkX representation of the skeletons.

Distance between nodes will be used as edge weights.

Returns
  • networkx.DiGraph
def get_segments(self, weight='weight', return_lengths=False): 
174    def get_segments(self,
175                     weight = 'weight',
176                     return_lengths = False):
177        """Generate a list of linear segments while maximizing segment lengths.
178
179        Parameters
180        ----------
181        weight :    'weight' | None, optional
182                    If ``"weight"`` use physical, geodesic length to determine
183                    segment length. If ``None`` use number of nodes (faster).
184        return_lengths : bool
185                    If True, also return lengths of segments according to ``weight``.
186
187        Returns
188        -------
189        segments :  list
190                    Segments as list of lists containing node IDs. List is
191                    sorted by segment lengths.
192        lengths :   list
193                    Length for each segment according to ``weight``. Only provided
194                    if `return_lengths` is True.
195
196        """
197        assert weight in ('weight', None), f'Unable to use weight "{weight}"'
198
199        # Get graph representation
200        G = self.get_graph()
201
202        # Get distances to root
203        dists = {}
204        for root in self.swc[self.swc.parent_id < 0].node_id.values:
205            dists.update(nx.shortest_path_length(G, target=root, weight=weight))
206
207        # Sort leaf nodes
208        endNodeIDs = self.leafs[self.leafs.parent_id >= 0].node_id.values
209        endNodeIDs = sorted(endNodeIDs, key=lambda x: dists.get(x, 0), reverse=True)
210
211        seen: set = set()
212        sequences = []
213        for nodeID in endNodeIDs:
214            sequence = [nodeID]
215            parents = list(G.successors(nodeID))
216            while True:
217                if not parents:
218                    break
219                parentID = parents[0]
220                sequence.append(parentID)
221                if parentID in seen:
222                    break
223                seen.add(parentID)
224                parents = list(G.successors(parentID))
225
226            if len(sequence) > 1:
227                sequences.append(sequence)
228
229        # Sort sequences by length
230        lengths = [dists[s[0]] - dists[s[-1]] for s in sequences]
231        sequences = [x for _, x in sorted(zip(lengths, sequences), reverse=True)]
232
233        if return_lengths:
234            return sequences, sorted(lengths, reverse=True)
235        else:
236            return sequences

Generate a list of linear segments while maximizing segment lengths.

Parameters
  • weight ('weight' | None, optional): If "weight" use physical, geodesic length to determine segment length. If None use number of nodes (faster).
  • return_lengths (bool): If True, also return lengths of segments according to weight.
Returns
  • segments (list): Segments as list of lists containing node IDs. List is sorted by segment lengths.
  • lengths (list): Length for each segment according to weight. Only provided if return_lengths is True.
def save_swc(self, filepath): 
238    def save_swc(self, filepath):
239        """Save skeleton in SWC format.
240
241        Parameters
242        ----------
243        filepath :      path-like
244                        Filepath to save SWC to.
245
246        """
247        header = dedent(f"""\
248        # SWC format file
249        # based on specifications at http://www.neuronland.org/NLMorphologyConverter/MorphologyFormats/SWC/Spec.html
250        # Created on {datetime.date.today()} using skeletor (https://github.com/navis-org/skeletor)
251        # PointNo Label X Y Z Radius Parent
252        # Labels:
253        # 0 = undefined, 1 = soma, 5 = fork point, 6 = end point
254        """)
255
256        # Make copy of SWC table
257        swc = self.swc.copy()
258
259        # Set all labels to undefined
260        swc['label'] = 0
261        swc.loc[~swc.node_id.isin(swc.parent_id.values), 'label'] = 6
262        n_childs = swc.groupby('parent_id').size()
263        bp = n_childs[n_childs > 1].index.values
264        swc.loc[swc.node_id.isin(bp), 'label'] = 5
265
266        # Add radius if missing
267        if 'radius' not in swc.columns:
268            swc['radius'] = 0
269
270        # Get things in order
271        swc = swc[['node_id', 'label', 'x', 'y', 'z', 'radius', 'parent_id']]
272
273        # Adjust column titles
274        swc.columns = ['PointNo', 'Label', 'X', 'Y', 'Z', 'Radius', 'Parent']
275
276        with open(filepath, 'w') as file:
277            # Write header
278            file.write(header)
279
280            # Write data
281            writer = csv.writer(file, delimiter=' ')
282            writer.writerows(swc.astype(str).values)

Save skeleton in SWC format.

Parameters
  • filepath (path-like): Filepath to save SWC to.
def scene(self, mesh=False, **kwargs): 
284    def scene(self, mesh=False, **kwargs):
285        """Return a Scene object containing the skeleton.
286
287        Returns
288        -------
289        scene :     trimesh.scene.scene.Scene
290                    Contains the skeleton and optionally the mesh.
291
292        """
293        if mesh:
294            if isinstance(self.mesh, type(None)):
295                raise ValueError('Skeleton has no mesh.')
296
297            self.mesh.visual.face_colors = [100, 100, 100, 100]
298
299            # Note the copy(): without it the transform in show() changes
300            # the original meshes
301            sc = tm.Scene([self.mesh.copy(), self.skeleton.copy()], **kwargs)
302        else:
303            sc = tm.Scene(self.skeleton.copy(), **kwargs)
304
305        return sc

Return a Scene object containing the skeleton.

Returns
  • scene (trimesh.scene.scene.Scene): Contains the skeleton and optionally the mesh.
def show(self, mesh=False, **kwargs): 
307    def show(self, mesh=False, **kwargs):
308        """Render the skeleton in an opengl window. Requires pyglet.
309
310        Parameters
311        ----------
312        mesh :      bool
313                    If True, will render transparent mesh on top of the
314                    skeleton.
315
316        Returns
317        --------
318        scene :     trimesh.scene.Scene
319                    Scene with skeleton in it.
320
321        """
322        scene = self.scene(mesh=mesh)
323
324        # I encountered some issues if object space is big and the easiest
325        # way to work around this is to apply a transform such that the
326        # coordinates have -5 to +5 bounds
327        fac = 5 / np.fabs(self.skeleton.bounds).max()
328        scene.apply_transform(np.diag([fac, fac, fac, 1]))
329
330        return scene.show(**kwargs)

Render the skeleton in an opengl window. Requires pyglet.

Parameters
  • mesh (bool): If True, will render transparent mesh on top of the skeleton.
Returns
  • scene (trimesh.scene.Scene): Scene with skeleton in it.
def example_mesh(): 
43def example_mesh():
44    """Load and return example mesh.
45
46    The example mesh is a fruit fly neuron (an olfactory projection neuron of
47    the DA1 glomerulus) segmented from an EM image data set. It is part of the
48    Janelia hemibrain data set (see [here](https://neuprint.janelia.org)) [1].
49
50    References
51    ----------
52    [1] Louis K. Scheffer et al., eLife. 2020. doi: 10.7554/eLife.57443
53    A connectome and analysis of the adult Drosophila central brain
54
55    Returns
56    -------
57    trimesh.Trimesh
58
59    Examples
60    --------
61    >>> import skeletor as sk
62    >>> # Load this example mesh
63    >>> mesh = sk.example_mesh()
64
65    """
66    return tm.load_mesh(obj_path)

Load and return example mesh.

The example mesh is a fruit fly neuron (an olfactory projection neuron of the DA1 glomerulus) segmented from an EM image data set. It is part of the Janelia hemibrain data set (see here) [1].

References

[1] Louis K. Scheffer et al., eLife. 2020. doi: 10.7554/eLife.57443 A connectome and analysis of the adult Drosophila central brain

Returns
  • trimesh.Trimesh
Examples
>>> import skeletor as sk
>>> # Load this example mesh
>>> mesh = sk.example_mesh()