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.

Hi, I am interesting in this code.
I have a question is the convex hall calculated as:

R,
Yaniv

I am sorry, I don't understand what your question is. My code uses a convex hull routine to compute a power triangulation. Note the a convex hull routine can also be used to compute a Delaunay triangulation, and thus a Voronoi diagram.

Hi marmakoide
I am getting the following error running your example code (Spyder python 3.7):

voronoi_cell_map = get_voronoi_cells(S, V, tri_list)
Traceback (most recent call last):

File "", line 1, in
voronoi_cell_map = get_voronoi_cells(S, V, tri_list)

File "E:/py_code/laguerre-voronoi-2d.py", line 157, in get_voronoi_cells
return { i : order_segment_list(segment_list) for i, segment_list in voronoi_cell_map.iteritems() }

AttributeError: 'dict' object has no attribute 'iteritems'

I am guessing that the error relates to the migration of dict.iteritems to dict.items between python 2 and 3. However, even updating the method to items does not seem to do the trick.
Thanks

Hi mbex9ts3,

I updated the code snippet so that it works with Python 3. I tested it, It Works On My Computer (tm).
Thank you to notify me about this issue.

Cheers,
Marmakoide

Thanks Marmakoide! I tested the code using Python 2.7 on another machine and it was working fine. I will give it a try with Py 3.

Incidentally, I was trying to pick apart the storage medium ('voronoi_cell_map') and was wondering how to pull out the individual polygons for each weighted vertex (S, R1,...,n) or the voronoi edge list (i.e. the output from voronoi(x,y) in Matlab). I am a bit new to Python (I use Matlab mostly) and am not very familiar with the 'dict' datatype tbh.

Any advice would be much appreciated!

The function get_voronoi_cells returns a map associating the index i of an input point to a list of edges that defines the corresponding voronoi cell. Each edge is defined as (i, j), (A, U, tmin, tmax), where i and j are the edge vertices, other parameters define the segment parameters.

Hi marmakoide,

Could you tell me if there is a way to get this to compute finite closed polygons circumscribing the circles?

@danvip10 I'm not sure of what do you mean by "finite closed polygons circumscribing the circles".

You can use that diagram to find a set of closed polygons enclosing all the circles, each polygon containing one circle not contained by any other polygon. You just compute that diagram with 4 additional circles with radius 0 around the set of input circles. This will enforce having closed polygons for the input circles.

Hi marmakoide,

I have a question regarding the boundary edges. The end points of edges at the boundary seems to go to infinity. If I want to restrict the cells within [0,1] domain, how should I modify the code? Thanks in advance.

@xchhuang I can suggest two approaches

1. Simple idea, complicated and expensive to implement

   • Compute the bounding box of your points
   • Clip every cell with the bounding box using a polygon clipping algorithm

2. Complicated idea, simple and cheap to implement

   • Compute the bounding box of your points
   • Scale the bounding box a bit, say, by a factor 2
   • Compute the diagram by including the 4 corners of the bounding box, give them a null weight
   • Ignore all the 4 cells associated with the 4 corners

Method 2 does not give you any guaranties about the shape of the cells, but it guaranties you that no
cell vertices is at infinity.

Hi @marmakoide,
I would like to use your code for further work in the geospatial computing field. I am planning to make my work open source. I wanted to ask your permission for using your code and also could I request you to put a license on this Gist, so it is clear when I use your work.
Thank you very much
Sunayana

@sunyana You can use this code as you wish, even for closed source work, as long as I'm free of any responsibility. In about 24 hours, I'll put a MIT license on this repository, to make this as official as it can be.

@sunyana You can use this code as you wish, even for closed source work, as long as I'm free of any responsibility. In about 24 hours, I'll put a MIT license on this repository, to make this as official as it can be.

@marmakoide : Thank you very much :)!

Hi @marmakoide,

I'm a current fourth year university student trying to use your code for a geospatial research project I'm working on and had a few questions regarding how your code works (and more so the properties of a Laguerre-Voronoi diagram).

1. How are the bisecting lines determined when two circles do no overlap? Initially I believed it to be the midpoint between the outer edges of the two circles but have found that not to be the case.

2. Is there a source you used to learn more about the Power Diagram and the algorithm needed to create this?

Thanks,
Jordan

@21jordanmorris

I found various course materials on the topic with Google. My previous experience implementing computational geometry stuffs helped me to understand those materials. I can't quote from memory what exact materials I used, but the slides bellow are a good start, with answer to your first question on slide 55.

http://www.ens-lyon.fr/LIP/Arenaire/ERGeoAlgo/JDB-ens-lyon-I.pdf

Hi @marmakoide,

Can 2 or more circle centres belong to a Voronoi cell?

Thanks, Bryan

@bbpajarito Yes, if circle A fully contains circle B, then they will be in the same Voronoi cell

Hi, thanks for your wonderful code! I find two minor issues:

1. If three points are on a line, for example, if the point set is [(0,0), (0,1), (0,2), (1,1)], it will return a triangle made by these three points and cause a numerical issue that gives NaN afterward. I think a quick co-linear checking when constructing tri_list will fix this.
2. If two infinite-long Voronoi edges start from a same point, the display will not display correctly since on line 190, they will have the same 'edge' and will be skipped. One of the case is when point set is [(0, 0), (0, 1), (1, 0), (1, 1)]. I think some naming rules for Voronoi edge can solve this quickly.

Is there a way to modify this so that the voronoi edges are medial axis between two circles ?

@jayakvenu I am not sure to understand. Voronoi edges ARE medial axis between two circles

@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!