creating a closed (voxelizable) mesh from delaunay-3d triangulation of 3D point cloud

[nwvann]

Hello,

I'm trying to combine two 3D pixel masks (separated by some distance) into a single mask that envelopes both features. I'm doing this by applying the delaunay-3D functionality to the ensemble of points defining both masks. The issue is that the resulting surface derived from the delaunay-3D algorithm is not "closed" so voxelizing it produces non-sensible results. I tried to fix the surface via pymeshfix, but the "repaired" surface looks nothing like the original. Here's what I have tried:

points3d.txt

import pyvista as pv
import numpy as np

#load in xyz coordinates
points3d=np.loadtxt('points3d.txt')

# set up the pyvista point cloud structure
cloud = pv.PolyData(points3d) 

#plot of points inside each mask, separated by some distance
cloud.plot(screenshot='pointcloud.png')

I apply the delaunay-3d triangulation...

#set alpha parameter
alpha=5

#run delaunay 3D triangulation
mesh = cloud.delaunay_3d(alpha=alpha)

#convert unstructured grid to polydata
surface=mesh.extract_surface()

surface.plot(screenshot="surface.png")

That results in a surface that makes sense to me. However, when I voxelize it, it's not closed, forcing me to use the "check_surface=False" parameterization to get it to run. It produces strange results:

#I want to fetch all points inside surface; surface is not closed so I have to add "check_surface=False"
voxels = pv.voxelize(surface,check_surface=False)
voxels.plot(screenshot='voxels.png')

I try to repair the surface via pymeshfix so that doesn't happen, but the "repaired" surface looks nothing like the original:

#I tried to repair the surface using pymeshfix.
import pymeshfix as mf

meshfix = mf.MeshFix(surface)
meshfix.repair(verbose=True)
meshfix.plot(screenshot='meshfix.png')

Any ideas how to fix this and get a new set of points within the new larger surface enveloping both original sets of points?

[akaszynski]

Played around with this today and had a tough time as well getting a manifold mesh from the input cloud. There's two problems that I can see:

1. The output from delaunay_3d isn't close enough to a manifold surface for meshfix to be able to handle.
2. meshfix isn't a good candidate for cleaning up meshes that aren't "mostly there". It seems well suited for repairing meshes that are nearly manifold instead of not even close, which is what we get from delaunay_3d.

Both scipy and open3d also generated poor surfaces. Here's my attempt with open3d

import open3d as o3d
pcd = o3d.geometry.PointCloud()
pcd.points = o3d.utility.Vector3dVector(points3d)
pcd.estimate_normals()

radius = 2
bpa_mesh = o3d.geometry.TriangleMesh.create_from_point_cloud_ball_pivoting(pcd, o3d.utility.DoubleVector([radius, radius * 2]))
print(bpa_mesh)

tri = np.asarray(bpa_mesh.triangles)
faces = np.empty((tri.shape[0], 4), dtype=np.int64)
faces[:, 0] = 3
faces[:, 1:] = tri
mesh = pv.PolyData(points3d, faces)

Given this is similar to #170, I think this is a feature that should be improved or at least documented within pyvista noting that delaunay_3d does not usually generate manifold meshes.

[banesullivan]

These data appear to have originated from a regular grid. You can see a uniform grid-like structure to the points. I used to work with data like these a ton and created a nifty filter in PVGeo (PVGeo is essentially an extension package to PyVista) for creating voxelized versions of the gridded point clouds. Effectively, the algorithm evaluates the coordinates that are unique to the grid on a rotated local coordinate system then builds out voxels around every node. I would use this filter with some caution though... it is sometimes unstable but usually works for recovering voxel models from regularly gridded points like yours.

See:

• https://pvgeo.org/examples/filters-general/voxelize-points.html#voxelize-points
• https://pvgeo.org/content/PVGeo/filters/voxelize.html#module-PVGeo.filters.voxelize
• The rotation to find the grid's local coordinate system needs to be reworked and tends to yield bad results. See Leverage PCA for RotationTool Estimation Algorithm OpenGeoVis/PVGeo#54

import pyvista as pv
import numpy as np
import PVGeo

points = np.loadtxt("points3d.txt")
pc = pv.PolyData(points)

# Use with caution
grid = PVGeo.filters.VoxelizePoints().apply(pc)

p = pv.Plotter()
p.add_mesh(grid, opacity=0.5, show_edges=True)
p.add_mesh(pc, point_size=5, color="red")
p.show_grid()
p.show()

[nwvann]

Many thanks to @akaszynski and @banesullivan for taking the time to reply. It's appreciated!

@banesullivan these points did indeed originate from a regular grid. I already had a volumetric dataset that I thresholded to extract the points inside some iso-threshold, which is what I sent in as points3d. The issue is that I now need to voxelize the "connected" surface shown in the second plot, but as @akaszynski pointed out it looks like delaunay-3D doesn't generate manifold surfaces, so I think this may be beyond the scope of pyvista. But again, thank you!

[JacobBumgarner]

Many thanks to @akaszynski and @banesullivan for taking the time to reply. It's appreciated!

@banesullivan these points did indeed originate from a regular grid. I already had a volumetric dataset that I thresholded to extract the points inside some iso-threshold, which is what I sent in as points3d. The issue is that I now need to voxelize the "connected" surface shown in the second plot, but as @akaszynski pointed out it looks like delaunay-3D doesn't generate manifold surfaces, so I think this may be beyond the scope of pyvista. But again, thank you!

Did you ever find a solution to this problem?

Thanks