How to prune the skeleton

[mikami520]

Hi, I am working on the skeleton of a linear structure. Since it has a lot of small branches, I want to prune them and select the largest one. Do you have any recommendation or give me some help about how to do it? The vtk file is attached and you should convert to Trimesh object firstly. By the way, is there a way to output skeleton as a vtk or stl or 2D imagefile? Since I tried to use trimesh.Path3D.exportto output a stl file, but when I open the ParaView to load the path, I cannot see it. Do you have any idea about this issue? I am looking forward to your response, thanks!

pred_ETV2_178.vtk.zip

[schlegelp]

OK, first re the export: Trimesh appears to only support dxf, dict and svg as output formats for Path3D. See file_type in the docstring:

>>> help(trimesh.Path3D.export)
Signature: s.skeleton.export(file_obj=None, file_type=None, **kwargs)
Docstring:
Export the path to a file object or return data.

Parameters
---------------
file_obj : None, str, or file object
  File object or string to export to
file_type : None or str
  Type of file: dxf, dict, svg

Returns
---------------
exported : bytes or str
  Exported as specified type

I work mostly with neurons (see below) and in that world, SWC is the most common format for skeletons. I'm not familiar with the vtk format but I would perhaps look into meshio and check if it's capable of dealing with meshes that have only vertices and edges (but no faces).

Now for the post-processing: although this might change in the future, skeletor's capabilities in this department are fairly limited. If I was you, I'd use navis (which is what I wrote skeletor for). navis is nominally meant for skeletons of neurons but at the end of the day it doesn't matter what the skeleton represents. Here are some quick pointers from your example #27:

import navis
# Create a neuron from your skeleton
n = navis.TreeNeuron(skel, soma=None)
# Example 1: Remove all twigs under a certain size (adjust `size` parameter!)
no_twigs = navis.prune_twigs(n, size=100)
# Example 2: keep only the longest linear section in your skeleton
long1 = navis.longest_neurite(n, from_root=False)
# Example 3: keep only the two longest linear section sin your skeleton
long2 = navis.longest_neurite(n, n=2, from_root=False)

# Turn any of the above back into a skeletor.Skeleton
# (note that processing might have broken the mesh vertices <->  skeleton nodes mapping
skel2 = sk.Skeleton(no_twigs.nodes[['node_id', 'parent_id', 'x', 'y', 'z', 'radius']])

[mikami520]

OK, first re the export: Trimesh appears to only support dxf, dict and svg as output formats for Path3D. See file_type in the docstring:

>>> help(trimesh.Path3D.export)
Signature: s.skeleton.export(file_obj=None, file_type=None, **kwargs)
Docstring:
Export the path to a file object or return data.

Parameters
---------------
file_obj : None, str, or file object
  File object or string to export to
file_type : None or str
  Type of file: dxf, dict, svg

Returns
---------------
exported : bytes or str
  Exported as specified type

I work mostly with neurons (see below) and in that world, SWC is the most common format for skeletons. I'm not familiar with the vtk format but I would perhaps look into meshio and check if it's capable of dealing with meshes that have only vertices and edges (but no faces).

Now for the post-processing: although this might change in the future, skeletor's capabilities in this department are fairly limited. If I was you, I'd use navis (which is what I wrote skeletor for). navis is nominally meant for skeletons of neurons but at the end of the day it doesn't matter what the skeleton represents. Here are some quick pointers from your example #27:

import navis
# Create a neuron from your skeleton
n = navis.TreeNeuron(skel, soma=None)
# Example 1: Remove all twigs under a certain size (adjust `size` parameter!)
no_twigs = navis.prune_twigs(n, size=100)
# Example 2: keep only the longest linear section in your skeleton
long1 = navis.longest_neurite(n, from_root=False)
# Example 3: keep only the two longest linear section sin your skeleton
long2 = navis.longest_neurite(n, n=2, from_root=False)

# Turn any of the above back into a skeletor.Skeleton
# (note that processing might have broken the mesh vertices <->  skeleton nodes mapping
skel2 = sk.Skeleton(no_twigs.nodes[['node_id', 'parent_id', 'x', 'y', 'z', 'radius']])

@schlegelp I tried to convert long1 to the Skeleton object and show the pruning result. However, it said that the index out of the bounds. Do you know the reasons?

[schlegelp]

Ah, I had forgotten that skeletor expects sequentially numbered node IDs while navis doesn't care much about that. During pruning we dropped a bunch of nodes and navis doesn't renumber them.

We have to add an intermediate step at which we convert the more forgiving navis node table into a SWC table that skeletor is happy with;

# This renumbers nodes
swc = navis.io.swc_io.make_swc_table(long1)
# We also need to rename some columns
swc = swc.rename({'PointNo': 'node_id', 'Parent': 'parent_id', 'X': 'x', 'Y': 'y', 'Z': 'z', 'Radius': 'radius'}, axis=1).drop('Label', axis=1)
# Skeletor excepts node IDs to start with 0, but navis starts at 1 for SWC 
swc['node_id'] -= 1
swc.loc[swc.parent_id > 0, 'parent_id'] -= 1 
# Create the skeletor.Skeleton
skel2 = sk.Skeleton(swc)

For what it's worth: navis has an independent built-in plotting system too. So if all you're after is visualization, you don't need to go back to skeletor.

# Plot only the skeleton
navis.plot3d(long1)

# From terminal only: to clear the 3d plot
navis.clear3d()

# From terminal only: to close the 3d plot
navis.close3d()

# Plot skeleton with mesh 
navis.plot3d([long1,  skel.mesh])

[mikami520]

Ah, I had forgotten that skeletor expects sequentially numbered node IDs while navis doesn't care much about that. During pruning we dropped a bunch of nodes and navis doesn't renumber them.

We have to add an intermediate step at which we convert the more forgiving navis node table into a SWC table that skeletor is happy with;

# This renumbers nodes
swc = navis.io.swc_io.make_swc_table(long1)
# We also need to rename some columns
swc = swc.rename({'PointNo': 'node_id', 'Parent': 'parent_id', 'X': 'x', 'Y': 'y', 'Z': 'z', 'Radius': 'radius'}, axis=1).drop('Label', axis=1)
# Skeletor excepts node IDs to start with 0, but navis starts at 1 for SWC 
swc['node_id'] -= 1
swc.loc[swc.parent_id > 0, 'parent_id'] -= 1 
# Create the skeletor.Skeleton
skel2 = sk.Skeleton(swc)

For what it's worth: navis has an independent built-in plotting system too. So if all you're after is visualization, you don't need to go back to skeletor.

# Plot only the skeleton
navis.plot3d(long1)

# From terminal only: to clear the 3d plot
navis.clear3d()

# From terminal only: to close the 3d plot
navis.close3d()

# Plot skeleton with mesh 
navis.plot3d([long1,  skel.mesh])

Thanks for your help. The following part for me is to compute the closest distance between pruning skeleton and original mesh so that I can construct a new mapping from new skeleton to original mesh, since you mentioned before that mapping from skeleton to mesh will be affected by using navis