diff --git a/CHANGELOG.md b/CHANGELOG.md
new file mode 100644
--- /dev/null
+++ b/CHANGELOG.md
@@ -0,0 +1,38 @@
+# Changelog for `pure-cdt`
+
+All notable changes to this project will be documented in this file.
+
+The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),
+and this project adheres to the
+[Haskell Package Versioning Policy](https://pvp.haskell.org/).
+
+## 0.2.0.0 - 2026-09-04
+
+First public release. Version 0.1.0.0 existed only in the repository, under the name
+`triangulation`, and was never uploaded.
+
+- `triangulate`: Delaunay triangulation of a point set by divide and conquer, with
+  halves of at least 256 points evaluated in parallel through
+  `Control.Parallel.Strategies`.
+- `constrainedTriangulate`: triangulation of a polygon with holes. The boundary and
+  hole edges are forced in, triangles outside the region are dropped, and the boundary
+  may be non-convex. An input point on a polygon edge subdivides it.
+- `refine`: Delaunay refinement (Ruppert's algorithm) to a minimum angle and a maximum
+  triangle area, preserving the boundary and hole edges. `refineWithBudget` returns the
+  mesh reached when the insertion budget runs out.
+- `toMesh`, `boundaryEdges`: the triangulation as node coordinates and
+  counter-clockwise index triples, with the boundary edges identified — the form finite
+  element code and mesh file formats expect.
+- Exact geometric predicates: orientation and in-circle are decided by a
+  floating-point filter backed by `Rational` arithmetic, and the algorithm runs on
+  symbolically perturbed points (Simulation of Simplicity), so collinear and
+  cocircular input — grids, regular polygons, concentric shapes — is triangulated
+  correctly.
+- No partial functions and no exceptions in the library: a point set that spans no
+  triangle comes back as `Nothing`.
+- `Point`, `Edge`, `Triangle` and `Polygon` derive `Eq`, `Ord`, `Hashable` and
+  `NFData`. `Edge` and `Triangle` are abstract and normalise their vertices, so
+  equality and hashing do not depend on construction order.
+- Tests on `tasty`: QuickCheck properties for the triangulation (Delaunay condition,
+  planarity, triangle count, convex hull), for the region mesh and for the refinement,
+  next to unit tests for each geometric primitive. Benchmarks on `tasty-bench`.
diff --git a/LICENSE b/LICENSE
new file mode 100644
--- /dev/null
+++ b/LICENSE
@@ -0,0 +1,26 @@
+Copyright 2025-2026 Alex Elyukov
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are met:
+
+1. Redistributions of source code must retain the above copyright notice,
+   this list of conditions and the following disclaimer.
+
+2. Redistributions in binary form must reproduce the above copyright notice,
+   this list of conditions and the following disclaimer in the documentation
+   and/or other materials provided with the distribution.
+
+3. Neither the name of the copyright holder nor the names of its contributors
+   may be used to endorse or promote products derived from this software
+   without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
+ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
+WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR
+ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
+(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
+LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
+ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
diff --git a/README.md b/README.md
new file mode 100644
--- /dev/null
+++ b/README.md
@@ -0,0 +1,131 @@
+# pure-cdt
+
+[![CI](https://github.com/alexelyukov/triangulation/actions/workflows/ci.yml/badge.svg)](https://github.com/alexelyukov/triangulation/actions/workflows/ci.yml)
+[![Hackage](https://img.shields.io/hackage/v/pure-cdt.svg)](https://hackage.haskell.org/package/pure-cdt)
+
+Constrained Delaunay triangulation and mesh refinement in 2D, written in Haskell with
+no bindings to a C or C++ library.
+
+![a gear with a keyed bore, triangulated and refined to a 28° minimum angle](https://raw.githubusercontent.com/alexelyukov/triangulation/master/assets/gear.png)
+
+*A gear with a keyed bore: `constrainedTriangulate` on the two outlines alone, then
+`refine` to a 28° minimum angle and a maximum triangle area. The only input points are
+the polygon vertices.*
+
+## Usage
+
+```haskell
+import Data.List.NonEmpty (NonEmpty (..))
+import Triangulation
+
+-- Delaunay triangulation of a point set (Nothing if it spans no triangle).
+Just t = triangulate [Point 0 0, Point 4 0, Point 2 3, Point 2 1]
+triangles t  -- :: [Triangle]
+hull t       -- :: Polygon, the convex hull
+
+-- A polygon with holes: the boundary and hole edges are forced into the
+-- triangulation, triangles outside the boundary or inside a hole are dropped.
+-- The third argument holds interior points, if there are any.
+boundary = Polygon (Point 0 0 :| [Point 0 10, Point 10 10, Point 10 0])
+hole     = Polygon (Point 4 4 :| [Point 4 6, Point 6 6, Point 6 4])
+Just ts = constrainedTriangulate boundary [hole] []
+
+-- Mesh refinement (Ruppert's algorithm): insert vertices until no angle is
+-- below 25° and no triangle is larger than 0.5, keeping the polygon edges.
+Just fine = refine defaultQuality {minAngle = 25, maxArea = Just 0.5} ts
+
+-- The indexed form finite element code wants: nodes, and triangles as
+-- counter-clockwise triples of node indices; plus the boundary as index pairs.
+mesh = toMesh fine
+nodes mesh         -- :: Vector Point
+elements mesh      -- :: Vector (Int, Int, Int)
+boundaryEdges mesh -- :: [(Int, Int)], the edges that belong to a single triangle
+```
+
+Two conventions to know before the first call: the y axis points up, and polygons are
+wound clockwise, so that walking along a boundary the interior is on the right. `Edge`
+and `Triangle` are abstract; build them with `mkEdge`/`mkTriangle` and take them apart
+with the read-only patterns `Edge a b`/`Triangle a b c`.
+
+## What this package is for
+
+It does not compete on speed. A mature Delaunay implementation in C or C++ is an order
+of magnitude or two faster, and nothing here will close that gap. What it offers
+instead:
+
+* **A pure API.** A handful of ordinary functions that can be mapped over a list or
+  used inside a QuickCheck property. A point set that spans no triangle comes back as
+  `Nothing`: there are no partial functions and no exceptions in the library.
+* **No toolchain beyond GHC.** The dependencies are `base`, `deepseq`, `hashable`,
+  `parallel`, `random`, `unordered-containers` and `vector`, so the package builds
+  wherever GHC does, including under Nix, when cross-compiling, and on the JavaScript
+  and WebAssembly backends.
+* **Degenerate input handled, not merely tolerated.** The predicates are exact and the
+  algorithm runs on symbolically perturbed points (Edelsbrunner and Mücke's Simulation
+  of Simplicity), so lattices, regular polygons and concentric shapes are triangulated
+  correctly instead of coming out with holes or overlaps. They cost about three times
+  as long as points in general position, and that is a constant factor, not a worse
+  complexity.
+* **A small readable core.** The code reads as geometry rather than as bookkeeping of
+  array indices.
+
+Making it faster would mean giving up the last of those. Passing indices into an array
+instead of `Point` values through the algorithm is likely worth a factor of two or
+three, and that trade has been declined deliberately.
+
+## Performance
+
+Wall-clock time, median of three runs on an Intel Core i7-14700KF (eight performance
+cores plus twelve efficiency cores), GHC 9.10.3:
+
+| Points | Triangles | 1 core | 4 cores | 8 cores |
+| ---: | ---: | ---: | ---: | ---: |
+| 50 000 | 99 967 | 0.41 s | 0.17 s (2.4×) | 0.13 s (3.2×) |
+| 100 000 | 199 965 | 0.86 s | 0.37 s (2.3×) | 0.30 s (2.9×) |
+| 500 000 | 999 966 | 5.11 s | 2.57 s (2.0×) | 2.15 s (2.4×) |
+| 1 000 000 | 1 999 968 | 11.4 s | 5.80 s (2.0×) | 5.16 s (2.2×) |
+
+Halves of a large point set are triangulated in parallel and then merged. The speedup
+plateaus between eight and twelve threads and falls away after that, because the merges
+along the spine of the recursion are sequential and the parallel collector has more
+capabilities to synchronise.
+
+## Building
+
+```
+stack build          # library and the example executable
+stack test           # tasty: unit tests and QuickCheck properties
+stack bench          # tasty-bench: triangulate, refine, toMesh
+stack exec pure-cdt-examples   # renders the examples to assets/*.png
+
+# wall-clock scaling, and the same on a regular lattice
+stack bench pure-cdt:pure-cdt-scaling --ba '50000 500000 +RTS -N8'
+stack bench pure-cdt:pure-cdt-scaling --ba 'lattice 100000 +RTS -N8'
+```
+
+`tasty-bench` reports CPU time, which grows with the number of cores and so says
+nothing about a parallel speedup; the table above comes from `pure-cdt-scaling`, which
+measures the clock. The example renderer sits behind the cabal flag `examples`, off by
+default, so that the library carries no image dependencies; `stack.yaml` turns it on
+for development, and plain cabal takes `-f examples`.
+
+Formatting is `fourmolu`, linting is `hlint`; both run in CI.
+
+## Limitations
+
+* Polygons must be simple and must not cross each other. A hole may touch the boundary
+  at a vertex, in which case the corner is cut off.
+* Input whose points all lie on one line has no triangulation, and `triangulate`
+  returns `Nothing` for it.
+* Ruppert's refinement is guaranteed to terminate for angle bounds up to about 20.7°
+  when no two polygon edges meet at less than 60°, and in practice works up to about
+  30°. A demanding bound next to a small input angle can exhaust the insertion budget;
+  `refine` then returns `Nothing`, and `refineWithBudget` returns the mesh it reached.
+* `refine` works on the triangles of a region, so an interior constraint edge with
+  triangles on both sides is not preserved. It splits segments at points computed in
+  floating point, so a vertex it adds to the boundary lies on the original polygon edge
+  only up to rounding.
+
+## License
+
+BSD-3-Clause.
diff --git a/Setup.hs b/Setup.hs
new file mode 100644
--- /dev/null
+++ b/Setup.hs
@@ -0,0 +1,2 @@
+import Distribution.Simple
+main = defaultMain
diff --git a/app/Drawer.hs b/app/Drawer.hs
new file mode 100644
--- /dev/null
+++ b/app/Drawer.hs
@@ -0,0 +1,78 @@
+-- | PNG rendering of triangulations. The library's y axis points up; here it
+-- is flipped to pixel rows.
+module Drawer (
+  drawTriangulation,
+) where
+
+import Codec.Picture (Image, PixelRGBA8 (..), writePng)
+import Data.Foldable (traverse_)
+import Graphics.Rasterific (
+  Cap (CapRound),
+  Drawing,
+  Join (JoinRound),
+  V2 (..),
+  circle,
+  fill,
+  line,
+  renderDrawing,
+  stroke,
+  withTexture,
+ )
+import Graphics.Rasterific.Texture (uniformTexture)
+import Triangulation.Geometry (
+  Edge (..),
+  Point (..),
+  Polygon (..),
+  Triangle (..),
+  cyclicPairs,
+  mkEdge,
+ )
+
+drawTriangulation :: FilePath -> [Triangle] -> [[Point]] -> [Polygon] -> IO ()
+drawTriangulation path triangles polylines polygons =
+  writePng path $ drawBackground $ do
+    traverse_ drawPolyline polylines
+    traverse_ drawPolygon polygons
+    traverse_ drawTriangle triangles
+
+-- | Side of the square canvas, in pixels.
+canvasSize :: Int
+canvasSize = 2000
+
+drawBackground :: Drawing PixelRGBA8 () -> Image PixelRGBA8
+drawBackground =
+  let backgroundColor = PixelRGBA8 255 255 255 255
+      drawColor = PixelRGBA8 0x00 0x00 0x00 255
+   in renderDrawing canvasSize canvasSize backgroundColor . withTexture (uniformTexture drawColor)
+
+drawTriangle :: Triangle -> Drawing PixelRGBA8 ()
+drawTriangle (Triangle p1 p2 p3) = traverse_ drawEdge [mkEdge p1 p2, mkEdge p2 p3, mkEdge p3 p1]
+
+drawEdge :: Edge -> Drawing PixelRGBA8 ()
+drawEdge (Edge p0 p1) = do
+  drawSegment 1 (p0, p1)
+  drawPoint p0
+  drawPoint p1
+
+drawPolygon :: Polygon -> Drawing PixelRGBA8 ()
+drawPolygon (Polygon points) = do
+  traverse_ (drawSegment 3) (cyclicPairs points)
+  traverse_ drawPoint points
+
+drawPolyline :: [Point] -> Drawing PixelRGBA8 ()
+drawPolyline points = do
+  traverse_ (drawSegment 3) (zip points (drop 1 points))
+  traverse_ drawPoint points
+
+drawSegment :: Float -> (Point, Point) -> Drawing PixelRGBA8 ()
+drawSegment width (p0, p1) =
+  stroke width JoinRound (CapRound, CapRound) $ line (preparePoint p0) (preparePoint p1)
+
+drawPoint :: Point -> Drawing PixelRGBA8 ()
+drawPoint pc =
+  let pointColor = PixelRGBA8 0xFF 0x00 0x00 255
+   in withTexture (uniformTexture pointColor) $ fill $ circle (preparePoint pc) 3
+
+-- | Map a library point (y up) to a pixel position (y down).
+preparePoint :: Point -> V2 Float
+preparePoint (Point x y) = V2 (realToFrac x) (realToFrac (fromIntegral canvasSize - y))
diff --git a/app/Examples/Circle.hs b/app/Examples/Circle.hs
new file mode 100644
--- /dev/null
+++ b/app/Examples/Circle.hs
@@ -0,0 +1,37 @@
+module Examples.Circle (
+  drawCircle,
+  drawTorus,
+) where
+
+import Data.List.NonEmpty (NonEmpty (..))
+import Data.List.NonEmpty qualified as NE
+import Drawer (drawTriangulation)
+import Examples.Common (canvas, orFail)
+import Triangulation (Point (..), Polygon (..), constrainedTriangulate, vertices)
+import Triangulation.Generator (generatePointsWithDistance)
+
+drawCircle :: FilePath -> IO ()
+drawCircle path = do
+  triangles <- orFail $ constrainedTriangulate circle [] points
+  drawTriangulation path triangles [] [circle]
+  where
+    circle = circleOf (Point 1000 1000) 800 150
+    points = generatePointsWithDistance 2 8000 20 canvas (vertices circle)
+
+drawTorus :: FilePath -> IO ()
+drawTorus path = do
+  triangles <- orFail $ constrainedTriangulate outer [inner] points
+  drawTriangulation path triangles [] [outer, inner]
+  where
+    outer = circleOf (Point 1000 1000) 800 100
+    inner = circleOf (Point 1000 1000) 400 50
+    points = generatePointsWithDistance 1 4000 20 canvas (vertices outer ++ vertices inner)
+
+-- | A regular @n@-gon inscribed in the circle, wound clockwise.
+circleOf :: Point -> Double -> Int -> Polygon
+circleOf (Point x0 y0) radius n =
+  let angleStep = 2 * pi / fromIntegral n
+      angles = NE.map ((angleStep *) . fromIntegral) (0 :| [1 .. n - 1])
+   in Polygon
+        $ NE.reverse
+        $ NE.map (\angle -> Point (x0 + radius * cos angle) (y0 + radius * sin angle)) angles
diff --git a/app/Examples/Common.hs b/app/Examples/Common.hs
new file mode 100644
--- /dev/null
+++ b/app/Examples/Common.hs
@@ -0,0 +1,16 @@
+-- | Shared pieces of the examples.
+module Examples.Common (
+  canvas,
+  orFail,
+) where
+
+import Triangulation.Generator (Rectangle (..))
+import Triangulation.Geometry (Point (..))
+
+-- | The area random points are drawn from; matches the canvas of "Drawer".
+canvas :: Rectangle
+canvas = Rectangle {minCorner = Point 0 0, maxCorner = Point 2000 2000}
+
+-- | Unwrap a triangulation result, aborting the program if it failed.
+orFail :: Maybe a -> IO a
+orFail = maybe (fail "triangulation failed: fewer than three points or inconsistent hulls") pure
diff --git a/app/Examples/Gear.hs b/app/Examples/Gear.hs
new file mode 100644
--- /dev/null
+++ b/app/Examples/Gear.hs
@@ -0,0 +1,79 @@
+-- | A gear with a keyed bore, meshed for finite elements: the outline and the
+-- hole are the only input, the density comes from the refinement.
+module Examples.Gear (
+  drawGear,
+) where
+
+import Data.List.NonEmpty qualified as NE
+import Drawer (drawTriangulation)
+import Examples.Common (orFail)
+import Triangulation (
+  Point (..),
+  Polygon (..),
+  Quality (..),
+  constrainedTriangulate,
+  defaultQuality,
+  refine,
+ )
+
+drawGear :: FilePath -> IO ()
+drawGear path = do
+  coarse <- orFail $ constrainedTriangulate gear [bore] []
+  triangles <- orFail $ refine defaultQuality {minAngle = 28, maxArea = Just 1200} coarse
+  drawTriangulation path triangles [] [gear, bore]
+
+centre :: Point
+centre = Point 1000 1000
+
+-- | The outline: 16 trapezoidal teeth on a root circle.
+gear :: Polygon
+gear = clockwise (concatMap tooth [0 .. teeth - 1])
+  where
+    teeth = 16 :: Int
+    (root, tip) = (700, 850)
+    pitch = 2 * pi / fromIntegral teeth
+    tooth i =
+      let a0 = fromIntegral i * pitch
+          at = polar
+       in [ at root a0
+          , at root (a0 + 0.1 * pitch)
+          , at root (a0 + 0.2 * pitch)
+          , at tip (a0 + 0.3 * pitch)
+          , at tip (a0 + 0.4 * pitch)
+          , at tip (a0 + 0.5 * pitch)
+          , at tip (a0 + 0.6 * pitch)
+          , at tip (a0 + 0.7 * pitch)
+          , at root (a0 + 0.8 * pitch)
+          , at root (a0 + 0.9 * pitch)
+          ]
+
+-- | The hole: a circle with a rectangular keyway cut into it at the top.
+bore :: Polygon
+bore = clockwise (concatMap arcPoint [0 .. n - 1])
+  where
+    n = 48 :: Int
+    radius = 260
+    (halfWidth, depth) = (45, 320)
+    step = 2 * pi / fromIntegral n
+    keywayEdge = asin (halfWidth / radius) -- half-angle of the keyway opening
+    arcPoint i
+      | a > pi / 2 - keywayEdge && a < pi / 2 + keywayEdge = []
+      | a <= pi / 2 && a + step > pi / 2 - keywayEdge =
+          [ polar radius a
+          , Point (1000 + halfWidth) (1000 + rim)
+          , Point (1000 + halfWidth) (1000 + depth)
+          , Point (1000 - halfWidth) (1000 + depth)
+          , Point (1000 - halfWidth) (1000 + rim)
+          ]
+      | otherwise = [polar radius a]
+      where
+        a = fromIntegral i * step
+    rim = sqrt (radius * radius - halfWidth * halfWidth)
+
+-- | A point at the given radius and angle around the centre.
+polar :: Double -> Double -> Point
+polar r a = Point (px centre + r * cos a) (py centre + r * sin a)
+
+-- | Points listed counter-clockwise (increasing angle) as a clockwise polygon.
+clockwise :: [Point] -> Polygon
+clockwise = Polygon . NE.reverse . NE.fromList
diff --git a/app/Examples/Simple.hs b/app/Examples/Simple.hs
new file mode 100644
--- /dev/null
+++ b/app/Examples/Simple.hs
@@ -0,0 +1,13 @@
+module Examples.Simple (
+  drawSimple,
+) where
+
+import Drawer (drawTriangulation)
+import Examples.Common (canvas, orFail)
+import Triangulation (triangles, triangulate)
+import Triangulation.Generator (generatePoints)
+
+drawSimple :: FilePath -> IO ()
+drawSimple path = do
+  triangulation <- orFail $ triangulate (generatePoints 2 8 canvas)
+  drawTriangulation path (triangles triangulation) [] []
diff --git a/app/Examples/Solenoid.hs b/app/Examples/Solenoid.hs
new file mode 100644
--- /dev/null
+++ b/app/Examples/Solenoid.hs
@@ -0,0 +1,63 @@
+module Examples.Solenoid (
+  drawSolenoid,
+  drawRefinedSolenoid,
+) where
+
+import Data.List.NonEmpty (NonEmpty (..))
+import Drawer (drawTriangulation)
+import Examples.Common (canvas, orFail)
+import Triangulation (
+  Point (..),
+  Polygon (..),
+  Quality (..),
+  constrainedTriangulate,
+  defaultQuality,
+  refine,
+  vertices,
+ )
+import Triangulation.Generator (Rectangle (..), generatePointsWithDistance)
+
+drawSolenoid :: FilePath -> IO ()
+drawSolenoid path = do
+  triangles <- orFail $ constrainedTriangulate outer [innerLeft, innerRight] points
+  drawTriangulation path triangles [] [outer, innerLeft, innerRight]
+  where
+    outer = rectangleOf Rectangle {minCorner = Point 100 200, maxCorner = Point 1900 1800} (30, 20)
+    innerLeft = rectangleOf Rectangle {minCorner = Point 300 400, maxCorner = Point 900 1600} (10, 20)
+    innerRight = rectangleOf Rectangle {minCorner = Point 1100 400, maxCorner = Point 1700 1600} (10, 20)
+    points = generatePointsWithDistance 1 4000 20 canvas (concatMap vertices [outer, innerLeft, innerRight])
+
+-- | The boundary of a rectangle as a clockwise polygon, with @numH@ extra
+-- points on each horizontal side and @numV@ on each vertical side.
+rectangleOf :: Rectangle -> (Int, Int) -> Polygon
+rectangleOf Rectangle {minCorner = Point x1 y1, maxCorner = Point x2 y2} (numH, numV) =
+  let lowerLeft = Point x1 y1
+      lowerRight = Point x2 y1
+      upperRight = Point x2 y2
+      upperLeft = Point x1 y2
+      sideRight = pointsBetween upperRight lowerRight numV ++ [lowerRight]
+      sideBottom = pointsBetween lowerRight lowerLeft numH ++ [lowerLeft]
+      sideLeft = pointsBetween lowerLeft upperLeft numV ++ [upperLeft]
+      sideTop = pointsBetween upperLeft upperRight numH
+   in Polygon (upperRight :| (sideRight ++ sideBottom ++ sideLeft ++ sideTop))
+
+-- | @n@ equally spaced points strictly between two points.
+pointsBetween :: Point -> Point -> Int -> [Point]
+pointsBetween (Point x1 y1) (Point x2 y2) n =
+  let dx = (x2 - x1) / fromIntegral (n + 1)
+      dy = (y2 - y1) / fromIntegral (n + 1)
+   in [Point (x1 + dx * fromIntegral i) (y1 + dy * fromIntegral i) | i <- [1 .. n]]
+
+-- | The same section, meshed for finite elements: the polygons with only a
+-- few points on each side, refined to a 28° angle bound and a maximum
+-- triangle area, so that the density comes from the refinement rather than
+-- from random points.
+drawRefinedSolenoid :: FilePath -> IO ()
+drawRefinedSolenoid path = do
+  coarse <- orFail $ constrainedTriangulate outer [innerLeft, innerRight] []
+  triangles <- orFail $ refine defaultQuality {minAngle = 28, maxArea = Just 1500} coarse
+  drawTriangulation path triangles [] [outer, innerLeft, innerRight]
+  where
+    outer = rectangleOf Rectangle {minCorner = Point 100 200, maxCorner = Point 1900 1800} (3, 2)
+    innerLeft = rectangleOf Rectangle {minCorner = Point 300 400, maxCorner = Point 900 1600} (1, 2)
+    innerRight = rectangleOf Rectangle {minCorner = Point 1100 400, maxCorner = Point 1700 1600} (1, 2)
diff --git a/app/Main.hs b/app/Main.hs
new file mode 100644
--- /dev/null
+++ b/app/Main.hs
@@ -0,0 +1,18 @@
+module Main (
+  main,
+) where
+
+import Examples.Circle (drawCircle, drawTorus)
+import Examples.Gear (drawGear)
+import Examples.Simple (drawSimple)
+import Examples.Solenoid (drawRefinedSolenoid, drawSolenoid)
+
+main :: IO ()
+main = do
+  drawGear "assets/gear.png"
+  drawSimple "assets/triangulation_simple.png"
+  drawCircle "assets/triangulation_circle.png"
+  drawTorus "assets/triangulation_torus.png"
+  drawSolenoid "assets/triangulation_solenoid.png"
+  drawRefinedSolenoid "assets/triangulation_refined.png"
+  return ()
diff --git a/bench/Main.hs b/bench/Main.hs
new file mode 100644
--- /dev/null
+++ b/bench/Main.hs
@@ -0,0 +1,39 @@
+module Main (main) where
+
+import Data.List.NonEmpty (NonEmpty (..))
+import Data.Maybe (fromMaybe)
+import Test.Tasty.Bench (bench, bgroup, defaultMain, nf)
+import Triangulation (
+  Point (..),
+  Polygon (..),
+  Quality (..),
+  constrainedTriangulate,
+  defaultQuality,
+  refine,
+  toMesh,
+  triangles,
+  triangulate,
+ )
+import Triangulation.Generator (Rectangle (..), generatePoints)
+
+main :: IO ()
+main =
+  defaultMain
+    [ bgroup
+        "triangulate"
+        [ bench (show n ++ " points") $ nf (fmap (length . triangles) . triangulate) (points n)
+        | n <- [1000, 8000, 50000]
+        ]
+    , bgroup
+        "refine a square with a hole"
+        [ bench ("20°, area " ++ show area) $
+            nf (fmap length . refine defaultQuality {maxArea = Just area}) coarse
+        | area <- [4000, 1000]
+        ]
+    , bench "toMesh of 8000 points" $ nf (fmap (toMesh . triangles) . triangulate) (points 8000)
+    ]
+  where
+    points n = generatePoints 7 n Rectangle {minCorner = Point 0 0, maxCorner = Point 2000 2000}
+    canvas = Polygon (Point 0 0 :| [Point 0 2000, Point 2000 2000, Point 2000 0])
+    hole = Polygon (Point 600 500 :| [Point 600 1300, Point 1400 1300, Point 1400 500])
+    coarse = fromMaybe [] (constrainedTriangulate canvas [hole] (points 200))
diff --git a/bench/Scale.hs b/bench/Scale.hs
new file mode 100644
--- /dev/null
+++ b/bench/Scale.hs
@@ -0,0 +1,97 @@
+-- | How the triangulation scales with the number of cores, measured in
+-- wall-clock time.
+--
+-- The @pure-cdt-bench@ benchmarks run on @tasty-bench@, which reports CPU
+-- time; that number grows with the number of cores and says nothing about a
+-- parallel speedup. This program times the triangulation by the clock
+-- instead. Each run gets a freshly generated point set, so no work is shared
+-- between runs, and generating the points is not counted.
+--
+-- > stack bench pure-cdt:pure-cdt-scaling --ba '100000 500000 +RTS -N8'
+--
+-- With no sizes given it triangulates 50 000 points. Pass @+RTS -N\<k\>@ to
+-- choose how many cores to use.
+--
+-- A leading @lattice@ places the points on a regular lattice instead of
+-- drawing them at random, which is the worst case for degeneracy and the
+-- likeliest input of a finite element model: every row, column and diagonal
+-- is collinear and the corners of every cell are cocircular.
+--
+-- > stack bench pure-cdt:pure-cdt-scaling --ba 'lattice 100000 +RTS -N8'
+module Main (main) where
+
+import Control.DeepSeq (force)
+import Control.Exception (evaluate)
+import Data.List (sort)
+import GHC.Clock (getMonotonicTime)
+import GHC.Conc (getNumCapabilities)
+import System.Environment (getArgs)
+import Text.Printf (printf)
+import Triangulation (Point (..), triangles, triangulate)
+import Triangulation.Generator (Rectangle (..), generatePoints)
+
+repetitions :: Int
+repetitions = 3
+
+main :: IO ()
+main = do
+  args <- getArgs
+  cores <- getNumCapabilities
+  let (shape, rest) = case args of
+        "lattice" : more -> (Lattice, more)
+        more -> (Random, more)
+      sizes = case map read rest of
+        [] -> [50000]
+        given -> given
+  printf
+    "%s points, cores: %d, %d runs per size, median reported\n"
+    (if shape == Lattice then "lattice" else "random")
+    cores
+    repetitions
+  mapM_ (report shape) sizes
+
+-- | Where the points come from.
+data Shape = Random | Lattice
+  deriving (Eq)
+
+report :: Shape -> Int -> IO ()
+report shape n = do
+  results <- mapM (timeOne shape n) [1 .. repetitions]
+  let times = sort (map fst results)
+  case (median times, map snd results) of
+    (Just middle, count : _) ->
+      printf
+        "%8d points -> %8d triangles: %7.0f ms  (runs: %s)\n"
+        n
+        count
+        (middle * 1000)
+        (unwords (map (printf "%.0f" . (* 1000)) times))
+    _ -> printf "%8d points: no result\n" n
+
+-- | Generate a point set from the seed, then time triangulating it.
+timeOne :: Shape -> Int -> Int -> IO (Double, Int)
+timeOne shape n seed = do
+  points <- evaluate (force (pointsOf shape n seed))
+  before <- getMonotonicTime
+  result <- evaluate (force (triangulate points))
+  count <- evaluate (maybe 0 (length . triangles) result)
+  after <- getMonotonicTime
+  pure (after - before, count)
+
+median :: [Double] -> Maybe Double
+median xs = case drop (length xs `div` 2) xs of
+  middle : _ -> Just middle
+  [] -> Nothing
+
+-- | @n@ points of the given shape. A lattice ignores the seed: there is only
+-- one lattice of a given size, and its regularity is the point.
+pointsOf :: Shape -> Int -> Int -> [Point]
+pointsOf Random n seed = generatePoints seed n canvas
+pointsOf Lattice n _ =
+  take n [Point (fromIntegral i * step) (fromIntegral j * step) | i <- [0 .. side], j <- [0 .. side]]
+  where
+    side = ceiling (sqrt (fromIntegral n :: Double)) :: Int
+    step = 300 :: Double
+
+canvas :: Rectangle
+canvas = Rectangle {minCorner = Point 0 0, maxCorner = Point 100000 100000}
diff --git a/pure-cdt.cabal b/pure-cdt.cabal
new file mode 100644
--- /dev/null
+++ b/pure-cdt.cabal
@@ -0,0 +1,162 @@
+cabal-version: 2.2
+
+-- This file has been generated from package.yaml by hpack version 0.39.6.
+--
+-- see: https://github.com/sol/hpack
+
+name:           pure-cdt
+version:        0.2.0.0
+synopsis:       Constrained Delaunay triangulation and mesh refinement in 2D, in pure Haskell
+description:    Delaunay triangulation of a set of points in the plane, computed by divide
+                and conquer and parallelised with @Control.Parallel.Strategies@; constrained
+                triangulation of a polygon with holes, where the boundary edges are forced
+                into the triangulation and the triangles outside the region are dropped;
+                Delaunay refinement (Ruppert's algorithm) to a minimum angle and a maximum
+                triangle area; and an indexed mesh (node coordinates, index triples,
+                boundary edges) for finite element code.
+                .
+                Everything is written in Haskell, with a pure API and no bindings to a C or
+                C++ library, so the package builds wherever GHC does.
+                .
+                The geometric predicates (orientation, in-circle, segment crossing) are
+                exact: a floating-point filter decides the easy cases and @Rational@
+                arithmetic the rest, so collinear and cocircular inputs are handled
+                correctly.
+                .
+                Start with the "Triangulation" module. Please see the README on GitHub at
+                <https://github.com/alexelyukov/triangulation#readme>.
+category:       Geometry
+homepage:       https://github.com/alexelyukov/triangulation#readme
+bug-reports:    https://github.com/alexelyukov/triangulation/issues
+author:         Alex Elyukov
+maintainer:     alexelyukov@gmail.com
+copyright:      2025-2026 Alex Elyukov
+license:        BSD-3-Clause
+license-file:   LICENSE
+build-type:     Simple
+tested-with:
+    GHC == 9.6.7 || == 9.10.3 || == 9.14.1
+extra-doc-files:
+    README.md
+    CHANGELOG.md
+
+source-repository head
+  type: git
+  location: https://github.com/alexelyukov/triangulation
+
+flag examples
+  description: Build the executable that renders the example pictures
+  manual: True
+  default: False
+
+library
+  exposed-modules:
+      Triangulation
+      Triangulation.Check
+      Triangulation.Constrained
+      Triangulation.Flip
+      Triangulation.Generator
+      Triangulation.Geometry
+      Triangulation.Geometry.Edge
+      Triangulation.Geometry.Exact
+      Triangulation.Geometry.Point
+      Triangulation.Geometry.Polygon
+      Triangulation.Geometry.Ring
+      Triangulation.Geometry.Triangle
+      Triangulation.Leaf
+      Triangulation.Merge
+      Triangulation.Mesh
+      Triangulation.Parallel
+      Triangulation.Refine
+      Triangulation.Repair
+      Triangulation.Store
+      Triangulation.Types
+  other-modules:
+      Paths_pure_cdt
+  autogen-modules:
+      Paths_pure_cdt
+  hs-source-dirs:
+      src
+  ghc-options: -Wall -Wcompat -Widentities -Wincomplete-record-updates -Wincomplete-uni-patterns -Wmissing-export-lists -Wmissing-home-modules -Wpartial-fields -Wredundant-constraints -Wunused-packages
+  build-depends:
+      base >=4.18 && <5
+    , deepseq >=1.4 && <1.6
+    , hashable >=1.4 && <1.6
+    , parallel >=3.2 && <3.4
+    , random >=1.2 && <1.4
+    , unordered-containers >=0.2.19 && <0.3
+    , vector ==0.13.*
+  default-language: GHC2021
+
+executable pure-cdt-examples
+  main-is: Main.hs
+  other-modules:
+      Drawer
+      Examples.Circle
+      Examples.Common
+      Examples.Gear
+      Examples.Simple
+      Examples.Solenoid
+      Paths_pure_cdt
+  autogen-modules:
+      Paths_pure_cdt
+  hs-source-dirs:
+      app
+  ghc-options: -Wall -Wcompat -Widentities -Wincomplete-record-updates -Wincomplete-uni-patterns -Wmissing-export-lists -Wmissing-home-modules -Wpartial-fields -Wredundant-constraints -Wunused-packages -threaded -rtsopts -with-rtsopts=-maxN8
+  build-depends:
+      JuicyPixels ==3.3.*
+    , Rasterific ==0.7.*
+    , base >=4.18 && <5
+    , pure-cdt
+  default-language: GHC2021
+  if !flag(examples)
+    buildable: False
+
+test-suite pure-cdt-test
+  type: exitcode-stdio-1.0
+  main-is: Main.hs
+  other-modules:
+      Test.Geometry
+      Test.Mesh
+      Test.Refine
+      Test.Ring
+      Test.Triangulation
+      Paths_pure_cdt
+  autogen-modules:
+      Paths_pure_cdt
+  hs-source-dirs:
+      test
+  ghc-options: -Wall -Wcompat -Widentities -Wincomplete-record-updates -Wincomplete-uni-patterns -Wmissing-export-lists -Wmissing-home-modules -Wpartial-fields -Wredundant-constraints -Wunused-packages -threaded -rtsopts -with-rtsopts=-maxN8
+  build-depends:
+      base >=4.18 && <5
+    , hashable >=1.4 && <1.6
+    , pure-cdt
+    , tasty >=1.4 && <1.6
+    , tasty-hunit ==0.10.*
+    , tasty-quickcheck >=0.10 && <0.12
+    , vector ==0.13.*
+  default-language: GHC2021
+
+benchmark pure-cdt-bench
+  type: exitcode-stdio-1.0
+  main-is: Main.hs
+  hs-source-dirs:
+      bench
+  ghc-options: -Wall -Wcompat -Widentities -Wincomplete-record-updates -Wincomplete-uni-patterns -Wmissing-export-lists -Wmissing-home-modules -Wpartial-fields -Wredundant-constraints -Wunused-packages -threaded -rtsopts -with-rtsopts=-maxN8
+  build-depends:
+      base >=4.18 && <5
+    , pure-cdt
+    , tasty-bench >=0.3 && <0.5
+  default-language: GHC2021
+
+benchmark pure-cdt-scaling
+  type: exitcode-stdio-1.0
+  main-is: Scale.hs
+  hs-source-dirs:
+      bench
+  ghc-options: -Wall -Wcompat -Widentities -Wincomplete-record-updates -Wincomplete-uni-patterns -Wmissing-export-lists -Wmissing-home-modules -Wpartial-fields -Wredundant-constraints -Wunused-packages -threaded -rtsopts
+  build-depends:
+      base >=4.18 && <5
+    , deepseq >=1.4 && <1.6
+    , pure-cdt
+  default-language: GHC2021
diff --git a/src/Triangulation.hs b/src/Triangulation.hs
new file mode 100644
--- /dev/null
+++ b/src/Triangulation.hs
@@ -0,0 +1,34 @@
+-- | Constrained Delaunay triangulation and mesh refinement in 2D.
+--
+-- * 'triangulate' builds the Delaunay triangulation of a point set by
+--   divide and conquer.
+-- * 'constrainedTriangulate' triangulates a polygon with holes.
+-- * 'refine' inserts vertices until every triangle meets an angle and area
+--   bound (Ruppert's algorithm), keeping the boundary and hole edges.
+-- * 'toMesh' numbers the vertices and gives the triangles as index triples,
+--   the form finite element code and mesh file formats expect.
+--
+-- Coordinates follow the mathematical convention (y axis up); polygons are
+-- wound clockwise.
+module Triangulation (
+  Triangulation (..),
+  triangles,
+  triangulate,
+  constrainedTriangulate,
+  Quality (..),
+  defaultQuality,
+  refine,
+  refineWithBudget,
+  Mesh (..),
+  toMesh,
+  fromMesh,
+  boundaryEdges,
+  module Triangulation.Geometry,
+) where
+
+import Triangulation.Constrained (constrainedTriangulate)
+import Triangulation.Geometry
+import Triangulation.Mesh (Mesh (..), boundaryEdges, fromMesh, toMesh)
+import Triangulation.Parallel (triangulate)
+import Triangulation.Refine (Quality (..), defaultQuality, refine, refineWithBudget)
+import Triangulation.Types (Triangulation (..), triangles)
diff --git a/src/Triangulation/Check.hs b/src/Triangulation/Check.hs
new file mode 100644
--- /dev/null
+++ b/src/Triangulation/Check.hs
@@ -0,0 +1,43 @@
+-- | Validity checks on sets of triangles. Quadratic; meant for tests and debugging.
+module Triangulation.Check (
+  isDelaunay,
+  isLocallyDelaunay,
+  hasNoIntersections,
+) where
+
+import Data.HashSet qualified as HS
+import Data.Maybe (isNothing)
+import Triangulation.Geometry.Edge (Edge (..), intersection)
+import Triangulation.Geometry.Triangle (Triangle (..), isOutsideCircumcircle, triangleEdges)
+import Triangulation.Store qualified as Store
+
+-- | No vertex lies strictly inside the circumcircle of any triangle: the
+-- (unconstrained) Delaunay property.
+isDelaunay :: [Triangle] -> Bool
+isDelaunay ts =
+  let points = HS.toList . HS.fromList $ concatMap (\(Triangle p1 p2 p3) -> [p1, p2, p3]) ts
+      others (Triangle p1 p2 p3) = filter (`notElem` [p1, p2, p3]) points
+   in all (\t -> all (`isOutsideCircumcircle` t) (others t)) ts
+
+-- | Every edge shared by two triangles is locally Delaunay: the apex of each
+-- triangle lies outside (or on) the circumcircle of the other. Together with
+-- the constraint that segments are never flipped this characterises a
+-- constrained Delaunay triangulation.
+isLocallyDelaunay :: [Triangle] -> Bool
+isLocallyDelaunay ts = all locallyDelaunay (Store.edges store)
+  where
+    store = foldr Store.insert Store.empty ts
+    locallyDelaunay e = case Store.trianglesOn e store of
+      [t1, t2] -> isOutsideCircumcircle (apex e t2) t1 && isOutsideCircumcircle (apex e t1) t2
+      _ -> True
+    apex (Edge a b) (Triangle x y z) = case filter (`notElem` [a, b]) [x, y, z] of
+      p : _ -> p
+      [] -> x
+
+-- | No two edges cross.
+hasNoIntersections :: [Triangle] -> Bool
+hasNoIntersections = go
+  where
+    go [] = True
+    go (t : ts) = all (notCrossing t) ts && go ts
+    notCrossing t1 t2 = and [isNothing (intersection e1 e2) | e1 <- triangleEdges t1, e2 <- triangleEdges t2]
diff --git a/src/Triangulation/Constrained.hs b/src/Triangulation/Constrained.hs
new file mode 100644
--- /dev/null
+++ b/src/Triangulation/Constrained.hs
@@ -0,0 +1,92 @@
+-- | Constrained triangulation: forcing given edges into a Delaunay
+-- triangulation, and triangulating a polygon with holes.
+module Triangulation.Constrained (
+  constrainedTriangulate,
+  forceEdges,
+) where
+
+import Data.HashSet qualified as HS
+import Data.List qualified as List
+import Triangulation.Flip (legalize)
+import Triangulation.Geometry.Edge (Edge (..), intersection, mkEdge)
+import Triangulation.Geometry.Point (
+  Orientation (..),
+  Point (..),
+  manhattanDistance,
+  orientation,
+  turn,
+ )
+import Triangulation.Geometry.Polygon (Polygon, polygonEdges, vertices)
+import Triangulation.Geometry.Triangle (
+  Triangle,
+  isValidCandidate,
+  mkTriangle,
+  triangleEdges,
+  trianglesInside,
+  trianglesOutside,
+ )
+import Triangulation.Parallel (triangulate)
+import Triangulation.Repair (repairDegeneracies)
+import Triangulation.Store (Store)
+import Triangulation.Store qualified as Store
+import Triangulation.Types (Triangulation (..))
+
+-- | Triangulate the region inside the boundary polygon and outside the holes.
+-- The vertices of the polygons are always part of the triangulation; the
+-- given points (which may repeat them) are added. Polygon edges are forced
+-- into the triangulation and triangles outside the region are dropped.
+constrainedTriangulate :: Polygon -> [Polygon] -> [Point] -> Maybe [Triangle]
+constrainedTriangulate boundary holes points = do
+  triangulation <- triangulate allPoints
+  let store = triangleStore triangulation
+      edges = concatMap (subdivide allPoints) (concatMap polygonEdges (boundary : holes))
+      missingEdges = filter (not . (`Store.member` store)) edges
+      forced = forceEdges store missingEdges edges
+      Triangulation _ repaired = repairDegeneracies edges (Triangulation (hull triangulation) forced)
+      constrained = Store.triangles repaired
+  pure $ List.foldl' (flip trianglesOutside) (trianglesInside boundary constrained) holes
+  where
+    allPoints = HS.toList . HS.fromList $ points ++ concatMap vertices (boundary : holes)
+
+-- | The chain of edges the polygon edge becomes when the points lying on it
+-- (strictly between its endpoints) are made vertices: an edge with a vertex
+-- on it cannot exist in a triangulation, so the constraint is the chain.
+subdivide :: [Point] -> Edge -> [Edge]
+subdivide points (Edge a b) = zipWith mkEdge chain (drop 1 chain)
+  where
+    chain = a : List.sortOn (manhattanDistance a) (filter (liesBetween a b) points) ++ [b]
+    liesBetween u v p =
+      p /= u
+        && p /= v
+        && orientation u v p == Collinear
+        && min (px u) (px v) <= px p
+        && px p <= max (px u) (px v)
+        && min (py u) (py v) <= py p
+        && py p <= max (py u) (py v)
+
+-- | Force the given edges into the triangulation: remove every triangle
+-- crossed by an edge and re-triangulate the two resulting pockets. The
+-- restricted edges are never flipped away afterwards.
+forceEdges :: Store -> [Edge] -> [Edge] -> Store
+forceEdges store [] _ = store
+forceEdges store (edge@(Edge p1 p2) : edges) restrictedEdges =
+  let crossedEdges = [(e, p) | e <- Store.edges store, Just p <- [intersection edge e]]
+      edgesPoints = concatMap (\(Edge a b, _) -> [a, b]) (List.sortOn snd crossedEdges)
+      pointsOn side = dedupeConsecutive $ filter (\p -> turn p1 p2 p == side) edgesPoints
+      pockets = [p1 : p2 : reverse (pointsOn Clockwise), p2 : p1 : pointsOn CounterClockwise]
+      deletingTriangles = HS.toList $ HS.fromList (concatMap (\(e, _) -> Store.trianglesOn e store) crossedEdges)
+      store' = List.foldl' (flip Store.delete) store deletingTriangles
+      store'' = List.foldl' (\acc points -> fillPocket acc points restrictedEdges) store' pockets
+   in forceEdges store'' edges restrictedEdges
+
+dedupeConsecutive :: Eq a => [a] -> [a]
+dedupeConsecutive = concatMap (take 1) . List.group
+
+fillPocket :: Store -> [Point] -> [Edge] -> Store
+fillPocket store (p1 : p2 : p3 : ps) restrictedEdges
+  | isValidCandidate p1 p2 p3 ps =
+      let triangle = mkTriangle p1 p2 p3
+          store' = legalize (Store.insert triangle store) (triangleEdges triangle) restrictedEdges
+       in fillPocket store' (p1 : p3 : ps) restrictedEdges
+  | otherwise = fillPocket store (p2 : p3 : ps ++ [p1]) restrictedEdges
+fillPocket store _ _ = store
diff --git a/src/Triangulation/Flip.hs b/src/Triangulation/Flip.hs
new file mode 100644
--- /dev/null
+++ b/src/Triangulation/Flip.hs
@@ -0,0 +1,77 @@
+-- | Delaunay edge flips: restoring the empty-circumcircle property locally by
+-- swapping the diagonal of a convex quadrilateral.
+module Triangulation.Flip (
+  trianglesOnQuadrilateral,
+  legalize,
+) where
+
+import Data.List qualified as List
+import Triangulation.Geometry.Edge (Edge (..))
+import Triangulation.Geometry.Point (Orientation (..), Point, orientation)
+import Triangulation.Geometry.Polygon (Hull4 (..), hullOf4)
+import Triangulation.Geometry.Triangle (
+  Triangle (..),
+  isOutsideCircumcircle,
+  mkTriangle,
+  triangleEdges,
+ )
+import Triangulation.Store (Store)
+import Triangulation.Store qualified as Store
+
+-- | The two triangles of a convex quadrilateral (vertices in boundary order),
+-- split along the diagonal that satisfies the Delaunay condition.
+--
+-- When three of the four vertices are collinear only one diagonal gives two
+-- triangles of non-zero area, and that one is chosen regardless of the
+-- in-circle test.
+trianglesOnQuadrilateral :: Point -> Point -> Point -> Point -> (Triangle, Triangle)
+trianglesOnQuadrilateral p1 p2 p3 p4
+  | hasFlat acrossP2P4 = acrossP1P3
+  | hasFlat acrossP1P3 = acrossP2P4
+  | isOutsideCircumcircle p1 (mkTriangle p2 p3 p4) = acrossP2P4
+  | otherwise = acrossP1P3
+  where
+    acrossP2P4 = (mkTriangle p1 p2 p4, mkTriangle p2 p3 p4)
+    acrossP1P3 = (mkTriangle p1 p2 p3, mkTriangle p1 p3 p4)
+    hasFlat (t1, t2) = isFlat t1 || isFlat t2
+
+-- | Restore the Delaunay condition around the given edges by flipping the
+-- diagonal of every non-Delaunay pair of adjacent triangles, propagating to
+-- the edges of the new triangles. Restricted edges are never flipped.
+legalize :: Store -> [Edge] -> [Edge] -> Store
+legalize store [] _ = store
+legalize store (edge : es) restrictedEdges
+  | edge `elem` restrictedEdges = legalize store es restrictedEdges
+  | otherwise = case Store.trianglesOn edge store of
+      [tr1, tr2]
+        | Edge a b <- edge
+        , c <- apex edge tr1
+        , d <- apex edge tr2
+        , Quadrilateral p1 p2 p3 p4 <- hullOf4 a b c d
+        , (new1, new2) <- trianglesOnQuadrilateral p1 p2 p3 p4
+        , new1 `notElem` [tr1, tr2]
+        , not (isFlat new1 || isFlat new2) ->
+            let newEdges = dedupe (triangleEdges new1 ++ triangleEdges new2) ++ es
+                store' = List.foldl' (flip Store.delete) store [tr1, tr2]
+                store'' = List.foldl' (flip Store.insert) store' [new1, new2]
+             in legalize store'' newEdges restrictedEdges
+      _ -> legalize store es restrictedEdges
+
+-- | The vertex of the triangle that is not an endpoint of the edge. The two
+-- endpoints and the two apexes are the four distinct points of the pair of
+-- triangles sharing the edge, which is what the flip needs; building a set to
+-- find them, as this used to, costs an allocation per flip.
+apex :: Edge -> Triangle -> Point
+apex (Edge u v) (Triangle p1 p2 p3)
+  | p1 /= u && p1 /= v = p1
+  | p2 /= u && p2 /= v = p2
+  | otherwise = p3
+
+-- | The distinct elements of a list of six edges at most; cheaper than a set
+-- for that size, and unlike the previous code it does not walk the whole
+-- pending queue on every flip.
+dedupe :: [Edge] -> [Edge]
+dedupe = foldr (\e rest -> e : filter (/= e) rest) []
+
+isFlat :: Triangle -> Bool
+isFlat (Triangle a b c) = orientation a b c == Collinear
diff --git a/src/Triangulation/Generator.hs b/src/Triangulation/Generator.hs
new file mode 100644
--- /dev/null
+++ b/src/Triangulation/Generator.hs
@@ -0,0 +1,53 @@
+{-# LANGUAGE DerivingStrategies #-}
+
+-- | Deterministic pseudo-random point sets.
+module Triangulation.Generator (
+  Rectangle (..),
+  generatePoints,
+  generatePointsWithDistance,
+) where
+
+import System.Random (mkStdGen, randoms)
+import Triangulation.Geometry.Point (Point (..), manhattanDistance)
+
+-- | An axis-aligned rectangle given by its two extreme corners.
+data Rectangle = Rectangle
+  { minCorner :: !Point
+  -- ^ smallest x and y (lower-left with the y axis up)
+  , maxCorner :: !Point
+  -- ^ largest x and y (upper-right with the y axis up)
+  }
+  deriving stock (Eq, Show)
+
+-- | @n@ pseudo-random points in the rectangle, determined by the seed.
+generatePoints :: Int -> Int -> Rectangle -> [Point]
+generatePoints seed n r = take n $ randomPointsInRectangle seed r
+
+-- | @n@ random points in the rectangle, added to the existing points so that
+-- all points are at least @distance@ apart (in 'manhattanDistance').
+generatePointsWithDistance :: Int -> Int -> Double -> Rectangle -> [Point] -> [Point]
+generatePointsWithDistance seed n distance r =
+  addPointsWithDistance n distance (randomPointsInRectangle seed r)
+
+-- | One stream of numbers, taken two at a time: @split@ is deprecated in
+-- @random-1.3@ and its replacement does not exist in @random-1.2@.
+randomPointsInRectangle :: Int -> Rectangle -> [Point]
+randomPointsInRectangle seed r = pairUp (randoms (mkStdGen seed))
+  where
+    pairUp (x : y : rest) = Point (scaleX r x) (scaleY r y) : pairUp rest
+    pairUp _ = []
+
+addPointsWithDistance :: Int -> Double -> [Point] -> [Point] -> [Point]
+addPointsWithDistance n distance = go n
+  where
+    go 0 _ out = out
+    go _ [] out = out
+    go k (c : cs) out
+      | all ((>= distance) . manhattanDistance c) out = go (k - 1) cs (c : out)
+      | otherwise = go k cs out
+
+scaleX :: Rectangle -> Double -> Double
+scaleX Rectangle {minCorner = Point x1 _, maxCorner = Point x2 _} x = x1 + x * (x2 - x1)
+
+scaleY :: Rectangle -> Double -> Double
+scaleY Rectangle {minCorner = Point _ y1, maxCorner = Point _ y2} y = y1 + y * (y2 - y1)
diff --git a/src/Triangulation/Geometry.hs b/src/Triangulation/Geometry.hs
new file mode 100644
--- /dev/null
+++ b/src/Triangulation/Geometry.hs
@@ -0,0 +1,16 @@
+-- | Geometric primitives of the library, re-exported from one place.
+module Triangulation.Geometry (
+  module Triangulation.Geometry.Point,
+  module Triangulation.Geometry.Edge,
+  module Triangulation.Geometry.Triangle,
+  module Triangulation.Geometry.Polygon,
+  module Triangulation.Geometry.Ring,
+  module Triangulation.Geometry.Exact,
+) where
+
+import Triangulation.Geometry.Edge
+import Triangulation.Geometry.Exact
+import Triangulation.Geometry.Point
+import Triangulation.Geometry.Polygon
+import Triangulation.Geometry.Ring
+import Triangulation.Geometry.Triangle
diff --git a/src/Triangulation/Geometry/Edge.hs b/src/Triangulation/Geometry/Edge.hs
new file mode 100644
--- /dev/null
+++ b/src/Triangulation/Geometry/Edge.hs
@@ -0,0 +1,66 @@
+{-# LANGUAGE DeriveAnyClass #-}
+{-# LANGUAGE DerivingStrategies #-}
+{-# LANGUAGE PatternSynonyms #-}
+
+-- | Undirected segments between two points.
+--
+-- The endpoints are stored in ascending order, so the derived 'Eq', 'Ord' and
+-- 'Hashable' instances are lawful and @mkEdge a b == mkEdge b a@. Build edges
+-- with 'mkEdge'; take them apart with the read-only v'Edge' pattern.
+module Triangulation.Geometry.Edge (
+  Edge (Edge),
+  mkEdge,
+  intersection,
+) where
+
+import Control.DeepSeq (NFData)
+import Data.Hashable (Hashable, hashWithSalt)
+import GHC.Generics (Generic)
+import Triangulation.Geometry.Point (Orientation (..), Point (..), orientation)
+
+-- | An undirected segment; see the module header.
+data Edge = UnsafeEdge {-# UNPACK #-} !Point {-# UNPACK #-} !Point
+  deriving stock (Eq, Ord, Generic)
+  deriving anyclass (NFData)
+
+-- | Written out rather than derived through 'Generic', which dominated the
+-- profile. The endpoints are kept sorted, so equal edges hash equally.
+instance Hashable Edge where
+  hashWithSalt salt (UnsafeEdge a b) = salt `hashWithSalt` a `hashWithSalt` b
+  {-# INLINE hashWithSalt #-}
+
+instance Show Edge where
+  showsPrec d (UnsafeEdge a b) =
+    showParen (d > 10) $ showString "Edge " . showsPrec 11 a . showChar ' ' . showsPrec 11 b
+
+-- | Smart constructor: orders the endpoints.
+mkEdge :: Point -> Point -> Edge
+mkEdge a b
+  | a <= b = UnsafeEdge a b
+  | otherwise = UnsafeEdge b a
+
+-- | The endpoints, in ascending order.
+pattern Edge :: Point -> Point -> Edge
+pattern Edge a b <- UnsafeEdge a b
+
+{-# COMPLETE Edge #-}
+
+-- | The point where two edges cross, if they cross strictly inside both of
+-- them: each edge has the other's endpoints strictly on opposite sides.
+-- Edges sharing an endpoint, or merely touching one, never count as
+-- intersecting. The decision is exact; the returned point is computed in
+-- floating point.
+intersection :: Edge -> Edge -> Maybe Point
+intersection (Edge p1@(Point x1 y1) p2@(Point x2 y2)) (Edge p3@(Point x3 y3) p4@(Point x4 y4))
+  | haveCommonPoint = Nothing
+  | separates p1 p2 p3 p4 && separates p3 p4 p1 p2 = Just crossing
+  | otherwise = Nothing
+  where
+    haveCommonPoint = p1 == p3 || p1 == p4 || p2 == p3 || p2 == p4
+    separates a b c d = case (orientation a b c, orientation a b d) of
+      (Clockwise, CounterClockwise) -> True
+      (CounterClockwise, Clockwise) -> True
+      _ -> False
+    det = (x2 - x1) * (y4 - y3) - (y2 - y1) * (x4 - x3) -- non-zero: the edges are not parallel
+    t = ((x3 - x1) * (y4 - y3) - (y3 - y1) * (x4 - x3)) / det
+    crossing = Point (x1 + t * (x2 - x1)) (y1 + t * (y2 - y1))
diff --git a/src/Triangulation/Geometry/Exact.hs b/src/Triangulation/Geometry/Exact.hs
new file mode 100644
--- /dev/null
+++ b/src/Triangulation/Geometry/Exact.hs
@@ -0,0 +1,87 @@
+-- | Exact geometric predicates on 'Double' coordinates.
+--
+-- Each predicate is the sign of a determinant. The determinant is first
+-- evaluated in floating point together with a bound on its rounding error
+-- (the static filters of Shewchuk's /Adaptive Precision Floating-Point
+-- Arithmetic and Fast Robust Geometric Predicates/); when the magnitude of the
+-- result exceeds the bound, its sign is certain. Otherwise the determinant is
+-- recomputed in 'Rational' arithmetic, which is exact because every 'Double'
+-- converts to a 'Rational' without loss. The exact branch is a lazy thunk and
+-- is only evaluated in the rare near-degenerate cases.
+--
+-- The predicates take raw coordinates so that this module sits below the
+-- geometric types.
+module Triangulation.Geometry.Exact (
+  Coordinates,
+  orientationSign,
+  inCircleSign,
+) where
+
+-- | A point as an @(x, y)@ pair.
+type Coordinates = (Double, Double)
+
+-- | Half the machine epsilon of 'Double': the relative rounding error of one operation.
+epsilon :: Double
+epsilon = encodeFloat 1 (-53)
+
+-- | Sign of @(a - c) × (b - c)@: 'GT' when @a -> b -> c@ turns counter-clockwise
+-- (y axis up), 'LT' when clockwise, 'EQ' when the points are collinear.
+--
+-- The exact fallback is a separate function rather than an argument, so that
+-- the common case does not build a 'Rational' thunk it never looks at.
+orientationSign :: Coordinates -> Coordinates -> Coordinates -> Ordering
+orientationSign a@(ax, ay) b@(bx, by) c@(cx, cy)
+  | approximate > errorBound = GT
+  | approximate < negate errorBound = LT
+  | otherwise = exactOrientation a b c
+  where
+    detLeft = (ax - cx) * (by - cy)
+    detRight = (ay - cy) * (bx - cx)
+    approximate = detLeft - detRight
+    errorBound = (3 + 16 * epsilon) * epsilon * (abs detLeft + abs detRight)
+
+-- | The orientation determinant in exact arithmetic. Every 'Double' converts
+-- to a 'Rational' without loss, so the sign is the true one.
+exactOrientation :: Coordinates -> Coordinates -> Coordinates -> Ordering
+{-# NOINLINE exactOrientation #-}
+exactOrientation (ax, ay) (bx, by) (cx, cy) =
+  compare ((r ax - r cx) * (r by - r cy) - (r ay - r cy) * (r bx - r cx)) 0
+  where
+    r = toRational
+
+-- | Sign of the in-circle determinant: for @a@, @b@, @c@ in counter-clockwise
+-- order, 'GT' when @d@ lies strictly inside their circumcircle, 'LT' when
+-- strictly outside, 'EQ' when on it. The sign flips for a clockwise triple.
+inCircleSign :: Coordinates -> Coordinates -> Coordinates -> Coordinates -> Ordering
+inCircleSign a@(ax, ay) b@(bx, by) c@(cx, cy) d@(dx, dy)
+  | det > errorBound = GT
+  | det < negate errorBound = LT
+  | otherwise = exactInCircle a b c d
+  where
+    (adx, ady) = (ax - dx, ay - dy)
+    (bdx, bdy) = (bx - dx, by - dy)
+    (cdx, cdy) = (cx - dx, cy - dy)
+    (bdxcdy, cdxbdy, alift) = (bdx * cdy, cdx * bdy, adx * adx + ady * ady)
+    (cdxady, adxcdy, blift) = (cdx * ady, adx * cdy, bdx * bdx + bdy * bdy)
+    (adxbdy, bdxady, clift) = (adx * bdy, bdx * ady, cdx * cdx + cdy * cdy)
+    det = alift * (bdxcdy - cdxbdy) + blift * (cdxady - adxcdy) + clift * (adxbdy - bdxady)
+    permanent =
+      (abs bdxcdy + abs cdxbdy) * alift
+        + (abs cdxady + abs adxcdy) * blift
+        + (abs adxbdy + abs bdxady) * clift
+    errorBound = (10 + 96 * epsilon) * epsilon * permanent
+
+-- | The in-circle determinant in exact arithmetic; see 'exactOrientation'.
+exactInCircle :: Coordinates -> Coordinates -> Coordinates -> Coordinates -> Ordering
+{-# NOINLINE exactInCircle #-}
+exactInCircle (ax, ay) (bx, by) (cx, cy) (dx, dy) = compare determinant 0
+  where
+    (ax', ay', bx', by', cx', cy', dx', dy') = (r ax, r ay, r bx, r by, r cx, r cy, r dx, r dy)
+    (adx, ady) = (ax' - dx', ay' - dy')
+    (bdx, bdy) = (bx' - dx', by' - dy')
+    (cdx, cdy) = (cx' - dx', cy' - dy')
+    determinant =
+      (adx * adx + ady * ady) * (bdx * cdy - cdx * bdy)
+        + (bdx * bdx + bdy * bdy) * (cdx * ady - adx * cdy)
+        + (cdx * cdx + cdy * cdy) * (adx * bdy - bdx * ady)
+    r = toRational
diff --git a/src/Triangulation/Geometry/Point.hs b/src/Triangulation/Geometry/Point.hs
new file mode 100644
--- /dev/null
+++ b/src/Triangulation/Geometry/Point.hs
@@ -0,0 +1,154 @@
+{-# LANGUAGE DeriveAnyClass #-}
+{-# LANGUAGE DerivingStrategies #-}
+
+-- | Points in the Euclidean plane and predicates on them: orientation of a
+-- triple, angles, distances, extreme points of a set.
+--
+-- The library uses the mathematical convention: the y axis points up.
+-- Rendering code flips it when mapping to pixel rows.
+module Triangulation.Geometry.Point (
+  Point (..),
+  coordinates,
+  Orientation (..),
+  manhattanDistance,
+  orientation,
+  turn,
+  cosSquaredAngle,
+  bottomRight,
+  leftTop,
+  topLeft,
+  rightBottom,
+) where
+
+import Control.DeepSeq (NFData)
+import Data.Bits (shiftR, xor)
+import Data.Foldable (minimumBy)
+import Data.Hashable (Hashable, hashWithSalt)
+import Data.List.NonEmpty (NonEmpty)
+import Data.Ord (Down (..), comparing)
+import Data.Word (Word64)
+import GHC.Float (castDoubleToWord64)
+import GHC.Generics (Generic)
+import Triangulation.Geometry.Exact (orientationSign)
+
+-- | A point of the plane.
+data Point = Point {px :: !Double, py :: !Double}
+  deriving stock (Eq, Ord, Show, Generic)
+  deriving anyclass (NFData)
+
+-- | A hash of the two coordinates, written out rather than derived.
+--
+-- 'Data.Hashable's own instance for 'Double' runs a strong mixing function per
+-- field, and a profile taken after optimisation showed that mixing to be about
+-- half the running time of a triangulation: every lookup in the triangle store
+-- hashes points. The geometry only needs a hash good enough to spread keys
+-- across a 'Data.HashMap.Strict.HashMap', so the two bit patterns are combined
+-- and put through the MurmurHash3 finaliser once. Equal points have equal
+-- coordinates, so the instance is lawful.
+instance Hashable Point where
+  hashWithSalt salt (Point x y) =
+    fromIntegral (finalise (castDoubleToWord64 x * goldenGamma `xor` castDoubleToWord64 y))
+      `xor` salt
+  {-# INLINE hashWithSalt #-}
+
+-- | The odd multiplier of the golden ratio, as used by @splitmix@.
+goldenGamma :: Word64
+goldenGamma = 0x9E3779B97F4A7C15
+
+-- | The finaliser of MurmurHash3: two multiplications, three shifts, three
+-- exclusive ors, and no memory traffic.
+finalise :: Word64 -> Word64
+finalise w0 =
+  let w1 = (w0 `xor` (w0 `shiftR` 33)) * 0xFF51AFD7ED558CCD
+      w2 = (w1 `xor` (w1 `shiftR` 33)) * 0xC4CEB9FE1A85EC53
+   in w2 `xor` (w2 `shiftR` 33)
+{-# INLINE finalise #-}
+
+-- | The direction of the turn made at the middle point of an ordered triple.
+data Orientation = Clockwise | CounterClockwise | Collinear
+  deriving stock (Eq, Show)
+
+-- | The point as an @(x, y)@ pair, for the exact predicates.
+coordinates :: Point -> (Double, Double)
+coordinates (Point x y) = (x, y)
+
+-- | Manhattan (L1) distance.
+manhattanDistance :: Point -> Point -> Double
+manhattanDistance (Point x1 y1) (Point x2 y2) = abs (x2 - x1) + abs (y2 - y1)
+
+-- | Direction of the turn @a -> b -> c@, with the y axis pointing up.
+-- Exact: see "Triangulation.Geometry.Exact".
+orientation :: Point -> Point -> Point -> Orientation
+orientation a b c = case orientationSign (coordinates a) (coordinates b) (coordinates c) of
+  LT -> Clockwise
+  GT -> CounterClockwise
+  EQ -> Collinear
+
+-- | The orientation of the /symbolically perturbed/ points: 'Collinear' only
+-- when two of the points coincide.
+--
+-- This is Edelsbrunner and Mücke's Simulation of Simplicity. Every point is
+-- imagined displaced by an infinitesimal amount that depends only on its rank
+-- (its position in the 'Ord' order): a lower rank gets a larger displacement,
+-- and the y coordinate a larger one than x. Ties in 'orientation' are then
+-- broken by the first non-zero term of the perturbed determinant, which for
+-- three points ranked @i < j < k@ is, in decreasing significance,
+-- @x_k - x_j@, @y_j - y_k@, @x_i - x_k@ and finally a constant. Because the
+-- displacement is a fixed function of the points, every decision the
+-- algorithm makes is consistent with one and the same perturbed point set,
+-- which is what lets the convex-hull and merge code assume general position.
+turn :: Point -> Point -> Point -> Orientation
+turn a b c = case orientation a b c of
+  Collinear
+    | a == b || b == c || a == c -> Collinear
+    | otherwise -> perturbed
+  o -> o
+  where
+    perturbed =
+      let (evenPermutation, Point xi _, Point xj yj, Point xk yk) = rankSorted a b c
+          firstNonZero = case filter (/= 0) [xk - xj, yj - yk, xi - xk] of
+            t : _ -> compare t 0
+            [] -> GT
+          sign = if evenPermutation then firstNonZero else flipOrdering firstNonZero
+       in case sign of
+            GT -> CounterClockwise
+            LT -> Clockwise
+            EQ -> Collinear
+    flipOrdering LT = GT
+    flipOrdering GT = LT
+    flipOrdering EQ = EQ
+
+-- | The three points in increasing 'Ord' order, and whether that reordering is
+-- an even permutation of the arguments.
+rankSorted :: Point -> Point -> Point -> (Bool, Point, Point, Point)
+rankSorted a b c
+  | a <= b && b <= c = (True, a, b, c)
+  | a <= c && c <= b = (False, a, c, b)
+  | b <= a && a <= c = (False, b, a, c)
+  | b <= c && c <= a = (True, b, c, a)
+  | c <= a && a <= b = (True, c, a, b)
+  | otherwise = (False, c, b, a)
+
+-- | Squared cosine of the angle at @o@ between the rays @o -> a@ and @o -> b@.
+cosSquaredAngle :: Point -> Point -> Point -> Double
+cosSquaredAngle (Point x0 y0) (Point x1 y1) (Point x2 y2) =
+  let (dx1, dy1) = (x1 - x0, y1 - y0)
+      (dx2, dy2) = (x2 - x0, y2 - y0)
+      dot = dx1 * dx2 + dy1 * dy2
+   in dot * dot / ((dx1 * dx1 + dy1 * dy1) * (dx2 * dx2 + dy2 * dy2))
+
+-- | The rightmost of the bottom points.
+bottomRight :: NonEmpty Point -> Point
+bottomRight = minimumBy (comparing (\(Point x y) -> (y, Down x)))
+
+-- | The topmost of the left points.
+leftTop :: NonEmpty Point -> Point
+leftTop = minimumBy (comparing (\(Point x y) -> (x, Down y)))
+
+-- | The leftmost of the top points.
+topLeft :: NonEmpty Point -> Point
+topLeft = minimumBy (comparing (\(Point x y) -> (Down y, x)))
+
+-- | The bottommost of the right points.
+rightBottom :: NonEmpty Point -> Point
+rightBottom = minimumBy (comparing (\(Point x y) -> (Down x, y)))
diff --git a/src/Triangulation/Geometry/Polygon.hs b/src/Triangulation/Geometry/Polygon.hs
new file mode 100644
--- /dev/null
+++ b/src/Triangulation/Geometry/Polygon.hs
@@ -0,0 +1,182 @@
+{-# LANGUAGE DerivingStrategies #-}
+
+-- | Simple polygons, convex hulls of three and four points, and the polygon
+-- side of merging two triangulations: sliding a bridge between two hulls
+-- until it is tangent to both.
+--
+-- Throughout the library polygons are wound /clockwise/ (with the y axis
+-- pointing up): walking along the boundary, the interior is on the right.
+module Triangulation.Geometry.Polygon (
+  Polygon (..),
+  vertices,
+  polygonEdges,
+  Bridge,
+  Hull4 (..),
+  hullOf3,
+  hullOf4,
+  hull4Polygon,
+  tangents,
+  isConvex,
+  isPointInPolygon,
+) where
+
+import Control.DeepSeq (NFData)
+import Data.List.NonEmpty (NonEmpty (..))
+import Data.List.NonEmpty qualified as NE
+import Triangulation.Geometry.Edge (Edge, mkEdge)
+import Triangulation.Geometry.Point (Orientation (..), Point (..), turn)
+import Triangulation.Geometry.Ring (
+  arc,
+  cyclicPairs,
+  cyclicTriples,
+  predecessor,
+  splitLoop,
+  successor,
+ )
+
+-- | A simple polygon, wound clockwise; see the module header.
+newtype Polygon = Polygon (NonEmpty Point)
+  deriving stock (Eq, Show)
+  deriving newtype (NFData)
+
+-- | The vertices, in boundary order.
+vertices :: Polygon -> [Point]
+vertices (Polygon points) = NE.toList points
+
+-- | The boundary edges, including the one closing the ring.
+polygonEdges :: Polygon -> [Edge]
+polygonEdges (Polygon points) = map (uncurry mkEdge) (cyclicPairs points)
+
+-- hulls of 3 and 4 points
+
+-- | Convex hull of three points, wound clockwise. A collinear triple is kept
+-- as a degenerate polygon.
+hullOf3 :: Point -> Point -> Point -> Polygon
+hullOf3 a b c =
+  let (pivot, r1, r2) = sort3By compare a b c
+   in Polygon $ case hullOrder pivot r1 r2 of
+        GT -> pivot :| [r2, r1]
+        _ -> pivot :| [r1, r2]
+
+-- | Convex hull of four points.
+data Hull4
+  = -- | all four points are hull vertices, wound clockwise
+    Quadrilateral Point Point Point Point
+  | -- | a hull triangle, wound clockwise, and the point inside it
+    TriangleWithInner Point Point Point Point
+
+-- | Convex hull of four points, telling a convex quadrilateral from a
+-- triangle with a point inside.
+hullOf4 :: Point -> Point -> Point -> Point -> Hull4
+hullOf4 a b c d =
+  let (p1, r1, r2, r3) = sort4By compare a b c d
+      (p2, p3, p4) = sort3By (hullOrder p1) r1 r2 r3
+   in case (turn p2 p3 p4, turn p3 p4 p1) of
+        (CounterClockwise, _) -> TriangleWithInner p1 p2 p4 p3
+        (_, Clockwise) -> Quadrilateral p1 p2 p3 p4
+        _ -> TriangleWithInner p1 p2 p3 p4 -- rare case: points on the same line
+
+-- | The hull as a polygon, dropping an inner point.
+hull4Polygon :: Hull4 -> Polygon
+hull4Polygon (Quadrilateral p1 p2 p3 p4) = Polygon (p1 :| [p2, p3, p4])
+hull4Polygon (TriangleWithInner p1 p2 p3 _) = Polygon (p1 :| [p2, p3])
+
+-- | Angular order around a pivot: @a@ precedes @b@ when @b@ lies clockwise of @pivot -> a@.
+hullOrder :: Point -> Point -> Point -> Ordering
+hullOrder pivot a b = case turn pivot a b of
+  Clockwise -> LT
+  CounterClockwise -> GT
+  Collinear -> EQ -- unreachable: 'turn' resolves collinear triples
+
+sort3By :: (a -> a -> Ordering) -> a -> a -> a -> (a, a, a)
+sort3By cmp a b c
+  | cmp c lo /= GT = (c, lo, hi)
+  | cmp c hi /= GT = (lo, c, hi)
+  | otherwise = (lo, hi, c)
+  where
+    (lo, hi) = if cmp a b == GT then (b, a) else (a, b)
+
+sort4By :: (a -> a -> Ordering) -> a -> a -> a -> a -> (a, a, a, a)
+sort4By cmp a b c d
+  | cmp a x /= GT = (a, x, y, z)
+  | cmp a y /= GT = (x, a, y, z)
+  | cmp a z /= GT = (x, y, a, z)
+  | otherwise = (x, y, z, a)
+  where
+    (x, y, z) = sort3By cmp b c d
+
+-- tangents
+
+-- | A pair of vertices joining two polygons being merged:
+-- @(vertex of the first polygon, vertex of the second polygon)@.
+type Bridge = (Point, Point)
+
+-- | Starting from an initial bridge (used as both the bottom and the top
+-- one), slide the bridge endpoints along the two hulls until both bridges are
+-- tangent to both polygons. Returns the bottom and top tangents and the merged
+-- hull. 'Nothing' if a bridge endpoint is not a vertex of its polygon.
+tangents :: Polygon -> Polygon -> Bridge -> Maybe (Bridge, Bridge, Polygon)
+tangents (Polygon points1) (Polygon points2) start = go start start
+  where
+    go bottom@(bl, br) top@(tl, tr) = case slide bottom top of
+      Just (bottom', top') -> go bottom' top'
+      Nothing -> do
+        cut1 <- arc tl bl points1
+        cut2 <- arc br tr points2
+        pure (bottom, top, removeLoops (Polygon (cut1 <> cut2)) bottom top)
+
+    -- One step: move whichever bridge endpoint still sees its hull neighbour
+    -- on the wrong side. 'Nothing' when both bridges are tangent.
+    slide (bl, br) (tl, tr) = do
+      blPred <- predecessor bl points1
+      brSucc <- successor br points2
+      tlSucc <- successor tl points1
+      trPred <- predecessor tr points2
+      pick blPred brSucc tlSucc trPred
+      where
+        pick blPred brSucc tlSucc trPred
+          | turn blPred bl br == CounterClockwise = Just ((blPred, br), (tl, tr))
+          | turn bl br brSucc == CounterClockwise = Just ((bl, brSucc), (tl, tr))
+          | turn tr tl tlSucc == CounterClockwise = Just ((bl, br), (tlSucc, tr))
+          | turn trPred tr tl == CounterClockwise = Just ((bl, br), (tl, trPred))
+          | otherwise = Nothing
+
+-- | When a bridge degenerates to a single vertex on one side, the merged ring
+-- visits that vertex twice; keep the larger of the two loops.
+removeLoops :: Polygon -> Bridge -> Bridge -> Polygon
+removeLoops polygon (bl, br) (tl, tr)
+  | bl == tl = largerLoop polygon bl
+  | br == tr = largerLoop polygon br
+  | otherwise = polygon
+
+largerLoop :: Polygon -> Point -> Polygon
+largerLoop (Polygon points) point =
+  let (outer, inner) = splitLoop point points
+   in case NE.nonEmpty inner of
+        Just inner' | doubledArea (Polygon outer) <= doubledArea (Polygon inner') -> Polygon inner'
+        _ -> Polygon outer
+
+-- | Twice the area (shoelace formula); enough for comparisons.
+doubledArea :: Polygon -> Double
+doubledArea (Polygon points) =
+  abs $ sum [x1 * y2 - y1 * x2 | (Point x1 y1, Point x2 y2) <- cyclicPairs points]
+
+-- checking
+
+-- | Whether the polygon is convex (no counter-clockwise turn along the boundary).
+isConvex :: Polygon -> Bool
+isConvex (Polygon (_ :| [_])) = True
+isConvex (Polygon points) =
+  all (\(p1, p2, p3) -> turn p1 p2 p3 /= CounterClockwise) (cyclicTriples points)
+
+-- | Whether the point lies strictly inside the polygon, which need not be
+-- convex: the number of polygon edges crossed by a ray from the point to the
+-- right is odd (the even–odd rule). A vertex of the polygon does not count as
+-- inside; for a point exactly on an edge the answer is not specified.
+isPointInPolygon :: Polygon -> Point -> Bool
+isPointInPolygon (Polygon points) point@(Point x y) =
+  point `notElem` points && odd (length (filter crossesRay (cyclicPairs points)))
+  where
+    -- half-open in y, so that a ray through a vertex is counted once
+    crossesRay (Point x1 y1, Point x2 y2) =
+      (y1 > y) /= (y2 > y) && x < x1 + (y - y1) * (x2 - x1) / (y2 - y1)
diff --git a/src/Triangulation/Geometry/Ring.hs b/src/Triangulation/Geometry/Ring.hs
new file mode 100644
--- /dev/null
+++ b/src/Triangulation/Geometry/Ring.hs
@@ -0,0 +1,56 @@
+-- | A non-empty list viewed as a ring: the successor of the last element is the
+-- first one. Polygon vertex lists are rings.
+module Triangulation.Geometry.Ring (
+  cyclicPairs,
+  cyclicTriples,
+  successor,
+  predecessor,
+  arc,
+  splitLoop,
+) where
+
+import Data.List.NonEmpty (NonEmpty (..))
+import Data.List.NonEmpty qualified as NE
+import Data.Tuple (swap)
+
+-- | Every element paired with its successor.
+cyclicPairs :: NonEmpty a -> [(a, a)]
+cyclicPairs (x :| xs) = zip (x : xs) (xs ++ [x])
+
+-- | Every element with its neighbours: @(predecessor, element, successor)@.
+cyclicTriples :: NonEmpty a -> [(a, a, a)]
+cyclicTriples ring =
+  let xs = NE.toList ring
+      n = length xs
+   in zip3 (drop (n - 1) xs ++ take (n - 1) xs) xs (drop 1 xs ++ take 1 xs)
+
+-- | Successor of the first occurrence of an element.
+successor :: Eq a => a -> NonEmpty a -> Maybe a
+successor x = lookup x . cyclicPairs
+
+-- | Predecessor of the first occurrence of an element.
+predecessor :: Eq a => a -> NonEmpty a -> Maybe a
+predecessor x = lookup x . map swap . cyclicPairs
+
+-- | The arc from @from@ to @to@, walking forward and wrapping around, both ends
+-- included. When @from == to@ the arc goes all the way round and ends at
+-- @from@ again; the same happens when @to@ is absent. 'Nothing' when @from@
+-- is not on the ring.
+arc :: Eq a => a -> a -> NonEmpty a -> Maybe (NonEmpty a)
+arc from to ring = case break (== from) (NE.toList ring) of
+  (_, []) -> Nothing
+  (before, x : after) -> Just (x :| takeThrough (== to) (after ++ before ++ [x]))
+
+takeThrough :: (a -> Bool) -> [a] -> [a]
+takeThrough p xs = let (prefix, rest) = break p xs in prefix ++ take 1 rest
+
+-- | Split a ring that visits @x@ twice into the outer ring, which visits @x@
+-- once, and the inner loop between the two visits, which starts with @x@.
+-- A ring visiting @x@ at most once is returned unchanged with an empty loop.
+splitLoop :: Eq a => a -> NonEmpty a -> (NonEmpty a, [a])
+splitLoop x ring = case break (== x) (NE.toList ring) of
+  (_, []) -> (ring, [])
+  (before, _ : rest) ->
+    case break (== x) (reverse rest) of
+      (_, []) -> (ring, [])
+      (afterRev, _ : innerRev) -> (NE.prependList before (x :| reverse afterRev), x : reverse innerRev)
diff --git a/src/Triangulation/Geometry/Triangle.hs b/src/Triangulation/Geometry/Triangle.hs
new file mode 100644
--- /dev/null
+++ b/src/Triangulation/Geometry/Triangle.hs
@@ -0,0 +1,150 @@
+{-# LANGUAGE DeriveAnyClass #-}
+{-# LANGUAGE DerivingStrategies #-}
+{-# LANGUAGE PatternSynonyms #-}
+
+-- | Triangles and predicates on them: the Delaunay in-circle test, point
+-- containment, position relative to a polygon.
+--
+-- Vertices are stored in ascending order, so two triangles on the same three
+-- points are equal regardless of the order they were built from, and the
+-- derived 'Eq', 'Ord' and 'Hashable' instances are lawful. Build triangles
+-- with 'mkTriangle'; take them apart with the read-only v'Triangle' pattern.
+module Triangulation.Geometry.Triangle (
+  Triangle (Triangle),
+  mkTriangle,
+  triangleEdges,
+  triangleArea,
+  smallestAngle,
+  circumcenter,
+  isOutsideCircumcircle,
+  isPointInTriangle,
+  isValidCandidate,
+  isTriangleInPolygon,
+  trianglesInside,
+  trianglesOutside,
+) where
+
+import Control.DeepSeq (NFData)
+import Data.Hashable (Hashable, hashWithSalt)
+import Data.List (sort)
+import GHC.Generics (Generic)
+import Triangulation.Geometry.Edge (Edge, mkEdge)
+import Triangulation.Geometry.Exact (inCircleSign)
+import Triangulation.Geometry.Point (Orientation (..), Point (..), coordinates, orientation, turn)
+import Triangulation.Geometry.Polygon (Polygon, isPointInPolygon)
+
+-- | A triangle; see the module header.
+data Triangle = UnsafeTriangle {-# UNPACK #-} !Point {-# UNPACK #-} !Point {-# UNPACK #-} !Point
+  deriving stock (Eq, Ord, Generic)
+  deriving anyclass (NFData)
+
+-- | Written out rather than derived through 'Generic', which dominated the
+-- profile. The vertices are kept sorted, so equal triangles hash equally.
+instance Hashable Triangle where
+  hashWithSalt salt (UnsafeTriangle a b c) =
+    salt `hashWithSalt` a `hashWithSalt` b `hashWithSalt` c
+  {-# INLINE hashWithSalt #-}
+
+instance Show Triangle where
+  showsPrec d (UnsafeTriangle a b c) =
+    showParen (d > 10) $
+      showString "Triangle "
+        . showsPrec 11 a
+        . showChar ' '
+        . showsPrec 11 b
+        . showChar ' '
+        . showsPrec 11 c
+
+-- | Smart constructor: sorts the vertices.
+mkTriangle :: Point -> Point -> Point -> Triangle
+mkTriangle a b c
+  | c <= lo = UnsafeTriangle c lo hi
+  | c <= hi = UnsafeTriangle lo c hi
+  | otherwise = UnsafeTriangle lo hi c
+  where
+    lo = min a b
+    hi = max a b
+
+-- | The vertices, in ascending order.
+pattern Triangle :: Point -> Point -> Point -> Triangle
+pattern Triangle a b c <- UnsafeTriangle a b c
+
+{-# COMPLETE Triangle #-}
+
+-- | The three edges.
+triangleEdges :: Triangle -> [Edge]
+triangleEdges (Triangle p1 p2 p3) = [mkEdge p1 p2, mkEdge p2 p3, mkEdge p3 p1]
+
+-- | The (unsigned) area.
+triangleArea :: Triangle -> Double
+triangleArea (Triangle (Point x1 y1) (Point x2 y2) (Point x3 y3)) =
+  abs ((x2 - x1) * (y3 - y1) - (x3 - x1) * (y2 - y1)) / 2
+
+-- | The smallest of the three angles, in degrees; 0 for a flat triangle.
+smallestAngle :: Triangle -> Double
+smallestAngle (Triangle a b c)
+  | orientation a b c == Collinear = 0
+  | otherwise = case sort [squaredDistance a b, squaredDistance b c, squaredDistance c a] of
+      -- the smallest angle is opposite the shortest side; law of cosines
+      [s1, s2, s3] -> acos (max (-1) (min 1 ((s2 + s3 - s1) / (2 * sqrt (s2 * s3))))) * 180 / pi
+      _ -> 0
+
+squaredDistance :: Point -> Point -> Double
+squaredDistance (Point x1 y1) (Point x2 y2) = (x2 - x1) * (x2 - x1) + (y2 - y1) * (y2 - y1)
+
+-- | The centre of the circumcircle; 'Nothing' for a flat triangle. Computed
+-- in floating point, relative to the first vertex.
+circumcenter :: Triangle -> Maybe Point
+circumcenter (Triangle a@(Point ax ay) b@(Point bx by) c@(Point cx cy))
+  | orientation a b c == Collinear = Nothing
+  | otherwise = Just (Point (ax + ux) (ay + uy))
+  where
+    (bx', by') = (bx - ax, by - ay)
+    (cx', cy') = (cx - ax, cy - ay)
+    d = 2 * (bx' * cy' - by' * cx')
+    (b2, c2) = (bx' * bx' + by' * by', cx' * cx' + cy' * cy')
+    ux = (cy' * b2 - by' * c2) / d
+    uy = (bx' * c2 - cx' * b2) / d
+
+-- | Whether the point lies outside of, or exactly on, the circumcircle of the
+-- triangle: the Delaunay in-circle test. Exact (see
+-- "Triangulation.Geometry.Exact"); a degenerate triangle has no circumcircle
+-- and every point counts as outside.
+isOutsideCircumcircle :: Point -> Triangle -> Bool
+isOutsideCircumcircle p (Triangle a b c) = case orientation a b c of
+  CounterClockwise -> inCircle /= GT
+  Clockwise -> inCircle /= LT
+  Collinear -> True
+  where
+    inCircle = inCircleSign (coordinates a) (coordinates b) (coordinates c) (coordinates p)
+
+-- | Whether the point lies inside the triangle, in the symbolically perturbed
+-- sense of 'turn': a point exactly on an edge is consistently assigned to one
+-- side of it. Vertices are excluded.
+isPointInTriangle :: Triangle -> Point -> Bool
+isPointInTriangle (Triangle a b c) p =
+  p `notElem` [a, b, c] && not (Clockwise `elem` turns && CounterClockwise `elem` turns)
+  where
+    turns = [turn a b p, turn b c p, turn c a p]
+
+-- | A candidate triangle @p1 p2 p3@ is accepted when the walk @p1 -> p2 -> p3@
+-- turns clockwise in the perturbed sense of 'turn' (the merge front advances
+-- on that side) and no other point lies inside it. A collinear triple passes
+-- as a zero-area triangle of the perturbed points; "Triangulation.Repair"
+-- removes those afterwards.
+isValidCandidate :: Point -> Point -> Point -> [Point] -> Bool
+isValidCandidate p1 p2 p3 points =
+  turn p1 p2 p3 == Clockwise && not (any (isPointInTriangle (mkTriangle p1 p2 p3)) points)
+
+-- | Whether the centroid of the triangle lies inside the polygon.
+isTriangleInPolygon :: Triangle -> Polygon -> Bool
+isTriangleInPolygon (Triangle (Point x1 y1) (Point x2 y2) (Point x3 y3)) polygon =
+  isPointInPolygon polygon (Point ((x1 + x2 + x3) / 3) ((y1 + y2 + y3) / 3))
+
+-- | The triangles whose centroid lies inside the polygon.
+trianglesInside :: Polygon -> [Triangle] -> [Triangle]
+trianglesInside polygon = filter (`isTriangleInPolygon` polygon)
+
+-- | The triangles whose centroid lies outside the polygon.
+trianglesOutside :: Polygon -> [Triangle] -> [Triangle]
+trianglesOutside polygon = filter (not . (`isTriangleInPolygon` polygon))
diff --git a/src/Triangulation/Leaf.hs b/src/Triangulation/Leaf.hs
new file mode 100644
--- /dev/null
+++ b/src/Triangulation/Leaf.hs
@@ -0,0 +1,31 @@
+-- | Triangulations of the smallest point sets, the leaves of the
+-- divide-and-conquer tree.
+module Triangulation.Leaf (
+  triangulateLeaf,
+) where
+
+import Data.List.NonEmpty (NonEmpty (..))
+import Triangulation.Flip (trianglesOnQuadrilateral)
+import Triangulation.Geometry.Point (Point)
+import Triangulation.Geometry.Polygon (Hull4 (..), Polygon (..), hull4Polygon, hullOf3, hullOf4)
+import Triangulation.Geometry.Triangle (Triangle, mkTriangle)
+import Triangulation.Types (Triangulation, fromTriangles)
+
+-- | Triangulation of 3 or 4 points; 'Nothing' for any other number.
+--
+-- Two points give a degenerate triangulation: a two-vertex hull and no
+-- triangles. It only arises when a five-point set is split 2 + 3, and the
+-- merge fills the gap between the segment and the triangle.
+triangulateLeaf :: [Point] -> Maybe Triangulation
+triangulateLeaf [a, b] = Just $ fromTriangles (Polygon (a :| [b])) []
+triangulateLeaf [a, b, c] = Just $ fromTriangles (hullOf3 a b c) [mkTriangle a b c]
+triangulateLeaf [a, b, c, d] =
+  let hull4 = hullOf4 a b c d
+   in Just $ fromTriangles (hull4Polygon hull4) (trianglesOfHull4 hull4)
+triangulateLeaf _ = Nothing
+
+trianglesOfHull4 :: Hull4 -> [Triangle]
+trianglesOfHull4 (Quadrilateral p1 p2 p3 p4) =
+  let (t1, t2) = trianglesOnQuadrilateral p1 p2 p3 p4 in [t1, t2]
+trianglesOfHull4 (TriangleWithInner p1 p2 p3 inner) =
+  [mkTriangle p1 p2 inner, mkTriangle p2 p3 inner, mkTriangle p3 p1 inner]
diff --git a/src/Triangulation/Merge.hs b/src/Triangulation/Merge.hs
new file mode 100644
--- /dev/null
+++ b/src/Triangulation/Merge.hs
@@ -0,0 +1,81 @@
+-- | Merging two triangulations of point sets separated by an axis-parallel
+-- line: join the hulls with their common tangents and fill the gap between
+-- them with triangles.
+module Triangulation.Merge (
+  mergeTriangulations,
+) where
+
+import Data.List.NonEmpty qualified as NE
+import Triangulation.Flip (legalize)
+import Triangulation.Geometry.Edge (Edge)
+import Triangulation.Geometry.Point (
+  Point,
+  bottomRight,
+  cosSquaredAngle,
+  leftTop,
+  rightBottom,
+  topLeft,
+ )
+import Triangulation.Geometry.Polygon (Polygon (..), tangents)
+import Triangulation.Geometry.Ring (arc)
+import Triangulation.Geometry.Triangle (Triangle (..), isValidCandidate, mkTriangle, triangleEdges)
+import Triangulation.Store (Store)
+import Triangulation.Store qualified as Store
+import Triangulation.Types (Axis (..), Triangulation (..))
+
+-- | Merge two triangulations whose point sets were split along the given axis:
+-- the first one holds the points with smaller x (for 'X') or larger y (for 'Y').
+-- 'Nothing' only if the hulls are inconsistent (a bridge endpoint is not a hull vertex).
+mergeTriangulations :: Triangulation -> Triangulation -> Axis -> Maybe Triangulation
+mergeTriangulations (Triangulation hull1@(Polygon points1) store1) (Triangulation hull2@(Polygon points2) store2) axis = do
+  let bridge = case axis of
+        X -> (rightBottom points1, leftTop points2)
+        Y -> (bottomRight points1, topLeft points2)
+  ((bl, br), (tl, tr), mergedHull) <- tangents hull1 hull2 bridge
+  remains1 <- arc bl tl points1
+  remains2 <- arc tr br points2
+  let unitedStore = Store.union store1 store2
+  pure $
+    Triangulation
+      mergedHull
+      (fillGap unitedStore (NE.toList remains1) (reverse (NE.toList remains2)) [])
+
+-- | Fill the gap between two hull chains (left and right, sharing no points)
+-- with triangles, advancing along whichever chain gives the better triangle.
+fillGap :: Store -> [Point] -> [Point] -> [Edge] -> Store
+fillGap store lefts rights restrictedEdges = case (lefts, rights) of
+  ([], _) -> store
+  (_, []) -> store
+  ([_], [_]) -> store
+  (left : lefts'@(nextLeft : _), [right]) ->
+    addCandidate store (left, right, nextLeft) lefts' [right] restrictedEdges
+  ([left], right : rights'@(nextRight : _)) ->
+    addCandidate store (left, right, nextRight) [left] rights' restrictedEdges
+  (left : lefts'@(nextLeft : _), right : rights'@(nextRight : _)) ->
+    let leftIsValid = isValidCandidate left right nextLeft rights'
+        rightIsValid = isValidCandidate left right nextRight lefts'
+        (leftCandidate, rightCandidate) = (mkTriangle left right nextLeft, mkTriangle left right nextRight)
+        advanceLeft = (leftCandidate, lefts', rights)
+        advanceRight = (rightCandidate, lefts, rights')
+        (triangle, lefts'', rights'')
+          | leftIsValid && not rightIsValid = advanceLeft
+          | rightIsValid && not leftIsValid = advanceRight
+          | minCosSquared leftCandidate >= minCosSquared rightCandidate = advanceLeft
+          | otherwise = advanceRight
+        store' = legalize (Store.insert triangle store) (triangleEdges triangle) restrictedEdges
+     in fillGap store' lefts'' rights'' restrictedEdges
+
+-- | Add the candidate triangle if it is valid, then continue filling the gap.
+addCandidate :: Store -> (Point, Point, Point) -> [Point] -> [Point] -> [Edge] -> Store
+addCandidate store (p1, p2, p3) lefts rights restrictedEdges
+  | isValidCandidate p1 p2 p3 [] =
+      let triangle = mkTriangle p1 p2 p3
+          store' = legalize (Store.insert triangle store) (triangleEdges triangle) restrictedEdges
+       in fillGap store' lefts rights restrictedEdges
+  | otherwise = fillGap store lefts rights restrictedEdges
+
+-- | The smallest squared cosine over the three angles; larger means the
+-- triangle is closer to right-angled, which the merge prefers.
+minCosSquared :: Triangle -> Double
+minCosSquared (Triangle p1 p2 p3) =
+  minimum [cosSquaredAngle p1 p2 p3, cosSquaredAngle p2 p1 p3, cosSquaredAngle p3 p1 p2]
diff --git a/src/Triangulation/Mesh.hs b/src/Triangulation/Mesh.hs
new file mode 100644
--- /dev/null
+++ b/src/Triangulation/Mesh.hs
@@ -0,0 +1,74 @@
+{-# LANGUAGE DeriveAnyClass #-}
+{-# LANGUAGE DerivingStrategies #-}
+
+-- | The indexed form of a triangulation: an array of nodes and the triangles
+-- as triples of node indices. This is the shape finite element assembly,
+-- mesh file formats and graphics APIs expect, and it is where a triangulation
+-- built from 'Point's by coordinate turns into
+-- one addressed by integers.
+--
+-- The output is canonical: nodes are numbered in ascending 'Ord' order of the
+-- points, triangles are listed in ascending order of their index triples, and
+-- every triple is wound counter-clockwise (with the y axis up, so its signed
+-- area is positive). Two runs on the same triangles give the same mesh.
+module Triangulation.Mesh (
+  Mesh (..),
+  toMesh,
+  fromMesh,
+  boundaryEdges,
+) where
+
+import Control.DeepSeq (NFData)
+import Data.HashMap.Strict qualified as HM
+import Data.List (sort)
+import Data.List.NonEmpty qualified as NE
+import Data.Maybe (mapMaybe)
+import Data.Vector (Vector)
+import Data.Vector qualified as V
+import GHC.Generics (Generic)
+import Triangulation.Geometry.Point (Orientation (..), Point, orientation)
+import Triangulation.Geometry.Triangle (Triangle (..), mkTriangle)
+
+-- | A triangulation with its vertices numbered.
+data Mesh = Mesh
+  { nodes :: !(Vector Point)
+  -- ^ the distinct vertices, in ascending order
+  , elements :: !(Vector (Int, Int, Int))
+  -- ^ the triangles as indices into 'nodes', each wound counter-clockwise
+  }
+  deriving stock (Eq, Show, Generic)
+  deriving anyclass (NFData)
+
+-- | Number the vertices of the triangles and express the triangles through
+-- the numbers. Duplicate triangles are kept once. A zero-area triangle,
+-- which has no winding, keeps its vertices in ascending order.
+toMesh :: [Triangle] -> Mesh
+toMesh triangles = Mesh (V.fromList points) (V.fromList indexed)
+  where
+    points = distinct $ concatMap (\(Triangle a b c) -> [a, b, c]) triangles
+    index = HM.fromList (zip points [0 ..])
+    -- every vertex is indexed, so the lookups never fail
+    indexed = distinct $ mapMaybe counterClockwise triangles
+    counterClockwise (Triangle a b c) = case orientation a b c of
+      Clockwise -> (,,) <$> at a <*> at c <*> at b
+      _ -> (,,) <$> at a <*> at b <*> at c
+    at p = HM.lookup p index
+
+-- | Sorted, without repetitions.
+distinct :: Ord a => [a] -> [a]
+distinct xs = [x | x NE.:| _ <- NE.group (sort xs)]
+
+-- | The triangles back as geometry.
+fromMesh :: Mesh -> [Triangle]
+fromMesh (Mesh points triangles) =
+  [mkTriangle (points V.! i) (points V.! j) (points V.! k) | (i, j, k) <- V.toList triangles]
+
+-- | The edges that belong to exactly one triangle: the outer boundary and the
+-- boundaries of the holes. Each edge is a pair of node indices with the
+-- smaller one first; the list is in ascending order.
+boundaryEdges :: Mesh -> [(Int, Int)]
+boundaryEdges (Mesh _ triangles) =
+  [e | e NE.:| [] <- NE.group . sort $ concatMap edgesOf (V.toList triangles)]
+  where
+    edgesOf (i, j, k) = [ordered i j, ordered j k, ordered k i]
+    ordered a b = (min a b, max a b)
diff --git a/src/Triangulation/Parallel.hs b/src/Triangulation/Parallel.hs
new file mode 100644
--- /dev/null
+++ b/src/Triangulation/Parallel.hs
@@ -0,0 +1,138 @@
+-- | Divide-and-conquer driver: split the points along the longer side of their
+-- bounding box, triangulate the halves (in parallel while they are large),
+-- merge.
+module Triangulation.Parallel (
+  triangulate,
+) where
+
+import Control.Parallel.Strategies (rdeepseq, rparWith, runEval)
+import Data.List (sortBy)
+import Data.List qualified as List
+import Triangulation.Geometry.Point (Point (..))
+import Triangulation.Leaf (triangulateLeaf)
+import Triangulation.Merge (mergeTriangulations)
+import Triangulation.Repair (repairDegeneracies)
+import Triangulation.Store qualified as Store
+import Triangulation.Types (Axis (..), Triangulation (..))
+
+-- | Delaunay triangulation of a point set; 'Nothing' when the points do not
+-- span a triangle: fewer than three of them, or all of them on one line.
+triangulate :: [Point] -> Maybe Triangulation
+triangulate points
+  | length points < 3 = Nothing
+  | otherwise = do
+      triangulation <- repairDegeneracies [] <$> go points
+      if Store.null (triangleStore triangulation) then Nothing else Just triangulation
+
+-- | Point sets at least this large have their halves evaluated in parallel.
+parallelThreshold :: Int
+parallelThreshold = 256
+
+go :: [Point] -> Maybe Triangulation
+go points
+  | Just leaf <- triangulateLeaf points = Just leaf
+  | otherwise = do
+      let Extent n minX maxX minY maxY = extentOf points
+          axis = if maxX - minX > maxY - minY then X else Y
+          (lefts, rights) = split n axis points
+          (left, right)
+            | n >= parallelThreshold = runEval $ do
+                -- one spark for the left half; the right one is evaluated here,
+                -- which is what sparking both used to end up doing anyway, only
+                -- after paying for a spark that then fizzled
+                l <- rparWith rdeepseq (go lefts)
+                r <- rdeepseq (go rights)
+                pure (l, r)
+            | otherwise = (go lefts, go rights)
+      l <- left
+      r <- right
+      mergeTriangulations l r axis
+
+-- | Split into two parts along the axis. Halves, except that 9–11 points go
+-- 3 + rest (halving would leave a part of 5, which cannot be split into two
+-- leaves), and 5 points go 2 + 3 (see 'triangulateLeaf').
+split :: Int -> Axis -> [Point] -> ([Point], [Point])
+split n axis points
+  | n > selectionThreshold = selectSmallest (pointsOrder axis) leftSize n points
+  | otherwise = splitAt leftSize (sortBy (pointsOrder axis) points)
+  where
+    leftSize = if n `elem` [9, 10, 11] then 3 else n `div` 2
+
+-- | Sets larger than this are partitioned by selecting the median; smaller
+-- ones are sorted outright. Sorting a short list is cheap, and it keeps the
+-- points reaching a leaf in the order they have always arrived in, which the
+-- two-point leaf of a five-point set depends on.
+selectionThreshold :: Int
+selectionThreshold = 32
+
+-- | The @k@ smallest of @n@ points by the given order, and the rest.
+--
+-- Neither part comes back ordered, and neither needs to be: the divide step
+-- only requires that the two halves be separated by a line, and each half is
+-- partitioned again along its own axis. Selecting the median takes a pass per
+-- level of the selection instead of the @n log n@ of a full sort, and the
+-- sorting used to be about two fifths of the running time.
+selectSmallest :: (Point -> Point -> Ordering) -> Int -> Int -> [Point] -> ([Point], [Point])
+selectSmallest order = select
+  where
+    select k n points
+      | k <= 0 = ([], points)
+      | k >= n = (points, [])
+      | otherwise = case drop (n `div` 2) points of
+          [] -> (points, [])
+          pivot : _ ->
+            let (smaller, equal, larger) = partitionAround pivot points
+                belowCount = length smaller
+                equalCount = length equal
+             in case compare k belowCount of
+                  LT -> keepLeft (select k belowCount smaller) equal larger
+                  EQ -> (smaller, equal ++ larger)
+                  GT
+                    | k <= belowCount + equalCount ->
+                        let (taken, left) = splitAt (k - belowCount) equal
+                         in (smaller ++ taken, left ++ larger)
+                    | otherwise ->
+                        keepRight smaller equal (select (k - belowCount - equalCount) (n - belowCount - equalCount) larger)
+    keepLeft (chosen, rest) equal larger = (chosen, rest ++ equal ++ larger)
+    keepRight smaller equal (chosen, rest) = (smaller ++ equal ++ chosen, rest)
+    partitionAround pivot = foldr step ([], [], [])
+      where
+        step p (below, same, above) = case order p pivot of
+          LT -> (p : below, same, above)
+          EQ -> (below, p : same, above)
+          GT -> (below, same, p : above)
+
+-- | Along 'X' points go left to right; along 'Y' they go top to bottom.
+--
+-- Ties are broken the way the symbolic perturbation behind
+-- 'Triangulation.Geometry.Point.turn' would break
+-- them (a lower rank is displaced further, so among points with the same x
+-- the one with the smaller y is further right, and among points with the same
+-- y the one with the smaller x is further up). The two halves are then
+-- separated by a straight line in the perturbed plane, which the merge relies
+-- on; splitting ties arbitrarily would put points of one half between points
+-- of the other.
+-- Written out rather than through @comparing@ on a tuple: the comparator runs
+-- once per comparison of every sort at every level of the recursion, and the
+-- tuple and the 'Data.Ord.Down' wrapper were allocated every time.
+pointsOrder :: Axis -> Point -> Point -> Ordering
+pointsOrder X (Point ax ay) (Point bx by) = case compare ax bx of
+  EQ -> compare by ay
+  unequal -> unequal
+pointsOrder Y (Point ax ay) (Point bx by) = case compare by ay of
+  EQ -> compare ax bx
+  unequal -> unequal
+{-# INLINE pointsOrder #-}
+
+-- | How many points there are and the corners of their bounding box.
+data Extent = Extent !Int !Double !Double !Double !Double
+
+-- | One pass for the count and all four extremes; the axis and the split size
+-- both come from it. Asking for the length and then for each extreme
+-- separately walked the list five times per node of the recursion.
+extentOf :: [Point] -> Extent
+extentOf = List.foldl' step (Extent 0 inf (-inf) inf (-inf))
+  where
+    step (Extent n minX maxX minY maxY) (Point x y) =
+      Extent (n + 1) (min minX x) (max maxX x) (min minY y) (max maxY y)
+    inf = 1 / 0
diff --git a/src/Triangulation/Refine.hs b/src/Triangulation/Refine.hs
new file mode 100644
--- /dev/null
+++ b/src/Triangulation/Refine.hs
@@ -0,0 +1,310 @@
+{-# LANGUAGE DeriveAnyClass #-}
+{-# LANGUAGE DerivingStrategies #-}
+
+-- | Delaunay refinement: Ruppert's algorithm.
+--
+-- Vertices are inserted into a constrained Delaunay triangulation until every
+-- triangle meets the quality bounds. The /segments/ — the edges that belong to
+-- a single triangle, that is the boundary of the region and of its holes —
+-- are preserved as chains: a segment is only ever split at a point on it.
+--
+-- The two rules, applied until neither fires:
+--
+-- 1. A segment /encroached/ upon by a vertex (one strictly inside its
+--    diametral circle) is split.
+-- 2. A poor triangle (an angle below the bound, or an area above it) has its
+--    circumcenter inserted — unless the circumcenter would encroach upon
+--    segments, in which case those segments are split instead.
+--
+-- Every insertion is a Bowyer–Watson step: the triangles whose circumcircles
+-- contain the new vertex are removed, without crossing a segment, and the
+-- cavity is re-triangulated as a fan around the vertex, so the triangulation
+-- stays constrained Delaunay.
+--
+-- Segments meeting at an input vertex at an angle below 60° are handled
+-- the way Shewchuk's Triangle does: a segment with one such endpoint is
+-- split at a power-of-two distance from it (concentric shells), and a
+-- triangle whose smallest angle is formed by two segments is not asked to
+-- meet the angle bound, since no refinement can fix an input angle.
+module Triangulation.Refine (
+  Quality (..),
+  defaultQuality,
+  refine,
+  refineWithBudget,
+) where
+
+import Control.DeepSeq (NFData)
+import Data.HashMap.Strict qualified as HM
+import Data.HashSet (HashSet)
+import Data.HashSet qualified as HS
+import Data.List (sortOn)
+import GHC.Generics (Generic)
+import Triangulation.Geometry.Edge (Edge (..), mkEdge)
+import Triangulation.Geometry.Point (Orientation (..), Point (..), orientation)
+import Triangulation.Geometry.Triangle (
+  Triangle (..),
+  circumcenter,
+  isOutsideCircumcircle,
+  mkTriangle,
+  smallestAngle,
+  triangleArea,
+  triangleEdges,
+ )
+import Triangulation.Store (Store)
+import Triangulation.Store qualified as Store
+
+-- | What the refined triangulation must satisfy.
+data Quality = Quality
+  { minAngle :: !Double
+  -- ^ lower bound on every angle, in degrees; @0@ imposes none. Ruppert's
+  -- algorithm is guaranteed to terminate up to about 20.7°, and usually does
+  -- up to about 30°.
+  , maxArea :: !(Maybe Double)
+  -- ^ upper bound on the area of a triangle
+  , maxInsertions :: !Int
+  -- ^ how many vertices may be inserted before giving up
+  }
+  deriving stock (Eq, Show, Generic)
+  deriving anyclass (NFData)
+
+-- | A 20° angle bound, no area bound, at most 100 000 insertions.
+defaultQuality :: Quality
+defaultQuality = Quality {minAngle = 20, maxArea = Nothing, maxInsertions = 100000}
+
+-- | Refine the triangles of a region — typically the output of
+-- 'Triangulation.Constrained.constrainedTriangulate' — until every triangle
+-- meets the bounds. 'Nothing' when the insertion budget runs out, which is
+-- what happens when the bounds cannot be met.
+refine :: Quality -> [Triangle] -> Maybe [Triangle]
+refine quality triangles = case refineWithBudget quality triangles of
+  (True, refined) -> Just refined
+  (False, _) -> Nothing
+
+-- | Like 'refine', but never gives up empty-handed. The 'Bool' says whether
+-- the bounds were met; when the budget ran out it is 'False' and the
+-- triangles are the mesh reached so far — a valid constrained Delaunay
+-- triangulation of the region that does not meet the bounds.
+refineWithBudget :: Quality -> [Triangle] -> (Bool, [Triangle])
+refineWithBudget quality triangles
+  | null triangles = (True, [])
+  | otherwise = let (done, r) = loop quality initial in (done, Store.triangles (store r))
+  where
+    store0 = foldr Store.insert Store.empty triangles
+    segments0 = HS.fromList [e | e <- Store.edges store0, [_] <- [Store.trianglesOn e store0]]
+    initial =
+      Refinement
+        { store = store0
+        , count = length triangles
+        , segments = segments0
+        , acute = acuteVertices segments0
+        , budget = maxInsertions quality
+        , pendingSegments = HS.toList segments0
+        , pendingTriangles = triangles
+        }
+
+-- | The state of the refinement.
+data Refinement = Refinement
+  { store :: !Store
+  , count :: !Int
+  -- ^ number of triangles in the store
+  , segments :: !(HashSet Edge)
+  , acute :: !(HashSet Point)
+  -- ^ input vertices where two segments meet at less than 60°
+  , budget :: !Int
+  , pendingSegments :: ![Edge]
+  -- ^ segments to test for encroachment
+  , pendingTriangles :: ![Triangle]
+  -- ^ triangles to test for quality
+  }
+
+-- | Split encroached segments first; when none is left, fix poor triangles.
+-- 'False' when the budget ran out, with the state at that point.
+loop :: Quality -> Refinement -> (Bool, Refinement)
+loop quality r = case pendingSegments r of
+  s : rest
+    | isSegment s r && encroachedByApex s r -> continue (splitSegment s r {pendingSegments = rest})
+    | otherwise -> loop quality r {pendingSegments = rest}
+  [] -> case pendingTriangles r of
+    t : rest
+      | Store.hasTriangle t (store r) && isPoor quality r t ->
+          continue (fixTriangle t r {pendingTriangles = rest})
+      | otherwise -> loop quality r {pendingTriangles = rest}
+    [] -> (True, r)
+  where
+    continue = maybe (False, r) (loop quality)
+
+isSegment :: Edge -> Refinement -> Bool
+isSegment e r = HS.member e (segments r)
+
+-- | Whether the point lies strictly inside the diametral circle of the edge.
+encroaches :: Point -> Edge -> Bool
+encroaches (Point x y) (Edge (Point x1 y1) (Point x2 y2)) = (x - x1) * (x - x2) + (y - y1) * (y - y2) < 0
+
+-- | A segment encroached upon by any vertex is encroached upon by the apex of
+-- a triangle on it (Shewchuk), so only the apexes need testing.
+encroachedByApex :: Edge -> Refinement -> Bool
+encroachedByApex e r = any (`encroaches` e) [apex e t | t <- Store.trianglesOn e (store r)]
+
+-- | The vertex of the triangle not on the edge.
+apex :: Edge -> Triangle -> Point
+apex (Edge a b) (Triangle x y z) = case filter (`notElem` [a, b]) [x, y, z] of
+  p : _ -> p
+  [] -> x -- unreachable: an edge of a triangle has exactly two of its vertices
+
+-- | Input vertices where two segments meet at an angle below 60°.
+acuteVertices :: HashSet Edge -> HashSet Point
+acuteVertices segs = HS.fromList [v | (v, ns) <- HM.toList neighbours, hasSmallAngle v ns]
+  where
+    neighbours = HM.fromListWith (++) (concat [[(a, [b]), (b, [a])] | Edge a b <- HS.toList segs])
+    hasSmallAngle v ns = or [cosine v n1 n2 > 0.5 | (n1, i) <- zip ns [0 :: Int ..], (n2, j) <- zip ns [0 ..], i < j]
+    cosine (Point x0 y0) (Point x1 y1) (Point x2 y2) =
+      let (dx1, dy1, dx2, dy2) = (x1 - x0, y1 - y0, x2 - x0, y2 - y0)
+       in (dx1 * dx2 + dy1 * dy2) / sqrt ((dx1 * dx1 + dy1 * dy1) * (dx2 * dx2 + dy2 * dy2))
+
+-- | Whether the triangle violates a bound. The angle bound is waived when the
+-- smallest angle lies between two segments: it is an input angle.
+isPoor :: Quality -> Refinement -> Triangle -> Bool
+isPoor quality r t = tooSmallAngle || tooLarge
+  where
+    tooSmallAngle = minAngle quality > 0 && smallestAngle t < minAngle quality && not inputAngle
+    tooLarge = any (triangleArea t >) (maxArea quality)
+    inputAngle = all (`isSegment` r) (longerEdges t)
+    -- the smallest angle is opposite the shortest edge, between the two others
+    longerEdges = drop 1 . sortOn edgeLength . triangleEdges
+    edgeLength (Edge (Point x1 y1) (Point x2 y2)) = (x2 - x1) * (x2 - x1) + (y2 - y1) * (y2 - y1)
+
+-- | Split a segment at its midpoint, or at a power-of-two distance from an
+-- acute endpoint, and insert the new vertex.
+splitSegment :: Edge -> Refinement -> Maybe Refinement
+splitSegment e@(Edge a b) r
+  | p == a || p == b = Just r -- too short to split in floating point
+  | otherwise = case Store.trianglesOn e (store r) of
+      t : _ -> insertOnSegment p e t r
+      [] -> Just r
+  where
+    p
+      | HS.member b (acute r) = shell b a
+      | HS.member a (acute r) = shell a b
+      | otherwise = Point ((px a + px b) / 2) ((py a + py b) / 2)
+    -- the split point at a power-of-two distance from @from@, between a third
+    -- and two thirds of the way, so that splits on different segments meeting
+    -- at the acute vertex land on common concentric circles
+    shell from to =
+      let len = sqrt ((px to - px from) ^ (2 :: Int) + (py to - py from) ^ (2 :: Int))
+          d = 2 ^^ (floor (logBase 2 (2 * len / 3)) :: Int)
+          f = d / len
+       in Point (px from + f * (px to - px from)) (py from + f * (py to - py from))
+
+-- | Insert the circumcenter of a poor triangle, or split the segments it
+-- would encroach upon.
+fixTriangle :: Triangle -> Refinement -> Maybe Refinement
+fixTriangle t r = case circumcenter t of
+  Nothing -> Just r
+  Just c -> case locate c t r of
+    -- The circumcenter lies outside the region: as in Shewchuk's Triangle,
+    -- the segment the walk left through is split instead.
+    Exited e -> splitEncroached c [e]
+    Lost -> Just r
+    Found tc
+      | isVertex c tc -> Just r
+      | otherwise ->
+          let cav = cavity c tc r
+              boundary = cavityBoundary cav
+           in case [e | e <- boundary, isSegment e r, encroaches c e] of
+                [] -> insertInCavity c [] cav boundary r
+                encroached -> splitEncroached c encroached
+  where
+    -- Split the segments and look at the triangle again afterwards; if
+    -- nothing could be split (the segments are too short to halve), leave it.
+    splitEncroached _ encroached = do
+      r' <- foldr (\e acc -> acc >>= splitSegment e) (Just r) encroached
+      pure $ if budget r' < budget r then requeue r' else r'
+    requeue r' = r' {pendingTriangles = pendingTriangles r' ++ [t]}
+    isVertex p (Triangle x y z) = p `elem` [x, y, z]
+
+-- | Where a walk from a triangle towards a point ends.
+data Location
+  = -- | the triangle containing the point (possibly on its boundary)
+    Found Triangle
+  | -- | the walk left the region through this segment
+    Exited Edge
+  | -- | the walk did not terminate in a reasonable number of steps
+    Lost
+
+-- | Walk from the triangle towards the point, crossing at each step an edge
+-- that separates the point from the current triangle.
+locate :: Point -> Triangle -> Refinement -> Location
+locate p start r = go (0 :: Int) Nothing start
+  where
+    go steps cameFrom t
+      | steps > count r = Lost
+      | inside t = Found t
+      | otherwise = case [e | (e@(Edge u v), w) <- edgesWithApex t, Just e /= cameFrom, separates u v w] of
+          [] -> Lost
+          e : _
+            | isSegment e r -> Exited e
+            | otherwise -> case filter (/= t) (Store.trianglesOn e (store r)) of
+                n : _ -> go (steps + 1) (Just e) n
+                [] -> Exited e
+    inside (Triangle a b c) = all sameSide [(a, b, c), (b, c, a), (c, a, b)]
+    sameSide (u, v, w) = let o = orientation u v p in o == Collinear || o == orientation u v w
+    separates u v w = let o = orientation u v p in o /= Collinear && o /= orientation u v w
+    edgesWithApex (Triangle a b c) = [(mkEdge a b, c), (mkEdge b c, a), (mkEdge c a, b)]
+
+-- | The triangles whose circumcircles contain the point, reached from the
+-- containing triangle without crossing a segment.
+cavity :: Point -> Triangle -> Refinement -> HashSet Triangle
+cavity p start r = go (HS.singleton start) [start]
+  where
+    go seen [] = seen
+    go seen (t : queue) =
+      let fresh =
+            [ n
+            | e <- triangleEdges t
+            , not (isSegment e r)
+            , n <- Store.trianglesOn e (store r)
+            , n /= t
+            , not (HS.member n seen)
+            , not (isOutsideCircumcircle p n)
+            ]
+       in go (foldr HS.insert seen fresh) (fresh ++ queue)
+
+-- | The edges of the cavity that belong to exactly one of its triangles.
+cavityBoundary :: HashSet Triangle -> [Edge]
+cavityBoundary cav =
+  [ e
+  | (e, 1 :: Int) <- HM.toList (HM.fromListWith (+) [(e, 1) | t <- HS.toList cav, e <- triangleEdges t])
+  ]
+
+-- | Insert a vertex on the given segment, which has the given triangle on it.
+-- The segment is split into two at the vertex; it is not tested for
+-- collinearity, because the split point is computed in floating point and
+-- need not lie exactly on the line.
+insertOnSegment :: Point -> Edge -> Triangle -> Refinement -> Maybe Refinement
+insertOnSegment p e t r =
+  let cav = cavity p t r
+   in insertInCavity p [e] cav (cavityBoundary cav) r
+
+-- | Replace the cavity by a fan of triangles around the new vertex. The given
+-- edges, and any boundary edge exactly collinear with the vertex, are the
+-- segments the vertex lies on: each is split into two segments instead of
+-- becoming a flat triangle.
+insertInCavity :: Point -> [Edge] -> HashSet Triangle -> [Edge] -> Refinement -> Maybe Refinement
+insertInCavity p splitting cav boundary r
+  | budget r <= 0 = Nothing
+  | otherwise =
+      Just
+        r
+          { store = foldr Store.insert (foldr Store.delete (store r) (HS.toList cav)) newTriangles
+          , count = count r - HS.size cav + length newTriangles
+          , segments = foldr HS.insert (foldr HS.delete (segments r) split) newSegments
+          , budget = budget r - 1
+          , pendingSegments = newSegments ++ filter (`isSegment` r) fan ++ pendingSegments r
+          , pendingTriangles = newTriangles ++ pendingTriangles r
+          }
+  where
+    isSplit e@(Edge a b) = e `elem` splitting || orientation a b p == Collinear
+    split = filter isSplit boundary
+    fan = filter (not . isSplit) boundary
+    newTriangles = [mkTriangle a b p | Edge a b <- fan]
+    newSegments = concat [[mkEdge a p, mkEdge p b] | e@(Edge a b) <- split, isSegment e r]
diff --git a/src/Triangulation/Repair.hs b/src/Triangulation/Repair.hs
new file mode 100644
--- /dev/null
+++ b/src/Triangulation/Repair.hs
@@ -0,0 +1,94 @@
+-- | Repair of degeneracies left behind by collinear input.
+--
+-- Three collinear points that land in the same leaf of the divide-and-conquer
+-- tree become a zero-area triangle, and the merge builds on it: its long edge
+-- ends up spanning the middle vertex, which later blocks constrained edges.
+-- Once the whole triangulation is assembled these spots are easy to find and
+-- fix locally, so that is done here rather than in every step that could
+-- produce one.
+module Triangulation.Repair (
+  repairDegeneracies,
+) where
+
+import Data.HashSet qualified as HS
+import Data.List (findIndex)
+import Data.List qualified as List
+import Data.List.NonEmpty qualified as NE
+import Triangulation.Flip (legalize)
+import Triangulation.Geometry.Edge (Edge, mkEdge)
+import Triangulation.Geometry.Point (Orientation (..), Point (..), orientation)
+import Triangulation.Geometry.Polygon (Polygon (..))
+import Triangulation.Geometry.Ring (cyclicPairs)
+import Triangulation.Geometry.Triangle (Triangle (..), mkTriangle, triangleEdges)
+import Triangulation.Store (Store)
+import Triangulation.Store qualified as Store
+import Triangulation.Types (Triangulation (..))
+
+-- | Remove every zero-area triangle, re-triangulating its neighbourhood, and
+-- put the vertices this leaves on the hull boundary into the hull ring. The
+-- given edges are constraints that the Delaunay flips must not remove.
+-- Input in general position has no flat triangles and passes through
+-- untouched.
+repairDegeneracies :: [Edge] -> Triangulation -> Triangulation
+repairDegeneracies restrictedEdges triangulation@(Triangulation polygon store)
+  | not (Store.anyTriangle isFlat store) = triangulation
+  | otherwise =
+      let store' = removeFlats restrictedEdges store
+          points = HS.toList . HS.fromList $ concatMap (\(Triangle a b c) -> [a, b, c]) (Store.triangles store')
+       in Triangulation (foldr insertOnBoundary polygon points) store'
+
+-- | Repair flat triangles until none is left.
+--
+-- One scan of the store collects every flat triangle there is and all of them
+-- are repaired before the store is scanned again, because repairing one can
+-- leave a new one behind. Taking the first flat triangle of a fresh scan each
+-- time, as this used to, costs a scan per repair; on input like a lattice,
+-- where a large share of the triangles is flat, that made the whole
+-- triangulation quadratic.
+removeFlats :: [Edge] -> Store -> Store
+removeFlats restrictedEdges = rounds
+  where
+    rounds store = case filter isFlat (Store.triangles store) of
+      [] -> store
+      flats -> rounds (List.foldl' repair store flats)
+    repair store flat
+      | Store.hasTriangle flat store = repairOne restrictedEdges flat store
+      | otherwise = store -- an earlier repair in this round removed it
+
+isFlat :: Triangle -> Bool
+isFlat (Triangle a b c) = orientation a b c == Collinear
+
+-- | The vertices of a triangle are kept sorted, so in a flat triangle @a b c@
+-- the point @b@ lies between @a@ and @c@ and the edge @a c@ spans it. The
+-- triangle on the other side of that edge, @a c d@, is split at @b@ into
+-- @a b d@ and @b c d@; without such a triangle the flat one is simply dropped.
+repairOne :: [Edge] -> Triangle -> Store -> Store
+repairOne restrictedEdges flat@(Triangle a b c) store =
+  let store1 = Store.delete flat store
+   in case Store.trianglesOn (mkEdge a c) store1 of
+        [neighbour@(Triangle x y z)]
+          | [d] <- filter (`notElem` [a, c]) [x, y, z] ->
+              let replacements = filter (not . isFlat) [mkTriangle a b d, mkTriangle b c d]
+                  store2 = foldr Store.insert (Store.delete neighbour store1) replacements
+               in legalize store2 (concatMap triangleEdges replacements) restrictedEdges
+        _ -> store1
+
+-- | Insert a vertex lying strictly inside a hull edge between that edge's
+-- endpoints; a vertex already on the ring, or not on any edge, is left alone.
+insertOnBoundary :: Point -> Polygon -> Polygon
+insertOnBoundary p polygon@(Polygon ring)
+  | p `elem` ring = polygon
+  | otherwise = case findIndex (\(u, v) -> liesBetween u v p) (cyclicPairs ring) of
+      Just i ->
+        let (front, back) = splitAt (i + 1) (NE.toList ring) in Polygon (NE.fromList (front ++ p : back))
+      Nothing -> polygon
+
+liesBetween :: Point -> Point -> Point -> Bool
+liesBetween u v p =
+  p /= u
+    && p /= v
+    && orientation u v p == Collinear
+    && min (px u) (px v) <= px p
+    && px p <= max (px u) (px v)
+    && min (py u) (py v) <= py p
+    && py p <= max (py u) (py v)
diff --git a/src/Triangulation/Store.hs b/src/Triangulation/Store.hs
new file mode 100644
--- /dev/null
+++ b/src/Triangulation/Store.hs
@@ -0,0 +1,124 @@
+{-# LANGUAGE DeriveAnyClass #-}
+{-# LANGUAGE DerivingStrategies #-}
+
+-- | The triangles of a triangulation, indexed by edge: every edge maps to the
+-- (at most two) triangles it belongs to.
+--
+-- Meant to be imported qualified:
+--
+-- > import Triangulation.Store (Store)
+-- > import Triangulation.Store qualified as Store
+module Triangulation.Store (
+  Store (..),
+  empty,
+  null,
+  anyTriangle,
+  union,
+  member,
+  hasTriangle,
+  trianglesOn,
+  edges,
+  triangles,
+  insert,
+  delete,
+) where
+
+import Control.DeepSeq (NFData)
+import Data.HashMap.Strict qualified as HM
+import Data.HashSet qualified as HS
+import Data.List qualified as List
+import GHC.Generics (Generic)
+import Triangulation.Geometry.Edge (Edge)
+import Triangulation.Geometry.Triangle (Triangle, triangleEdges)
+import Prelude hiding (null)
+
+-- | Triangles indexed by their edges.
+-- | The triangles on one edge are kept as a short list rather than a
+-- 'Data.HashSet.HashSet': an edge belongs to one or two triangles, and a hash
+-- set of that size costs an allocation per insertion for nothing.
+newtype Store = Store (HM.HashMap Edge Adjacent)
+  deriving stock (Show)
+  deriving newtype (NFData)
+
+-- | The triangles on one edge. In a planar triangulation an edge belongs to
+-- one or two of them, and spelling that out keeps the common case free of
+-- cons cells and of a set per edge. 'Several' exists only so that the type is
+-- total: nothing in the library builds it.
+data Adjacent
+  = One !Triangle
+  | Two !Triangle !Triangle
+  | Several ![Triangle]
+  deriving stock (Show, Generic)
+  deriving anyclass (NFData)
+
+adjacentList :: Adjacent -> [Triangle]
+adjacentList (One t) = [t]
+adjacentList (Two t u) = [t, u]
+adjacentList (Several ts) = ts
+{-# INLINE adjacentList #-}
+
+-- | No triangles.
+empty :: Store
+empty = Store HM.empty
+
+-- | Whether the store holds no triangles.
+null :: Store -> Bool
+null (Store store) = HM.null store
+
+-- | Whether some triangle satisfies the predicate. Cheaper than filtering
+-- 'triangles': no set of all triangles is built.
+anyTriangle :: (Triangle -> Bool) -> Store -> Bool
+anyTriangle p (Store store) = any (any p . adjacentList) (HM.elems store)
+
+-- | The triangles of both stores; meant for stores on disjoint point sets.
+union :: Store -> Store -> Store
+union (Store store1) (Store store2) = Store $ HM.union store1 store2
+
+-- | Whether the edge belongs to some triangle of the store.
+member :: Edge -> Store -> Bool
+member edge (Store store) = HM.member edge store
+
+-- | Whether the triangle is in the store.
+hasTriangle :: Triangle -> Store -> Bool
+hasTriangle triangle (Store store) = case triangleEdges triangle of
+  edge : _ -> any (elem triangle . adjacentList) (HM.lookup edge store)
+  [] -> False
+
+-- | The triangles the edge belongs to.
+trianglesOn :: Edge -> Store -> [Triangle]
+trianglesOn edge (Store store) = maybe [] adjacentList (HM.lookup edge store)
+
+-- | Every edge of every triangle.
+edges :: Store -> [Edge]
+edges (Store store) = HM.keys store
+
+-- | Every triangle, once.
+triangles :: Store -> [Triangle]
+triangles (Store store) = HS.toList . HS.fromList . concatMap adjacentList $ HM.elems store
+
+-- | Add a triangle under each of its edges.
+insert :: Triangle -> Store -> Store
+insert triangle (Store store) = Store $ List.foldl' insertOn store (triangleEdges triangle)
+  where
+    insertOn store' edge = HM.insertWith addUnique edge (One triangle) store'
+    addUnique _ present = case present of
+      One t | t /= triangle -> Two t triangle
+      Two t u | t /= triangle && u /= triangle -> Several [triangle, t, u]
+      Several ts | triangle `notElem` ts -> Several (triangle : ts)
+      _ -> present
+
+-- | Remove a triangle from each of its edges; edges left without triangles disappear.
+delete :: Triangle -> Store -> Store
+delete triangle (Store store) = Store $ List.foldl' deleteOn store (triangleEdges triangle)
+  where
+    deleteOn store' edge = HM.update remaining edge store'
+    remaining present = case present of
+      One t | t == triangle -> Nothing
+      Two t u | t == triangle -> Just (One u)
+      Two t u | u == triangle -> Just (One t)
+      Several ts -> case filter (/= triangle) ts of
+        [] -> Nothing
+        [t] -> Just (One t)
+        [t, u] -> Just (Two t u)
+        rest -> Just (Several rest)
+      _ -> Just present
diff --git a/src/Triangulation/Types.hs b/src/Triangulation/Types.hs
new file mode 100644
--- /dev/null
+++ b/src/Triangulation/Types.hs
@@ -0,0 +1,40 @@
+{-# LANGUAGE DeriveAnyClass #-}
+{-# LANGUAGE DerivingStrategies #-}
+
+-- | The triangulation itself: the convex hull plus the triangles.
+module Triangulation.Types (
+  Triangulation (..),
+  triangles,
+  fromTriangles,
+  Axis (..),
+) where
+
+import Control.DeepSeq (NFData)
+import GHC.Generics (Generic)
+import Triangulation.Geometry.Polygon (Polygon)
+import Triangulation.Geometry.Triangle (Triangle)
+import Triangulation.Store (Store)
+import Triangulation.Store qualified as Store
+
+-- | A triangulation of a point set.
+data Triangulation = Triangulation
+  { hull :: !Polygon
+  -- ^ the convex hull of the points
+  , triangleStore :: !Store
+  -- ^ the triangles, indexed by edge
+  }
+  deriving stock (Show, Generic)
+  deriving anyclass (NFData)
+
+-- | All triangles.
+triangles :: Triangulation -> [Triangle]
+triangles = Store.triangles . triangleStore
+
+-- | Build a triangulation from its hull and triangles.
+fromTriangles :: Polygon -> [Triangle] -> Triangulation
+fromTriangles polygon = Triangulation polygon . foldr Store.insert Store.empty
+
+-- | A coordinate axis; the direction in which a point set is split before
+-- the halves are triangulated and merged.
+data Axis = X | Y
+  deriving stock (Eq, Show)
diff --git a/test/Main.hs b/test/Main.hs
new file mode 100644
--- /dev/null
+++ b/test/Main.hs
@@ -0,0 +1,20 @@
+module Main (main) where
+
+import Test.Geometry qualified
+import Test.Mesh qualified
+import Test.Refine qualified
+import Test.Ring qualified
+import Test.Tasty (defaultMain, testGroup)
+import Test.Triangulation qualified
+
+main :: IO ()
+main =
+  defaultMain $
+    testGroup
+      "triangulation"
+      [ Test.Ring.tests
+      , Test.Geometry.tests
+      , Test.Triangulation.tests
+      , Test.Mesh.tests
+      , Test.Refine.tests
+      ]
diff --git a/test/Test/Geometry.hs b/test/Test/Geometry.hs
new file mode 100644
--- /dev/null
+++ b/test/Test/Geometry.hs
@@ -0,0 +1,218 @@
+{-# OPTIONS_GHC -Wno-orphans #-}
+
+module Test.Geometry (tests, Points (..)) where
+
+import Data.Hashable (hash)
+import Data.List (permutations)
+import Data.List.NonEmpty (NonEmpty (..))
+import Test.Tasty (TestTree, testGroup)
+import Test.Tasty.HUnit (assertBool, testCase, (@?=))
+import Test.Tasty.QuickCheck (
+  Arbitrary (..),
+  Gen,
+  choose,
+  chooseInt,
+  shrinkList,
+  testProperty,
+  vectorOf,
+  (===),
+ )
+import Triangulation.Geometry
+
+-- | Points inside the 2000×2000 canvas used by the examples.
+instance Arbitrary Point where
+  arbitrary = Point <$> choose (0, 2000) <*> choose (0, 2000)
+
+-- | A point set of at least three points; shrinking never goes below three.
+newtype Points = Points [Point]
+  deriving (Show)
+
+instance Arbitrary Points where
+  arbitrary = do
+    n <- chooseInt (3, 60)
+    Points <$> vectorOf n (arbitrary :: Gen Point)
+  shrink (Points ps) = [Points ps' | ps' <- shrinkList (const []) ps, length ps' >= 3]
+
+tests :: TestTree
+tests =
+  testGroup
+    "Geometry"
+    [ orientationTests
+    , exactnessTests
+    , edgeTests
+    , triangleTests
+    , polygonTests
+    ]
+
+-- | Triples that are collinear or nearly so, where naive floating-point
+-- evaluation gets the sign wrong: @c@ is on the line through @a@ and @b@ up
+-- to a perturbation of a few ulps.
+newtype NearlyCollinear = NearlyCollinear (Point, Point, Point)
+  deriving (Show)
+
+instance Arbitrary NearlyCollinear where
+  arbitrary = do
+    a <- arbitrary
+    b <- arbitrary
+    t <- choose (-1, 2 :: Double)
+    ulps <- chooseInt (-3, 3)
+    let Point cx cy = Point (px a + t * (px b - px a)) (py a + t * (py b - py a))
+    pure $ NearlyCollinear (a, b, Point (nudge ulps cx) cy)
+    where
+      nudge n x = x + fromIntegral n * (x * 2.220446049250313e-16)
+
+-- | Points on a circle (up to rounding) around a centre, plus one more nearly on it.
+newtype NearlyCocircular = NearlyCocircular (Point, Point, Point, Point)
+  deriving (Show)
+
+instance Arbitrary NearlyCocircular where
+  arbitrary = do
+    centre <- arbitrary
+    radius <- choose (1, 1000)
+    t1 <- choose (0, 2 * pi)
+    t2 <- choose (0, 2 * pi)
+    t3 <- choose (0, 2 * pi)
+    t4 <- choose (0, 2 * pi)
+    ulps <- chooseInt (-3, 3)
+    let onCircle t = Point (px centre + radius * cos t) (py centre + radius * sin t)
+        Point dx dy = onCircle t4
+    pure $
+      NearlyCocircular
+        ( onCircle t1
+        , onCircle t2
+        , onCircle t3
+        , Point (dx + fromIntegral ulps * dx * 2.220446049250313e-16) dy
+        )
+
+exactnessTests :: TestTree
+exactnessTests =
+  testGroup
+    "exact predicates"
+    [ testProperty "orientation agrees with rational arithmetic on random triples" $
+        \a b c -> orientation a b c === referenceOrientation a b c
+    , testProperty "orientation agrees with rational arithmetic on nearly collinear triples" $
+        \(NearlyCollinear (a, b, c)) -> orientation a b c === referenceOrientation a b c
+    , testProperty "in-circle agrees with rational arithmetic on nearly cocircular points" $
+        \(NearlyCocircular (a, b, c, d)) ->
+          inCircleSign (coordinates a) (coordinates b) (coordinates c) (coordinates d)
+            === referenceInCircle a b c d
+    ]
+  where
+    r = toRational
+    referenceOrientation (Point ax ay) (Point bx by) (Point cx cy) =
+      case compare ((r ax - r cx) * (r by - r cy) - (r ay - r cy) * (r bx - r cx)) 0 of
+        LT -> Clockwise
+        GT -> CounterClockwise
+        EQ -> Collinear
+    referenceInCircle (Point ax ay) (Point bx by) (Point cx cy) (Point dx dy) =
+      let (adx, ady, bdx, bdy, cdx, cdy) = (r ax - r dx, r ay - r dy, r bx - r dx, r by - r dy, r cx - r dx, r cy - r dy)
+       in compare
+            ( (adx * adx + ady * ady) * (bdx * cdy - cdx * bdy)
+                + (bdx * bdx + bdy * bdy) * (cdx * ady - adx * cdy)
+                + (cdx * cdx + cdy * cdy) * (adx * bdy - bdx * ady)
+            )
+            0
+
+orientationTests :: TestTree
+orientationTests =
+  testGroup
+    "orientation"
+    [ testCase "left turn is counter-clockwise (y axis up)" $
+        orientation (Point 0 0) (Point 1 0) (Point 1 1) @?= CounterClockwise
+    , testCase "right turn is clockwise" $
+        orientation (Point 0 0) (Point 1 0) (Point 1 (-1)) @?= Clockwise
+    , testCase "points on a line are collinear" $
+        orientation (Point 0 0) (Point 1 0) (Point 2 0) @?= Collinear
+    , testCase "turn: collinear with the third point closer counts as clockwise" $
+        turn (Point 0 0) (Point 2 0) (Point 1 0) @?= Clockwise
+    , testCase "turn: collinear with the third point farther counts as counter-clockwise" $
+        turn (Point 0 0) (Point 1 0) (Point 2 0) @?= CounterClockwise
+    , testProperty "turn never reports collinear" $
+        \a b c -> turn a b c /= Collinear
+    , testProperty "reversing the walk flips the orientation" $
+        \a b c -> orientation a b c === flipOrientation (orientation c b a)
+    ]
+  where
+    flipOrientation Clockwise = CounterClockwise
+    flipOrientation CounterClockwise = Clockwise
+    flipOrientation Collinear = Collinear
+
+edgeTests :: TestTree
+edgeTests =
+  testGroup
+    "Edge"
+    [ testProperty "mkEdge is symmetric" $
+        \a b -> mkEdge a b === mkEdge b a
+    , testProperty "equal edges hash equally" $
+        \a b -> hash (mkEdge a b) === hash (mkEdge b a)
+    , testCase "crossing segments intersect at the crossing" $
+        intersection (mkEdge (Point 0 0) (Point 2 2)) (mkEdge (Point 0 2) (Point 2 0)) @?= Just (Point 1 1)
+    , testCase "parallel segments do not intersect" $
+        intersection (mkEdge (Point 0 0) (Point 1 0)) (mkEdge (Point 0 1) (Point 1 1)) @?= Nothing
+    , testCase "segments sharing an endpoint do not intersect" $
+        intersection (mkEdge (Point 0 0) (Point 1 1)) (mkEdge (Point 0 0) (Point 1 0)) @?= Nothing
+    , testCase "segments on crossing lines but apart do not intersect" $
+        intersection (mkEdge (Point 0 0) (Point 1 1)) (mkEdge (Point 2 0) (Point 3 1)) @?= Nothing
+    , testCase "a segment touching the other at an endpoint does not intersect" $
+        intersection (mkEdge (Point 0 0) (Point 2 0)) (mkEdge (Point 1 0) (Point 1 1)) @?= Nothing
+    ]
+
+triangleTests :: TestTree
+triangleTests =
+  testGroup
+    "Triangle"
+    [ testProperty "mkTriangle ignores the order of the vertices" $
+        \a b c -> all (== mkTriangle a b c) [mkTriangle x y z | [x, y, z] <- permutations [a, b, c]]
+    , testCase "a point inside the circumcircle" $
+        isOutsideCircumcircle (Point 1 1) rightTriangle @?= False
+    , testCase "a point far away is outside the circumcircle" $
+        isOutsideCircumcircle (Point 5 5) rightTriangle @?= True
+    , testCase "a point on the circumcircle counts as outside" $
+        isOutsideCircumcircle (Point 2 2) rightTriangle @?= True
+    , testCase "isPointInTriangle: inside" $
+        isPointInTriangle (mkTriangle (Point 0 0) (Point 4 0) (Point 0 4)) (Point 1 1) @?= True
+    , testCase "isPointInTriangle: outside" $
+        isPointInTriangle (mkTriangle (Point 0 0) (Point 4 0) (Point 0 4)) (Point 3 3) @?= False
+    , testCase "isPointInTriangle: a vertex is not inside" $
+        isPointInTriangle (mkTriangle (Point 0 0) (Point 4 0) (Point 0 4)) (Point 4 0) @?= False
+    , testCase "triangles are sorted into and out of a polygon by their centroid" $
+        let ts =
+              [ mkTriangle (Point 1 1) (Point 2 1) (Point 1 2)
+              , mkTriangle (Point 10 10) (Point 11 10) (Point 10 11)
+              ]
+         in (trianglesInside square ts, trianglesOutside square ts) @?= splitAt 1 ts
+    ]
+  where
+    rightTriangle = mkTriangle (Point 0 0) (Point 2 0) (Point 0 2)
+
+-- | The unit square scaled to 4, wound clockwise.
+square :: Polygon
+square = Polygon (Point 0 0 :| [Point 0 4, Point 4 4, Point 4 0])
+
+polygonTests :: TestTree
+polygonTests =
+  testGroup
+    "Polygon"
+    [ testCase "hullOf3 starts at the lower-left point and winds clockwise" $
+        vertices (hullOf3 (Point 2 0) (Point 1 1) (Point 0 0)) @?= [Point 0 0, Point 1 1, Point 2 0]
+    , testProperty "hullOf3 is convex and clockwise" $
+        \a b c -> isConvex (hullOf3 a b c)
+    , testCase "hullOf4 of a square keeps all four corners" $
+        case hullOf4 (Point 2 2) (Point 0 0) (Point 2 0) (Point 0 2) of
+          Quadrilateral p1 p2 p3 p4 -> [p1, p2, p3, p4] @?= [Point 0 0, Point 0 2, Point 2 2, Point 2 0]
+          TriangleWithInner {} -> assertBool "expected a quadrilateral" False
+    , testCase "hullOf4 finds the inner point" $
+        case hullOf4 (Point 0 0) (Point 2 1) (Point 4 0) (Point 2 4) of
+          TriangleWithInner _ _ _ inner -> inner @?= Point 2 1
+          Quadrilateral {} -> assertBool "expected a triangle with an inner point" False
+    , testProperty "hullOf4 is convex and clockwise" $
+        \a b c d -> isConvex (hull4Polygon (hullOf4 a b c d))
+    , testCase "isPointInPolygon: inside" $ isPointInPolygon square (Point 1 1) @?= True
+    , testCase "isPointInPolygon: outside" $ isPointInPolygon square (Point 5 5) @?= False
+    , testCase "isPointInPolygon: a vertex is not inside" $ isPointInPolygon square (Point 4 4) @?= False
+    , testCase "a square is convex" $ isConvex square @?= True
+    , testCase "a polygon with a reflex vertex is not convex" $
+        isConvex (Polygon (Point 0 0 :| [Point 0 3, Point 3 3, Point 1 2])) @?= False
+    , testCase "polygonEdges closes the ring" $
+        length (polygonEdges square) @?= 4
+    ]
diff --git a/test/Test/Mesh.hs b/test/Test/Mesh.hs
new file mode 100644
--- /dev/null
+++ b/test/Test/Mesh.hs
@@ -0,0 +1,67 @@
+module Test.Mesh (tests) where
+
+import Data.List (nub, sort)
+import Data.List.NonEmpty (NonEmpty (..))
+import Data.Vector qualified as V
+import Test.Geometry (Points (..))
+import Test.Tasty (TestTree, localOption, testGroup)
+import Test.Tasty.HUnit (testCase, (@?=))
+import Test.Tasty.QuickCheck (Property, QuickCheckTests (..), counterexample, testProperty)
+import Triangulation
+
+tests :: TestTree
+tests =
+  localOption (QuickCheckTests 300) $
+    testGroup
+      "Mesh"
+      [ testProperty "nodes are the distinct vertices in ascending order" $
+          meshHolds $ \ts (Mesh ns _) ->
+            V.toList ns == sort (nub (concatMap (\(Triangle a b c) -> [a, b, c]) ts))
+      , testProperty "elements index the nodes, ascending, no repeats" $
+          meshHolds $ \_ (Mesh ns es) ->
+            let n = V.length ns
+                triples = V.toList es
+             in all (\(i, j, k) -> all (\x -> 0 <= x && x < n) [i, j, k]) triples
+                  && and (zipWith (<) triples (drop 1 triples))
+      , testProperty "every element is wound counter-clockwise" $
+          meshHolds $ \_ m@(Mesh ns _) ->
+            all
+              (\(i, j, k) -> orientation (ns V.! i) (ns V.! j) (ns V.! k) == CounterClockwise)
+              (V.toList (elements m))
+      , testProperty "fromMesh . toMesh is the identity on the set of triangles" $
+          meshHolds $
+            \ts m -> sort (fromMesh m) == sort (nub ts)
+      , testProperty "the boundary of a Delaunay triangulation is its hull" $
+          \(Points points) -> case triangulate points of
+            Nothing -> counterexample "triangulate returned Nothing" False
+            Just t ->
+              let m = toMesh (triangles t)
+               in counterexample (show m) $
+                    edgesOf m (boundaryEdges m) == sort (polygonEdges (hull t))
+      , testProperty "the boundary of a polygon with a hole is both polygons" $
+          \(Points points) -> case constrainedTriangulate canvas [hole] points of
+            Nothing -> counterexample "constrainedTriangulate returned Nothing" False
+            Just ts ->
+              let m = toMesh ts
+               in counterexample (show m) $
+                    edgesOf m (boundaryEdges m) == sort (polygonEdges canvas ++ polygonEdges hole)
+      , testCase "a single triangle" $
+          toMesh [mkTriangle (Point 0 0) (Point 1 0) (Point 0 1)]
+            @?= Mesh (V.fromList [Point 0 0, Point 0 1, Point 1 0]) (V.fromList [(0, 2, 1)])
+      , testCase "no triangles" $
+          toMesh [] @?= Mesh V.empty V.empty
+      ]
+  where
+    canvas = Polygon (Point 0 0 :| [Point 0 2000, Point 2000 2000, Point 2000 0])
+    hole = Polygon (Point 600 500 :| [Point 600 1300, Point 1400 1300, Point 1400 500])
+    edgesOf (Mesh ns _) = sort . map (\(i, j) -> mkEdge (ns V.! i) (ns V.! j))
+
+-- | The property must hold for the mesh of every Delaunay triangulation of a
+-- random point set.
+meshHolds :: ([Triangle] -> Mesh -> Bool) -> Points -> Property
+meshHolds p (Points points) = case triangulate points of
+  Nothing -> counterexample "triangulate returned Nothing" False
+  Just t ->
+    let ts = triangles t
+        m = toMesh ts
+     in counterexample (show m) (p ts m)
diff --git a/test/Test/Refine.hs b/test/Test/Refine.hs
new file mode 100644
--- /dev/null
+++ b/test/Test/Refine.hs
@@ -0,0 +1,221 @@
+module Test.Refine (tests) where
+
+import Data.List (sort)
+import Data.List.NonEmpty (NonEmpty (..))
+import Data.Maybe (isNothing)
+import Data.Vector qualified as V
+import Test.Geometry (Points (..))
+import Test.Tasty (TestTree, localOption, testGroup)
+import Test.Tasty.HUnit (assertBool, testCase, (@?=))
+import Test.Tasty.QuickCheck (Property, QuickCheckTests (..), counterexample, testProperty, (===))
+import Triangulation
+import Triangulation.Check (hasNoIntersections, isLocallyDelaunay)
+import Triangulation.Generator (Rectangle (..), generatePoints)
+
+tests :: TestTree
+tests =
+  testGroup
+    "Refine"
+    [ localOption (QuickCheckTests 100) properties
+    , shapes
+    , edgeCases
+    ]
+
+properties :: TestTree
+properties =
+  testGroup
+    "random interior points in a square with a hole"
+    [ testProperty "20°: every angle is at least 20° and the mesh is valid" $
+        refined defaultQuality [canvas, hole] $
+          \ts -> wellRefined defaultQuality [canvas, hole] ts
+    , testProperty "20° and an area bound: no triangle is larger than the bound" $
+        let q = defaultQuality {maxArea = Just 20000}
+         in refined q [canvas, hole] $ \ts -> wellRefined q [canvas, hole] ts
+    , testProperty "the input points survive refinement" $
+        \(Points points) -> case constrainedTriangulate canvas [hole] points >>= refine defaultQuality of
+          Nothing -> counterexample "Nothing" False
+          Just ts ->
+            let ns = V.toList (nodes (toMesh ts))
+                kept = filter (\p -> strictlyInside canvas p && not (strictlyInside hole p)) points
+             in counterexample (show ns) $ all (`elem` ns) (kept ++ vertices canvas ++ vertices hole)
+    , testProperty "refining twice changes nothing" $
+        \(Points points) -> case constrainedTriangulate canvas [hole] points >>= refine defaultQuality of
+          Nothing -> counterexample "Nothing" False
+          Just ts -> fmap sort (refine defaultQuality ts) === Just (sort ts)
+    , testProperty "the Delaunay triangulation of a point set (hull as boundary) refines too" $
+        -- the corners keep the hull free of acute angles, which Ruppert's algorithm cannot fix
+        \(Points points) -> case triangulate (vertices canvas ++ points) of
+          Nothing -> counterexample "triangulate returned Nothing" False
+          Just t -> case refine defaultQuality (triangles t) of
+            Nothing -> counterexample "refine returned Nothing" False
+            Just ts -> counterexample (show ts) (wellRefined defaultQuality [hull t] ts)
+    ]
+
+-- | The refinement of the constrained triangulation of random points must
+-- succeed and satisfy the predicate.
+refined :: Quality -> [Polygon] -> ([Triangle] -> Bool) -> Points -> Property
+refined quality (boundary : holes) p (Points points) =
+  case constrainedTriangulate boundary holes points of
+    Nothing -> counterexample "constrainedTriangulate returned Nothing" False
+    Just coarse -> case refine quality coarse of
+      Nothing -> counterexample "refine returned Nothing" False
+      Just ts -> counterexample (show ts) (p ts)
+refined _ [] _ _ = counterexample "no boundary" False
+
+shapes :: TestTree
+shapes =
+  testGroup
+    "particular shapes"
+    [ testCase "an L-shaped region at 30°" $
+        check (Just q30) [lShape] (constrainedTriangulate lShape [] (grid 300) >>= refine q30)
+    , testCase "a square with a hole at 30° and an area bound" $
+        let q = q30 {maxArea = Just 50000}
+         in check (Just q) [canvas, hole] (constrainedTriangulate canvas [hole] (grid 500) >>= refine q)
+    , testCase "a wedge with a 17° input angle terminates and keeps its boundary" $
+        case constrainedTriangulate wedge [] [] >>= refine defaultQuality of
+          Nothing -> assertBool "refine returned Nothing" False
+          Just ts -> do
+            problems Nothing [wedge] ts @?= []
+            -- every angle away from the wedge apex meets the bound
+            assertBool "angles away from the apex" $
+              all (\t -> smallestAngle t >= 20 - tolerance || touches (Point 0 0) t) ts
+    , testCase "a wedge with a 17° input angle at 30° with an area bound" $
+        check Nothing [wedge] (constrainedTriangulate wedge [] [] >>= refine q30 {maxArea = Just 5000})
+    , testCase "a thin rectangle at 30°" $
+        let thin = Polygon (Point 0 0 :| [Point 0 100, Point 2000 100, Point 2000 0])
+         in check (Just q30) [thin] (constrainedTriangulate thin [] [] >>= refine q30)
+    , testCase "a hexagon with a triangular hole"
+        $ check (Just defaultQuality) [hexagon, triangleHole]
+        $ constrainedTriangulate hexagon [triangleHole] [] >>= refine defaultQuality
+    , testCase "1000 random points at 25°" $
+        let q = defaultQuality {minAngle = 25}
+            points = generatePoints 3 1000 Rectangle {minCorner = Point 0 0, maxCorner = Point 2000 2000}
+         in check (Just q) [canvas, hole] (constrainedTriangulate canvas [hole] points >>= refine q)
+    ]
+  where
+    q30 = defaultQuality {minAngle = 30}
+    grid h = [Point x y | x <- [h, 2 * h .. 2000 - h], y <- [h, 2 * h .. 2000 - h]]
+    check _ _ Nothing = assertBool "refine returned Nothing" False
+    check quality polygons (Just ts) = problems quality polygons ts @?= []
+
+edgeCases :: TestTree
+edgeCases =
+  testGroup
+    "edge cases"
+    [ testCase "no triangles" $ refine defaultQuality [] @?= Just []
+    , testCase "a triangle that already meets the bound is returned as is" $
+        let t = mkTriangle (Point 0 0) (Point 100 0) (Point 50 80)
+         in refine defaultQuality [t] @?= Just [t]
+    , testCase "no budget and a poor triangle give Nothing" $
+        let t = mkTriangle (Point 0 0) (Point 1000 0) (Point 500 10)
+         in assertBool "expected Nothing" $ isNothing (refine defaultQuality {maxInsertions = 0} [t])
+    , testCase "no budget and a good triangle succeed" $
+        let t = mkTriangle (Point 0 0) (Point 100 0) (Point 50 80)
+         in refine defaultQuality {maxInsertions = 0} [t] @?= Just [t]
+    , testCase "an impossible bound exhausts the budget rather than looping"
+        $
+        -- a square cut by its diagonal: two 45° angles at interior edges, which no refinement lifts to 59°
+        assertBool "expected Nothing"
+        $ isNothing (constrainedTriangulate canvas [] [] >>= refine (Quality 59 Nothing 200))
+    , testCase "an input angle is not asked to meet the bound" $
+        let t = mkTriangle (Point 0 0) (Point 100 0) (Point 50 80)
+         in refine (Quality 59 Nothing 200) [t] @?= Just [t]
+    , testCase "no bounds at all: only encroached segments are split" $
+        -- the apex sees the long edge at 118°, so it encroaches upon it
+        let t = mkTriangle (Point 0 0) (Point 1000 0) (Point 500 300)
+         in case refine (Quality 0 Nothing 100) [t] of
+              Nothing -> assertBool "Nothing" False
+              Just ts -> do
+                problems Nothing [Polygon (Point 0 0 :| [Point 500 300, Point 1000 0])] ts @?= []
+                assertBool "something was split" (length ts > 1)
+    ]
+
+-- Shapes; polygons are wound clockwise (y up: interior on the right).
+
+canvas, hole, lShape, wedge, hexagon, triangleHole :: Polygon
+canvas = Polygon (Point 0 0 :| [Point 0 2000, Point 2000 2000, Point 2000 0])
+hole = Polygon (Point 600 500 :| [Point 600 1300, Point 1400 1300, Point 1400 500])
+lShape =
+  Polygon
+    (Point 0 0 :| [Point 0 2000, Point 1000 2000, Point 1000 1000, Point 2000 1000, Point 2000 0])
+wedge = Polygon (Point 0 0 :| [Point 0 1000, Point 300 1000])
+hexagon =
+  Polygon
+    (Point 1000 0 :| [Point 134 500, Point 134 1500, Point 1000 2000, Point 1866 1500, Point 1866 500])
+triangleHole = Polygon (Point 700 700 :| [Point 1000 1300, Point 1300 700])
+
+tolerance :: Double
+tolerance = 1e-9
+
+-- | The mesh satisfies the bounds and is a valid triangulation of the region.
+wellRefined :: Quality -> [Polygon] -> [Triangle] -> Bool
+wellRefined quality polygons ts = null (problems (Just quality) polygons ts)
+
+-- | What is wrong with the mesh: bounds not met (when a quality is given), or
+-- not a constrained Delaunay triangulation of exactly the region bounded by
+-- the polygons (the first is the outer boundary, the rest are holes).
+problems :: Maybe Quality -> [Polygon] -> [Triangle] -> [String]
+problems quality polygons ts =
+  [ name
+  | (name, ok) <-
+      [
+        ( "an angle below the bound"
+        , all (\q -> all (\t -> smallestAngle t >= minAngle q - tolerance) ts) quality
+        )
+      ,
+        ( "an area above the bound"
+        , all (\a -> all (\t -> triangleArea t <= a) ts) (quality >>= maxArea)
+        )
+      , ("no triangles", not (null ts))
+      , ("not locally Delaunay", isLocallyDelaunay ts)
+      , ("crossing edges", hasNoIntersections ts)
+      , ("a flat triangle", all (\(Triangle a b c) -> orientation a b c /= Collinear) ts)
+      , ("a boundary edge off the polygons", boundaryOnPolygons)
+      , ("a polygon vertex missing", all (all (`elem` V.toList (nodes m)) . vertices) polygons)
+      , ("Euler's formula", eulerWithHoles)
+      , ("an encroached segment", noEncroachedSegments m)
+      ]
+  , not ok
+  ]
+  where
+    m = toMesh ts
+    boundary = boundaryEdges m
+    -- every boundary edge of the mesh lies on an edge of one of the polygons
+    boundaryOnPolygons = all onSomePolygonEdge boundary
+    onSomePolygonEdge (i, j) =
+      let (p, q) = (nodes m V.! i, nodes m V.! j)
+       in any (any (\(Edge u v) -> between u v p && between u v q) . polygonEdges) polygons
+    -- split points are computed in floating point, so "on the edge" is up to rounding
+    between u v p =
+      let (dx, dy) = (px v - px u, py v - py u)
+          len2 = dx * dx + dy * dy
+          t = ((px p - px u) * dx + (py p - py u) * dy) / len2
+          cross = dx * (py p - py u) - dy * (px p - px u)
+       in -1e-9 <= t && t <= 1 + 1e-9 && cross * cross <= 1e-18 * len2 * len2
+    -- for a triangulated region with V vertices, B boundary vertices and H holes: T = 2V - B + 2H - 2
+    eulerWithHoles =
+      let v = V.length (nodes m)
+          b = length boundary -- boundary loops have as many edges as vertices
+          h = length polygons - 1
+       in V.length (elements m) == 2 * v - b + 2 * h - 2
+
+-- | No vertex lies strictly inside the diametral circle of a boundary edge.
+noEncroachedSegments :: Mesh -> Bool
+noEncroachedSegments m =
+  and
+    [ not (encroaches p (nodes m V.! i) (nodes m V.! j))
+    | (i, j) <- boundaryEdges m
+    , p <- V.toList (nodes m)
+    ]
+  where
+    encroaches (Point x y) (Point x1 y1) (Point x2 y2) = (x - x1) * (x - x2) + (y - y1) * (y - y2) < 0
+
+touches :: Point -> Triangle -> Bool
+touches p (Triangle a b c) = p `elem` [a, b, c]
+
+-- | Strictly inside an axis-aligned rectangle given as a polygon.
+strictlyInside :: Polygon -> Point -> Bool
+strictlyInside polygon (Point x y) =
+  let xs = map px (vertices polygon)
+      ys = map py (vertices polygon)
+   in minimum xs < x && x < maximum xs && minimum ys < y && y < maximum ys
diff --git a/test/Test/Ring.hs b/test/Test/Ring.hs
new file mode 100644
--- /dev/null
+++ b/test/Test/Ring.hs
@@ -0,0 +1,50 @@
+module Test.Ring (tests) where
+
+import Data.List.NonEmpty (NonEmpty (..))
+import Test.Tasty (TestTree, testGroup)
+import Test.Tasty.HUnit (testCase, (@?=))
+import Triangulation.Geometry.Ring (
+  arc,
+  cyclicPairs,
+  cyclicTriples,
+  predecessor,
+  splitLoop,
+  successor,
+ )
+
+ring :: NonEmpty Int
+ring = 1 :| [2, 3, 4]
+
+tests :: TestTree
+tests =
+  testGroup
+    "Ring"
+    [ testCase "cyclicPairs wraps around" $
+        cyclicPairs ring @?= [(1, 2), (2, 3), (3, 4), (4, 1)]
+    , testCase "cyclicPairs of a singleton pairs it with itself" $
+        cyclicPairs (7 :| [] :: NonEmpty Int) @?= [(7, 7)]
+    , testCase "cyclicTriples gives predecessor and successor" $
+        cyclicTriples ring @?= [(4, 1, 2), (1, 2, 3), (2, 3, 4), (3, 4, 1)]
+    , testCase "successor of the last element is the first" $
+        successor 4 ring @?= Just 1
+    , testCase "predecessor of the first element is the last" $
+        predecessor 1 ring @?= Just 4
+    , testCase "neighbours of an absent element" $
+        (successor 9 ring, predecessor 9 ring) @?= (Nothing, Nothing)
+    , testCase "arc forward" $
+        arc 2 4 ring @?= Just (2 :| [3, 4])
+    , testCase "arc wrapping around" $
+        arc 4 2 ring @?= Just (4 :| [1, 2])
+    , testCase "arc from a point to itself goes all the way round" $
+        arc 3 3 ring @?= Just (3 :| [4, 1, 2, 3])
+    , testCase "arc from an absent point" $
+        arc 9 1 ring @?= Nothing
+    , testCase "arc to an absent point goes all the way round" $
+        arc 1 9 ring @?= Just (1 :| [2, 3, 4, 1])
+    , testCase "splitLoop separates the inner loop" $
+        splitLoop 2 (1 :| [2, 5, 6, 2, 3] :: NonEmpty Int) @?= (1 :| [2, 3], [2, 5, 6])
+    , testCase "splitLoop without a loop" $
+        splitLoop 2 ring @?= (ring, [])
+    , testCase "splitLoop on an absent element" $
+        splitLoop 9 ring @?= (ring, [])
+    ]
diff --git a/test/Test/Triangulation.hs b/test/Test/Triangulation.hs
new file mode 100644
--- /dev/null
+++ b/test/Test/Triangulation.hs
@@ -0,0 +1,314 @@
+module Test.Triangulation (tests) where
+
+import Data.List (nub, sort)
+import Data.List.NonEmpty (NonEmpty (..))
+import Data.Maybe (isNothing)
+import Data.Vector qualified as V
+import Test.Geometry (Points (..))
+import Test.Tasty (TestTree, localOption, testGroup)
+import Test.Tasty.HUnit (assertBool, testCase, (@?=))
+import Test.Tasty.QuickCheck (
+  Arbitrary (..),
+  Property,
+  QuickCheckTests (..),
+  chooseInt,
+  counterexample,
+  property,
+  shrinkList,
+  testProperty,
+  vectorOf,
+ )
+import Triangulation
+import Triangulation.Check (hasNoIntersections, isDelaunay, isLocallyDelaunay)
+import Triangulation.Constrained (forceEdges)
+import Triangulation.Generator (Rectangle (..), generatePoints)
+import Triangulation.Store qualified as Store
+
+tests :: TestTree
+tests =
+  testGroup
+    "Triangulation"
+    [ localOption (QuickCheckTests 300) properties
+    , localOption (QuickCheckTests 300) collinearHeavy
+    , edgeCases
+    , largeSets
+    , constrained
+    , localOption (QuickCheckTests 300) constrainedRegion
+    ]
+
+-- | Every triangulation of a random point set must satisfy the predicate.
+holds :: String -> (Triangulation -> Bool) -> Points -> Property
+holds what p (Points points) = case triangulate points of
+  Nothing -> counterexample "triangulate returned Nothing" False
+  Just t -> counterexample (what ++ " violated by " ++ show t) (p t)
+
+properties :: TestTree
+properties =
+  testGroup
+    "properties of random point sets"
+    [ testProperty "Delaunay condition" $ holds "Delaunay condition" (isDelaunay . triangles)
+    , testProperty "no triangle has zero area" $ holds "non-degeneracy" (not . any isFlat . triangles)
+    , testProperty "no edge spans another vertex" $ holds "no edge spans a vertex" noEdgeSpansAVertex
+    , testProperty "no two edges cross" $ holds "planarity" (hasNoIntersections . triangles)
+    , testProperty "Euler's formula: 2n - h - 2 triangles" $ holds "Euler's formula" hasEulerTriangleCount
+    , testProperty "the hull is convex" $ holds "convexity" (isConvex . hull)
+    , testProperty "the hull contains every point" $
+        \ps@(Points points) -> holds "hull containment" (hullContains points) ps
+    ]
+
+-- | Points drawn from a coarse grid, so that collinear and cocircular triples
+-- are the rule rather than the exception.
+collinearHeavy :: TestTree
+collinearHeavy =
+  testGroup
+    "collinear-heavy input (grid points)"
+    [ testProperty "Delaunay condition" $ gridHolds (isDelaunay . triangles)
+    , testProperty "no two edges cross" $ gridHolds (hasNoIntersections . triangles)
+    , testProperty "no triangle has zero area" $ gridHolds (not . any isFlat . triangles)
+    , testProperty "no edge spans another vertex" $ gridHolds noEdgeSpansAVertex
+    , testProperty "Euler's formula" $ gridHolds hasEulerTriangleCount
+    , testProperty "the hull contains every point" $
+        \(GridPoints points) -> gridHolds (hullContains points) (GridPoints points)
+    ]
+  where
+    -- all points on one line: no triangulation exists and Nothing is the right answer
+    gridHolds p (GridPoints points) = case triangulate points of
+      Nothing -> counterexample "triangulate returned Nothing" (allCollinear points)
+      Just t -> counterexample (show (triangles t)) (p t)
+    allCollinear (a : b : rest) = all (\c -> orientation a b c == Collinear) rest
+    allCollinear _ = True
+
+-- | At least three distinct points with coordinates that are multiples of 100;
+-- shrinking drops points but never below three.
+newtype GridPoints = GridPoints [Point]
+  deriving (Show)
+
+instance Arbitrary GridPoints where
+  arbitrary = do
+    n <- chooseInt (3, 12)
+    points <- vectorOf (3 * n) (Point <$> coordinate <*> coordinate)
+    pure (GridPoints (take n (nub points)))
+    where
+      coordinate = (* 100) . fromIntegral <$> chooseInt (0, 20)
+  shrink (GridPoints ps) = [GridPoints ps' | ps' <- shrinkList (const []) ps, length ps' >= 3]
+
+isFlat :: Triangle -> Bool
+isFlat (Triangle a b c) = orientation a b c == Collinear
+
+-- | No edge of the triangulation has a vertex strictly inside it. Overlapping
+-- collinear edges used to slip past 'hasNoIntersections', which only looks for
+-- proper crossings.
+noEdgeSpansAVertex :: Triangulation -> Bool
+noEdgeSpansAVertex triangulation =
+  not $ or [spans a b p | Edge a b <- Store.edges (triangleStore triangulation), p <- points]
+  where
+    points = nub $ concatMap (\(Triangle a b c) -> [a, b, c]) (triangles triangulation)
+    spans a b p =
+      p /= a
+        && p /= b
+        && orientation a b p == Collinear
+        && min (px a) (px b) <= px p
+        && px p <= max (px a) (px b)
+        && min (py a) (py b) <= py p
+        && py p <= max (py a) (py b)
+
+edgeCases :: TestTree
+edgeCases =
+  testGroup
+    "edge cases"
+    [ testCase "fewer than three points give Nothing"
+        $ assertBool "expected Nothing"
+        $ all (isNothing . triangulate) [[], [Point 1 1], [Point 1 1, Point 2 2]]
+    , testCase "five points (split 2 + 3) and 9-11 points (split 3 + rest)" $
+        mapM_
+          (\n -> assertBool ("size " ++ show n) (all wellFormed (deterministicSets n 200)))
+          ([5 .. 12] :: [Int])
+    ]
+
+largeSets :: TestTree
+largeSets =
+  testGroup
+    "large deterministic sets"
+    [ testCase "1000 points, three seeds"
+        $ assertBool "well-formed and Delaunay"
+        $ all
+          (any (\t -> hasEulerTriangleCount t && isDelaunay (triangles t)) . triangulate)
+          (deterministicSets 1000 3)
+    , -- A lattice is the worst case for degeneracy and the most likely input of
+      -- a finite element model: every row, column and diagonal is collinear and
+      -- the corners of every cell are cocircular. Repairing the flat triangles
+      -- it produces used to rescan the whole store per repair, which made a
+      -- lattice of this size quadratic.
+      testCase "a 40x40 lattice" $ case triangulate lattice of
+        Nothing -> assertBool "expected a triangulation" False
+        Just t -> do
+          let ts = triangles t
+          assertBool "no flat triangle" (not (any isFlat ts))
+          assertBool "locally Delaunay" (isLocallyDelaunay ts)
+          assertBool "no edge spans a vertex" (noEdgeSpansAVertex t)
+          assertBool "Euler's formula" (hasEulerTriangleCount t)
+          -- not 'isConvex': the hull ring of a lattice carries the boundary
+          -- points that lie between its corners, and under the symbolic
+          -- perturbation those are not collinear, so a strictly convex ring is
+          -- the wrong expectation here. The property tests on grid points omit
+          -- the same check for the same reason.
+          assertBool "the hull contains every point" (hullContains lattice t)
+          length (nub (concatMap (\(Triangle a b c) -> [a, b, c]) ts)) @?= length lattice
+    ]
+  where
+    lattice =
+      [ Point (fromIntegral i * 300) (fromIntegral j * 300)
+      | i <- [0 .. 39 :: Int]
+      , j <- [0 .. 39 :: Int]
+      ]
+
+constrained :: TestTree
+constrained =
+  testGroup
+    "constrained edges"
+    [ testProperty "a forced edge is present and nothing crosses" $
+        \(Points points) -> case (points, triangulate points) of
+          (a : b : _, Just (Triangulation _ store)) ->
+            let edge = mkEdge a b
+                store' = forceEdges store [edge | not (Store.member edge store)] [edge]
+             in property $ Store.member edge store' && hasNoIntersections (Store.triangles store')
+          _ -> counterexample "triangulate returned Nothing" False
+    , testProperty "constrainedTriangulate keeps the boundary and empties the hole" $
+        \(Points points) -> constrainedHolds hole points
+    , testProperty "constrainedTriangulate with the diagonal hole (collinear input)" $
+        \(Points points) -> constrainedHolds diagonalHole points
+    , -- Regression: hole corners on the canvas diagonals make many triples collinear;
+      -- a zero-area candidate triangle used to slip in and swallow a hole edge.
+      testCase "hole with corners on the canvas diagonals (collinear input)"
+        $ assertBool "boundary and hole edges present, hole empty"
+        $ constrainedOk
+          diagonalHole
+          [ Point 823.1104570397463 1295.569301574011
+          , Point 1705.4929691922034 740.4819744884217
+          , Point 97.24513811803548 174.33490747208324
+          ]
+    ]
+  where
+    canvas = Polygon (Point 0 0 :| [Point 0 2000, Point 2000 2000, Point 2000 0])
+    -- corners off the canvas diagonals, so the input is in general position
+    hole = Polygon (Point 600 500 :| [Point 600 1300, Point 1400 1300, Point 1400 500])
+    diagonalHole = Polygon (Point 500 500 :| [Point 500 1500, Point 1500 1500, Point 1500 500])
+    constrainedHolds h points = case constrainedTriangulate canvas [h] points of
+      Nothing -> counterexample "constrainedTriangulate returned Nothing" False
+      Just ts -> counterexample (show ts) (constrainedOkWith h ts)
+    constrainedOk h points = any (constrainedOkWith h) (constrainedTriangulate canvas [h] points)
+    constrainedOkWith h ts =
+      let edges = concatMap triangleEdges ts
+       in all (`elem` edges) (polygonEdges canvas ++ polygonEdges h)
+            && all (`isTriangleInPolygon` canvas) ts
+            && not (any (`isTriangleInPolygon` h) ts)
+
+-- | The constrained triangulation of a square with a hole, for interior points
+-- in general position and for lattice points (which land on the polygon edges,
+-- inside the hole, and on top of the corners).
+constrainedRegion :: TestTree
+constrainedRegion =
+  testGroup
+    "constrainedTriangulate as a mesh of the region"
+    [ testProperty "random points: valid region mesh" $ \(Points points) -> regionHolds hole points
+    , testProperty "lattice points: valid region mesh" $ \(GridPoints points) -> regionHolds latticeHole points
+    , testProperty "lattice points on the boundary subdivide it" $ \(GridPoints points) ->
+        case constrainedTriangulate canvas [latticeHole] points of
+          Nothing -> counterexample "Nothing" False
+          Just ts ->
+            let m = toMesh ts
+                onBoundary =
+                  [ p
+                  | p <- nub points
+                  , any (\(Edge u v) -> between u v p) (polygonEdges canvas ++ polygonEdges latticeHole)
+                  ]
+             in counterexample (show m) $ all (`elem` V.toList (nodes m)) onBoundary
+    , testCase "an empty square: two triangles" $
+        fmap length (constrainedTriangulate canvas [] []) @?= Just 2
+    , testCase "a square with a hole and no interior points: eight triangles" $
+        fmap length (constrainedTriangulate canvas [hole] []) @?= Just 8
+    , testCase "points inside the hole disappear" $
+        fmap
+          (length . nodes . toMesh)
+          (constrainedTriangulate canvas [hole] [Point 1000 900, Point 900 1000])
+          @?= Just 8
+    , -- the corner itself belongs to no triangle of the region and disappears
+      testCase "a hole touching the boundary at a corner cuts the corner off" $
+        let cornerHole = Polygon (Point 0 0 :| [Point 0 500, Point 500 500, Point 500 0])
+         in case constrainedTriangulate canvas [cornerHole] [] of
+              Nothing -> assertBool "Nothing" False
+              Just ts -> do
+                length ts @?= 4
+                assertBool "locally Delaunay and planar" (isLocallyDelaunay ts && hasNoIntersections ts)
+                sort (V.toList (nodes (toMesh ts)))
+                  @?= sort [Point 0 500, Point 0 2000, Point 2000 2000, Point 2000 0, Point 500 0, Point 500 500]
+                length (boundaryEdges (toMesh ts)) @?= 6
+    , testCase "two holes" $
+        let hole2 = Polygon (Point 1500 1500 :| [Point 1500 1900, Point 1900 1900, Point 1900 1500])
+         in assertBool "valid" $
+              any
+                (regionOkWith [hole, hole2])
+                (constrainedTriangulate canvas [hole, hole2] (concat (deterministicSets 40 1)))
+    ]
+  where
+    canvas = Polygon (Point 0 0 :| [Point 0 2000, Point 2000 2000, Point 2000 0])
+    hole = Polygon (Point 600 500 :| [Point 600 1300, Point 1400 1300, Point 1400 500])
+    latticeHole = Polygon (Point 600 600 :| [Point 600 1400, Point 1400 1400, Point 1400 600])
+    regionHolds h points = case constrainedTriangulate canvas [h] points of
+      Nothing -> counterexample "constrainedTriangulate returned Nothing" False
+      Just ts -> counterexample (show ts) (regionOk h ts)
+    regionOk h = regionOkWith [h]
+    regionOkWith holes ts =
+      let m = toMesh ts
+          polygons = canvas : holes
+          boundary = boundaryEdges m
+          onPolygonEdge (i, j) =
+            let (p, q) = (nodes m V.! i, nodes m V.! j)
+             in any (any (\(Edge u v) -> between u v p && between u v q) . polygonEdges) polygons
+          nodeCount = V.length (nodes m)
+          b = length boundary
+       in not (null ts)
+            && isLocallyDelaunay ts
+            && hasNoIntersections ts
+            && not (any isFlat ts)
+            && all onPolygonEdge boundary
+            && all (all (`elem` V.toList (nodes m)) . vertices) polygons
+            && all (`isTriangleInPolygon` canvas) ts
+            && not (any (\t -> any (isTriangleInPolygon t) holes) ts)
+            && length ts == 2 * nodeCount - b + 2 * length holes - 2
+    between u v p =
+      orientation u v p == Collinear
+        && min (px u) (px v) <= px p
+        && px p <= max (px u) (px v)
+        && min (py u) (py v) <= py p
+        && py p <= max (py u) (py v)
+
+-- | For n points in general position with h of them on the convex hull, any
+-- triangulation has exactly 2n - h - 2 triangles. This guards against a
+-- triangulation that silently drops triangles: the Delaunay and intersection
+-- checks are vacuously true on an empty triangulation.
+hasEulerTriangleCount :: Triangulation -> Bool
+hasEulerTriangleCount triangulation@(Triangulation (Polygon hullPoints) _) =
+  let ts = triangles triangulation
+      n = length (nub (concatMap (\(Triangle p1 p2 p3) -> [p1, p2, p3]) ts))
+      h = length hullPoints
+   in length ts == 2 * n - h - 2
+
+hullContains :: [Point] -> Triangulation -> Bool
+hullContains points (Triangulation polygon@(Polygon hullPoints) _) =
+  all (\p -> p `elem` hullPoints || isPointInPolygon polygon p) points
+
+wellFormed :: [Point] -> Bool
+wellFormed points = case triangulate points of
+  Nothing -> False
+  Just t ->
+    isDelaunay (triangles t)
+      && hasNoIntersections (triangles t)
+      && hasEulerTriangleCount t
+      && isConvex (hull t)
+
+deterministicSets :: Int -> Int -> [[Point]]
+deterministicSets size seeds =
+  [ generatePoints seed size Rectangle {minCorner = Point 0 0, maxCorner = Point 2000 2000}
+  | seed <- [1 .. seeds]
+  ]