How to create a Delaunay triangulation with SciPy

A scattered coordinate set does not show which samples share a local face. Delaunay triangulation turns those coordinates into a non-overlapping simplex mesh for spatial lookup, interpolation, and mesh preparation.

SciPy's scipy.spatial.Delaunay accepts a floating-point NumPy array and delegates tessellation to Qhull. Its simplices rows contain indexes into the original point array; each row has three counterclockwise vertex indexes for two-dimensional input.

The five input points include a rectangular boundary and one interior sample, so the triangle areas can be checked against the boundary area. Two inside queries and one outside query also expose the -1 result that must be guarded before simplices is indexed.

1. Create delaunay_demo.py with the SciPy import and two-dimensional point array.

   $ cat > delaunay_demo.py <<'PY'
   import numpy as np
   from scipy.spatial import Delaunay
   
   points = np.array([
       [0.0, 0.0],
       [2.0, 0.0],
       [2.0, 2.0],
       [0.0, 2.0],
       [0.9, 0.7],
   ])
   PY

   The first four rows form the boundary, and the last row adds an interior vertex.

2. Append the triangulation and stable vertex-index report to delaunay_demo.py.

   $ cat >> delaunay_demo.py <<'PY'
   
   triangulation = Delaunay(points)
   triangles = points[triangulation.simplices]
   
   vertex_sets = sorted(
       sorted(int(index) for index in simplex)
       for simplex in triangulation.simplices
   )
   print("triangle vertex indexes:")
   for vertices in vertex_sets:
       print(f"  {vertices}")
   PY

   SciPy may return simplex rows in a different order when floating-point rounding differs. Sorting the report keeps its presentation stable without changing triangulation.simplices or assuming a fixed simplex number.

3. Append the triangle-area coverage check to delaunay_demo.py.

   $ cat >> delaunay_demo.py <<'PY'
   
   edge_a = triangles[:, 1] - triangles[:, 0]
   edge_b = triangles[:, 2] - triangles[:, 0]
   triangle_areas = 0.5 * np.abs(
       edge_a[:, 0] * edge_b[:, 1]
       - edge_a[:, 1] * edge_b[:, 0]
   )
   
   boundary = points[:4]
   next_boundary = np.roll(boundary, -1, axis=0)
   boundary_area = 0.5 * abs(
       np.sum(boundary[:, 0] * next_boundary[:, 1])
       - np.sum(boundary[:, 1] * next_boundary[:, 0])
   )
   
   np.testing.assert_allclose(triangle_areas.sum(), boundary_area)
   print(f"triangles: {len(triangles)}")
   print(f"triangulated area: {triangle_areas.sum():.3f}")
   print(f"boundary area: {boundary_area:.3f}")
   PY

   Duplicate, collinear, or nearly coplanar input can omit vertices or produce degenerate simplices. The triangulation.coplanar attribute records omitted points, while a Qhull option such as QJ perturbs input and may create zero-area triangles.

4. Append the guarded point-location check to delaunay_demo.py.

   $ cat >> delaunay_demo.py <<'PY'
   
   query_points = np.array([
       [0.2, 0.2],
       [1.7, 0.4],
       [2.3, 0.5],
   ])
   query_simplexes = triangulation.find_simplex(query_points)
   
   assert np.all(query_simplexes[:2] >= 0)
   assert query_simplexes[2] == -1
   
   print("point lookup:")
   for point, simplex_index in zip(query_points, query_simplexes):
       label = (
           "outside triangulation"
           if simplex_index == -1
           else f"inside vertices {sorted(triangulation.simplices[simplex_index].tolist())}"
       )
       print(f"  ({point[0]:.1f}, {point[1]:.1f}) -> {label}")
   PY

   An outside lookup returns -1, which Python would otherwise interpret as the last row of triangulation.simplices.

5. Run delaunay_demo.py to print the complete geometric report.

   $ python3 delaunay_demo.py
   triangle vertex indexes:
     [0, 1, 4]
     [0, 3, 4]
     [1, 2, 4]
     [2, 3, 4]
   triangles: 4
   triangulated area: 4.000
   boundary area: 4.000
   point lookup:
     (0.2, 0.2) -> inside vertices [0, 3, 4]
     (1.7, 0.4) -> inside vertices [1, 2, 4]
     (2.3, 0.5) -> outside triangulation