Voronoi diagram is a diagram for set $P$ of $N$ points from $\mathbb R^n $ space (generally to any set $X$). It is constructed by splitting $\mathbb R^n$ space into $N$ subspaces $A_i$, in such a way that arbitrary point $x\in\mathbb R^n$ belongs to $A_i$, only if point $p_i\in P$ is the nearest point to $x$ for given metric $d$. From the mathematical perspective one can define each subspace in the following way:

\[A_i = \{ x\in \mathbb R^n : \forall_{i\neq j} d(x,p_i) \le d(x,p_j) \}, \\]

where $d(x,y)$ is given metric. One can use an arbitrary metric. However, most common ones are Euclidean (standard) and Manhattan (sometimes referred to as taxicab) metrics. The first defined by function: $d(x,y) = \sqrt{\displaystyle\sum_{i=1}^{n}(x_i-y_i)^2} $ and the other one by: $d(x,y) = \displaystyle\sum_{i=1}^n |x_i-y_i|$. Voronoi diagrams are very useful. There can be found in many aspects of our lives: in natural sciences, engineering, geometry, computer games or even in medicine. The main goal of the post is to introduced to the reader some basics concepts and thought regarding the Voronoi diagram. I prepared some basics scripts in gnuplot and python to present to you some methods for solving this kind of problem. Let’s begin with the simple brute-force algorithm in python. From now on we will be focus only on $n=2$ case. Below a simple program is presented. Here we assume that we work in a square area: $(x,y)\in \mathbb R^2, 0 \le x,y \le $L. we generate 20 points and we sample points in a square with L/samples steps and for each we check, which point from data is the nearest.

#!/bin/python
import random
import sys

L = 10.
samples = 1000.
points = 20

data = []
for i in range(points):
    x = L * random.random()
    y = L * random.random()
    data.append([x,y])

def d(x,y):
    #euclidan
    #return (x[0]-y[0])**2 + (x[1]-y[1])**2
    #metro metric
    return abs(x[0]-y[0]) + abs(x[1]-y[1])

def min_d(x,P):
    MIN = []
    for i in P:
        MIN.append(d(x,i))
    return MIN.index(min(MIN))

for i in range(int(samples)):
    for j in range(int(samples)):
        px = i * L/samples
        py = j * L/samples
        print px,py,min_d([px,py],data)

After the procedure, one can visualize the results using gnuplot program. Our output has the following form: $x$, $y$, $i$, where $i$ is the index of subspace $A_i$.

#!/bin/gnuplot

unset key
unset colorbox
unset tics
set border 0
set size ratio 1
set xr[0:10]
set yr[0:10]
filename = '<python voronoi.py'

p filename u 1:2:3 w image

where voronoi.py is the filename of the previous python script. Here you have an example and comparison of these two different metrics on the same random points:

Fig. 1: An example Voronoi diagram for Euclidean (left) / Manhattan (right) metric.

As you can see the result with the Manhattan metric is similar to a polyline. Each area’s border is a polyline, in contrast to the Euclidean metric, where borders between neighbors’ areas are just straight lines. In many blogs, websites, people explain Voronoi diagrams by using some graphical programs. By generating cones, one can find the Voronoi diagram (for Euclidean metric) by simply searching for the intersection between them. However, I have never seen a solution to the problem, which was fully designed in gnuplot. Gnuplot is plotting and also a graphical program. Most of the community presents an example only using the cones, but similar understanding can be considered for other geometric objects, for example for pyramids to solve the problem in Manhattan metric.

Fig. 2: Cone/pyramid landscape - graphical method for finding Voronoi diagram.

The problem comes down to simply drawing the selected geometrical shapes. Then one should only map these objects into $xy$ plain. Here a gnuplot example is presented:

#!/bin/gnuplot

# square dim
L = 10
# number of points
N = 20
# radius
r = 3

# euclidean metric
# d(x,y) = sqrt(x**2+y**2)
# manhattan metric
d(x,y) = abs(x)+abs(y)
f(x,y,x0,y0,r) = (d(x-x0,y-y0)<r)? -d(x-x0,y-y0) : 1/0

# random seed
sys = rand(int(`date +%s.%N`))

# some settings
set border 0
unset tics
unset colorbox
set view 0,0
set isosamples 100
set samples 1000
set view equal xy
set surface
set hidden3d
unset key

set ur[0:L]
set vr[0:L]

array A[N];array B[N];
do for [i=1:N] { 
	A[i] = rand(0)*L
	B[i] = rand(0)*L 
}

set pm3d depthorder 
sp for [i=1:N]  '++' u 1:2:(f(x,y,A[i],B[i],r)):(i) w pm3d

Fig. 3: Maping cone (left) / pyramid (right) landscape into plain - Voronoi diagram.

In the previous code snippet, there is another variable, which can be changed: distance radius r. Changing r value, creation of the intersection between objects can be seen. In the next animations, this phenomenon was presented.

Fig. 4: Radial growth of Voronoi diagram: Euklidean (left), Manhattan (right) metric.

In this post, the simplest solutions to these kinds of problems were presented. There are many more (much more sophisticated) algorithms. Probably the famous one: Fortune’s algorithm.

AW