module Geometry.Trajectory.PathSimplifier.VisvalingamWhyatt (
simplifyTrajectoryVW
, simplifyTrajectoryVWBy
) where
import Control.Monad
import Control.Monad.ST
import Data.Foldable
import Data.Heap (Entry (..), Heap)
import qualified Data.Heap as H
import Data.Maybe
import Data.Vector (Vector, (!?))
import qualified Data.Vector as V
import qualified Data.Vector.Mutable as VM
import Geometry.Core
data Triangle = Triangle
{ Triangle -> Int
_left :: Int
, Triangle -> Int
_center :: Int
, Triangle -> Int
_right :: Int
} deriving (Triangle -> Triangle -> Bool
(Triangle -> Triangle -> Bool)
-> (Triangle -> Triangle -> Bool) -> Eq Triangle
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: Triangle -> Triangle -> Bool
== :: Triangle -> Triangle -> Bool
$c/= :: Triangle -> Triangle -> Bool
/= :: Triangle -> Triangle -> Bool
Eq, Eq Triangle
Eq Triangle
-> (Triangle -> Triangle -> Ordering)
-> (Triangle -> Triangle -> Bool)
-> (Triangle -> Triangle -> Bool)
-> (Triangle -> Triangle -> Bool)
-> (Triangle -> Triangle -> Bool)
-> (Triangle -> Triangle -> Triangle)
-> (Triangle -> Triangle -> Triangle)
-> Ord Triangle
Triangle -> Triangle -> Bool
Triangle -> Triangle -> Ordering
Triangle -> Triangle -> Triangle
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 :: Triangle -> Triangle -> Ordering
compare :: Triangle -> Triangle -> Ordering
$c< :: Triangle -> Triangle -> Bool
< :: Triangle -> Triangle -> Bool
$c<= :: Triangle -> Triangle -> Bool
<= :: Triangle -> Triangle -> Bool
$c> :: Triangle -> Triangle -> Bool
> :: Triangle -> Triangle -> Bool
$c>= :: Triangle -> Triangle -> Bool
>= :: Triangle -> Triangle -> Bool
$cmax :: Triangle -> Triangle -> Triangle
max :: Triangle -> Triangle -> Triangle
$cmin :: Triangle -> Triangle -> Triangle
min :: Triangle -> Triangle -> Triangle
Ord, Int -> Triangle -> ShowS
[Triangle] -> ShowS
Triangle -> String
(Int -> Triangle -> ShowS)
-> (Triangle -> String) -> ([Triangle] -> ShowS) -> Show Triangle
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> Triangle -> ShowS
showsPrec :: Int -> Triangle -> ShowS
$cshow :: Triangle -> String
show :: Triangle -> String
$cshowList :: [Triangle] -> ShowS
showList :: [Triangle] -> ShowS
Show)
triangleArea :: Vec2 -> Vec2 -> Vec2 -> Double
triangleArea :: Vec2 -> Vec2 -> Vec2 -> Double
triangleArea Vec2
left Vec2
center Vec2
right = Double -> Double
forall a. Num a => a -> a
abs (Vec2 -> Vec2 -> Double
cross (Vec2
left Vec2 -> Vec2 -> Vec2
forall v. VectorSpace v => v -> v -> v
-. Vec2
center) (Vec2
right Vec2 -> Vec2 -> Vec2
forall v. VectorSpace v => v -> v -> v
-. Vec2
center))
mkTriangleAreaPQ :: Vector Vec2 -> Heap (Entry Double Triangle)
mkTriangleAreaPQ :: Vector Vec2 -> Heap (Entry Double Triangle)
mkTriangleAreaPQ Vector Vec2
vec
| Vector Vec2 -> Bool
forall a. Eq a => Vector a -> Bool
isCyclic Vector Vec2
vec = Vector Vec2 -> Heap (Entry Double Triangle)
mkTriangleAreaPQCyclic Vector Vec2
vec
| Bool
otherwise = Vector Vec2 -> Heap (Entry Double Triangle)
mkTriangleAreaPQLinear Vector Vec2
vec
mkTriangleAreaPQCyclic :: Vector Vec2 -> Heap (Entry Double Triangle)
mkTriangleAreaPQCyclic :: Vector Vec2 -> Heap (Entry Double Triangle)
mkTriangleAreaPQCyclic Vector Vec2
vec = Vector (Entry Double Triangle) -> Heap (Entry Double Triangle)
forall a. Ord a => Vector a -> Heap a
heapFromVector (Vector (Entry Double Triangle) -> Heap (Entry Double Triangle))
-> Vector (Entry Double Triangle) -> Heap (Entry Double Triangle)
forall a b. (a -> b) -> a -> b
$ (Int -> Vec2 -> Vec2 -> Vec2 -> Entry Double Triangle)
-> Vector Vec2
-> Vector Vec2
-> Vector Vec2
-> Vector (Entry Double Triangle)
forall a b c d.
(Int -> a -> b -> c -> d)
-> Vector a -> Vector b -> Vector c -> Vector d
V.izipWith3
(\Int
i Vec2
a Vec2
b Vec2
c ->
let l :: Int
l = Vector Vec2 -> Int
forall a. Vector a -> Int
V.length Vector Vec2
vec
in Double -> Triangle -> Entry Double Triangle
forall p a. p -> a -> Entry p a
Entry
(Vec2 -> Vec2 -> Vec2 -> Double
triangleArea Vec2
a Vec2
b Vec2
c)
Triangle
{ _left :: Int
_left = Int
i Int -> Int -> Int
forall a. Integral a => a -> a -> a
`mod` Int
l
, _center :: Int
_center = (Int
iInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1) Int -> Int -> Int
forall a. Integral a => a -> a -> a
`mod` Int
l
, _right :: Int
_right = (Int
iInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
2) Int -> Int -> Int
forall a. Integral a => a -> a -> a
`mod` Int
l })
Vector Vec2
vec
(Vector Vec2
vec Vector Vec2 -> Vector Vec2 -> Vector Vec2
forall a. Semigroup a => a -> a -> a
<> Int -> Vector Vec2 -> Vector Vec2
forall a. Int -> Vector a -> Vector a
V.take Int
2 Vector Vec2
vec)
(Vector Vec2
vec Vector Vec2 -> Vector Vec2 -> Vector Vec2
forall a. Semigroup a => a -> a -> a
<> Int -> Vector Vec2 -> Vector Vec2
forall a. Int -> Vector a -> Vector a
V.take Int
2 Vector Vec2
vec)
isCyclic :: Eq a => Vector a -> Bool
isCyclic :: forall a. Eq a => Vector a -> Bool
isCyclic Vector a
vec = Vector a -> a
forall a. Vector a -> a
V.head Vector a
vec a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== Vector a -> a
forall a. Vector a -> a
V.last Vector a
vec
mkTriangleAreaPQLinear :: Vector Vec2 -> Heap (Entry Double Triangle)
mkTriangleAreaPQLinear :: Vector Vec2 -> Heap (Entry Double Triangle)
mkTriangleAreaPQLinear Vector Vec2
vec = Vector (Entry Double Triangle) -> Heap (Entry Double Triangle)
forall a. Ord a => Vector a -> Heap a
heapFromVector (Vector (Entry Double Triangle) -> Heap (Entry Double Triangle))
-> Vector (Entry Double Triangle) -> Heap (Entry Double Triangle)
forall a b. (a -> b) -> a -> b
$ (Int -> Vec2 -> Vec2 -> Vec2 -> Entry Double Triangle)
-> Vector Vec2
-> Vector Vec2
-> Vector Vec2
-> Vector (Entry Double Triangle)
forall a b c d.
(Int -> a -> b -> c -> d)
-> Vector a -> Vector b -> Vector c -> Vector d
V.izipWith3
(\Int
i Vec2
a Vec2
b Vec2
c -> Double -> Triangle -> Entry Double Triangle
forall p a. p -> a -> Entry p a
Entry
(Vec2 -> Vec2 -> Vec2 -> Double
triangleArea Vec2
a Vec2
b Vec2
c)
Triangle
{ _left :: Int
_left = Int
i
, _center :: Int
_center = Int
iInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1
, _right :: Int
_right = Int
iInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
2 })
Vector Vec2
vec
(Int -> Vector Vec2 -> Vector Vec2
forall a. Int -> Vector a -> Vector a
V.drop Int
1 Vector Vec2
vec)
(Int -> Vector Vec2 -> Vector Vec2
forall a. Int -> Vector a -> Vector a
V.drop Int
2 Vector Vec2
vec)
heapFromVector :: Ord a => Vector a -> Heap a
heapFromVector :: forall a. Ord a => Vector a -> Heap a
heapFromVector = (Heap a -> a -> Heap a) -> Heap a -> Vector a -> Heap a
forall a b. (a -> b -> a) -> a -> Vector b -> a
V.foldl' (\Heap a
heap a
entry -> a -> Heap a -> Heap a
forall a. Ord a => a -> Heap a -> Heap a
H.insert a
entry Heap a
heap) Heap a
forall a. Monoid a => a
mempty
vwSimplifyIndices :: Double -> Vector Vec2 -> Vector Int
vwSimplifyIndices :: Double -> Vector Vec2 -> Vector Int
vwSimplifyIndices Double
minArea Vector Vec2
inputPoints = (forall s. ST s (Vector Int)) -> Vector Int
forall a. (forall s. ST s a) -> a
runST ((forall s. ST s (Vector Int)) -> Vector Int)
-> (forall s. ST s (Vector Int)) -> Vector Int
forall a b. (a -> b) -> a -> b
$ do
let
inputLength :: Int
inputLength = Vector Vec2 -> Int
forall a. Vector a -> Int
V.length Vector Vec2
inputPoints
ix :: Int -> Int
ix Int
i
| Vector Vec2 -> Bool
forall a. Eq a => Vector a -> Bool
isCyclic Vector Vec2
inputPoints = Int
i Int -> Int -> Int
forall a. Integral a => a -> a -> a
`mod` Int
inputLength
| Bool
otherwise = Int
i
MVector s (Int, Int)
adjacentMut <- Vector (Int, Int) -> ST s (MVector (PrimState (ST s)) (Int, Int))
forall (m :: * -> *) a.
PrimMonad m =>
Vector a -> m (MVector (PrimState m) a)
V.thaw (Int -> (Int -> (Int, Int)) -> Vector (Int, Int)
forall a. Int -> (Int -> a) -> Vector a
V.generate Int
inputLength (\Int
i -> (Int -> Int
ix (Int
iInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1), Int -> Int
ix (Int
iInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1))))
let vwLoop :: Heap (Entry Double Triangle) -> ST s ()
vwLoop Heap (Entry Double Triangle)
heap = case Heap (Entry Double Triangle)
-> Maybe (Entry Double Triangle, Heap (Entry Double Triangle))
forall a. Heap a -> Maybe (a, Heap a)
H.uncons Heap (Entry Double Triangle)
heap of
Maybe (Entry Double Triangle, Heap (Entry Double Triangle))
Nothing -> () -> ST s ()
forall a. a -> ST s a
forall (f :: * -> *) a. Applicative f => a -> f a
pure ()
Just (Entry Double
area Triangle
triangle, Heap (Entry Double Triangle)
restOfHeap)
| Double
area Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>= Double
minArea -> Heap (Entry Double Triangle) -> ST s ()
vwLoop Heap (Entry Double Triangle)
restOfHeap
| Bool
otherwise -> do
{ (Int
left, Int
right) <- MVector (PrimState (ST s)) (Int, Int) -> Int -> ST s (Int, Int)
forall (m :: * -> *) a.
PrimMonad m =>
MVector (PrimState m) a -> Int -> m a
VM.read MVector s (Int, Int)
MVector (PrimState (ST s)) (Int, Int)
adjacentMut (Triangle -> Int
_center Triangle
triangle)
; if Int
left Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Triangle -> Int
_left Triangle
triangle Bool -> Bool -> Bool
|| Int
right Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Triangle -> Int
_right Triangle
triangle
then
Heap (Entry Double Triangle) -> ST s ()
vwLoop Heap (Entry Double Triangle)
restOfHeap
else do
{ (Int
ll, Int
rr) <- MVector s (Int, Int) -> Int -> Triangle -> Int -> ST s (Int, Int)
forall s.
STVector s (Int, Int) -> Int -> Triangle -> Int -> ST s (Int, Int)
deleteTriangle MVector s (Int, Int)
adjacentMut Int
left Triangle
triangle Int
right
; let
newTriangles :: [Triangle]
newTriangles = [Int -> Int -> Int -> Triangle
Triangle Int
ll Int
left Int
right, Int -> Int -> Int -> Triangle
Triangle Int
left Int
right Int
rr]
newHeap :: Heap (Entry Double Triangle)
newHeap = Vector Vec2
-> [Triangle]
-> Heap (Entry Double Triangle)
-> Heap (Entry Double Triangle)
insertNewTriangles Vector Vec2
inputPoints [Triangle]
newTriangles Heap (Entry Double Triangle)
restOfHeap
; Heap (Entry Double Triangle) -> ST s ()
vwLoop Heap (Entry Double Triangle)
newHeap
}}
Heap (Entry Double Triangle) -> ST s ()
vwLoop (Vector Vec2 -> Heap (Entry Double Triangle)
mkTriangleAreaPQ Vector Vec2
inputPoints)
Vector (Int, Int)
adjacent <- MVector (PrimState (ST s)) (Int, Int) -> ST s (Vector (Int, Int))
forall (m :: * -> *) a.
PrimMonad m =>
MVector (PrimState m) a -> m (Vector a)
V.freeze MVector s (Int, Int)
MVector (PrimState (ST s)) (Int, Int)
adjacentMut
Vector Int -> ST s (Vector Int)
forall a. a -> ST s a
forall (f :: * -> *) a. Applicative f => a -> f a
pure (Vector Int -> ST s (Vector Int))
-> Vector Int -> ST s (Vector Int)
forall a b. (a -> b) -> a -> b
$ Vector (Maybe Int) -> Vector Int
forall a. Vector (Maybe a) -> Vector a
V.catMaybes (Vector (Maybe Int) -> Vector Int)
-> Vector (Maybe Int) -> Vector Int
forall a b. (a -> b) -> a -> b
$ (Int -> Vec2 -> (Int, Int) -> Maybe Int)
-> Vector Vec2 -> Vector (Int, Int) -> Vector (Maybe Int)
forall a b c.
(Int -> a -> b -> c) -> Vector a -> Vector b -> Vector c
V.izipWith
(\Int
i Vec2
_ (Int, Int)
adj ->
if (Int, Int)
adj (Int, Int) -> (Int, Int) -> Bool
forall a. Eq a => a -> a -> Bool
/= (Int
0,Int
0)
then Int -> Maybe Int
forall a. a -> Maybe a
Just Int
i
else Maybe Int
forall a. Maybe a
Nothing)
Vector Vec2
inputPoints
Vector (Int, Int)
adjacent
insertNewTriangles
:: Vector Vec2
-> [Triangle]
-> Heap (Entry Double Triangle)
-> Heap (Entry Double Triangle)
insertNewTriangles :: Vector Vec2
-> [Triangle]
-> Heap (Entry Double Triangle)
-> Heap (Entry Double Triangle)
insertNewTriangles Vector Vec2
inputPoints [Triangle]
current Heap (Entry Double Triangle)
restOfHeap =
(Heap (Entry Double Triangle)
-> Triangle -> Heap (Entry Double Triangle))
-> Heap (Entry Double Triangle)
-> [Triangle]
-> Heap (Entry Double Triangle)
forall b a. (b -> a -> b) -> b -> [a] -> b
forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl' (Vector Vec2
-> Heap (Entry Double Triangle)
-> Triangle
-> Heap (Entry Double Triangle)
insertNewTriangle Vector Vec2
inputPoints) Heap (Entry Double Triangle)
restOfHeap [Triangle]
current
insertNewTriangle
:: Vector Vec2
-> Heap (Entry Double Triangle)
-> Triangle
-> Heap (Entry Double Triangle)
insertNewTriangle :: Vector Vec2
-> Heap (Entry Double Triangle)
-> Triangle
-> Heap (Entry Double Triangle)
insertNewTriangle Vector Vec2
inputPoints Heap (Entry Double Triangle)
heap triangle :: Triangle
triangle@(Triangle Int
ai Int
bi Int
ci) = Heap (Entry Double Triangle)
-> Maybe (Heap (Entry Double Triangle))
-> Heap (Entry Double Triangle)
forall a. a -> Maybe a -> a
fromMaybe Heap (Entry Double Triangle)
heap (Maybe (Heap (Entry Double Triangle))
-> Heap (Entry Double Triangle))
-> Maybe (Heap (Entry Double Triangle))
-> Heap (Entry Double Triangle)
forall a b. (a -> b) -> a -> b
$ do
Bool -> Maybe ()
forall (f :: * -> *). Alternative f => Bool -> f ()
guard (Int
ai Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
bi)
Bool -> Maybe ()
forall (f :: * -> *). Alternative f => Bool -> f ()
guard (Int
ai Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
ci)
Bool -> Maybe ()
forall (f :: * -> *). Alternative f => Bool -> f ()
guard (Int
bi Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
ci)
Vec2
a <- Vector Vec2
inputPoints Vector Vec2 -> Int -> Maybe Vec2
forall a. Vector a -> Int -> Maybe a
!? Int
ai
Vec2
b <- Vector Vec2
inputPoints Vector Vec2 -> Int -> Maybe Vec2
forall a. Vector a -> Int -> Maybe a
!? Int
bi
Vec2
c <- Vector Vec2
inputPoints Vector Vec2 -> Int -> Maybe Vec2
forall a. Vector a -> Int -> Maybe a
!? Int
ci
let newArea :: Double
newArea = Vec2 -> Vec2 -> Vec2 -> Double
triangleArea Vec2
a Vec2
b Vec2
c
Heap (Entry Double Triangle)
-> Maybe (Heap (Entry Double Triangle))
forall a. a -> Maybe a
forall (f :: * -> *) a. Applicative f => a -> f a
pure (Entry Double Triangle
-> Heap (Entry Double Triangle) -> Heap (Entry Double Triangle)
forall a. Ord a => a -> Heap a -> Heap a
H.insert (Double -> Triangle -> Entry Double Triangle
forall p a. p -> a -> Entry p a
Entry Double
newArea Triangle
triangle) Heap (Entry Double Triangle)
heap)
deleteTriangle
:: VM.STVector s (Int, Int)
-> Int
-> Triangle
-> Int
-> ST s (Int, Int)
deleteTriangle :: forall s.
STVector s (Int, Int) -> Int -> Triangle -> Int -> ST s (Int, Int)
deleteTriangle STVector s (Int, Int)
adjacentMut Int
newLeft Triangle
center Int
newRight = do
(Int
ll, Int
_l) <- MVector (PrimState (ST s)) (Int, Int) -> Int -> ST s (Int, Int)
forall (m :: * -> *) a.
PrimMonad m =>
MVector (PrimState m) a -> Int -> m a
VM.read STVector s (Int, Int)
MVector (PrimState (ST s)) (Int, Int)
adjacentMut (Triangle -> Int
_left Triangle
center)
(Int
_r, Int
rr) <- MVector (PrimState (ST s)) (Int, Int) -> Int -> ST s (Int, Int)
forall (m :: * -> *) a.
PrimMonad m =>
MVector (PrimState m) a -> Int -> m a
VM.read STVector s (Int, Int)
MVector (PrimState (ST s)) (Int, Int)
adjacentMut (Triangle -> Int
_right Triangle
center)
MVector (PrimState (ST s)) (Int, Int)
-> Int -> (Int, Int) -> ST s ()
forall (m :: * -> *) a.
PrimMonad m =>
MVector (PrimState m) a -> Int -> a -> m ()
VM.write STVector s (Int, Int)
MVector (PrimState (ST s)) (Int, Int)
adjacentMut (Triangle -> Int
_left Triangle
center) (Int
ll, Int
newRight)
MVector (PrimState (ST s)) (Int, Int)
-> Int -> (Int, Int) -> ST s ()
forall (m :: * -> *) a.
PrimMonad m =>
MVector (PrimState m) a -> Int -> a -> m ()
VM.write STVector s (Int, Int)
MVector (PrimState (ST s)) (Int, Int)
adjacentMut (Triangle -> Int
_right Triangle
center) (Int
newLeft, Int
rr)
MVector (PrimState (ST s)) (Int, Int)
-> Int -> (Int, Int) -> ST s ()
forall (m :: * -> *) a.
PrimMonad m =>
MVector (PrimState m) a -> Int -> a -> m ()
VM.write STVector s (Int, Int)
MVector (PrimState (ST s)) (Int, Int)
adjacentMut (Triangle -> Int
_center Triangle
center) (Int
0, Int
0)
(Int, Int) -> ST s (Int, Int)
forall a. a -> ST s a
forall (f :: * -> *) a. Applicative f => a -> f a
pure (Int
ll, Int
rr)
simplifyTrajectoryVW
:: Sequential vector
=> Double
-> vector Vec2
-> Vector Vec2
simplifyTrajectoryVW :: forall (vector :: * -> *).
Sequential vector =>
Double -> vector Vec2 -> Vector Vec2
simplifyTrajectoryVW Double
epsilon vector Vec2
trajectorySequence =
let inputPoints :: Vector Vec2
inputPoints = vector Vec2 -> Vector Vec2
forall a. vector a -> Vector a
forall (f :: * -> *) a. Sequential f => f a -> Vector a
toVector vector Vec2
trajectorySequence
indices :: Vector Int
indices = Double -> Vector Vec2 -> Vector Int
vwSimplifyIndices Double
epsilon Vector Vec2
inputPoints
in Vector Vec2 -> Vector Int -> Vector Vec2
forall a. Vector a -> Vector Int -> Vector a
V.backpermute Vector Vec2
inputPoints Vector Int
indices
simplifyTrajectoryVWBy
:: (Ord a, Sequential vector)
=> Double
-> (a -> Vec2)
-> vector a
-> Vector a
simplifyTrajectoryVWBy :: forall a (vector :: * -> *).
(Ord a, Sequential vector) =>
Double -> (a -> Vec2) -> vector a -> Vector a
simplifyTrajectoryVWBy Double
epsilon a -> Vec2
vec2in vector a
trajectorySequence =
let inputPoints :: Vector a
inputPoints = vector a -> Vector a
forall a. vector a -> Vector a
forall (f :: * -> *) a. Sequential f => f a -> Vector a
toVector vector a
trajectorySequence
indices :: Vector Int
indices = Double -> Vector Vec2 -> Vector Int
vwSimplifyIndices Double
epsilon ((a -> Vec2) -> Vector a -> Vector Vec2
forall a b. (a -> b) -> Vector a -> Vector b
V.map a -> Vec2
vec2in Vector a
inputPoints)
in Vector a -> Vector Int -> Vector a
forall a. Vector a -> Vector Int -> Vector a
V.backpermute Vector a
inputPoints Vector Int
indices