module Geo.Computations.Trail
(
AvgMethod(..)
, Selected(..)
, PointGrouping
, TransformGrouping
, isSelected
, isNotSelected
, onSelected
, selLength
, totalDistance
, totalTime
, avgSpeeds
, slidingAverageSpeed
, closestDistance
, convexHull
, bezierCurveAt
, bezierCurve
, linearTime
, filterPoints
, betweenSpeeds
, restLocations
, spansTime
, everyNPoints
, intersectionOf
, invertSelection
, firstGrouping
, lastGrouping
, unionOf
, refineGrouping
, (/\), (\/)
, smoothRests
, smoothTrail
, bezierPoint
) where
import Text.Show.Functions ()
import Geo.Computations.Basic
import Geo.GPX.Conduit ()
import Control.Arrow (first)
import Control.Monad
import Data.Fixed (mod')
import Data.Function (on)
import Data.List as L
import Data.Maybe
import Data.Ord
import Data.Time
import Statistics.Function as F
import Statistics.Sample
import qualified Data.Vector.Unboxed as V
takeWhileEnd :: (a -> Bool) -> [a] -> (Maybe a, [a],[a])
takeWhileEnd p xs = go xs Nothing
where
go [] e = (e, [], [])
go (a:as) e
| p a = let (e',xs2,zs) = go as (Just a) in (e',a:xs2,zs)
| otherwise = (e,[], a:as)
data AvgMethod c
= AvgMean
| AvgHarmonicMean
| AvgGeometricMean
| AvgMedian
| AvgEndPoints
| AvgMinOf [AvgMethod c]
| AvgWith ([c] -> Speed)
avgSpeeds :: NominalDiffTime -> Trail Point -> [(UTCTime, Speed)]
avgSpeeds = slidingAverageSpeed AvgHarmonicMean
slidingAverageSpeed :: AvgMethod Point -> NominalDiffTime -> Trail Point -> [(UTCTime, Speed)]
slidingAverageSpeed _ _ [] = []
slidingAverageSpeed m minTime xs =
let pts = map unSelect (spansTime minTime xs)
spds = map (getAvg m) pts
times = map getAvgTimes pts
in concatMap maybeToList $ zipWith (\t s -> fmap (,s) t) times spds
where
getTimeDiff a b = on (liftM2 diffUTCTime) pntTime a b
getAvg _ [] = 0
getAvg _ [_] = 0
getAvg m2 cs =
let ss = getSpeedsV cs
in case m2 of
AvgMean -> mean ss
AvgHarmonicMean -> harmonicMean ss
AvgGeometricMean -> geometricMean ss
AvgMedian ->
let ss' = F.sort $ getSpeedsV cs
len = V.length ss'
mid = len `div` 2
in if V.length ss' < 3
then mean ss'
else if odd len then ss' V.! mid else mean (V.slice mid 2 ss')
AvgEndPoints -> fromMaybe 0 $ speed (head cs) (last cs)
AvgMinOf as -> minimum $ map (flip getAvg cs) as
AvgWith f -> f cs
getAvgTimes [] = Nothing
getAvgTimes [x] = pntTime x
getAvgTimes ps = getAvgTime (head ps) (last ps)
getAvgTime a b = liftM2 addUTCTime (getTimeDiff b a) (pntTime a)
getSpeedsV = V.fromList . getSpeeds
getSpeeds zs = concatMap maybeToList $ zipWith speed zs (drop 1 zs)
type PointGrouping c = Trail c -> [Selected (Trail c)]
type TransformGrouping c = [Selected (Trail c)] -> [Selected (Trail c)]
data Selected a = Select {unSelect :: a} | NotSelect {unSelect :: a}
deriving (Eq, Ord, Show)
isSelected :: Selected a -> Bool
isSelected (Select _) = True
isSelected _ = False
isNotSelected :: Selected a -> Bool
isNotSelected = not . isSelected
selLength :: Selected [a] -> Int
selLength = length . unSelect
onSelected :: (a -> b) -> (a -> b) -> Selected a -> b
onSelected f _ (Select a) = f a
onSelected _ g (NotSelect a) = g a
instance Functor Selected where
fmap f (Select x) = Select $ f x
fmap f (NotSelect x) = NotSelect $ f x
dropExact :: Int -> [Selected [a]] -> [Selected [a]]
dropExact _ [] = []
dropExact i (x:xs) =
case compare (selLength x) i of
EQ -> xs
LT -> dropExact (i selLength x) xs
GT -> fmap (drop i) x : xs
betweenSpeeds :: Double -> Double -> PointGrouping Point
betweenSpeeds low hi ps =
let spds = concatMap maybeToList $ zipWith speed ps (drop 1 ps)
psSpds = [(p,s) | p <- ps, s <- maybeToList (listToMaybe spds) ++ spds]
inRange x = x >= low && x <= hi
chunk [] = []
chunk xs@(x:_) =
let op p = if inRange (snd x) then first Select . span p else first NotSelect . break p
(r,rest) = op (inRange . snd) xs
in r : chunk rest
in map (fmap (map fst)) $ chunk psSpds
restLocations :: Distance -> NominalDiffTime -> PointGrouping Point
restLocations d s ps =
let consToFirst x [] = [NotSelect [x]]
consToFirst x (a:as) = (fmap (x:) a) : as
go [] [] = []
go [] nonRests = [NotSelect $ reverse nonRests]
go (a:as) nonRests =
case takeWhileEnd ((<=) d . distance a) as of
(Just l, close, far) ->
case (pntTime a, pntTime l) of
(Just t1, Just t2) ->
let diff = diffUTCTime t2 t1
in if diff >= s then NotSelect (reverse nonRests) : Select (a:close) : go far [] else go as (a:nonRests)
_ -> consToFirst a $ go as nonRests
_ -> consToFirst a $ go as nonRests
in go ps []
spansTime :: NominalDiffTime -> PointGrouping Point
spansTime n ps =
let times = mkTimePair ps
chunk [] = []
chunk xs@(x:_) =
let (good,rest) = span ((<= addUTCTime n (snd x)) . snd) xs
in if null good then [xs] else good : chunk rest
in map (Select . map fst) $ chunk times
intersectionOf :: [PointGrouping Point] -> PointGrouping Point
intersectionOf gs ps =
let groupings = map ($ ps) gs
chunk _ [] = []
chunk ggs xs =
let minLen = max 1 . minimum . concatMap (take 1) $ map (map selLength) ggs
sel = if all isSelected (concatMap (take 1) ggs) then Select else NotSelect
(c,rest) = splitAt minLen xs
in sel c : chunk (filter (not . null) $ map (dropExact minLen) ggs) rest
in chunk groupings ps
unionOf :: [PointGrouping Point] -> PointGrouping Point
unionOf gs ps =
let groupings = map ($ ps) gs
chunk _ [] = []
chunk ggs xs =
let getSegs = concatMap (take 1)
segs = getSegs ggs
len =
if any isSelected segs
then max 1 . maximum . getSegs . map (map selLength) . map (filter isSelected) $ ggs
else max 1 . minimum . getSegs . map (map selLength) $ ggs
sel = if any isSelected segs then Select else NotSelect
(c,rest) = splitAt len xs
in sel c : chunk (filter (not . null) $ map (dropExact len) ggs) rest
in chunk groupings ps
(/\) :: [Selected (Trail a)] -> TransformGrouping a
(/\) _ [] = []
(/\) [] _ = []
(/\) xsL@(Select x:_) ysL@(Select y:_) =
let z = if length x < length y then x else y
xs' = selListDrop (length z) xsL
ys' = selListDrop (length z) ysL
in Select z : (xs' /\ ys')
(/\) xs (NotSelect y:ys) = NotSelect y : (selListDrop (length y) xs /\ ys)
(/\) (NotSelect x:xs) ys = NotSelect x : (xs /\ selListDrop (length x) ys)
(\/) :: [Selected (Trail a)] -> TransformGrouping a
(\/) xs [] = xs
(\/) [] ys = ys
(\/) (Select x:xs) (Select y : ys) =
let xLen = length x
yLen = length y
in if xLen < yLen
then (Select y :) (selListDrop (yLen xLen) xs \/ ys)
else (Select x :) (xs \/ selListDrop (xLen yLen) ys)
(\/) (Select x:_) ys = Select x : selListDrop (length x) ys
(\/) xs (Select y:_) = Select y : selListDrop (length y) xs
(\/) xsL@(NotSelect x:xs) ysL@(NotSelect y:ys) =
let xLen = length x
yLen = length y
in if xLen < yLen
then (NotSelect x:) (xs \/ selListDrop xLen ysL)
else (NotSelect y:) (selListDrop yLen xsL \/ ys)
selListDrop :: Int -> [Selected [a]] -> [Selected [a]]
selListDrop 0 xs = xs
selListDrop _ [] = []
selListDrop n (x:xs) =
let x' = drop n (unSelect x)
in fmap (const x') x : selListDrop (n (selLength x length x')) xs
invertSelection :: TransformGrouping a
invertSelection = map (onSelected NotSelect Select)
firstGrouping :: TransformGrouping a
firstGrouping ps = take 1 ps ++ map (NotSelect . unSelect) (drop 1 ps)
lastGrouping :: TransformGrouping a
lastGrouping ps = let ps' = reverse ps in reverse $ take 1 ps' ++ map (NotSelect . unSelect) (drop 1 ps')
everyNPoints :: Int -> PointGrouping a
everyNPoints n ps
| n <= 0 = [NotSelect ps]
| otherwise = go ps
where
go [] = []
go xs = let (h,t) = splitAt n xs in Select h : go t
refineGrouping :: PointGrouping a -> TransformGrouping a
refineGrouping b = concatMap (onSelected b (\x -> [NotSelect x]))
filterPoints :: PointGrouping a -> Trail a -> Trail a
filterPoints g = concatMap unSelect . filter isSelected . g
mkTimePair :: Trail Point -> [(Point,UTCTime)]
mkTimePair xs =
let timesM = map (\x-> fmap (x,) $ pntTime x) xs
in concatMap maybeToList timesM
bezierCurveAt :: [UTCTime] -> Trail Point -> Trail Point
bezierCurveAt _ [] = []
bezierCurveAt selectedTimes xs =
let timesDef = mkTimePair xs
end = last timesDef
top = head timesDef
tTime = diffUTCTime (snd end) (snd top)
times = if null selectedTimes then map snd timesDef else selectedTimes
diffTimes = [diffUTCTime t (snd top) / tTime | t <- times]
queryTimes = map realToFrac diffTimes
in if tTime <= 0 || any (\x -> x < 0 || x > 1) queryTimes
then xs
else
if null timesDef || any (\x -> x < 0 || x > 1) queryTimes
then xs
else let curvePoints = (map (bezierPoint xs) queryTimes)
newTimes = [addUTCTime t (snd top) | t <- diffTimes]
in zipWith (\t p -> p { pntTime = Just t}) newTimes curvePoints
bezierPoint :: [Point] -> Double -> Point
bezierPoint pnts t = go pnts
where
go [] = error "GPS Package: Can not create a bezier point from an empty list"
go [p] = p
go ps = interpolate (go (init ps)) (go (tail ps)) t
bezierCurve :: [Selected (Trail Point)] -> Trail Point
bezierCurve = concatMap (onSelected (bezierCurveAt []) Prelude.id)
linearTime :: [Point] -> [Point]
linearTime [] = []
linearTime (p:ps) = go (pntTime p) ps
where
go _ [] = []
go t (x:xs) = if pntTime x < t then go t xs else x : go (pntTime x) xs
closestDistance :: Trail Point -> Trail Point -> Maybe Distance
closestDistance as bs = listToMaybe $ L.sort [distance a b | a <- as, b <- bs]
totalDistance :: [Point] -> Distance
totalDistance as = sum $ zipWith distance as (drop 1 as)
totalTime :: Trail Point -> NominalDiffTime
totalTime [] = 0
totalTime xs@(x:_) = fromMaybe 0 $ liftM2 diffUTCTime (pntTime x) (pntTime $ last xs)
convexHull :: [Point] -> [Point]
convexHull lst =
let frst = southMost lst
in case frst of
Nothing -> []
Just f ->
let sorted = L.sortBy (comparing (eastZeroHeading f)) (filter (/= f) lst)
in case sorted of
(a:b:cs) -> grahamScan (b:a:f:[]) cs
cs -> f : cs
where
grahamScan [] _ = []
grahamScan ps [] = ps
grahamScan (x:[]) _ = [x]
grahamScan (p2:p1:ps) (x:xs) =
case turn p1 p2 x of
LeftTurn -> grahamScan (x:p2:p1:ps) xs
Straight -> grahamScan (x:p2:p1:ps) xs
_ -> grahamScan (p1:ps) (x:xs)
eastZeroHeading :: Point -> Point -> Heading
eastZeroHeading s = (`mod'` (2*pi)) . (+ pi/2) . heading s
data Turn = LeftTurn | RightTurn | Straight deriving (Eq, Ord, Show, Read, Enum)
turn :: Point -> Point -> Point -> Turn
turn a b c =
let h1 = eastZeroHeading a b
h2 = eastZeroHeading b c
d = h2 h1
in if d >= 0 && d < pi then LeftTurn else RightTurn
southMost :: [Point] -> Maybe Point
southMost [] = Nothing
southMost cs = Just . minimumBy (comparing pntLat) $ cs
smoothRests :: Trail Point -> Trail Point
smoothRests = bezierCurve . refineGrouping (everyNPoints 8) . restLocations 30 60
smoothTrail :: Trail Point -> Trail Point
smoothTrail = gSmoothSome 7
gSmoothSome :: Int -> Trail Point -> Trail Point
gSmoothSome n = bezierCurve . everyNPoints n