How to save the skeleton without the mesh?

[Ameliecc]

Hi~ I would like to save the skeleton without the mesh, however I don't know how to do it. Would you please give me some advice? Following are the code I use. Thank you for your kind help.

import open3d as o3d
import matplotlib.pyplot as plt
import numpy as np
import skeletor as sk
import trimesh
import pandas as pd

# reading a mesh to trimesh
mesh_name = '400.stl'
mesh = trimesh.load(mesh_name)
mesh_o3d = o3d.io.read_triangle_mesh(mesh_name)

fixed = sk.pre.fix_mesh(mesh, remove_disconnected=5, inplace=False)

# Contract the mesh
cont = sk.pre.contract(fixed, iter_lim=15)  

# Extract the skeleton from the contracted mesh
skel = sk.skeletonize.by_vertex_clusters(cont, sampling_dist=1)  

# # Clean up the skeleton
skel_clean = sk.post.clean_up(skel, mesh) 
# # Add/update radii
# # swc['radius'] = sk.post.radii(swc, mesh, method='knn', n=5, aggregate='mean')
sk.post.radii(skel_clean, method='knn') 

skel.show(mesh=False)
skel.swc.head()
skel.save_swc('save_swc')
print(skel.swc.head())

[schlegelp]

Hi. Looks like you asked almost the same question in #4.

You seem to have already found the .save_swc() method. SWC is a commonly used format to store skeletons (see format specs here) and the only one currently implemented.

They way to use it would be to run this:

skel.save_swc('skeleton.swc')

This will create a skeleton.swc file in the current directory. These files are effectively just comma-separated tables, i.e. you can open them in a text editor if you like.

Whether the SWC format is helpful in your case really depends on what you want to do with the file. I might have other/better suggestions if you tell me a bit more e.g. about what other software you want to load it into.

[Ameliecc]

Hi, thank you very much for your help, and I'm really sorry for taking so long to reply to your message. My goal is to save the calculated skeleton lines. As you can see the green point (skeleton) in the middle of the figure. And I wonder if there is a solution.

400.zip

[schlegelp]

No problem. The question is not whether there is a way to save the skeleton - there are plenty. The question is rather what you want to do with them downstream because that will determine the format.

[Ameliecc]

I would like to optimize the skeleton line and calculate the length of the extracted skeleton line. However, when i use " skel.save_swc('skeleton.swc')", it saved the mesh only. As depicted in the picture. I looked for instructions, but I couldn't find a way to save the skel alone. I tried using two point clouds to differentiate, but found that kel.save_swc('save_swc.txt') is different from the point cloud in the last saved pcd file. Therefore, I would like to ask how to save the calculated skeleton alone. Thanks a lot.

[Ameliecc]

Attached is the code.

import open3d as o3d
import matplotlib.pyplot as plt
import numpy as np
import skeletor as sk
import trimesh
import pandas as pd

# reading a mesh to trimesh
mesh_name = '400.stl'
mesh = trimesh.load(mesh_name)
mesh_o3d = o3d.io.read_triangle_mesh(mesh_name)

fixed = sk.pre.fix_mesh(mesh, remove_disconnected=5, inplace=False)

# Contract the mesh
cont = sk.pre.contract(fixed, iter_lim=20)   

# Extract the skeleton from the contracted mesh
skel = sk.skeletonize.by_vertex_clusters(cont, sampling_dist=1)    # 生成骨架

# # Clean up the skeleton
skel_clean = sk.post.clean_up(skel, mesh)   
# # Add/update radii
# # swc['radius'] = sk.post.radii(swc, mesh, method='knn', n=5, aggregate='mean')
sk.post.radii(skel_clean, method='knn', n=5, aggregate='mean') 

skel.show(mesh=False)

skel.swc.head()
skel.save_swc('save_swc.txt')
print(skel.swc.head())

# # visualizing scattered points
max_val = max(max(skel_clean.swc.x), max(skel_clean.swc.y), max(skel_clean.swc.z))
fig = plt.figure()
ax = fig.gca(projection='3d')
plt.axis('off')
plt.xlim([-max_val, max_val])
plt.ylim([-max_val, max_val])
ax.set_zlim(-max_val, max_val)
ax.scatter(skel_clean.swc.x, skel_clean.swc.y, skel_clean.swc.z, s=10)
plt.show()
# #
# # # visualizing using open3d
xyz = np.zeros((skel_clean.swc.x.size, 3))
xyz[:, 0] = skel_clean.swc.x
xyz[:, 1] = skel_clean.swc.y
xyz[:, 2] = skel_clean.swc.z
pcd = o3d.geometry.PointCloud()
pcd.points = o3d.utility.Vector3dVector(xyz)
mesh_o3d.compute_vertex_normals()
o3d.visualization.draw_geometries([mesh_o3d, pcd])
o3d.visualization.draw_geometries([pcd])
o3d.io.write_point_cloud("400_output.pcd", pcd)
# #
# # # exporting csv file
export_file_path = mesh_name[:-4] + '.csv'
print(export_file_path)
# skel.to_csv(export_file_path, index=False, header=True)

# skel.show(mesh=True)

[schlegelp]

I would like to optimize the skeleton line and calculate the length of the extracted skeleton line.

That's fair enough. I might have recommended you use navis for this. It's principally made for analysing neurons but since neurons are also just directed acyclic graphs it works just as well. Importantly, it plays nicely with skeletor (since I wrote skeletor for navis) and lets you process your skeleton:

>>> import skeletor as sk
>>> mesh = sk.example_mesh()
>>> fixed = sk.pre.fix_mesh(mesh, remove_disconnected=5, inplace=False)
>>> skel = sk.skeletonize.by_wavefront(fixed, waves=1, step_size=1)

 # Turn skeleton into a navis Neuron

>>> import navis
>>> n = navis.TreeNeuron(skel) 

# Some examples of what you can do with it now:

>>> n.cable_length
193260.085048636

>>> s = navis.smooth_skeleton(n, window=5)  # Smooth
>>> s.cable_length
209428.27956084194

>>> pr = navis.prune_twigs(n, size=1000)  # prune small twigs
>>> pr.cable_length
120032.67049704675

>>> ln = navis.longest_neurite(n, from_root=False)  # Extract the longest continuous segment
>>> ln.cable_length
55611.49568331024

However, when i use " skel.save_swc('skeleton.swc')", it saved the mesh only. As depicted in the picture.

I'm not sure what I'm seeing in the plot you've attached but it doesn't seem to match the green points in the earlier rendering. The bottom line is this: .save_swc definitely doesn't save the mesh, it saves the skeleton. If the points in the SWC file don't match your expectations then something went wrong with the skeletonization.

If you want to get your skeleton into open3d you could do this:

>>> import open3d as o3d
>>> ls = o3d.geometry.LineSet(o3d.utility.Vector3dVector(skel.vertices),
...                           o3d.utility.Vector2iVector(skel.edges))
>>> open3d.visualization.draw_geometries([ls])

[Ameliecc]

Thank you for your help. ❤