navis-org/skeletor

I've created a simple code for using Skeletor without using cloudvolume

TzabarDolev opened this issue · 4 comments

Hey,
I was unable to work with your usage example of the code. Cloudvolume is not a particularly user-friendly library and I wanted a local solution using a stl file. I used trimesh instead of cloudvolume and loaded a file and worked on it directly. I add a short code I built for a skeletal visualization using trimesh and a stl file.

Thank you so much for implementing this article, it has helped me a lot.
Tzabar

https://github.com/TzabarDolev/Mesh_skeltonization/blob/main/skeletor_test.py

Hi @TzabarDolev. Thanks for example! I should probably put something more basic in the README. In the meantime, I'll leave this issue open.

Hi, how to save the skeleton without the mesh?

I run the following code, however, the mesh seems strange, could you tell me what the problem is?
Thanks a lot.
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=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')  # 通过 k 最近邻或光线投射提取半径

skel.show()

skel.swc.head()
skel.save_swc('save_swc')
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])
# #
# # # 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)

After "skel.show()"
image

The mesh:
image

Hi. I took the liberty to format your code (see this link) - makes reading it easier.

I see you've found the save to SWC method - that's a commonly used format for skeletons.

As to why your mesh isn't skeletonizing correctly: it's really hard to tell without having the mesh in front of me but at first glance I suspect that the contraction didn't work (well) and that perhaps your sampling_dist for the skeletonization step is too low.

First, the mesh contraction is unfortunately finicky and really depends on your mesh. So the first thing I would do is to inspect the contracted mesh (cont in your code) and see if it indeed is somewhat contracted.

Second, if the mesh contraction looks useful I would try to play with the sampling_dist parameter - this depends on the scale of your (contracted) mesh. To get a feel for where to start it might be worth checking the minimum/mean edge lengths:

>>> # From the example mesh
>>> cont.edges_unique_length.mean()
95.78825249307354
>>> cont.edges_unique_length.min()
1.0634794638400474e-05

For the example mesh the contraction works fairy well which is why some the edges end up being effectively zero. So here, I'd start with sampling distance of around 1 and see where it gets me.

Lastly, if all things fail I would perhaps try the wavefront skeletonization from the example. It's more robust and if it still produces garbage there might be something wrong with your mesh (looks like there are some funky faces in the screenshot).

Hope this helps. I'm also happy to have a crack at this if you would share your 400.stl here.