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

README.md

2d Laguerre-Voronoi diagrams

This code sample demonstrates how to compute a Laguerre-Voronoi diagram (also known as power diagram) in 2d.

Power diagrams have a wonderful property : they decompose the union of (overlapping) circles into clipped circles that don't overlap. The cells have a simple geometry, just straight lines.

It works as following

• The circles centers are considered as points associated with a positive weight : the circle's radius
• 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 power triangulation
• The power diagram is the dual of the power 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

laguerre-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)



# --- Delaunay triangulation --------------------------------------------------

def get_triangle_normal(A, B, C):
	return normalized(numpy.cross(A, B) + numpy.cross(B, C) + numpy.cross(C, A))



def get_power_circumcenter(A, B, C):
	N = get_triangle_normal(A, B, C)
	return (-.5 / N[2]) * N[:2]



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_power_triangulation(S, R):
	# Compute the lifted weighted points
	S_norm = numpy.sum(S ** 2, axis = 1) - R ** 2
	S_lifted = numpy.concatenate([S, S_norm[:,None]], axis = 1)

	# Special case for 3 points
	if S.shape[0] == 3:
		if is_ccw_triangle(S[0], S[1], S[2]):
			return [[0, 1, 2]], numpy.array([get_power_circumcenter(*S_lifted)])
		else:
			return [[0, 2, 1]], numpy.array([get_power_circumcenter(*S_lifted)])

	# Compute the convex hull of the lifted weighted points
	hull = ConvexHull(S_lifted)
	
	# Extract the Delaunay triangulation from the lower hull
	tri_list = 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 the Voronoi points
	V = numpy.array([get_power_circumcenter(*S_lifted[tri]) for tri in tri_list])

	# Job done
	return tri_list, V



# --- 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 -1, for minus infinity
     * tmax is the parameter for the right end of the segment. Can be -1, for infinity
     * Therefore, the endpoints are [A + tmin * U, A + tmax * U]
'''
def get_voronoi_cells(S, V, tri_list):
	# Keep track of which circles are included in the triangulation
	vertices_set = frozenset(itertools.chain(*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_map = { i : [] for i in vertices_set }

	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)):
		# Finite Voronoi edge
			edge = tuple(sorted((u, v)))
			if len(edge_map[edge]) == 2:
				j, k = edge_map[edge]
				if k == i:
					j, k = k, j
				
				# Compute the segment parameters
				U = V[k] - V[j]
				U_norm = norm2(U)				

				# Add the segment
				voronoi_cell_map[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_map[u].append(((edge_map[edge][0], -1), (D,  W, 0, None)))				
				voronoi_cell_map[v].append(((-1, edge_map[edge][0]), (D, -W, None, 0)))				

	# Order the segments
	def order_segment_list(segment_list):
		# 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 { i : order_segment_list(segment_list) for i, segment_list in voronoi_cell_map.items() }



# --- Plot all the things -----------------------------------------------------

def display(S, R, tri_list, voronoi_cell_map):
	# Setup
	fig, ax = plot.subplots()
	plot.axis('equal')
	plot.axis('off')	

	# Set min/max display size, as Matplotlib does it wrong
	min_corner = numpy.amin(S, axis = 0) - numpy.max(R)
	max_corner = numpy.amax(S, axis = 0) + numpy.max(R)
	plot.xlim((min_corner[0], max_corner[0]))
	plot.ylim((min_corner[1], max_corner[1]))

	# Plot the samples
	for Si, Ri in zip(S, R):
		ax.add_artist(plot.Circle(Si, Ri, fill = True, alpha = .4, lw = 0., color = '#8080f0', zorder = 1))

	# Plot the power 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 = '.9')
	line_list.set_zorder(0)
	ax.add_collection(line_list)

	# Plot the Voronoi cells
	edge_map = { }
	for segment_list in voronoi_cell_map.values():
		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 contains circles center, R contains circles radius
	sample_count = 32
	S = 5 * disc_uniform_pick(sample_count)
	R = .8 * numpy.random.random(sample_count) + .2

	# Compute the power triangulation of the circles
	tri_list, V = get_power_triangulation(S, R)

	# Compute the Voronoi cells
	voronoi_cell_map = get_voronoi_cells(S, V, tri_list)

	# Display the result
	display(S, R, tri_list, voronoi_cell_map)



if __name__ == '__main__':
	main()

LICENSE.md

MIT License

Copyright (c) 2021 Devert Alexandre

Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

@marmakoide In the context of the below image, you can see that the distance to the bigger circle(green) is smaller than distance to the smaller circle(dark blue). Is there any way in which I can have both the distances equal when constructing the edge

s.

That's would be a very different kind of diagram, not a Laguerre Voronoi diagram. I don't see a simple change of the code that would fulfill such a property. That sounds like some neat geometry problem.

Hi @marmakoide, thanks for this excellent gist. I have a question re: the function get_power_circumcenter(A, B, C).

1. As far as I understand, this is supposed to compute the Voronoi vertices as the circumcenters of Laguerre-Delauany triangles. Am I understanding this correctly?
2. Also, if 1 is yes, does this routine work out to be the same as the standard circumcenter calculation (e.g. https://en.wikipedia.org/wiki/Circumcircle#Higher_dimensions)? It's amazing how much shorter your approach is..
3. What would break in this gist if I tried to compute the power diagram of N-dimensional spheres with this code? Alternatively, how do I go about computing N-dimensional spheres?

@lkm1321

1. Yes
2. I use a N+1 dimensions convex hull to get a N dimensions triangulation, it's a trick called "lifting", and yes it's very neat because N-d convex hull are available and well-tested, it handles several types of triangulations and you get the circumcenters.
3. Not a lot of coding, but a lot of thinking :) Managing the graph structure of the Laguerre-Voronoi diagram is the hard part.

A couple of years ago, I implemented 3d Laguerre-Voronoi diagrams on Euclidean space, and spherical Laguerre-Voronoi diagrams on 2-spheres. Maybe I should post it ?

@marmakoide Yes! Personally 3D Euclidean would be very helpful. I'm sure others will find the 2-sphere version very helpful as well. Thanks!

Hi @marmakoide From your paper, how did you get the resulting image in figure 4d? Figure 4d is exactly what I want the cells to look like when they are touching each other.

@sofijadimitrijevic1 Hi, which paper are you referring to ?

@marmakoide My bad, you did not write the paper but were cited in it since they used your code. Here is the paper, I am interested in figure 4d in the paper: https://doi.org/10.3390/ijgi10030157

@marmakoide Alright looks to me like this code didn't come from you actually, I think the citation that I'm interested in is by Royer. The paper is called Mesh Generation with Alpha Complexes. Sorry about this! I think the model you built is super cool and insightful!