# Voronoi Diagrams

Some time ago I was looking for an algorithms that can generate a ‘map like’ like pictures –  e.g. tessellation of a plane into set of more or less random polygons.    I found Voronoi diagrams –   which give very nice pictures and have many useful properties.
Most common case of Voronoi diagram is known in  an Euclidean plane,   where we have a set of points (call seeds) then Voronoi  diagram splits the plane into areas – called Voronoi cells – around each seed,   where inside each area any point is closer to that seed then to any other.   Areas are then convex polygons (for Euclidean metric). This definition is best illustrated on the picture below – the Voronoi diagram for 100 random points in range 0-100 – Voronoi cells are marked by red lines, blue points are seeds: To play with Voronoi diagrams  we can use python – recent version of scipy  already contains function for calculating Voronoi diagram (see later).   But we can use duality of Voronoi diagram to Delaunay triangulation. Both scipy and matplotlib   contain   functions for Delaunay triangulation. From a given triangulation we can get Voronoi diagram, if we connect centers  of  circumcircles  of each triangle through edges of adjacent triangles.  Also there is a pure python implementation of Fortune algorithm available on web here.

Below is the code for calculation Voronoi diagram using Delaunay triangulation function in matplotlib and using numpy for numerical calculations:

```#Inspired by code https://github.com/rougier/neural-networks/blob/master/voronoi.py from Nicolas P. Rougier

import sys
import numpy as np
import matplotlib.tri
import matplotlib.pyplot as plt
import time

def circumcircle2(T):
P1,P2,P3=T[:,0], T[:,1], T[:,2]
b = P2 - P1
c = P3 - P1
d=2*(b[:,0]*c[:,1]-b[:,1]*c[:,0])
center_x=(c[:,1]*(np.square(b[:,0])+np.square(b[:,1]))- b[:,1]*(np.square(c[:,0])+np.square(c[:,1])))/d + P1[:,0]
center_y=(b[:,0]*(np.square(c[:,0])+np.square(c[:,1]))- c[:,0]*(np.square(b[:,0])+np.square(b[:,1])))/d + P1[:,1]
return np.array((center_x, center_y)).T

def check_outside(point, bbox):
point=np.round(point, 4)
return point<bbox or point>bbox or point< bbox or point>bbox

def move_point(start, end, bbox):
vector=end-start
c=calc_shift(start, vector, bbox)
if c>0 and c<1:
start=start+c*vector
return start

def calc_shift(point, vector, bbox):
c=sys.float_info.max
for l,m in enumerate(bbox):
a=(float(m)-point[l%2])/vector[l%2]
if  a>0 and  not check_outside(point+a*vector, bbox):
if abs(a)<abs(c):
c=a
return c if c<sys.float_info.max else None

def voronoi2(P, bbox=None):
if not isinstance(P, np.ndarray):
P=np.array(P)
if not bbox:
xmin=P[:,0].min()
xmax=P[:,0].max()
ymin=P[:,1].min()
ymax=P[:,1].max()
xrange=(xmax-xmin) * 0.3333333
yrange=(ymax-ymin) * 0.3333333
bbox=(xmin-xrange, ymin-yrange, xmax+xrange, ymax+yrange)
bbox=np.round(bbox,4)

D = matplotlib.tri.Triangulation(P[:,0],P[:,1])
T = D.triangles
n = T.shape
C = circumcircle2(P[T])

segments = []
for i in range(n):
for j in range(3):
k = D.neighbors[i][j]
if k != -1:
#cut segment to part in bbox
start,end=C[i], C[k]
if check_outside(start, bbox):
start=move_point(start,end, bbox)
if  start is None:
continue
if check_outside(end, bbox):
end=move_point(end,start, bbox)
if  end is None:
continue
segments.append( [start, end] )
else:
#ignore center outside of bbox
if check_outside(C[i], bbox) :
continue
first, second, third=P[T[i,j]], P[T[i,(j+1)%3]], P[T[i,(j+2)%3]]
edge=np.array([first, second])
vector=np.array([[0,1], [-1,0]]).dot(edge-edge)
line=lambda p: (p-first)*(second-first)/(second-first)  -p + first
orientation=np.sign(line(third))*np.sign( line(first+vector))
if orientation>0:
vector=-orientation*vector
c=calc_shift(C[i], vector, bbox)
if c is not None:
segments.append([C[i],C[i]+c*vector])
return segments

if __name__=='__main__':
points=np.random.rand(100,2)*100
lines=voronoi2(points, (-20,-20, 120, 120))
plt.scatter(points[:,0], points[:,1], color="blue")
lines = matplotlib.collections.LineCollection(lines, color='red')
plt.axis((-20,120, -20,120))
plt.show()```

I was inspired by code mentioned on first line comment,  I refactor  it to leverage numpy matrix computations and also added feature to limit diagram to a bounding box.

Pictures below illustrateshow Voronoi diagram is created from triangulation –   we have  a set of seed points (blue),   then Delaunay triangulation is calculated (yellow lines), for each triangle we calculate circumcircle  –   center of circle is green cross and circle itself is grey.   Finally we connect circle centers (red lines)  and for border triangles we create line from the circumcircle center in direction perpendicular to the outer triangle edge and line extends to infinity ( to bounding box boundary in our case).    As already said  there is also function scipy.spatial.Voronoi, which returns directly vertices and edges of Voronoi diagram. Generally this function is faster then our example above (app. 2 – 3 times) – below is code that also calculates missing edges (leading to infinity)  and  limits picture to a bounding box ( similar to previous example):

```import sys
import numpy as np

def check_outside(point, bbox):
point=np.round(point, 4)
return point<bbox or point>bbox or point< bbox or point>bbox

def move_point(start, end, bbox):
vector=end-start
c=calc_shift(start, vector, bbox)
if c>0 and c<1:
start=start+c*vector
return start

def calc_shift(point, vector, bbox):
c=sys.float_info.max
for l,m in enumerate(bbox):
a=(float(m)-point[l%2])/vector[l%2]
if  a>0 and  not check_outside(point+a*vector, bbox):
if abs(a)<abs(c):
c=a
return c if c<sys.float_info.max else None

from scipy.spatial import Voronoi
def voronoi3(P, bbox=None):
P=np.asarray(P)
if not bbox:
xmin=P[:,0].min()
xmax=P[:,0].max()
ymin=P[:,1].min()
ymax=P[:,1].max()
xrange=(xmax-xmin) * 0.3333333
yrange=(ymax-ymin) * 0.3333333
bbox=(xmin-xrange, ymin-yrange, xmax+xrange, ymax+yrange)
bbox=np.round(bbox,4)
vor=Voronoi(P)
center = vor.points.mean(axis=0)
vs=vor.vertices
segments=[]
for i,(istart,iend) in enumerate(vor.ridge_vertices):
if istart<0 or iend<=0:
start=vs[istart] if istart>=0 else vs[iend]
if check_outside(start, bbox) :
continue
first,second = vor.ridge_points[i]
first,second = vor.points[first], vor.points[second]
edge= second - first
vector=np.array([-edge, edge])
midpoint= (second+first)/2
orientation=np.sign(np.dot(midpoint-center, vector))
vector=orientation*vector
c=calc_shift(start, vector, bbox)
if c is not None:
segments.append([start,start+c*vector])
else:
start,end=vs[istart], vs[iend]
if check_outside(start, bbox):
start=move_point(start,end, bbox)
if  start is None:
continue
if check_outside(end, bbox):
end=move_point(end,start, bbox)
if  end is None:
continue
segments.append( [start, end] )

return segments```

Voronoi diagrams have many different application , some of them in geography – for instance below is a map of Czech Republic, overlayed with Voronoi diagram, where seeds are breweries in that country. 