Euclidean Voronoi labelling on the multidimensional grid

Given a set of k points (sites) S ⊂ R^d, and a finite d-dimensional grid G_n^d=(I, n)^d, we describe an efficient algorithm computing the mapping D : G_n^d → {l...,k} such that D(_p) is the index of the site closest to p under the Euclidean norm. This mapping is called the Voronoi labelling of the grid. The algorithm traverses the points of G on a "locality-preserving" space-filling curve, exploiting the correlation between the value of D on adjacent grid points to reduce computation time. The advantage of this algorithm over existing ones is that it works in any dimension, and is general-purpose, in the sense that it is easily modified to accomodate many variants of the problem, such as non-integral sites, and computation on a subset of G. The runtimes of our algorithm in two dimensions compare with those of the existing algorithms, and are better than the few existing algorithms operating in higher dimensions.