.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "api/gallery/geometry/triangulation/plot_02_delaunay_triangulation.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_api_gallery_geometry_triangulation_plot_02_delaunay_triangulation.py: Delaunay triangulation: empty circumcircles and fat triangles =================================================================== Triangulates a random point set with :func:`~mathematicskit.geometry.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. .. GENERATED FROM PYTHON SOURCE LINES 14-43 .. code-block:: Python 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 .. GENERATED FROM PYTHON SOURCE LINES 44-46 The empty-circle property ----------------------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 46-58 .. code-block:: Python 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}") .. rst-class:: sphx-glr-script-out .. code-block:: none 22 triangles; points strictly inside a circumcircle: 0 .. GENERATED FROM PYTHON SOURCE LINES 59-64 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. .. GENERATED FROM PYTHON SOURCE LINES 64-74 .. code-block:: Python 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}") .. rst-class:: sphx-glr-script-out .. code-block:: none Delaunay uses diagonal [1, 3] smallest angle: Delaunay 25.23 deg, other diagonal 8.53 deg both cover the same area: 7.2000 .. GENERATED FROM PYTHON SOURCE LINES 75-77 The triangulation and its empty circumcircles ----------------------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 77-97 .. code-block:: Python 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() .. image-sg:: /api/gallery/geometry/triangulation/images/sphx_glr_plot_02_delaunay_triangulation_001.png :alt: Every circumcircle is empty, Delaunay diagonal: min angle 25.2 deg, other diagonal: min angle 8.5 deg :srcset: /api/gallery/geometry/triangulation/images/sphx_glr_plot_02_delaunay_triangulation_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 0.065 seconds) .. _sphx_glr_download_api_gallery_geometry_triangulation_plot_02_delaunay_triangulation.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_02_delaunay_triangulation.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_02_delaunay_triangulation.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: plot_02_delaunay_triangulation.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_