Compute and display a Voronoi diagram, only relying on a 3d convex hull routine. The Voronoi cells are guaranted to be consistently oriented.

README.md

2d Voronoi diagrams

This code sample demonstrates how to compute a Voronoi diagram in 2d.

It works as following

• Transform the points (called here lifting the points) from 2d to 3d.
• Compute the convex hull of the lifted points
• The lower enveloppe of the convex hull gives us the Delaunay triangulation
• The Voronoi diagram is the dual of the Delaunay triangulation.

The complexity of this algorithm is O(n log(n)), with most of the heavy lifting done by the convex hull routine.

Prerequisites

To run this sample, you will need

• Python 2.7 or above
• Numpy
• Scipy
• Matplotlib, 2.0 or above

voronoi-2d.py

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