-- SPDX-License-Identifier: MPL-2.0

{- | Deterministic Bowyer-Watson Delaunay triangulation (at most 256 unique
points), edge-flip constrained Delaunay triangulation, polygon ear clipping,
and bounded Voronoi cells. Triangles use original zero-based point indices
and counterclockwise orientation. Collinear clouds return no triangles.
Constraints must be noncrossing input-vertex segments with no interior input
vertices; they do not specify holes or discard the outer convex hull. Exact
Rational orientation/incircle signs avoid epsilon topology decisions.
-}
module HMath.Geometry.Triangulation (
    IndexedTriangle,
    delaunay,
    constrainedDelaunay,
    voronoi,
    earClip,
) where

import Control.Monad (foldM)
import Data.List (delete, nub, sort)
import HMath.Error (MathError (..))
import HMath.Geometry.Internal
import HMath.Geometry.Primitives (Intersection2 (..), segmentIntersection)

-- | Counterclockwise indices into the caller's input point list.
type IndexedTriangle = (Int, Int, Int)

triEdges :: IndexedTriangle -> [(Int, Int)]
triEdges :: IndexedTriangle -> [(Int, Int)]
triEdges (Int
a, Int
b, Int
c) = [(Int
a, Int
b), (Int
b, Int
c), (Int
c, Int
a)]
canonical :: (Int, Int) -> (Int, Int)
canonical :: (Int, Int) -> (Int, Int)
canonical (Int
a, Int
b) = if Int
a Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
b then (Int
a, Int
b) else (Int
b, Int
a)
orientTri :: [Point2] -> IndexedTriangle -> IndexedTriangle
orientTri :: [Point2] -> IndexedTriangle -> IndexedTriangle
orientTri [Point2]
ps (Int
a, Int
b, Int
c) = if Point2 -> Point2 -> Point2 -> Ordering
orient ([Point2]
ps [Point2] -> Int -> Point2
forall a. HasCallStack => [a] -> Int -> a
!! Int
a) ([Point2]
ps [Point2] -> Int -> Point2
forall a. HasCallStack => [a] -> Int -> a
!! Int
b) ([Point2]
ps [Point2] -> Int -> Point2
forall a. HasCallStack => [a] -> Int -> a
!! Int
c) Ordering -> Ordering -> Bool
forall a. Eq a => a -> a -> Bool
== Ordering
LT then (Int
a, Int
c, Int
b) else (Int
a, Int
b, Int
c)
insideCircle :: Point2 -> Point2 -> Point2 -> Point2 -> Ordering
insideCircle :: Point2 -> Point2 -> Point2 -> Point2 -> Ordering
insideCircle (P2 Double
ax Double
ay) (P2 Double
bx Double
by) (P2 Double
cx Double
cy) (P2 Double
px Double
py) =
    let a :: Rational
a = Double -> Rational
forall a. Real a => a -> Rational
toRational Double
ax Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
- Double -> Rational
forall a. Real a => a -> Rational
toRational Double
px
        b :: Rational
b = Double -> Rational
forall a. Real a => a -> Rational
toRational Double
ay Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
- Double -> Rational
forall a. Real a => a -> Rational
toRational Double
py
        c :: Rational
c = Double -> Rational
forall a. Real a => a -> Rational
toRational Double
bx Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
- Double -> Rational
forall a. Real a => a -> Rational
toRational Double
px
        d :: Rational
d = Double -> Rational
forall a. Real a => a -> Rational
toRational Double
by Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
- Double -> Rational
forall a. Real a => a -> Rational
toRational Double
py
        e :: Rational
e = Double -> Rational
forall a. Real a => a -> Rational
toRational Double
cx Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
- Double -> Rational
forall a. Real a => a -> Rational
toRational Double
px
        f :: Rational
f = Double -> Rational
forall a. Real a => a -> Rational
toRational Double
cy Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
- Double -> Rational
forall a. Real a => a -> Rational
toRational Double
py
        determinant :: Rational
determinant =
            (Rational
a Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
* Rational
a Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
+ Rational
b Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
* Rational
b) Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
* (Rational
c Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
* Rational
f Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
- Rational
d Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
* Rational
e)
                Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
- (Rational
c Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
* Rational
c Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
+ Rational
d Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
* Rational
d) Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
* (Rational
a Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
* Rational
f Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
- Rational
b Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
* Rational
e)
                Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
+ (Rational
e Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
* Rational
e Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
+ Rational
f Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
* Rational
f) Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
* (Rational
a Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
* Rational
d Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
- Rational
b Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
* Rational
c)
        orientationSign :: Rational
orientationSign = Point2 -> Point2 -> Point2 -> Rational
orientExact (Double -> Double -> Point2
P2 Double
ax Double
ay) (Double -> Double -> Point2
P2 Double
bx Double
by) (Double -> Double -> Point2
P2 Double
cx Double
cy)
     in Rational -> Rational -> Ordering
forall a. Ord a => a -> a -> Ordering
compare (if Rational
orientationSign Rational -> Rational -> Bool
forall a. Ord a => a -> a -> Bool
> Rational
0 then Rational
determinant else -Rational
determinant) Rational
0
validateCloud :: [Point2] -> Either MathError ()
validateCloud :: [Point2] -> Either MathError ()
validateCloud [Point2]
ps = do
    Bool -> String -> Either MathError ()
require ([Point2] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Point2]
ps Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
256) String
"triangulation reference budget is 256 points"
    (Point2 -> Either MathError ()) -> [Point2] -> Either MathError ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
(a -> m b) -> t a -> m ()
mapM_ Point2 -> Either MathError ()
check2 [Point2]
ps
    Bool -> String -> Either MathError ()
require ([Point2] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length ([Point2] -> [Point2]
forall a. Eq a => [a] -> [a]
nub [Point2]
ps) Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== [Point2] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Point2]
ps) String
"duplicate triangulation points"

{- | Incremental cavity replacement. Cocircular insertion uses a strict
in-circle predicate, giving a deterministic (not unique) valid diagonal.
-}
delaunay :: [Point2] -> Either MathError [IndexedTriangle]
delaunay :: [Point2] -> Either MathError [IndexedTriangle]
delaunay [Point2]
ps = do
    [Point2] -> Either MathError ()
validateCloud [Point2]
ps
    if [Point2] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Point2]
ps Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
3 Bool -> Bool -> Bool
|| (Point2 -> Bool) -> [Point2] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all (\Point2
p -> Point2 -> Point2 -> Point2 -> Ordering
orient (Point2 -> [Point2] -> Point2
forall a. a -> [a] -> a
firstOr (Double -> Double -> Point2
P2 Double
0 Double
0) [Point2]
ps) ([Point2]
ps [Point2] -> Int -> Point2
forall a. HasCallStack => [a] -> Int -> a
!! Int
1) Point2
p Ordering -> Ordering -> Bool
forall a. Eq a => a -> a -> Bool
== Ordering
EQ) (Int -> [Point2] -> [Point2]
forall a. Int -> [a] -> [a]
drop Int
2 [Point2]
ps)
        then [IndexedTriangle] -> Either MathError [IndexedTriangle]
forall a b. b -> Either a b
Right []
        else do
            let xs :: [Double]
xs = [Double
x | P2 Double
x Double
_ <- [Point2]
ps]
                ys :: [Double]
ys = [Double
y | P2 Double
_ Double
y <- [Point2]
ps]
                cx :: Double
cx = ([Double] -> Double
forall a. Ord a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
minimum [Double]
xs Double -> Double -> Double
forall a. Num a => a -> a -> a
+ [Double] -> Double
forall a. Ord a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum [Double]
xs) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2
                cy :: Double
cy = ([Double] -> Double
forall a. Ord a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
minimum [Double]
ys Double -> Double -> Double
forall a. Num a => a -> a -> a
+ [Double] -> Double
forall a. Ord a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum [Double]
ys) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2
                span' :: Double
span' = Double -> Double -> Double
forall a. Ord a => a -> a -> a
max ([Double] -> Double
forall a. Ord a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum [Double]
xs Double -> Double -> Double
forall a. Num a => a -> a -> a
- [Double] -> Double
forall a. Ord a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
minimum [Double]
xs) ([Double] -> Double
forall a. Ord a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum [Double]
ys Double -> Double -> Double
forall a. Num a => a -> a -> a
- [Double] -> Double
forall a. Ord a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
minimum [Double]
ys)
                s :: Double
s = Double -> Double -> Double
forall a. Ord a => a -> a -> a
max Double
span' Double
1e-100
                n :: Int
n = [Point2] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Point2]
ps
                allPoints :: [Point2]
allPoints = [Point2]
ps [Point2] -> [Point2] -> [Point2]
forall a. [a] -> [a] -> [a]
++ [Double -> Double -> Point2
P2 (Double
cx Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
32 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
s) (Double
cy Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
16 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
s), Double -> Double -> Point2
P2 (Double
cx Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
32 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
s) (Double
cy Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
16 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
s), Double -> Double -> Point2
P2 Double
cx (Double
cy Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
32 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
s)]
                insert :: [IndexedTriangle] -> Int -> [IndexedTriangle]
insert [IndexedTriangle]
active Int
i =
                    let p :: Point2
p = [Point2]
allPoints [Point2] -> Int -> Point2
forall a. HasCallStack => [a] -> Int -> a
!! Int
i
                        bad :: [IndexedTriangle]
bad =
                            [IndexedTriangle
t | t :: IndexedTriangle
t@(Int
a, Int
b, Int
c) <- [IndexedTriangle]
active, Point2 -> Point2 -> Point2 -> Point2 -> Ordering
insideCircle ([Point2]
allPoints [Point2] -> Int -> Point2
forall a. HasCallStack => [a] -> Int -> a
!! Int
a) ([Point2]
allPoints [Point2] -> Int -> Point2
forall a. HasCallStack => [a] -> Int -> a
!! Int
b) ([Point2]
allPoints [Point2] -> Int -> Point2
forall a. HasCallStack => [a] -> Int -> a
!! Int
c) Point2
p Ordering -> Ordering -> Bool
forall a. Eq a => a -> a -> Bool
== Ordering
GT]
                        allEdges :: [(Int, Int)]
allEdges = (IndexedTriangle -> [(Int, Int)])
-> [IndexedTriangle] -> [(Int, Int)]
forall (t :: * -> *) a b. Foldable t => (a -> [b]) -> t a -> [b]
concatMap IndexedTriangle -> [(Int, Int)]
triEdges [IndexedTriangle]
bad
                        boundary :: [(Int, Int)]
boundary = [(Int, Int)
e | (Int, Int)
e <- [(Int, Int)]
allEdges, [(Int, Int)] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length (((Int, Int) -> Bool) -> [(Int, Int)] -> [(Int, Int)]
forall a. (a -> Bool) -> [a] -> [a]
filter (((Int, Int) -> (Int, Int) -> Bool
forall a. Eq a => a -> a -> Bool
== (Int, Int) -> (Int, Int)
canonical (Int, Int)
e) ((Int, Int) -> Bool)
-> ((Int, Int) -> (Int, Int)) -> (Int, Int) -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Int, Int) -> (Int, Int)
canonical) [(Int, Int)]
allEdges) Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
1]
                        created :: [IndexedTriangle]
created =
                            [ [Point2] -> IndexedTriangle -> IndexedTriangle
orientTri [Point2]
allPoints (Int
a, Int
b, Int
i) | (Int
a, Int
b) <- [(Int, Int)]
boundary, Point2 -> Point2 -> Point2 -> Ordering
orient ([Point2]
allPoints [Point2] -> Int -> Point2
forall a. HasCallStack => [a] -> Int -> a
!! Int
a) ([Point2]
allPoints [Point2] -> Int -> Point2
forall a. HasCallStack => [a] -> Int -> a
!! Int
b) Point2
p Ordering -> Ordering -> Bool
forall a. Eq a => a -> a -> Bool
/= Ordering
EQ
                            ]
                     in (IndexedTriangle -> Bool) -> [IndexedTriangle] -> [IndexedTriangle]
forall a. (a -> Bool) -> [a] -> [a]
filter (IndexedTriangle -> [IndexedTriangle] -> Bool
forall (t :: * -> *) a. (Foldable t, Eq a) => a -> t a -> Bool
`notElem` [IndexedTriangle]
bad) [IndexedTriangle]
active [IndexedTriangle] -> [IndexedTriangle] -> [IndexedTriangle]
forall a. [a] -> [a] -> [a]
++ [IndexedTriangle]
created
                triangles :: [IndexedTriangle]
triangles = ([IndexedTriangle] -> Int -> [IndexedTriangle])
-> [IndexedTriangle] -> [Int] -> [IndexedTriangle]
forall b a. (b -> a -> b) -> b -> [a] -> b
forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl [IndexedTriangle] -> Int -> [IndexedTriangle]
insert [(Int
n, Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1, Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
2)] [Int
0 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]
            [IndexedTriangle] -> Either MathError [IndexedTriangle]
forall a b. b -> Either a b
Right ([IndexedTriangle] -> [IndexedTriangle]
forall a. Ord a => [a] -> [a]
sort [IndexedTriangle
t | t :: IndexedTriangle
t@(Int
a, Int
b, Int
c) <- [IndexedTriangle]
triangles, Int
a Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
n, Int
b Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
n, Int
c Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
n])

{- | Insert noncrossing constrained edges by convex quadrilateral flips, then
Lawson-flip every unconstrained illegal edge. Work is bounded explicitly.
-}
constrainedDelaunay :: [Point2] -> [(Int, Int)] -> Either MathError [IndexedTriangle]
constrainedDelaunay :: [Point2] -> [(Int, Int)] -> Either MathError [IndexedTriangle]
constrainedDelaunay [Point2]
ps [(Int, Int)]
constraints = do
    [Point2] -> Either MathError ()
validateCloud [Point2]
ps
    let n :: Int
n = [Point2] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Point2]
ps; cs :: [(Int, Int)]
cs = [(Int, Int)] -> [(Int, Int)]
forall a. Eq a => [a] -> [a]
nub (((Int, Int) -> (Int, Int)) -> [(Int, Int)] -> [(Int, Int)]
forall a b. (a -> b) -> [a] -> [b]
map (Int, Int) -> (Int, Int)
canonical [(Int, Int)]
constraints)
    Bool -> String -> Either MathError ()
require
        ([(Int, Int)] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [(Int, Int)]
cs Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
512 Bool -> Bool -> Bool
&& ((Int, Int) -> Bool) -> [(Int, Int)] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all (\(Int
a, Int
b) -> Int
a Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
0 Bool -> Bool -> Bool
&& Int
b Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
n Bool -> Bool -> Bool
&& Int
a Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
b) [(Int, Int)]
cs)
        String
"invalid constraint indices/budget"
    Bool -> String -> Either MathError ()
require
        ( [Bool] -> Bool
forall (t :: * -> *). Foldable t => t Bool -> Bool
and
            [Bool -> Bool
not (Point2 -> Point2 -> Point2 -> Bool
onSegment ([Point2]
ps [Point2] -> Int -> Point2
forall a. HasCallStack => [a] -> Int -> a
!! Int
a) ([Point2]
ps [Point2] -> Int -> Point2
forall a. HasCallStack => [a] -> Int -> a
!! Int
b) ([Point2]
ps [Point2] -> Int -> Point2
forall a. HasCallStack => [a] -> Int -> a
!! Int
i)) | (Int
a, Int
b) <- [(Int, Int)]
cs, Int
i <- [Int
0 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1], Int
i Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
a, Int
i Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
b]
        )
        String
"constraint contains an interior input point"
    checks <-
        (((Int, Int), (Int, Int)) -> Either MathError Intersection2)
-> [((Int, Int), (Int, Int))] -> Either MathError [Intersection2]
forall (t :: * -> *) (m :: * -> *) a b.
(Traversable t, Monad m) =>
(a -> m b) -> t a -> m (t b)
forall (m :: * -> *) a b. Monad m => (a -> m b) -> [a] -> m [b]
mapM
            (\((Int
a, Int
b), (Int
c, Int
d)) -> Point2
-> Point2 -> Point2 -> Point2 -> Either MathError Intersection2
segmentIntersection ([Point2]
ps [Point2] -> Int -> Point2
forall a. HasCallStack => [a] -> Int -> a
!! Int
a) ([Point2]
ps [Point2] -> Int -> Point2
forall a. HasCallStack => [a] -> Int -> a
!! Int
b) ([Point2]
ps [Point2] -> Int -> Point2
forall a. HasCallStack => [a] -> Int -> a
!! Int
c) ([Point2]
ps [Point2] -> Int -> Point2
forall a. HasCallStack => [a] -> Int -> a
!! Int
d))
            [((Int, Int)
x, (Int, Int)
y) | (Int, Int)
x <- [(Int, Int)]
cs, (Int, Int)
y <- [(Int, Int)]
cs, (Int, Int)
x (Int, Int) -> (Int, Int) -> Bool
forall a. Ord a => a -> a -> Bool
< (Int, Int)
y, Bool -> Bool
not ((Int, Int) -> (Int, Int) -> Bool
forall {a}. Eq a => (a, a) -> (a, a) -> Bool
shared (Int, Int)
x (Int, Int)
y)]
    require (all (== NoIntersection) checks) "constraints intersect"
    initial <- delaunay ps
    (inserted, _) <- foldM insert (initial, []) cs
    legalize inserted cs 0
  where
    shared :: (a, a) -> (a, a) -> Bool
shared (a
a, a
b) (a
c, a
d) = a
a a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
c Bool -> Bool -> Bool
|| a
a a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
d Bool -> Bool -> Bool
|| a
b a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
c Bool -> Bool -> Bool
|| a
b a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
d
    uniqueEdges :: t IndexedTriangle -> [(Int, Int)]
uniqueEdges t IndexedTriangle
ts = [(Int, Int)] -> [(Int, Int)]
forall a. Eq a => [a] -> [a]
nub (((Int, Int) -> (Int, Int)) -> [(Int, Int)] -> [(Int, Int)]
forall a b. (a -> b) -> [a] -> [b]
map (Int, Int) -> (Int, Int)
canonical ((IndexedTriangle -> [(Int, Int)])
-> t IndexedTriangle -> [(Int, Int)]
forall (t :: * -> *) a b. Foldable t => (a -> [b]) -> t a -> [b]
concatMap IndexedTriangle -> [(Int, Int)]
triEdges t IndexedTriangle
ts))
    adjacent :: (Int, Int) -> [IndexedTriangle] -> [IndexedTriangle]
adjacent (Int, Int)
edge [IndexedTriangle]
ts = [IndexedTriangle
t | IndexedTriangle
t <- [IndexedTriangle]
ts, (Int, Int)
edge (Int, Int) -> [(Int, Int)] -> Bool
forall a. Eq a => a -> [a] -> Bool
forall (t :: * -> *) a. (Foldable t, Eq a) => a -> t a -> Bool
`elem` ((Int, Int) -> (Int, Int)) -> [(Int, Int)] -> [(Int, Int)]
forall a b. (a -> b) -> [a] -> [b]
map (Int, Int) -> (Int, Int)
canonical (IndexedTriangle -> [(Int, Int)]
triEdges IndexedTriangle
t)]
    opposite :: (a, a) -> (a, a, a) -> a
opposite (a
a, a
b) (a
x, a
y, a
z) = a -> [a] -> a
forall a. a -> [a] -> a
firstOr a
a ((a -> Bool) -> [a] -> [a]
forall a. (a -> Bool) -> [a] -> [a]
filter (\a
p -> a
p a -> a -> Bool
forall a. Eq a => a -> a -> Bool
/= a
a Bool -> Bool -> Bool
&& a
p a -> a -> Bool
forall a. Eq a => a -> a -> Bool
/= a
b) [a
x, a
y, a
z])
    flipEdge :: (Int, Int) -> [IndexedTriangle] -> Maybe [IndexedTriangle]
flipEdge edge :: (Int, Int)
edge@(Int
a, Int
b) [IndexedTriangle]
ts = case (Int, Int) -> [IndexedTriangle] -> [IndexedTriangle]
adjacent (Int, Int)
edge [IndexedTriangle]
ts of
        [IndexedTriangle
u, IndexedTriangle
v] ->
            let c :: Int
c = (Int, Int) -> IndexedTriangle -> Int
forall {a}. Eq a => (a, a) -> (a, a, a) -> a
opposite (Int, Int)
edge IndexedTriangle
u; d :: Int
d = (Int, Int) -> IndexedTriangle -> Int
forall {a}. Eq a => (a, a) -> (a, a, a) -> a
opposite (Int, Int)
edge IndexedTriangle
v
             in if Point2 -> Point2 -> Point2 -> Ordering
orient ([Point2]
ps [Point2] -> Int -> Point2
forall a. HasCallStack => [a] -> Int -> a
!! Int
c) ([Point2]
ps [Point2] -> Int -> Point2
forall a. HasCallStack => [a] -> Int -> a
!! Int
d) ([Point2]
ps [Point2] -> Int -> Point2
forall a. HasCallStack => [a] -> Int -> a
!! Int
a) Ordering -> Ordering -> Bool
forall a. Eq a => a -> a -> Bool
== Point2 -> Point2 -> Point2 -> Ordering
orient ([Point2]
ps [Point2] -> Int -> Point2
forall a. HasCallStack => [a] -> Int -> a
!! Int
c) ([Point2]
ps [Point2] -> Int -> Point2
forall a. HasCallStack => [a] -> Int -> a
!! Int
d) ([Point2]
ps [Point2] -> Int -> Point2
forall a. HasCallStack => [a] -> Int -> a
!! Int
b)
                    Bool -> Bool -> Bool
|| Point2 -> Point2 -> Point2 -> Ordering
orient ([Point2]
ps [Point2] -> Int -> Point2
forall a. HasCallStack => [a] -> Int -> a
!! Int
c) ([Point2]
ps [Point2] -> Int -> Point2
forall a. HasCallStack => [a] -> Int -> a
!! Int
d) ([Point2]
ps [Point2] -> Int -> Point2
forall a. HasCallStack => [a] -> Int -> a
!! Int
a) Ordering -> Ordering -> Bool
forall a. Eq a => a -> a -> Bool
== Ordering
EQ
                    Bool -> Bool -> Bool
|| Point2 -> Point2 -> Point2 -> Ordering
orient ([Point2]
ps [Point2] -> Int -> Point2
forall a. HasCallStack => [a] -> Int -> a
!! Int
c) ([Point2]
ps [Point2] -> Int -> Point2
forall a. HasCallStack => [a] -> Int -> a
!! Int
d) ([Point2]
ps [Point2] -> Int -> Point2
forall a. HasCallStack => [a] -> Int -> a
!! Int
b) Ordering -> Ordering -> Bool
forall a. Eq a => a -> a -> Bool
== Ordering
EQ
                    then Maybe [IndexedTriangle]
forall a. Maybe a
Nothing
                    else [IndexedTriangle] -> Maybe [IndexedTriangle]
forall a. a -> Maybe a
Just ((IndexedTriangle -> Bool) -> [IndexedTriangle] -> [IndexedTriangle]
forall a. (a -> Bool) -> [a] -> [a]
filter (\IndexedTriangle
t -> IndexedTriangle
t IndexedTriangle -> IndexedTriangle -> Bool
forall a. Eq a => a -> a -> Bool
/= IndexedTriangle
u Bool -> Bool -> Bool
&& IndexedTriangle
t IndexedTriangle -> IndexedTriangle -> Bool
forall a. Eq a => a -> a -> Bool
/= IndexedTriangle
v) [IndexedTriangle]
ts [IndexedTriangle] -> [IndexedTriangle] -> [IndexedTriangle]
forall a. [a] -> [a] -> [a]
++ [[Point2] -> IndexedTriangle -> IndexedTriangle
orientTri [Point2]
ps (Int
c, Int
d, Int
a), [Point2] -> IndexedTriangle -> IndexedTriangle
orientTri [Point2]
ps (Int
d, Int
c, Int
b)])
        [IndexedTriangle]
_ -> Maybe [IndexedTriangle]
forall a. Maybe a
Nothing
    insert :: ([IndexedTriangle], [(Int, Int)])
-> (Int, Int) -> Either MathError ([IndexedTriangle], [(Int, Int)])
insert ([IndexedTriangle]
ts, [(Int, Int)]
locked) (Int, Int)
edge = do
        out <- [IndexedTriangle] -> Int -> Either MathError [IndexedTriangle]
seek [IndexedTriangle]
ts Int
0
        Right (out, edge : locked)
      where
        seek :: [IndexedTriangle] -> Int -> Either MathError [IndexedTriangle]
seek [IndexedTriangle]
current Int
budget
            | (Int, Int)
edge (Int, Int) -> [(Int, Int)] -> Bool
forall a. Eq a => a -> [a] -> Bool
forall (t :: * -> *) a. (Foldable t, Eq a) => a -> t a -> Bool
`elem` [IndexedTriangle] -> [(Int, Int)]
forall {t :: * -> *}.
Foldable t =>
t IndexedTriangle -> [(Int, Int)]
uniqueEdges [IndexedTriangle]
current = [IndexedTriangle] -> Either MathError [IndexedTriangle]
forall a b. b -> Either a b
Right [IndexedTriangle]
current
            | Int
budget Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
200000 = MathError -> Either MathError [IndexedTriangle]
forall a b. a -> Either a b
Left (String -> Int -> MathError
NonConvergence String
"constrained edge insertion" Int
budget)
            | Bool
otherwise =
                let (Int
a, Int
b) = (Int, Int)
edge
                    crossing :: [(Int, Int)]
crossing =
                        [ (Int, Int)
e
                        | e :: (Int, Int)
e@(Int
c, Int
d) <- [IndexedTriangle] -> [(Int, Int)]
forall {t :: * -> *}.
Foldable t =>
t IndexedTriangle -> [(Int, Int)]
uniqueEdges [IndexedTriangle]
current
                        , (Int, Int)
e (Int, Int) -> [(Int, Int)] -> Bool
forall (t :: * -> *) a. (Foldable t, Eq a) => a -> t a -> Bool
`notElem` [(Int, Int)]
locked
                        , Bool -> Bool
not ((Int, Int) -> (Int, Int) -> Bool
forall {a}. Eq a => (a, a) -> (a, a) -> Bool
shared (Int, Int)
edge (Int, Int)
e)
                        , Point2 -> Point2 -> Point2 -> Ordering
orient ([Point2]
ps [Point2] -> Int -> Point2
forall a. HasCallStack => [a] -> Int -> a
!! Int
a) ([Point2]
ps [Point2] -> Int -> Point2
forall a. HasCallStack => [a] -> Int -> a
!! Int
b) ([Point2]
ps [Point2] -> Int -> Point2
forall a. HasCallStack => [a] -> Int -> a
!! Int
c) Ordering -> Ordering -> Bool
forall a. Eq a => a -> a -> Bool
/= Point2 -> Point2 -> Point2 -> Ordering
orient ([Point2]
ps [Point2] -> Int -> Point2
forall a. HasCallStack => [a] -> Int -> a
!! Int
a) ([Point2]
ps [Point2] -> Int -> Point2
forall a. HasCallStack => [a] -> Int -> a
!! Int
b) ([Point2]
ps [Point2] -> Int -> Point2
forall a. HasCallStack => [a] -> Int -> a
!! Int
d)
                        , Point2 -> Point2 -> Point2 -> Ordering
orient ([Point2]
ps [Point2] -> Int -> Point2
forall a. HasCallStack => [a] -> Int -> a
!! Int
c) ([Point2]
ps [Point2] -> Int -> Point2
forall a. HasCallStack => [a] -> Int -> a
!! Int
d) ([Point2]
ps [Point2] -> Int -> Point2
forall a. HasCallStack => [a] -> Int -> a
!! Int
a) Ordering -> Ordering -> Bool
forall a. Eq a => a -> a -> Bool
/= Point2 -> Point2 -> Point2 -> Ordering
orient ([Point2]
ps [Point2] -> Int -> Point2
forall a. HasCallStack => [a] -> Int -> a
!! Int
c) ([Point2]
ps [Point2] -> Int -> Point2
forall a. HasCallStack => [a] -> Int -> a
!! Int
d) ([Point2]
ps [Point2] -> Int -> Point2
forall a. HasCallStack => [a] -> Int -> a
!! Int
b)
                        ]
                    candidates :: [[IndexedTriangle]]
candidates = [[IndexedTriangle]
out | (Int, Int)
e <- [(Int, Int)]
crossing, Just [IndexedTriangle]
out <- [(Int, Int) -> [IndexedTriangle] -> Maybe [IndexedTriangle]
flipEdge (Int, Int)
e [IndexedTriangle]
current]]
                 in case [[IndexedTriangle]]
candidates of
                        [] -> MathError -> Either MathError [IndexedTriangle]
forall a b. a -> Either a b
Left (String -> MathError
DomainError String
"constraint cannot be inserted into this triangulation")
                        [IndexedTriangle]
out : [[IndexedTriangle]]
_ -> [IndexedTriangle] -> Int -> Either MathError [IndexedTriangle]
seek [IndexedTriangle]
out (Int
budget Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1)
    legalize :: [IndexedTriangle]
-> t (Int, Int) -> Int -> Either MathError [IndexedTriangle]
legalize [IndexedTriangle]
ts t (Int, Int)
locked Int
budget
        | Int
budget Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
200000 = MathError -> Either MathError [IndexedTriangle]
forall a b. a -> Either a b
Left (String -> Int -> MathError
NonConvergence String
"constrained Delaunay legalization" Int
budget)
        | Bool
otherwise =
            let illegal :: [[IndexedTriangle]]
illegal =
                    [ [IndexedTriangle]
out
                    | (Int, Int)
edge <- [IndexedTriangle] -> [(Int, Int)]
forall {t :: * -> *}.
Foldable t =>
t IndexedTriangle -> [(Int, Int)]
uniqueEdges [IndexedTriangle]
ts
                    , (Int, Int)
edge (Int, Int) -> t (Int, Int) -> Bool
forall (t :: * -> *) a. (Foldable t, Eq a) => a -> t a -> Bool
`notElem` t (Int, Int)
locked
                    , [(Int
a, Int
b, Int
c), IndexedTriangle
v] <- [(Int, Int) -> [IndexedTriangle] -> [IndexedTriangle]
adjacent (Int, Int)
edge [IndexedTriangle]
ts]
                    , let d :: Int
d = (Int, Int) -> IndexedTriangle -> Int
forall {a}. Eq a => (a, a) -> (a, a, a) -> a
opposite (Int, Int)
edge IndexedTriangle
v
                    , Point2 -> Point2 -> Point2 -> Point2 -> Ordering
insideCircle ([Point2]
ps [Point2] -> Int -> Point2
forall a. HasCallStack => [a] -> Int -> a
!! Int
a) ([Point2]
ps [Point2] -> Int -> Point2
forall a. HasCallStack => [a] -> Int -> a
!! Int
b) ([Point2]
ps [Point2] -> Int -> Point2
forall a. HasCallStack => [a] -> Int -> a
!! Int
c) ([Point2]
ps [Point2] -> Int -> Point2
forall a. HasCallStack => [a] -> Int -> a
!! Int
d) Ordering -> Ordering -> Bool
forall a. Eq a => a -> a -> Bool
== Ordering
GT
                    , Just [IndexedTriangle]
out <- [(Int, Int) -> [IndexedTriangle] -> Maybe [IndexedTriangle]
flipEdge (Int, Int)
edge [IndexedTriangle]
ts]
                    ]
             in case [[IndexedTriangle]]
illegal of [] -> [IndexedTriangle] -> Either MathError [IndexedTriangle]
forall a b. b -> Either a b
Right [IndexedTriangle]
ts; [IndexedTriangle]
out : [[IndexedTriangle]]
_ -> [IndexedTriangle]
-> t (Int, Int) -> Int -> Either MathError [IndexedTriangle]
legalize [IndexedTriangle]
out t (Int, Int)
locked (Int
budget Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1)

{- | Voronoi cells clipped to an ordered nonzero rectangle, one polygon ring per
site (empty only if the site lies outside the clipping rectangle). Direct
pairwise half-plane clipping is O(n^3), limited to 128 sites. Bisector
classification and intersections use exact Rational arithmetic; final vertices
are rounded once to Double.
-}
voronoi :: Rectangle -> [Point2] -> Either MathError [[Point2]]
voronoi :: Rectangle -> [Point2] -> Either MathError [[Point2]]
voronoi (Rectangle lo :: Point2
lo@(P2 Double
x0 Double
y0) hi :: Point2
hi@(P2 Double
x1 Double
y1)) [Point2]
sites = do
    [Point2] -> Either MathError ()
validateCloud [Point2]
sites
    (Point2 -> Either MathError ()) -> [Point2] -> Either MathError ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
(a -> m b) -> t a -> m ()
mapM_ Point2 -> Either MathError ()
check2 [Point2
lo, Point2
hi]
    Bool -> String -> Either MathError ()
require ([Point2] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Point2]
sites Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
128 Bool -> Bool -> Bool
&& Double
x0 Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
x1 Bool -> Bool -> Bool
&& Double
y0 Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
y1) String
"invalid Voronoi rectangle/site budget"
    (Point2 -> Either MathError [Point2])
-> [Point2] -> Either MathError [[Point2]]
forall (t :: * -> *) (m :: * -> *) a b.
(Traversable t, Monad m) =>
(a -> m b) -> t a -> m (t b)
forall (m :: * -> *) a b. Monad m => (a -> m b) -> [a] -> m [b]
mapM Point2 -> Either MathError [Point2]
cell [Point2]
sites
  where
    rational :: Point2 -> (Rational, Rational)
rational (P2 Double
x Double
y) = (Double -> Rational
forall a. Real a => a -> Rational
toRational Double
x, Double -> Rational
forall a. Real a => a -> Rational
toRational Double
y)
    initial :: [(Rational, Rational)]
initial = (Point2 -> (Rational, Rational))
-> [Point2] -> [(Rational, Rational)]
forall a b. (a -> b) -> [a] -> [b]
map Point2 -> (Rational, Rational)
rational [Double -> Double -> Point2
P2 Double
x0 Double
y0, Double -> Double -> Point2
P2 Double
x1 Double
y0, Double -> Double -> Point2
P2 Double
x1 Double
y1, Double -> Double -> Point2
P2 Double
x0 Double
y1]
    cell :: Point2 -> Either MathError [Point2]
cell Point2
site =
        ((Rational, Rational) -> Either MathError Point2)
-> [(Rational, Rational)] -> Either MathError [Point2]
forall (t :: * -> *) (m :: * -> *) a b.
(Traversable t, Monad m) =>
(a -> m b) -> t a -> m (t b)
forall (m :: * -> *) a b. Monad m => (a -> m b) -> [a] -> m [b]
mapM
            (\(Rational
x, Rational
y) -> Double -> Double -> Either MathError Point2
point2 (Rational -> Double
forall a. Fractional a => Rational -> a
fromRational Rational
x) (Rational -> Double
forall a. Fractional a => Rational -> a
fromRational Rational
y))
            (([(Rational, Rational)]
 -> (Rational, Rational) -> [(Rational, Rational)])
-> [(Rational, Rational)]
-> [(Rational, Rational)]
-> [(Rational, Rational)]
forall b a. (b -> a -> b) -> b -> [a] -> b
forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl ((Rational, Rational)
-> [(Rational, Rational)]
-> (Rational, Rational)
-> [(Rational, Rational)]
forall {a}.
(Ord a, Fractional a) =>
(a, a) -> [(a, a)] -> (a, a) -> [(a, a)]
clip (Point2 -> (Rational, Rational)
rational Point2
site)) [(Rational, Rational)]
initial ((Point2 -> (Rational, Rational))
-> [Point2] -> [(Rational, Rational)]
forall a b. (a -> b) -> [a] -> [b]
map Point2 -> (Rational, Rational)
rational ((Point2 -> Bool) -> [Point2] -> [Point2]
forall a. (a -> Bool) -> [a] -> [a]
filter (Point2 -> Point2 -> Bool
forall a. Eq a => a -> a -> Bool
/= Point2
site) [Point2]
sites)))
    clip :: (a, a) -> [(a, a)] -> (a, a) -> [(a, a)]
clip (a, a)
_ [] (a, a)
_ = []
    clip (a
sx, a
sy) [(a, a)]
input (a
ox, a
oy) =
        -- Shift the bisector to the first site. Exact rational accumulation
        -- prevents subtracting two huge squared world-coordinate norms and
        -- rounds intersections only once, when the final cell is emitted.
        let nx :: a
nx = a
ox a -> a -> a
forall a. Num a => a -> a -> a
- a
sx
            ny :: a
ny = a
oy a -> a -> a
forall a. Num a => a -> a -> a
- a
sy
            offset :: a
offset = (a
nx a -> a -> a
forall a. Num a => a -> a -> a
* a
nx a -> a -> a
forall a. Num a => a -> a -> a
+ a
ny a -> a -> a
forall a. Num a => a -> a -> a
* a
ny) a -> a -> a
forall a. Fractional a => a -> a -> a
/ a
2
            value :: (a, a) -> a
value (a
px, a
py) = (a
px a -> a -> a
forall a. Num a => a -> a -> a
- a
sx) a -> a -> a
forall a. Num a => a -> a -> a
* a
nx a -> a -> a
forall a. Num a => a -> a -> a
+ (a
py a -> a -> a
forall a. Num a => a -> a -> a
- a
sy) a -> a -> a
forall a. Num a => a -> a -> a
* a
ny a -> a -> a
forall a. Num a => a -> a -> a
- a
offset
            crossing :: (a, a) -> (a, a) -> (a, a)
crossing a :: (a, a)
a@(a
ax, a
ay) b :: (a, a)
b@(a
bx, a
by) =
                let t :: a
t = (a, a) -> a
value (a, a)
a a -> a -> a
forall a. Fractional a => a -> a -> a
/ ((a, a) -> a
value (a, a)
a a -> a -> a
forall a. Num a => a -> a -> a
- (a, a) -> a
value (a, a)
b) in (a
ax a -> a -> a
forall a. Num a => a -> a -> a
+ a
t a -> a -> a
forall a. Num a => a -> a -> a
* (a
bx a -> a -> a
forall a. Num a => a -> a -> a
- a
ax), a
ay a -> a -> a
forall a. Num a => a -> a -> a
+ a
t a -> a -> a
forall a. Num a => a -> a -> a
* (a
by a -> a -> a
forall a. Num a => a -> a -> a
- a
ay))
         in [[(a, a)]] -> [(a, a)]
forall (t :: * -> *) a. Foldable t => t [a] -> [a]
concat
                [ case ((a, a) -> a
value (a, a)
a a -> a -> Bool
forall a. Ord a => a -> a -> Bool
<= a
0, (a, a) -> a
value (a, a)
b a -> a -> Bool
forall a. Ord a => a -> a -> Bool
<= a
0) of
                    (Bool
True, Bool
True) -> [(a, a)
b]
                    (Bool
True, Bool
False) -> [(a, a) -> (a, a) -> (a, a)
crossing (a, a)
a (a, a)
b]
                    (Bool
False, Bool
True) -> [(a, a) -> (a, a) -> (a, a)
crossing (a, a)
a (a, a)
b, (a, a)
b]
                    (Bool, Bool)
_ -> []
                | ((a, a)
a, (a, a)
b) <- [(a, a)] -> [((a, a), (a, a))]
forall a. [a] -> [(a, a)]
edges [(a, a)]
input
                ]

{- | Ear clipping of a simple polygon. Collinear boundary vertices are removed
from the active ring; returned triangles retain original indices.
-}
earClip :: Polygon -> Either MathError [IndexedTriangle]
earClip :: Polygon -> Either MathError [IndexedTriangle]
earClip Polygon
poly = do
    points <- Polygon -> Either MathError [Point2]
checkPolygon Polygon
poly
    let indices = if [Point2] -> Double
signedArea [Point2]
points Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
0 then [Int] -> [Int]
forall a. [a] -> [a]
reverse [Int
0 .. [Point2] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Point2]
points Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1] else [Int
0 .. [Point2] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Point2]
points Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]
    clip points indices []
  where
    clip :: [Point2]
-> [Int] -> [IndexedTriangle] -> Either MathError [IndexedTriangle]
clip [Point2]
_ [] [IndexedTriangle]
acc = [IndexedTriangle] -> Either MathError [IndexedTriangle]
forall a b. b -> Either a b
Right ([IndexedTriangle] -> [IndexedTriangle]
forall a. [a] -> [a]
reverse [IndexedTriangle]
acc)
    clip [Point2]
_ [Int
_] [IndexedTriangle]
acc = [IndexedTriangle] -> Either MathError [IndexedTriangle]
forall a b. b -> Either a b
Right ([IndexedTriangle] -> [IndexedTriangle]
forall a. [a] -> [a]
reverse [IndexedTriangle]
acc)
    clip [Point2]
_ [Int
_, Int
_] [IndexedTriangle]
acc = [IndexedTriangle] -> Either MathError [IndexedTriangle]
forall a b. b -> Either a b
Right ([IndexedTriangle] -> [IndexedTriangle]
forall a. [a] -> [a]
reverse [IndexedTriangle]
acc)
    clip [Point2]
_ [Int
a, Int
b, Int
c] [IndexedTriangle]
acc = [IndexedTriangle] -> Either MathError [IndexedTriangle]
forall a b. b -> Either a b
Right ([IndexedTriangle] -> [IndexedTriangle]
forall a. [a] -> [a]
reverse ((Int
a, Int
b, Int
c) IndexedTriangle -> [IndexedTriangle] -> [IndexedTriangle]
forall a. a -> [a] -> [a]
: [IndexedTriangle]
acc))
    clip [Point2]
ps [Int]
ring [IndexedTriangle]
acc =
        let triples :: [IndexedTriangle]
triples = [Int] -> [Int] -> [Int] -> [IndexedTriangle]
forall a b c. [a] -> [b] -> [c] -> [(a, b, c)]
zip3 ([Int] -> Int
forall a. HasCallStack => [a] -> a
last [Int]
ring Int -> [Int] -> [Int]
forall a. a -> [a] -> [a]
: [Int] -> [Int]
forall a. HasCallStack => [a] -> [a]
init [Int]
ring) [Int]
ring (Int -> [Int] -> [Int]
forall a. Int -> [a] -> [a]
drop Int
1 [Int]
ring [Int] -> [Int] -> [Int]
forall a. [a] -> [a] -> [a]
++ [Int -> [Int] -> Int
forall a. a -> [a] -> a
firstOr Int
0 [Int]
ring])
            collinear :: [Int]
collinear = [Int
b | (Int
a, Int
b, Int
c) <- [IndexedTriangle]
triples, Point2 -> Point2 -> Point2 -> Ordering
orient ([Point2]
ps [Point2] -> Int -> Point2
forall a. HasCallStack => [a] -> Int -> a
!! Int
a) ([Point2]
ps [Point2] -> Int -> Point2
forall a. HasCallStack => [a] -> Int -> a
!! Int
b) ([Point2]
ps [Point2] -> Int -> Point2
forall a. HasCallStack => [a] -> Int -> a
!! Int
c) Ordering -> Ordering -> Bool
forall a. Eq a => a -> a -> Bool
== Ordering
EQ]
            inside :: Point2 -> Point2 -> Point2 -> Point2 -> Bool
inside Point2
a Point2
b Point2
c Point2
p = (Ordering -> Bool) -> [Ordering] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all (Ordering -> Ordering -> Bool
forall a. Eq a => a -> a -> Bool
/= Ordering
LT) [Point2 -> Point2 -> Point2 -> Ordering
orient Point2
a Point2
b Point2
p, Point2 -> Point2 -> Point2 -> Ordering
orient Point2
b Point2
c Point2
p, Point2 -> Point2 -> Point2 -> Ordering
orient Point2
c Point2
a Point2
p]
            ears :: [IndexedTriangle]
ears =
                [ (Int
a, Int
b, Int
c)
                | (Int
a, Int
b, Int
c) <- [IndexedTriangle]
triples
                , Point2 -> Point2 -> Point2 -> Ordering
orient ([Point2]
ps [Point2] -> Int -> Point2
forall a. HasCallStack => [a] -> Int -> a
!! Int
a) ([Point2]
ps [Point2] -> Int -> Point2
forall a. HasCallStack => [a] -> Int -> a
!! Int
b) ([Point2]
ps [Point2] -> Int -> Point2
forall a. HasCallStack => [a] -> Int -> a
!! Int
c) Ordering -> Ordering -> Bool
forall a. Eq a => a -> a -> Bool
== Ordering
GT
                , (Int -> Bool) -> [Int] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all (\Int
i -> Int
i Int -> [Int] -> Bool
forall a. Eq a => a -> [a] -> Bool
forall (t :: * -> *) a. (Foldable t, Eq a) => a -> t a -> Bool
`elem` [Int
a, Int
b, Int
c] Bool -> Bool -> Bool
|| Bool -> Bool
not (Point2 -> Point2 -> Point2 -> Point2 -> Bool
inside ([Point2]
ps [Point2] -> Int -> Point2
forall a. HasCallStack => [a] -> Int -> a
!! Int
a) ([Point2]
ps [Point2] -> Int -> Point2
forall a. HasCallStack => [a] -> Int -> a
!! Int
b) ([Point2]
ps [Point2] -> Int -> Point2
forall a. HasCallStack => [a] -> Int -> a
!! Int
c) ([Point2]
ps [Point2] -> Int -> Point2
forall a. HasCallStack => [a] -> Int -> a
!! Int
i))) [Int]
ring
                ]
         in case [Int]
collinear of
                Int
b : [Int]
_ -> [Point2]
-> [Int] -> [IndexedTriangle] -> Either MathError [IndexedTriangle]
clip [Point2]
ps (Int -> [Int] -> [Int]
forall a. Eq a => a -> [a] -> [a]
delete Int
b [Int]
ring) [IndexedTriangle]
acc
                [] -> case [IndexedTriangle]
ears of
                    [] -> MathError -> Either MathError [IndexedTriangle]
forall a b. a -> Either a b
Left (String -> Int -> MathError
NonConvergence String
"ear clipping" ([Int] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Int]
ring))
                    t :: IndexedTriangle
t@(Int
_, Int
b, Int
_) : [IndexedTriangle]
_ -> [Point2]
-> [Int] -> [IndexedTriangle] -> Either MathError [IndexedTriangle]
clip [Point2]
ps (Int -> [Int] -> [Int]
forall a. Eq a => a -> [a] -> [a]
delete Int
b [Int]
ring) (IndexedTriangle
t IndexedTriangle -> [IndexedTriangle] -> [IndexedTriangle]
forall a. a -> [a] -> [a]
: [IndexedTriangle]
acc)