# moonlight-triangulation
> Part of **Moonlight**, the sheaf-theoretic computation layer beneath
> [Melusine](https://bluerose.blue) and Pale Meridian.
`moonlight-triangulation` carries Delaunay and constrained Delaunay
triangulations as a lawful finite-set algebra under canonical observation: a
mesh represents its site set, a join returns a valid Delaunay representative,
and the result is a triangulation again — so the operations close, compose, and
fold. Operations return typed obstructions where the finite arena cannot
represent a result.
Delaunay triangulation, constrained Delaunay (CDT), the Voronoi dual,
natural-neighbour interpolation, Ruppert refinement, walk point location,
convex hull, exact Shewchuk predicates, incremental insertion and removal,
versioned binary serialization.
## Operations
| Operation | Role | What it returns |
| --- | --- | --- |
| `union` | join | A valid Delaunay representative of both site sets; overlapping annotations glue through `JoinSemilattice`. |
| `unions` | balanced fold | The same join over a list, associated for the tournament rather than the left. |
| `siteRelation` | support order | Exact equality, proper-subset, disjointness or partial-overlap classification. |
| `intersection` | geometry-only meet | The sites both meshes hold. |
| `intersectionWith` | annotated meet | Shared sites with a left-then-right annotation combiner. |
| `difference` | relative complement | The left's sites and annotations, less the right's support. |
| `symmetricDifference` | exclusive or | The sites and annotations carried by exactly one mesh. |
| `canonicalize` | physical normal form | Construction-independent dense numbering, derived explicitly when required. |
| `refine` | quality | Steiner insertion, composed *after* any of the above — never a second kind of mesh. |
| `constrainedDelaunay` | boundary | Sites plus segments, after which interiority is computable. |
### Choose the operation
| What you need | Use | Do not substitute |
| --- | --- | --- |
| Build a geometry-only mesh from coordinates | `delaunayGeometry` | Do not manufacture unit annotations and erase them afterward. |
| Build from separate coordinates and annotations | `delaunayFromCoordinates` | Do not fabricate a `HasPosition` instance for an annotation that does not own geometry. |
| Combine two unconstrained meshes | `union` | Do not concatenate vertices and rebuild manually; `union` owns overlap annotations and the measured schedule. |
| Combine many unconstrained meshes | `unions` | Do not left-fold `union`; `unions` owns balanced association. |
| Classify coordinate support without constructing a mesh | `siteRelation` | Do not compare vertex counts or resident `Eq`; neither answers support order. |
| Keep coordinates present in both geometry-only meshes | `intersection` | Use `intersectionWith` instead when annotations must survive or be recomputed. |
| Keep shared coordinates and choose their annotation | `intersectionWith combine` | For a left annotation-preserving mask, use `intersectionWith const left mask`. |
| Remove a coordinate mask from a mesh | `difference source mask` | Do not remove by stale `VertexId`; removal compacts arenas while coordinates remain stable. |
| Keep coordinates present in exactly one mesh | `symmetricDifference` | Do not spell it as two differences plus `union`; the direct operation owns the persistent toggle schedule. |
| Build the first constrained mesh | `constrainedDelaunay` | Do not build unconstrained and treat a rendered outline as topology. |
| Canonically combine arbitrary constrained meshes | `unionConstrainedWith` or geometry-only `unionConstrained` | Do not use unconstrained `union`; constrained union owns complete conflict witnesses and constraint recovery. |
| Join two strictly x-separated constrained meshes | `joinSeparatedConstrainedWith` | Use the general constrained union when separation is not proved. |
| Append one constrained section to an authoritative constrained base | `extendConstrainedWith` | Do not use symmetric constrained union when base identity and already-solved constraints must remain resident. |
| Insert one site persistently | `BulkLoad.insert` or `BulkLoad.insertAt` | Do not open a manual transaction; the singleton entry owns the dense/copy-on-write crossover. |
| Insert a vector of sites | `BulkLoad.insertMany` | Do not fold singleton insertion; one batch thaws and publishes once. |
| Compose insertions and removals | one `Session.withSession` | Do not publish every intermediate mesh. Use coordinate-keyed removal after the first compaction. |
| Refine the whole mesh | `refine` | Do not manufacture a domain witness merely to reach the local API. |
| Refine a proved face section without changing protected faces | `refineWithinDomain` | Supply its exact permitted faces and interface edges; a guessed boundary is a typed refusal, not a hint. |
| Require construction-independent numbering | `canonicalize` at the observation boundary | Do not canonicalize every intermediate value; it is intentionally global work. |
## Foreign bindings
The `moonlight-triangulation-c` foreign-library component exposes opaque
immutable geometry meshes, batch insertion, finite-set operations, and dense
vertex and triangle projections through one versioned C ABI. Python,
TypeScript, and Rust bindings live beneath [`bindings/`](bindings/README.md)
and descend through that same ABI; none restates the triangulation engine or
introduces a second mesh representation.
`union a a` is `a`; commutativity and associativity hold after explicit
`canonicalize`; and a join adds no sites: the result carries
`|A| + |B| − |A ∩ B|` of them. Structural `Eq` remains exact resident equality
for caches and serialization rather than secretly rebuilding the mesh.
### Persistent publication schedules
The algebraic result is independent of the execution schedule. These are the
currently measured publication choices; `canonicalize` remains the explicit
global observation when construction-independent numbering is required.
| Operation context | Publication schedule |
| --- | --- |
| `siteRelation left right` | Index the smaller support in one transient exact open-addressed section, scan the other operand, then discard the index; no support maps or published cache state. |
| `difference mesh empty`, `symmetricDifference mesh empty`, `symmetricDifference empty mesh` | Return the surviving representative verbatim. |
| `difference left right` with `size right <= (size left - size right) / 128` | Remove the right support through one local copy-on-write session; return `left` verbatim when the supports are disjoint. Larger masks rebuild. |
| Geometry-only `intersection left right` | Return an existing operand for equality or subset, return the empty mesh for disjoint supports, and locally remove the smaller complement when it is at most `overlap / 128`. Other partial overlaps rebuild. |
| `intersectionWith` | Rebuild, because the annotation combiner may rewrite every surviving payload even when topology changes locally. |
| `symmetricDifference left right` with `small <= (large - small) / 128` | Toggle the smaller operand through one local session. Comparable operands and small-output, large-input cases rebuild. |
| Comparable `symmetricDifference left right` | Partition through one transient exact index and a matched bitset, then rebuild only the exclusive output section. |
| `BulkLoad.insert` / `BulkLoad.insertAt` | Dense publication below 10,000 resident sites; copy-on-write publication at 10,000 and above. A sequence still belongs in one `Session`. |
| `extendConstrainedWith` | Copy-on-write only for a base of at least 200,000 sites, at most 128 incoming sites, exactly one incoming segment, and a pre-thaw corridor with no resident intersection. Every unmeasured or resident-corridor case stays dense. |
| `refineWithinDomain` | Dense publication. The local transaction candidate preserved semantics but did not improve wall time, so it was removed. |
On the retained one-million-site / five-thousand-site witnesses, raw
`difference` takes 0.169 s. Near-full `intersection` takes 0.116 s and 0.100 s
across the cold and pre-forced contexts. A five-thousand-site
`symmetricDifference` result takes 0.121 s and allocates 244 MB. Exact
`siteRelation` takes 0.083 s and 0.115 s across those contexts. Singleton
insertion takes 5.326 ms. These are raw publication measurements;
`canonicalize` is measured separately when construction-independent numbering
is required.
```haskell
union :: JoinSemilattice annotation
=> Triangulation 'Unconstrained annotation () () ()
-> Triangulation 'Unconstrained annotation () () ()
-> Either BuildError (Triangulation 'Unconstrained annotation () () ())
unions :: JoinSemilattice annotation
=> [Triangulation 'Unconstrained annotation () () ()]
-> Either BuildError (Triangulation 'Unconstrained annotation () () ())
siteRelation
:: Triangulation leftMode leftAnnotation leftDirected leftUndirected leftFace
-> Triangulation rightMode rightAnnotation rightDirected rightUndirected rightFace
-> SiteRelation
intersection
:: Triangulation 'Unconstrained () () () ()
-> Triangulation 'Unconstrained () () () ()
-> Either BuildError (Triangulation 'Unconstrained () () () ())
intersectionWith
:: (leftAnnotation -> rightAnnotation -> annotation)
-> Triangulation 'Unconstrained leftAnnotation () () ()
-> Triangulation 'Unconstrained rightAnnotation () () ()
-> Either BuildError (Triangulation 'Unconstrained annotation () () ())
difference
:: Triangulation 'Unconstrained leftAnnotation () () ()
-> Triangulation 'Unconstrained rightAnnotation () () ()
-> Either BuildError (Triangulation 'Unconstrained leftAnnotation () () ())
symmetricDifference
:: Triangulation 'Unconstrained annotation () () ()
-> Triangulation 'Unconstrained annotation () () ()
-> Either BuildError (Triangulation 'Unconstrained annotation () () ())
```
## Constraints
* No `Semigroup`, no `Monoid`, no `<>`. A join the finite arena cannot
represent is a `Left`; a partial class instance would lie about totality.
* `union` and `unions` glue overlapping annotations through
`JoinSemilattice`. `intersectionWith` supplies the corresponding explicit
overlap combiner. `difference` preserves the left annotation and
`symmetricDifference` preserves the annotation of whichever exclusive site
survives. Plain `intersection` remains the geometry-only specialization.
Annotation-preserving restriction is derived rather than stored as a second
operator: `intersectionWith const left mask`.
* `refine` composes after an operation. It is not a mode, a flag, or a second
triangulation type.
* Mutation lives in `ST` behind `Moonlight.Triangulation.Internal.Mutable` and
never escapes.
* `serialize` is absent from the facade by design; import
`Moonlight.Triangulation.Serialization` to reach it.
* Coordinates are binary64 throughout the public surface. `vertexPoints` and
`innerFaceVertexTriples` project dense vectors directly from the authoritative
DCEL; callers do not reconstruct a sibling mesh DTO.
* One half-edge mesh underneath. `Cdt` is a mode index on it, not a second
structure.
## Use
A triangulation is a value of its site set: construction returns `Either` with
typed obstructions, `union` returns the same typed refusal when the finite arena
cannot represent its result, and mesh quality is a
composition over the result rather than a second operation. This program
compiles against the facade alone; `fromList` is `GHC.Exts`, so nothing beyond
`base` and `moonlight-triangulation` is in scope.
```haskell
module Main where
import GHC.Exts (fromList)
import Moonlight.Triangulation
main :: IO ()
main = do
let build :: [(Double, Double)] -> IO (DelaunayTriangulation ())
build coords =
either (fail . show) pure (delaunayGeometry (fromList [Point x y | (x, y) <- coords]))
a <- build [(x, y) | x <- [0 .. 9], y <- [0 .. 9]]
b <- build [(x + 6, y) | x <- [0 .. 9], y <- [0 .. 9]]
-- valid Delaunay representative of the union of sites
joined <- either (fail . show) pure (union a b)
-- meet and differences return typed obstructions, never partial results
met <- either (fail . show) pure (intersection a b)
-- mesh quality is a composition, not a second operation
parameters <- either (fail . show) pure (withMinimumAngle 14 defaultRefinementParameters)
result <- either (fail . show) pure (refine (const ()) parameters joined)
let healed = refinedTriangulation result
print (numVertices joined, numVertices met, numVertices healed, refinementComplete result)
```
Output: `(160,40,160,True)` — the union carries 160 sites, the meet 40, and
refinement inserts nothing because every triangle already clears the 14° bar;
`refinementComplete` states that sufficiency. Where operands leave concavities,
the same composition places Steiner sites exactly at the slivers it dissolves.
`delaunayGeometry` constructs the unit-annotated carrier directly. Payloads can
instead enter with their coordinates through `delaunay`, or independently
through `delaunayFromCoordinates`; annotation-preserving restriction is
available through `intersectionWith`, `difference`, and `symmetricDifference`. Constraints enter through
`constrainedDelaunay`; bulk incremental work names the machine-room module
directly (`Moonlight.Triangulation.BulkLoad` for `insertMany`,
`Moonlight.Triangulation.Session` for the owned editing transaction —
`withSession`: thaw once, insert and remove freely, publish once).
### Interior without hull-fill
A point set has no boundary, so `delaunay` necessarily meshes the convex hull:
concavities and holes are spanned by faces that belong to the hull, not to any
intended region. The region becomes real the moment its boundary is authored:
`constrainedDelaunay` takes the sites plus constraint segments as input-index
pairs, and interiority is then computable — `facesAtEvenBarrierDepth` runs a
0–1 BFS from the outer face and returns every face at even barrier depth,
which is exactly the outside (depth 0) plus anything nested behind a second
loop. The interior is the complement. `refine` consumes the same parity
through `refineExcludeOuterFaces`, so Steiner sites respect the boundary too.
```haskell
module Main where
import GHC.Exts (fromList)
import Moonlight.Triangulation
ring :: Double -> Int -> [(Double, Double)]
ring radius n =
[ (radius * cos t, radius * sin t)
| k <- [0 .. n - 1]
, let t = 2 * pi * fromIntegral k / fromIntegral n
]
main :: IO ()
main = do
let outer = ring 4 32
inner = ring 2 16
middle = ring 3 24
pts = fromList [Point x y | (x, y) <- outer <> inner <> middle]
loop base count = [(base + k, base + (k + 1) `mod` count) | k <- [0 .. count - 1]]
constraints = fromList (loop 0 32 <> loop 32 16)
annulus <- case constrainedDelaunay unitElementDefaults pts constraints of
Left err -> fail (show err)
Right result -> pure (buildTriangulation result)
let outside =
[ fromIntegral raw :: Int
| FaceId raw <- facesAtEvenBarrierDepth annulus (isConstraintEdge annulus)
]
interior =
[ FaceId (fromIntegral k)
| k <- [0 .. numFaces annulus - 1]
, not (k `elem` outside)
]
print (numFaces annulus, length interior)
```
Output: `(111,97)` — the annulus band is the 97 interior faces; the 14
excluded faces are the hole's hull-fill and the outer region. Rendering the
interior list draws the ring with a genuine void: no crop, no edge-length
heuristic, the engine's own verdict. `faceVertices` walks each interior face
for display, and a nested loop flips parity again, so islands inside holes
come back automatically.
### Deriving the boundary
When no boundary is known, the mesh itself carries one. Every Delaunay face
has a circumradius: faces inside a sampled region sit near the local pitch,
while faces spanning concavities and voids circumscribe them and blow up.
Keeping the faces below a threshold derived from the data — a multiple of the
median circumradius, so no authored constant — is the alpha-complex, and its
boundary falls out as the edges with exactly one kept side. Those edges are
already edges of the triangulation, so feeding them to `constrainedDelaunay`
as index pairs recovers without conflicts, and the parity machinery above
takes over from there.
```haskell
module Main where
import Data.List (sort, span)
import GHC.Exts (fromList)
import Moonlight.Triangulation
ring :: Double -> Int -> [(Double, Double)]
ring radius n =
[ (radius * cos t, radius * sin t)
| k <- [0 .. n - 1]
, let t = 2 * pi * fromIntegral k / fromIntegral n
]
main :: IO ()
main = do
let pts = ring 4 32 <> ring 2 16 <> ring 3 24
mesh <-
either
(fail . show)
pure
(delaunayGeometry (fromList [Point x y | (x, y) <- pts]))
let corners f =
[(x, y) | v <- faceVertices mesh f, let Point x y = vertexPoint mesh v]
circumradius (ax, ay) (bx, by) (cx, cy) =
let dab = sqrt ((bx - ax) ** 2 + (by - ay) ** 2)
dbc = sqrt ((cx - bx) ** 2 + (cy - by) ** 2)
dca = sqrt ((ax - cx) ** 2 + (ay - cy) ** 2)
area2 = abs ((bx - ax) * (cy - ay) - (cx - ax) * (by - ay))
in if area2 == 0 then 1 / 0 else dab * dbc * dca / (2 * area2)
radii =
[ (k, circumradius p q s)
| k <- [0 .. numFaces mesh - 1]
, [p, q, s] <- [corners (FaceId (fromIntegral k))]
]
sorted = sort (map snd radii)
medianR = case drop (length sorted `div` 2) sorted of
m : _ -> m
[] -> 1
kept = [FaceId (fromIntegral k) | (k, r) <- radii, r <= 1.35 * medianR]
norm (p, q) = if p <= q then (p, q) else (q, p)
edges f = case corners f of
[p, q, s] -> [norm (p, q), norm (q, s), norm (s, p)]
_ -> []
runs xs = case xs of
[] -> []
x : rest -> let (same, more) = span (== x) rest in (x, 1 + length same) : runs more
boundary = [e | (e, n) <- runs (sort (concatMap edges kept)), n == 1]
print (length radii, length kept, length boundary)
```
Output: `(110,96,48)` — 110 finite faces, 96 in the alpha-complex, and the
derived boundary is exactly the 48 ring edges: both loops recovered from the
point set alone. The derived boundary walks through data sites, so it is as
jagged as the sampling; an authored boundary stays smooth at any pitch and
wins where the generator is known.
## The growing city
A triangulation is a value so that a large one can be extended without being
rebuilt or locked: a city mesh a district is added to, and then another.
`union` rebuilds last. `planPair` classifies the pair first — an empty or
repeated operand returns the other verbatim, separable operands merge along
their seam, a subset inserts into its superset, skewed sizes insert the smaller
into the larger.
```haskell
insertionIsCheaper addition base = addition <= 64 || addition <= base `quot` 8
```
A district is small against a city, so it inserts through a local copy-on-write
transaction. Publication scales with the insertion work and pages dirtied
rather than copying and renumbering the whole city.
The mesh being read is never the mesh being written. A render or pathfinding
thread queries the published value while the next is built: no lock, no
defensive copy, no interval in which the world and the mesh disagree. Mutation
offers only the choice between stalling in place and answering from a clone
that does not yet know what was built, and both are visible from the frame.
The same fact makes a failed union a `Left` over an untouched operand,
reverting an edit a selection, and a preview a second value rather than a copy.
`Moonlight.Triangulation.Parallel` reads the operation the other way: regions
triangulated apart, folded up a balanced tournament, each node's sides
concurrent.
## Representation
Structure-of-arrays over paged, copy-on-write storage. Half-edge twins are index
complements (`e xor 1`), so traversal is arithmetic rather than indirection.
Local insertion transactions copy only the pages they touch; mutation is
confined to `ST` behind `Moonlight.Triangulation.Internal.Mutable` and never
escapes.
## Predicates
`Moonlight.Triangulation.Internal.Dyadic` carries the exact layer: mantissas are
decoded and aligned to a common exponent, and determinants are evaluated over
`Integer`. The floating approximation is trusted only inside Shewchuk's error
bounds — `(3 + 16u)u` for orientation, `(10 + 96u)u` for incircle — and falls
through to the exact evaluation otherwise.
## Public sublibraries
Depend on the sublibrary you actually use, not on the facade. The footprint is
the interface: each row states what the token costs you, and a build that pulls
more than the row says is a bug in this table or in the cabal.
| Sublibrary | Surface |
| --- | --- |
| `core` | Exact-arithmetic scalars (`Scalar`, `LineSideInfo`) over paged, copy-on-write storage (`Internal.Dyadic`, `.Paged`, `.BoxedPaged`, `.PageDirectory`, `.Growable`, `.PackedIndex`, `.FaceQueue`). Depends only on `base`/`containers`/`deepseq`/`vector` — no other sublibrary. |
| `dcel` | The finite half-edge mesh and its whole read surface: `Types`, `Math`, `Interop`, `Dcel`, `Payload`, `JoinSemilattice`, `Handles` with its iterator family, `PointLocation`, `Validation`, `FloodFillIterator`, `IntersectionIterator`. Adds `primitive` and `vector-algorithms` over `core`. |
| `build` | Everything that *constructs*: `BulkLoad`, `Session`, `Removal`, typed `SetAlgebra`, `Cdt` constraint recovery, and `Refinement`. Adds no external dependency over `core` and `dcel`. |
| `parallel` | Concurrent evaluation of the pure union plan. Adds `async` at this effect boundary rather than below it. |
| `serialize` | `Serialization` — the versioned binary envelope, and the only sublibrary that costs you `binary`, `bytestring`, and `transformers`. Deliberately absent from the facade. |
| `dual` | The Voronoi dual and what reads it: `Voronoi`, `Voronoi.Handles`, `Interpolation` (natural-neighbour), `HintGenerator` (Delaunay hierarchy hints). Sits over `core`, `dcel` and `build`. |
| facade (`moonlight-triangulation`) | `Moonlight.Triangulation` alone: the equational surface, one export list stating a theory. It re-exports selected names from `core`/`dcel`/`build`/`dual` — and *not* `serialize`. A caller who wants more than the theory names the machine-room module directly. |
## Surface
* `Moonlight.Triangulation` — the apex facade.
* `.BulkLoad` — circle-sweep construction and incremental insertion.
* `.Cdt` — constraint recovery by conflict strip, with requeue on re-intersection.
* `.PointLocation`, `.HintGenerator` — walk location and Delaunay hierarchy hints.
* `.IntersectionIterator`, `.FloodFillIterator` — ordered line traversal and barrier fill.
* `.Voronoi`, `.Interpolation` — dual cells and natural-neighbour interpolation.
* `.Refinement` — Ruppert-style angle and area refinement.
* `.Validation` — structural and Delaunay-property audits.
* `.Serialization` (sublibrary `serialize`) — versioned binary envelope.