Getting coordinates of surface nodes using pyvista

Question

I'm wondering if anyone could help me figure out how to apply pyvista to extract the surface nodes of a 3D object. For example, suppose I have a collection of points that builds out a sphere, including 'interior' and 'surface' points:

import numpy as np
import matplotlib.pyplot as plt

N = 50
max_rad = 1
thetavec = np.linspace(0,np.pi,N)
phivec = np.linspace(0,2*np.pi,2*N)
[th, ph] = np.meshgrid(thetavec,phivec)
R = np.random.rand(*th.shape) * max_rad

x = R*np.sin(th)*np.cos(ph)
y = R*np.sin(th)*np.sin(ph)
z = R*np.cos(th)

ax = plt.axes(projection='3d')
ax.plot(x.flatten(), y.flatten(), z.flatten(), '*')

Now I'd like to apply pyvista's extract_surface to locate the 'nodes' that live on the surface, together with their coordinates. That is, I'd like for extract_surface to return an array or dataframe of the coordinates of the surface points. I've tried to build a polydata object just with the vertices above (see link and section 'Initialize with just vertices')

Any help is much appreciated. Thanks!

Answer 1

Since you've confirmed in a comment that you're looking for a convex hull, you can do this using the delaunay_3d() filter. The output of the triangulation is an UnstructuredGrid that contains a grid of tetrahedra that fills the convex hull of you mesh. Calling extract_surface() on this space-filling mesh will give you the actual exterior, i.e. the convex hull:

import numpy as np
import pyvista as pv

# your example data
N = 50
max_rad = 1
thetavec = np.linspace(0,np.pi,N)
phivec = np.linspace(0,2*np.pi,2*N)
[th, ph] = np.meshgrid(thetavec,phivec)
R = np.random.rand(*th.shape) * max_rad

x = R*np.sin(th)*np.cos(ph)
y = R*np.sin(th)*np.sin(ph)
z = R*np.cos(th)

# create a PyVista point cloud (in a PolyData)
points = np.array([x, y, z]).reshape(3, -1).T  # shape (n_points, 3)
cloud = pv.PolyData(points)

# extract surface by Delaunay triangulation to get the convex hull
convex_hull = cloud.delaunay_3d().extract_surface()  # contains faces
surface_points = convex_hull.cast_to_pointset()  # only points

# check what we've got
surface_points.plot(
    render_points_as_spheres=True,
    point_size=10,
    background='paleturquoise',
    scalar_bar_args={'color': 'black'},
)

(On older PyVista versions where PolyData.cast_to_pointset() is not available, one can convex_hull.extract_points(range(convex_hull.n_points))).

The result looks like this:

Playing around with this interactively it's obvious that it only contains points from the convex hull (i.e. it doesn't contain interior points).

Also note the colouring: the scalars used are called 'vtkOriginalPointIds' which are what you would actually expect if you tried to guess: it is the index of each point in the original point cloud. So we can use these scalars to extract the indices of the points making up the point cloud:

# grab original point indices
surface_point_inds = surface_points.point_data['vtkOriginalPointIds']
# confirm that the indices are correct
print(np.array_equal(surface_points.points, cloud.points[surface_point_inds, :]))
# True

Of course if you don't need to identify the surface points in the original point cloud then you can just use surface_points.points or even convex_hull.points to get a standalone array of convex hull point coordinates.