diff --git a/Data/GPS.hs b/Data/GPS.hs
--- a/Data/GPS.hs
+++ b/Data/GPS.hs
@@ -1,273 +1,12 @@
-{-# LANGUAGE MultiParamTypeClasses, FunctionalDependencies, FlexibleInstances, EmptyDataDecls, BangPatterns #-}
+{-# LANGUAGE MultiParamTypeClasses, FunctionalDependencies, FlexibleInstances, EmptyDataDecls, BangPatterns, TupleSections #-}
 -- |A basic GPS library with calculations for distance and speed along
 -- with helper functions for filtering/smoothing trails.  All distances
 -- are in meters and time is in seconds.  Speed is thus meters/second
 
 module Data.GPS
-       ( -- * Types
-         Distance
-       , Heading
-       , Speed
-       , Vector
-       , Trail
-         -- * Constants
-       , north
-       , south
-       , east
-       , west
-       , radiusOfEarth
-         -- * Coordinate Functions
-       , heading
-       , distance
-       , speed
-       , addVector
-       , getRadianPair
-       , getDMSPair
-       , divideArea
-         -- * Trail Functions
-       , totalDistance
-       , restLocations
-       , closestDistance
-       , filterByMaxSpeed
-       , convexHull
-         -- * Other helpers
-       , readGPX
-       , module Data.Geo.GPX
+       ( module Data.GPS.Core
+       , module Data.GPS.Trail
        ) where
 
-import Data.Function (on)
-import Data.Ord (comparing)
-import Data.List (sort, mapAccumL, minimumBy, maximumBy, sortBy)
-import Data.Geo.GPX hiding (none, cmt)
-
-import Text.XML.HXT.Arrow
-import Text.XML.XSD.DateTime(DateTime, toUTCTime)
-
-import Data.Time
-import Data.Maybe (listToMaybe)
-import Data.Fixed
-import Control.Monad
-
--- |Distances are expressed in meters
-type Distance = Double
-
--- |Angles are expressed in radians from North.
--- 	0	== North
--- 	pi/2	== West
--- 	pi 	== South
--- 	(3/2)pi	== East    == - (pi / 2)
-type Heading = Double
-
--- |Speed is hard coded as meters per second
-type Speed = Double
-
-type Vector = (Distance, Heading)
-type Trail a = [a]
-
-getUTCTime :: (Lat a, Lon a, Time a) => a -> Maybe UTCTime
-getUTCTime = fmap toUTCTime . time
-
-distance :: (Lat a, Lon a, Lat b, Lon b) => a -> b -> Distance
-distance a b =
-	let x  = sin lat1 * sin lat2 + cos lat1 * cos lat2 * cos (lon2 - lon1)
-	    x' = if x > 1 then 1 else x
-	in radiusOfEarth * acos x'
- where
-  (lat1, lon1) = getRadianPairD a
-  (lat2, lon2) = getRadianPairD b
-
--- | Find the total distance traveled
-totalDistance :: (Lat a, Lon a) => [a] -> Distance
-totalDistance as = sum $ zipWith distance as (drop 1 as)
-
--- | Direction two points aim toward (0 = North, pi/2 = West, pi = South, 3pi/2 = East)
-heading         :: (Lat a, Lon a, Lat b, Lon b) => a -> b -> Heading	-- ^ 0 = North, pi/2 = West...
-heading a b =
-	atan2	(sin (diffLon) * cos (lat2)) 
-		(cos(lat1) * sin (lat2) - sin(lat1) * cos lat2 * cos (diffLon))
- where
-  (lat1, lon1) = getRadianPairD a
-  (lat2, lon2) = getRadianPairD b
-  diffLon = lon1 - lon2
-
-getVector :: (Lat a, Lon a, Lat b, Lon b) => a -> b -> Vector
-getVector a b = (distance a b, heading a b)
-
--- | Speed in meters per second, only if a 'Time' was recorded for each waypoint.
-speed :: (Lat loc, Lon loc, Time loc, Lat b, Lon b, Time b) => loc -> b -> Maybe Speed
-speed a b = 
-  case (getUTCTime b, getUTCTime a) of
-    (Just x, Just y) -> 
-      let timeDiff = realToFrac (diffUTCTime x y)
-      in if timeDiff == 0 then Nothing else Just $ (distance a b) / timeDiff
-    _ -> Nothing
-
--- |Filter out all points that result in a speed greater than a given
--- value (the second point is dropped)
-filterByMaxSpeed :: (Lat loc, Lon loc, Time loc) => Speed -> Trail loc -> Trail loc
-filterByMaxSpeed mx xs =
-	let ss = zipWith speed xs (drop 1 xs)
-	    fs = filter ((< Just mx) . fst) (zip ss $ drop 1 xs)
-	in take 1 xs ++ map snd fs
-
-data TempTrail a = T (Trail a) a
-
--- |radius of the earth in meters
-radiusOfEarth :: Double
-radiusOfEarth = 6378700
-
--- |North is 0 radians
-north :: Heading
-north = 0
-
--- |South, being 180 degrees from North, is pi.
-south :: Heading
-south = pi
-
--- |East is 270 degrees from North
-east :: Heading
-east = (3 / 2) * pi
-
--- |West is 90 degrees (pi/2)
-west :: Heading
-west = pi / 2
-
-toDecimal = (*) (180 / pi)
-
-
--- |Given a vector and coordinate, computes a new coordinate.
--- Within some epsilon it should hold that if
---
--- 	@dest = addVector (dist,heading) start@
---
--- then
---
--- 	@heading == dmsHeading start dest@
--- 	
--- 	@dist    == distance start dest@
-addVector :: (Lat c, Lon c) => Vector -> c -> c
-addVector (d,h) p = setLon (longitudeType lon2) . setLat (latitudeType lat2) $ p
-  where
-	(lat,lon) = getRadianPairD p
-	lat2 = lat + (cos h) * (d / radiusOfEarth)
-	lon2 = lon + acos ( (cos (d/radiusOfEarth) - sin lat * sin lat2) / (cos lat * cos lat2))
-
-getRadianPairD :: (Lat c, Lon c) => c -> (Double,Double)
-getRadianPairD = (\(a,b) -> (realToFrac a, realToFrac b)) . getRadianPair
-
-getDMSPair :: (Lat c, Lon c) => c -> (LatitudeType, LongitudeType)
-getDMSPair c = (lat c, lon c)
-
--- |Provides a lat/lon pair of doubles in radians
-getRadianPair :: (Lat p, Lon p) => p -> (LatitudeType, LongitudeType)
-getRadianPair p = (toRadians (lat p), toRadians (lon p))
-
-toRadians :: Floating f => f -> f
-toRadians = (*) (pi / 180)
-
--- |Filters out any points that go backward in time (thus must not be valid if this is a trail)
-linearTime :: (Lat a, Lon a, Time a) => Trail a -> Trail a
-linearTime [] = []
-linearTime (p:ps) = go (getUTCTime p) ps
-  where
-  go _ [] = []
-  go t (p:ps) = if getUTCTime p < t then go t ps else p : go (getUTCTime p) ps
-
--- |Creates a list of trails all of which are within the given distance of each
--- other spanning atleast the given amount of time.
---
--- For example @restLocations 50 600@
--- would return lists of all points that are within 50 meters of each other and
--- span at least 10 minutes (600 seconds).
---
--- Note this gives points within fifty meters of the earliest point - wandering
--- in a rest area with a 50 meter radius could result in several rest points
--- ([a,b..]) or even none if the distance between individual points exceeds 50m.
-restLocations :: (Lat a, Lon a, Time a) => Distance -> NominalDiffTime -> Trail a -> [Trail a]
-restLocations d s xs = go xs
-  where
-  go [] = []
-  go (a:as) =
-	case takeWhileLast ((<=) d . distance a) as of
-		(Just l, close, far) ->
-			case (getUTCTime a, getUTCTime l) of
-				(Just t1, Just t2) ->
-					let diff = diffUTCTime t2 t1
-					in if diff >= s then (a:close) : go far else go as
-				_ -> go as
-		_ -> go as
-
-takeWhileLast :: (a -> Bool) -> [a] -> (Maybe a, [a], [a])
-takeWhileLast p [] = (Nothing, [], [])
-takeWhileLast p (x:xs)
-	| not (p x) = (Nothing, [], x:xs)
-	| otherwise = go x xs
-  where
-  go !a [] = (Just a, [], [])
-  go !a (b:bs)
-	| p b = let (c,d,f) = go b bs in (c, a:d, f)
-	| otherwise = (Just a, [], b:bs)
-
--- |Returns the closest distance between two trails (or Nothing if a trail is empty)
--- O( (n * m) * log (n * m) )
-closestDistance :: (Lat a, Lon a) => Trail a -> Trail a -> Maybe Distance
-closestDistance as bs = listToMaybe $ sort [distance a b | a <- as, b <- bs]
-
-
--- |@divideArea vDist hDist nw se@ divides an area into a grid of equally
--- spaced coordinates within the box drawn by the northwest point (nw) and
--- southeast point (se).  Because this uses floating point there might be a
--- different number of points in some rows (the last might be too far east based
--- on a heading from the se point).
-divideArea :: (Lat c, Lon c) => Distance -> Distance -> c -> c -> [[c]]
-divideArea vDist hDist nw se =
-	let (top,left)  = (lat nw, lon nw)
-	    (btm,right) = (lat se, lon se)
-	    columnOne = takeWhile ( (<= west) . heading se) . iterate (addVector (vDist, south)) $ nw
-	    buildRow  = takeWhile ((>= north) . heading se) . iterate (addVector (hDist, east))
-	in map buildRow columnOne
-
--- |Reads a GPX file (using the GPX library) by simply concatenating all the
--- tracks, segments, and points ('trkpts', 'trksegs', 'trks') into a single 'Trail'.
-readGPX :: FilePath -> IO (Trail WptType)
-readGPX = liftM (concatMap trkpts . concatMap trksegs . concatMap trks) . readGpxFile
-
--- | 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 :: (Eq c, Lat c, Lon c) => [c] -> [c]
-convexHull xs =
-	let first = southMost xs
-	in case first of
-		Nothing -> []
-		Just f  ->
-	    	     let sorted = sortBy (comparing (eastZeroHeading f)) (filter (/= f) xs)
-		     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 :: (Lat c, Lon c) => c -> c -> Heading
-eastZeroHeading s = (`mod'` (2*pi)) . (+ pi/2) . heading s
-
-data Turn = LeftTurn | RightTurn | Straight deriving (Eq, Ord, Show, Read, Enum)
-
-turn :: (Lat c, Lon c) => c -> c -> c -> 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 :: (Lat c) => [c] -> Maybe c
-southMost []  = Nothing
-southMost cs = Just . minimumBy (comparing lat) $ cs
+import Data.GPS.Core
+import Data.GPS.Trail
diff --git a/Data/GPS/Core.hs b/Data/GPS/Core.hs
new file mode 100644
--- /dev/null
+++ b/Data/GPS/Core.hs
@@ -0,0 +1,183 @@
+module Data.GPS.Core
+       ( -- * Types
+         Distance
+       , Heading
+       , Speed
+       , Vector
+       , Trail
+         -- * Constants
+       , north
+       , south
+       , east
+       , west
+       , radiusOfEarth
+         -- * Coordinate Functions
+       , heading
+       , distance
+       , speed
+       , getVector
+       , addVector
+       , getRadianPair
+       , getDMSPair
+       , divideArea
+       , interpolate
+         -- * IO helpers
+       , writeGPX
+       , readGPX
+       , readGPXSegments
+         -- * Utility
+       , getUTCTime
+       , module Data.Geo.GPX
+         ) where
+
+import Data.Time
+import Data.Maybe
+import Control.Monad
+import Text.XML.HXT.Arrow
+import Text.XML.XSD.DateTime(DateTime,toUTCTime)
+import Data.Geo.GPX
+
+-- |Distances are expressed in meters
+type Distance = Double
+
+-- |Angles are expressed in radians from North.
+-- 	0	== North
+-- 	pi/2	== West
+-- 	pi 	== South
+-- 	(3/2)pi	== East    == - (pi / 2)
+type Heading = Double
+
+-- |Speed is hard coded as meters per second
+type Speed = Double
+type Vector = (Distance, Heading)
+
+type Trail a = [a]
+
+getUTCTime :: (Time a) => a -> Maybe UTCTime
+getUTCTime = fmap toUTCTime . time
+
+acos' x = if x > 1 then acos 1 else if x < (-1) then acos (-1) else acos x
+
+distance :: (Lat a, Lon a, Lat b, Lon b) => a -> b -> Distance
+distance a b =
+	let x  = sin lat1 * sin lat2 + cos lat1 * cos lat2 * cos (lon2 - lon1)
+	in radiusOfEarth * acos' x
+ where
+  (lat1, lon1) = getRadianPairD a
+  (lat2, lon2) = getRadianPairD b
+
+-- | Direction two points aim toward (0 = North, pi/2 = West, pi = South, 3pi/2 = East)
+heading         :: (Lat a, Lon a, Lat b, Lon b) => a -> b -> Heading
+heading a b =
+	atan2	(sin (diffLon) * cos (lat2)) 
+		(cos(lat1) * sin (lat2) - sin(lat1) * cos lat2 * cos (diffLon))
+ where
+  (lat1, lon1) = getRadianPairD a
+  (lat2, lon2) = getRadianPairD b
+  diffLon = lon1 - lon2
+
+getVector :: (Lat a, Lon a, Lat b, Lon b) => a -> b -> Vector
+getVector a b = (distance a b, heading a b)
+
+-- |Given a vector and coordinate, computes a new coordinate.
+-- Within some epsilon it should hold that if
+--
+-- 	@dest = addVector (dist,heading) start@
+--
+-- then
+--
+-- 	@heading == heading start dest@
+-- 	
+-- 	@dist    == distance start dest@
+addVector :: (Lat c, Lon c) => Vector -> c -> c
+addVector (d,h) p = setLon (longitudeType $ toDegrees lon2) 
+                  . setLat (latitudeType $ toDegrees lat2) $ p
+  where
+	(lat,lon) = getRadianPairD p
+	lat2 = lat + (cos h) * (d / radiusOfEarth)
+	lon2 = lon - acos' ( (cos (d/radiusOfEarth) - sin lat * sin lat2) / (cos lat * cos lat2))
+
+-- | Speed in meters per second, only if a 'Time' was recorded for each waypoint.
+speed :: (Lat loc, Lon loc, Time loc, Lat b, Lon b, Time b) => loc -> b -> Maybe Speed
+speed a b = 
+  case (getUTCTime b, getUTCTime a) of
+    (Just x, Just y) -> 
+      let timeDiff = realToFrac (diffUTCTime x y)
+      in if timeDiff == 0 then Nothing else Just $ (distance a b) / timeDiff
+    _ -> Nothing
+
+-- |radius of the earth in meters
+radiusOfEarth :: Double
+radiusOfEarth = 6378700
+
+-- |North is 0 radians
+north :: Heading
+north = 0
+
+-- |South, being 180 degrees from North, is pi.
+south :: Heading
+south = pi
+
+-- |East is 270 degrees (3 pi / 2)
+east :: Heading
+east = (3 / 2) * pi
+
+-- |West is 90 degrees (pi/2)
+west :: Heading
+west = pi / 2
+
+toDegrees = (*) (180 / pi)
+
+getRadianPairD :: (Lat c, Lon c) => c -> (Double,Double)
+getRadianPairD = (\(a,b) -> (realToFrac a, realToFrac b)) . getRadianPair
+
+getDMSPair :: (Lat c, Lon c) => c -> (LatitudeType, LongitudeType)
+getDMSPair c = (lat c, lon c)
+
+-- |Provides a lat/lon pair of doubles in radians
+getRadianPair :: (Lat p, Lon p) => p -> (LatitudeType, LongitudeType)
+getRadianPair p = (toRadians (lat p), toRadians (lon p))
+
+toRadians :: Floating f => f -> f
+toRadians = (*) (pi / 180)
+
+-- | @interpolate c1 c2 w@ where 0 <= w <= 1 Gives a point on the line
+-- between c1 and c2 equal to @c1 when @w == 0@ (weighted linearly
+-- toward c2).
+interpolate :: (Lat a, Lon a) => a -> a -> Double -> a
+interpolate c1 c2 w
+  | w < 0 || w > 1 = error "Interpolate only works with a weight between zero and one"
+  | otherwise = 
+  let (h,d) = (heading c1 c2, distance c1 c2)
+      v = (d * w, h)
+  in addVector v c1
+     
+-- |@divideArea vDist hDist nw se@ divides an area into a grid of equally
+-- spaced coordinates within the box drawn by the northwest point (nw) and
+-- southeast point (se).  Because this uses floating point there might be a
+-- different number of points in some rows (the last might be too far east based
+-- on a heading from the se point).
+divideArea :: (Lat c, Lon c) => Distance -> Distance -> c -> c -> [[c]]
+divideArea vDist hDist nw se =
+	let (top,left)  = (lat nw, lon nw)
+	    (btm,right) = (lat se, lon se)
+	    columnOne = takeWhile ( (<= west) . heading se) . iterate (addVector (vDist, south)) $ nw
+	    buildRow  = takeWhile ((>= north) . heading se) . iterate (addVector (hDist, east))
+	in map buildRow columnOne
+
+-- |Reads a GPX file (using the GPX library) by simply concatenating all the
+-- tracks, segments, and points ('trkpts', 'trksegs', 'trks') into a single 'Trail'.
+readGPX :: FilePath -> IO (Trail WptType)
+readGPX = liftM (concatMap trkpts . concatMap trksegs . concatMap trks) . readGpxFile
+
+writeGPX :: FilePath -> Trail WptType -> IO ()
+writeGPX fp ps = writeGpxFile fp $ gpx $ gpxType "1.0" "Haskell GPS Package (via the GPX package)" Nothing [] [] [trkType Nothing Nothing Nothing Nothing [] Nothing Nothing Nothing [trksegType ps Nothing]] Nothing
+
+-- writeGpxFile should go in the GPX package
+writeGpxFile :: FilePath -> Gpx -> IO ()
+writeGpxFile fp gpx = runX_ (constA gpx >>> xpickleDocument (xpickle :: PU Gpx) [] fp)
+
+runX_ t = runX t >> return ()
+
+readGPXSegments :: FilePath -> IO [Trail WptType]
+readGPXSegments = liftM (map (concatMap trkpts) . map trksegs . concatMap trks) . readGpxFile
diff --git a/Data/GPS/Trail.hs b/Data/GPS/Trail.hs
new file mode 100644
--- /dev/null
+++ b/Data/GPS/Trail.hs
@@ -0,0 +1,441 @@
+{-# LANGUAGE TupleSections #-}
+module Data.GPS.Trail
+       ( -- * Types         
+         AvgMethod(..)
+       , Selected(..)
+         -- * Utility Functions
+       , isSelected
+       , isNotSelected
+       , onSelected
+       , selLength
+         -- * Trail Functions
+         -- ** Queries
+       , totalDistance
+       , 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
+--       , smoothSegments
+--       , smoothPath
+         ) where
+
+import Text.Show.Functions ()
+import Data.GPS.Core hiding (fix)
+
+import Control.Arrow (first, second)
+import Control.Monad
+import Data.Fixed (mod')
+import Data.Function (on,fix)
+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',xs,zs) = go as (Just a) in (e',a:xs,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 :: (Lat a, Lon a, Time a) => NominalDiffTime -> Trail a -> [(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 :: (Lat a, Lon a, Time a) => 
+                       AvgMethod a -> NominalDiffTime -> Trail a -> [(UTCTime, Speed)]
+slidingAverageSpeed _ _ [] = []
+slidingAverageSpeed m n (x:xs) =
+    let avg = getAvg (x:xs') m
+        avgTime = getAvgTime x (fromMaybe x e)
+    in case avgTime of
+	Nothing -> []
+        Just t  -> (t,avg) : slidingAverageSpeed m n xs
+  where
+  (e,xs',rest) = takeWhileEnd (\c -> getTimeDiff c x <= Just n) xs
+  getTimeDiff a b = on (liftM2 diffUTCTime) getUTCTime a b
+  
+  --  getAvg :: [] -> AvgMethod -> Speed
+  getAvg cs AvgMean =
+    let ss = getSpeedsV cs
+    in case m 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 . join . fmap (speed x) $ e
+        AvgMinOf as -> minimum $ map (getAvg cs) as
+        AvgWith f -> f cs
+  getAvgTime a b = liftM2 addUTCTime (getTimeDiff b a) (getUTCTime a)
+  getSpeedsV = V.fromList . getSpeeds
+  getSpeeds zs = concatMap (maybeToList . uncurry speed) $ zip zs (drop 1 zs)
+
+type TrailTransformation c = [Selected (Trail c)] -> Trail c
+
+-- | 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)]
+
+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 i [] = []
+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 :: (Lat a, Lon a, Time a) => Double -> Double -> PointGrouping a
+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 xs
+  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 :: (Lat a, Lon a, Time a) => Distance -> NominalDiffTime -> PointGrouping a
+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 (getUTCTime a, getUTCTime 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 :: (Lat a, Lon a, Time a) => NominalDiffTime -> PointGrouping a
+spansTime n ps =
+  let times  = mkTimePair ps
+      chunk [] = []
+      chunk (x:xs) =
+        let (good,rest) = span ((<= addUTCTime n (snd x)) . snd) xs in good : chunk rest
+  in map (Select . map fst) $ chunk times
+
+-- | intersects the given groupings
+intersectionOf :: (Lat a, Lon a, Time a) => [PointGrouping a] -> PointGrouping a
+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 :: (Lat a, Lon a, Time a) => [PointGrouping a] -> PointGrouping a
+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
+(/\) xs [] = []
+(/\) [] ys = []
+(/\) xsL@(Select x:xs) ysL@(Select y:ys) =
+  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:xs) ys = Select x : selListDrop (length x) ys
+(\/) xs (Select y:ys) = 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)
+      
+selListTake :: Int -> [Selected [a]] -> [Selected [a]]
+selListTake 0 _ = []
+selListTake n [] = []
+selListTake n (x:xs) =
+  let x' = take n (unSelect x)
+  in fmap (const x') x : selListTake (n - length x') xs
+
+selListDrop :: Int -> [Selected [a]] -> [Selected [a]]
+selListDrop 0 xs = xs
+selListDrop n [] = []
+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 :: (Time a) => Trail a -> [(a,UTCTime)]
+mkTimePair xs =
+  let timesM = map (\x-> fmap (x,) $ getUTCTime 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 :: (Lat a, Lon a, Time a) => [UTCTime] -> Trail a -> Trail a
+bezierCurveAt selectedTimes xs = 
+  let timesDef = mkTimePair xs
+      end = last timesDef
+      top = head timesDef
+      totalTime  = diffUTCTime (snd end) (snd top)
+      times = if null selectedTimes then map snd timesDef else selectedTimes
+      queryTimes = [fromTo (diffUTCTime (snd end) t / totalTime) | t <- times]
+      fromTo = fromRational . toRational
+  in if null timesDef || totalTime == 0 || any (\x -> x < 0 || x > 1) queryTimes
+      then xs
+      else map (bezierPoint xs) queryTimes
+
+bezierPoint :: (Lat a, Lon a) => [a] -> Double -> a
+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 ::  (Lat a, Lon a, Time a) => [Selected (Trail a)] -> Trail a
+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 :: (Lon a, Lat a, Time a) => [a] -> [a]
+linearTime [] = []
+linearTime (p:ps) = go (getUTCTime p) ps
+  where
+  go _ [] = []
+  go t (p:ps) = if getUTCTime p < t then go t ps else p : go (getUTCTime p) ps
+
+-- |Returns the closest distance between two trails (or Nothing if a
+-- trail is empty).  Inefficient implementation:
+-- O( (n * m) * log (n * m) )
+closestDistance :: (Lat a, Lon a) => Trail a -> Trail a -> Maybe Distance
+closestDistance as bs = listToMaybe $ L.sort [distance a b | a <- as, b <- bs]
+
+-- | Find the total distance traveled
+totalDistance :: (Lat a, Lon a) => [a] -> Distance
+totalDistance as = sum $ zipWith distance as (drop 1 as)
+
+-- | 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 :: (Eq c, Lat c, Lon c) => [c] -> [c]
+convexHull xs =
+	let first = southMost xs
+	in case first of
+		Nothing -> []
+		Just f  ->
+	    	     let sorted = sortBy (comparing (eastZeroHeading f)) (filter (/= f) xs)
+		     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 :: (Lat c, Lon c) => c -> c -> Heading
+eastZeroHeading s = (`mod'` (2*pi)) . (+ pi/2) . heading s
+
+data Turn = LeftTurn | RightTurn | Straight deriving (Eq, Ord, Show, Read, Enum)
+
+turn :: (Lat c, Lon c) => c -> c -> c -> 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 :: (Lat c) => [c] -> Maybe c
+southMost []  = Nothing
+southMost cs = Just . minimumBy (comparing lat) $ 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.
+------------------------------------------
+
+smoothRests :: (Lat a, Lon a, Time a) => Trail a -> Trail a
+smoothRests = bezierCurve . refineGrouping (everyNPoints 8) . restLocations 30 60
+
+smoothSegments :: (Lat a, Lon a, Time a) => Trail a -> Trail a
+smoothSegments ps = 
+  let ps' = bezierCurve . everyNPoints 5 $ ps
+      (h,t) = splitAt 2 ps'
+      ps'' = bezierCurve . everyNPoints 5 $ t
+  in h ++ ps''
+
+smoothPath :: (Lat a, Lon a, Time a) => Trail a -> Trail a
+smoothPath ps = undefined
+  
+slidingWindow :: Int -> Int -> ([a] -> [a]) -> [a] -> [a]
+slidingWindow width step f trail = go trail
+  where
+    go [] = []
+    go xs =
+      let (window,e) = splitAt width xs
+          (hT,eT) = splitAt step (f window)
+      in hT ++ go (eT ++ e)
diff --git a/Setup.hs b/Setup.hs
--- a/Setup.hs
+++ b/Setup.hs
diff --git a/gps.cabal b/gps.cabal
--- a/gps.cabal
+++ b/gps.cabal
@@ -1,10 +1,10 @@
 name:		gps
-version:	0.7
+version:	0.8
 license:	BSD3
 license-file:	LICENSE
 author:		Thomas DuBuisson <thomas.dubuisson@gmail.com>
 maintainer:	Thomas DuBuisson
-description:	Useful for manipulating GPS coordinages (in various forms), building paths, and performing basic computations
+description:	Useful for manipulating GPS coordinages (in various forms), building paths, and performing basic computations.  NOTE: Version range 0.8.* won't strictly follow PVP - I will be adding additional functions in minor releases 0.8.x.
 synopsis:	For manipulating GPS coordinates and trails.
 category:	Data
 stability:	stable
@@ -16,7 +16,13 @@
 Library
   Build-Depends: base >= 3 && < 6,
                    pretty >= 1.0 , prettyclass >= 1.0,
-                   time >= 1.1 && < 1.3, GPX == 0.4.*, hxt >= 8.5, xsd == 0.3.*
+                   time >= 1.1, GPX == 0.4.*, hxt >= 8.5, xsd == 0.3.*,
+                   vector >= 0.7, statistics >= 0.9
   hs-source-dirs:
   exposed-modules: Data.GPS
+  other-modules: Data.GPS.Core Data.GPS.Trail
   ghc-options:
+
+source-repository head
+  type:     git
+  location: git@github.com:TomMD/gps.git