module Geography.VectorTile.Geometry
(
Geometry(..)
, Point(..)
, LineString(..)
, Polygon(..)
, area
, surveyor
, distance
, Command(..)
, commands
, uncommands
, zig
, unzig
) where
import Control.DeepSeq (NFData)
import Control.Monad.Trans.State.Lazy
import Data.Bits
import Data.Foldable (foldlM)
import Data.Int
import Data.Monoid
import Data.Text (Text)
import qualified Data.Text as T
import qualified Data.Vector as V
import qualified Data.Vector.Unboxed as U
import Data.Word
import GHC.Generics (Generic)
import Geography.VectorTile.Util
import Text.Printf.TH
type Point = (Int,Int)
pattern Point :: Int -> Int -> (Int, Int)
pattern Point{x, y} = (x, y)
instance Monoid Point where
mempty = Point 0 0
(Point a b) `mappend` (Point a' b') = Point (a + a') (b + b')
newtype LineString = LineString { lsPoints :: U.Vector Point } deriving (Eq,Show,Generic)
instance NFData LineString
data Polygon = Polygon { polyPoints :: U.Vector Point
, inner :: V.Vector Polygon } deriving (Eq,Show,Generic)
instance NFData Polygon
area :: Polygon -> Float
area p = surveyor (polyPoints p) + sum (V.map area $ inner p)
surveyor :: U.Vector Point -> Float
surveyor v = (/ 2) . fromIntegral . U.sum $ U.zipWith3 (\xn yn yp -> xn * (yn yp)) xs yns yps
where v' = U.init v
xs = U.map x v'
yns = U.map y . U.tail $ U.snoc v' (U.head v')
yps = U.map y . U.init $ U.cons (U.last v') v'
distance :: Point -> Point -> Float
distance p1 p2 = sqrt . fromIntegral $ dx ^ 2 + dy ^ 2
where dx = x p1 x p2
dy = y p1 y p2
class Geometry g where
fromCommands :: [Command] -> Either Text (V.Vector g)
toCommands :: V.Vector g -> [Command]
instance Geometry Point where
fromCommands (MoveTo ps : []) = Right . U.convert $ evalState (U.mapM expand ps) (0,0)
fromCommands (c:_) = Left $ [st|Invalid command found in Point feature: %s|] (show c)
fromCommands [] = Left "No points given!"
toCommands ps = [MoveTo $ evalState (U.mapM collapse $ U.convert ps) (0,0)]
instance Geometry LineString where
fromCommands cs = evalState (f cs) (0,0)
where f (MoveTo p : LineTo ps : rs) = fmap . V.cons <$> ls <*> f rs
where ls = LineString <$> U.mapM expand (p <> ps)
f [] = pure $ Right V.empty
f _ = pure $ Left "LineString decode: Invalid command sequence given."
toCommands ls = concat $ evalState (mapM f ls) (0,0)
where f (LineString ps) = do
curr <- get
let (h,t) = (U.head ps, U.tail ps)
put h
l <- U.mapM collapse t
pure [MoveTo $ U.singleton (x h x curr, y h y curr), LineTo l]
instance Geometry Polygon where
fromCommands cs = do
ps <- evalState (f cs) (0,0)
let (h,t) = (V.head ps, V.tail ps)
(ps',p') = runState (foldlM g V.empty t) h
pure $ V.snoc ps' p'
where f (MoveTo p : LineTo ps : ClosePath : rs) = do
curr <- get
let h = U.head p
here = (x h + x curr, y h + y curr)
po <- flip U.snoc here <$> U.mapM expand (U.cons h ps)
fmap (V.cons (Polygon po V.empty)) <$> f rs
f [] = pure $ Right V.empty
f _ = pure . Left $ [st|Polygon decode: Invalid command sequence given: %s|] (show cs)
g acc p | area p > 0 = do
curr <- get
put p
pure $ V.snoc acc curr
| otherwise = do
modify (\s -> s { inner = V.snoc (inner s) p })
pure acc
toCommands ps = concat $ evalState (mapM f ps) (0,0)
where f (Polygon p i) = do
curr <- get
let (h,t) = (U.head p, U.tail $ U.init p)
put h
l <- U.mapM collapse t
let cs = [MoveTo $ U.singleton (x h x curr, y h y curr), LineTo l, ClosePath]
concat . V.cons cs <$> mapM f i
data Command = MoveTo (U.Vector (Int,Int))
| LineTo (U.Vector (Int,Int))
| ClosePath deriving (Eq,Show)
zig :: Int -> Word32
zig n = fromIntegral $ shift n 1 `xor` shift n (63)
unzig :: Word32 -> Int
unzig n = fromIntegral (fromIntegral unzigged :: Int32)
where unzigged = shift n (1) `xor` negate (n .&. 1)
parseCmd :: Word32 -> Either T.Text (Int,Int)
parseCmd n = case (cmd,count) of
(1,m) -> Right $ both fromIntegral (1,m)
(2,m) -> Right $ both fromIntegral (2,m)
(7,1) -> Right (7,1)
(7,m) -> Left $ "ClosePath was given a parameter count: " <> T.pack (show m)
(m,_) -> Left $ [st|Invalid command integer %d found in: %X|] m n
where cmd = n .&. 7
count = shift n (3)
unparseCmd :: (Int,Int) -> Word32
unparseCmd (cmd,count) = fromIntegral $ (cmd .&. 7) .|. shift count 3
commands :: [Word32] -> Either T.Text [Command]
commands [] = Right []
commands (n:ns) = parseCmd n >>= f
where f (1,count) = do
mts <- MoveTo . U.fromList . map (both unzig) <$> pairs (take (count * 2) ns)
(mts :) <$> commands (drop (count * 2) ns)
f (2,count) = do
mts <- LineTo . U.fromList . map (both unzig) <$> pairs (take (count * 2) ns)
(mts :) <$> commands (drop (count * 2) ns)
f (7,_) = (ClosePath :) <$> commands ns
f _ = Left "Sentinel: You should never see this."
uncommands :: [Command] -> [Word32]
uncommands = U.toList . U.concat . map f
where f (MoveTo ps) = (U.cons $ unparseCmd (1, U.length ps)) $ params ps
f (LineTo ls) = (U.cons $ unparseCmd (2, U.length ls)) $ params ls
f ClosePath = U.singleton $ unparseCmd (7,1)
params :: U.Vector (Int,Int) -> U.Vector Word32
params = U.foldr (\(a,b) acc -> U.cons (zig a) $ U.cons (zig b) acc) U.empty
expand :: (Int,Int) -> State (Int,Int) Point
expand p = do
curr <- get
let here = (x p + x curr, y p + y curr)
put here
pure here
collapse :: Point -> State (Int,Int) (Int,Int)
collapse p = do
curr <- get
let diff = (x p x curr, y p y curr)
put p
pure diff