-- | Nice API for Delaunator’s technical output.
module Geometry.Algorithms.Delaunay.Internal.Delaunator.Api where



import           Control.DeepSeq
import           Control.Monad
import           Control.Monad.ST
import qualified Data.Map               as M
import           Data.Vector            (Vector, (!))
import qualified Data.Vector            as V
import qualified Data.Vector.Mutable    as VM
import           Geometry.Core
import           Numerics.Interpolation
import           Util

import qualified Geometry.Algorithms.Delaunay.Internal.Delaunator.Raw as D

import Draw
import Graphics.Rendering.Cairo as C




-- | Abstract data type supporting many efficient Delaunay and Voronoi properties.
data DelaunayTriangulation = Triangulation
    { DelaunayTriangulation -> Vector Polygon
_triangles             :: Vector Polygon
    , DelaunayTriangulation -> Vector Line
_edges                 :: Vector Line
    , DelaunayTriangulation -> Vector Vec2
_voronoiCorners        :: Vector Vec2
    , DelaunayTriangulation -> Vector (Either Line Ray)
_voronoiEdges          :: Vector (Either Line Ray)
    , DelaunayTriangulation -> Vector VoronoiPolygon
_voronoiCells          :: Vector VoronoiPolygon
    , DelaunayTriangulation -> Vector Vec2
_convexHull            :: Vector Vec2
    , DelaunayTriangulation -> Vec2 -> Int -> Int
_findClosestInputPoint :: Vec2 -> Int -> Int
    }

-- | Create a 'DelaunayTriangulation' from a set of points.
delaunayTriangulation :: Sequential vector => vector Vec2 -> DelaunayTriangulation
delaunayTriangulation :: forall (vector :: * -> *).
Sequential vector =>
vector Vec2 -> DelaunayTriangulation
delaunayTriangulation vector Vec2
points' =
    let points :: Vector Vec2
points = vector Vec2 -> Vector Vec2
forall a. vector a -> Vector a
forall (f :: * -> *) a. Sequential f => f a -> Vector a
toVector vector Vec2
points'
        raw :: TriangulationRaw
raw = HasCallStack => Vector Vec2 -> TriangulationRaw
Vector Vec2 -> TriangulationRaw
D.triangulate Vector Vec2
points
        (Vector Polygon
triPolygons, Vector Vec2
triCircumcenters) = Vector Vec2 -> TriangulationRaw -> (Vector Polygon, Vector Vec2)
triangles Vector Vec2
points TriangulationRaw
raw
        extRays :: Vector ExtRays
extRays = Vector Vec2 -> TriangulationRaw -> Vector ExtRays
exteriorRays Vector Vec2
points TriangulationRaw
raw
        inedges :: Map Int Int
inedges = TriangulationRaw -> Map Int Int
bulidInedgesLookup TriangulationRaw
raw
        hullIndex :: Vector Int
hullIndex = Vector Vec2 -> TriangulationRaw -> Vector Int
createHullIndex Vector Vec2
points TriangulationRaw
raw
        findTriangle :: Vec2 -> Int -> Int
findTriangle Vec2
needle Int
i0 = Vector Vec2
-> Map Int Int
-> Vector Int
-> TriangulationRaw
-> Vec2
-> Int
-> Int
findClosestInputPointIndex Vector Vec2
points Map Int Int
inedges Vector Int
hullIndex TriangulationRaw
raw Vec2
needle Int
i0
    in Triangulation
        { _triangles :: Vector Polygon
_triangles = Vector Polygon
triPolygons
        , _edges :: Vector Line
_edges = Vector Vec2 -> TriangulationRaw -> Vector Line
edges Vector Vec2
points TriangulationRaw
raw
        , _findClosestInputPoint :: Vec2 -> Int -> Int
_findClosestInputPoint = Vec2 -> Int -> Int
findTriangle
        , _voronoiCorners :: Vector Vec2
_voronoiCorners = Vector Vec2
triCircumcenters
        , _voronoiEdges :: Vector (Either Line Ray)
_voronoiEdges = Vector Vec2
-> TriangulationRaw -> Vector ExtRays -> Vector (Either Line Ray)
voronoiEdges Vector Vec2
triCircumcenters TriangulationRaw
raw Vector ExtRays
extRays
        , _voronoiCells :: Vector VoronoiPolygon
_voronoiCells = Vector Vec2
-> Vector Vec2
-> Map Int Int
-> TriangulationRaw
-> Vector VoronoiPolygon
voronoiCells Vector Vec2
points Vector Vec2
triCircumcenters Map Int Int
inedges TriangulationRaw
raw
        , _convexHull :: Vector Vec2
_convexHull = Vector Vec2 -> TriangulationRaw -> Vector Vec2
convexHullViaDelaunay Vector Vec2
points TriangulationRaw
raw
        }

instance NFData DelaunayTriangulation where
    rnf :: DelaunayTriangulation -> ()
rnf Triangulation
        { _triangles :: DelaunayTriangulation -> Vector Polygon
_triangles = Vector Polygon
x1
        , _edges :: DelaunayTriangulation -> Vector Line
_edges = Vector Line
x2
        , _findClosestInputPoint :: DelaunayTriangulation -> Vec2 -> Int -> Int
_findClosestInputPoint = Vec2 -> Int -> Int
_
        , _voronoiCorners :: DelaunayTriangulation -> Vector Vec2
_voronoiCorners = Vector Vec2
x3
        , _voronoiEdges :: DelaunayTriangulation -> Vector (Either Line Ray)
_voronoiEdges = Vector (Either Line Ray)
x4
        , _voronoiCells :: DelaunayTriangulation -> Vector VoronoiPolygon
_voronoiCells = Vector VoronoiPolygon
x5
        , _convexHull :: DelaunayTriangulation -> Vector Vec2
_convexHull = Vector Vec2
x6
        } = (Vector Polygon, Vector Line, Vector Vec2,
 Vector (Either Line Ray), Vector VoronoiPolygon, Vector Vec2)
-> ()
forall a. NFData a => a -> ()
rnf (Vector Polygon
x1, Vector Line
x2, Vector Vec2
x3, Vector (Either Line Ray)
x4, Vector VoronoiPolygon
x5, Vector Vec2
x6)

triangles :: Vector Vec2 -> D.TriangulationRaw -> (Vector Polygon, Vector Vec2)
triangles :: Vector Vec2 -> TriangulationRaw -> (Vector Polygon, Vector Vec2)
triangles Vector Vec2
points TriangulationRaw
triangulation =
    let triangleIxs :: Vector Int
triangleIxs = TriangulationRaw -> Vector Int
D._triangles TriangulationRaw
triangulation
        corners :: Vector Vec2
corners = Vector Vec2 -> Vector Int -> Vector Vec2
forall a. Vector a -> Vector Int -> Vector a
V.backpermute Vector Vec2
points Vector Int
triangleIxs
    in ( (Vec2 -> Vec2 -> Vec2 -> Polygon) -> Vector Vec2 -> Vector Polygon
forall a b. (a -> a -> a -> b) -> Vector a -> Vector b
mapChunksOf3 (\Vec2
x Vec2
y Vec2
z -> [Vec2] -> Polygon
Polygon [Vec2
x,Vec2
y,Vec2
z]) Vector Vec2
corners
       , (Vec2 -> Vec2 -> Vec2 -> Vec2) -> Vector Vec2 -> Vector Vec2
forall a b. (a -> a -> a -> b) -> Vector a -> Vector b
mapChunksOf3 (\Vec2
x Vec2
y Vec2
z -> Vec2 -> Vec2 -> Vec2 -> Vec2
D.circumcenter Vec2
x Vec2
y Vec2
z) Vector Vec2
corners)

-- | @mapChunksOf3 f [a,b,c,  i,j,k,  p,q,r] = [f a b c,  f i j k,  f p q r]@
mapChunksOf3 :: (a -> a -> a -> b) -> Vector a -> Vector b
mapChunksOf3 :: forall a b. (a -> a -> a -> b) -> Vector a -> Vector b
mapChunksOf3 a -> a -> a -> b
f Vector a
vec = (forall s. ST s (MVector s b)) -> Vector b
forall a. (forall s. ST s (MVector s a)) -> Vector a
V.create ((forall s. ST s (MVector s b)) -> Vector b)
-> (forall s. ST s (MVector s b)) -> Vector b
forall a b. (a -> b) -> a -> b
$ do
    let len :: Int
len = Vector a -> Int
forall a. Vector a -> Int
V.length Vector a
vec Int -> Int -> Int
forall a. Integral a => a -> a -> a
`div` Int
3
    MVector s b
result <- Int -> ST s (MVector (PrimState (ST s)) b)
forall (m :: * -> *) a.
PrimMonad m =>
Int -> m (MVector (PrimState m) a)
VM.new Int
len
    MVector (PrimState (ST s)) b -> (Int -> b -> ST s ()) -> ST s ()
forall (m :: * -> *) a b.
PrimMonad m =>
MVector (PrimState m) a -> (Int -> a -> m b) -> m ()
VM.iforM_ MVector s b
MVector (PrimState (ST s)) b
result ((Int -> b -> ST s ()) -> ST s ())
-> (Int -> b -> ST s ()) -> ST s ()
forall a b. (a -> b) -> a -> b
$ \Int
i b
_x -> do
        let x :: a
x = Vector a
vec Vector a -> Int -> a
forall a. Vector a -> Int -> a
! (Int
3Int -> Int -> Int
forall a. Num a => a -> a -> a
*Int
iInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
0)
            y :: a
y = Vector a
vec Vector a -> Int -> a
forall a. Vector a -> Int -> a
! (Int
3Int -> Int -> Int
forall a. Num a => a -> a -> a
*Int
iInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1)
            z :: a
z = Vector a
vec Vector a -> Int -> a
forall a. Vector a -> Int -> a
! (Int
3Int -> Int -> Int
forall a. Num a => a -> a -> a
*Int
iInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
2)
        MVector (PrimState (ST s)) b -> Int -> b -> ST s ()
forall (m :: * -> *) a.
PrimMonad m =>
MVector (PrimState m) a -> Int -> a -> m ()
VM.write MVector s b
MVector (PrimState (ST s)) b
result Int
i (a -> a -> a -> b
f a
x a
y a
z)
    MVector s b -> ST s (MVector s b)
forall a. a -> ST s a
forall (f :: * -> *) a. Applicative f => a -> f a
pure MVector s b
result

edges :: Vector Vec2 -> D.TriangulationRaw -> Vector Line
edges :: Vector Vec2 -> TriangulationRaw -> Vector Line
edges Vector Vec2
points TriangulationRaw
triangulation = do
    let triangleIxs :: Vector Int
triangleIxs = TriangulationRaw -> Vector Int
D._triangles TriangulationRaw
triangulation
        halfedges :: Vector Int
halfedges = TriangulationRaw -> Vector Int
D._halfedges TriangulationRaw
triangulation
        numHalfedges :: Int
numHalfedges = Vector Int -> Int
forall a. Vector a -> Int
V.length Vector Int
halfedges
    Int
e <- Int -> Int -> Vector Int
forall a. Num a => a -> Int -> Vector a
V.enumFromN Int
0 Int
numHalfedges

    -- We arbitrarily select the larger of the two edges here. Note that this also
    -- covers the pair-less case, in which the opposite halfedge has index -1.
    Bool -> Vector ()
forall (f :: * -> *). Alternative f => Bool -> f ()
guard (Int
e Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Vector Int
halfedgesVector Int -> Int -> Int
forall a. Vector a -> Int -> a
!Int
e)
    let p :: Vec2
p = Vector Vec2
pointsVector Vec2 -> Int -> Vec2
forall a. Vector a -> Int -> a
!(Vector Int
triangleIxsVector Int -> Int -> Int
forall a. Vector a -> Int -> a
!Int
e)
        q :: Vec2
q = Vector Vec2
pointsVector Vec2 -> Int -> Vec2
forall a. Vector a -> Int -> a
!(Vector Int
triangleIxsVector Int -> Int -> Int
forall a. Vector a -> Int -> a
!Int -> Int
D.nextHalfedge Int
e)
    Line -> Vector Line
forall a. a -> Vector a
forall (f :: * -> *) a. Applicative f => a -> f a
pure (Vec2 -> Vec2 -> Line
Line Vec2
p Vec2
q)

convexHullViaDelaunay :: Vector Vec2 -> D.TriangulationRaw -> Vector Vec2
convexHullViaDelaunay :: Vector Vec2 -> TriangulationRaw -> Vector Vec2
convexHullViaDelaunay Vector Vec2
points TriangulationRaw
triangulation = Vector Vec2 -> Vector Int -> Vector Vec2
forall a. Vector a -> Vector Int -> Vector a
V.backpermute Vector Vec2
points (TriangulationRaw -> Vector Int
D._convexHull TriangulationRaw
triangulation)

-- ^ Given a single edge, what’s the index of the start of the triangle?
triangleOfEdge :: Int -> Int
triangleOfEdge :: Int -> Int
triangleOfEdge Int
e = Int -> Int -> Int
forall a. Integral a => a -> a -> a
div Int
e Int
3

voronoiEdges :: Vector Vec2 -> D.TriangulationRaw -> Vector ExtRays -> Vector (Either Line Ray)
voronoiEdges :: Vector Vec2
-> TriangulationRaw -> Vector ExtRays -> Vector (Either Line Ray)
voronoiEdges Vector Vec2
circumcenters TriangulationRaw
triangulation Vector ExtRays
extRays = do
    let halfedges :: Vector Int
halfedges = TriangulationRaw -> Vector Int
D._halfedges TriangulationRaw
triangulation
        numHalfedges :: Int
numHalfedges = Vector Int -> Int
forall a. Vector a -> Int
V.length Vector Int
halfedges

    Int
e <- Int -> Int -> Vector Int
forall a. Num a => a -> Int -> Vector a
V.enumFromN Int
0 Int
numHalfedges
    let e' :: Int
e' = Vector Int
halfedgesVector Int -> Int -> Int
forall a. Vector a -> Int -> a
!Int
e
        pStart :: Vec2
pStart = Vector Vec2
circumcenters Vector Vec2 -> Int -> Vec2
forall a. Vector a -> Int -> a
! Int -> Int
triangleOfEdge Int
e
    if
        | Int
e' Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
D.tEMPTY ->
            let ExtRays Vec2
_inDir Vec2
outDir = Vector ExtRays
extRays Vector ExtRays -> Int -> ExtRays
forall a. Vector a -> Int -> a
! (TriangulationRaw -> Vector Int
D._triangles TriangulationRaw
triangulation Vector Int -> Int -> Int
forall a. Vector a -> Int -> a
! Int
e)
                -- I don’t know why it’s outDir and not inDir, but I’m quite happy
                -- it’s consistent in my tests. I would have expected more random
                -- behavior. Lucky me!
            in Either Line Ray -> Vector (Either Line Ray)
forall a. a -> Vector a
forall (f :: * -> *) a. Applicative f => a -> f a
pure (Ray -> Either Line Ray
forall a b. b -> Either a b
Right (Vec2 -> Vec2 -> Ray
Ray Vec2
pStart Vec2
outDir))
        | Int
e Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
e' ->
            let pEnd :: Vec2
pEnd = Vector Vec2
circumcenters Vector Vec2 -> Int -> Vec2
forall a. Vector a -> Int -> a
! Int -> Int
triangleOfEdge Int
e'
            in Either Line Ray -> Vector (Either Line Ray)
forall a. a -> Vector a
forall (f :: * -> *) a. Applicative f => a -> f a
pure (Line -> Either Line Ray
forall a b. a -> Either a b
Left (Vec2 -> Vec2 -> Line
Line Vec2
pStart Vec2
pEnd))
        | Bool
otherwise -> Vector (Either Line Ray)
forall a. Monoid a => a
mempty

-- | All edges around a point. The point is specified by an incoming edge.
edgesAroundPoint
    :: D.TriangulationRaw
    -> Int -- ^ Incoming (!) edge to the center point
    -> [Int]
edgesAroundPoint :: TriangulationRaw -> Int -> [Int]
edgesAroundPoint TriangulationRaw
delaunay Int
start = Int -> [Int]
loop Int
start
  where
    loop :: Int -> [Int]
loop Int
incoming =
        let outgoing :: Int
outgoing = Int -> Int
D.nextHalfedge Int
incoming
            incoming' :: Int
incoming' = TriangulationRaw -> Vector Int
D._halfedges TriangulationRaw
delaunay Vector Int -> Int -> Int
forall a. Vector a -> Int -> a
! Int
outgoing
        in if Int
incoming' Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
D.tEMPTY Bool -> Bool -> Bool
&& Int
incoming' Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
start
            then Int
incoming Int -> [Int] -> [Int]
forall a. a -> [a] -> [a]
: Int -> [Int]
loop Int
incoming'
            else [Int
incoming]

-- | A Voronoi Cell can either be an ordinary (finite) polygon,
-- or one that extends to infinity for boundary polygons.
data VoronoiPolygon
    = VoronoiFinite !Polygon -- ^ Ordinary polygon
    | VoronoiInfinite !Vec2 [Vec2] !Vec2
        -- ^ The polygon consists of a list of finite points, and extends to
        -- infinity at the beginning\/end in the direction of the first\/last
        -- argument. For example, the bottom\/right quadrant (in screen coordinates)
        -- would be @'VoronoiInfinite' ('Vec2' 0 1) ['Vec2' 0 0] ('Vec2' 1 0)@.
    deriving (VoronoiPolygon -> VoronoiPolygon -> Bool
(VoronoiPolygon -> VoronoiPolygon -> Bool)
-> (VoronoiPolygon -> VoronoiPolygon -> Bool) -> Eq VoronoiPolygon
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: VoronoiPolygon -> VoronoiPolygon -> Bool
== :: VoronoiPolygon -> VoronoiPolygon -> Bool
$c/= :: VoronoiPolygon -> VoronoiPolygon -> Bool
/= :: VoronoiPolygon -> VoronoiPolygon -> Bool
Eq, Eq VoronoiPolygon
Eq VoronoiPolygon
-> (VoronoiPolygon -> VoronoiPolygon -> Ordering)
-> (VoronoiPolygon -> VoronoiPolygon -> Bool)
-> (VoronoiPolygon -> VoronoiPolygon -> Bool)
-> (VoronoiPolygon -> VoronoiPolygon -> Bool)
-> (VoronoiPolygon -> VoronoiPolygon -> Bool)
-> (VoronoiPolygon -> VoronoiPolygon -> VoronoiPolygon)
-> (VoronoiPolygon -> VoronoiPolygon -> VoronoiPolygon)
-> Ord VoronoiPolygon
VoronoiPolygon -> VoronoiPolygon -> Bool
VoronoiPolygon -> VoronoiPolygon -> Ordering
VoronoiPolygon -> VoronoiPolygon -> VoronoiPolygon
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 :: VoronoiPolygon -> VoronoiPolygon -> Ordering
compare :: VoronoiPolygon -> VoronoiPolygon -> Ordering
$c< :: VoronoiPolygon -> VoronoiPolygon -> Bool
< :: VoronoiPolygon -> VoronoiPolygon -> Bool
$c<= :: VoronoiPolygon -> VoronoiPolygon -> Bool
<= :: VoronoiPolygon -> VoronoiPolygon -> Bool
$c> :: VoronoiPolygon -> VoronoiPolygon -> Bool
> :: VoronoiPolygon -> VoronoiPolygon -> Bool
$c>= :: VoronoiPolygon -> VoronoiPolygon -> Bool
>= :: VoronoiPolygon -> VoronoiPolygon -> Bool
$cmax :: VoronoiPolygon -> VoronoiPolygon -> VoronoiPolygon
max :: VoronoiPolygon -> VoronoiPolygon -> VoronoiPolygon
$cmin :: VoronoiPolygon -> VoronoiPolygon -> VoronoiPolygon
min :: VoronoiPolygon -> VoronoiPolygon -> VoronoiPolygon
Ord, Int -> VoronoiPolygon -> ShowS
[VoronoiPolygon] -> ShowS
VoronoiPolygon -> String
(Int -> VoronoiPolygon -> ShowS)
-> (VoronoiPolygon -> String)
-> ([VoronoiPolygon] -> ShowS)
-> Show VoronoiPolygon
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> VoronoiPolygon -> ShowS
showsPrec :: Int -> VoronoiPolygon -> ShowS
$cshow :: VoronoiPolygon -> String
show :: VoronoiPolygon -> String
$cshowList :: [VoronoiPolygon] -> ShowS
showList :: [VoronoiPolygon] -> ShowS
Show)

instance Sketch VoronoiPolygon where
    sketch :: VoronoiPolygon -> Render ()
sketch (VoronoiFinite Polygon
polygon) = Polygon -> Render ()
forall a. Sketch a => a -> Render ()
sketch Polygon
polygon
    sketch (VoronoiInfinite Vec2
_ [] Vec2
_) = () -> Render ()
forall a. a -> Render a
forall (f :: * -> *) a. Applicative f => a -> f a
pure ()
    sketch (VoronoiInfinite Vec2
dirIn [Vec2]
points Vec2
dirOut) = do
        let sketchRay :: Vec2 -> Vec2 -> Render ()
sketchRay Vec2
start Vec2
dir = Render () -> Render ()
forall a. Render a -> Render a
cairoScope (Render () -> Render ()) -> Render () -> Render ()
forall a b. (a -> b) -> a -> b
$ [Double] -> Double -> Render ()
setDash [Double
3,Double
3] Double
0 Render () -> Render () -> Render ()
forall a b. Render a -> Render b -> Render b
forall (m :: * -> *) a b. Monad m => m a -> m b -> m b
>> Line -> Render ()
forall a. Sketch a => a -> Render ()
sketch ((Double -> Double) -> Line -> Line
resizeLine (Double -> Double -> Double
forall a b. a -> b -> a
const Double
30) (Vec2 -> Vec2 -> Line
Line Vec2
start (Vec2
start Vec2 -> Vec2 -> Vec2
forall v. VectorSpace v => v -> v -> v
+. Vec2
dir)))
        Vec2 -> Vec2 -> Render ()
sketchRay ([Vec2] -> Vec2
forall a. HasCallStack => [a] -> a
head [Vec2]
points) Vec2
dirIn
        Polyline -> Render ()
forall a. Sketch a => a -> Render ()
sketch ([Vec2] -> Polyline
Polyline [Vec2]
points)
        Vec2 -> Vec2 -> Render ()
sketchRay ([Vec2] -> Vec2
forall a. HasCallStack => [a] -> a
head [Vec2]
points) Vec2
dirOut

-- | Construct a single Voronoi polygon.
voronoiPolygon
    :: Vector Vec2 -- ^ Circumcenters
    -> D.TriangulationRaw
    -> Vector ExtRays -- ^ Exterior rays
    -> Int -- ^ Index of the point itself
    -> Int -- ^ Index of an incoming edge towards the point
    -> VoronoiPolygon
voronoiPolygon :: Vector Vec2
-> TriangulationRaw
-> Vector ExtRays
-> Int
-> Int
-> VoronoiPolygon
voronoiPolygon Vector Vec2
circumcenters TriangulationRaw
delaunay Vector ExtRays
extRays Int
p Int
e =
    let cellEdges :: [Int]
cellEdges = TriangulationRaw -> Int -> [Int]
edgesAroundPoint TriangulationRaw
delaunay Int
e
        cellTriangles :: [Int]
cellTriangles = (Int -> Int) -> [Int] -> [Int]
forall a b. (a -> b) -> [a] -> [b]
map Int -> Int
triangleOfEdge [Int]
cellEdges
        vertices :: [Vec2]
vertices = (Int -> Vec2) -> [Int] -> [Vec2]
forall a b. (a -> b) -> [a] -> [b]
map (Vector Vec2
circumcentersVector Vec2 -> Int -> Vec2
forall a. Vector a -> Int -> a
!) [Int]
cellTriangles

    in case Vector ExtRays
extRaysVector ExtRays -> Int -> ExtRays
forall a. Vector a -> Int -> a
!Int
p of
        ExtRays Vec2
dirIn Vec2
dirOut -> Vec2 -> [Vec2] -> Vec2 -> VoronoiPolygon
VoronoiInfinite Vec2
dirIn [Vec2]
vertices Vec2
dirOut
        ExtRays
NoExtRays -> Polygon -> VoronoiPolygon
VoronoiFinite ([Vec2] -> Polygon
Polygon [Vec2]
vertices)

instance NFData VoronoiPolygon where
    rnf :: VoronoiPolygon -> ()
rnf VoronoiFinite{} = ()
    rnf (VoronoiInfinite Vec2
_in [Vec2]
ps Vec2
_out) = [Vec2] -> ()
forall a. NFData a => a -> ()
rnf [Vec2]
ps

-- | Map of point index to an incoming halfedge ID. Originates on the hull for hull
-- points possible, required for reconstructing edge polygons correctly.
bulidInedgesLookup :: D.TriangulationRaw -> M.Map Int Int
bulidInedgesLookup :: TriangulationRaw -> Map Int Int
bulidInedgesLookup TriangulationRaw
delaunay = (Map Int Int -> Int -> Int -> Map Int Int)
-> Map Int Int -> Vector Int -> Map Int Int
forall a b. (a -> Int -> b -> a) -> a -> Vector b -> a
V.ifoldl' Map Int Int -> Int -> Int -> Map Int Int
forall {p}. Map Int Int -> Int -> p -> Map Int Int
addToIndex Map Int Int
forall k a. Map k a
M.empty (TriangulationRaw -> Vector Int
D._triangles TriangulationRaw
delaunay)
  where
    addToIndex :: Map Int Int -> Int -> p -> Map Int Int
addToIndex Map Int Int
acc Int
e p
_t =
        let endpoint :: Int
endpoint = TriangulationRaw -> Vector Int
D._triangles TriangulationRaw
delaunay Vector Int -> Int -> Int
forall a. Vector a -> Int -> a
! Int -> Int
D.nextHalfedge Int
e
            hasSiblingHalfedge :: Bool
hasSiblingHalfedge = TriangulationRaw -> Vector Int
D._halfedges TriangulationRaw
delaunay Vector Int -> Int -> Int
forall a. Vector a -> Int -> a
! Int
e Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
D.tEMPTY
            seen :: Bool
seen = Int -> Map Int Int -> Bool
forall k a. Ord k => k -> Map k a -> Bool
M.member Int
endpoint Map Int Int
acc
        in if Bool -> Bool
not Bool
seen Bool -> Bool -> Bool
|| Bool -> Bool
not Bool
hasSiblingHalfedge
            then Int -> Int -> Map Int Int -> Map Int Int
forall k a. Ord k => k -> a -> Map k a -> Map k a
M.insert Int
endpoint Int
e Map Int Int
acc
            else Map Int Int
acc

voronoiCells :: Vector Vec2 -> Vector Vec2 -> M.Map Int Int -> D.TriangulationRaw -> Vector VoronoiPolygon
voronoiCells :: Vector Vec2
-> Vector Vec2
-> Map Int Int
-> TriangulationRaw
-> Vector VoronoiPolygon
voronoiCells Vector Vec2
points Vector Vec2
circumcenters Map Int Int
inedges TriangulationRaw
delaunay =
    let extRays :: Vector ExtRays
extRays = Vector Vec2 -> TriangulationRaw -> Vector ExtRays
exteriorRays Vector Vec2
points TriangulationRaw
delaunay
    in ((Int -> Vec2 -> VoronoiPolygon)
 -> Vector Vec2 -> Vector VoronoiPolygon)
-> Vector Vec2
-> (Int -> Vec2 -> VoronoiPolygon)
-> Vector VoronoiPolygon
forall a b c. (a -> b -> c) -> b -> a -> c
flip (Int -> Vec2 -> VoronoiPolygon)
-> Vector Vec2 -> Vector VoronoiPolygon
forall a b. (Int -> a -> b) -> Vector a -> Vector b
V.imap Vector Vec2
points ((Int -> Vec2 -> VoronoiPolygon) -> Vector VoronoiPolygon)
-> (Int -> Vec2 -> VoronoiPolygon) -> Vector VoronoiPolygon
forall a b. (a -> b) -> a -> b
$ \Int
pIx Vec2
_pCoord ->
        let incoming :: Int
incoming = Map Int Int
inedges Map Int Int -> Int -> Int
forall k a. Ord k => Map k a -> k -> a
M.! Int
pIx
        in Vector Vec2
-> TriangulationRaw
-> Vector ExtRays
-> Int
-> Int
-> VoronoiPolygon
voronoiPolygon Vector Vec2
circumcenters TriangulationRaw
delaunay Vector ExtRays
extRays Int
pIx Int
incoming

data ExtRays = NoExtRays | ExtRays !Vec2 !Vec2
    deriving (ExtRays -> ExtRays -> Bool
(ExtRays -> ExtRays -> Bool)
-> (ExtRays -> ExtRays -> Bool) -> Eq ExtRays
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: ExtRays -> ExtRays -> Bool
== :: ExtRays -> ExtRays -> Bool
$c/= :: ExtRays -> ExtRays -> Bool
/= :: ExtRays -> ExtRays -> Bool
Eq, Eq ExtRays
Eq ExtRays
-> (ExtRays -> ExtRays -> Ordering)
-> (ExtRays -> ExtRays -> Bool)
-> (ExtRays -> ExtRays -> Bool)
-> (ExtRays -> ExtRays -> Bool)
-> (ExtRays -> ExtRays -> Bool)
-> (ExtRays -> ExtRays -> ExtRays)
-> (ExtRays -> ExtRays -> ExtRays)
-> Ord ExtRays
ExtRays -> ExtRays -> Bool
ExtRays -> ExtRays -> Ordering
ExtRays -> ExtRays -> ExtRays
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 :: ExtRays -> ExtRays -> Ordering
compare :: ExtRays -> ExtRays -> Ordering
$c< :: ExtRays -> ExtRays -> Bool
< :: ExtRays -> ExtRays -> Bool
$c<= :: ExtRays -> ExtRays -> Bool
<= :: ExtRays -> ExtRays -> Bool
$c> :: ExtRays -> ExtRays -> Bool
> :: ExtRays -> ExtRays -> Bool
$c>= :: ExtRays -> ExtRays -> Bool
>= :: ExtRays -> ExtRays -> Bool
$cmax :: ExtRays -> ExtRays -> ExtRays
max :: ExtRays -> ExtRays -> ExtRays
$cmin :: ExtRays -> ExtRays -> ExtRays
min :: ExtRays -> ExtRays -> ExtRays
Ord, Int -> ExtRays -> ShowS
[ExtRays] -> ShowS
ExtRays -> String
(Int -> ExtRays -> ShowS)
-> (ExtRays -> String) -> ([ExtRays] -> ShowS) -> Show ExtRays
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> ExtRays -> ShowS
showsPrec :: Int -> ExtRays -> ShowS
$cshow :: ExtRays -> String
show :: ExtRays -> String
$cshowList :: [ExtRays] -> ShowS
showList :: [ExtRays] -> ShowS
Show)

instance NFData ExtRays where
    rnf :: ExtRays -> ()
rnf ExtRays
NoExtRays = ()
    rnf ExtRays{} = ()

-- | Each point on the Delaunay hull defines two rays:
--     1. The incoming edge (in hull traversal order), rotated by 90° outwards
--     2. The outgoing edge, rotated 90° outwards
--
-- We traverse the hull in order, and create a vector mapping point indices
-- to the rays originating from the incoming/outgoing edge.
--
-- The result vector has the structure /point -> (incoming, outgoing)/.
exteriorRays :: Vector Vec2 -> D.TriangulationRaw -> Vector ExtRays
exteriorRays :: Vector Vec2 -> TriangulationRaw -> Vector ExtRays
exteriorRays Vector Vec2
points TriangulationRaw
delaunay = (forall s. ST s (Vector ExtRays)) -> Vector ExtRays
forall a. (forall s. ST s a) -> a
runST ((forall s. ST s (Vector ExtRays)) -> Vector ExtRays)
-> (forall s. ST s (Vector ExtRays)) -> Vector ExtRays
forall a b. (a -> b) -> a -> b
$ do
    let hull :: Vector Int
hull = TriangulationRaw -> Vector Int
D._convexHull TriangulationRaw
delaunay
    MVector (PrimState (ST s)) (Maybe Vec2)
inRays <- Int -> Maybe Vec2 -> ST s (MVector (PrimState (ST s)) (Maybe Vec2))
forall (m :: * -> *) a.
PrimMonad m =>
Int -> a -> m (MVector (PrimState m) a)
VM.replicate (Vector Vec2 -> Int
forall a. Vector a -> Int
V.length Vector Vec2
points) Maybe Vec2
forall a. Maybe a
Nothing
    MVector (PrimState (ST s)) (Maybe Vec2)
outRays <- Int -> Maybe Vec2 -> ST s (MVector (PrimState (ST s)) (Maybe Vec2))
forall (m :: * -> *) a.
PrimMonad m =>
Int -> a -> m (MVector (PrimState m) a)
VM.replicate (Vector Vec2 -> Int
forall a. Vector a -> Int
V.length Vector Vec2
points) Maybe Vec2
forall a. Maybe a
Nothing
    let recordRays :: Int -> Int -> ST s ()
recordRays = \Int
pStart Int
pEnd -> do
            let vecStart :: Vec2
vecStart = Vector Vec2
pointsVector Vec2 -> Int -> Vec2
forall a. Vector a -> Int -> a
!Int
pStart
                vecEnd :: Vec2
vecEnd = Vector Vec2
pointsVector Vec2 -> Int -> Vec2
forall a. Vector a -> Int -> a
!Int
pEnd

                rayDir :: Vec2
rayDir = Vec2 -> Vec2
rotate90 (Vec2
vecEnd Vec2 -> Vec2 -> Vec2
forall v. VectorSpace v => v -> v -> v
-. Vec2
vecStart)

            -- Record as outgoing ray for pStart
            MVector (PrimState (ST s)) (Maybe Vec2)
-> Int -> Maybe Vec2 -> ST s ()
forall (m :: * -> *) a.
PrimMonad m =>
MVector (PrimState m) a -> Int -> a -> m ()
VM.write MVector (PrimState (ST s)) (Maybe Vec2)
outRays Int
pStart (Vec2 -> Maybe Vec2
forall a. a -> Maybe a
Just Vec2
rayDir)

            -- Record as incoming ray for pEnd
            MVector (PrimState (ST s)) (Maybe Vec2)
-> Int -> Maybe Vec2 -> ST s ()
forall (m :: * -> *) a.
PrimMonad m =>
MVector (PrimState m) a -> Int -> a -> m ()
VM.write MVector (PrimState (ST s)) (Maybe Vec2)
inRays Int
pEnd (Vec2 -> Maybe Vec2
forall a. a -> Maybe a
Just Vec2
rayDir)

    (Int -> Int -> ST s ()) -> Vector Int -> Vector Int -> ST s ()
forall (m :: * -> *) a b c.
Monad m =>
(a -> b -> m c) -> Vector a -> Vector b -> m ()
V.zipWithM_ Int -> Int -> ST s ()
recordRays Vector Int
hull (Vector Int -> Vector Int
forall a. Vector a -> Vector a
V.tail Vector Int
hull)
    ()
_ <- Int -> Int -> ST s ()
recordRays (Vector Int -> Int
forall a. Vector a -> a
V.last Vector Int
hull) (Vector Int -> Int
forall a. Vector a -> a
V.head Vector Int
hull) -- zip omits the cyclic pair, so we do it manually

    (Vector (Maybe Vec2)
a, Vector (Maybe Vec2)
b) <- (,) (Vector (Maybe Vec2)
 -> Vector (Maybe Vec2)
 -> (Vector (Maybe Vec2), Vector (Maybe Vec2)))
-> ST s (Vector (Maybe Vec2))
-> ST
     s
     (Vector (Maybe Vec2) -> (Vector (Maybe Vec2), Vector (Maybe Vec2)))
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> MVector (PrimState (ST s)) (Maybe Vec2)
-> ST s (Vector (Maybe Vec2))
forall (m :: * -> *) a.
PrimMonad m =>
MVector (PrimState m) a -> m (Vector a)
V.unsafeFreeze MVector (PrimState (ST s)) (Maybe Vec2)
inRays ST
  s
  (Vector (Maybe Vec2) -> (Vector (Maybe Vec2), Vector (Maybe Vec2)))
-> ST s (Vector (Maybe Vec2))
-> ST s (Vector (Maybe Vec2), Vector (Maybe Vec2))
forall a b. ST s (a -> b) -> ST s a -> ST s b
forall (f :: * -> *) a b. Applicative f => f (a -> b) -> f a -> f b
<*> MVector (PrimState (ST s)) (Maybe Vec2)
-> ST s (Vector (Maybe Vec2))
forall (m :: * -> *) a.
PrimMonad m =>
MVector (PrimState m) a -> m (Vector a)
V.unsafeFreeze MVector (PrimState (ST s)) (Maybe Vec2)
outRays
    Vector ExtRays -> ST s (Vector ExtRays)
forall a. a -> ST s a
forall (f :: * -> *) a. Applicative f => a -> f a
pure ((Maybe Vec2 -> Maybe Vec2 -> ExtRays)
-> Vector (Maybe Vec2) -> Vector (Maybe Vec2) -> Vector ExtRays
forall a b c. (a -> b -> c) -> Vector a -> Vector b -> Vector c
V.zipWith
        (\Maybe Vec2
x Maybe Vec2
y -> case (Maybe Vec2
x,Maybe Vec2
y) of
            (Just Vec2
inDir, Just Vec2
outDir) -> Vec2 -> Vec2 -> ExtRays
ExtRays Vec2
inDir Vec2
outDir
            (Maybe Vec2
Nothing, Maybe Vec2
Nothing) -> ExtRays
NoExtRays
            (Maybe Vec2, Maybe Vec2)
other -> String -> String -> ExtRays
forall a. String -> String -> a
bugError String
"exteriorRays" (String
"Bad external ray pair: " String -> ShowS
forall a. [a] -> [a] -> [a]
++ (Maybe Vec2, Maybe Vec2) -> String
forall a. Show a => a -> String
show (Maybe Vec2, Maybe Vec2)
other)
        )
        Vector (Maybe Vec2)
a
        Vector (Maybe Vec2)
b)

-- | Rotate a 'Vec2' by 90°
rotate90 :: Vec2 -> Vec2
rotate90 :: Vec2 -> Vec2
rotate90 (Vec2 Double
x Double
y) = Double -> Double -> Vec2
Vec2 (-Double
y) Double
x

-- | A ray is a line that extends to infinity on one side.
data Ray = Ray !Vec2 !Vec2 -- ^ Starting point and direction (!)
    deriving (Ray -> Ray -> Bool
(Ray -> Ray -> Bool) -> (Ray -> Ray -> Bool) -> Eq Ray
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: Ray -> Ray -> Bool
== :: Ray -> Ray -> Bool
$c/= :: Ray -> Ray -> Bool
/= :: Ray -> Ray -> Bool
Eq, Eq Ray
Eq Ray
-> (Ray -> Ray -> Ordering)
-> (Ray -> Ray -> Bool)
-> (Ray -> Ray -> Bool)
-> (Ray -> Ray -> Bool)
-> (Ray -> Ray -> Bool)
-> (Ray -> Ray -> Ray)
-> (Ray -> Ray -> Ray)
-> Ord Ray
Ray -> Ray -> Bool
Ray -> Ray -> Ordering
Ray -> Ray -> Ray
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 :: Ray -> Ray -> Ordering
compare :: Ray -> Ray -> Ordering
$c< :: Ray -> Ray -> Bool
< :: Ray -> Ray -> Bool
$c<= :: Ray -> Ray -> Bool
<= :: Ray -> Ray -> Bool
$c> :: Ray -> Ray -> Bool
> :: Ray -> Ray -> Bool
$c>= :: Ray -> Ray -> Bool
>= :: Ray -> Ray -> Bool
$cmax :: Ray -> Ray -> Ray
max :: Ray -> Ray -> Ray
$cmin :: Ray -> Ray -> Ray
min :: Ray -> Ray -> Ray
Ord, Int -> Ray -> ShowS
[Ray] -> ShowS
Ray -> String
(Int -> Ray -> ShowS)
-> (Ray -> String) -> ([Ray] -> ShowS) -> Show Ray
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> Ray -> ShowS
showsPrec :: Int -> Ray -> ShowS
$cshow :: Ray -> String
show :: Ray -> String
$cshowList :: [Ray] -> ShowS
showList :: [Ray] -> ShowS
Show)

instance NFData Ray where
    rnf :: Ray -> ()
rnf Ray
_ = () -- Already strict

instance Sketch Ray where
    sketch :: Ray -> Render ()
sketch (Ray Vec2
start Vec2
dir) = Line -> Render ()
forall a. Sketch a => a -> Render ()
sketch ((Double -> Double) -> Line -> Line
resizeLine (Double -> Double -> Double
forall a b. a -> b -> a
const Double
100000) (Vec2 -> Vec2 -> Line
Line Vec2
start (Vec2
start Vec2 -> Vec2 -> Vec2
forall v. VectorSpace v => v -> v -> v
+. Vec2
dir)))

-- | Reverse lookup table for the hull. @'D._convexHull'!i@ yields the i-th point’s
-- ID on the hull, this index gives us the number i, given a point ID.
createHullIndex :: Vector Vec2 -> D.TriangulationRaw -> Vector Int
createHullIndex :: Vector Vec2 -> TriangulationRaw -> Vector Int
createHullIndex Vector Vec2
points TriangulationRaw
raw = (forall s. ST s (MVector s Int)) -> Vector Int
forall a. (forall s. ST s (MVector s a)) -> Vector a
V.create ((forall s. ST s (MVector s Int)) -> Vector Int)
-> (forall s. ST s (MVector s Int)) -> Vector Int
forall a b. (a -> b) -> a -> b
$ do
    let hull :: Vector Int
hull = TriangulationRaw -> Vector Int
D._convexHull TriangulationRaw
raw
    MVector s Int
hullIndex <- Int -> Int -> ST s (MVector (PrimState (ST s)) Int)
forall (m :: * -> *) a.
PrimMonad m =>
Int -> a -> m (MVector (PrimState m) a)
VM.replicate (Vector Vec2 -> Int
forall a. Vector a -> Int
V.length Vector Vec2
points) (-Int
1)
    Vector Int -> (Int -> Int -> ST s ()) -> ST s ()
forall (m :: * -> *) a b.
Monad m =>
Vector a -> (Int -> a -> m b) -> m ()
V.iforM_ Vector Int
hull ((Int -> Int -> ST s ()) -> ST s ())
-> (Int -> Int -> ST s ()) -> ST s ()
forall a b. (a -> b) -> a -> b
$ \Int
i Int
hull_i -> MVector (PrimState (ST s)) Int -> Int -> Int -> ST s ()
forall (m :: * -> *) a.
PrimMonad m =>
MVector (PrimState m) a -> Int -> a -> m ()
VM.write MVector s Int
MVector (PrimState (ST s)) Int
hullIndex Int
hull_i Int
i
    MVector s Int -> ST s (MVector s Int)
forall a. a -> ST s a
forall (f :: * -> *) a. Applicative f => a -> f a
pure MVector s Int
hullIndex

-- | Find the input point closest to the needle. Search starts at specified point i.
findClosestInputPointIndex
    :: Vector Vec2   -- ^ Input points
    -> M.Map Int Int -- ^ Incoming edges table
    -> Vector Int    -- ^ Hull index, see 'createHullIndex'
    -> D.TriangulationRaw
    -> Vec2          -- ^ Needle: which input point is closest to this?
    -> Int           -- ^ Start the search at this index. 0 searches from the beginning.
    -> Int           -- ^ Index of the closest point
findClosestInputPointIndex :: Vector Vec2
-> Map Int Int
-> Vector Int
-> TriangulationRaw
-> Vec2
-> Int
-> Int
findClosestInputPointIndex Vector Vec2
points Map Int Int
inedges Vector Int
hullIndex TriangulationRaw
tri Vec2
needle Int
i0 = Int -> Int
loopFind Int
i0
  where
    loopFind :: Int -> Int
loopFind Int
i =
        let c :: Int
c = Int -> Int
step Int
i
        in if Int
c Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
0 Bool -> Bool -> Bool
&& Int
c Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
i Bool -> Bool -> Bool
&& Int
c Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
i0
            then Int -> Int
loopFind Int
c
            else Int
c

    -- The idea of one step is to look at outgoing edges of a point, and following
    -- the one that leads us closer to the needle.
    step :: Int -> Int
step Int
j =
        let c :: Int
c = Int
j
            dc :: Double
dc = Vec2 -> Double
normSquare (Vec2
needle Vec2 -> Vec2 -> Vec2
forall v. VectorSpace v => v -> v -> v
-. Vector Vec2
pointsVector Vec2 -> Int -> Vec2
forall a. Vector a -> Int -> a
!Int
j)
            e0 :: Int
e0 = Map Int Int
inedges Map Int Int -> Int -> Int
forall k a. Ord k => Map k a -> k -> a
M.! Int
j
            e :: Int
e = Int
e0
        in Int -> Int -> Double -> Int -> Int -> Int
loopStep Int
j Int
c Double
dc Int
e0 Int
e

    loopStep
        :: Int
        -> Int    -- c:  Start of search (candidate)
        -> Double -- dc: Distance² from candidate to needle
        -> Int    -- e0: inedge the search has started
        -> Int    -- e:  inedge we’re currently searching
        -> Int    -- Better candidate after the step
    loopStep :: Int -> Int -> Double -> Int -> Int -> Int
loopStep Int
j Int
c Double
dc Int
e0 Int
e =
        let t :: Int
t = TriangulationRaw -> Vector Int
D._triangles TriangulationRaw
tri Vector Int -> Int -> Int
forall a. Vector a -> Int -> a
! Int
e
            dt :: Double
dt = Vec2 -> Double
normSquare (Vec2
needle Vec2 -> Vec2 -> Vec2
forall v. VectorSpace v => v -> v -> v
-. Vector Vec2
pointsVector Vec2 -> Int -> Vec2
forall a. Vector a -> Int -> a
!Int
t)
            (Double
dc', Int
c') | Double
dt Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
dc   = (Double
dt, Int
t)
                      | Bool
otherwise = (Double
dc, Int
c)
            e' :: Int
e' = TriangulationRaw -> Vector Int
D._halfedges TriangulationRaw
tri Vector Int -> Int -> Int
forall a. Vector a -> Int -> a
! Int -> Int
D.nextHalfedge Int
e
        in if Int
e' Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
D.tEMPTY
            then -- The next edge has no partner: we’re on the hull
                let e'' :: Int
e'' = TriangulationRaw -> Vector Int
D._convexHull TriangulationRaw
tri Vector Int -> Int -> Int
forall a. Vector a -> Int -> a
! (((Vector Int
hullIndexVector Int -> Int -> Int
forall a. Vector a -> Int -> a
!Int
j) Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) Int -> Int -> Int
forall a. Integral a => a -> a -> a
`mod` Vector Int -> Int
forall a. Vector a -> Int
V.length (TriangulationRaw -> Vector Int
D._convexHull TriangulationRaw
tri))
                in if Int
e'' Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
t Bool -> Bool -> Bool
&& Vec2 -> Double
normSquare (Vec2
needle Vec2 -> Vec2 -> Vec2
forall v. VectorSpace v => v -> v -> v
-. Vector Vec2
pointsVector Vec2 -> Int -> Vec2
forall a. Vector a -> Int -> a
!Int
e'') Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
dc'
                    then Int
e''
                    else Int
c'
            else -- We’re not on the hull
                if Int
e' Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
e0
                    then Int -> Int -> Double -> Int -> Int -> Int
loopStep Int
j Int
c' Double
dc' Int
e0 Int
e'
                    else Int
c'

-- | Voronoi stippling. Unfortunately not fast enough for stippling images
-- pixel-wise, hence not exposed from the user-facing module.
stipple
    :: Double -- ^ \(\omega\) convergence speed parameter, see 'lloydRelaxation'.
    -> Int -- ^ Width of the input data. x values will be picked in the integer range \([0\ldots w)\).
    -> Int -- ^ Height of the input data. x values will be picked in the integer range \([0\ldots h)\).
    -> (Int -> Int -> Double) -- ^ How much weight does the input have at \(f(x,y)\)?
    -> Vector Vec2 -- ^ Input points, chosen as last parameter for 'iterate' convenience.
    -> Int -- ^ Number of generations
    -> Vector Vec2
stipple :: Double
-> Int
-> Int
-> (Int -> Int -> Double)
-> Vector Vec2
-> Int
-> Vector Vec2
stipple Double
omega Int
width Int
height Int -> Int -> Double
f Vector Vec2
points Int
n = (Vector Vec2 -> Vector Vec2) -> Vector Vec2 -> [Vector Vec2]
forall a. (a -> a) -> a -> [a]
iterate (Double
-> Int
-> Int
-> (Int -> Int -> Double)
-> Vector Vec2
-> Vector Vec2
stippleStep Double
omega Int
width Int
height Int -> Int -> Double
f) Vector Vec2
points [Vector Vec2] -> Int -> Vector Vec2
forall a. HasCallStack => [a] -> Int -> a
!! Int
n

-- | One step towards Voronoi stippling. See 'stipple' for argument docs.
stippleStep
    :: Double
    -> Int
    -> Int
    -> (Int -> Int -> Double)
    -> Vector Vec2
    -> Vector Vec2
stippleStep :: Double
-> Int
-> Int
-> (Int -> Int -> Double)
-> Vector Vec2
-> Vector Vec2
stippleStep Double
omega Int
width Int
height Int -> Int -> Double
f Vector Vec2
points = (forall s. ST s (Vector Vec2)) -> Vector Vec2
forall a. (forall s. ST s a) -> a
runST ((forall s. ST s (Vector Vec2)) -> Vector Vec2)
-> (forall s. ST s (Vector Vec2)) -> Vector Vec2
forall a b. (a -> b) -> a -> b
$ do
    let numPoints :: Int
numPoints = Vector Vec2 -> Int
forall a. Vector a -> Int
V.length Vector Vec2
points
        tri :: DelaunayTriangulation
tri = Vector Vec2 -> DelaunayTriangulation
forall (vector :: * -> *).
Sequential vector =>
vector Vec2 -> DelaunayTriangulation
delaunayTriangulation Vector Vec2
points

    MVector (PrimState (ST s)) Vec2
centroidsMut <- Int -> Vec2 -> ST s (MVector (PrimState (ST s)) Vec2)
forall (m :: * -> *) a.
PrimMonad m =>
Int -> a -> m (MVector (PrimState m) a)
VM.replicate Int
numPoints Vec2
forall v. VectorSpace v => v
zero
    MVector (PrimState (ST s)) Double
weightsMut <- Int -> Double -> ST s (MVector (PrimState (ST s)) Double)
forall (m :: * -> *) a.
PrimMonad m =>
Int -> a -> m (MVector (PrimState m) a)
VM.replicate Int
numPoints Double
0
    let loopY :: Int -> Int -> ST s ()
loopY Int
_ Int
y | Int
y Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
height = () -> ST s ()
forall a. a -> ST s a
forall (f :: * -> *) a. Applicative f => a -> f a
pure ()
        loopY Int
i Int
y = Int -> Int -> Int -> ST s ()
loopX Int
i Int
0 Int
y

        loopX :: Int -> Int -> Int -> ST s ()
loopX Int
i Int
x Int
y | Int
x Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
width = Int -> Int -> ST s ()
loopY Int
i (Int
yInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1)
        loopX Int
i Int
x Int
y = do
            let w :: Double
w = Int -> Int -> Double
f Int
x Int
y
                pixelTopLeft :: Vec2
pixelTopLeft = Double -> Double -> Vec2
Vec2 (Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
x) (Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
y)
                pixelCenter :: Vec2
pixelCenter = Vec2
pixelTopLeft Vec2 -> Vec2 -> Vec2
forall v. VectorSpace v => v -> v -> v
+. Double -> Double -> Vec2
Vec2 Double
0.5 Double
0.5
                i' :: Int
i' = DelaunayTriangulation -> Vec2 -> Int -> Int
_findClosestInputPoint DelaunayTriangulation
tri Vec2
pixelTopLeft Int
i
            MVector (PrimState (ST s)) Double
-> (Double -> Double) -> Int -> ST s ()
forall (m :: * -> *) a.
PrimMonad m =>
MVector (PrimState m) a -> (a -> a) -> Int -> m ()
VM.modify MVector (PrimState (ST s)) Double
weightsMut (Double -> Double -> Double
forall a. Num a => a -> a -> a
+Double
w) Int
i'
            MVector (PrimState (ST s)) Vec2 -> (Vec2 -> Vec2) -> Int -> ST s ()
forall (m :: * -> *) a.
PrimMonad m =>
MVector (PrimState m) a -> (a -> a) -> Int -> m ()
VM.modify MVector (PrimState (ST s)) Vec2
centroidsMut (Vec2 -> Vec2 -> Vec2
forall v. VectorSpace v => v -> v -> v
+. (Double
w Double -> Vec2 -> Vec2
forall v. VectorSpace v => Double -> v -> v
*. Vec2
pixelCenter)) Int
i'
            Int -> Int -> Int -> ST s ()
loopX Int
i' (Int
xInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1) Int
y
    Int -> Int -> ST s ()
loopY Int
0 Int
0

    Vector Vec2
centroids <- MVector (PrimState (ST s)) Vec2 -> ST s (Vector Vec2)
forall (m :: * -> *) a.
PrimMonad m =>
MVector (PrimState m) a -> m (Vector a)
V.unsafeFreeze MVector (PrimState (ST s)) Vec2
centroidsMut
    Vector Double
weights <- MVector (PrimState (ST s)) Double -> ST s (Vector Double)
forall (m :: * -> *) a.
PrimMonad m =>
MVector (PrimState m) a -> m (Vector a)
V.unsafeFreeze MVector (PrimState (ST s)) Double
weightsMut

    Vector Vec2 -> ST s (Vector Vec2)
forall a. a -> ST s a
forall (f :: * -> *) a. Applicative f => a -> f a
pure (Vector Vec2 -> ST s (Vector Vec2))
-> Vector Vec2 -> ST s (Vector Vec2)
forall a b. (a -> b) -> a -> b
$ (Vec2 -> Double -> Vec2 -> Vec2)
-> Vector Vec2 -> Vector Double -> Vector Vec2 -> Vector Vec2
forall a b c d.
(a -> b -> c -> d) -> Vector a -> Vector b -> Vector c -> Vector d
V.zipWith3
        (\Vec2
c Double
w Vec2
p ->
            let p' :: Vec2
p' | Double
w Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0 = Vec2
c Vec2 -> Double -> Vec2
forall v. VectorSpace v => v -> Double -> v
/. Double
w
                   | Bool
otherwise = Vec2
p
            in (Double, Double) -> (Vec2, Vec2) -> Double -> Vec2
forall vec.
VectorSpace vec =>
(Double, Double) -> (vec, vec) -> Double -> vec
lerp (Double
0,Double
1) (Vec2
p,Vec2
p') Double
omega
        )
        Vector Vec2
centroids
        Vector Double
weights
        Vector Vec2
points