A quite complex vascular tree with strange radii

[Evave1204]

Hi, I was trying to skeletonize a vascular tree(.stl) and I tried quite a lot skeleton methods but all got either a strange branch shape or incorrect radius for the thickest branch.

[schlegelp]

Mhmm, those vasculatures should skeletonize pretty well. Would you mind sharing the mesh (you might have to zip it for GitHub to accept the upload) so I can take a closer look?

[Evave1204]

Can I send the data by email?

[schlegelp]

I'd prefer here as my inbox is overflowing as is but email is Ok too.

[Evave1204]

Hi, I was trying to upload the compressed zip file, but the size is still too big which has 109MB. Is the simplified file okay for you to process?

[Evave1204]

Uploading final_model(simple_0.01).stl.zip…

[schlegelp]

Sorry, you seem to have posted a link to this issue rather than a file.

[Evave1204]

Hi, I ve found out why the radii is incorrect(because I reconstructed a multi-nodes tree to re-define the root point, and incorrectly assign the wrong radius value to the corresponding edge), but I still have a question about what's the relationship between skel.radius and skel.edges ?

[schlegelp]

The radii in skel.radius match the skel.vertices, i.e. there is a radius for every node in the skeleton.

If you need a radius for every edge rather than every node, you could just take the mean over the radii of the two nodes making up the edge:

>>> node_radii = skel.radius[0] # see * note below
>>> len(skel.vertices), len(node_radii)
(1103, 1103)
>>> edge_radii = np.mean(node_radii[skel.edges], axis=1)
>>> len(skel.edges), len(edge_radii)
(1102, 1102)

*I just realised that the radius property has a trailing , in the return statement which means it currently returns a tuple of (radius, ). Not sure how I missed this but will be fixed in the next release.

[schlegelp]

I also just added a Skeleton.reroot() method. Not released yet but if you re-install skeletor from Github you can take it for a spin.

[Evave1204]

Thanks you for ur solution!! But I still have some questions below:

1. If I have a tree structure like my data, is there any method which can help me to quickly find out where the root is?
2. Which skeletonize method you recommend to skeletonize this tree structure(I want lines to be as smooth as possible)
3. Which method could help me to get more exact radius?

[schlegelp]

Re 1

Do you mean find out where the root of the skeleton is (either by vertex ID or by xyz location), or the base of your vasculature?

Re 2+3

It's hard to tell without the mesh (looks like you posted it but then removed it again?) or even having seen a screenshot with the radii you got. In any event: I'd try the wave front method.

[Evave1204]

1.The base of vasculature
2/3. Yes I posted it, but it is the my lab data so i dont want to remain it here, I'm posting it again in attachment. If you have downloaded it please tell me.

[schlegelp]

I downloaded it. Already had a play and as expected the wavefront method is doing best:

import trimesh as tm
import skeletor as sk
m = tm.load_mesh('final_model(simple_0.01).stl')
fixed = sk.pre.fix_mesh(m, remove_disconnected=10, inplace=False)
skel = sk.skeletonize.by_wavefront(fixed)

The first screenshot shows the mesh with the skeleton in yellow. It looks OK but there are issues: for example the skeleton seems to have been pulled out of the mesh at the branch points (red circles). That may be fixable by using a higher res mesh but I doubt it will go away entirely.

In addition the radii seem to consistently be too small and the error is bigger in the backbone (see second screenshot).

These are issues I have come across before and the reason seems to be that the wave (by design) "jumps" from vertex to vertex instead of traveling a fix distance. That has two implications:

1. The resolution of the skeleton is bounded by the resolution of the mesh: higher res mesh = higher res skeleton and vice versa
2. The wave can distort and form elongated shapes instead of nice round circles which causes e.g. radii to be off. This seems to be less of an issue with thin, spindly meshes (like neurons) - yours is pretty big, which may be why it's causing more trouble here.

I have worked on a version of the wavefront skeletonization that does travel a fixed distance and is independent of the resolution of the mesh. It's computationally a lot more expensive (still only a few minutes on your mesh) but it does produce nicer results. If you want to try it out you need to install skeletor from the better_wavefront branch:

skel = sk.skeletonize.by_wavefront_exact(fixed, step_size=25)

Overall that new method looks very encouraging. Unfortunately, there are still some issues that causes the topology (i.e. the connections between nodes) to be wrong. I will need to look into it again.

[Evave1204]

Thanks for your suggestions and new methods!!

I have tried the new method which is sk.skeletonize.by_wavefront_exact(fixed, step_size=25), but as your said, the topology of the tree looks not pretty much good, please let me know if there is any method that could optimize the topology.

[Evave1204]

Hi, I noticed u said the resolution does matter, and I have two questions below:

1. I am wondering if the resolution makes great effect on the result of radius? Because the document I sent to you has been simplified to 0.01.
2. I am trying to understand how the wavefront works, and my question is if we set the seed point at the root, will the wavefront work better on returning more exact radius?

[Evave1204]

Hi, I have tried your methods but there are some issues with the tree topology.
My code:

fixed = sk.pre.fix_mesh(simple_mesh, remove_disconnected=10, inplace=False)
skel = sk.skeletonize.by_wavefront_exact(fixed, step_size=25)
skel.show(mesh=True)

and the result is:

which not looks like the skeleton you given.

And I wanna know how did you show the mesh with the ball(presenting the radii).
Thanks for your time!!!

[schlegelp]

Your skeleton looks much worse than what I remember from my first trials. Is it possible that the mesh you're using has lower resolution? Looks pretty coarse in your screenshot.

In any event: the by_wavefront_exact method is still WIP and what's currently on Github clearly isn't ready for primetime yet. I have some uncommitted changes that improve things a lot (faster, better topology) but it's still a bit finicky. I'll be working on this on and off - I can't give you an ETA but will keep you posted on this channel.

And I wanna know how did you show the mesh with the ball(presenting the radii).

I use octarine:

import octarine as oc 

# Open a viewer 
# (if you're in an iPython terminal you may have to start a gui event loop via e.g. `%gui qt6`)
v = oc.Viewer()

# Add the skeleton points
v.add_points(skel.vertices, size=skel.radius, size_space='world')

You can also visualize the skeleton itself but for that you currently have to go through navis and the octarine-navis plugin:

>>> import navis
>>> v.add_neurons(navis.TreeNeuron(skel))

I should probably add a standalone octarine-skeletor plugin.

[Evave1204]

1. I did recheck the mesh what I used is the correct resolution which is 0.01, but the topology still look like the screenshot I posted above.

m = tm.load('/Users/hcj/科研/new_project/Data/Skeleton Data/final model/simple_0.01/final_model(simple_0.01).stl')
fixed = sk.pre.fix_mesh(m, remove_disconnected=10, inplace=False)
skel = sk.skeletonize.by_wavefront_exact(fixed, step_size=25)
skel.show(mesh=True)

2. For the installation of octarine, pip install octarine3d[all] on the documentation should be pip install 'octarine3d[all]', otherwise it does not work. After installing it, there is an error making import octarine as oc do not work:

File ~/Library/Python/3.9/lib/python/site-packages/octarine/plugins.py:11, in register_plugins()
      [9](https://file+.vscode-resource.vscode-cdn.net/Users/hcj/%E7%A7%91%E7%A0%94/new_project/~/Library/Python/3.9/lib/python/site-packages/octarine/plugins.py:9) """Register a plugin."""
     [10](https://file+.vscode-resource.vscode-cdn.net/Users/hcj/%E7%A7%91%E7%A0%94/new_project/~/Library/Python/3.9/lib/python/site-packages/octarine/plugins.py:10) # Find all plugins that defined an entry point
---> [11](https://file+.vscode-resource.vscode-cdn.net/Users/hcj/%E7%A7%91%E7%A0%94/new_project/~/Library/Python/3.9/lib/python/site-packages/octarine/plugins.py:11) discovered_plugins = entry_points(group="octarine.plugins")
     [13](https://file+.vscode-resource.vscode-cdn.net/Users/hcj/%E7%A7%91%E7%A0%94/new_project/~/Library/Python/3.9/lib/python/site-packages/octarine/plugins.py:13) # Go over each of the plugins
     [14](https://file+.vscode-resource.vscode-cdn.net/Users/hcj/%E7%A7%91%E7%A0%94/new_project/~/Library/Python/3.9/lib/python/site-packages/octarine/plugins.py:14) for plugin in discovered_plugins:
     [15](https://file+.vscode-resource.vscode-cdn.net/Users/hcj/%E7%A7%91%E7%A0%94/new_project/~/Library/Python/3.9/lib/python/site-packages/octarine/plugins.py:15)     # Import the module

TypeError: entry_points() got an unexpected keyword argument 'group'

which might be caused by the old version of importlib.metadata

[schlegelp]

1. The topology issues currently in by_wavefront_exact have nothing to do with the resolution of the mesh.
2. Thanks for reporting it. I opened importlib.metadata versions causing issues schlegelp/octarine#8 in octarine and it should be fixed with new version 0.2.1

[Evave1204]

I have tried to use octarine to visualize the tree, and I found the vertices of the skeleton by wavefront_exact are absolutely correct, but they might have wrong edges, like you said above.

So I think the reason causing this case might be:
point 1, point 2, point 3 are in the same mesh, and the actual connection relationship is p1 -> p2 -> p3, but if these three points form a Equilateral Triangle, the connection will be p1 -> p2, p1 -> p3 ? (I am not sure)

[Evave1204]

Hi, I have looked into your algorithm and tried to optimize it by comparing the distance of center1 and center2 with the sum of their radius to judge if they can be the relationship of parent of child. But I am confused why the all positions of centers are correct(shown by octarine) but connection is wrong(shown by skeleton.show())

I am grateful if you can give me some tips

[Evave1204]

1. Hi, I got some success on optimizing the skeletonize_exact algorithm, what I did is replace the spanning tree which used to record distance with distances of centers. the code is below:

        this_step_track_centers = centers[is_this_step]
        prev_step_track_centers = centers[is_prev_step]
        this_step_track_verts = np.unique(
            data[is_this_step, 2]
        )  # All track-vertices in this step
        prev_step_track_verts = np.unique(
            data[is_prev_step, 2]
        )  # All track-vertices in the previous step

        # Get distances between current and previous track vertices
        # Note to self: this step takes up about 60% of the time at the moment
        # There should be a way to speed this up.
        # vertices = [10, 11, 20, 21, 30, 31] | this_step_track_verts = [20, 21] | prev_step_track_verts = [10, 11]
        # According to the connection relationship of the tree, get the min distance from 'this_step_track_verts' to # 'prev_step_track_verts' 
        # min_dist(20 -> 10), min_dist(20 -> 11), min_dist(21 -> 10), min_dist(21 -> 11), ...
        # -> d = [[3.0, 2.0],[1.5, 2.5]]
        d = np.linalg.norm(
        this_step_track_centers[:, np.newaxis, :] - prev_step_track_centers[np.newaxis, :, :],
        axis=2
        )
        # Get the closest previous track-vertex for each current track-vertex
        # d.argmin(axis=1) returns the index of min value of each element 
        # d = [[3.0, 2.0],[1.5, 2.5]] | min([3.0, 2.0]) = 2.0 = index(1) | min([1.5, 2.5]) = 1.5 = index(0) | d.argmin(axis=1) = [1, 0]
        # prev_step_track_verts = [10, 11] | d.argmin(axis=1) = [1, 0] | cloest_prev_track_verts        # = prev_step_track_verts[d.argmin(axis=1) = [11, 10]
        closest_prev_track_verts = prev_step_track_verts[d.argmin(axis=1)]
        # 将this_step_track_verts 和 cloest_prev_track_verts中的一一对应
        # this_step_track_verts = [20, 21] | cloest_prev_track_verts = [11, 10] | i = 2 | 
        # is_this_step = [False,False,True,True,False,False] 代表的当前[21, 20]对应的是index 2 3
        for this, prev in zip(this_step_track_verts, closest_prev_track_verts):
            # parents[2] = step_id_map[(1, 11)] = 1 | parents[3] = step_id_map[(1, 10)] = 1
            parents[is_this_step & (data[:, 2] == this)] = step_id_map[(i - 1, prev)]
    return centers, radii, parents

This optimized algorithm could give a better connection between nodes

Could you give me some feedback please because the skeletonized mesh has not been perfect yet.