|
import itertools |
|
import numpy |
|
from scipy.spatial import ConvexHull |
|
|
|
from matplotlib.collections import LineCollection |
|
from matplotlib import pyplot as plot |
|
|
|
|
|
|
|
# --- Misc. geometry code ----------------------------------------------------- |
|
|
|
''' |
|
Pick N points uniformly from the unit disc |
|
This sampling algorithm does not use rejection sampling. |
|
''' |
|
def disc_uniform_pick(N): |
|
angle = (2 * numpy.pi) * numpy.random.random(N) |
|
out = numpy.stack([numpy.cos(angle), numpy.sin(angle)], axis = 1) |
|
out *= numpy.sqrt(numpy.random.random(N))[:,None] |
|
return out |
|
|
|
|
|
|
|
def norm2(X): |
|
return numpy.sqrt(numpy.sum(X ** 2)) |
|
|
|
|
|
|
|
def normalized(X): |
|
return X / norm2(X) |
|
|
|
|
|
|
|
def get_circumcenter(A, B, C): |
|
U, V = B - A, C - A |
|
U_norm, V_norm = norm2(U), norm2(V) |
|
U /= U_norm |
|
V /= V_norm |
|
UdV = numpy.dot(U, V) |
|
|
|
D = numpy.array([U_norm - UdV * V_norm, V_norm - UdV * U_norm]) |
|
D /= 2. * (1. - UdV ** 2) |
|
|
|
center = D[0] * U + D[1] * V |
|
|
|
return A + center |
|
|
|
|
|
|
|
# --- Delaunay triangulation -------------------------------------------------- |
|
|
|
def is_ccw_triangle(A, B, C): |
|
M = numpy.concatenate([numpy.stack([A, B, C]), numpy.ones((3, 1))], axis = 1) |
|
return numpy.linalg.det(M) > 0 |
|
|
|
|
|
|
|
def get_delaunay_triangulation(S): |
|
# Special case for 3 points |
|
if S.shape[0] == 3: |
|
if is_ccw_triangle(S[0], S[1], S[2]): |
|
return [[0, 1, 2]] |
|
else: |
|
return [[0, 2, 1]] |
|
|
|
# Compute the lifted points |
|
S_norm = numpy.sum(S ** 2, axis = 1) |
|
S_lifted = numpy.concatenate([S, S_norm[:,None]], axis = 1) |
|
|
|
# Compute the convex hull of the lifted points |
|
hull = ConvexHull(S_lifted) |
|
|
|
# Extract the Delaunay triangulation from the lower hull, all triangles are ccw |
|
return tuple([a, b, c] if is_ccw_triangle(S[a], S[b], S[c]) else [a, c, b] for (a, b, c), eq in zip(hull.simplices, hull.equations) if eq[2] <= 0) |
|
|
|
|
|
|
|
# --- Compute Voronoi cells --------------------------------------------------- |
|
|
|
''' |
|
Compute the segments and half-lines that delimits each Voronoi cell |
|
* The segments are oriented so that they are in CCW order |
|
* Each cell is a list of (i, j), (A, U, tmin, tmax) where |
|
* i, j are the indices of two ends of the segment. Segments end points are |
|
the circumcenters. If i or j is set to None, then it's an infinite end |
|
* A is the origin of the segment |
|
* U is the direction of the segment, as a unit vector |
|
* tmin is the parameter for the left end of the segment. Can be None, for minus infinity |
|
* tmax is the parameter for the right end of the segment. Can be None, for infinity |
|
* Therefore, the endpoints are [A + tmin * U, A + tmax * U] |
|
''' |
|
def get_voronoi_cells(S, tri_list): |
|
# Compute the circumcenter of each Delaunay triangle |
|
V = numpy.array([get_circumcenter(*S[tri]) for tri in tri_list]) |
|
|
|
# Keep track of which edge separate which triangles |
|
edge_map = { } |
|
for i, tri in enumerate(tri_list): |
|
for edge in itertools.combinations(tri, 2): |
|
edge = tuple(sorted(edge)) |
|
if edge in edge_map: |
|
edge_map[edge].append(i) |
|
else: |
|
edge_map[edge] = [i] |
|
|
|
# For each triangle |
|
voronoi_cell_list = [[] for i in range(S.shape[0])] |
|
|
|
for i, (a, b, c) in enumerate(tri_list): |
|
# For each edge of the triangle |
|
for u, v, w in ((a, b, c), (b, c, a), (c, a, b)): |
|
edge = tuple(sorted((u, v))) |
|
|
|
# Finite Voronoi edge |
|
if len(edge_map[edge]) == 2: |
|
j, k = edge_map[edge] |
|
if k == i: |
|
j, k = k, j |
|
|
|
# Compute the segment [Vj, Vk] parameters |
|
U = V[k] - V[j] |
|
U_norm = norm2(U) |
|
|
|
# Add the segment |
|
voronoi_cell_list[u].append(((j, k), (V[j], U / U_norm, 0, U_norm))) |
|
else: |
|
# Infinite Voronoi edge |
|
# Compute the segment parameters |
|
A, B, C, D = S[u], S[v], S[w], V[i] |
|
U = normalized(B - A) |
|
I = A + numpy.dot(D - A, U) * U |
|
W = normalized(I - D) |
|
if numpy.dot(W, I - C) < 0: |
|
W = -W |
|
|
|
# Add the segment |
|
voronoi_cell_list[u].append(((edge_map[edge][0], None), (D, W, 0, None))) |
|
voronoi_cell_list[v].append(((None, edge_map[edge][0]), (D, -W, None, 0))) |
|
|
|
# Order in-place the segments in ccw order |
|
def order_segment_list(segment_list): |
|
if len(segment_list) > 0: |
|
# Pick the first element |
|
first = min((seg[0][0], i) for i, seg in enumerate(segment_list))[1] |
|
|
|
# In-place ordering |
|
segment_list[0], segment_list[first] = segment_list[first], segment_list[0] |
|
for i in range(len(segment_list) - 1): |
|
for j in range(i + 1, len(segment_list)): |
|
if segment_list[i][0][1] == segment_list[j][0][0]: |
|
segment_list[i+1], segment_list[j] = segment_list[j], segment_list[i+1] |
|
break |
|
|
|
# Job done |
|
return segment_list |
|
|
|
# Job done |
|
return [order_segment_list(segment_list) for segment_list in voronoi_cell_list] |
|
|
|
|
|
|
|
# --- Plot all the things ----------------------------------------------------- |
|
|
|
def display(S, tri_list, voronoi_cell_list): |
|
# Setup |
|
fig, ax = plot.subplots() |
|
plot.axis('equal') |
|
plot.axis('off') |
|
|
|
# Plot the samples |
|
plot.scatter(S[:,0], S[:,1], lw = 0., c = 'k', zorder = 1) |
|
|
|
# Plot the Delaunay triangulation |
|
edge_set = frozenset(tuple(sorted(edge)) for tri in tri_list for edge in itertools.combinations(tri, 2)) |
|
line_list = LineCollection([(S[i], S[j]) for i, j in edge_set], lw = 1., colors = '.8') |
|
line_list.set_zorder(0) |
|
ax.add_collection(line_list) |
|
|
|
# Plot the Voronoi cells |
|
edge_map = { } |
|
for segment_list in voronoi_cell_list: |
|
for edge, (A, U, tmin, tmax) in segment_list: |
|
edge = tuple(sorted(edge)) |
|
if edge not in edge_map: |
|
if tmax is None: |
|
tmax = 10 |
|
if tmin is None: |
|
tmin = -10 |
|
|
|
edge_map[edge] = (A + tmin * U, A + tmax * U) |
|
|
|
line_list = LineCollection(edge_map.values(), lw = 1., colors = 'k') |
|
line_list.set_zorder(0) |
|
ax.add_collection(line_list) |
|
|
|
# Job done |
|
plot.show() |
|
|
|
|
|
|
|
# --- Main entry point -------------------------------------------------------- |
|
|
|
def main(): |
|
# Generate samples |
|
S = 10. * disc_uniform_pick(256) |
|
|
|
# Compute the Delaunay triangulation of the points |
|
tri_list = get_delaunay_triangulation(S) |
|
|
|
# Compute the Voronoi cells |
|
voronoi_cell_list = get_voronoi_cells(S, tri_list) |
|
|
|
# Display the result |
|
display(S, tri_list, voronoi_cell_list) |
|
|
|
|
|
|
|
if __name__ == '__main__': |
|
main() |
Hi,
thank you for the code but when I run it I got:
File "C:/Users/yanivv/laguerre-voronoi-2d.py", line 144, in order_segment_list
first = min((seg[0][0], 0) for i, seg in enumerate(segment_list))[1]
TypeError: '<' not supported between instances of 'NoneType' and 'int'
can you assist?
R,
Yaniv