# 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 |
| --- | --- | --- |
| 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. |

`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 Double ())
      build coords = case delaunay unitElementDefaults (fromList [Point x y | (x, y) <- coords]) of
        Left err -> fail (show err)
        Right result -> pure (mapVertices (const ()) (buildTriangulation result))
  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.

Payloads annotate geometry through `mapVertices`; the input `Point`s arrive as
their own vertex payload. The example deliberately erases them with
`mapVertices (const ())`, while 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 <- case delaunay unitElementDefaults (fromList [Point x y | (x, y) <- pts]) of
    Left err -> fail (show err)
    Right result -> pure (mapVertices (const () :: Point -> ()) (buildTriangulation result))
  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.