moonlight-triangulation

Part of Moonlight, the sheaf-theoretic computation layer beneath Melusine 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

Choose the operation

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.

On the retained one-million-site / five-thousand-site witnesses, raw difference fell from 4.997 s to 0.169 s. Near-full intersection now takes 0.116 s from a cold persistent-index context and 0.100 s after that cache has already been forced, down from 3.475 s before indexed descent. A five-thousand-site symmetricDifference result fell from 4.022 s / 2.19 GB allocated to 0.121 s / 244 MB. Exact siteRelation now takes 0.083 s cold and 0.115 s after forcing the persistent cache, down from 3.417 s before indexed descent. Singleton insertion remains 5.326 ms from 48.037 ms. These are raw publication measurements; paying for canonicalize afterward is deliberately reported separately rather than smuggled into the local-operation claim.

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.

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

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.

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.

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.

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.