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.
Advertisement
Answer
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.)