2D Voronoi Diagrams using a Method by Byers
Bruno Guerrieri
Florida A&M University
Tallahassee, FL 32312
Given N points (sites) in the plane (a 1x1 square) we would like to generate a tessellation of that domain by assigning every point in the plane to its nearest site. The process, as presented by Byers [1], is to select a particular site, call it C, and determine its Voronoi polygon. This is achieved by a) generating all segments from C to the remaining sites, b) considering all bisecting lines of each of these segments (in other words, the normal to the segment through the midpoint), c) determining all intersections points of these 'bisectors", d) eliminate those intersection points not in the 1x1 square, and e) determine which ones of the remaining intersection points are "closest" to C. Their convex hull is the Voronoi cell we are looking for. We then repeat the process for each of the remaining sites. We are mainly interested in the discovery process and, therefore, are not always concerned with efficiency or storage minimization.
Procedures
Randomly Selected Sites
Procedure for randomly generating N points (sites) in a 1x1 square. We will try to have the selected sites not"fall" too close to any of the edges.
Remove Elements from a List
Once a point has been found to belong to the convex hull, we want to eliminate it from further consideration. We also will be removing all the "interior" points that have been "passed by".
Is Intersection Point a Candidate?
It is a candidate if there is no bisecting line strictly between that intersection point and the site (xc,yc).
Elevation
The formula we use for the arctan is one that returns angles in interval (-Pi,Pi) (Consult [2]). "a" is the list of candidates and "qg" is the observation point (site whose Voronoi cell we are trying to determine) (xc,yc)
Voronoi
This is the condensed version of the process. It will be used in Main 2. In Main 1 (see below), one finds all the ingredients found in procedure Voronoi, shown one step at a time.
Main 1
Initially, let us select a certain number of random points that will become the sites. Let us record their coordinates and plot them as well. To determine the Voronoi diagram, we are first going to determine the equations of all the "bisecting" lines. Then we will determine all the points of intersection of these different lines and then determine which of these intersection points belong to the Voronoi diagram. N represents the number of sites. We will maintain a copy of the original list of sites in "original_list'. We will also have a "running_list" of sites which will dynamically change through the process.
At this moment, we will proceed slowly. Once the sequence of steps to follow has proven successful to generate a particular Voronoi cell, we will aggregate all commands into a single procedure.. We are going to select one of our original sites (the k-th one), call it (xc,yc), as the one under scrutiny at the moment. In the graphics we will identify it by circling it in red. We wish to determine its Voronoi cell. Let us find the equations for all the "bisecting" lines about (xc,yc).
Next, we calculate the equations of the lines that are perpendicular bisectors between (xc,yc), the chosen point, and each of its neighbors (as well as the four boundary lines). Once that is done, we will calculate all the possible intersection points. We are doing so by solving systems of 2x2 equations in matrix form.
Let us convert this mass of matrices into a list of ordered pairs representing the intersection points.
We need to also include the points of intersection that are the result of intersection between bisecting lines and boundary lines.
We have a problem here if some of the calculations returns no intersection points!!!!!!
We are going to filter out all the points of intersection that are outside of the [0,1] x [0,1] square and keep the remaining candidates in "filtered_pairs".
Let us show graphically the situation as it is so far. We will "draw" the sites as well as the bisecting lines and the relevant intersection points that will help use determine the Voronoi cell associated with (xc,yc) (circled in red). Please, recall that the bisecting lines you see are only the ones involving (xc,yc) with the rest of the sites. By then, the eye can pretty much figure out the edges of the polygonal cell about (xc,yc).
What needs to be done now? We need to find out which intersection points, in the filtered list, need to be retained because they are on the "same" side as (xc,yc) with respect to the bisecting lines. This is the crux of a second filtering process. An intersection point will be retained if there is no "bisecting" line strictly between it and (xc,yc). We say strictly because the intersection point under scrutiny will be itself on a bisecting line.
We would like to, once more, represent the situation graphically, circling the vertives of the polygonal cell forming the Voronoi cell. Depending on the point (xc,yc) being analyzed, it may tun out that the cell is fully contained within the square while, at other times, the boundary will "intervene". We will, therefore have to analyze subcases since the polygonal contours will be sometime "closed" and sometime "open".
Here is a little technicality. We need to "order" the vertices in a clockwise manner so as to draw the "correct" (neigbhor to neighbor) polygonal contour.We will achieve this by calculating the polar angle for each of the verticeses and use these values to order them.
What do we want to do now. It seems that whenever the cell has an edge on the boundary (whether along an edge or in a corner) all we need to do is rotate the list "corner" to the left until the first pair containing a 0 in the x-position or the y-position is on the far right. We think the second edge point will have naturally migrated to the head of the pack!
Main 2
We have automatized the process. We call the procedure Voronoi to generate all Voronoi cells, one at a time, and plot the final result. All this was developed on Maple purely for pedagogical results and within the context of a Research Experience for Undergraduates. . Neither storage nor efficiency of processes have been optimized. Note that, if you wish to run Main 2 for a different value of N, then remove the "comment indicator" # on the next line. Choosing a large value for N (say 20) will result in delays that grow exponentially with N.
To display Voronoi cell 3, look at the relevant frame.
Bibliography
[1] Byers J.A.,"Dirichlet Tessellation of Bark Beetle Spatial Attack Points", Journal of Animal Ecology, 61(1992) p. 759-768.
[2] Baker J.A. and Sluis W., "Arg", Mathematics Magazine, 70 (1997), no. 4, p. 291.
