--------------------------------------------------------------------------------
-- |
-- Module      :  Algorithms.Geometry.DelaunayTriangulation.Naive
-- Copyright   :  (C) Frank Staals
-- License     :  see the LICENSE file
-- Maintainer  :  Frank Staals
--------------------------------------------------------------------------------
module Algorithms.Geometry.DelaunayTriangulation.Naive where

import           Algorithms.Geometry.DelaunayTriangulation.Types
import           Control.Lens
import           Control.Monad (forM_)
import qualified Data.CircularList as C
import           Data.Ext
import qualified Data.Foldable as F
import           Data.Function (on)
import           Data.Geometry
import           Data.Geometry.Ball (disk, insideBall)
import qualified Data.List as L
import qualified Data.List.NonEmpty as NonEmpty
import qualified Data.Map as M
import qualified Data.Vector as V
import qualified Data.Vector.Mutable as MV

--------------------------------------------------------------------------------

-- | Naive \( O(n^4) \) time implementation of the delaunay triangulation. Simply
-- tries each triple (p,q,r) and tests if it is delaunay, i.e. if there are no
-- other points in the circle defined by p, q, and r.
--
-- pre: the input is a *SET*, i.e. contains no duplicate points. (If the
-- input does contain duplicate points, the implementation throws them away)
delaunayTriangulation     :: (Ord r, Fractional r)
                          => NonEmpty.NonEmpty (Point 2 r :+ p) -> Triangulation p r
delaunayTriangulation :: NonEmpty (Point 2 r :+ p) -> Triangulation p r
delaunayTriangulation NonEmpty (Point 2 r :+ p)
pts = Map (Point 2 r) VertexID
-> Vector (Point 2 r :+ p)
-> Vector (CList VertexID)
-> Triangulation p r
forall p r.
Map (Point 2 r) VertexID
-> Vector (Point 2 r :+ p)
-> Vector (CList VertexID)
-> Triangulation p r
Triangulation Map (Point 2 r) VertexID
ptIds Vector (Point 2 r :+ p)
ptsV Vector (CList VertexID)
adjV
  where
    ptsV :: Vector (Point 2 r :+ p)
ptsV   = [Point 2 r :+ p] -> Vector (Point 2 r :+ p)
forall a. [a] -> Vector a
V.fromList ([Point 2 r :+ p] -> Vector (Point 2 r :+ p))
-> (NonEmpty (Point 2 r :+ p) -> [Point 2 r :+ p])
-> NonEmpty (Point 2 r :+ p)
-> Vector (Point 2 r :+ p)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. NonEmpty (Point 2 r :+ p) -> [Point 2 r :+ p]
forall (t :: * -> *) a. Foldable t => t a -> [a]
F.toList (NonEmpty (Point 2 r :+ p) -> [Point 2 r :+ p])
-> (NonEmpty (Point 2 r :+ p) -> NonEmpty (Point 2 r :+ p))
-> NonEmpty (Point 2 r :+ p)
-> [Point 2 r :+ p]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. ((Point 2 r :+ p) -> (Point 2 r :+ p) -> Bool)
-> NonEmpty (Point 2 r :+ p) -> NonEmpty (Point 2 r :+ p)
forall a. (a -> a -> Bool) -> NonEmpty a -> NonEmpty a
NonEmpty.nubBy (Point 2 r -> Point 2 r -> Bool
forall a. Eq a => a -> a -> Bool
(==) (Point 2 r -> Point 2 r -> Bool)
-> ((Point 2 r :+ p) -> Point 2 r)
-> (Point 2 r :+ p)
-> (Point 2 r :+ p)
-> Bool
forall b c a. (b -> b -> c) -> (a -> b) -> a -> a -> c
`on` ((Point 2 r :+ p)
-> Getting (Point 2 r) (Point 2 r :+ p) (Point 2 r) -> Point 2 r
forall s a. s -> Getting a s a -> a
^.Getting (Point 2 r) (Point 2 r :+ p) (Point 2 r)
forall core extra core'.
Lens (core :+ extra) (core' :+ extra) core core'
core)) (NonEmpty (Point 2 r :+ p) -> Vector (Point 2 r :+ p))
-> NonEmpty (Point 2 r :+ p) -> Vector (Point 2 r :+ p)
forall a b. (a -> b) -> a -> b
$ NonEmpty (Point 2 r :+ p)
pts
    ptIds :: Map (Point 2 r) VertexID
ptIds  = [(Point 2 r, VertexID)] -> Map (Point 2 r) VertexID
forall k a. Ord k => [(k, a)] -> Map k a
M.fromList ([(Point 2 r, VertexID)] -> Map (Point 2 r) VertexID)
-> [(Point 2 r, VertexID)] -> Map (Point 2 r) VertexID
forall a b. (a -> b) -> a -> b
$ [Point 2 r] -> [VertexID] -> [(Point 2 r, VertexID)]
forall a b. [a] -> [b] -> [(a, b)]
zip (((Point 2 r :+ p) -> Point 2 r) -> [Point 2 r :+ p] -> [Point 2 r]
forall a b. (a -> b) -> [a] -> [b]
map ((Point 2 r :+ p)
-> Getting (Point 2 r) (Point 2 r :+ p) (Point 2 r) -> Point 2 r
forall s a. s -> Getting a s a -> a
^.Getting (Point 2 r) (Point 2 r :+ p) (Point 2 r)
forall core extra core'.
Lens (core :+ extra) (core' :+ extra) core core'
core) ([Point 2 r :+ p] -> [Point 2 r])
-> (Vector (Point 2 r :+ p) -> [Point 2 r :+ p])
-> Vector (Point 2 r :+ p)
-> [Point 2 r]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Vector (Point 2 r :+ p) -> [Point 2 r :+ p]
forall a. Vector a -> [a]
V.toList (Vector (Point 2 r :+ p) -> [Point 2 r])
-> Vector (Point 2 r :+ p) -> [Point 2 r]
forall a b. (a -> b) -> a -> b
$ Vector (Point 2 r :+ p)
ptsV) [VertexID
0..]
    adjV :: Vector (CList VertexID)
adjV   = Mapping p r -> [(VertexID, VertexID)] -> Vector (CList VertexID)
forall r p.
(Num r, Ord r) =>
Mapping p r -> [(VertexID, VertexID)] -> Vector (CList VertexID)
toAdjLists (Map (Point 2 r) VertexID
ptIds,Vector (Point 2 r :+ p)
ptsV) ([(VertexID, VertexID)] -> Vector (CList VertexID))
-> ([(VertexID, VertexID, VertexID)] -> [(VertexID, VertexID)])
-> [(VertexID, VertexID, VertexID)]
-> Vector (CList VertexID)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [(VertexID, VertexID, VertexID)] -> [(VertexID, VertexID)]
extractEdges ([(VertexID, VertexID, VertexID)] -> Vector (CList VertexID))
-> [(VertexID, VertexID, VertexID)] -> Vector (CList VertexID)
forall a b. (a -> b) -> a -> b
$ [(VertexID, VertexID, VertexID)]
fs
    n :: VertexID
n      = Vector (Point 2 r :+ p) -> VertexID
forall a. Vector a -> VertexID
V.length Vector (Point 2 r :+ p)
ptsV VertexID -> VertexID -> VertexID
forall a. Num a => a -> a -> a
- VertexID
1

    -- construct the list of faces/triangles in the delaunay triangulation
    fs :: [(VertexID, VertexID, VertexID)]
fs = [ (VertexID
p,VertexID
q,VertexID
r)
         | VertexID
p <- [VertexID
0..VertexID
n], VertexID
q <- [VertexID
p..VertexID
n], VertexID
r <- [VertexID
q..VertexID
n], Mapping p r -> VertexID -> VertexID -> VertexID -> Bool
forall r p.
(Fractional r, Ord r) =>
Mapping p r -> VertexID -> VertexID -> VertexID -> Bool
isDelaunay (Map (Point 2 r) VertexID
ptIds,Vector (Point 2 r :+ p)
ptsV) VertexID
p VertexID
q VertexID
r
         ]

-- | Given a list of edges, as vertexId pairs, construct a vector with the
-- adjacency lists, each in CW sorted order.
toAdjLists             :: (Num r, Ord r) => Mapping p r -> [(VertexID,VertexID)]
                       -> V.Vector (C.CList VertexID)
toAdjLists :: Mapping p r -> [(VertexID, VertexID)] -> Vector (CList VertexID)
toAdjLists m :: Mapping p r
m@(Map (Point 2 r) VertexID
_,Vector (Point 2 r :+ p)
ptsV) [(VertexID, VertexID)]
es = (VertexID -> [VertexID] -> CList VertexID)
-> Vector [VertexID] -> Vector (CList VertexID)
forall a b. (VertexID -> a -> b) -> Vector a -> Vector b
V.imap VertexID -> [VertexID] -> CList VertexID
toCList (Vector [VertexID] -> Vector (CList VertexID))
-> Vector [VertexID] -> Vector (CList VertexID)
forall a b. (a -> b) -> a -> b
$ (forall s. ST s (MVector s [VertexID])) -> Vector [VertexID]
forall a. (forall s. ST s (MVector s a)) -> Vector a
V.create ((forall s. ST s (MVector s [VertexID])) -> Vector [VertexID])
-> (forall s. ST s (MVector s [VertexID])) -> Vector [VertexID]
forall a b. (a -> b) -> a -> b
$ do
    MVector s [VertexID]
v <- VertexID
-> [VertexID] -> ST s (MVector (PrimState (ST s)) [VertexID])
forall (m :: * -> *) a.
PrimMonad m =>
VertexID -> a -> m (MVector (PrimState m) a)
MV.replicate (Vector (Point 2 r :+ p) -> VertexID
forall a. Vector a -> VertexID
V.length Vector (Point 2 r :+ p)
ptsV) []
    [(VertexID, VertexID)]
-> ((VertexID, VertexID) -> ST s ()) -> ST s ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
t a -> (a -> m b) -> m ()
forM_ [(VertexID, VertexID)]
es (((VertexID, VertexID) -> ST s ()) -> ST s ())
-> ((VertexID, VertexID) -> ST s ()) -> ST s ()
forall a b. (a -> b) -> a -> b
$ \(VertexID
i,VertexID
j) -> do
      MVector (PrimState (ST s)) [VertexID]
-> VertexID -> VertexID -> ST s ()
forall (m :: * -> *) a.
PrimMonad m =>
MVector (PrimState m) [a] -> VertexID -> a -> m ()
addAt MVector s [VertexID]
MVector (PrimState (ST s)) [VertexID]
v VertexID
i VertexID
j
      MVector (PrimState (ST s)) [VertexID]
-> VertexID -> VertexID -> ST s ()
forall (m :: * -> *) a.
PrimMonad m =>
MVector (PrimState m) [a] -> VertexID -> a -> m ()
addAt MVector s [VertexID]
MVector (PrimState (ST s)) [VertexID]
v VertexID
j VertexID
i
    MVector s [VertexID] -> ST s (MVector s [VertexID])
forall (f :: * -> *) a. Applicative f => a -> f a
pure MVector s [VertexID]
v
  where
    updateAt :: MVector (PrimState m) a -> VertexID -> (a -> a) -> m ()
updateAt MVector (PrimState m) a
v VertexID
i a -> a
f = MVector (PrimState m) a -> VertexID -> m a
forall (m :: * -> *) a.
PrimMonad m =>
MVector (PrimState m) a -> VertexID -> m a
MV.read MVector (PrimState m) a
v VertexID
i m a -> (a -> m ()) -> m ()
forall (m :: * -> *) a b. Monad m => m a -> (a -> m b) -> m b
>>= \a
x -> MVector (PrimState m) a -> VertexID -> a -> m ()
forall (m :: * -> *) a.
PrimMonad m =>
MVector (PrimState m) a -> VertexID -> a -> m ()
MV.write MVector (PrimState m) a
v VertexID
i (a -> a
f a
x)
    addAt :: MVector (PrimState m) [a] -> VertexID -> a -> m ()
addAt    MVector (PrimState m) [a]
v VertexID
i a
j = MVector (PrimState m) [a] -> VertexID -> ([a] -> [a]) -> m ()
forall (m :: * -> *) a.
PrimMonad m =>
MVector (PrimState m) a -> VertexID -> (a -> a) -> m ()
updateAt MVector (PrimState m) [a]
v VertexID
i (a
ja -> [a] -> [a]
forall a. a -> [a] -> [a]
:)

    -- convert to a CList, sorted in CCW order around point u
    toCList :: VertexID -> [VertexID] -> CList VertexID
toCList VertexID
u = [VertexID] -> CList VertexID
forall a. [a] -> CList a
C.fromList ([VertexID] -> CList VertexID)
-> ([VertexID] -> [VertexID]) -> [VertexID] -> CList VertexID
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Mapping p r -> VertexID -> [VertexID] -> [VertexID]
forall r p.
(Num r, Ord r) =>
Mapping p r -> VertexID -> [VertexID] -> [VertexID]
sortAroundMapping Mapping p r
m VertexID
u

-- | Given a particular point u and a list of points vs, sort the points vs in
-- CW order around u.
-- running time: \( O(m log m) \), where m=|vs| is the number of vertices to sort.
sortAroundMapping               :: (Num r, Ord r)
                          => Mapping p r -> VertexID -> [VertexID] -> [VertexID]
sortAroundMapping :: Mapping p r -> VertexID -> [VertexID] -> [VertexID]
sortAroundMapping (Map (Point 2 r) VertexID
_,Vector (Point 2 r :+ p)
ptsV) VertexID
u [VertexID]
vs = [VertexID] -> [VertexID]
forall a. [a] -> [a]
reverse ([VertexID] -> [VertexID])
-> ([Point 2 r :+ VertexID] -> [VertexID])
-> [Point 2 r :+ VertexID]
-> [VertexID]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. ((Point 2 r :+ VertexID) -> VertexID)
-> [Point 2 r :+ VertexID] -> [VertexID]
forall a b. (a -> b) -> [a] -> [b]
map ((Point 2 r :+ VertexID)
-> Getting VertexID (Point 2 r :+ VertexID) VertexID -> VertexID
forall s a. s -> Getting a s a -> a
^.Getting VertexID (Point 2 r :+ VertexID) VertexID
forall core extra extra'.
Lens (core :+ extra) (core :+ extra') extra extra'
extra) ([Point 2 r :+ VertexID] -> [VertexID])
-> [Point 2 r :+ VertexID] -> [VertexID]
forall a b. (a -> b) -> a -> b
$ (Point 2 r :+ VertexID)
-> [Point 2 r :+ VertexID] -> [Point 2 r :+ VertexID]
forall r q p.
(Ord r, Num r) =>
(Point 2 r :+ q) -> [Point 2 r :+ p] -> [Point 2 r :+ p]
sortAround' (VertexID -> Point 2 r :+ VertexID
f VertexID
u) ((VertexID -> Point 2 r :+ VertexID)
-> [VertexID] -> [Point 2 r :+ VertexID]
forall a b. (a -> b) -> [a] -> [b]
map VertexID -> Point 2 r :+ VertexID
f [VertexID]
vs)
  where
    f :: VertexID -> Point 2 r :+ VertexID
f VertexID
v = (Vector (Point 2 r :+ p)
ptsV Vector (Point 2 r :+ p) -> VertexID -> Point 2 r :+ p
forall a. Vector a -> VertexID -> a
V.! VertexID
v)(Point 2 r :+ p)
-> ((Point 2 r :+ p) -> Point 2 r :+ VertexID)
-> Point 2 r :+ VertexID
forall a b. a -> (a -> b) -> b
&(p -> Identity VertexID)
-> (Point 2 r :+ p) -> Identity (Point 2 r :+ VertexID)
forall core extra extra'.
Lens (core :+ extra) (core :+ extra') extra extra'
extra ((p -> Identity VertexID)
 -> (Point 2 r :+ p) -> Identity (Point 2 r :+ VertexID))
-> VertexID -> (Point 2 r :+ p) -> Point 2 r :+ VertexID
forall s t a b. ASetter s t a b -> b -> s -> t
.~ VertexID
v

-- | Given a list of faces, construct a list of edges
extractEdges :: [(VertexID,VertexID,VertexID)] -> [(VertexID,VertexID)]
extractEdges :: [(VertexID, VertexID, VertexID)] -> [(VertexID, VertexID)]
extractEdges = ([(VertexID, VertexID)] -> (VertexID, VertexID))
-> [[(VertexID, VertexID)]] -> [(VertexID, VertexID)]
forall a b. (a -> b) -> [a] -> [b]
map [(VertexID, VertexID)] -> (VertexID, VertexID)
forall a. [a] -> a
L.head ([[(VertexID, VertexID)]] -> [(VertexID, VertexID)])
-> ([(VertexID, VertexID, VertexID)] -> [[(VertexID, VertexID)]])
-> [(VertexID, VertexID, VertexID)]
-> [(VertexID, VertexID)]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [(VertexID, VertexID)] -> [[(VertexID, VertexID)]]
forall a. Eq a => [a] -> [[a]]
L.group ([(VertexID, VertexID)] -> [[(VertexID, VertexID)]])
-> ([(VertexID, VertexID, VertexID)] -> [(VertexID, VertexID)])
-> [(VertexID, VertexID, VertexID)]
-> [[(VertexID, VertexID)]]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [(VertexID, VertexID)] -> [(VertexID, VertexID)]
forall a. Ord a => [a] -> [a]
L.sort
               ([(VertexID, VertexID)] -> [(VertexID, VertexID)])
-> ([(VertexID, VertexID, VertexID)] -> [(VertexID, VertexID)])
-> [(VertexID, VertexID, VertexID)]
-> [(VertexID, VertexID)]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. ((VertexID, VertexID, VertexID) -> [(VertexID, VertexID)])
-> [(VertexID, VertexID, VertexID)] -> [(VertexID, VertexID)]
forall (t :: * -> *) a b. Foldable t => (a -> [b]) -> t a -> [b]
concatMap (\(VertexID
p,VertexID
q,VertexID
r) -> [(VertexID
p,VertexID
q), (VertexID
q,VertexID
r), (VertexID
p,VertexID
r)])
               -- we encounter every edge twice. To get rid of the duplicates
               -- we sort, group, and take the head of the lists


-- | \( O(n) \) Test if the given three points form a triangle in the delaunay triangulation.
isDelaunay                :: (Fractional r, Ord r)
                          => Mapping p r -> VertexID -> VertexID -> VertexID -> Bool
isDelaunay :: Mapping p r -> VertexID -> VertexID -> VertexID -> Bool
isDelaunay (Map (Point 2 r) VertexID
_,Vector (Point 2 r :+ p)
ptsV) VertexID
p VertexID
q VertexID
r = case Point 2 r -> Point 2 r -> Point 2 r -> Maybe (Disk () r)
forall r.
(Eq r, Fractional r) =>
Point 2 r -> Point 2 r -> Point 2 r -> Maybe (Disk () r)
disk (VertexID -> Point 2 r
pt VertexID
p) (VertexID -> Point 2 r
pt VertexID
q) (VertexID -> Point 2 r
pt VertexID
r) of
    Maybe (Disk () r)
Nothing -> Bool
False -- if the points are colinear, we interpret this as: all
                     -- pts in the plane are in the circle.
    Just Disk () r
d  -> Bool -> Bool
not (Bool -> Bool) -> Bool -> Bool
forall a b. (a -> b) -> a -> b
$ (Point 2 r -> Bool) -> [Point 2 r] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
any (Point 2 r -> Disk () r -> Bool
forall (d :: Nat) r p.
(Arity d, Ord r, Num r) =>
Point d r -> Ball d p r -> Bool
`insideBall` Disk () r
d)
      [VertexID -> Point 2 r
pt VertexID
i | VertexID
i <- [VertexID
0..(Vector (Point 2 r :+ p) -> VertexID
forall a. Vector a -> VertexID
V.length Vector (Point 2 r :+ p)
ptsV VertexID -> VertexID -> VertexID
forall a. Num a => a -> a -> a
- VertexID
1)], VertexID
i VertexID -> VertexID -> Bool
forall a. Eq a => a -> a -> Bool
/= VertexID
p, VertexID
i VertexID -> VertexID -> Bool
forall a. Eq a => a -> a -> Bool
/= VertexID
q, VertexID
i VertexID -> VertexID -> Bool
forall a. Eq a => a -> a -> Bool
/= VertexID
r]
   where
     pt :: VertexID -> Point 2 r
pt VertexID
i = (Vector (Point 2 r :+ p)
ptsV Vector (Point 2 r :+ p) -> VertexID -> Point 2 r :+ p
forall a. Vector a -> VertexID -> a
V.! VertexID
i)(Point 2 r :+ p)
-> Getting (Point 2 r) (Point 2 r :+ p) (Point 2 r) -> Point 2 r
forall s a. s -> Getting a s a -> a
^.Getting (Point 2 r) (Point 2 r :+ p) (Point 2 r)
forall core extra core'.
Lens (core :+ extra) (core' :+ extra) core core'
core