-- | Port of Delaunator, a very fast Delaunay triangulation algorithm.
-- * Original Javascript source: https://github.com/mapbox/delaunator
-- * Rust port used as second reference: https://github.com/mourner/delaunator-rs
--
-- This module aimed to be as dumb and close to the Rust source implementation as
-- possible. Refactor as long as the tests pass.
--
-- Users of this internal module are responsible for implementing a nice API.
module Geometry.Algorithms.Delaunay.Internal.Delaunator.Raw where



import           Control.DeepSeq
import           Control.Monad
import           Control.Monad.ST
import           Data.Function
import           Data.Ord
import           Data.STRef
import           Data.Vector                (Vector, (!))
import qualified Data.Vector                as V
import qualified Data.Vector.Algorithms.Tim as VM
import           Data.Vector.Mutable        (STVector)
import qualified Data.Vector.Mutable        as VM
import           GHC.Stack                  (HasCallStack)
import           Geometry.Core



-- ^ Squared distance between two points
distSquare :: Vec2 -> Vec2 -> Double
distSquare :: Vec2 -> Vec2 -> Double
distSquare Vec2
x Vec2
y = Vec2 -> Double
normSquare (Vec2
x Vec2 -> Vec2 -> Vec2
forall v. VectorSpace v => v -> v -> v
-. Vec2
y)

-- | Orientation in screen coordinates (y downward).
data Orientation = Clockwise | Counterclockwise | Degenerate
    deriving (Orientation -> Orientation -> Bool
(Orientation -> Orientation -> Bool)
-> (Orientation -> Orientation -> Bool) -> Eq Orientation
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: Orientation -> Orientation -> Bool
== :: Orientation -> Orientation -> Bool
$c/= :: Orientation -> Orientation -> Bool
/= :: Orientation -> Orientation -> Bool
Eq, Eq Orientation
Eq Orientation
-> (Orientation -> Orientation -> Ordering)
-> (Orientation -> Orientation -> Bool)
-> (Orientation -> Orientation -> Bool)
-> (Orientation -> Orientation -> Bool)
-> (Orientation -> Orientation -> Bool)
-> (Orientation -> Orientation -> Orientation)
-> (Orientation -> Orientation -> Orientation)
-> Ord Orientation
Orientation -> Orientation -> Bool
Orientation -> Orientation -> Ordering
Orientation -> Orientation -> Orientation
forall a.
Eq a
-> (a -> a -> Ordering)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> a)
-> (a -> a -> a)
-> Ord a
$ccompare :: Orientation -> Orientation -> Ordering
compare :: Orientation -> Orientation -> Ordering
$c< :: Orientation -> Orientation -> Bool
< :: Orientation -> Orientation -> Bool
$c<= :: Orientation -> Orientation -> Bool
<= :: Orientation -> Orientation -> Bool
$c> :: Orientation -> Orientation -> Bool
> :: Orientation -> Orientation -> Bool
$c>= :: Orientation -> Orientation -> Bool
>= :: Orientation -> Orientation -> Bool
$cmax :: Orientation -> Orientation -> Orientation
max :: Orientation -> Orientation -> Orientation
$cmin :: Orientation -> Orientation -> Orientation
min :: Orientation -> Orientation -> Orientation
Ord, Int -> Orientation -> ShowS
[Orientation] -> ShowS
Orientation -> String
(Int -> Orientation -> ShowS)
-> (Orientation -> String)
-> ([Orientation] -> ShowS)
-> Show Orientation
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> Orientation -> ShowS
showsPrec :: Int -> Orientation -> ShowS
$cshow :: Orientation -> String
show :: Orientation -> String
$cshowList :: [Orientation] -> ShowS
showList :: [Orientation] -> ShowS
Show)

-- | Orientation of a triangle in screen coordinates (y downward).
--
-- >>> (a,b,c) = (Vec2 0 0, Vec2 100 0, Vec2 0 100)
-- >>> orientation a b c
-- Clockwise
-- >>> orientation b a c
-- Counterclockwise
orientation :: Vec2 -> Vec2 -> Vec2 -> Orientation
orientation :: Vec2 -> Vec2 -> Vec2 -> Orientation
orientation Vec2
x Vec2
q Vec2
r =
    let lineQ :: Vec2
lineQ = Vec2
q Vec2 -> Vec2 -> Vec2
forall v. VectorSpace v => v -> v -> v
-. Vec2
x
        lineR :: Vec2
lineR = Vec2
r Vec2 -> Vec2 -> Vec2
forall v. VectorSpace v => v -> v -> v
-. Vec2
x
    in case Double -> Double -> Ordering
forall a. Ord a => a -> a -> Ordering
compare (Vec2 -> Vec2 -> Double
cross Vec2
lineQ Vec2
lineR) Double
0 of
        Ordering
GT -> Orientation
Clockwise
        Ordering
EQ -> Orientation
Degenerate
        Ordering
LT -> Orientation
Counterclockwise

-- | Offset of the circumcenter of the triangle (a,b,c) from a.
--
-- >>> (a,b,c) = (Vec2 0 0, Vec2 100 0, Vec2 0 100)
-- >>> circumdelta a b c
-- Vec2 50.0 50.0
circumdelta :: Vec2 -> Vec2 -> Vec2 -> Vec2
circumdelta :: Vec2 -> Vec2 -> Vec2 -> Vec2
circumdelta Vec2
a Vec2
b Vec2
c =
    let Vec2 Double
dx Double
dy = Vec2
b Vec2 -> Vec2 -> Vec2
forall v. VectorSpace v => v -> v -> v
-. Vec2
a
        Vec2 Double
ex Double
ey = Vec2
c Vec2 -> Vec2 -> Vec2
forall v. VectorSpace v => v -> v -> v
-. Vec2
a

        bl :: Double
bl = Double
dxDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
dx Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
dyDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
dy;
        cl :: Double
cl = Double
exDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
ex Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
eyDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
ey;
        d :: Double
d = Double
0.5 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
dxDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
ey Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
dyDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
ex);

        x :: Double
x = (Double
eyDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
bl Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
dyDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
cl)Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
d;
        y :: Double
y = (Double
dxDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
cl Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
exDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
bl)Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
d;
    in Double -> Double -> Vec2
Vec2 Double
x Double
y

-- | Square of the circumradius’ norm of the triangle (a,b,c).
--
-- >>> (a,b,c) = (Vec2 0 0, Vec2 100 0, Vec2 0 100)
-- >>> circumradiusSquare a b c
-- 5000.0
circumradiusSquare :: Vec2 -> Vec2 -> Vec2 -> Double
circumradiusSquare :: Vec2 -> Vec2 -> Vec2 -> Double
circumradiusSquare Vec2
a Vec2
b Vec2
c = Vec2 -> Double
normSquare (Vec2 -> Vec2 -> Vec2 -> Vec2
circumdelta Vec2
a Vec2
b Vec2
c)

-- | Coordinate of the circumcenter of the triangle (a,b,c).
--
-- >>> (a,b,c) = (Vec2 0 0, Vec2 100 0, Vec2 0 100)
-- >>> circumcenter a b c
-- Vec2 50.0 50.0
circumcenter :: Vec2 -> Vec2 -> Vec2 -> Vec2
circumcenter :: Vec2 -> Vec2 -> Vec2 -> Vec2
circumcenter Vec2
a Vec2
b Vec2
c = Vec2
a Vec2 -> Vec2 -> Vec2
forall v. VectorSpace v => v -> v -> v
+. Vec2 -> Vec2 -> Vec2 -> Vec2
circumdelta Vec2
a Vec2
b Vec2
c

-- | Check whether a point is inside the circumcircle of a triangle. The triangle
-- must be oriented in counter-clockwise orientation in screen coordinates.
--
-- >>> abc@(a,b,c) = (Vec2 0 0, Vec2 0 100, Vec2 100 0)
--
-- >>> orientation a b c
-- Counterclockwise
--
-- >>> inCircle abc (Vec2 25 25)
-- True
--
-- >>> inCircle abc (Vec2 150 150)
-- False
inCircle
    :: (Vec2, Vec2, Vec2) -- ^ Triangle’s corners
    -> Vec2 -- ^ Point
    -> Bool
inCircle :: (Vec2, Vec2, Vec2) -> Vec2 -> Bool
inCircle (Vec2
a, Vec2
b, Vec2
c) Vec2
p =
    let d :: Vec2
d@(Vec2 Double
dx Double
dy) = Vec2
a Vec2 -> Vec2 -> Vec2
forall v. VectorSpace v => v -> v -> v
-. Vec2
p
        e :: Vec2
e@(Vec2 Double
ex Double
ey) = Vec2
b Vec2 -> Vec2 -> Vec2
forall v. VectorSpace v => v -> v -> v
-. Vec2
p
        f :: Vec2
f@(Vec2 Double
fx Double
fy) = Vec2
c Vec2 -> Vec2 -> Vec2
forall v. VectorSpace v => v -> v -> v
-. Vec2
p

        ap' :: Double
ap' = Vec2 -> Double
normSquare Vec2
d -- ap shadowed by Control.Monad.ap
        bp :: Double
bp  = Vec2 -> Double
normSquare Vec2
e
        cp :: Double
cp  = Vec2 -> Double
normSquare Vec2
f

        val :: Double
val = Double
dxDouble -> Double -> Double
forall a. Num a => a -> a -> a
*(Double
eyDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
cp Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
bpDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
fy) Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
dyDouble -> Double -> Double
forall a. Num a => a -> a -> a
*(Double
exDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
cp Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
bpDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
fx) Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
ap'Double -> Double -> Double
forall a. Num a => a -> a -> a
*(Double
exDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
fy Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
eyDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
fx)

    in Double
val Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
0

-- | Are both axes at most 'epsilon' apart?
nearlyEquals :: Vec2 -> Vec2 -> Bool
nearlyEquals :: Vec2 -> Vec2 -> Bool
nearlyEquals Vec2
a Vec2
b =
    let Vec2 Double
dx Double
dy = Vec2
a Vec2 -> Vec2 -> Vec2
forall v. VectorSpace v => v -> v -> v
-. Vec2
b
    in Double -> Double
forall a. Num a => a -> a
abs Double
dx Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
epsilon Bool -> Bool -> Bool
&& Double -> Double
forall a. Num a => a -> a
abs Double
dy Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
epsilon
  where
    epsilon :: Double
    epsilon :: Double
epsilon = Double
2Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
ieee_float64_epsilon

    ieee_float64_epsilon :: Double
    -- Got this by printing f64::EPSILON in a Rust repl.it.
    -- Couldn’t find it in Haskell’s Base.
    ieee_float64_epsilon :: Double
ieee_float64_epsilon = Double
2.2204460492503131e-16

-- | Represents the area outside of the triangulation. Halfedges on the convex hull
-- (which don't have an adjacent halfedge) will have this value.
tEMPTY :: Int
tEMPTY :: Int
tEMPTY = -Int
1

-- | Given one halfedge, go to the next one. This allows walking around a single
-- triangle.
nextHalfedge :: Int -> Int
nextHalfedge :: Int -> Int
nextHalfedge Int
i = if Int -> Int -> Int
forall a. Integral a => a -> a -> a
mod Int
i Int
3 Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
2 then Int
iInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
2 else Int
iInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1
{-# INLINE nextHalfedge #-}

-- | Inverse of 'nextHalfedge'.
prevHalfedge :: Int -> Int
prevHalfedge :: Int -> Int
prevHalfedge Int
i = if Int -> Int -> Int
forall a. Integral a => a -> a -> a
mod Int
i Int
3 Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0 then Int
iInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
2 else Int
iInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1
{-# INLINE prevHalfedge #-}

data TriangulationST s = TriangulationST
    { forall s. TriangulationST s -> STVector s Int
__triangles :: STVector s Int
    -- ^ A vector of point indices where each triple represents a Delaunay
    -- triangle. All triangles are directed counter-clockwise (in screen
    -- coordinates, i.e. y pointing downwards).
    --
    --   The i-th triangle is points[3i], points[3i+1], points[3i+2].

    , forall s. TriangulationST s -> STRef s Int
__trianglesLen :: STRef s Int
    -- ^ Number of valid entries in '__triangles'. n/3 is the number of triangles
    --   in the triangulation.

    , forall s. TriangulationST s -> STVector s Int
__halfedges :: STVector s Int
    -- ^ A vector of adjacent halfedge indices that allows traversing the triangulation graph.
    --
    -- `i`-th half-edge in the array corresponds to vertex `triangles[i]`
    -- the half-edge is coming from. `halfedges[i]` is the index of a twin half-edge
    -- in an adjacent triangle (or `tEMPTY` for outer halfedges on the convex hull).
    }

-- | Initialize a new mutable triangulation struct.
triangulation_new
    :: HasCallStack
    => Int -- ^ Number of points
    -> ST s (TriangulationST s)
triangulation_new :: forall s. HasCallStack => Int -> ST s (TriangulationST s)
triangulation_new Int
n = do
    let maxTriangles :: Int
maxTriangles = if Int
n Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
2 then Int
2Int -> Int -> Int
forall a. Num a => a -> a -> a
*Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
5 else Int
0
    STVector s Int
triangles <- Int -> ST s (MVector (PrimState (ST s)) Int)
forall (m :: * -> *) a.
PrimMonad m =>
Int -> m (MVector (PrimState m) a)
VM.new (Int
maxTrianglesInt -> Int -> Int
forall a. Num a => a -> a -> a
*Int
3)
    STRef s Int
trianglesLen <- Int -> ST s (STRef s Int)
forall a s. a -> ST s (STRef s a)
newSTRef Int
0
    STVector s Int
halfedges <- Int -> ST s (MVector (PrimState (ST s)) Int)
forall (m :: * -> *) a.
PrimMonad m =>
Int -> m (MVector (PrimState m) a)
VM.new (Int
maxTrianglesInt -> Int -> Int
forall a. Num a => a -> a -> a
*Int
3)
    TriangulationST s -> ST s (TriangulationST s)
forall a. a -> ST s a
forall (f :: * -> *) a. Applicative f => a -> f a
pure TriangulationST
        { __triangles :: STVector s Int
__triangles = STVector s Int
triangles
        , __trianglesLen :: STRef s Int
__trianglesLen = STRef s Int
trianglesLen
        , __halfedges :: STVector s Int
__halfedges = STVector s Int
halfedges
        }

-- | Add a new triangle to the triangulation; report the old (!) size (why, Rust source?!).
triangulation_add_triangle
    :: HasCallStack
    => TriangulationST s
    -> Int -- ^ Corner i0
    -> Int -- ^ Corner i1
    -> Int -- ^ Corner i2
    -> Int -- ^ Halfedge a
    -> Int -- ^ Halfedge b
    -> Int -- ^ Halfedge c
    -> ST s Int
triangulation_add_triangle :: forall s.
HasCallStack =>
TriangulationST s
-> Int -> Int -> Int -> Int -> Int -> Int -> ST s Int
triangulation_add_triangle TriangulationST s
tgl Int
i0 Int
i1 Int
i2 Int
a Int
b Int
c = do
    Int
t <- STRef s Int -> ST s Int
forall s a. STRef s a -> ST s a
readSTRef (TriangulationST s -> STRef s Int
forall s. TriangulationST s -> STRef s Int
__trianglesLen TriangulationST s
tgl)

    MVector (PrimState (ST s)) Int -> Int -> Int -> ST s ()
forall (m :: * -> *) a.
PrimMonad m =>
MVector (PrimState m) a -> Int -> a -> m ()
VM.write (TriangulationST s -> STVector s Int
forall s. TriangulationST s -> STVector s Int
__triangles TriangulationST s
tgl)  Int
t    Int
i0
    MVector (PrimState (ST s)) Int -> Int -> Int -> ST s ()
forall (m :: * -> *) a.
PrimMonad m =>
MVector (PrimState m) a -> Int -> a -> m ()
VM.write (TriangulationST s -> STVector s Int
forall s. TriangulationST s -> STVector s Int
__triangles TriangulationST s
tgl) (Int
tInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1) Int
i1
    MVector (PrimState (ST s)) Int -> Int -> Int -> ST s ()
forall (m :: * -> *) a.
PrimMonad m =>
MVector (PrimState m) a -> Int -> a -> m ()
VM.write (TriangulationST s -> STVector s Int
forall s. TriangulationST s -> STVector s Int
__triangles TriangulationST s
tgl) (Int
tInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
2) Int
i2
    STRef s Int -> (Int -> Int) -> ST s ()
forall s a. STRef s a -> (a -> a) -> ST s ()
modifySTRef' (TriangulationST s -> STRef s Int
forall s. TriangulationST s -> STRef s Int
__trianglesLen TriangulationST s
tgl) (Int -> Int -> Int
forall a. Num a => a -> a -> a
+Int
3)

    MVector (PrimState (ST s)) Int -> Int -> Int -> ST s ()
forall (m :: * -> *) a.
PrimMonad m =>
MVector (PrimState m) a -> Int -> a -> m ()
VM.write (TriangulationST s -> STVector s Int
forall s. TriangulationST s -> STVector s Int
__halfedges TriangulationST s
tgl)  Int
t    Int
a
    MVector (PrimState (ST s)) Int -> Int -> Int -> ST s ()
forall (m :: * -> *) a.
PrimMonad m =>
MVector (PrimState m) a -> Int -> a -> m ()
VM.write (TriangulationST s -> STVector s Int
forall s. TriangulationST s -> STVector s Int
__halfedges TriangulationST s
tgl) (Int
tInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1) Int
b
    MVector (PrimState (ST s)) Int -> Int -> Int -> ST s ()
forall (m :: * -> *) a.
PrimMonad m =>
MVector (PrimState m) a -> Int -> a -> m ()
VM.write (TriangulationST s -> STVector s Int
forall s. TriangulationST s -> STVector s Int
__halfedges TriangulationST s
tgl) (Int
tInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
2) Int
c

    Bool -> ST s () -> ST s ()
forall (f :: * -> *). Applicative f => Bool -> f () -> f ()
when (Int
a Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
tEMPTY) (MVector (PrimState (ST s)) Int -> Int -> Int -> ST s ()
forall (m :: * -> *) a.
PrimMonad m =>
MVector (PrimState m) a -> Int -> a -> m ()
VM.write (TriangulationST s -> STVector s Int
forall s. TriangulationST s -> STVector s Int
__halfedges TriangulationST s
tgl) Int
a  Int
t   )
    Bool -> ST s () -> ST s ()
forall (f :: * -> *). Applicative f => Bool -> f () -> f ()
when (Int
b Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
tEMPTY) (MVector (PrimState (ST s)) Int -> Int -> Int -> ST s ()
forall (m :: * -> *) a.
PrimMonad m =>
MVector (PrimState m) a -> Int -> a -> m ()
VM.write (TriangulationST s -> STVector s Int
forall s. TriangulationST s -> STVector s Int
__halfedges TriangulationST s
tgl) Int
b (Int
tInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1))
    Bool -> ST s () -> ST s ()
forall (f :: * -> *). Applicative f => Bool -> f () -> f ()
when (Int
c Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
tEMPTY) (MVector (PrimState (ST s)) Int -> Int -> Int -> ST s ()
forall (m :: * -> *) a.
PrimMonad m =>
MVector (PrimState m) a -> Int -> a -> m ()
VM.write (TriangulationST s -> STVector s Int
forall s. TriangulationST s -> STVector s Int
__halfedges TriangulationST s
tgl) Int
c (Int
tInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
2))

    Int -> ST s Int
forall a. a -> ST s a
forall (f :: * -> *) a. Applicative f => a -> f a
pure Int
t

-- | If the pair of triangles doesn't satisfy the Delaunay condition (p1 is inside
-- the circumcircle of [p0, pl, pr]), flip them, then do the same check/flip
-- recursively for the new pair of triangles.
--
-- >          pl                    pl
-- >         /||\                  /  \
-- >      al/ || \bl            al/    \a
-- >       /  ||  \              /      \
-- >      /  a||b  \    flip    /___ar___\
-- >    p0\   ||   /p1   =>   p0\---bl---/p1
-- >       \  ||  /              \      /
-- >      ar\ || /br             b\    /br
-- >         \||/                  \  /
-- >          pr                    pr
triangulation_legalize
    :: HasCallStack
    => TriangulationST s -- ^ Triangulation that needs legalization
    -> Int             -- ^ ID of the halfedge to potentially flip
    -> Vector Vec2       -- ^ Delaunay input points
    -> HullST s
    -> ST s Int -- ^ Halfedge adjacent to a. I don’t fully understand this value yet.
triangulation_legalize :: forall s.
HasCallStack =>
TriangulationST s -> Int -> Vector Vec2 -> HullST s -> ST s Int
triangulation_legalize TriangulationST s
tgl !Int
a Vector Vec2
points HullST s
hull = do
    Int
b <- MVector (PrimState (ST s)) Int -> Int -> ST s Int
forall (m :: * -> *) a.
PrimMonad m =>
MVector (PrimState m) a -> Int -> m a
VM.read (TriangulationST s -> STVector s Int
forall s. TriangulationST s -> STVector s Int
__halfedges TriangulationST s
tgl) Int
a
    let ar :: Int
ar = Int -> Int
prevHalfedge Int
a

    -- b is a’s opposite halfedge. If it’s not there (i.e. a is on the hull), we won’t need to flip anything.
    case Int
b Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
tEMPTY of
        Bool
True -> Int -> ST s Int
forall a. a -> ST s a
forall (f :: * -> *) a. Applicative f => a -> f a
pure Int
ar
        Bool
False -> do {
    ; do
        let al :: Int
al = Int -> Int
nextHalfedge Int
a
            bl :: Int
bl = Int -> Int
prevHalfedge Int
b

        Int
p0 <- MVector (PrimState (ST s)) Int -> Int -> ST s Int
forall (m :: * -> *) a.
PrimMonad m =>
MVector (PrimState m) a -> Int -> m a
VM.read (TriangulationST s -> STVector s Int
forall s. TriangulationST s -> STVector s Int
__triangles TriangulationST s
tgl) Int
ar
        Int
pr <- MVector (PrimState (ST s)) Int -> Int -> ST s Int
forall (m :: * -> *) a.
PrimMonad m =>
MVector (PrimState m) a -> Int -> m a
VM.read (TriangulationST s -> STVector s Int
forall s. TriangulationST s -> STVector s Int
__triangles TriangulationST s
tgl) Int
a
        Int
pl <- MVector (PrimState (ST s)) Int -> Int -> ST s Int
forall (m :: * -> *) a.
PrimMonad m =>
MVector (PrimState m) a -> Int -> m a
VM.read (TriangulationST s -> STVector s Int
forall s. TriangulationST s -> STVector s Int
__triangles TriangulationST s
tgl) Int
al
        Int
p1 <- MVector (PrimState (ST s)) Int -> Int -> ST s Int
forall (m :: * -> *) a.
PrimMonad m =>
MVector (PrimState m) a -> Int -> m a
VM.read (TriangulationST s -> STVector s Int
forall s. TriangulationST s -> STVector s Int
__triangles TriangulationST s
tgl) Int
bl

        -- Delaunay condition: No point may be inside the circumcircle of another triangle.
        let illegal :: Bool
illegal = (Vec2, Vec2, Vec2) -> Vec2 -> Bool
inCircle (Vector Vec2
pointsVector Vec2 -> Int -> Vec2
forall a. Vector a -> Int -> a
!Int
p0, Vector Vec2
pointsVector Vec2 -> Int -> Vec2
forall a. Vector a -> Int -> a
!Int
pr, Vector Vec2
pointsVector Vec2 -> Int -> Vec2
forall a. Vector a -> Int -> a
!Int
pl) (Vector Vec2
pointsVector Vec2 -> Int -> Vec2
forall a. Vector a -> Int -> a
!Int
p1)
        case Bool
illegal of
            Bool
False -> Int -> ST s Int
forall a. a -> ST s a
forall (f :: * -> *) a. Applicative f => a -> f a
pure Int
ar
            Bool
True -> do {
    ; do
        MVector (PrimState (ST s)) Int -> Int -> Int -> ST s ()
forall (m :: * -> *) a.
PrimMonad m =>
MVector (PrimState m) a -> Int -> a -> m ()
VM.write (TriangulationST s -> STVector s Int
forall s. TriangulationST s -> STVector s Int
__triangles TriangulationST s
tgl) Int
a Int
p1
        MVector (PrimState (ST s)) Int -> Int -> Int -> ST s ()
forall (m :: * -> *) a.
PrimMonad m =>
MVector (PrimState m) a -> Int -> a -> m ()
VM.write (TriangulationST s -> STVector s Int
forall s. TriangulationST s -> STVector s Int
__triangles TriangulationST s
tgl) Int
b Int
p0

        Int
hbl <- MVector (PrimState (ST s)) Int -> Int -> ST s Int
forall (m :: * -> *) a.
PrimMonad m =>
MVector (PrimState m) a -> Int -> m a
VM.read (TriangulationST s -> STVector s Int
forall s. TriangulationST s -> STVector s Int
__halfedges TriangulationST s
tgl) Int
bl
        Int
har <- MVector (PrimState (ST s)) Int -> Int -> ST s Int
forall (m :: * -> *) a.
PrimMonad m =>
MVector (PrimState m) a -> Int -> m a
VM.read (TriangulationST s -> STVector s Int
forall s. TriangulationST s -> STVector s Int
__halfedges TriangulationST s
tgl) Int
ar

        -- // edge swapped on the other side of the hull (rare); fix the halfedge reference
        Bool -> ST s () -> ST s ()
forall (f :: * -> *). Applicative f => Bool -> f () -> f ()
when (Int
hbl Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
tEMPTY) (ST s () -> ST s ()) -> ST s () -> ST s ()
forall a b. (a -> b) -> a -> b
$ do
            STRef s Int
eRef <- Int -> ST s (STRef s Int)
forall a s. a -> ST s (STRef s a)
newSTRef (Int -> ST s (STRef s Int)) -> ST s Int -> ST s (STRef s Int)
forall (m :: * -> *) a b. Monad m => (a -> m b) -> m a -> m b
=<< STRef s Int -> ST s Int
forall s a. STRef s a -> ST s a
readSTRef (HullST s -> STRef s Int
forall s. HullST s -> STRef s Int
_start HullST s
hull)
            (ST s () -> ST s ()) -> ST s ()
forall a. (a -> a) -> a
fix ((ST s () -> ST s ()) -> ST s ())
-> (ST s () -> ST s ()) -> ST s ()
forall a b. (a -> b) -> a -> b
$ \ST s ()
loop -> do
                Int
e <- STRef s Int -> ST s Int
forall s a. STRef s a -> ST s a
readSTRef STRef s Int
eRef
                Int
hull_tri_e <- MVector (PrimState (ST s)) Int -> Int -> ST s Int
forall (m :: * -> *) a.
PrimMonad m =>
MVector (PrimState m) a -> Int -> m a
VM.read (HullST s -> STVector s Int
forall s. HullST s -> STVector s Int
_tri HullST s
hull) Int
e
                if Int
hull_tri_e Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
bl
                    then do
                        MVector (PrimState (ST s)) Int -> Int -> Int -> ST s ()
forall (m :: * -> *) a.
PrimMonad m =>
MVector (PrimState m) a -> Int -> a -> m ()
VM.write (HullST s -> STVector s Int
forall s. HullST s -> STVector s Int
_tri HullST s
hull) Int
e Int
a
                        () -> ST s ()
forall a. a -> ST s a
forall (f :: * -> *) a. Applicative f => a -> f a
pure () -- break
                    else do
                        Int
hull_prev_e <- MVector (PrimState (ST s)) Int -> Int -> ST s Int
forall (m :: * -> *) a.
PrimMonad m =>
MVector (PrimState m) a -> Int -> m a
VM.read (HullST s -> STVector s Int
forall s. HullST s -> STVector s Int
_prev HullST s
hull) Int
e
                        let e' :: Int
e' = Int
hull_prev_e
                        STRef s Int -> Int -> ST s ()
forall s a. STRef s a -> a -> ST s ()
writeSTRef STRef s Int
eRef Int
e'
                        Int
hull_start <- STRef s Int -> ST s Int
forall s a. STRef s a -> ST s a
readSTRef (HullST s -> STRef s Int
forall s. HullST s -> STRef s Int
_start HullST s
hull)
                        Bool -> ST s () -> ST s ()
forall (f :: * -> *). Applicative f => Bool -> f () -> f ()
when (Int
e' Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
hull_start) ST s ()
loop

        MVector (PrimState (ST s)) Int -> Int -> Int -> ST s ()
forall (m :: * -> *) a.
PrimMonad m =>
MVector (PrimState m) a -> Int -> a -> m ()
VM.write (TriangulationST s -> STVector s Int
forall s. TriangulationST s -> STVector s Int
__halfedges TriangulationST s
tgl) Int
a Int
hbl
        MVector (PrimState (ST s)) Int -> Int -> Int -> ST s ()
forall (m :: * -> *) a.
PrimMonad m =>
MVector (PrimState m) a -> Int -> a -> m ()
VM.write (TriangulationST s -> STVector s Int
forall s. TriangulationST s -> STVector s Int
__halfedges TriangulationST s
tgl) Int
b Int
har
        MVector (PrimState (ST s)) Int -> Int -> Int -> ST s ()
forall (m :: * -> *) a.
PrimMonad m =>
MVector (PrimState m) a -> Int -> a -> m ()
VM.write (TriangulationST s -> STVector s Int
forall s. TriangulationST s -> STVector s Int
__halfedges TriangulationST s
tgl) Int
ar Int
bl

        Bool -> ST s () -> ST s ()
forall (f :: * -> *). Applicative f => Bool -> f () -> f ()
when (Int
hbl Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
tEMPTY) (MVector (PrimState (ST s)) Int -> Int -> Int -> ST s ()
forall (m :: * -> *) a.
PrimMonad m =>
MVector (PrimState m) a -> Int -> a -> m ()
VM.write (TriangulationST s -> STVector s Int
forall s. TriangulationST s -> STVector s Int
__halfedges TriangulationST s
tgl) Int
hbl Int
a)
        Bool -> ST s () -> ST s ()
forall (f :: * -> *). Applicative f => Bool -> f () -> f ()
when (Int
har Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
tEMPTY) (MVector (PrimState (ST s)) Int -> Int -> Int -> ST s ()
forall (m :: * -> *) a.
PrimMonad m =>
MVector (PrimState m) a -> Int -> a -> m ()
VM.write (TriangulationST s -> STVector s Int
forall s. TriangulationST s -> STVector s Int
__halfedges TriangulationST s
tgl) Int
har Int
b)
        Bool -> ST s () -> ST s ()
forall (f :: * -> *). Applicative f => Bool -> f () -> f ()
when (Int
bl Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
tEMPTY) (MVector (PrimState (ST s)) Int -> Int -> Int -> ST s ()
forall (m :: * -> *) a.
PrimMonad m =>
MVector (PrimState m) a -> Int -> a -> m ()
VM.write (TriangulationST s -> STVector s Int
forall s. TriangulationST s -> STVector s Int
__halfedges TriangulationST s
tgl) Int
bl Int
ar)

        let br :: Int
br = Int -> Int
nextHalfedge Int
b
        Int
_ <- TriangulationST s -> Int -> Vector Vec2 -> HullST s -> ST s Int
forall s.
HasCallStack =>
TriangulationST s -> Int -> Vector Vec2 -> HullST s -> ST s Int
triangulation_legalize TriangulationST s
tgl Int
a Vector Vec2
points HullST s
hull
        TriangulationST s -> Int -> Vector Vec2 -> HullST s -> ST s Int
forall s.
HasCallStack =>
TriangulationST s -> Int -> Vector Vec2 -> HullST s -> ST s Int
triangulation_legalize TriangulationST s
tgl Int
br Vector Vec2
points HullST s
hull
    }}

data HullST s = Hull
    { forall s. HullST s -> STVector s Int
_prev :: STVector s Int -- ^ Edge to previous edge
    , forall s. HullST s -> STVector s Int
_next :: STVector s Int -- ^ Edge to next edge
    , forall s. HullST s -> STVector s Int
_tri :: STVector s Int  -- ^ Edge to adjacent halfedge
    , forall s. HullST s -> STVector s Int
_hash :: STVector s Int -- ^ angular edge hash
    , forall s. HullST s -> Int
_hashLen :: Int
    , forall s. HullST s -> STRef s Int
_start :: STRef s Int
    , forall s. HullST s -> Vec2
_center :: Vec2
    }

hull_new
    :: HasCallStack
    => Int       -- ^ Number of points.
                 --   (Redundant since we’ve also got the input points vector,
                 --   but the Rust source does it this way.)
    -> Vec2      -- ^ Circumcenter of the initial triangle
    -> Int       -- ^ First corner of the initial triangle
    -> Int       -- ^ Second corner of the initial triangle
    -> Int       -- ^ Third corner of the initial triangle
    -> Vector Vec2 -- ^ Input points
    -> ST s (HullST s)
hull_new :: forall s.
HasCallStack =>
Int -> Vec2 -> Int -> Int -> Int -> Vector Vec2 -> ST s (HullST s)
hull_new Int
n Vec2
center Int
i0 Int
i1 Int
i2 Vector Vec2
points = do
    let hash_len :: Int
        hash_len :: Int
hash_len = Double -> Int
forall b. Integral b => Double -> b
forall a b. (RealFrac a, Integral b) => a -> b
floor (Double -> Double
forall a. Floating a => a -> a
sqrt (Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n))

    STVector s Int
prev <- Int -> Int -> ST s (MVector (PrimState (ST s)) Int)
forall (m :: * -> *) a.
PrimMonad m =>
Int -> a -> m (MVector (PrimState m) a)
VM.replicate Int
n Int
0
    STVector s Int
next <- Int -> Int -> ST s (MVector (PrimState (ST s)) Int)
forall (m :: * -> *) a.
PrimMonad m =>
Int -> a -> m (MVector (PrimState m) a)
VM.replicate Int
n Int
0
    STVector s Int
tri <- Int -> Int -> ST s (MVector (PrimState (ST s)) Int)
forall (m :: * -> *) a.
PrimMonad m =>
Int -> a -> m (MVector (PrimState m) a)
VM.replicate Int
n Int
0
    STVector s Int
hash <- Int -> Int -> ST s (MVector (PrimState (ST s)) Int)
forall (m :: * -> *) a.
PrimMonad m =>
Int -> a -> m (MVector (PrimState m) a)
VM.replicate Int
hash_len Int
tEMPTY
    STRef s Int
i0Ref <- Int -> ST s (STRef s Int)
forall a s. a -> ST s (STRef s a)
newSTRef Int
i0

    let hull :: HullST s
hull = Hull
            { _prev :: STVector s Int
_prev = STVector s Int
prev
            , _next :: STVector s Int
_next = STVector s Int
next
            , _tri :: STVector s Int
_tri = STVector s Int
tri
            , _hash :: STVector s Int
_hash = STVector s Int
hash
            , _hashLen :: Int
_hashLen = Int
hash_len
            , _start :: STRef s Int
_start = STRef s Int
i0Ref
            , _center :: Vec2
_center = Vec2
center
            }

    -- Forward direction: i0 -> i1 -> i2 -> i0
    MVector (PrimState (ST s)) Int -> Int -> Int -> ST s ()
forall (m :: * -> *) a.
PrimMonad m =>
MVector (PrimState m) a -> Int -> a -> m ()
VM.write STVector s Int
MVector (PrimState (ST s)) Int
next Int
i0 Int
i1
    MVector (PrimState (ST s)) Int -> Int -> Int -> ST s ()
forall (m :: * -> *) a.
PrimMonad m =>
MVector (PrimState m) a -> Int -> a -> m ()
VM.write STVector s Int
MVector (PrimState (ST s)) Int
next Int
i1 Int
i2
    MVector (PrimState (ST s)) Int -> Int -> Int -> ST s ()
forall (m :: * -> *) a.
PrimMonad m =>
MVector (PrimState m) a -> Int -> a -> m ()
VM.write STVector s Int
MVector (PrimState (ST s)) Int
next Int
i2 Int
i0

    -- Backwards direction: i0 -> i2 -> i1 -> i0
    MVector (PrimState (ST s)) Int -> Int -> Int -> ST s ()
forall (m :: * -> *) a.
PrimMonad m =>
MVector (PrimState m) a -> Int -> a -> m ()
VM.write STVector s Int
MVector (PrimState (ST s)) Int
prev Int
i0 Int
i2
    MVector (PrimState (ST s)) Int -> Int -> Int -> ST s ()
forall (m :: * -> *) a.
PrimMonad m =>
MVector (PrimState m) a -> Int -> a -> m ()
VM.write STVector s Int
MVector (PrimState (ST s)) Int
prev Int
i2 Int
i1
    MVector (PrimState (ST s)) Int -> Int -> Int -> ST s ()
forall (m :: * -> *) a.
PrimMonad m =>
MVector (PrimState m) a -> Int -> a -> m ()
VM.write STVector s Int
MVector (PrimState (ST s)) Int
prev Int
i1 Int
i0

    MVector (PrimState (ST s)) Int -> Int -> Int -> ST s ()
forall (m :: * -> *) a.
PrimMonad m =>
MVector (PrimState m) a -> Int -> a -> m ()
VM.write STVector s Int
MVector (PrimState (ST s)) Int
tri Int
i0 Int
0
    MVector (PrimState (ST s)) Int -> Int -> Int -> ST s ()
forall (m :: * -> *) a.
PrimMonad m =>
MVector (PrimState m) a -> Int -> a -> m ()
VM.write STVector s Int
MVector (PrimState (ST s)) Int
tri Int
i1 Int
1
    MVector (PrimState (ST s)) Int -> Int -> Int -> ST s ()
forall (m :: * -> *) a.
PrimMonad m =>
MVector (PrimState m) a -> Int -> a -> m ()
VM.write STVector s Int
MVector (PrimState (ST s)) Int
tri Int
i2 Int
2

    HullST s -> Vec2 -> Int -> ST s ()
forall s. HasCallStack => HullST s -> Vec2 -> Int -> ST s ()
hull_hash_edge HullST s
hull (Vector Vec2
pointsVector Vec2 -> Int -> Vec2
forall a. Vector a -> Int -> a
!Int
i0) Int
i0
    HullST s -> Vec2 -> Int -> ST s ()
forall s. HasCallStack => HullST s -> Vec2 -> Int -> ST s ()
hull_hash_edge HullST s
hull (Vector Vec2
pointsVector Vec2 -> Int -> Vec2
forall a. Vector a -> Int -> a
!Int
i1) Int
i1
    HullST s -> Vec2 -> Int -> ST s ()
forall s. HasCallStack => HullST s -> Vec2 -> Int -> ST s ()
hull_hash_edge HullST s
hull (Vector Vec2
pointsVector Vec2 -> Int -> Vec2
forall a. Vector a -> Int -> a
!Int
i2) Int
i2

    HullST s -> ST s (HullST s)
forall a. a -> ST s a
forall (f :: * -> *) a. Applicative f => a -> f a
pure HullST s
hull

hull_hash_key
    :: HasCallStack
    => HullST s
    -> Vec2
    -> ST s Int
hull_hash_key :: forall s. HasCallStack => HullST s -> Vec2 -> ST s Int
hull_hash_key HullST s
hull Vec2
p = do
    let a :: Double
a = Vec2 -> Double
pseudoAngle (Vec2
p Vec2 -> Vec2 -> Vec2
forall v. VectorSpace v => v -> v -> v
-. HullST s -> Vec2
forall s. HullST s -> Vec2
_center HullST s
hull)
        len :: Int
len = HullST s -> Int
forall s. HullST s -> Int
_hashLen HullST s
hull
    Int -> ST s Int
forall a. a -> ST s a
forall (f :: * -> *) a. Applicative f => a -> f a
pure ((Double -> Int
forall b. Integral b => Double -> b
forall a b. (RealFrac a, Integral b) => a -> b
floor((Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
len :: Double) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
a) :: Int) Int -> Int -> Int
forall a. Integral a => a -> a -> a
`mod` Int
len)

hull_hash_edge
    :: HasCallStack
    => HullST s
    -> Vec2
    -> Int
    -> ST s ()
hull_hash_edge :: forall s. HasCallStack => HullST s -> Vec2 -> Int -> ST s ()
hull_hash_edge HullST s
hull Vec2
p Int
i = do
    Int
key <- HullST s -> Vec2 -> ST s Int
forall s. HasCallStack => HullST s -> Vec2 -> ST s Int
hull_hash_key HullST s
hull Vec2
p
    MVector (PrimState (ST s)) Int -> Int -> Int -> ST s ()
forall (m :: * -> *) a.
PrimMonad m =>
MVector (PrimState m) a -> Int -> a -> m ()
VM.write (HullST s -> STVector s Int
forall s. HullST s -> STVector s Int
_hash HullST s
hull) Int
key Int
i

hull_find_visible_edge
    :: HasCallStack
    => HullST s
    -> Vec2        -- ^ Newly inserted point
    -> Vector Vec2 -- ^ Input points
    -> ST s (Int, Bool)
hull_find_visible_edge :: forall s.
HasCallStack =>
HullST s -> Vec2 -> Vector Vec2 -> ST s (Int, Bool)
hull_find_visible_edge HullST s
hull Vec2
p Vector Vec2
points = do
    STRef s Int
startRef <- Int -> ST s (STRef s Int)
forall a s. a -> ST s (STRef s a)
newSTRef Int
0
    Int
key <- HullST s -> Vec2 -> ST s Int
forall s. HasCallStack => HullST s -> Vec2 -> ST s Int
hull_hash_key HullST s
hull Vec2
p
    let len :: Int
len = HullST s -> Int
forall s. HullST s -> Int
_hashLen HullST s
hull

    -- // find a visible edge on the convex hull using edge hash
    do  let loop :: Int -> ST s ()
loop Int
j | Int
j Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
len = () -> ST s ()
forall a. a -> ST s a
forall (f :: * -> *) a. Applicative f => a -> f a
pure () -- end of loop
            loop Int
j = do
                Int
start <- MVector (PrimState (ST s)) Int -> Int -> ST s Int
forall (m :: * -> *) a.
PrimMonad m =>
MVector (PrimState m) a -> Int -> m a
VM.read (HullST s -> STVector s Int
forall s. HullST s -> STVector s Int
_hash HullST s
hull) ((Int
keyInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
j) Int -> Int -> Int
forall a. Integral a => a -> a -> a
`mod` Int
len)
                STRef s Int -> Int -> ST s ()
forall s a. STRef s a -> a -> ST s ()
writeSTRef STRef s Int
startRef Int
start
                if Int
start Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
tEMPTY
                    then do
                        Int
next <- MVector (PrimState (ST s)) Int -> Int -> ST s Int
forall (m :: * -> *) a.
PrimMonad m =>
MVector (PrimState m) a -> Int -> m a
VM.read (HullST s -> STVector s Int
forall s. HullST s -> STVector s Int
_next HullST s
hull) Int
start
                        if Int
next Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
tEMPTY
                            then () -> ST s ()
forall a. a -> ST s a
forall (f :: * -> *) a. Applicative f => a -> f a
pure () -- break: we found a good start
                            else Int -> ST s ()
loop (Int
jInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1)
                    else Int -> ST s ()
loop (Int
jInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1)
        Int -> ST s ()
loop Int
0

    Int
start <- MVector (PrimState (ST s)) Int -> Int -> ST s Int
forall (m :: * -> *) a.
PrimMonad m =>
MVector (PrimState m) a -> Int -> m a
VM.read (HullST s -> STVector s Int
forall s. HullST s -> STVector s Int
_prev HullST s
hull) (Int -> ST s Int) -> ST s Int -> ST s Int
forall (m :: * -> *) a b. Monad m => (a -> m b) -> m a -> m b
=<< STRef s Int -> ST s Int
forall s a. STRef s a -> ST s a
readSTRef STRef s Int
startRef
    STRef s Int
eRef <- Int -> ST s (STRef s Int)
forall a s. a -> ST s (STRef s a)
newSTRef Int
start

    (ST s (Int, Bool) -> ST s (Int, Bool)) -> ST s (Int, Bool)
forall a. (a -> a) -> a
fix ((ST s (Int, Bool) -> ST s (Int, Bool)) -> ST s (Int, Bool))
-> (ST s (Int, Bool) -> ST s (Int, Bool)) -> ST s (Int, Bool)
forall a b. (a -> b) -> a -> b
$ \ST s (Int, Bool)
loop -> do
        Int
e <- STRef s Int -> ST s Int
forall s a. STRef s a -> ST s a
readSTRef STRef s Int
eRef
        Int
next_e <- MVector (PrimState (ST s)) Int -> Int -> ST s Int
forall (m :: * -> *) a.
PrimMonad m =>
MVector (PrimState m) a -> Int -> m a
VM.read (HullST s -> STVector s Int
forall s. HullST s -> STVector s Int
_next HullST s
hull) Int
e
        if Vec2 -> Vec2 -> Vec2 -> Orientation
orientation Vec2
p (Vector Vec2
pointsVector Vec2 -> Int -> Vec2
forall a. Vector a -> Int -> a
!Int
e) (Vector Vec2
pointsVector Vec2 -> Int -> Vec2
forall a. Vector a -> Int -> a
!Int
next_e) Orientation -> Orientation -> Bool
forall a. Eq a => a -> a -> Bool
/= Orientation
Clockwise -- Rust source: <=0, equivalent to degenerate or counterclockwise
            then do
                STRef s Int -> Int -> ST s ()
forall s a. STRef s a -> a -> ST s ()
writeSTRef STRef s Int
eRef Int
next_e
                if Int
next_e Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
start
                    then (Int, Bool) -> ST s (Int, Bool)
forall a. a -> ST s a
forall (f :: * -> *) a. Applicative f => a -> f a
pure (Int
tEMPTY, Bool
False)
                    else ST s (Int, Bool)
loop
            else (Int, Bool) -> ST s (Int, Bool)
forall a. a -> ST s a
forall (f :: * -> *) a. Applicative f => a -> f a
pure (Int
e, Int
e Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
start)

-- | The seed triangle is the first triangle, located at the center of the input
-- points. Other triangles will be grown around the seed.
findSeedTriangle :: HasCallStack => Vector Vec2 -> Maybe (Int, Int, Int)
findSeedTriangle :: HasCallStack => Vector Vec2 -> Maybe (Int, Int, Int)
findSeedTriangle Vector Vec2
points
    -- Make sure there is some initial triangle at all. Fails when the input starts
    -- with two identical points even though the rest might be legal, too bad for that user.
    | Vector Vec2 -> Int
forall a. Vector a -> Int
V.length Vector Vec2
points Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
3 = Maybe (Int, Int, Int)
forall a. Maybe a
Nothing
    | Double -> Bool
forall a. RealFloat a => a -> Bool
isInfinite (Vec2 -> Vec2 -> Vec2 -> Double
circumradiusSquare (Vector Vec2
pointsVector Vec2 -> Int -> Vec2
forall a. Vector a -> Int -> a
!Int
0) (Vector Vec2
pointsVector Vec2 -> Int -> Vec2
forall a. Vector a -> Int -> a
!Int
1) (Vector Vec2
pointsVector Vec2 -> Int -> Vec2
forall a. Vector a -> Int -> a
!Int
2)) = Maybe (Int, Int, Int)
forall a. Maybe a
Nothing
findSeedTriangle Vector Vec2
points = (Int, Int, Int) -> Maybe (Int, Int, Int)
forall a. a -> Maybe a
Just ((Int, Int, Int) -> Maybe (Int, Int, Int))
-> (Int, Int, Int) -> Maybe (Int, Int, Int)
forall a b. (a -> b) -> a -> b
$

        -- Pick the first point close to the center
    let bboxCenter :: Vec2
bboxCenter = Vector Vec2 -> Vec2
forall a. HasBoundingBox a => a -> Vec2
boundingBoxCenter Vector Vec2
points
        i0 :: Int
i0 = (Vec2 -> Vec2 -> Ordering) -> Vector Vec2 -> Int
forall a. (a -> a -> Ordering) -> Vector a -> Int
V.minIndexBy (\Vec2
p Vec2
q -> (Vec2 -> Double) -> Vec2 -> Vec2 -> Ordering
forall a b. Ord a => (b -> a) -> b -> b -> Ordering
comparing (Vec2 -> Vec2 -> Double
distSquare Vec2
bboxCenter) Vec2
p Vec2
q) Vector Vec2
points
        p0 :: Vec2
p0 = Vector Vec2
pointsVector Vec2 -> Int -> Vec2
forall a. Vector a -> Int -> a
!Int
i0

        -- Find a second point close to the first
        (Int
i1,Vec2
p1) = ((Int, Vec2) -> (Int, Vec2) -> Ordering)
-> Vector (Int, Vec2) -> (Int, Vec2)
forall a. (a -> a -> Ordering) -> Vector a -> a
V.minimumBy (\(Int
_,Vec2
p) (Int
_,Vec2
q) -> (Vec2 -> Double) -> Vec2 -> Vec2 -> Ordering
forall a b. Ord a => (b -> a) -> b -> b -> Ordering
comparing (Vec2 -> Vec2 -> Double
distSquare Vec2
p0) Vec2
p Vec2
q) (Vector (Int, Vec2) -> (Int, Vec2))
-> (Vector (Int, Vec2) -> Vector (Int, Vec2))
-> Vector (Int, Vec2)
-> (Int, Vec2)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. ((Int, Vec2) -> Bool) -> Vector (Int, Vec2) -> Vector (Int, Vec2)
forall a. (a -> Bool) -> Vector a -> Vector a
V.filter (\(Int
i,Vec2
_) -> Int
i Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
i0) (Vector (Int, Vec2) -> (Int, Vec2))
-> Vector (Int, Vec2) -> (Int, Vec2)
forall a b. (a -> b) -> a -> b
$ Vector Vec2 -> Vector (Int, Vec2)
forall a. Vector a -> Vector (Int, a)
V.indexed Vector Vec2
points

        -- Find the third point which forms the smallest circumcircle with the first two
        (Int
i2,Vec2
p2) = ((Int, Vec2) -> (Int, Vec2) -> Ordering)
-> Vector (Int, Vec2) -> (Int, Vec2)
forall a. (a -> a -> Ordering) -> Vector a -> a
V.minimumBy (\(Int
_,Vec2
p) (Int
_,Vec2
q) -> (Vec2 -> Double) -> Vec2 -> Vec2 -> Ordering
forall a b. Ord a => (b -> a) -> b -> b -> Ordering
comparing (Vec2 -> Vec2 -> Vec2 -> Double
circumradiusSquare Vec2
p0 Vec2
p1) Vec2
p Vec2
q) (Vector (Int, Vec2) -> (Int, Vec2))
-> (Vector Vec2 -> Vector (Int, Vec2))
-> Vector Vec2
-> (Int, Vec2)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. ((Int, Vec2) -> Bool) -> Vector (Int, Vec2) -> Vector (Int, Vec2)
forall a. (a -> Bool) -> Vector a -> Vector a
V.filter (\(Int
i,Vec2
_) -> Int
i Int -> [Int] -> Bool
forall (t :: * -> *) a. (Foldable t, Eq a) => a -> t a -> Bool
`notElem` [Int
i0,Int
i1]) (Vector (Int, Vec2) -> Vector (Int, Vec2))
-> (Vector Vec2 -> Vector (Int, Vec2))
-> Vector Vec2
-> Vector (Int, Vec2)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Vector Vec2 -> Vector (Int, Vec2)
forall a. Vector a -> Vector (Int, a)
V.indexed (Vector Vec2 -> (Int, Vec2)) -> Vector Vec2 -> (Int, Vec2)
forall a b. (a -> b) -> a -> b
$ Vector Vec2
points

    in if Vec2 -> Vec2 -> Vec2 -> Orientation
orientation Vec2
p0 Vec2
p1 Vec2
p2 Orientation -> Orientation -> Bool
forall a. Eq a => a -> a -> Bool
== Orientation
Clockwise
        then (Int
i0, Int
i2, Int
i1) -- Flip two points to ensure counterclockwise orientation
        else (Int
i0, Int
i1, Int
i2)

-- | Inplace sort of the input by distance.
--
-- __Rust source:__ @sortf@
sortf :: STVector s (a, Double) -> ST s ()
sortf :: forall s a. STVector s (a, Double) -> ST s ()
sortf = Comparison (a, Double)
-> MVector (PrimState (ST s)) (a, Double) -> ST s ()
forall (m :: * -> *) (v :: * -> * -> *) e.
(PrimMonad m, MVector v e) =>
Comparison e -> v (PrimState m) e -> m ()
VM.sortBy (((a, Double) -> Double) -> Comparison (a, Double)
forall a b. Ord a => (b -> a) -> b -> b -> Ordering
comparing (\(a
_, Double
d) -> Double
d))

-- /// Order collinear points by dx (or dy if all x are identical) and return the list as a hull
-- fn handle_collinear_points(points: &[Point]) -> Triangulation {
--     let Point { x, y } = points.first().cloned().unwrap_or_default();
--
--     let mut dist: Vec<_> = points
--         .iter()
--         .enumerate()
--         .map(|(i, p)| {
--             let mut d = p.x - x;
--             if d == 0.0 {
--                 d = p.y - y;
--             }
--             (i, d)
--         })
--         .collect();
--     sortf(&mut dist);
--
--     let mut triangulation = Triangulation::new(0);
--     let mut d0 = f64::NEG_INFINITY;
--     for (i, distance) in dist {
--         if distance > d0 {
--             triangulation.hull.push(i);
--             d0 = distance;
--         }
--     }
--
--     triangulation
-- }

-- | Triangulate a set of 2D points. Returns the triangulation for the input
-- points. For the degenerated case when all points are collinear, returns an empty
-- triangulation where all points are in the hull.
--
-- __Rust source:__ @triangulate@
triangulateMut :: HasCallStack => Vector Vec2 -> ST s (TriangulationST s, HullST s)
triangulateMut :: forall s.
HasCallStack =>
Vector Vec2 -> ST s (TriangulationST s, HullST s)
triangulateMut Vector Vec2
points = do
    case HasCallStack => Vector Vec2 -> Maybe (Int, Int, Int)
Vector Vec2 -> Maybe (Int, Int, Int)
findSeedTriangle Vector Vec2
points of
        Maybe (Int, Int, Int)
Nothing -> String -> ST s (TriangulationST s, HullST s)
forall a. HasCallStack => String -> a
error String
"Can’t find a seed triangle, and handle_collinear_points is not implemented" -- TODO!
        Just (Int, Int, Int)
seed_triangle -> (Int, Int, Int) -> ST s (TriangulationST s, HullST s)
forall s. (Int, Int, Int) -> ST s (TriangulationST s, HullST s)
triangulate_for_real (Int, Int, Int)
seed_triangle
  where
    triangulate_for_real :: (Int, Int, Int) -> ST s (TriangulationST s, HullST s)
    triangulate_for_real :: forall s. (Int, Int, Int) -> ST s (TriangulationST s, HullST s)
triangulate_for_real (Int, Int, Int)
seed_triangle = do
        let numPoints :: Int
numPoints = Vector Vec2 -> Int
forall a. Vector a -> Int
V.length Vector Vec2
points
            (Int
i0, Int
i1, Int
i2) = (Int, Int, Int)
seed_triangle
            center :: Vec2
center = Vec2 -> Vec2 -> Vec2 -> Vec2
circumcenter (Vector Vec2
pointsVector Vec2 -> Int -> Vec2
forall a. Vector a -> Int -> a
!Int
i0) (Vector Vec2
pointsVector Vec2 -> Int -> Vec2
forall a. Vector a -> Int -> a
!Int
i1) (Vector Vec2
pointsVector Vec2 -> Int -> Vec2
forall a. Vector a -> Int -> a
!Int
i2)

        TriangulationST s
tgl <- Int -> ST s (TriangulationST s)
forall s. HasCallStack => Int -> ST s (TriangulationST s)
triangulation_new Int
numPoints
        Int
_ <- TriangulationST s
-> Int -> Int -> Int -> Int -> Int -> Int -> ST s Int
forall s.
HasCallStack =>
TriangulationST s
-> Int -> Int -> Int -> Int -> Int -> Int -> ST s Int
triangulation_add_triangle TriangulationST s
tgl Int
i0 Int
i1 Int
i2 Int
tEMPTY Int
tEMPTY Int
tEMPTY

        -- // sort the points by distance from the seed triangle circumcenter
        let dists :: Vector (Int, Double)
dists = (forall s. MVector s (Int, Double) -> ST s ())
-> Vector (Int, Double) -> Vector (Int, Double)
forall a.
(forall s. MVector s a -> ST s ()) -> Vector a -> Vector a
V.modify STVector s (Int, Double) -> ST s ()
forall s. MVector s (Int, Double) -> ST s ()
forall s a. STVector s (a, Double) -> ST s ()
sortf ((Int -> Vec2 -> (Int, Double))
-> Vector Vec2 -> Vector (Int, Double)
forall a b. (Int -> a -> b) -> Vector a -> Vector b
V.imap (\Int
i Vec2
p -> (Int
i, Vec2 -> Vec2 -> Double
distSquare Vec2
center Vec2
p)) Vector Vec2
points)

        HullST s
hull <- Int -> Vec2 -> Int -> Int -> Int -> Vector Vec2 -> ST s (HullST s)
forall s.
HasCallStack =>
Int -> Vec2 -> Int -> Int -> Int -> Vector Vec2 -> ST s (HullST s)
hull_new Int
numPoints Vec2
center Int
i0 Int
i1 Int
i2 Vector Vec2
points

        -- Those indices, argh!
        -- k: the loop runs over the k-th closest input point to the seed triangle center.
        -- i: index of the current subject point, as in »we are looking at point[i]«
        Vector (Int, Double)
-> (Int -> (Int, Double) -> ST s ()) -> ST s ()
forall (m :: * -> *) a b.
Monad m =>
Vector a -> (Int -> a -> m b) -> m ()
V.iforM_ Vector (Int, Double)
dists ((Int -> (Int, Double) -> ST s ()) -> ST s ())
-> (Int -> (Int, Double) -> ST s ()) -> ST s ()
forall a b. (a -> b) -> a -> b
$ \Int
k (Int
i, Double
_) -> do
            let p :: Vec2
p = Vector Vec2
pointsVector Vec2 -> Int -> Vec2
forall a. Vector a -> Int -> a
!Int
i
                isNearDuplicate :: Bool
isNearDuplicate =
                    let (Int
previousI, Double
_previousDist) = Vector (Int, Double)
distsVector (Int, Double) -> Int -> (Int, Double)
forall a. Vector a -> Int -> a
!(Int
kInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1)
                        previousP :: Vec2
previousP = Vector Vec2
pointsVector Vec2 -> Int -> Vec2
forall a. Vector a -> Int -> a
!Int
previousI
                    in Int
k Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
0 Bool -> Bool -> Bool
&& Vec2 -> Vec2 -> Bool
nearlyEquals Vec2
p Vec2
previousP -- k>0 so we don’t underrun in the k-1 lookup
                isInSeedTriangle :: Bool
isInSeedTriangle = Int
i Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
i0 Bool -> Bool -> Bool
|| Int
i Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
i1 Bool -> Bool -> Bool
|| Int
i Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
i2

            -- // skip near-duplicates
            -- // skip seed triangle points
            if Bool
isNearDuplicate Bool -> Bool -> Bool
|| Bool
isInSeedTriangle
                then () -> ST s ()
forall a. a -> ST s a
forall (f :: * -> *) a. Applicative f => a -> f a
pure () -- continue
                else do

                    -- // find a visible edge on the convex hull using edge hash
                    (Int
e0, Bool
walk_back) <- HullST s -> Vec2 -> Vector Vec2 -> ST s (Int, Bool)
forall s.
HasCallStack =>
HullST s -> Vec2 -> Vector Vec2 -> ST s (Int, Bool)
hull_find_visible_edge HullST s
hull Vec2
p Vector Vec2
points

                    if Int
e0 Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
tEMPTY
                        then () -> ST s ()
forall a. a -> ST s a
forall (f :: * -> *) a. Applicative f => a -> f a
pure () -- // likely a near-duplicate point; skip it
                        else HullST s
-> TriangulationST s -> Int -> Bool -> Int -> Vec2 -> ST s ()
forall s.
HullST s
-> TriangulationST s -> Int -> Bool -> Int -> Vec2 -> ST s ()
addPoint HullST s
hull TriangulationST s
tgl Int
e0 Bool
walk_back Int
i Vec2
p

        (TriangulationST s, HullST s) -> ST s (TriangulationST s, HullST s)
forall a. a -> ST s a
forall (f :: * -> *) a. Applicative f => a -> f a
pure (TriangulationST s
tgl, HullST s
hull)

    addPoint
        :: HullST s
        -> TriangulationST s
        -> Int  -- ^ Visible edge from p
        -> Bool -- ^ walk_back
        -> Int  -- ^ Index of p
        -> Vec2 -- ^ Coordinates of p
        -> ST s ()
    addPoint :: forall s.
HullST s
-> TriangulationST s -> Int -> Bool -> Int -> Vec2 -> ST s ()
addPoint HullST s
hull TriangulationST s
tgl Int
e0 Bool
walk_back Int
i Vec2
p = do

        STRef s Int
eRef <- Int -> ST s (STRef s Int)
forall a s. a -> ST s (STRef s a)
newSTRef Int
e0
        -- // add the first triangle from the point
        do
            Int
e <- STRef s Int -> ST s Int
forall s a. STRef s a -> ST s a
readSTRef STRef s Int
eRef
            Int
t <- do
                Int
hull_next_e <- MVector (PrimState (ST s)) Int -> Int -> ST s Int
forall (m :: * -> *) a.
PrimMonad m =>
MVector (PrimState m) a -> Int -> m a
VM.read (HullST s -> STVector s Int
forall s. HullST s -> STVector s Int
_next HullST s
hull) Int
e
                Int
hull_tri_e <- MVector (PrimState (ST s)) Int -> Int -> ST s Int
forall (m :: * -> *) a.
PrimMonad m =>
MVector (PrimState m) a -> Int -> m a
VM.read (HullST s -> STVector s Int
forall s. HullST s -> STVector s Int
_tri HullST s
hull) Int
e
                TriangulationST s
-> Int -> Int -> Int -> Int -> Int -> Int -> ST s Int
forall s.
HasCallStack =>
TriangulationST s
-> Int -> Int -> Int -> Int -> Int -> Int -> ST s Int
triangulation_add_triangle TriangulationST s
tgl Int
e Int
i Int
hull_next_e Int
tEMPTY Int
tEMPTY Int
hull_tri_e

            -- // recursively flip triangles from the point until they satisfy the Delaunay condition
            MVector (PrimState (ST s)) Int -> Int -> Int -> ST s ()
forall (m :: * -> *) a.
PrimMonad m =>
MVector (PrimState m) a -> Int -> a -> m ()
VM.write (HullST s -> STVector s Int
forall s. HullST s -> STVector s Int
_tri HullST s
hull) Int
i (Int -> ST s ()) -> ST s Int -> ST s ()
forall (m :: * -> *) a b. Monad m => (a -> m b) -> m a -> m b
=<< TriangulationST s -> Int -> Vector Vec2 -> HullST s -> ST s Int
forall s.
HasCallStack =>
TriangulationST s -> Int -> Vector Vec2 -> HullST s -> ST s Int
triangulation_legalize TriangulationST s
tgl (Int
tInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
2) Vector Vec2
points HullST s
hull
            MVector (PrimState (ST s)) Int -> Int -> Int -> ST s ()
forall (m :: * -> *) a.
PrimMonad m =>
MVector (PrimState m) a -> Int -> a -> m ()
VM.write (HullST s -> STVector s Int
forall s. HullST s -> STVector s Int
_tri HullST s
hull) Int
e Int
t -- // keep track of boundary triangles on the hull

        -- // walk forward through the hull, adding more triangles and flipping recursively
        STRef s Int
nRef <- do
            Int
e <- STRef s Int -> ST s Int
forall s a. STRef s a -> ST s a
readSTRef STRef s Int
eRef
            Int
hull_next_e <- MVector (PrimState (ST s)) Int -> Int -> ST s Int
forall (m :: * -> *) a.
PrimMonad m =>
MVector (PrimState m) a -> Int -> m a
VM.read (HullST s -> STVector s Int
forall s. HullST s -> STVector s Int
_next HullST s
hull) Int
e
            Int -> ST s (STRef s Int)
forall a s. a -> ST s (STRef s a)
newSTRef Int
hull_next_e
        (ST s () -> ST s ()) -> ST s ()
forall a. (a -> a) -> a
fix ((ST s () -> ST s ()) -> ST s ())
-> (ST s () -> ST s ()) -> ST s ()
forall a b. (a -> b) -> a -> b
$ \ST s ()
loop -> do
            Int
n <- STRef s Int -> ST s Int
forall s a. STRef s a -> ST s a
readSTRef STRef s Int
nRef
            Int
q <- MVector (PrimState (ST s)) Int -> Int -> ST s Int
forall (m :: * -> *) a.
PrimMonad m =>
MVector (PrimState m) a -> Int -> m a
VM.read (HullST s -> STVector s Int
forall s. HullST s -> STVector s Int
_next HullST s
hull) Int
n
            if Vec2 -> Vec2 -> Vec2 -> Orientation
orientation Vec2
p (Vector Vec2
pointsVector Vec2 -> Int -> Vec2
forall a. Vector a -> Int -> a
!Int
n) (Vector Vec2
pointsVector Vec2 -> Int -> Vec2
forall a. Vector a -> Int -> a
!Int
q) Orientation -> Orientation -> Bool
forall a. Eq a => a -> a -> Bool
/= Orientation
Clockwise -- Rust source: <=0, equivalent to degenerate-or-counterclockwise
                then () -> ST s ()
forall a. a -> ST s a
forall (f :: * -> *) a. Applicative f => a -> f a
pure () -- break
                else do
                    Int
hull_tri_i <- MVector (PrimState (ST s)) Int -> Int -> ST s Int
forall (m :: * -> *) a.
PrimMonad m =>
MVector (PrimState m) a -> Int -> m a
VM.read (HullST s -> STVector s Int
forall s. HullST s -> STVector s Int
_tri HullST s
hull) Int
i
                    Int
hull_tri_n <- MVector (PrimState (ST s)) Int -> Int -> ST s Int
forall (m :: * -> *) a.
PrimMonad m =>
MVector (PrimState m) a -> Int -> m a
VM.read (HullST s -> STVector s Int
forall s. HullST s -> STVector s Int
_tri HullST s
hull) Int
n
                    Int
t <- TriangulationST s
-> Int -> Int -> Int -> Int -> Int -> Int -> ST s Int
forall s.
HasCallStack =>
TriangulationST s
-> Int -> Int -> Int -> Int -> Int -> Int -> ST s Int
triangulation_add_triangle TriangulationST s
tgl Int
n Int
i Int
q Int
hull_tri_i Int
tEMPTY Int
hull_tri_n
                    MVector (PrimState (ST s)) Int -> Int -> Int -> ST s ()
forall (m :: * -> *) a.
PrimMonad m =>
MVector (PrimState m) a -> Int -> a -> m ()
VM.write (HullST s -> STVector s Int
forall s. HullST s -> STVector s Int
_tri HullST s
hull) Int
i (Int -> ST s ()) -> ST s Int -> ST s ()
forall (m :: * -> *) a b. Monad m => (a -> m b) -> m a -> m b
=<< TriangulationST s -> Int -> Vector Vec2 -> HullST s -> ST s Int
forall s.
HasCallStack =>
TriangulationST s -> Int -> Vector Vec2 -> HullST s -> ST s Int
triangulation_legalize TriangulationST s
tgl (Int
tInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
2) Vector Vec2
points HullST s
hull
                    MVector (PrimState (ST s)) Int -> Int -> Int -> ST s ()
forall (m :: * -> *) a.
PrimMonad m =>
MVector (PrimState m) a -> Int -> a -> m ()
VM.write (HullST s -> STVector s Int
forall s. HullST s -> STVector s Int
_next HullST s
hull) Int
n Int
tEMPTY -- // mark as removed
                    STRef s Int -> Int -> ST s ()
forall s a. STRef s a -> a -> ST s ()
writeSTRef STRef s Int
nRef Int
q
                    ST s ()
loop

        -- // walk backward from the other side, adding more triangles and flipping
        Bool -> ST s () -> ST s ()
forall (f :: * -> *). Applicative f => Bool -> f () -> f ()
when Bool
walk_back (ST s () -> ST s ()) -> ST s () -> ST s ()
forall a b. (a -> b) -> a -> b
$ do
            (ST s () -> ST s ()) -> ST s ()
forall a. (a -> a) -> a
fix ((ST s () -> ST s ()) -> ST s ())
-> (ST s () -> ST s ()) -> ST s ()
forall a b. (a -> b) -> a -> b
$ \ST s ()
loop -> do
                Int
e <- STRef s Int -> ST s Int
forall s a. STRef s a -> ST s a
readSTRef STRef s Int
eRef
                Int
q <- MVector (PrimState (ST s)) Int -> Int -> ST s Int
forall (m :: * -> *) a.
PrimMonad m =>
MVector (PrimState m) a -> Int -> m a
VM.read (HullST s -> STVector s Int
forall s. HullST s -> STVector s Int
_prev HullST s
hull) Int
e
                if Vec2 -> Vec2 -> Vec2 -> Orientation
orientation Vec2
p (Vector Vec2
pointsVector Vec2 -> Int -> Vec2
forall a. Vector a -> Int -> a
!Int
q) (Vector Vec2
pointsVector Vec2 -> Int -> Vec2
forall a. Vector a -> Int -> a
!Int
e) Orientation -> Orientation -> Bool
forall a. Eq a => a -> a -> Bool
/= Orientation
Clockwise -- Rust source: <=0, equivalent to degenerate-or-counterclockwise
                    then () -> ST s ()
forall a. a -> ST s a
forall (f :: * -> *) a. Applicative f => a -> f a
pure () -- break
                    else do
                        Int
hull_tri_e <- MVector (PrimState (ST s)) Int -> Int -> ST s Int
forall (m :: * -> *) a.
PrimMonad m =>
MVector (PrimState m) a -> Int -> m a
VM.read (HullST s -> STVector s Int
forall s. HullST s -> STVector s Int
_tri HullST s
hull) Int
e
                        Int
hull_tri_q <- MVector (PrimState (ST s)) Int -> Int -> ST s Int
forall (m :: * -> *) a.
PrimMonad m =>
MVector (PrimState m) a -> Int -> m a
VM.read (HullST s -> STVector s Int
forall s. HullST s -> STVector s Int
_tri HullST s
hull) Int
q
                        Int
t <- TriangulationST s
-> Int -> Int -> Int -> Int -> Int -> Int -> ST s Int
forall s.
HasCallStack =>
TriangulationST s
-> Int -> Int -> Int -> Int -> Int -> Int -> ST s Int
triangulation_add_triangle TriangulationST s
tgl Int
q Int
i Int
e Int
tEMPTY Int
hull_tri_e Int
hull_tri_q
                        Int
_ <- TriangulationST s -> Int -> Vector Vec2 -> HullST s -> ST s Int
forall s.
HasCallStack =>
TriangulationST s -> Int -> Vector Vec2 -> HullST s -> ST s Int
triangulation_legalize TriangulationST s
tgl (Int
tInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
2) Vector Vec2
points HullST s
hull
                        MVector (PrimState (ST s)) Int -> Int -> Int -> ST s ()
forall (m :: * -> *) a.
PrimMonad m =>
MVector (PrimState m) a -> Int -> a -> m ()
VM.write (HullST s -> STVector s Int
forall s. HullST s -> STVector s Int
_tri HullST s
hull) Int
q Int
t
                        MVector (PrimState (ST s)) Int -> Int -> Int -> ST s ()
forall (m :: * -> *) a.
PrimMonad m =>
MVector (PrimState m) a -> Int -> a -> m ()
VM.write (HullST s -> STVector s Int
forall s. HullST s -> STVector s Int
_next HullST s
hull) Int
e Int
tEMPTY -- // mark as removed
                        STRef s Int -> Int -> ST s ()
forall s a. STRef s a -> a -> ST s ()
writeSTRef STRef s Int
eRef Int
q
                        ST s ()
loop

        -- // update the hull indices
        Int
e <- STRef s Int -> ST s Int
forall s a. STRef s a -> ST s a
readSTRef STRef s Int
eRef
        Int
n <- STRef s Int -> ST s Int
forall s a. STRef s a -> ST s a
readSTRef STRef s Int
nRef
        MVector (PrimState (ST s)) Int -> Int -> Int -> ST s ()
forall (m :: * -> *) a.
PrimMonad m =>
MVector (PrimState m) a -> Int -> a -> m ()
VM.write (HullST s -> STVector s Int
forall s. HullST s -> STVector s Int
_prev HullST s
hull) Int
i Int
e
        MVector (PrimState (ST s)) Int -> Int -> Int -> ST s ()
forall (m :: * -> *) a.
PrimMonad m =>
MVector (PrimState m) a -> Int -> a -> m ()
VM.write (HullST s -> STVector s Int
forall s. HullST s -> STVector s Int
_next HullST s
hull) Int
i Int
n
        MVector (PrimState (ST s)) Int -> Int -> Int -> ST s ()
forall (m :: * -> *) a.
PrimMonad m =>
MVector (PrimState m) a -> Int -> a -> m ()
VM.write (HullST s -> STVector s Int
forall s. HullST s -> STVector s Int
_prev HullST s
hull) Int
n Int
i
        MVector (PrimState (ST s)) Int -> Int -> Int -> ST s ()
forall (m :: * -> *) a.
PrimMonad m =>
MVector (PrimState m) a -> Int -> a -> m ()
VM.write (HullST s -> STVector s Int
forall s. HullST s -> STVector s Int
_next HullST s
hull) Int
e Int
i
        STRef s Int -> Int -> ST s ()
forall s a. STRef s a -> a -> ST s ()
writeSTRef (HullST s -> STRef s Int
forall s. HullST s -> STRef s Int
_start HullST s
hull) Int
e

        -- // expose hull as a vector of point indices
        -- ==> This is done in 'freezeHull' outside of this function.

        -- // save the two new edges in the hash table
        HullST s -> Vec2 -> Int -> ST s ()
forall s. HasCallStack => HullST s -> Vec2 -> Int -> ST s ()
hull_hash_edge HullST s
hull Vec2
p Int
i
        HullST s -> Vec2 -> Int -> ST s ()
forall s. HasCallStack => HullST s -> Vec2 -> Int -> ST s ()
hull_hash_edge HullST s
hull (Vector Vec2
pointsVector Vec2 -> Int -> Vec2
forall a. Vector a -> Int -> a
!Int
e) Int
e

-- | Create the list of hull points based on the previously calculated hull data.
freezeHull :: HullST s -> ST s (Vector Int)
freezeHull :: forall s. HullST s -> ST s (Vector Int)
freezeHull HullST s
hull = do
    Int
hull_start <- STRef s Int -> ST s Int
forall s a. STRef s a -> ST s a
readSTRef (HullST s -> STRef s Int
forall s. HullST s -> STRef s Int
_start HullST s
hull)
    ((Int, Bool) -> ST s (Maybe (Int, (Int, Bool))))
-> (Int, Bool) -> ST s (Vector Int)
forall (m :: * -> *) b a.
Monad m =>
(b -> m (Maybe (a, b))) -> b -> m (Vector a)
V.unfoldrM
        -- We add an 'isFirst' flag here so we can distinguish the cases
        --    1. We start with e == hull_start but want to continue
        --    2. We have e == hull_start after circling the hull and want to terminate
        (\(Int
e, Bool
isFirst) -> do
            Int
next_e <- MVector (PrimState (ST s)) Int -> Int -> ST s Int
forall (m :: * -> *) a.
PrimMonad m =>
MVector (PrimState m) a -> Int -> m a
VM.read (HullST s -> STVector s Int
forall s. HullST s -> STVector s Int
_next HullST s
hull) Int
e -- prev/next switch hull orientation here. Pick whatever you like.
            if Bool -> Bool
not Bool
isFirst Bool -> Bool -> Bool
&& Int
e Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
hull_start
                then Maybe (Int, (Int, Bool)) -> ST s (Maybe (Int, (Int, Bool)))
forall a. a -> ST s a
forall (f :: * -> *) a. Applicative f => a -> f a
pure Maybe (Int, (Int, Bool))
forall a. Maybe a
Nothing
                else Maybe (Int, (Int, Bool)) -> ST s (Maybe (Int, (Int, Bool)))
forall a. a -> ST s a
forall (f :: * -> *) a. Applicative f => a -> f a
pure ((Int, (Int, Bool)) -> Maybe (Int, (Int, Bool))
forall a. a -> Maybe a
Just (Int
e, (Int
next_e, Bool
False)))
        )
        (Int
hull_start, Bool
True)

-- | User-facing raw Delaunay triangulation data. Build the API based on this and
-- 'triangulate'. See 'TriangulationST' for documentation of the fields.
data TriangulationRaw = TriangulationRaw
    { TriangulationRaw -> Vector Int
_triangles :: Vector Int
    , TriangulationRaw -> Vector Int
_halfedges :: Vector Int
    , TriangulationRaw -> Vector Int
_convexHull :: Vector Int
    } deriving (TriangulationRaw -> TriangulationRaw -> Bool
(TriangulationRaw -> TriangulationRaw -> Bool)
-> (TriangulationRaw -> TriangulationRaw -> Bool)
-> Eq TriangulationRaw
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: TriangulationRaw -> TriangulationRaw -> Bool
== :: TriangulationRaw -> TriangulationRaw -> Bool
$c/= :: TriangulationRaw -> TriangulationRaw -> Bool
/= :: TriangulationRaw -> TriangulationRaw -> Bool
Eq, Eq TriangulationRaw
Eq TriangulationRaw
-> (TriangulationRaw -> TriangulationRaw -> Ordering)
-> (TriangulationRaw -> TriangulationRaw -> Bool)
-> (TriangulationRaw -> TriangulationRaw -> Bool)
-> (TriangulationRaw -> TriangulationRaw -> Bool)
-> (TriangulationRaw -> TriangulationRaw -> Bool)
-> (TriangulationRaw -> TriangulationRaw -> TriangulationRaw)
-> (TriangulationRaw -> TriangulationRaw -> TriangulationRaw)
-> Ord TriangulationRaw
TriangulationRaw -> TriangulationRaw -> Bool
TriangulationRaw -> TriangulationRaw -> Ordering
TriangulationRaw -> TriangulationRaw -> TriangulationRaw
forall a.
Eq a
-> (a -> a -> Ordering)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> a)
-> (a -> a -> a)
-> Ord a
$ccompare :: TriangulationRaw -> TriangulationRaw -> Ordering
compare :: TriangulationRaw -> TriangulationRaw -> Ordering
$c< :: TriangulationRaw -> TriangulationRaw -> Bool
< :: TriangulationRaw -> TriangulationRaw -> Bool
$c<= :: TriangulationRaw -> TriangulationRaw -> Bool
<= :: TriangulationRaw -> TriangulationRaw -> Bool
$c> :: TriangulationRaw -> TriangulationRaw -> Bool
> :: TriangulationRaw -> TriangulationRaw -> Bool
$c>= :: TriangulationRaw -> TriangulationRaw -> Bool
>= :: TriangulationRaw -> TriangulationRaw -> Bool
$cmax :: TriangulationRaw -> TriangulationRaw -> TriangulationRaw
max :: TriangulationRaw -> TriangulationRaw -> TriangulationRaw
$cmin :: TriangulationRaw -> TriangulationRaw -> TriangulationRaw
min :: TriangulationRaw -> TriangulationRaw -> TriangulationRaw
Ord, Int -> TriangulationRaw -> ShowS
[TriangulationRaw] -> ShowS
TriangulationRaw -> String
(Int -> TriangulationRaw -> ShowS)
-> (TriangulationRaw -> String)
-> ([TriangulationRaw] -> ShowS)
-> Show TriangulationRaw
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> TriangulationRaw -> ShowS
showsPrec :: Int -> TriangulationRaw -> ShowS
$cshow :: TriangulationRaw -> String
show :: TriangulationRaw -> String
$cshowList :: [TriangulationRaw] -> ShowS
showList :: [TriangulationRaw] -> ShowS
Show)

instance NFData TriangulationRaw where
    rnf :: TriangulationRaw -> ()
rnf (TriangulationRaw Vector Int
a Vector Int
b Vector Int
c) = (Vector Int, Vector Int, Vector Int) -> ()
forall a. NFData a => a -> ()
rnf (Vector Int
a,Vector Int
b,Vector Int
c)

-- | NB: unsafeFreeze is safe here, since this module encapsulates the mutable
-- triangulation along with 'triangulation_new'.
freezeTriangulation :: HasCallStack => TriangulationST s -> HullST s -> ST s TriangulationRaw
freezeTriangulation :: forall s.
HasCallStack =>
TriangulationST s -> HullST s -> ST s TriangulationRaw
freezeTriangulation TriangulationST s
tglMut HullST s
hullMut = do
    Int
trianglesLen <- STRef s Int -> ST s Int
forall s a. STRef s a -> ST s a
readSTRef (TriangulationST s -> STRef s Int
forall s. TriangulationST s -> STRef s Int
__trianglesLen TriangulationST s
tglMut)
    Vector Int
triangles <- MVector (PrimState (ST s)) Int -> ST s (Vector Int)
forall (m :: * -> *) a.
PrimMonad m =>
MVector (PrimState m) a -> m (Vector a)
V.unsafeFreeze (Int -> MVector s Int -> MVector s Int
forall s a. Int -> MVector s a -> MVector s a
VM.take Int
trianglesLen (TriangulationST s -> MVector s Int
forall s. TriangulationST s -> STVector s Int
__triangles TriangulationST s
tglMut))
    Vector Int
halfedges <- MVector (PrimState (ST s)) Int -> ST s (Vector Int)
forall (m :: * -> *) a.
PrimMonad m =>
MVector (PrimState m) a -> m (Vector a)
V.unsafeFreeze (Int -> MVector s Int -> MVector s Int
forall s a. Int -> MVector s a -> MVector s a
VM.take Int
trianglesLen (TriangulationST s -> MVector s Int
forall s. TriangulationST s -> STVector s Int
__halfedges TriangulationST s
tglMut))

    Vector Int
hull <- HullST s -> ST s (Vector Int)
forall s. HullST s -> ST s (Vector Int)
freezeHull HullST s
hullMut
    TriangulationRaw -> ST s TriangulationRaw
forall a. a -> ST s a
forall (f :: * -> *) a. Applicative f => a -> f a
pure TriangulationRaw
        { _triangles :: Vector Int
_triangles = Vector Int
triangles
        , _halfedges :: Vector Int
_halfedges = Vector Int
halfedges
        , _convexHull :: Vector Int
_convexHull = Vector Int
hull
        }

-- | Main result from the raw internal module: a frozen-but-raw triangulation, to
-- be boxed nicely as vectors of polygons etc.
triangulate :: HasCallStack => Vector Vec2 -> TriangulationRaw
triangulate :: HasCallStack => Vector Vec2 -> TriangulationRaw
triangulate Vector Vec2
points = (forall s. ST s TriangulationRaw) -> TriangulationRaw
forall a. (forall s. ST s a) -> a
runST ((forall s. ST s TriangulationRaw) -> TriangulationRaw)
-> (forall s. ST s TriangulationRaw) -> TriangulationRaw
forall a b. (a -> b) -> a -> b
$ do
    (TriangulationST s
tglMut, HullST s
hullMut) <- Vector Vec2 -> ST s (TriangulationST s, HullST s)
forall s.
HasCallStack =>
Vector Vec2 -> ST s (TriangulationST s, HullST s)
triangulateMut Vector Vec2
points
    TriangulationST s -> HullST s -> ST s TriangulationRaw
forall s.
HasCallStack =>
TriangulationST s -> HullST s -> ST s TriangulationRaw
freezeTriangulation TriangulationST s
tglMut HullST s
hullMut