Breakpoint Convergence in Fortune's Algorithm

Question

I am implementing Fortune's sweepline algorithm for computing Voronoi diagrams. My primary reference is "Computational Geometry: Algorithms and Applications" by de Berg et al., and while their coverage of the topic is very clear, they pass over several small but important details that I have been having trouble working out myself. I've searched the web for help, but other websites either give an even higher overview than the textbook, or give the exact same pseudocode provided by the book.

I need a way to determine whether a pair of breakpoints determined by a triple of arcs on the beach line converges or diverges, in order to detect upcoming circle events. It seems that to make a decision I would need knowledge about the shape of the Voronoi cell edges that the breakpoints trace out as Fortune's algorithm progresses. For example, if I could find the slope of the edge traced by a breakpoint I could calculate where the two lines formed by the breakpoints and their respective slopes intersect, and decide whether they converge based on that result. However, I have no idea how to get information on the slopes, only the current position of the breakpoints.

The only information I have to work with is the x,y location of the three sites and the current y-coordinate of the sweepline (I'm using a horizontal sweepline).

Actually, I do have one idea for determining convergence. Given two sites, the breakpoint between the two sections of the beachline they define is governed only by the current position of the sweep line. I thought about recording the position of the two breakpoints, temporarily advancing the sweep line a small amount, and recording their new positions. Because edges in a normal Voronoi diagram do not curve, if the distance between the new pair of breakpoints is less than the distance between the old pair, then the breakpoints converge; otherwise, they diverge. But this seems both dangerous (I have no idea if it always works) and ugly. Surely there must be a better way.

Any ideas would be appreciated, and pseudocode (in a C#-like syntax if possible) especially so. Also I am aware that there are computational geometry libraries that I could use to get Voronoi diagrams, but this is a personal learning exercise, so I want to implement all parts of the algorithm myself.

Answer 1

Welcome Drake. I implemented it by checking whether the breakpoints physically converge on the circle center in a 'fictitious' increment of the sweepline position. This actually complicates itself a bit because in certain cases the circle center can be almost or precisely at the sweepline position, so the sweepline increment needs to be proportional to the difference between the current sweepline position and the circle center generated as you recommend.

Say:

1. currentSweeplineY = 1.0f; circleCenterY = 0.5f (and we are moving downwards, i.e. in the decreasing y direction).

2. Set sweepYIncrement = (circleCenterY - currentSweepLineY) / 10.0f (the 10.0f divisor is arbitrarily chosen).

3. Check new breakpoint positions at new sweepline position.

4. Check distance to circle center.

5. If both distances are less than current distances, the breakpoints converge.

I know this is probably very expensive, since you have to calculate the breakpoint positions multiple times, but I'm confident it takes care of all possible cases.

Anyway, I'm finding serious issues with floating point precision error elsewhere in the algorithm. Definitely not as straightforward as I thought initially.

Answer 2

So this is rather embarassing, but after sleeping on the problem the answer seems obvious. I'm writing this to hopefully help students in the future with the same question as me.

The Voronoi edge between two sites perpendicularly bisects the (imaginary) line segment connecting the sites. You could derive the slope of the edge by taking the perpendicular of the slope of the connecting line segment, and then performing a line intersection test on the two edges, but there is an even easier way.

As long as three sites are not collinear, then the edges that perpendicularly bisect the segments between the sites are also tangent to the circle whose edge contains all three sites. Therefore the breakpoints defined by a triple of Voronoi sites converge if the center of the circle defined by the three sites is in front of the middle site, where "in front" and "behind" depend on the coordinate system and sweepline alignment you have chosen.

In my case, I have a horizontal sweepline that I am moving from minimum y to maximum y, so the breakpoints converge if the y-coordinate of the center of the circle is greater than the y-coordinate of the middle site, and diverge otherwise.

Edit: Kristian D'Amato rightfully points out that the algorithm above misses some convergence cases. The final algorithm I ended up using is below. Of course, I'm not confident that its 100% correct, but it seems to work for all the cases I tried it out on.

Given left, middle, right sites
    if they are collinear, return false
    center = ComputeCircleCenterDefinedBy3Points(left, middle, right)
    return IsRightOfLine(left, middle, center) && IsRightOfLine(middle, right, center)

IsRightOfLine(start, end, point)
    ((end.X - start.X) * (point.Y - start.Y) - (end.Y - start.Y) * (point.X - start.X)) <= 0

Answer 3

If the sites are ordered clockwise around the center of the circle, the arc is converging. If they are ordered counterclockwise around the center of the circle, the arc is diverging. (or vice versa, depending on your implementation). Testing for cw or ccw falls out of the code you use to find the center of the circle.

Here's a snippet of C# code for computing the circumcenter d of points a,b,c:

        Vector2 ba = b - a;
        Vector2 ca = c - a;     
        float baLength = (ba.x * ba.x) + (ba.y * ba.y);
        float caLength = (ca.x * ca.x) + (ca.y * ca.y); 
        float denominator = 2f * (ba.x * ca.y - ba.y * ca.x);   
        if (denominator <= 0f ) { // Equals 0 for colinear points.  Less than zero if points are ccw and arc is diverging.
            return false;  // Don't use this circle event!
        };
        d.x = a.x + (ca.y * baLength - ba.y * caLength) / denominator ;
        d.y = a.y + (ba.x * caLength - ca.x * baLength) / denominator ;