Generating cloud of points around a mesh

[brianjimenez]

First of all, congratulations on this amazing piece of software. I have a very specific problem to solve, hope you could give some insights on using PyVista for solving it.

Description

The goal is to generate a uniform cloud of points over a surface. I am generating the surface using a cloud of points which are atomic coordinates of a molecule. For that purpose, I've been playing with delaunay_3d and decimate filters, please see:

points = np.load("coord.npy")
point_cloud = pv.PolyData(points)
# Sanity check
#print(np.allclose(points, point_cloud.points))

surf = point_cloud.delaunay_3d(alpha=7.0, offset=4.0)

shell = surf.extract_geometry().triangulate()

decimated = shell.decimate(0.8).extract_surface().clean()

This generates a nice polygonal surface. Then, using cell centers and normals for extruding the points:

Complete code here:

import numpy as np
import pyvista as pv

points = np.load("coord.npy")
point_cloud = pv.PolyData(points)
# Sanity check
#print(np.allclose(points, point_cloud.points))

surf = point_cloud.delaunay_3d(alpha=7.0, offset=4.0)

shell = surf.extract_geometry().triangulate()

decimated = shell.decimate(0.8).extract_surface().clean()

decimated.compute_normals(cell_normals=True, point_normals=False, inplace=True)

centers = decimated.cell_centers()
centers.translate(decimated['Normals']*10.0)

p = pv.Plotter(notebook=False, shape=(1,3))
p.add_mesh(point_cloud, point_size=8.0, render_points_as_spheres=True, show_edges=True)
p.subplot(0,1)
p.add_mesh(shell)
p.subplot(0,2)
p.add_mesh(decimated)
p.add_mesh(centers, color="r", point_size=8.0, render_points_as_spheres=True)
p.link_views()
p.show()

I think extruding the polygonal surface and then using the vertices instead of cell centers might be a better description of the surface (some regions might not be uniformly sampled), but I don't know it that is possible at the moment with PyVista.

The previous algorithm was a ray-tracing approach from uniformly-distributed points over a sphere, but that approach seems to be limited in corner cases (think for example in cavities).

Example Data

Here you may find the code and NumPy data used: example.zip

Thank you in advance for your help!

[akaszynski]

Here's a slightly different way of trying to get a uniform cloud of points around your initial point cloud. It works by:

• Expanding the initial mesh using vtkWarpVector
• Try to make a uniform surface out of that mesh.
• Simply plot the uniform surface points.

The biggest (and most frustrating) limitation right now is with delaunay_3d. The surface isn't manifold, and it causes havoc with anything downstream (especially with pyacvd's uniform remesher). If you can figure out how to get delanuay_3d working (including trying out some external packages), I'd use that mesh, uniformly remesh it, and expand it in small increments (to ensure the point normals are updated in between steps).

The solution below doesn't work because of delaunay_3d...

import numpy as np
import pyvista as pv

points = np.load("coord.npy")
point_cloud = pv.PolyData(points)

surf = point_cloud.delaunay_3d(alpha=7.0, offset=4.0)
shell = surf.extract_geometry().triangulate()
shell = shell.compute_normals()

# warp each point by the normal vectors
import vtk
warp = vtk.vtkWarpVector()
warp.SetInputData(shell)
warp.SetInputArrayToProcess(0, 0, 0,
                            vtk.vtkDataObject.FIELD_ASSOCIATION_POINTS,
                            vtk.vtkDataSetAttributes.NORMALS)
warp.SetScaleFactor(10)
warp.Update()
expanded_mesh = pv.wrap(warp.GetOutput())

import pyacvd

clus = pyacvd.Clustering(expanded_mesh)
clus.subdivide(3)
clus.cluster(5000)
uniform = clus.create_mesh()
# out.plot()

p = pv.Plotter(notebook=False, shape=(1,2))
p.add_mesh(point_cloud, point_size=8.0, render_points_as_spheres=True, show_edges=True)
p.subplot(0,1)
p.add_mesh(shell)
# p.subplot(0,2)
# p.add_mesh(decimated)
p.add_points(np.asarray(uniform.points), color="r",
             point_size=8.0, render_points_as_spheres=True)
# p.add_mesh(centers, color="b", point_size=8.0, render_points_as_spheres=True)
# p.add_mesh(expanded_mesh, smooth_shading=True)
p.link_views()
p.show_bounds()
p.show()

[brianjimenez]

Thanks a lot @akaszynski for spending some time on the problem and the proposed solution. Indeed, it's a nice approach and probably doable if, instead of using the atom coordinates as a simple point cloud, I could generate a set of spheres/cubes from the atom coordinates. Kind of overkilling, but cubes (voxels) might be combined in a union and then remeshed as a surface. Does it make sense?

[brianjimenez]

I've been playing a bit with @akaszynski code, but I think there is a problem with the warp expansion using the normal. Please see in a more globular set of coordinates:

Data used: 1ppe.zip

Code is very similar:

import numpy as np
import pyvista as pv
import vtk
import pyacvd

points = np.load("coord.npy")

# Use PyVista
point_cloud = pv.PolyData(points)

# Use Delaunay 3D algorithm to get a first mesh
surf = point_cloud.delaunay_3d(alpha=8.0, offset=10.0)
shell_raw = surf.extract_geometry().triangulate()

# Decimate surface to clean artifacts
shell = shell_raw.decimate(0.8).extract_surface().clean()
shell = shell.compute_normals()

# Warp each point by the normal vectors
warp = vtk.vtkWarpVector()
warp.SetInputData(shell)
warp.SetInputArrayToProcess(0, 0, 0,
                            vtk.vtkDataObject.FIELD_ASSOCIATION_POINTS,
                            vtk.vtkDataSetAttributes.NORMALS)
warp.SetScaleFactor(20)
warp.Update()
expanded_mesh = pv.wrap(warp.GetOutput())

# Uniform mesh after remeshing
clus = pyacvd.Clustering(expanded_mesh)
clus.subdivide(3)
clus.cluster(400)
uniform = clus.create_mesh()

# Plot structures, uniform mesh and points
p = pv.Plotter(notebook=False, shape=(1,4))
# First plot is atoms and swarm centers
p.add_mesh(point_cloud, point_size=8.0, render_points_as_spheres=True, show_edges=True)
p.add_points(np.asarray(uniform.points), color="r", point_size=8.0, render_points_as_spheres=True)
# Second plot is the uniform mesh generated as a surface
p.subplot(0,1)
p.add_mesh(uniform, smooth_shading=True)
# Third plot is decimated surface
p.subplot(0,2)
p.add_mesh(shell)
# Fourth is the expanded mesh
p.subplot(0,3)
p.add_mesh(expanded_mesh)
p.link_views()
p.show()

I would expect in plot 4 a version 3 inflated, but it looks like there are some points which are incorrect. Please, could you point me to the possible error in the code? Thanks a lot!

[whoIsClownHere]

Hi everyone, I also have such a problem right now, please, could you explain what the red dots mean?