hgeometry-0.13: src/Algorithms/Geometry/PolyLineSimplification/ImaiIri.hs

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

import           Algorithms.Graph.BFS (bfs')
import           Control.Lens
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.List.NonEmpty (NonEmpty(..))
import qualified Data.List.NonEmpty as NonEmpty
import qualified Data.Sequence as Seq
import           Data.Tree
import qualified Data.Vector as V
import           Witherable

-- import Data.RealNumber.Rational
-- type R = RealNumber 5

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

-- | Line simplification with the Imai-Iri alogrithm. Given a distance
-- value eps and 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) \) time.
simplify     :: (Ord r, Fractional r, Arity d)
             => r -> PolyLine d p r -> PolyLine d p r
simplify eps = simplifyWith $ \shortcut subPoly -> all (closeTo shortcut) (subPoly^.points)
  where
    closeTo seg (p :+ _) = sqDistanceToSeg p seg  <= epsSq
    epsSq = eps*eps

-- | Given a function that tests if the shortcut is valid, compute a
-- simplification using the Imai-Iri algorithm.
--
-- Running time: \( O(Tn^2 \) time, where \(T\) is the time to
-- evaluate the predicate.
simplifyWith            :: (LineSegment d p r -> PolyLine  d p r -> Bool)
                        -> PolyLine d p r -> PolyLine d p r
simplifyWith isValid pl = pl&points %~ (LSeq.promise @2 . extract path)
  where
    g    = mkGraph isValid pl
    spt  = bfs' 0 g
    path = case pathsTo (pl^.points.to F.length - 1) spt of
             []      -> error "no path found?"
             (pth:_) -> pth

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

type Graph = V.Vector [Int]

-- | Constructs the shortcut graph
mkGraph         :: (LineSegment d p r -> PolyLine d p r -> Bool) -> PolyLine d p r -> Graph
mkGraph isValid = flip V.snoc [] . V.imap f . V.fromList . F.toList . allPrefixes
  where
    f i subPl = catMaybes
              $ zipWith isValid' [i+1..] . F.toList . allSuffixes $ subPl

    isValid' j subPoly = let shortcut = ClosedLineSegment (subPoly^.start) (subPoly^.end)
                         in if isValid shortcut subPoly then Just j else Nothing

-- | Generates all prefixes of the polyline; i.e. all contiguous
-- polylines all starting at the original starting point.
allPrefixes    :: PolyLine d p r -> Seq.Seq (PolyLine d p r)
allPrefixes pl = mapMaybe mkPolyLine . Seq.tails . LSeq.toSeq $ pl^.points

mkPolyLine :: Seq.Seq (Point d r :+ p) -> Maybe (PolyLine d p r)
mkPolyLine = fmap PolyLine . LSeq.eval @2 . LSeq.fromSeq

-- | Generates all suffixes of the polyline.
allSuffixes :: PolyLine d p r -> Seq.Seq (PolyLine d p r)
allSuffixes pl = mapMaybe mkPolyLine . Seq.drop 2 . Seq.inits . LSeq.toSeq $ pl^.points






-- | Get all paths to the particular element in the tree.
pathsTo   :: Eq a => a -> Tree a -> [NonEmpty a]
pathsTo x = findPaths (== x)

-- | All paths to the nodes satisfying the predicate.
findPaths   :: (a -> Bool) -> Tree a -> [NonEmpty a]
findPaths p = go
  where
    go (Node x chs) = case foldMap go chs of
                        []    | p x       -> [x:|[]]
                              | otherwise -> []
                        paths | p x       -> (x:|[]) : map (x NonEmpty.<|) paths
                              | otherwise ->           map (x NonEmpty.<|) paths




-- | Given a non-empty list of indices, and some LSeq, extract the elemnets
-- on those indices.
--
-- running time: \(O(n)\)
extract    :: NonEmpty Int -> LSeq.LSeq n a -> LSeq.LSeq 0 a
extract is = LSeq.fromList . extract' (F.toList is) 0 . F.toList

extract'                                 :: [Int] -> Int -> [a] -> [a]
extract' []         _ _                  = []
extract' (_:_)      _ []                 = []
extract' is'@(i:is) j (x:xs) | i == j    = x : extract' is (j+1) xs
                             | otherwise = extract' is' (j+1) xs

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


-- tr :: Tree Int
-- tr = Node 0 [Node 1 [], Node 2 [Node 3 [], Node 2 [], Node 4 [Node 5 []]]]

-- poly :: PolyLine 2 Int R
-- poly = case fromPoints [origin :+ 0, Point2 1 1 :+ 1, Point2 2 2 :+ 2, Point2 3 3 :+ 3] of
--          Just p -> p

-- test = Seq.fromList [0..5]

-- myTree :: Tree Int
-- myTree = Node {rootLabel = 0, subForest = [Node {rootLabel = 1, subForest = []}
--                                        ,Node {rootLabel = 2, subForest = []}
--                                        ,Node {rootLabel = 3, subForest = []}]
--            }