## Mathematical Recreations

 MathRec Home Interactive Solution to this Problem Last Month's Problem Introductory Problem Older Problemsand Links Educational Resources Steve's Personal Page

#### Generalized Wigner-Seitz Cells (February 2002)

Note: The solution to this problem has been implemented as an interactive web page. Feel free to try it out at any time.

Given a set of distinct points inside the unit circle, construct generalized Wigner-Seitz cells defined such that the generalized Wigner-Seitz cell around each point is the area inside the unit circle which is closer to that point than to any other point.

Develop and algorithm to divide the unit circle into such cells for an arbitrary set of such points. The result should be a set of vertices where three or more cells meet, and the edges connecting those vertices. That is, the result is a planar graph with specified coordinates for the nodes

Do the same for an arbitrary set of points on the surface of the unit sphere.

Along the way you should get answers to the following:

• Given a set of three points (xi, yi) in a plane, what is the point equidistant from these points?
• Given two points (xi, yi) inside a unit circle, what points on the circle are equidistant from these points?
• Given a set of three points (θ, φ) on a sphere, what points on the sphere are equidistant from these points?

I make extensive use of vector notation in this solution. You may want to look at my summary of vector operations.

The edges between the cells are the perpendicular bisectors of the lines between the points. Once you realize this, then the problem is to find the intersections of those perpendicular bisectors. Unfortunately, this is slightly difficult.

First, let's illustrate what we are trying to achieve. If we look at a circle with nineteen points arranged for optimum packing, and divide the circle up into generalized Wigner-Seitz cells around each point, we get: Wigner-Seitz division of a circle with 19 points

Note that there are six places in the figure where four cells meet at a common point. Since our objective is to classify the topology of this arrangment of points, it makes a difference which diagonal is shorter, even if it is very close. Take a look at the following two figures:  Nearly Degenerate Vertices

These two arrangments are very close to each other, but they have different topologies. If you transform the first configuration smoothly into the second, then the cell boundaries meet in a common point at some time during the transformation. The challenge for this algorithm is to deal with the situations that are degenerate, or nearly so.

Graphically, we can distinguish between these two cases and the truly degenerate case by adding line segments connecting points whose Wigner-Seitz cells are neighbors. In the two graphs below, for example, the neighbor links indicate that the upper left and lower right cells are neighbors in the first graph, but the four cells meet at a common vertex in the second graph.  Use of neighbor links to distinguish nearly degenerate cases

There is a secondary challenge of dealing with the boundary. In particular, of dealing with the case where two (or more) cell boundaries meet on the unit circle.

I have developed an algorithm for finding generalized Wigner-Seitz cells around points in the unit circle, and it can be generalized easily for any convex boundary. I have not completed a solution for the case of circles on a sphere. I encourage you to work through it and submit your results. I have some partial results, however, which you can find below.

I've divided the algorithm into four pieces:

#### Define a polygonal boundary

The key to dealing with the degenerate situations is to make sure that you don't get bitten by roundoff errors. In addition, if you do have roundoff errors, you have to at least get a consistent result when you are done. The circular boundary poses a particularly nasty challenge this way. After some extensive work, I realized that the easiest way to deal with this problem is to put it off until the cell boundaries have been determined, then to trim the cells afterward.

This might seem to be unnecessary, but let me take just a moment to highlight some of the challenges that are avoided by this approach.

First, a linear cell boundary is determined by two points. The boundary of a cell with the unit circle is determined by two ordered points. That is, if A and B are two points on the unit circle, they can be connected by the counterclockwise arc from A to B or by the counterclookwise arc from B to A. This is fine until A and B are close enough to each other that the numerical inaccuracies of the computer get in the way. What happens when the arc from A to B is supposed to be very small, but the position of A is off by just a bit, so that the counterclockwise arc from A to B is just a bit less than 2π instead? It's important to make sure that this never happens or that it is always caught if it does.

This can be avoided if the points on the unit circle are saved as one-dimensional angles instead of as a two-dimensional Cartesian coordinates. But then you have to deal with conversion between the two systems, and you have to make sure that there are no traps with that process.

The third annoyance which is avoided by using a polygonal boundary, instead of a circular one, is that a line (a perpendicular bisector between two points in particular) can intersect an arc in zero, one or two points. In contrast, it can intersect a line segment in zero or one point.

Well, the whole issue can be avoided by defining a convex polyonal boundary outside of the unit circle. The choice of a different outer boundary won't affect any of the cell boundaries between points, and the intersections with the polygonal outer boundary can be handled in almost exactly the same way as the boundaries between cells.

I defined a square boundary at x = ±1.0625 and y = ±1.0625, but ±1.001 or any other value would work just fine. As would any other convex shape. For systems with a large number of points, it is advantageous to define a boundary which tracks the circle fairly closely. This will allow the process of defining the boundary around each cell to terminate sooner.

#### Find the generalized Wigner-Seitz cell for each point

Let A be the point for which we want to find the Wigner-Seitz cell. We start with an initial cell (defined by the outer boundary), and lop off parts of that cell that are closer to one of the other points than they are to A. We can test every point, but it's slightly more efficient to sort the other points by distance from A and then test them in order starting with the closest ones. Then we can stop as soon as we've tested all points which are close enough to further trim the cell.

At any particular step, we have a convex polygonal cell, and we have a perpendicular bisector between A and some other point B. The question is, does the perpendicular bisector of AB intersect any of the segments of the polygonal cell? If it does, then some of the current cell is closer to B than it is to A.

If we look at each vertex of the cell, we can determine whether that vertex is closer to A than it is to B. I did this by defining the vector k such that

k  =  (ayby, bxax)  =  z × (ba)

where the coordinates of A are (axay); the coordinates of B are (bxby); z is (0,0,1), the unit vector in the z direction; a is the positional vector (axay, 0), and similarly for b. From this point on, I will not make any distiction between a point P and its positional vector p.

We can also define c to be the midpoint of AB. Then we find that the perpendicular bisector l is defined by

l  =  c + u k.

A point p is on the outside of the perpendicular bisector (that is, the same side as B) if the z-component of (p – c) × k is positive. We can test each of the vertices of the current cell. The cell needs to be trimmed if and only if one or more vertices are on the outside of the perpendicular bisector.

Now, we should acknowledge the beast that is roundoff error. By classifying the vertices first, then considerng the line segments, we avoid some unlikely, but potentially troublesome numerical errors.

Consider, for example, what happens at a vertex where four cells meet. Since we are looking at the boundaries with the closest points first, we will already have defined the two sides of the cell which meet at that vertex, when we consider the perpendicular bisector between A and the third point. This perpendicular bisector should pass directly through the vertex that we established as the intersection between the other two sides.

But numerical equality is a tricky thing in most systems. We want to avoid a situation where our program decides that the perpendicular bisector intersects one side, but not the other. For that reason, test each vertex once, and define an intersection with one of the existing boundary segments as a situation where one vertex is inside the perpendicular bisector (closer to A) and the other vertex is on or outside the perpendicular bisector.

You may be wondering why I used a vector approach to decide whether a vertex was on the inside or outside of the perpendicular bisector. There is an obvious alternative—to use the Pythagorean distance formula and simply compare the distances directly. The reason that I didn't do it that way was that I found it very useful to keep track of the direction of the cell boundary. That is, I consider the cell boundary to go counterclockwise around A. If we keep this information, then we know that the first intersection (going counterclockwise) occurs when the first vertex of a boundary segment is inside the perpendicular bisector and the second vertex is on or outside the perpendicualar bisector.

So, we store the boundary of the cell as an ordered array of vertices, and we consider each perpendicular bisector to see whether it intersects the current cell at zero, one or two points. But, what if it intersects the cell at three or more points? This would be completely impossible without roundoff errors, but we should have an error handling mechanism if it occurs.

We are guaranteed to have at least one clockwise intersection if we have a counterclockwise intersection, and vice versa. My approach to handling this situation was to compare all candidate clockwise and counterclockwise intersections, and see where they fall along l. The first (most clockwise) intersection should have the smallest value of u, and the last intersection should have the largest value, and I assumed that everything in between could be cut out. I have since realized that in extreme cases where five or more cells meet at a common vertex, it might be possible to have nearly degenerate nodes in a configuration such that the counterclockwise path from the "first" intersection to the "last" intersection goes the long way around A, so that most of the cell is discarded.

I haven't revised this error handling yet, but I think the best revision in a case where there are more than two intersections is to identify the vertex with the minimum value of (p – c) × k. That vertex will be on the other side of the cell. Make sure it is not discarded.

The intersection of an existing cell boundary with the perpendicular bisector between A and B is the points which is equidistant from A, B, and the point on the other side of the cell boundary segment B′). Let's find the formula for that intersection.

Let the coordinates of A, B, and B′ be (x0, y0), (x1, y1), and (x2, y2). We can use any number of equations for the two lines and find the intersection. For example, the perpendicular bisector between A and B is given by

y  =  –
 x1 – x0 y1 – y0
(x  –
 x0 + x1 2
)  +
 x0 + x1 2

But a slightly easier approach is to use the relation

(pc) × k  =  0

where p is a point on the line; c is the midpoint of AB; k is (y0 – y1, x1 – x0) = z × (b – a); and 0 is the zero vector. This defines a plane in three dimensions where (x, y) is on the perpendicular bisector for any value of z. We find

x k1yy k1x  =  c1 × k1
x k2yy k2x  =  c2 × k2

where p = (x, y); ki = (kix, kiy); and ci × ki in this scalar context refers to the z-component of the cross product. Solving these simultaneous equations to find the intersection, we find

x (k1x k2yk2x k1y)  =  k1x (c2 × k2) – k2x (c1 × k1)
y (k1x k2yk2x k1y)  =  k1y (c2 × k2) – k2y (c1 × k1)

Let's simplify the denominator first:

 k1x k2y – k2x k1y =  (y0 – y1) (x2 – x0)  –  (y0 – y2) (x1 – x0) =  x0 y1 + x1 y2 + x2 y0 – x0 y2 – x1 y0 – x2 y1

Now, as we simplify the numerator, we find:

 k1x (c2 × k2) =  (y0 – y1) [ ½ (x0 + x2) (x2 – x0) – ½ (y0 + y2) (y0 – y2) ] =  (y0 – y1) ( x22 + y22 – x02 – y02 ) / 2 k2x (c1 × k1) =  (y0 – y2) ( x12 + y12 – x02 – y02 ) / 2

Now, we can combine terms to find

x  =
 ( x02 + y02 ) ( y1 – y2 ) + ( x12 + y12 ) ( y2 – y0 ) + ( x22 + y22 ) ( y0 – y1 ) 2 (x0 y1 + x1 y2 + x2 y0 – x0 y2 – x1 y0 – x2 y1 )

Similarly, for y we find

y  =
 ( x02 + y02 ) ( x2 – x1 ) + ( x12 + y12 ) ( x0 – x2 ) + ( x22 + y22 ) ( x1 – x0 ) 2 (x0 y1 + x1 y2 + x2 y0 – x0 y2 – x1 y0 – x2 y1 )

There is just one other trick, which is necessary to implement this algorithm. When the perpendicular bisector intersects a cell segment on the outer boundary, then I use the above formulas for x and y to find the new vertex, but I have to create phantom coordinates for (x2y2). For my simple choice of square boundaries, this is trivial, although the formula is different for each of the four sides of the square. For the right-hand side of the square (at +1.0625), the phantom coordinates are (x2y2) = (2.125 –x0y0).

#### Trim the cell to the actual boundary

After we have found the Wigner-Seitz cell inside the polygonal boundary, we can then trim it to the circular boundary. Simply check all cell segments which are boundaries with other cells (that is, they are not segments on the polygonal boundary), and replace any vertices which are outside of the circle with the corresponding vertex on the circle. Be careful to note that a boundary between two cells may lie entirely outside of the circle. Those segments must be discarded.

If the endpoints of the segment are p1 = (x1y1) and p2 = (x2y2), then we can use a parametric form for the line

p1 + u (p2p1)

and solve for u. We find that the intersection satisfies

[ x1 + u (x2  –  x1) ]2  +  [ y1 + u (y2  –  y1) ]2  =  1
u2 [ (x2  –  x1)2 + (y2  –  y1)2 ]  +  u ( 2x2 + 2y2 – 2x1 – 2y1 )  +  ( x12 + y12 – 1 )

This is quadratic in u and is easily solved using the quadratic formula. It is worth noting that this is guaranteed to have two roots for any x12 + y12 < 1. Note also that u = 0 corresponds to p = p1 and u = 1 corresponds to p = p2. When both vertices are outside of the circle, this simplifies the process of determining whether or not the link passes through the circle.

#### Reconcile the results from different cells with each other

Now just because we have an algorithm which gives a generalized Wigner-Seitz cell for each data point, it does not follow that these cells are consistent with each other. It can (and does) happen, for example, that the cell for A seems to neighbor the cell for B, but not the other way around. It is necessary to remove these inconsistencies before we can produce the final Wigner-Seitz graph.

In practice, this is surprisingly easy, although my algoritm probably is not rigorously unbreakable. One thing that helped is that I chose to keep the nodes in clockwise order for each cell. The edge going from A to B is the same as the edge going from B to A, except that the sense of clockwise and counterclockwise is switched.

I've numbered all of the vertices, and vertex numbers are associated with the clockwise and counterclockwise ends of each edge. The data points are numbered as well, and data point numbers are associated with the "inside" and "outside" of each edge. When I reconcile the individual cells with each other, I go through all of the edges and I compare each of the edges with its mirror image. That is, I compare the edge from cell A to cell B with the corresponding edge from cell B to A. If the vertex numbers do not match, then I know that I have two vertex numbers assigned to the same vertex, and I simply replace every occurence of one number with the other number.

In practice, this gives a completely self-consitent graph and correctly finds a number of degenerate vertices. I am aware, however, that it is theoretically possible to defeat this reconciliation process. (I haven't generated a graphic, so please tolerate the long text description.) Consider a situation where six cells (a, b, c, d, e, and f) meet at a common vertex. A counterclockwise path around the vertex would pass through those six cells in alphabetical order. Now suppose that the numerical program finds that cell a borders cell d as well as cells f and b, and that similarly each cell is found to border the opposite cell, as well as the adjacent cells. My algorithm will selfconsistently resolve this single degenerate vertex to two vertices.

It is not immediately clear to me how to resolve this inconsistency. I welcome responses. In the meantime, I think it's unlikely that my algorithm will fail with any actual set of points.

I have turned my algorithm into Perl code. In fact, I've turned it into a .cgi script, and you can play with it. This .cgi allows you to feed the script a set of up to twenty points and it will create the Wigner-Seitz graph. Try to break it (the algorithm, not the script—please don't hack me). I've tried to highlight the potential weaknesses.

#### Partial Results for Wigner-Seitz Division of the Sphere

The problem of defining an algorithm for Wigner-Seitz division of the Sphere is similar to the corresponding problem for division of the plane. There are two main differences, however. The first is that two "perpendicular bisectors" meet in two points, not one. The second is that two points do not determine a line segment.

I haven't worked through these issues well enough to develop a solution, but I'll put up what I have. Please send me what you have.

I've been representing points by their position vector in Cartesian coordinates. The vector is not always normalized, though, because that often adds computational load without any significant benefit. For example, the position vector for (θφ) = (π/4, π/3) is (x, y, z) = ( (√2)/4, (√6)/4, (√2)/2 ), but I might also use ( 1, √3, 2) to represent the same point on the sphere.

I use the normal to the plane of the great cirle to represent a line. The normal has direction, so it represents a direction along the great circle as well. The normal points in the direction from which the path along the great circle appears to be counterclockwise.

The perpendicular bisector between two points A and B is given by

n  =  ab

This is not a normalized vector. That is, the length of the vector is not equal to one. This represents a counterclockwise path as seen from A.

The intersection between two great circles n1 and n2 is given by

p  =  ± n1 × n2

Again, this vector is not normalized.

Two distinct points which are not antipodal (that is, on opposite sides of the sphere) determine a great circle, but as soon as you have three data points on the sphere, their cells are all wedges with antipodal nodes. It may be that the same schema that I used for the planar solution can be applied to the spherical solution. In that case, I stored each edge as two vertex indices and two data point indices. In this case, the index information would have to be supplemented with the coordinate information to reconstruct the edge, but the same process of deciding which vertices are inside of a new perpendicular bisector and which ones are outside would still work well.

Good luck. If you tackle this problem, be sure to let me know.

Send all responses to .

Thanks,
Steve