+from __future__ import print_function
+from __future__ import absolute_import
+from __future__ import division
+
+import random
+
+from compas.geometry import centroid_points
+from compas.geometry import distance_point_point
+from compas.geometry import add_vectors
+from compas.geometry import bounding_box
+
+from compas.geometry import is_point_in_polygon_xy
+from compas.geometry import is_point_in_triangle_xy
+from compas.geometry import is_point_in_circle_xy
+
+from compas.geometry import circle_from_points_xy
+
+
+[docs]def delaunay_from_points(points, boundary=None, holes=None, tiny=1e-12):
+
"""Computes the delaunay triangulation for a list of points.
+
+
Parameters
+
----------
+
points : sequence[[float, float, float] | :class:`compas.geometry.Point`]
+
XYZ coordinates of the original points.
+
boundary : sequence[[float, float, float] | :class:`compas.geometry.Point`] | :class:`compas.geometry.Polygon`, optional
+
List of ordered points describing the outer boundary.
+
holes : sequence[sequence[[float, float, float] | :class:`compas.geometry.Point`] | :class:`compas.geometry.Polygon`], optional
+
List of polygons (ordered points describing internal holes.
+
+
Returns
+
-------
+
list[[int, int, int]]
+
The faces of the triangulation.
+
Each face is a triplet of indices referring to the list of point coordinates.
+
+
Notes
+
-----
+
For more info, see [1]_.
+
+
References
+
----------
+
.. [1] Sloan, S. W., 1987 *A fast algorithm for constructing Delaunay triangulations in the plane*
+
Advances in Engineering Software 9(1): 34-55, 1978.
+
+
Examples
+
--------
+
>>>
+
+
"""
+
from compas.datastructures import Mesh
+
from compas.datastructures import trimesh_swap_edge
+
+
def super_triangle(coords, ccw=True):
+
centpt = centroid_points(coords)
+
bbpts = bounding_box(coords)
+
dis = distance_point_point(bbpts[0], bbpts[2])
+
dis = dis * 300
+
v1 = (0 * dis, 2 * dis, 0)
+
v2 = (1.73205 * dis, -1.0000000000001 * dis, 0) # due to numerical issues
+
v3 = (-1.73205 * dis, -1 * dis, 0)
+
pt1 = add_vectors(centpt, v1)
+
pt2 = add_vectors(centpt, v2)
+
pt3 = add_vectors(centpt, v3)
+
if ccw:
+
return pt1, pt3, pt2
+
return pt1, pt2, pt3
+
+
mesh = Mesh()
+
+
# to avoid numerical issues for perfectly structured point sets
+
points = [
+
(
+
point[0] + random.uniform(-tiny, tiny),
+
point[1] + random.uniform(-tiny, tiny),
+
0.0,
+
)
+
for point in points
+
]
+
+
# create super triangle
+
pt1, pt2, pt3 = super_triangle(points)
+
+
# add super triangle vertices to mesh
+
n = len(points)
+
super_keys = n, n + 1, n + 2
+
+
mesh.add_vertex(super_keys[0], {"x": pt1[0], "y": pt1[1], "z": pt1[2]})
+
mesh.add_vertex(super_keys[1], {"x": pt2[0], "y": pt2[1], "z": pt2[2]})
+
mesh.add_vertex(super_keys[2], {"x": pt3[0], "y": pt3[1], "z": pt3[2]})
+
+
mesh.add_face(super_keys)
+
+
# iterate over points
+
for key, point in enumerate(points):
+
# newtris should be intialised here
+
+
# check in which triangle this point falls
+
for fkey in list(mesh.faces()):
+
abc = mesh.face_coordinates(fkey)
+
+
if is_point_in_triangle_xy(point, abc, True):
+
# generate 3 new triangles (faces) and delete surrounding triangle
+
key, newtris = mesh.insert_vertex(fkey, key=key, xyz=point, return_fkeys=True)
+
break
+
+
while newtris:
+
fkey = newtris.pop()
+
+
face = mesh.face_vertices(fkey)
+
i = face.index(key)
+
u = face[i - 2]
+
v = face[i - 1]
+
+
nbr = mesh.halfedge[v][u]
+
+
if nbr is not None:
+
a, b, c = mesh.face_coordinates(nbr)
+
circle = circle_from_points_xy(a, b, c)
+
+
if is_point_in_circle_xy(point, circle):
+
fkey, nbr = trimesh_swap_edge(mesh, (u, v))
+
newtris.append(fkey)
+
newtris.append(nbr)
+
+
# Delete faces adjacent to supertriangle
+
for key in super_keys:
+
mesh.delete_vertex(key)
+
+
# Delete faces outside of boundary
+
if boundary:
+
for fkey in list(mesh.faces()):
+
centroid = mesh.face_centroid(fkey)
+
if not is_point_in_polygon_xy(centroid, boundary):
+
mesh.delete_face(fkey)
+
+
# Delete faces inside of inside boundaries
+
if holes:
+
for polygon in holes:
+
for fkey in list(mesh.faces()):
+
centroid = mesh.face_centroid(fkey)
+
if is_point_in_polygon_xy(centroid, polygon):
+
mesh.delete_face(fkey)
+
+
return [mesh.face_vertices(fkey) for fkey in mesh.faces()]
+
+
+# def voronoi_from_delaunay(delaunay):
+# """Construct the Voronoi dual of the triangulation of a set of points.
+
+# Parameters
+# ----------
+# delaunay : Mesh
+# A delaunay mesh.
+
+# Returns
+# -------
+# Mesh
+# The corresponding voronoi mesh.
+
+# Warnings
+# --------
+# This function does not work properly if all vertices of the delaunay
+# are on the boundary.
+
+# Examples
+# --------
+# .. plot::
+# :include-source:
+
+# from compas.datastructures import Mesh
+# from compas.datastructures import trimesh_remesh
+# from compas.geometry import delaunay_from_points
+# from compas.geometry import voronoi_from_delaunay
+# from compas.geometry import pointcloud_xy
+# from compas_plotters import MeshPlotter
+
+# points = pointcloud_xy(10, (0, 10))
+# faces = delaunay_from_points(points)
+# delaunay = Mesh.from_vertices_and_faces(points, faces)
+
+# trimesh_remesh(delaunay, 1.0, allow_boundary_split=True)
+
+# points = [delaunay.vertex_coordinates(key) for key in delaunay.vertices()]
+# faces = delaunay_from_points(points)
+# delaunay = Mesh.from_vertices_and_faces(points, faces)
+
+# voronoi = voronoi_from_delaunay(delaunay)
+
+# lines = []
+# for u, v in voronoi.edges():
+# lines.append({
+# 'start': voronoi.vertex_coordinates(u, 'xy'),
+# 'end' : voronoi.vertex_coordinates(v, 'xy'),
+# 'width': 1.0
+# })
+
+# plotter = MeshPlotter(delaunay)
+
+# plotter.draw_lines(lines)
+
+# plotter.draw_vertices(
+# radius=0.075,
+# facecolor={key: '#0092d2' for key in delaunay.vertices() if key not in delaunay.vertices_on_boundary()})
+
+# plotter.draw_edges(color='#cccccc')
+
+# plotter.show()
+
+# """
+# voronoi = mesh_dual(delaunay)
+
+# for key in voronoi.vertices():
+# a, b, c = delaunay.face_coordinates(key)
+# center, radius, normal = circle_from_points_xy(a, b, c)
+# voronoi.vertex[key]['x'] = center[0]
+# voronoi.vertex[key]['y'] = center[1]
+# voronoi.vertex[key]['z'] = center[2]
+
+# return voronoi
+