{-# LANGUAGE TupleSections #-} module Geo.Computations.Trail ( -- * Types AvgMethod(..) , Selected(..) , PointGrouping , TransformGrouping -- * Utility Functions , isSelected , isNotSelected , onSelected , selLength -- * Trail Functions -- ** Queries , totalDistance , totalTime , avgSpeeds , slidingAverageSpeed , closestDistance , convexHull -- ** Transformations , bezierCurveAt , bezierCurve , linearTime , filterPoints -- ** Grouping Methods , betweenSpeeds , restLocations , spansTime , everyNPoints -- ** Group Transformations , intersectionOf , invertSelection , firstGrouping , lastGrouping , unionOf , refineGrouping , (/\), (\/) -- ** Composite Operations (Higher Level) , smoothRests , smoothTrail -- * Misc , 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 -- ^ Obtain the 'mean' of the considered points | AvgHarmonicMean -- ^ Obtain the 'harmonicMean' | AvgGeometricMean -- ^ Obtain the 'geometricMean' | AvgMedian -- ^ Obtain the median of the considered points | AvgEndPoints -- ^ Compute the speed considering only the given endpoints | AvgMinOf [AvgMethod c] -- ^ Take the minimum of the speeds from the given methods | AvgWith ([c] -> Speed) -- | @avgSpeeds n points@ -- Average speed using a window of up to @n@ seconds and averaging by taking the -- Median ('AvgMedian'). avgSpeeds :: NominalDiffTime -> Trail Point -> [(UTCTime, Speed)] avgSpeeds = slidingAverageSpeed AvgHarmonicMean -- | @slidingAverageSpeed m n@ Average speed using a moving window of up to @n@ seconds -- and an 'AvgMethod' of @m@. 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 :: [] -> AvgMethod -> Speed 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) -- | A PointGrouping is a function that selects segments of a trail. -- -- Grouping point _does not_ result in deleted points. It is always true that: -- -- forall g :: PointGrouping c --> -- concatMap unSelect (g ts) == ts -- -- The purpose of grouping is usually for later processing. Any desire to drop -- points that didn't meet a particular grouping criteria can be filled with -- a composition with 'filter' (or directly via 'filterPoints'). type PointGrouping c = Trail c -> [Selected (Trail c)] -- | Given a selection of coordinates, transform the selected -- coordinates in some way (while leaving the non-selected -- coordinates unaffected). type TransformGrouping c = [Selected (Trail c)] -> [Selected (Trail c)] -- | When grouping points, lists of points are either marked as 'Select' or 'NotSelect'. 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 -- | Groups trail segments into contiguous points within the speed -- and all others outside of the speed. The "speed" from point p(i) -- to p(i+1) is associated with p(i) (execpt for the first speed -- value, which is associated with both the first and second point) 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 -- | A "rest point" means the coordinates remain within a given distance -- for at least a particular amount of time. 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 [] -- |chunking points into groups spanning at most the given time -- interval. 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 -- | intersects the given groupings intersectionOf :: [PointGrouping Point] -> PointGrouping Point intersectionOf gs ps = let groupings = map ($ ps) gs -- chunk :: [[Selected [pnts]]] -> pnts -> [pnts] chunk _ [] = [] chunk ggs xs = let minLen = max 1 . minimum . concatMap (take 1) $ map (map selLength) ggs -- FIXME this is all manner of broken 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 -- | Union all the groupings 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 -- | Intersection binary operator (/\) :: [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) -- | Union binary operator (\/) :: [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 -- |Inverts the selected/nonselected segments invertSelection :: TransformGrouping a invertSelection = map (onSelected NotSelect Select) -- |@firstGrouping f ps@ only the first segment remains 'Select'ed, and only -- if it was already selected by @f@. firstGrouping :: TransformGrouping a firstGrouping ps = take 1 ps ++ map (NotSelect . unSelect) (drop 1 ps) -- | Only the last segment, if any, is selected (note: the current -- implementation is inefficient, using 'reverse') lastGrouping :: TransformGrouping a lastGrouping ps = let ps' = reverse ps in reverse $ take 1 ps' ++ map (NotSelect . unSelect) (drop 1 ps') -- | chunk the trail into groups of N points 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 -- |For every selected group, refine the selection using the second -- grouping method. This differs from 'IntersectionOf' by restarting -- the second grouping algorithm at the beginning each group selected -- by the first algorithm. refineGrouping :: PointGrouping a -> TransformGrouping a refineGrouping b = concatMap (onSelected b (\x -> [NotSelect x])) -- |Remove all points that remain 'NotSelect'ed by the given grouping algorithm. filterPoints :: PointGrouping a -> Trail a -> Trail a filterPoints g = concatMap unSelect . filter isSelected . g -- Extract the time from each coordinate. If no time is available then -- the coordinate is dropped! mkTimePair :: Trail Point -> [(Point,UTCTime)] mkTimePair xs = let timesM = map (\x-> fmap (x,) $ pntTime x) xs in concatMap maybeToList timesM -- |Construct a bezier curve using the provided trail. Construct a -- new trail by sampling the given bezier curve at the given times. -- The current implementation assumes the times of the input -- coordinates are available and all equal (Ex: all points are 5 -- seconds apart), the results will be poor if this is not the case! 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 -- error "bezierCurveAt has a out-of-bound time!" 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 -- |Interpolate selected points onto a bezier curve. Note this gets -- exponentially more expensive with the length of the segement being -- transformed - it is not advisable to perform this operation on -- trail segements with more than ten points! bezierCurve :: [Selected (Trail Point)] -> Trail Point bezierCurve = concatMap (onSelected (bezierCurveAt []) Prelude.id) -- |Filters out any points that go backward in time (thus must not be -- valid if this is a trail) 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 -- |Returns the closest distance between two trails (or Nothing if a -- trail is empty). Inefficient implementation: -- O( (n * m) * log (n * m) ) closestDistance :: Trail Point -> Trail Point -> Maybe Distance closestDistance as bs = listToMaybe $ L.sort [distance a b | a <- as, b <- bs] -- | Find the total distance traveled 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) -- | Uses Grahams scan to compute the convex hull of the given points. -- This operation requires sorting of the points, so don't try it unless -- you have notably more memory than the list of points will consume. 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 -- | Find the southmost point southMost :: [Point] -> Maybe Point southMost [] = Nothing southMost cs = Just . minimumBy (comparing pntLat) $ cs ---------- COMPOSIT OPERATIONS --------------- -- These operations are simply implemented using the previously -- defined functions. They can serve either for concise use for novice -- users or as instructional examples. ------------------------------------------ -- | Smooth points with rest areas using a bezierCurve. -- -- Parameters: rest for 1 minute within 30 meters get smoothed -- in a bezier curve over every 8 points. smoothRests :: Trail Point -> Trail Point smoothRests = bezierCurve . refineGrouping (everyNPoints 8) . restLocations 30 60 -- |Smooth every 7 points using a bezier curve smoothTrail :: Trail Point -> Trail Point smoothTrail = gSmoothSome 7 -- |Smooth every n points using a bezier curve gSmoothSome :: Int -> Trail Point -> Trail Point gSmoothSome n = bezierCurve . everyNPoints n