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)

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

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.vertex_clusters() to generate a skeleton
  5. skeletor.post.clean_up() to clean up some potential issues with the skeleton
  6. 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
  • if the mesh consists of multiple disconnected pieces the skeleton will likewise be fragmented (i.e. will have multiple roots)

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)
 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
 88------
 89
 90See docstrings of the respective functions for details.
 91
 92A pipeline might look like this:
 93
 94 1. `skeletor.pre.fix_mesh()` to fix the mesh
 95 2. `skeletor.pre.simplify()` to simplify the mesh
 96 3. `skeletor.pre.contract()` to contract the mesh [1]
 97 4. `skeletor.skeletonize.vertex_clusters()` to generate a skeleton
 98 5. `skeletor.post.clean_up()` to clean up some potential issues with the skeleton
 99 6. `skeletor.post.radii()` to extract radii either by k-nearest neighbours or ray-casting
100
101In my experience there is no one-size-fits-all. You will have to play around to
102find the right approach and parameters to get nice skeletons for your meshes.
103If you need help just open an [issue](https://github.com/navis-org/skeletor/issues).
104
105Also check out the Gotchas below!
106
107# Examples
108
109First load the example mesh (a fruit fly neuron):
110
111```Python
112>>> import skeletor as sk
113>>> mesh = sk.example_mesh()
114>>> mesh
115<trimesh.Trimesh(vertices.shape=(6582, 3), faces.shape=(13772, 3))>
116```
117
118Next see if there is stuff to fix in the mesh (degenerate faces, duplicate
119vertices, etc.):
120
121```Python
122>>> fixed = sk.pre.fix_mesh(mesh, remove_disconnected=5, inplace=False)
123>>> fixed
124<trimesh.Trimesh(vertices.shape=(6213, 3), faces.shape=(12805, 3))>
125```
126
127Now for tubular meshes like this neuron, the "wave front" skeletonization method
128performs really well: it works by casting waves across the mesh and collapsing
129the resulting rings into a skeleton (kinda like when you throw a stone in a
130pond and track the expanding ripples).
131
132```Python
133>>> skel = sk.skeletonize.by_wavefront(fixed, waves=1, step_size=1)
134>>> skel
135<Skeleton(vertices=(1258, 3), edges=(1194, 2), method=wavefront)>
136```
137
138All skeletonization methods return a `Skeleton` object. These are just
139convenient objects to bundle the various outputs of the skeletonization.
140
141```Python
142>>> # x/y/z location of skeleton vertices (nodes)
143>>> skel.vertices
144array([[16744, 36720, 26407],
145       ...,
146       [22076, 23217, 24472]])
147>>> # child -> parent edges
148>>> skel.edges
149array([[  64,   31],
150       ...,
151       [1257, 1252]])
152>>> # Mapping for mesh to skeleton vertex indices
153>>> skel.mesh_map
154array([ 157,  158, 1062, ...,  525,  474,  547])
155>>> # SWC table
156>>> skel.swc.head()
157   node_id  parent_id             x             y             z    radius
1580        0         -1  16744.005859  36720.058594  26407.902344  0.000000
1591        1         -1   5602.751953  22266.756510  15799.991211  7.542587
1602        2         -1  16442.666667  14999.978516  10887.916016  5.333333
161```
162
163SWC is a commonly used format for saving skeletons. `Skeleton` objects
164have a method for quickly saving a correctly formatted SWC file:
165
166```Python
167>>> skel.save_swc('~/Documents/my_skeleton.swc')
168```
169
170If you installed `pyglet` (see above) you can also use `trimesh`'s plotting
171capabilities to inspect the results:
172
173```Python
174>>> skel.show(mesh=True)
175```
176
177<img src="https://github.com/navis-org/skeletor/raw/master/_static/example1.png" alt="skeletor_example" width="100%"/>
178
179That looks pretty good already but let's run some pro-forma postprocessing.
180
181```Python
182>>> sk.post.clean_up(skel, inplace=True)
183<Skeleton(vertices=(1071, 3), edges=(1070, 2))>
184```
185
186So that would be a full pipeline mesh to skeleton. Don't expect your own meshes
187to produce such nice results off the bat though. Chances are you will need to
188play around to find the right recipe. If you don't know where to start, I suggest
189you try out mesh contraction + vertex clustering first:
190
191```Python
192>>> import skeletor as sk
193>>> # Load the example mesh that ships with skeletor
194>>> mesh = sk.example_mesh()
195>>> # Alternatively use trimesh to load/construct your own mesh:
196>>> # import trimesh as tm
197>>> # mesh = tm.Trimesh(vertices, faces) 
198>>> # mesh = tm.load_mesh('some/mesh.obj')
199>>> # Run some general clean-up (see docstring for details)
200>>> fixed = sk.pre.fix_mesh(mesh, remove_disconnected=5, inplace=False)
201>>> # Contract mesh to 10% (0.1) of original volume
202>>> cont = sk.pre.contract(fixed, epsilon=0.1)
203>>> # Skeletonize
204>>> skel = sk.skeletonize.by_vertex_clusters(cont, sampling_dist=100)
205>>> # Replace contracted mesh with original for postprocessing and plotting
206>>> skel.mesh = fixed
207>>> # Add radii (vertex cluster method does not do that automatically)
208>>> sk.post.radii(skel, method='knn')
209>>> skel.show(mesh=True)
210```
211
212# Gotchas
213
214- while this is a general purpose library, my personal focus is on neurons and
215  this has certainly influenced things like default parameter values and certain
216  post-processing steps
217- meshes need to be triangular (we are using `trimesh`)
218- use `sk.pre.simplify` if your mesh is very complex (half a million vertices is
219  where things start getting sluggish)
220- a good mesh contraction is often half the battle
221- if the mesh consists of multiple disconnected pieces the skeleton will
222  likewise be fragmented (i.e. will have multiple roots)
223
224# Benchmarks
225
226<img src="https://github.com/navis-org/skeletor/raw/master/benchmarks/benchmark_2.png" alt="skeletor_benchmark" width="100%"/>
227
228[Benchmarks](https://github.com/navis-org/skeletor/blob/master/benchmarks/skeletor_benchmark.ipynb)
229were run on a 2018 MacBook Pro (2.2 GHz Core i7, 32Gb memory) with optional
230`fastremap` dependency installed. Note some of these functions (e.g.
231contraction and TEASAR/vertex cluster skeletonization) vary a lot in
232speed based on parameterization.
233
234# What about algorithm `X`?
235
236`skeletor` contains some algorithms that I found easy enough to implement
237and useful for my work with neurons. If you have some interesting paper/approach
238that could make a nice addition to `skeletor`, please get in touch on Github.
239Pull requests are always welcome!
240
241# References
242
243`[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.`
244
245The abstract and the paper can be found [here](http://visgraph.cse.ust.hk/projects/skeleton/).
246Also see [this](https://www.youtube.com/watch?v=-H7n59YQCRM&feature=youtu.be) YouTube video.
247
248Some of the code in skeletor was modified from the
249[Py_BL_MeshSkeletonization](https://github.com/aalavandhaann/Py_BL_MeshSkeletonization)
250addon created by #0K Srinivasan Ramachandran and published under GPL3.
251
252# Top-level functions and classes
253At top-level we only expose `example_mesh()` and the `Skeleton` class (which
254you probably won't ever need to touch manually). Everything else is neatly
255tucked away into submodules (see side-bar or above table).
256
257"""
258
259__version__ = "1.2.3"
260__version_vector__ = (1, 2, 3)
261
262from . import skeletonize
263from . import pre
264from . import post
265
266from .skeletonize.base import Skeleton
267from .data import example_mesh
268
269__docformat__ = "numpy"
270
271__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_weighted_edges_from(zip(nodes.node_id.values, nodes.parent_id.values, dists))
170
171        return G
172
173    def save_swc(self, filepath):
174        """Save skeleton in SWC format.
175
176        Parameters
177        ----------
178        filepath :      path-like
179                        Filepath to save SWC to.
180
181        """
182        header = dedent(f"""\
183        # SWC format file
184        # based on specifications at http://www.neuronland.org/NLMorphologyConverter/MorphologyFormats/SWC/Spec.html
185        # Created on {datetime.date.today()} using skeletor (https://github.com/navis-org/skeletor)
186        # PointNo Label X Y Z Radius Parent
187        # Labels:
188        # 0 = undefined, 1 = soma, 5 = fork point, 6 = end point
189        """)
190
191        # Make copy of SWC table
192        swc = self.swc.copy()
193
194        # Set all labels to undefined
195        swc['label'] = 0
196        swc.loc[~swc.node_id.isin(swc.parent_id.values), 'label'] = 6
197        n_childs = swc.groupby('parent_id').size()
198        bp = n_childs[n_childs > 1].index.values
199        swc.loc[swc.node_id.isin(bp), 'label'] = 5
200
201        # Add radius if missing
202        if 'radius' not in swc.columns:
203            swc['radius'] = 0
204
205        # Get things in order
206        swc = swc[['node_id', 'label', 'x', 'y', 'z', 'radius', 'parent_id']]
207
208        # Adjust column titles
209        swc.columns = ['PointNo', 'Label', 'X', 'Y', 'Z', 'Radius', 'Parent']
210
211        with open(filepath, 'w') as file:
212            # Write header
213            file.write(header)
214
215            # Write data
216            writer = csv.writer(file, delimiter=' ')
217            writer.writerows(swc.astype(str).values)
218
219    def scene(self, mesh=False, **kwargs):
220        """Return a Scene object containing the skeleton.
221
222        Returns
223        -------
224        scene :     trimesh.scene.scene.Scene
225                    Contains the skeleton and optionally the mesh.
226
227        """
228        if mesh:
229            if isinstance(self.mesh, type(None)):
230                raise ValueError('Skeleton has no mesh.')
231
232            self.mesh.visual.face_colors = [100, 100, 100, 100]
233
234            # Note the copy(): without it the transform in show() changes
235            # the original meshes
236            sc = tm.Scene([self.mesh.copy(), self.skeleton.copy()], **kwargs)
237        else:
238            sc = tm.Scene(self.skeleton.copy(), **kwargs)
239
240        return sc
241
242    def show(self, mesh=False, **kwargs):
243        """Render the skeleton in an opengl window. Requires pyglet.
244
245        Parameters
246        ----------
247        mesh :      bool
248                    If True, will render transparent mesh on top of the
249                    skeleton.
250
251        Returns
252        --------
253        scene :     trimesh.scene.Scene
254                    Scene with skeleton in it.
255
256        """
257        scene = self.scene(mesh=mesh)
258
259        # I encountered some issues if object space is big and the easiest
260        # way to work around this is to apply a transform such that the
261        # coordinates have -5 to +5 bounds
262        fac = 5 / np.fabs(self.skeleton.bounds).max()
263        scene.apply_transform(np.diag([fac, fac, fac, 1]))
264
265        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
edges

Return skeleton edges.

vertices

Return skeleton vertices (nodes).

radius

Return radii.

skeleton

Skeleton as trimesh Path3D.

skel_map

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

leafs

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_weighted_edges_from(zip(nodes.node_id.values, nodes.parent_id.values, dists))
170
171        return G

Generate networkX representation of the skeletons.

Distance between nodes will be used as edge weights.

Returns
  • networkx.DiGraph
def save_swc(self, filepath):
173    def save_swc(self, filepath):
174        """Save skeleton in SWC format.
175
176        Parameters
177        ----------
178        filepath :      path-like
179                        Filepath to save SWC to.
180
181        """
182        header = dedent(f"""\
183        # SWC format file
184        # based on specifications at http://www.neuronland.org/NLMorphologyConverter/MorphologyFormats/SWC/Spec.html
185        # Created on {datetime.date.today()} using skeletor (https://github.com/navis-org/skeletor)
186        # PointNo Label X Y Z Radius Parent
187        # Labels:
188        # 0 = undefined, 1 = soma, 5 = fork point, 6 = end point
189        """)
190
191        # Make copy of SWC table
192        swc = self.swc.copy()
193
194        # Set all labels to undefined
195        swc['label'] = 0
196        swc.loc[~swc.node_id.isin(swc.parent_id.values), 'label'] = 6
197        n_childs = swc.groupby('parent_id').size()
198        bp = n_childs[n_childs > 1].index.values
199        swc.loc[swc.node_id.isin(bp), 'label'] = 5
200
201        # Add radius if missing
202        if 'radius' not in swc.columns:
203            swc['radius'] = 0
204
205        # Get things in order
206        swc = swc[['node_id', 'label', 'x', 'y', 'z', 'radius', 'parent_id']]
207
208        # Adjust column titles
209        swc.columns = ['PointNo', 'Label', 'X', 'Y', 'Z', 'Radius', 'Parent']
210
211        with open(filepath, 'w') as file:
212            # Write header
213            file.write(header)
214
215            # Write data
216            writer = csv.writer(file, delimiter=' ')
217            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):
219    def scene(self, mesh=False, **kwargs):
220        """Return a Scene object containing the skeleton.
221
222        Returns
223        -------
224        scene :     trimesh.scene.scene.Scene
225                    Contains the skeleton and optionally the mesh.
226
227        """
228        if mesh:
229            if isinstance(self.mesh, type(None)):
230                raise ValueError('Skeleton has no mesh.')
231
232            self.mesh.visual.face_colors = [100, 100, 100, 100]
233
234            # Note the copy(): without it the transform in show() changes
235            # the original meshes
236            sc = tm.Scene([self.mesh.copy(), self.skeleton.copy()], **kwargs)
237        else:
238            sc = tm.Scene(self.skeleton.copy(), **kwargs)
239
240        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):
242    def show(self, mesh=False, **kwargs):
243        """Render the skeleton in an opengl window. Requires pyglet.
244
245        Parameters
246        ----------
247        mesh :      bool
248                    If True, will render transparent mesh on top of the
249                    skeleton.
250
251        Returns
252        --------
253        scene :     trimesh.scene.Scene
254                    Scene with skeleton in it.
255
256        """
257        scene = self.scene(mesh=mesh)
258
259        # I encountered some issues if object space is big and the easiest
260        # way to work around this is to apply a transform such that the
261        # coordinates have -5 to +5 bounds
262        fac = 5 / np.fabs(self.skeleton.bounds).max()
263        scene.apply_transform(np.diag([fac, fac, fac, 1]))
264
265        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()