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

import           Control.Lens hiding (only)
import           Data.Ext
import qualified Data.Foldable as F
import           Data.Geometry.LineSegment
import           Data.Geometry.Point
import           Data.Geometry.PolyLine
import           Data.Geometry.Vector
import qualified Data.LSeq as LSeq
import           Data.LSeq (LSeq, pattern (:|>))
import           Data.Ord (comparing)

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

-- | Line simplification with the well-known Douglas Peucker alogrithm. Given a distance
-- value eps adn a polyline pl, constructs a simplification of pl (i.e. with
-- vertices from pl) s.t. all other vertices are within dist eps to the
-- original polyline.
--
-- Running time: \( O(n^2) \) worst case, \( O(n log n) \) expected.
douglasPeucker         :: (Ord r, Fractional r, Arity d)
                       => r -> PolyLine d p r -> PolyLine d p r
douglasPeucker :: r -> PolyLine d p r -> PolyLine d p r
douglasPeucker r
eps PolyLine d p r
pl
    | r
dst r -> r -> Bool
forall a. Ord a => a -> a -> Bool
<= (r
epsr -> r -> r
forall a. Num a => a -> a -> a
*r
eps) = [Point d r :+ p] -> PolyLine d p r
forall (d :: Nat) r p. [Point d r :+ p] -> PolyLine d p r
fromPointsUnsafe [Point d r :+ p
a,Point d r :+ p
b] -- at least two points, so we are fine.
    | Bool
otherwise        = r -> PolyLine d p r -> PolyLine d p r
forall r (d :: Nat) p.
(Ord r, Fractional r, Arity d) =>
r -> PolyLine d p r -> PolyLine d p r
douglasPeucker r
eps PolyLine d p r
pref PolyLine d p r -> PolyLine d p r -> PolyLine d p r
forall (d :: Nat) p r.
PolyLine d p r -> PolyLine d p r -> PolyLine d p r
`merge` r -> PolyLine d p r -> PolyLine d p r
forall r (d :: Nat) p.
(Ord r, Fractional r, Arity d) =>
r -> PolyLine d p r -> PolyLine d p r
douglasPeucker r
eps PolyLine d p r
subf
  where
    pts :: LSeq 2 (Point d r :+ p)
pts         = PolyLine d p r
plPolyLine d p r
-> Getting
     (LSeq 2 (Point d r :+ p))
     (PolyLine d p r)
     (LSeq 2 (Point d r :+ p))
-> LSeq 2 (Point d r :+ p)
forall s a. s -> Getting a s a -> a
^.Getting
  (LSeq 2 (Point d r :+ p))
  (PolyLine d p r)
  (LSeq 2 (Point d r :+ p))
forall (d1 :: Nat) p1 r1 (d2 :: Nat) p2 r2.
Iso
  (PolyLine d1 p1 r1)
  (PolyLine d2 p2 r2)
  (LSeq 2 (Point d1 r1 :+ p1))
  (LSeq 2 (Point d2 r2 :+ p2))
points
    a :: Point d r :+ p
a           = LSeq (1 + 1) (Point d r :+ p) -> Point d r :+ p
forall (n :: Nat) a. LSeq (1 + n) a -> a
LSeq.head LSeq 2 (Point d r :+ p)
LSeq (1 + 1) (Point d r :+ p)
pts
    b :: Point d r :+ p
b           = LSeq (1 + 1) (Point d r :+ p) -> Point d r :+ p
forall (n :: Nat) a. LSeq (1 + n) a -> a
LSeq.last LSeq 2 (Point d r :+ p)
LSeq (1 + 1) (Point d r :+ p)
pts
    (Int
i,r
dst)     = LSeq 2 (Point d r :+ p) -> LineSegment d p r -> (Int, r)
forall r (d :: Nat) (n :: Nat) p.
(Ord r, Fractional r, Arity d) =>
LSeq n (Point d r :+ p) -> LineSegment d p r -> (Int, r)
maxDist LSeq 2 (Point d r :+ p)
pts ((Point d r :+ p) -> (Point d r :+ p) -> LineSegment d p r
forall (d :: Nat) r p.
(Point d r :+ p) -> (Point d r :+ p) -> LineSegment d p r
ClosedLineSegment Point d r :+ p
a Point d r :+ p
b)

    (PolyLine d p r
pref,PolyLine d p r
subf) = Int -> PolyLine d p r -> (PolyLine d p r, PolyLine d p r)
forall (d :: Nat) p r.
Int -> PolyLine d p r -> (PolyLine d p r, PolyLine d p r)
split Int
i PolyLine d p r
pl

--------------------------------------------------------------------------------
-- * Internal functions

-- | Concatenate the two polylines, dropping their shared vertex
merge          :: PolyLine d p r -> PolyLine d p r -> PolyLine d p r
merge :: PolyLine d p r -> PolyLine d p r -> PolyLine d p r
merge PolyLine d p r
pref PolyLine d p r
sub = LSeq 2 (Point d r :+ p) -> PolyLine d p r
forall (d :: Nat) p r. LSeq 2 (Point d r :+ p) -> PolyLine d p r
PolyLine (LSeq 2 (Point d r :+ p) -> PolyLine d p r)
-> LSeq 2 (Point d r :+ p) -> PolyLine d p r
forall a b. (a -> b) -> a -> b
$ LSeq 1 (Point d r :+ p)
pref' LSeq 1 (Point d r :+ p)
-> LSeq 2 (Point d r :+ p) -> LSeq 2 (Point d r :+ p)
forall (n :: Nat) a (m :: Nat). LSeq n a -> LSeq m a -> LSeq m a
`append` (PolyLine d p r
subPolyLine d p r
-> Getting
     (LSeq 2 (Point d r :+ p))
     (PolyLine d p r)
     (LSeq 2 (Point d r :+ p))
-> LSeq 2 (Point d r :+ p)
forall s a. s -> Getting a s a -> a
^.Getting
  (LSeq 2 (Point d r :+ p))
  (PolyLine d p r)
  (LSeq 2 (Point d r :+ p))
forall (d1 :: Nat) p1 r1 (d2 :: Nat) p2 r2.
Iso
  (PolyLine d1 p1 r1)
  (PolyLine d2 p2 r2)
  (LSeq 2 (Point d1 r1 :+ p1))
  (LSeq 2 (Point d2 r2 :+ p2))
points)
  where
    (LSeq 1 (Point d r :+ p)
pref' :|> Point d r :+ p
_) = PolyLine d p r
prefPolyLine d p r
-> Getting
     (LSeq 2 (Point d r :+ p))
     (PolyLine d p r)
     (LSeq 2 (Point d r :+ p))
-> LSeq 2 (Point d r :+ p)
forall s a. s -> Getting a s a -> a
^.Getting
  (LSeq 2 (Point d r :+ p))
  (PolyLine d p r)
  (LSeq 2 (Point d r :+ p))
forall (d1 :: Nat) p1 r1 (d2 :: Nat) p2 r2.
Iso
  (PolyLine d1 p1 r1)
  (PolyLine d2 p2 r2)
  (LSeq 2 (Point d1 r1 :+ p1))
  (LSeq 2 (Point d2 r2 :+ p2))
points
    append     :: LSeq.LSeq n a ->  LSeq.LSeq m a -> LSeq.LSeq m a
    append :: LSeq n a -> LSeq m a -> LSeq m a
append LSeq n a
a LSeq m a
b = LSeq (n + m) a -> LSeq m a
forall (m :: Nat) (n :: Nat) a. LSeq m a -> LSeq n a
LSeq.promise (LSeq (n + m) a -> LSeq m a) -> LSeq (n + m) a -> LSeq m a
forall a b. (a -> b) -> a -> b
$ LSeq n a -> LSeq m a -> LSeq (n + m) a
forall (n :: Nat) a (m :: Nat).
LSeq n a -> LSeq m a -> LSeq (n + m) a
LSeq.append LSeq n a
a LSeq m a
b
    -- our seq is actually even longer


-- | Split the polyline at the given vertex. Both polylines contain this vertex
split                  :: Int -> PolyLine d p r
                       -> (PolyLine d p r, PolyLine d p r)
split :: Int -> PolyLine d p r -> (PolyLine d p r, PolyLine d p r)
split Int
i (PolyLine LSeq 2 (Point d r :+ p)
pts) = (LSeq 0 (Point d r :+ p) -> PolyLine d p r)
-> (LSeq 0 (Point d r :+ p) -> PolyLine d p r)
-> (LSeq 0 (Point d r :+ p), LSeq 0 (Point d r :+ p))
-> (PolyLine d p r, PolyLine d p r)
forall (p :: * -> * -> *) a b c d.
Bifunctor p =>
(a -> b) -> (c -> d) -> p a c -> p b d
bimap LSeq 0 (Point d r :+ p) -> PolyLine d p r
forall (m :: Nat) (d :: Nat) r p.
LSeq m (Point d r :+ p) -> PolyLine d p r
f LSeq 0 (Point d r :+ p) -> PolyLine d p r
forall (m :: Nat) (d :: Nat) r p.
LSeq m (Point d r :+ p) -> PolyLine d p r
f (LSeq 0 (Point d r :+ p)
as,LSeq 0 (Point d r :+ p)
bs)
  where
    f :: LSeq m (Point d r :+ p) -> PolyLine d p r
f = LSeq 2 (Point d r :+ p) -> PolyLine d p r
forall (d :: Nat) p r. LSeq 2 (Point d r :+ p) -> PolyLine d p r
PolyLine (LSeq 2 (Point d r :+ p) -> PolyLine d p r)
-> (LSeq m (Point d r :+ p) -> LSeq 2 (Point d r :+ p))
-> LSeq m (Point d r :+ p)
-> PolyLine d p r
forall b c a. (b -> c) -> (a -> b) -> a -> c
. C 2 -> LSeq m (Point d r :+ p) -> LSeq 2 (Point d r :+ p)
forall (n :: Nat) (proxy :: Nat -> *) (m :: Nat) a.
KnownNat n =>
proxy n -> LSeq m a -> LSeq n a
LSeq.forceLSeq (C 2
forall (n :: Nat). C n
C  :: C 2)
    as :: LSeq 0 (Point d r :+ p)
as = Int -> LSeq 2 (Point d r :+ p) -> LSeq 0 (Point d r :+ p)
forall (n :: Nat) a. Int -> LSeq n a -> LSeq 0 a
LSeq.take (Int
iInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1) LSeq 2 (Point d r :+ p)
pts
    bs :: LSeq 0 (Point d r :+ p)
bs = Int -> LSeq 2 (Point d r :+ p) -> LSeq 0 (Point d r :+ p)
forall (n :: Nat) a. Int -> LSeq n a -> LSeq 0 a
LSeq.drop Int
i     LSeq 2 (Point d r :+ p)
pts

-- | Given a sequence of points, find the index of the point that has the
-- Furthest distance to the LineSegment. The result is the index of the point
-- and this distance.
maxDist       :: (Ord r, Fractional r, Arity d)
              => LSeq n (Point d r :+ p) -> LineSegment d p r -> (Int,r)
maxDist :: LSeq n (Point d r :+ p) -> LineSegment d p r -> (Int, r)
maxDist LSeq n (Point d r :+ p)
pts LineSegment d p r
s = ((Int, r) -> (Int, r) -> Ordering) -> LSeq n (Int, r) -> (Int, r)
forall (t :: * -> *) a.
Foldable t =>
(a -> a -> Ordering) -> t a -> a
F.maximumBy (((Int, r) -> r) -> (Int, r) -> (Int, r) -> Ordering
forall a b. Ord a => (b -> a) -> b -> b -> Ordering
comparing (Int, r) -> r
forall a b. (a, b) -> b
snd) (LSeq n (Int, r) -> (Int, r))
-> (LSeq n (Point d r :+ p) -> LSeq n (Int, r))
-> LSeq n (Point d r :+ p)
-> (Int, r)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Int -> (Point d r :+ p) -> (Int, r))
-> LSeq n (Point d r :+ p) -> LSeq n (Int, r)
forall a b (n :: Nat). (Int -> a -> b) -> LSeq n a -> LSeq n b
LSeq.mapWithIndex (\Int
i (Point d r
p :+ p
_) ->
                                                     (Int
i,Point d r -> LineSegment d p r -> r
forall (d :: Nat) r p.
(Arity d, Fractional r, Ord r) =>
Point d r -> LineSegment d p r -> r
sqDistanceToSeg Point d r
p LineSegment d p r
s)) (LSeq n (Point d r :+ p) -> (Int, r))
-> LSeq n (Point d r :+ p) -> (Int, r)
forall a b. (a -> b) -> a -> b
$ LSeq n (Point d r :+ p)
pts