diff --git a/colored_pentagon_reduction/example.py b/colored_pentagon_reduction/example.py
index 8ac83d1..cc0b1d1 100644
--- a/colored_pentagon_reduction/example.py
+++ b/colored_pentagon_reduction/example.py
@@ -5,10 +5,9 @@ import json
 from collections import defaultdict
 from pathlib import Path
 from typing import Any, cast, TypedDict, Literal
-import math
-
 from sage.all import graphs, Graph, save, load  # type: ignore[attr-defined]  # pylint: disable=no-name-in-module
 from lib.tutte_embedding import tutte_embedding
+from lib.planar_embedding import outer_face, get_embedding_from_pos
 
 DIR = Path(__file__).parent.parent
 PALETTE = ['red', 'blue', 'green', 'yellow']
@@ -86,50 +85,6 @@ def save_colored_graph(g: Graph, coloring: VertexColoring) -> tuple[Graph, Verte
     save(g_canon, str(out_dir / 'graph'))
     return g_canon, canonical_coloring, cid
 
-def outer_face(g: Graph) -> list[Any]:
-    """Return the vertices of the outer (unbounded) face of g using its planar embedding and positions."""
-    pos = g.get_pos()
-    embedding = g.get_embedding()
-    if not pos or not embedding:
-        raise Exception("Position and embedding required to find outer face")
-
-    visited: set[tuple[Any, Any]] = set()
-    faces: list[list[Any]] = []
-
-    for u in g.vertices():
-        for v in embedding[u]:
-            if (u, v) not in visited:
-                face: list[Any] = []
-                cu, cv = u, v
-                while (cu, cv) not in visited:
-                    visited.add((cu, cv))
-                    face.append(cu)
-                    neighbors = embedding[cv]
-                    cw = neighbors[(neighbors.index(cu) + 1) % len(neighbors)]
-                    cu, cv = cv, cw
-                faces.append(face)
-
-    def signed_area(face: list[Any]) -> float:
-        coords = [pos[v] for v in face]
-        n = len(coords)
-        return sum(
-            coords[i][0] * coords[(i + 1) % n][1] - coords[(i + 1) % n][0] * coords[i][1]
-            for i in range(n)
-        ) / 2
-
-    return min(faces, key=signed_area)
-
-def get_embedding_from_pos(g: Graph) -> dict[Any, list[Any]]:
-    """Compute a combinatorial planar embedding by sorting each vertex's neighbors by angle."""
-    pos_or_none = g.get_pos()
-    if pos_or_none is None:
-        raise Exception("Cannot get embedding without position")
-    pos: dict[Any, Any] = pos_or_none
-    return {
-        v: sorted(g.neighbors(v), key=lambda w: math.atan2(pos[w][1] - pos[v][1], pos[w][0] - pos[v][0]))
-        for v in g.vertices()
-    }
-
 def _neighbors_form_cycle(g: Graph, v: Any) -> bool:
     """Return True if the neighbors of v induce a cycle in g."""
     return bool(cast(Graph, g.subgraph(g.neighbors(v))).is_cycle())
diff --git a/lib/planar_embedding.py b/lib/planar_embedding.py
new file mode 100644
index 0000000..3d39e90
--- /dev/null
+++ b/lib/planar_embedding.py
@@ -0,0 +1,51 @@
+"""Utilities for working with planar embeddings of graphs."""
+import math
+from typing import Any
+
+from sage.all import Graph  # type: ignore
+
+
+def outer_face(g: Graph) -> list[Any]:
+    """Return the vertices of the outer (unbounded) face of g using its planar embedding and positions."""
+    pos = g.get_pos()
+    embedding = g.get_embedding()
+    if not pos or not embedding:
+        raise Exception("Position and embedding required to find outer face")
+
+    visited: set[tuple[Any, Any]] = set()
+    faces: list[list[Any]] = []
+
+    for u in g.vertices():
+        for v in embedding[u]:
+            if (u, v) not in visited:
+                face: list[Any] = []
+                cu, cv = u, v
+                while (cu, cv) not in visited:
+                    visited.add((cu, cv))
+                    face.append(cu)
+                    neighbors = embedding[cv]
+                    cw = neighbors[(neighbors.index(cu) + 1) % len(neighbors)]
+                    cu, cv = cv, cw
+                faces.append(face)
+
+    def signed_area(face: list[Any]) -> float:
+        coords = [pos[v] for v in face]
+        n = len(coords)
+        return sum(
+            coords[i][0] * coords[(i + 1) % n][1] - coords[(i + 1) % n][0] * coords[i][1]
+            for i in range(n)
+        ) / 2
+
+    return min(faces, key=signed_area)
+
+
+def get_embedding_from_pos(g: Graph) -> dict[Any, list[Any]]:
+    """Compute a combinatorial planar embedding by sorting each vertex's neighbors by angle."""
+    pos_or_none = g.get_pos()
+    if pos_or_none is None:
+        raise Exception("Cannot get embedding without position")
+    pos: dict[Any, Any] = pos_or_none
+    return {
+        v: sorted(g.neighbors(v), key=lambda w: math.atan2(pos[w][1] - pos[v][1], pos[w][0] - pos[v][0]))
+        for v in g.vertices()
+    }