What does a Voronoi diagram have to do with a flow network?

[simontatham@hachyderm.io]

What does a Voronoi diagram have to do with a flow network?

Some years ago I was trying to make 3D-printable models of Britain for a friend. (Both the whole of Great Britain as a single model, and also zooming in to model individual mountains.) Our source material was the Ordnance Survey OpenData, which provides 10m contour lines plus an outline for the coast itself.

The data was in the form of XML files, each covering a 10km × 10km square of the OS national grid, describing segments of whichever outlines intersect that square, each as a sequence of (x,y) points. (No clear specification of how the points are supposed to be connected; for simplicity I just assumed a straight line between each pair of points.)

The coast outline, unfortunately, was not a closed curve. Partly that was because the segments in adjacent grid-square files didn't always match up exactly – apparently nobody had checked the endpoints of line segments were consistent across grid-square boundaries. But also, I think they'd deliberately decided not to bother with rivers: the coast outline would just stop, on one side of a river mouth, and resume on the other. That saved them having to decide how much of the river itself to include, but also left a gap in the outline.

I needed a properly closed loop, in order to determine what was inside it and what was outside it. So I had to close all those gaps. But the data was too huge and detailed to go round the coast by hand with Inkscape or the like. It'd have taken forever. It had to be automated.

So, the problem: given a _mostly_ complete outline of a polygon, fill in the gaps, in some way that looks minimal, such as 'smallest total length of the line segments you add'.

The idea I hit on was to combine two algorithms I already knew about, but hadn't ever considered using together.

(1/5)

[simontatham@hachyderm.io]

Given a set S of points in the plane, the _Voronoi cell_ of a point P ∈ S is the region of the plane which is closer to P than to any other point of S. A Voronoi cell is a polygon; where two cells touch, say for points P and Q, the edge between them is the perpendicular bisector of the line PQ. If the point set S is finite, then on the edges of the diagram the Voronoi cells extend off to infinity.

The Voronoi diagram is what you get by dividing the whole plane into those Voronoi cells. Regarded as a graph in its own right, the edges of the diagram are all segments of perpendicular bisectors as I just described, and the vertices are points equidistant from three or more of the points of S.

A Voronoi diagram for n points can be built in O(n log n) time, using https://en.wikipedia.org/wiki/Fortune%27s_algorithm, which sweeps a line across the plane and maintains a sorted list of the partially built Voronoi cells along the current location of the line.

(2/5)

[simontatham@hachyderm.io]

A flow network is a graph, with each edge imagined as a pipe through which some kind of fluid can flow. Edges are annotated with _capacities_ – the maximum rate at which fluid can flow through that pipe. Some vertices are _sources_ (fluid is imagined to enter the graph there) and some are _sinks_ (fluid arriving there leaves the graph).

The most standard problem involving a flow network is: what's the maximum rate at which this set of pipes can move fluid from the source(s) to the sink(s)? And what does that optimal flow look like, in the sense of how much fluid goes along each individual pipe (and in which direction)?

A key theorem is the max-flow min-cut theorem. If you imagine partitioning the graph into two components, with all the sources on one side and all the sinks on the other, then the total capacity of the pipes crossing that divide is an upper bound on the maximum flow of the network. If you choose the partition for which the total capacity is smallest - the minimum cut - then that gives you the strongest bound on the maximum flow. The theorem says that in fact these quantities are equal: if the flow from sources to sinks can't exceed X, it's _because_ there is some cut across the graph with maximum total capacity X. A corollary is that in the best possible flow, the pipes across the minimum cut must all be operating at maximum capacity.

The https://en.wikipedia.org/wiki/Ford%E2%80%93Fulkerson_algorithm finds the maximum flow in such a network. But more: it _also_ finds the minimum cut, as a by-product of the algorithm. That is, when it terminates and tells you the flow can't be increased any more, it also tells you _why_ it can't be improved, by showing you the set of pipes that are the limiting factor and already at maximum capacity.

(3/5)

[simontatham@hachyderm.io]

My problem was to fill in the gaps in an almost-complete polygon. I had the idea of imagining that you tried to fill the polygon with fluid. It would leak out of the holes in the outline – out of each hole proportionally to its size. So if we can find the maximum flow from the interior of the polygon to the exterior, that tells us the total size of the holes. And the minimum cut is precisely the smallest amount of extra barrier you'd have to add to block up that flow.

The actual "size of a hole" for this purpose is the length of the line segment you'd have to add to block it, say between points P and Q. The flow of fluid is _across_ that line segment, say perpendicular to it. So I imagined the Voronoi diagram of the set of polygon vertices, with the edge separating P and Q's Voronoi cells being a pipe, with capacity equal to the distance PQ.

So my actual algorithm was:
1. Run Fortune's algorithm on the points of my partial outline, to build its Voronoi diagram.
2. Make a flow network out of the edges and vertices of the Voronoi cells. If two points P,Q are connected by the partial outline, the capacity of the edge between them is zero. Otherwise it equals the distance PQ.
3. Choose a pair of points by hand that you want to be inside and outside the final polygon. (For outside, it's convenient to choose the 'point at infinity', where all the edges of the unbounded Voronoi cells are imagined to meet.)
4. Run the Ford-Fulkerson algorithm to solve the max-flow/min-cut problem for that flow network.

The minimum cut output by the algorithm tells you which pairs of outline points to connect with an extra edge!

(4/5)

[simontatham@hachyderm.io]

I was really pleased with this daft idea. And it worked! I got a successful run of it on the OS OpenData's incomplete outline of the coastline of Great Britain.

With hindsight, there are several problems with it. What if an edge of the input near-polygon connects two polygon vertices whose Voronoi cells aren't adjacent? Similarly, what if the shortest line segment that would block up a gap is between two points in non-adjacent cells?

If I were trying to polish up this idea for general use, I'd probably want to do something more sophisticated. In principle, you can make a Voronoi diagram for things other than points: you could imagine a set of _line segments_ in the plane having Voronoi cells in the same sense. So I could have made the Voronoi diagram of the _edges_ of the provided partial outline, which would have at least solved problem 1: any two edges that connect at a vertex certainly would have adjacent Voronoi cells.

But that's much harder: it's not such a standard well-studied problem, and it's inherently more tricky, because the Voronoi cells of line segments are not polygons any more (their edges are curved).

So since I was doing this as a one-off, I did it the simple way and just hoped that this kind of problem wouldn't arise. (A coastline mostly doesn't double back on itself very sharply, so generally, adjacent points in the outline _can_ be expected to be in touching Voronoi cells.)

And it didn't – phew!

(5/5)

[darkling@mstdn.social]

@simontatham I love it. That's a really nice algorithm.

[klausman@mas.to]

@simontatham You can't post this without also providing photos of the 3d prints! :D

[simontatham@hachyderm.io]

@klausman alas, *this* particular project of filling in the OS OpenData outline of Great Britain was for a print that I never finished. Thinking up this weird algorithm was the most fun I got out of it.

But just so you're not 100% disappointed, here's a picture of the GB model I *did* successfully print a few years earlier, based on a much lower-resolution data set (GTOPO30).

[klausman@mas.to]

@simontatham Thank you!