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)
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:
skeletor.pre.fix_mesh()
to fix the meshskeletor.pre.simplify()
to simplify the meshskeletor.pre.contract()
to contract the mesh [1]skeletor.skeletonize.vertex_clusters()
to generate a skeletonskeletor.post.clean_up()
to clean up some potential issues with 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
- if the mesh consists of multiple disconnected pieces the skeleton will likewise be fragmented (i.e. will have multiple roots)
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) 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']
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.
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_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
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.
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.
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.
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()