Delaunay triangulation: empty circumcircles and fat triangles

Triangulates a random point set with delaunay_triangulation(), checks Delaunay’s defining property (no point lies inside any triangle’s circumcircle), and shows on four points that choosing the Delaunay diagonal maximizes the smallest angle and avoids needle-like triangles.

import matplotlib.pyplot as plt
import numpy as np

from mathematicskit.geometry import delaunay_triangulation, polygon_area
from mathematicskit.geometry.visualizers.plots import plot_triangulation


def circumcircle(a, b, c):
    """Center and radius of the circle through three points."""
    d = 2 * (a[0] * (b[1] - c[1]) + b[0] * (c[1] - a[1]) + c[0] * (a[1] - b[1]))
    sq = [p @ p for p in (a, b, c)]
    ux = (sq[0] * (b[1] - c[1]) + sq[1] * (c[1] - a[1]) + sq[2] * (a[1] - b[1])) / d
    uy = (sq[0] * (c[0] - b[0]) + sq[1] * (a[0] - c[0]) + sq[2] * (b[0] - a[0])) / d
    center = np.array([ux, uy])
    return center, np.linalg.norm(a - center)


def min_angle(points, simplices):
    """Smallest interior angle (degrees) over all triangles."""
    worst = 180.0
    for tri in simplices:
        p = points[tri]
        for k in range(3):
            u, v = p[(k + 1) % 3] - p[k], p[(k + 2) % 3] - p[k]
            cos = u @ v / (np.linalg.norm(u) * np.linalg.norm(v))
            worst = min(worst, np.degrees(np.arccos(np.clip(cos, -1, 1))))
    return worst

The empty-circle property

rng = np.random.default_rng(0)
points = rng.uniform(0, 10, size=(15, 2))
result = delaunay_triangulation(points)

violations = 0
for tri in result.simplices:
    center, radius = circumcircle(*points[tri])
    others = np.delete(points, tri, axis=0)
    violations += int(np.sum(np.linalg.norm(others - center, axis=1) < radius - 1e-9))
print(f"{result.simplices.shape[0]} triangles; points strictly inside a circumcircle: {violations}")

22 triangles; points strictly inside a circumcircle: 0

Delaunay maximizes the smallest angle

Four points can be triangulated in two ways, by either diagonal. The Delaunay choice avoids the thin pair of triangles and keeps the other diagonal’s endpoints out of each circumcircle.

quad = np.array([[0.0, 0.0], [4.0, -0.6], [8.0, 0.0], [4.0, 1.2]])
delaunay_quad = delaunay_triangulation(quad).simplices
diagonal = set(delaunay_quad[0]) & set(delaunay_quad[1])
print(f"Delaunay uses diagonal {sorted(int(i) for i in diagonal)}")
# the other triangulation uses the opposite diagonal
other_quad = np.array([[0, 1, 2], [0, 2, 3]]) if diagonal == {1, 3} else np.array([[0, 1, 3], [1, 2, 3]])
print(f"smallest angle: Delaunay {min_angle(quad, delaunay_quad):.2f} deg, other diagonal {min_angle(quad, other_quad):.2f} deg")
print(f"both cover the same area: {polygon_area(quad):.4f}")

Delaunay uses diagonal [1, 3]
smallest angle: Delaunay 25.23 deg, other diagonal 8.53 deg
both cover the same area: 7.2000

The triangulation and its empty circumcircles

fig, axes = plt.subplots(1, 3, figsize=(14, 5), gridspec_kw={"width_ratios": [1.2, 1, 1]})
plot_triangulation(result, ax=axes[0])
for tri in result.simplices:
    center, radius = circumcircle(*points[tri])
    axes[0].add_patch(plt.Circle(center, radius, fill=False, color="tab:orange", lw=0.6, alpha=0.6))
axes[0].set_xlim(-2, 12)
axes[0].set_ylim(-2, 12)
axes[0].set_title("Every circumcircle is empty")
for ax, simplices, label in [(axes[1], delaunay_quad, "Delaunay diagonal"), (axes[2], other_quad, "other diagonal")]:
    ax.triplot(quad[:, 0], quad[:, 1], simplices, color="gray")
    ax.scatter(quad[:, 0], quad[:, 1], color="steelblue", zorder=2)
    for tri in simplices:
        center, radius = circumcircle(*quad[tri])
        ax.add_patch(plt.Circle(center, radius, fill=False, color="tab:orange", lw=0.8))
    ax.set_xlim(-2, 10)
    ax.set_ylim(-7, 7)
    ax.set_aspect("equal")
    ax.set_title(f"{label}: min angle {min_angle(quad, simplices):.1f} deg")
plt.show()