de Fraysseix-Pach-Pollack algorithm in Python
Linear-time straight-line drawing on a (2n-4) by (n-2) grid using the shift method.
- File
de_fraysseix_pach_pollack.py, 173 lines- Level
- I0 not checked yet
- Toolchain
- Python 3.14, pytest 9
- Test framework
- pytest
- Run the tests
python3 -m pytest
# SPDX-License-Identifier: MIT
# Part of https://algorithm.tamnd.com/de-fraysseix-pach-pollack/. The authors are listed in prov/FQHSKN.jsonl.
"""Straight-line grid drawing of a maximal planar graph by the shift method.
The embedding is given as a rotation system: rot[v] lists the neighbours of v
in counterclockwise order. The outer face is the triangle (v1, v2, vn), listed
counterclockwise. The result puts v1 at (0, 0), v2 at (2n - 4, 0) and vn at
(n - 2, n - 2), with integer coordinates and no crossings.
"""
def canonical_ordering(rot, v1, v2, vn):
"""Return (order, left, right) for a maximal planar graph.
order is a canonical ordering v1, v2, ..., vn. For k >= 4, left[v] and
right[v] are the first and last neighbours of v on the contour of the
graph induced by the vertices before v.
"""
n = len(rot)
prev = [-1] * n
nxt = [-1] * n
on = [False] * n
chords = [0] * n
left = [-1] * n
right = [-1] * n
nxt[v1], prev[vn], nxt[vn], prev[v2] = vn, v1, v2, vn
on[v1] = on[v2] = on[vn] = True
removed_order = []
stack = []
v = vn
for k in range(n, 3, -1):
if k < n:
while True:
v = stack.pop()
if on[v] and chords[v] == 0 and v != v1 and v != v2:
break
a, b = prev[v], nxt[v]
left[v], right[v] = a, b
removed_order.append(v)
on[v] = False
r = rot[v]
d = len(r)
i = r.index(a)
new = []
while True:
i = (i + 1) % d
w = r[i]
if w == b:
break
new.append(w)
if not new:
# the edge ab was a chord and is now a contour edge
chords[a] -= 1
chords[b] -= 1
nxt[a], prev[b] = b, a
for x in (a, b):
if chords[x] == 0:
stack.append(x)
continue
chain = [a] + new + [b]
for j in range(1, len(chain)):
nxt[chain[j - 1]] = chain[j]
prev[chain[j]] = chain[j - 1]
for w in new:
# mark the new vertices one at a time, so that a chord
# between two new vertices is counted once
on[w] = True
for y in rot[w]:
if on[y] and y != w and y != prev[w] and y != nxt[w]:
chords[w] += 1
chords[y] += 1
for x in [a] + new + [b]:
if chords[x] == 0:
stack.append(x)
third = nxt[v1]
order = [v1, v2, third] + removed_order[::-1]
return order, left, right
def check_triangulation(rot, v1, v2, vn):
"""Raise ValueError unless rot is a plane triangulation with outer face
v1, v2, vn in counterclockwise order."""
n = len(rot)
if n < 3:
raise ValueError("need at least 3 vertices")
where = []
m2 = 0
for v, r in enumerate(rot):
if len(set(r)) != len(r) or v in r:
raise ValueError("not a simple graph")
where.append({w: i for i, w in enumerate(r)})
m2 += len(r)
for v, r in enumerate(rot):
for w in r:
if not 0 <= w < n or v not in where[w]:
raise ValueError("adjacency is not symmetric")
seen = set()
faces = 0
for v in range(n):
for w in rot[v]:
if (v, w) in seen:
continue
faces += 1
length = 0
a, b = v, w
while (a, b) not in seen:
seen.add((a, b))
length += 1
r = rot[b]
a, b = b, r[(where[b][a] - 1) % len(r)]
if length != 3:
raise ValueError("a face is not a triangle")
if n - m2 // 2 + faces != 2:
raise ValueError("the rotation system is not planar")
for a, b in ((v1, v2), (v2, vn), (vn, v1)):
if a not in where[b]:
raise ValueError("outer face is not a triangle")
for a, b, c in ((v1, vn, v2), (v2, v1, vn), (vn, v2, v1)):
r = rot[a]
if r[(where[a][b] + 1) % len(r)] != c:
raise ValueError("outer face is not v1, v2, vn counterclockwise")
def draw(rot, v1, v2, vn):
"""Return integer coordinates [(x, y)] of a crossing-free drawing."""
check_triangulation(rot, v1, v2, vn)
n = len(rot)
if n == 3:
pos = [None] * 3
pos[v1], pos[v2], pos[vn] = (0, 0), (2, 0), (1, 1)
return pos
order, lft, rgt = canonical_ordering(rot, v1, v2, vn) # @line L1
dx = [0] * n # x offset relative to the parent in the shift tree
y = [0] * n
nxt = [-1] * n # next vertex on the contour, or in a covered chain
under = [-1] * n # first vertex covered by v
v3 = order[2]
dx[v1], dx[v3], dx[v2] = 0, 1, 1
y[v3] = 1
nxt[v1], nxt[v3] = v3, v2 # @line L2
for v in order[3:]: # @line L3
wp, wq = lft[v], rgt[v] # @line L4
first = nxt[wp]
dx[first] += 1 # @line L5
dx[wq] += 1 # @line L6
delta = 0
w = first
while True:
delta += dx[w]
if w == wq:
break
w = nxt[w]
dx[v] = (delta - y[wp] + y[wq]) // 2 # @line L7
y[v] = (delta + y[wp] + y[wq]) // 2
dx[wq] = delta - dx[v] # @line L8
if first != wq:
dx[first] -= dx[v] # @line L9
w = first
while nxt[w] != wq:
w = nxt[w]
nxt[w] = -1
under[v] = first
nxt[wp], nxt[v] = v, wq # @line L10
x = [0] * n # @line L11
stack = [(v1, 0)]
while stack:
v, base = stack.pop()
x[v] = base + dx[v]
if nxt[v] != -1:
stack.append((nxt[v], x[v]))
if under[v] != -1:
stack.append((under[v], x[v]))
return [(x[v], y[v]) for v in range(n)] # @line L12
test_de_fraysseix_pach_pollack.py 266 lines
# SPDX-License-Identifier: MIT
# Part of https://algorithm.tamnd.com/de-fraysseix-pach-pollack/. The authors are listed in prov/FQHSKN.jsonl.
import itertools
import random
import pytest
from de_fraysseix_pach_pollack import canonical_ordering, draw
def rot_from_faces(n, inner):
"""Rotation system from the counterclockwise inner triangles.
The outer face is 0, 1, 2 (counterclockwise), so the outer boundary
walked with the interior on the right is 0, 2, 1.
"""
succ = [dict() for _ in range(n)]
for f in list(inner) + [(0, 2, 1)]:
for t in range(3):
x, y, z = f[t - 1], f[t], f[(t + 1) % 3]
succ[y][z] = x
rot = []
for v in range(n):
start = min(succ[v])
r = [start]
w = succ[v][start]
while w != start:
r.append(w)
w = succ[v][w]
rot.append(r)
return rot
def random_triangulation(n, rng, flips=0):
"""Rotation system of a random maximal planar graph with outer face 0, 1, 2.
Vertices are inserted into random inner faces, then random edge flips are
applied. Faces are kept as counterclockwise triangles.
"""
faces = [(0, 1, 2)]
for v in range(3, n):
i = rng.randrange(len(faces))
a, b, c = faces[i]
faces[i] = (a, b, v)
faces += [(b, c, v), (c, a, v)]
outer_edges = {frozenset((0, 1)), frozenset((1, 2)), frozenset((0, 2))}
for _ in range(flips):
dart_face = {}
for k, f in enumerate(faces):
for t in range(3):
dart_face[(f[t], f[(t + 1) % 3])] = k
k = rng.randrange(len(faces))
f = faces[k]
t = rng.randrange(3)
a, b, c = f[t], f[(t + 1) % 3], f[(t + 2) % 3]
if frozenset((a, b)) in outer_edges:
continue
g = faces[dart_face[(b, a)]]
d = [x for x in g if x not in (a, b)][0]
if (c, d) in dart_face or (d, c) in dart_face:
continue
faces[k] = (a, d, c)
faces[dart_face[(b, a)]] = (b, c, d)
return rot_from_faces(n, faces)
def oracle_drawing(rot, order):
"""The original O(n^2) shift method with explicit sets M(w)."""
n = len(rot)
v1, v2, v3 = order[0], order[1], order[2]
x = [0] * n
y = [0] * n
x[v2], x[v3], y[v3] = 2, 1, 1
contour = [v1, v3, v2]
move = {v1: {v1, v2, v3}, v3: {v3, v2}, v2: {v2}}
placed = {v1, v2, v3}
for v in order[3:]:
idx = [i for i, w in enumerate(contour) if w in rot[v]]
p, q = idx[0], idx[-1]
assert idx == list(range(p, q + 1)), "neighbours not contiguous"
assert p < q
for w in move[contour[p + 1]]:
x[w] += 1
for w in move[contour[q]]:
x[w] += 1
wp, wq = contour[p], contour[q]
x[v] = (x[wp] + x[wq] + y[wq] - y[wp]) // 2
y[v] = (x[wq] - x[wp] + y[wp] + y[wq]) // 2
for i in range(p + 1):
move[contour[i]].add(v)
move[v] = {v} | move[contour[p + 1]]
contour = contour[: p + 1] + [v] + contour[q:]
placed.add(v)
return [(x[v], y[v]) for v in range(n)]
def orient(a, b, c):
return (b[0] - a[0]) * (c[1] - a[1]) - (b[1] - a[1]) * (c[0] - a[0])
def on_segment(a, b, c):
return (min(a[0], b[0]) <= c[0] <= max(a[0], b[0])
and min(a[1], b[1]) <= c[1] <= max(a[1], b[1]))
def segments_meet(a, b, c, d):
d1, d2 = orient(c, d, a), orient(c, d, b)
d3, d4 = orient(a, b, c), orient(a, b, d)
if ((d1 > 0 > d2) or (d1 < 0 < d2)) and ((d3 > 0 > d4) or (d3 < 0 < d4)):
return True
if d1 == 0 and on_segment(c, d, a):
return True
if d2 == 0 and on_segment(c, d, b):
return True
if d3 == 0 and on_segment(a, b, c):
return True
if d4 == 0 and on_segment(a, b, d):
return True
return False
def check_plane(rot, pos):
n = len(rot)
assert len(set(pos)) == n, "two vertices share a point"
edges = [(u, v) for u in range(n) for v in rot[u] if u < v]
for i, (a, b) in enumerate(edges):
for c in range(n):
if c != a and c != b and orient(pos[a], pos[b], pos[c]) == 0 \
and on_segment(pos[a], pos[b], pos[c]):
raise AssertionError("vertex on an edge")
for (c, d) in edges[i + 1:]:
if len({a, b, c, d}) < 4:
continue
assert not segments_meet(pos[a], pos[b], pos[c], pos[d]), "crossing"
def check_rotation(rot, pos):
"""The drawing realizes the given rotation system (same cyclic orders)."""
import math
for v, r in enumerate(rot):
got = sorted(r, key=lambda w: math.atan2(pos[w][1] - pos[v][1],
pos[w][0] - pos[v][0]))
i = got.index(r[0])
assert got[i:] + got[:i] == r
def check_grid(pos, v1, v2, vn):
n = len(pos)
w, h = 2 * n - 4, n - 2
if n == 3:
w, h = 2, 1
for (x, y) in pos:
assert 0 <= x <= w and 0 <= y <= h
assert pos[v1] == (0, 0) and pos[v2] == (w, 0)
def test_triangle():
rot = [[1, 2], [2, 0], [0, 1]]
pos = draw(rot, 0, 1, 2)
assert pos == [(0, 0), (2, 0), (1, 1)]
def test_k4():
# triangle 0, 1, 2 counterclockwise with 3 inside
rot = [[1, 3, 2], [2, 3, 0], [0, 3, 1], [0, 1, 2]]
pos = draw(rot, 0, 1, 2)
assert pos == [(0, 0), (4, 0), (2, 2), (2, 1)]
check_plane(rot, pos)
check_rotation(rot, pos)
def test_octahedron():
import math
pts = [(0, 0), (10, 0), (5, 9), (5, 2), (6.5, 4), (3.5, 4)]
adj = [[] for _ in range(6)]
for i in range(6):
for j in range(i + 1, 6):
if (i, j) not in [(0, 4), (1, 5), (2, 3)]:
adj[i].append(j)
adj[j].append(i)
rot = [sorted(adj[u], key=lambda w: math.atan2(pts[w][1] - pts[u][1],
pts[w][0] - pts[u][0]))
for u in range(6)]
order, _, _ = canonical_ordering(rot, 0, 1, 2)
assert order == [0, 1, 3, 5, 4, 2]
pos = draw(rot, 0, 1, 2)
assert pos == [(0, 0), (8, 0), (4, 4), (5, 1), (4, 3), (3, 2)]
assert pos == oracle_drawing(rot, order)
check_plane(rot, pos)
check_rotation(rot, pos)
def test_random_against_oracle():
rng = random.Random(2024)
for trial in range(600):
n = rng.randint(3, 30)
rot = random_triangulation(n, rng, flips=rng.randint(0, 3 * n))
pos = draw(rot, 0, 1, 2)
check_grid(pos, 0, 1, 2)
if n > 3:
order, _, _ = canonical_ordering(rot, 0, 1, 2)
assert sorted(order) == list(range(n))
assert order[:2] == [0, 1] and order[-1] == 2
assert pos == oracle_drawing(rot, order)
if n <= 18:
check_plane(rot, pos)
check_rotation(rot, pos)
def test_large():
rng = random.Random(7)
n = 20000
rot = random_triangulation(n, rng)
pos = draw(rot, 0, 1, 2)
check_grid(pos, 0, 1, 2)
assert pos[2] == (n - 2, n - 2)
def test_deep_nesting():
# vertex v is inserted into the triangle (0, 1, v - 1): a long chain of
# nested triangles, which makes the contour and the shift tree long
n = 3000
faces = [(0, 1, 2)]
for v in range(3, n):
faces[0] = (0, 1, v)
faces += [(1, v - 1, v), (v - 1, 0, v)]
rot = rot_from_faces(n, faces)
pos = draw(rot, 0, 1, 2)
check_grid(pos, 0, 1, 2)
check_rotation(rot, pos)
def test_rejects_non_triangulations():
# a single vertex, a single edge, a path (a tree)
for rot in ([[]], [[1], [0]], [[1], [0, 2], [1]]):
with pytest.raises(ValueError):
draw(rot, 0, 1, len(rot) - 1)
# a star with three leaves: a tree, faces are not triangles
with pytest.raises(ValueError):
draw([[1, 2, 3], [0], [0], [0]], 1, 2, 3)
# K4 with the outer face listed clockwise
with pytest.raises(ValueError):
draw([[1, 3, 2], [2, 3, 0], [0, 3, 1], [0, 1, 2]], 0, 2, 1)
def test_k5_and_k33_have_no_plane_rotation():
# every rotation system of K5 fails, and so does every one of K3,3
def all_rotations(adj):
choices = []
for v, nb in enumerate(adj):
first, rest = nb[0], nb[1:]
choices.append([[first] + list(p) for p in itertools.permutations(rest)])
return itertools.product(*choices)
k5 = [[w for w in range(5) if w != v] for v in range(5)]
for rot in all_rotations(k5):
with pytest.raises(ValueError):
draw([list(r) for r in rot], 0, 1, 2)
k33 = [[3, 4, 5]] * 3 + [[0, 1, 2]] * 3
count = 0
for rot in all_rotations(k33):
count += 1
with pytest.raises(ValueError):
draw([list(r) for r in rot], 0, 3, 1)
assert count == 2 ** 6
What the next level needs
- the tests did not run in this build
Read the full page about de Fraysseix-Pach-Pollack algorithm