Planar Max Flow, in depth: dual shortest paths, recovering the flow from potentials, and beyond st-planar graphs

Cuts are dual paths

Fix a planar embedding. Its faces are the regions the drawing cuts the plane into, including the unbounded outer face. The dual graph has one node per face and one dual edge per primal edge, joining the two faces on either side; give each dual edge the capacity of the primal edge it crosses.

Assume s and t are both on the outer face. Draw an imaginary edge from t back to s through the outer face. That splits the outer face in two: call the part above the s-to-t boundary face A and the part below it face B. Any closed curve separating s from t must go from A to B through the drawing, and the edges it crosses form an s-t cut. Conversely, every minimal cut corresponds to such a curve. So the minimum cut capacity equals the shortest-path distance from A to B in the dual, and by max-flow min-cut that is the max-flow value.

Historically this is the oldest flow algorithm there is. Ford and Fulkerson's 1956 paper handled the st-planar case by repeatedly augmenting along the uppermost path; Itai and Shiloach made that O(n log n) in 1979, and Hassin observed in 1981 that dual shortest paths give both the value and the flow. With the linear-time planar shortest-path algorithm of Henzinger, Klein, Rao and Subramanian, the st-planar case is O(n).

From distances to a flow

Run Dijkstra in the dual from A and let d(f) be the distance to face f. For a primal edge traversed from u to v, let L be the face on its left and R the face on its right. Set the flow from u to v to d(R) - d(L). Three facts make this a maximum flow.

Capacity. L and R are joined by a dual edge of length c(e), so by the triangle inequality |d(R) - d(L)| ≤ c(e). The flow fits in either direction.

Conservation. Walk around an interior vertex v. Its edges, in rotation order, are separated by faces, and the right face of one outgoing edge is the left face of the next. The sum of d(R) - d(L) over v's edges telescopes to zero. The only vertices where it does not are s and t, because there the outer face is split into A and B by the imaginary edge.

Value. At s the telescoping sum leaves d(B) - d(A) = d(B), which is exactly the minimum cut. A feasible flow whose value equals a cut capacity is maximum.

Nothing iterates. One shortest-path computation produces the potentials, and each edge's flow is a subtraction. Integer capacities give integer distances and therefore an integral flow.

Implementation

The implementation below handles undirected graphs with integer or real capacities, a straight-line drawing given by coordinates, and s and t on the outer face. It sorts each vertex's neighbours by angle to obtain the rotation system, traces faces, finds the outer face as the one with negative signed area, splits it at s and t, and runs Dijkstra.

import heapq, math
from collections import defaultdict

def faces_of(pos, edges):
    nbrs = defaultdict(list)
    for u, v, _ in edges:
        nbrs[u].append(v); nbrs[v].append(u)
    for u in nbrs:                                   # counter-clockwise rotation
        nbrs[u].sort(key=lambda w: math.atan2(pos[w][1] - pos[u][1], pos[w][0] - pos[u][0]))
    face_of, faces = {}, []
    for u in nbrs:
        for v in nbrs[u]:
            if (u, v) in face_of: continue
            fid, walk, d = len(faces), [], (u, v)
            while d not in face_of:                  # face lies left of each dart
                face_of[d] = fid; walk.append(d)
                a, b = d; ring = nbrs[b]
                d = (b, ring[(ring.index(a) - 1) % len(ring)])
            faces.append(walk)
    return face_of, faces

def planar_max_flow(pos, edges, s, t):
    face_of, faces = faces_of(pos, edges)
    area = lambda w: sum(pos[a][0]*pos[b][1] - pos[b][0]*pos[a][1] for a, b in w)
    walk = faces[min(range(len(faces)), key=lambda f: area(faces[f]))]
    k = next(i for i, (a, _) in enumerate(walk) if a == s)
    A, B = len(faces), len(faces) + 1                # split the outer face at s and t
    side, cur = {}, A
    for d in walk[k:] + walk[:k]:
        if d[0] == t: cur = B
        side[d] = cur
    face = lambda d: side.get(d, face_of[d])
    adj = defaultdict(list)
    for u, v, cap in edges:
        f, g = face((u, v)), face((v, u))
        adj[f].append((g, cap)); adj[g].append((f, cap))
    dist, pq = {A: 0}, [(0, A)]
    while pq:
        d0, f = heapq.heappop(pq)
        if d0 > dist[f]: continue
        for g, cap in adj[f]:
            if d0 + cap < dist.get(g, math.inf):
                dist[g] = d0 + cap; heapq.heappush(pq, (d0 + cap, g))
    flow = {(u, v): dist[face((v, u))] - dist[face((u, v))] for u, v, _ in edges}
    return dist[B], flow

I compared this with a plain BFS augmenting-path solver on two thousand random layered grids, up to six by six with random diagonals, with s and t fanned out to the left and right columns. Every value matched, every edge flow stayed within capacity and every interior vertex balanced.

Worked example and measurement

Take a 3 by 2 grid with s joined to the left column and t to the right column, and integer capacities between 2 and 9. The function returns a value of 10. The figure shows the per-edge flows it computes; red edges are saturated.

Check conservation at a couple of vertices. Vertex (0,0) receives 7 from s and sends 2 right and 5 up: 7 in, 7 out. Vertex (1,1) receives 8 from the left and sends 4 right and 4 down to (1,0). Vertex (2,0) receives 6 from the left and 2 from above, and sends 8 to t. Note the edge (1,0)-(1,1): the code reported -4 on the dart from (1,0) to (1,1), meaning 4 units flow downward. Undirected edges carry signed flow, and the sign comes free from the potential difference.

On the same machine, a 200 by 200 grid (40,002 vertices and 91,936 edges, value 810) took 0.83 seconds with this code and 44.9 seconds with the augmenting-path solver, both in pure CPython. A tuned Dinic or push-relabel implementation would close much of that gap, but not the asymptotic one: here the work is one heap-based shortest path, O(n log n).

Recovering the cut, and reusing the dual

Applications often need the cut, not the flow: which pixels are foreground, or which links to reinforce. Dijkstra already found it: follow the predecessor pointers from B back to A, and the primal edges crossed by that dual path are a minimum cut, since their capacities sum to d(B). To name the vertex sides, take the source side as every vertex reachable from s through edges that are not saturated in the forward direction, exactly as in the residual-graph argument for general flows. In the example that search from s reaches every vertex except t, because the only saturated edges that block it are the two sink edges.

When many cuts are needed on one graph with different terminals on the outer face, keep the face structure and dual adjacency, and rerun only Dijkstra. Face tracing is the expensive, error-prone part; shortest paths are cheap. With small integer capacities, a bucket queue (Dial's algorithm) replaces the heap and makes each run linear in the number of edges plus the maximum distance.

Beyond the st-planar case

For directed st-planar graphs the same construction works with a small change: give each dual edge length c(e) when crossed in the direction that counts toward the cut and 0 in the other. The triangle inequalities then give 0 ≤ flow ≤ c(e) instead of the symmetric bound.

The super source trick in the example stays planar only because s connects to vertices on one boundary side. Joining a source to arbitrary interior pixels, as segmentation with scattered seeds does, usually destroys planarity. So does an 8-connected pixel grid, whose diagonals cross. That is why vision libraries mostly ship general solvers such as Boykov-Kolmogorov, and why planar graph-cut methods restrict the model to 4-connectivity with boundary terminals.

Failure modes

• Not actually planar. Angle-sorted rotations from coordinates are only a valid embedding if the straight-line drawing has no crossings. Test with a planarity test and compute a combinatorial embedding rather than trusting coordinates from a real-world map, where overpasses break planarity.
• s or t not on the outer face. The A/B split then does not exist. The listing fails loudly but unhelpfully: StopIteration if s is missing from the outer walk, KeyError on dist[B] if t is. Assert both explicitly with a clear message; if they share an inner face, re-embed so that face becomes the outer one.
• Degenerate geometry. Two neighbours at the same angle (collinear overlapping edges) or parallel edges give an ambiguous rotation. Merge parallel edges by summing capacity.
• Disconnected graphs. Faces of separate components nest and face tracing no longer yields one dual graph. Drop components without s or t first.
• Floating-point capacities. Flows are differences of large distances, so tiny rounding errors appear as slightly negative residuals and conservation that is off in the last digits. Compare with a tolerance, or scale capacities to integers.

Trade-offs

Use the dual method when the graph is genuinely planar, s and t share a face, and you solve many instances or very large ones: it is simple, has no iteration-count risk, and gives the cut and the flow together. Use a general solver when planarity is approximate, terminals are interior, or capacities change incrementally and you want warm starts, which augmenting-path and push-relabel methods support and a one-shot shortest path does not. Memory is modest too: by Euler's formula a connected planar graph with n vertices and m edges has m - n + 2 faces, so the dual has exactly m edges and at most 2n - 4 nodes, since a simple planar graph has at most 3n - 6 edges. The arbitrary-terminal planar algorithms in the table are mostly of theoretical interest; implementations are rare and constants are high.

What to do next

1. Run the code above on a small grid and verify conservation and capacity yourself.
2. Add an assertion that s and t are on the outer face, and a planarity check before trusting a coordinate embedding.
3. Swap Dijkstra for a bucket queue when capacities are small integers.
4. Extend it to directed edges with the c(e)/0 dual lengths and test it against a general solver.
5. Measure it against your current max-flow library on your real instances before switching.
6. Read the arbitrary-terminal results in the table only if your terminals are interior.