pure-cdt (empty) → 0.2.0.0

raw patch · 40 files changed

+3591/−0 lines, 40 filesdep +JuicyPixelsdep +Rasterificdep +basesetup-changed

Dependencies added: JuicyPixels, Rasterific, base, deepseq, hashable, parallel, pure-cdt, random, tasty, tasty-bench, tasty-hunit, tasty-quickcheck, unordered-containers, vector

Files

• CHANGELOG.md +38/−0
• LICENSE +26/−0
• README.md +131/−0
• Setup.hs +2/−0
• app/Drawer.hs +78/−0
• app/Examples/Circle.hs +37/−0
• app/Examples/Common.hs +16/−0
• app/Examples/Gear.hs +79/−0
• app/Examples/Simple.hs +13/−0
• app/Examples/Solenoid.hs +63/−0
• app/Main.hs +18/−0
• bench/Main.hs +39/−0
• bench/Scale.hs +97/−0
• pure-cdt.cabal +162/−0
• src/Triangulation.hs +34/−0
• src/Triangulation/Check.hs +43/−0
• src/Triangulation/Constrained.hs +92/−0
• src/Triangulation/Flip.hs +77/−0
• src/Triangulation/Generator.hs +53/−0
• src/Triangulation/Geometry.hs +16/−0
• src/Triangulation/Geometry/Edge.hs +66/−0
• src/Triangulation/Geometry/Exact.hs +87/−0
• src/Triangulation/Geometry/Point.hs +154/−0
• src/Triangulation/Geometry/Polygon.hs +182/−0
• src/Triangulation/Geometry/Ring.hs +56/−0
• src/Triangulation/Geometry/Triangle.hs +150/−0
• src/Triangulation/Leaf.hs +31/−0
• src/Triangulation/Merge.hs +81/−0
• src/Triangulation/Mesh.hs +74/−0
• src/Triangulation/Parallel.hs +138/−0
• src/Triangulation/Refine.hs +310/−0
• src/Triangulation/Repair.hs +94/−0
• src/Triangulation/Store.hs +124/−0
• src/Triangulation/Types.hs +40/−0
• test/Main.hs +20/−0
• test/Test/Geometry.hs +218/−0
• test/Test/Mesh.hs +67/−0
• test/Test/Refine.hs +221/−0
• test/Test/Ring.hs +50/−0
• test/Test/Triangulation.hs +314/−0

@@ -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`.

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

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

@@ -0,0 +1,2 @@+import Distribution.Simple+main = defaultMain

@@ -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))

@@ -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

@@ -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

@@ -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

@@ -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) [] []

@@ -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)

@@ -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 ()

@@ -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))

@@ -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}

@@ -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

@@ -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)

@@ -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]

@@ -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

@@ -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

@@ -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)

@@ -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

@@ -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))

@@ -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

@@ -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)))

@@ -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)

@@ -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)

@@ -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))

@@ -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]

@@ -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]

@@ -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)

@@ -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

@@ -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]

@@ -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)

@@ -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

@@ -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)

@@ -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+      ]

@@ -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+    ]

@@ -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)

@@ -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

@@ -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, [])+    ]

@@ -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]+  ]