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:
- Typically smaller (less vertices) than the mesh
- 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:
- Some pre-processing of the mesh (e.g. fixing some potential errors like degenerate faces, unreferenced vertices, etc.)
- The skeletonization itself
- 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:
skeletor.pre.fix_mesh()
to fix the meshskeletor.pre.simplify()
to simplify the meshskeletor.pre.contract()
to contract the mesh [1]skeletor.skeletonize.by_vertex_clusters()
to generate a skeletonskeletor.post.clean_up()
to clean up some potential issues with the skeletonskeletor.post.smooth()
to smooth the skeletonskeletor.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)
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
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']
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.
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.
88 @property 89 def vertices(self): 90 """Return skeleton vertices (nodes).""" 91 return self.swc[['x', 'y', 'z']].values
Return skeleton vertices (nodes).
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.
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.
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
.
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).
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.
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.
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
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. IfNone
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 ifreturn_lengths
is True.
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.
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.
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.
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()