Segmentation of mesh into smaller parts ?

[GreKro]

Hi everyone !

I'm wondering, is there the possibility to do segmentation of mesh in pyvista or other tools related to mesh processing. As segmentation, I mean dividing one solid to smaller ones.

I have a 3D models of the trees and I want to split its branches into separate files. Someone has some solution or idea which tools or algorithms would be useful to this task? On the below images, I painted what I mean, on the second picture you can see green lines which show an example of segment edges. Additionally, I attached test data to this post.

tree.zip

Best wishes,

[akaszynski]

Hello,

Thanks for the challenge! Turns out there are some features that should be added to pyvista, but for the moment I can get this to work with the following:

"""This example demonstrates how to extract meshes using ghost cells.

It works by:
1.  Displaying the entire mesh as a grid
2.  Having the user select new sub-meshes with a rectangular selector
3.  Removing those selected cells from the original grid.
4.  Storing the new selected meshes
5.  Finally converting those meshes back to PolyData

"""
import numpy as np
import pyvista as pv

mesh = pv.read('./tree.ply')
grid = mesh.cast_to_unstructured_grid()

pl = pv.Plotter()
pl.add_mesh(grid, color='w')

# picking example
picked = []
legend = []
def split_mesh(selected_grid):
    """Adds a new mesh to the plotter each time cells are picked, and
    removes them from the original mesh"""

    # if nothing selected
    if not selected_grid.n_cells:
        return

    # remove the picked cells from main grid
    ghost_cells = np.zeros(grid.n_cells, np.uint8)
    ghost_cells[selected_grid['orig_extract_id']] = 1
    grid.cell_arrays[vtk.vtkDataSetAttributes.GhostArrayName()] = ghost_cells
    grid.RemoveGhostCells()

    # add the selected mesh this to the main plotter
    color = np.random.random(3)
    legend.append(['picked mesh %d' % len(picked), color])
    pl.add_mesh(selected_grid, color=color)
    pl.add_legend(legend)

    # track the picked meshes and label them
    selected_grid['picked_index'] = np.ones(selected_grid.n_points)*len(picked)
    picked.append(selected_grid)

# enable cell picking with our custom callback
pl.enable_cell_picking(mesh=grid, callback=split_mesh, show=False)
pl.show()

Finally, convert these back to meshes and then plot and save them.

# convert these meshes back to surface meshes (PolyData)
meshes = []
for selected_grid in picked:
    meshes.append(selected_grid.extract_surface())

# plot final separated meshes for fun
pv.plot(meshes)

# save these meshes somewhere
for i, mesh in enumerate(meshes):
    mesh.save('mesh_%03d.ply' % i)

[GreKro]

Thanks for answering. I'm impressed with the ability to work interactively in pyvista. I'm looking for a solution that would help me do it without manual intervention. Your solution is great too. Maybe I will use it. Thanks once more time. Best wishes,

[akaszynski]

There is a solution that works without manual intervention, but it will just attempt to divide up the mesh into equal parts. You can also cut by planes, but you'll have to figure out the location of the planes as well.

Please close this issue if you think it's been solved. Otherwise, please let me know if there's something else to be addressed.

[GreKro]

Yep, indeed I need to figure out how to define the location of cutting planes. Thanks a lot for your help. This functionality of pyvista is great.

[banesullivan]

This is a really cool challenge! At the moment, I cannot think of a branching segmentation algorithm that would be able to detect when a branch occurs and label that. However, maybe you could design such an algorithm? Or maybe we could train an ML model to do it using PyVista's interactive selection tools to generate the training data set?

---

In the meantime, it may be possible to do something with the mesh's obbTree:

reference https://github.com/pyvista/pyvista/issues/555 for more info on obbTrees and the fact that we need to improve them

import pyvista as pv

tree = pv.read("tree.ply")

boxes = pv.PolyData()
obb = tree.obbTree
obb.GenerateRepresentation(5, boxes)

p = pv.Plotter(notebook=False)
p.add_mesh(boxes, show_edges=True, opacity=0.25, color="green")
p.add_mesh(tree)
p.show()

[GreKro]

@banesullivan I know that some people do it by using the Graph theory - Minimum Spanning Tree (MST) and Dijkstra algorithms. When you convert mesh into one direction 3D graph where the lowest point of the tree is the base node and you simplified this graph you can detect branch nodes (these nodes have more than two connections (?)). When you have coordinates of these nodes you can use it as cutting planes. I think this process is called skeletonization in your field. I tried to find some Python or R ready-to-use libraries to mesh skeletonization but I can't find it. Another option is to use some algorithms to create directly centerline of 3D models (skeleton) maybe here you know some algorithm?

If you will create these boxes smaller and get their centroids you get nodes to create MST. Is there the option to get box centroid coordinates? how did it work? Can I define the box dimensions? Another option is to build Octree and get boxes centroid coordinates. Is there an option to build an Octree in pyvista ?

I think there are simpler tools than ML to do it. Of course, this is interesting but to be honest, I suppose that will be time-consuming.

[GreKro]

Hi,

@banesullivan @akaszynski
If I have the coordinates of the places where the branches start, how can I divide mesh properly in pyvista ? can you help me with it ? below I write an example of coordinates for 3 points for my test tree dataset.

cuting_plane_coords1 = [34.08, 71.23, 125.62]
cuting_plane_coords2 = [33.07, 72.46, 134.05]
cuting_plane_coords3 = [33.90, 70.43, 136.72]

Thanks in advance for the response. Best wishes

[akaszynski]

Hi @GreKro,

Please look into the clip method for a mesh. For example:

import pyvista as pv

mesh = pv.Sphere()
clipped = mesh.clip(normal=[1, 0, 0], origin=[0, 0, 0])

You'll of course have to compute the origin and normal of your planes, but that's not too difficult since you already have a fully defined plane by three points.

[GreKro]

@akaszynski I tried this method but I don't know how to define properly normal direction when I try with different values the results are rather poor. Any idea? What you mean to compute the origin? my coordinates are not the origin of planes? Thanks in advance for the answer

[banesullivan]

You will need some sort of normal as well (you already have the origin points) and perform a cut by a clipping plane using the clip filter.

---

Also, I stumbled upon the Vascular Modeling Toolkit which has tools for exactly what you are trying to do. Check out these two tutorials

• http://www.vmtk.org/tutorials/Centerlines.html
• http://www.vmtk.org/tutorials/BranchSplitting.html

And some other useful stuff in the VTK ecosystem:

• https://discourse.slicer.org/t/extracting-data-from-vmtk-centerline/7205
• https://www.slicer.org/wiki/Modules:VMTK_in_3D_Slicer_Tutorial:_Coronary_Artery_Centerline_Extraction

If you have a way of generating centerlines, then we should be able to easily perform the tubular clipping to segment your tree meshes.

---

Since vmtk is easily installable via conda, I went ahead and played around with this. I installed everything with:

conda create -n vmtk_testing -c vmtk -c conda-forge python=3.6 itk vtk vmtk pyvista

Then I tried making the centerlines and splitting the meshes following those two tutorials. Unfortunately, I got errors and it couldn't really compute the centerlines. If you want to try, its:

vmtkcenterlines -ifile tree -ofile tree_centerlines.vtp

---

I am actually also working on generating the skeleton mesh in PyVista using a forcing function in the 3D mesh.... I'll try to wrap it up and post later

[GreKro]

@banesullivan Yes, recently I found vmtk library too and for now, I'm able to generate a skeleton of a tree (centerline with unique ids for branches). But I can't do branch clipping in vmtk, I have some errors and I can't figure out how to fix it. Now I tried to use coordinates from skeleton to clip branches in pyvista but I have no idea how to calculte properly normal for a planes. If you want to see a skeleton, I prepared txt file with xyz and few other vmtk centerline attributes, ids of branches contain CenterlineIds attribute.

pine_qsm_skeleton_cleaned.txt

[banesullivan]

Can you save those centerlines as a vtp file to preserve the line geometry?

[GreKro]

@banesullivan Yes, sure. I attached two files in vtp - first output of centerline function, second centerlines after brachextractor algorithm tree_cl.zip

[banesullivan]

Actually, check this out..

You could simply interpolate the skeleton's group ids onto the original mesh then split the mesh by a thresholding criteria

import pyvista as pv
import numpy as np

tree = pv.read("tree-segmentation/tree.ply")

import pandas as pd
df = pd.read_csv("tree-segmentation/pine_qsm_skeleton_cleaned.txt", delim_whitespace=True)
skeleton = pv.PolyData(df[['x', 'y', 'z']].values)
for k in ['MaximumInscribedSphereRadius', 'CenterlineIds',
       'TractIds', 'Blanking', 'GroupIds']:
    skeleton[k] = df[k].values
skeleton

interpolated = tree.interpolate(skeleton)

segments = pv.MultiBlock()
name = "CenterlineIds"
for v in np.unique(skeleton[name]):
    segment = interpolated.threshold([v-0.5, v+0.5], scalars=name)
    segments["{}{}".format(name,v)] = segment

segments.plot(multi_colors=True, notebook=False)

[akaszynski]

Nice solution @banesullivan!

[banesullivan]

I think that interpolation (since it uses a uniform Gaussian kernel) would be almost the same as clipping with tubes along the skeleton-like mentioned in one of the tutorials I shared above.

[GreKro]

Hi, @banesullivan I checked out your solution, It is promising... but when you look at segments this approach works not properly (at your screen look similar, you have many colors ( id values) on single branches, it should be one value per branch) Look at segment with index i.e. 2

I suppose it is related to the fixed radius (if I'm wrong, correct me). Is there the possibility to set for each point of skeleton different radius ? in this file which I prepared you have attribute 'MaximumInscribedSphereRadius' which could help to define proper radius for each point.

Best wishes,

[banesullivan]

So we can make tubes with varying radii, however these tubes aren't well formed becasue the skeleton mesh has really sharp turns.

Check out the following

import pyvista as pv
import numpy as np

tree = pv.read("tree.ply")
skeleton = pv.read("tree_cl/pine_qsm_CL_branchextractor.vtp")

# split up each group:
branches = pv.MultiBlock()
for idx in np.unique(skeleton["CenterlineIds"]):
    cell_ids = np.argwhere(skeleton["CenterlineIds"] == idx)
    branches["group{}".format(idx)] = skeleton.extract_cells(cell_ids).extract_geometry()

tubes = pv.MultiBlock()
for branch in branches:
    tube = branch.tube(scalars="MaximumInscribedSphereRadius", 
                       radius_factor=1.0)
    tubes.append(tube)

p = pv.Plotter(notebook=0, shape=(1,2))
p.add_mesh(skeleton, line_width=5)
p.subplot(0,1)
p.add_mesh(tubes, multi_colors=True, opacity=0.5)
p.add_mesh(tree, color=True,)
p.link_views()
p.show()

If we can smooth those poly lines, then the tubes would be better formed and we could then use the clip_surface filter to extract regions of the tree within those tubes. There are many ways to smooth those lines, but I'm not sure what would work best. I'll see what I can do when comping back to this

[GreKro]

@banesullivan Wow, this approach probably will work. Referring to the issue with the not well shape of the tube model, at first glance, I can see that the created tubes do not have a metric radius preserved on the scale. All the branches have the same radius. These branches from the upper part of the crown should be much more tiniest than the main stem. You can see the thickness ratio of branches on the raw tree model. I don't think so it is related to sharpy lines of the skeleton.

[akaszynski]

@banesullivan, how about we include this as a demo? tree.ply is already in the pyvista-data repo.

[banesullivan]

Hah, I didn't realize we already had these/similar data.

This would make an excellent "Advanced Demo" once we actually figure this out. I'll try to play around with it one night this week

[venugovh]

I am actually also working on generating the skeleton mesh in PyVista using a forcing function in the 3D mesh.... I'll try to wrap it up and post later

@banesullivan Does pyvista have a mesh skeletonization filter implemented? Could you please let me know?

I am trying to find methods to calculate the skeleton/centerline of an STL file and also find the maximum inscribed sphere radius within the STL file. Thank you!

A sample stl file will look like this u-bend with a radius of 0.25mm.

[koenterheegde0507]

grid.cell_arrays[vtk.vtkDataSetAttributes.GhostArrayName()] = ghost_cells i get the error: AttributeError: 'UnstructuredGrid' object has no attribute 'cell_arrays'. Is it deprecated? Can't find anything on it.

[banesullivan]

cell_arrays is cell_data now