pythonscipytriangulationdelaunaypyvista

PyVista mesh triangulation - vertex indices

I have a triangulated mesh that I generated with the Delaunay 3D function in PyVista. I would like to calculate the surface area of the mesh adding up the areas of all the triangles. Is there a way to obtain the indices of the simplices triangles from the delaunay result?

import pyvista as pv
cloud = pv.PolyData(points)
volume = cloud.delaunay_3d(alpha = 2.5)
shell = volume.extract_geometry()
shell.plot()

I know I can do it with Scipy but for whatever reason Scipy generates an incorrect mesh (and does not have attributes I can adjust in the Delaunay method):

from scipy.spatial import Delaunay
tri = Delaunay(points)
print(tri.simplices)


[[386 466 377 613]
 [159 386 377 613]
 [159 386 466 613]
 ...
 [696 709 695 691]
 [696 710 711 691]
 [696 697 711 691]]

My goal is to loop through the triangles and calculate the surface area of the mesh.

Solution

  • PyVista PolyData objects already have an area property that adds up the cell areas.

    For example, consider random points on the unit sphere:

    import numpy as np
import pyvista as pv

# random point cloud on a unit sphere
rng = np.random.default_rng()
N = 1000 # points
thetas = rng.random(N)*np.pi
phis = rng.random(N)*2*np.pi
points = np.array([
    np.sin(thetas)*np.cos(phis),
    np.sin(thetas)*np.sin(phis),
    np.cos(thetas),
]).T

# triangulate and compute total area
mesh = pv.PolyData(points)
triangulated = mesh.delaunay_3d().extract_surface()
print(triangulated.area, 4*np.pi)

    For me this printed

    12.47386243049973 12.566370614359172

    The first value is the total area of our triangulated point cloud on a sphere, and the second value is the exact surface of a perfect unit sphere. Looks good.

    Additionally, there are other similar attributes, for instance volume for watertight surfaces (which you have):

    >>> print(triangulated.volume, 4/3*np.pi)
4.127984347614561 4.1887902047863905

    In other words, you don't need the vertex indices if you're only looking for the total area.

    (If you do still want the vertices, the information is accessible via triangulated.faces. My recommendation would be to look at triangulated.faces.reshape(-1, 4)[:, 1:] which is a 2d array of shape (n_cells, 3), where each row corresponds to a given triangle and the three integers in the row are the indices of the three points forming the corresponding triangle.)