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.

Permutation

Did you solve the delauny triangulation? – Soso

No, this is my first attempt at computational geometry. I know they're duals of each other, but I'd like to implement the straightforward Fortune's algorithm first. – Permutation

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.

Renner

Very interesting. More of a logical than mathematical approach, but thinking it over it seems like this should work. I think that this part of the algorithm is a prime candidate for unit testing. – Permutation

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

Permutation

I too was a bit annoyed that nobody seemed to offer a quick and easy way of checking for convergence. Your question (and answer) was a great find. But I think there are certain cases that converge and which the algorithm as you described it does not capture. These are the cases in which the breakpoints actually travel backwards (in the sense you defined it; opposite to the sweepline movement). This can happen when three sites are arranged, going from left to right, in the opposite sense in which their parabola arcs appear; the breakpoints will still converge in the future, below the beachline. – Clein

@KristianD'Amato Thanks for reminding me of my answer here, I've updated at the bottom of the post with the improved algorithm I think is correct. – Permutation

Based on KristianD'Amato's full answer I should clarify that my IsRightOfLine function is specific to the coordinate system and sweepline orientation (specifically, sweepline goes from greater y to lesser y) I chose, and if you adapt the function for use in your own implementation you should check to make sure it works for your coordinate setup. – Permutation

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 ;

Illmannered

This is the only right answer. All others gives only too roundabout solutions. Having the determinant value we also can detect the case of collinear points, which leads to parallel edges. – Damn

@Orient I'm confused about the intuition behind this, possibly because I haven't fully looked into the intuition behind the determinant method for computing circumcircles. Any chance you could explain? In particular, I'm not fully sure of what clockwise vs. counterclockwise is supposed to mean, and have some examples that I'm not sure how to classify. – Howitzer

@Howitzer Me? It is Brian Upton's answer. – Damn

@Orient Alright. The reason I pinged you was the original author isn't active and you seemed to indicate you understood the answer in your comment. – Howitzer

@Howitzer I don't think I can explain better, then, say, wiki. Here is my implementation of the Fortune's algorithm. – Damn

@Orient No worries. Thanks for including the link to your work! – Howitzer