{-# LANGUAGE ApplicativeDo #-}
{-# LANGUAGE BangPatterns #-}
{-# LANGUAGE DataKinds #-}
{-# LANGUAGE DerivingStrategies #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE GeneralizedNewtypeDeriving #-}
{-# LANGUAGE Rank2Types #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE TypeFamilies #-}
{-# LANGUAGE TypeInType #-}
{-# LANGUAGE UnboxedTuples #-}
{-# LANGUAGE UndecidableInstances #-}
module Geography.MapAlgebra
(
Raster(..)
, lazy, strict
, RGBARaster(..)
, constant, fromFunction, fromVector, fromRGBA, fromGray
, grayscale
, Histogram(..), histogram, breaks
, invisible
, greenRed, spectrum, blueGreen, purpleYellow, brownBlue
, grayBrown, greenPurple, brownYellow, purpleGreen, purpleRed
, writeImage, writeImageAuto
, displayImage
, png
, Projection(..)
, reproject
, Sphere, LatLng, WebMercator
, Point(..)
, zipWith
, classify
, lmax, lmin
, lmean, lvariety, lmajority, lminority, lvariance
, fsum, fproduct, fmonoid, fmean
, fmax, fmin
, fmajority, fminority, fvariety
, fpercentage, fpercentile
, Line(..)
, flinkage, flength
, Corners(..), Surround(..)
, fpartition, fshape, ffrontage, farea
, Drain(..), Direction(..)
, direction, directions, drainage
, fvolume, fgradient, faspect, faspect', fdownstream, fupstream
, leftPseudo, tau
, D(..), DW(..), S(..), P(..), U(..), B(..), N(..)
, Ix2(..), Comp(..)
, Pixel(..), RGBA, Y
) where
import Control.Concurrent (getNumCapabilities)
import Control.DeepSeq (NFData(..), deepseq)
import Control.Monad.ST
import Data.Bits (testBit)
import Data.Bool (bool)
import qualified Data.ByteString.Lazy as BL
import Data.Foldable
import qualified Data.List as L
import Data.List.NonEmpty (NonEmpty(..))
import qualified Data.List.NonEmpty as NE
import qualified Data.Map.Strict as M
import Data.Massiv.Array hiding (zipWith)
import qualified Data.Massiv.Array as A
import Data.Massiv.Array.IO
import qualified Data.Massiv.Array.Manifest.Vector as A
import Data.Massiv.Array.Unsafe as A
import Data.Proxy (Proxy(..))
import Data.Semigroup
import qualified Data.Set as S
import Data.Typeable (Typeable)
import qualified Data.Vector.Generic as GV
import qualified Data.Vector.Storable as VS
import qualified Data.Vector.Storable.Mutable as VSM
import Data.Word
import GHC.TypeLits
import Graphics.ColorSpace (ColorSpace, Elevator, Pixel(..), RGBA, Y)
import qualified Numeric.LinearAlgebra as LA
import Prelude hiding (zipWith)
import qualified Prelude as P
import Text.Printf (printf)
data Point p = Point { x :: !Double, y :: !Double } deriving (Eq, Show)
class Projection p where
toSphere :: Point p -> Point Sphere
fromSphere :: Point Sphere -> Point p
reproject :: (Projection p, Projection r) => Point p -> Point r
reproject = fromSphere . toSphere
{-# INLINE reproject #-}
data Sphere
instance Projection Sphere where
toSphere = id
fromSphere = id
data LatLng
data WebMercator
newtype Raster u p (r :: Nat) (c :: Nat) a = Raster { _array :: Array u Ix2 a }
instance (Show a, Load (R u) Ix2 a, Extract u Ix2 a) => Show (Raster u p r c a) where
show (Raster a) = show . computeAs B $ extract' zeroIndex (Sz2 (P.min 10 r) (P.min 10 c)) a
where Sz (r :. c) = size a
instance (Eq a, Unbox a) => Eq (Raster U p r c a) where
Raster a == Raster b = a == b
{-# INLINE (==) #-}
instance (Eq a, Storable a) => Eq (Raster S p r c a) where
Raster a == Raster b = a == b
{-# INLINE (==) #-}
instance (Eq a, Prim a) => Eq (Raster P p r c a) where
Raster a == Raster b = a == b
{-# INLINE (==) #-}
instance (Eq a, NFData a) => Eq (Raster N p r c a) where
Raster a == Raster b = a == b
{-# INLINE (==) #-}
instance Eq a => Eq (Raster B p r c a) where
Raster a == Raster b = a == b
{-# INLINE (==) #-}
instance Eq a => Eq (Raster D p r c a) where
Raster a == Raster b = a == b
{-# INLINE (==) #-}
instance Functor (Raster DW p r c) where
fmap f (Raster a) = Raster $ fmap f a
{-# INLINE fmap #-}
instance Functor (Raster D p r c) where
fmap f (Raster a) = Raster $ fmap f a
{-# INLINE fmap #-}
instance Functor (Raster DI p r c) where
fmap f (Raster a) = Raster $ fmap f a
{-# INLINE fmap #-}
instance (KnownNat r, KnownNat c) => Applicative (Raster D p r c) where
pure = constant D Par
{-# INLINE pure #-}
fs <*> as = zipWith ($) fs as
{-# INLINE (<*>) #-}
instance Semigroup a => Semigroup (Raster D p r c a) where
a <> b = zipWith (<>) a b
{-# INLINE (<>) #-}
instance (Monoid a, KnownNat r, KnownNat c) => Monoid (Raster D p r c a) where
mempty = constant D Par mempty
{-# INLINE mempty #-}
a `mappend` b = zipWith mappend a b
{-# INLINE mappend #-}
instance (Num a, KnownNat r, KnownNat c) => Num (Raster D p r c a) where
a + b = zipWith (+) a b
{-# INLINE (+) #-}
a - b = zipWith (-) a b
{-# INLINE (-) #-}
a * b = zipWith (*) a b
{-# INLINE (*) #-}
abs = fmap abs
{-# INLINE abs #-}
signum = fmap signum
{-# INLINE signum #-}
fromInteger = constant D Par . fromInteger
{-# INLINE fromInteger #-}
instance (Fractional a, KnownNat r, KnownNat c) => Fractional (Raster D p r c a) where
a / b = zipWith (/) a b
{-# INLINE (/) #-}
fromRational = constant D Par . fromRational
{-# INLINE fromRational #-}
instance Foldable (Raster D p r c) where
foldMap f (Raster a) = foldMap f a
{-# INLINE foldMap #-}
sum (Raster a) = A.sum a
{-# INLINE sum #-}
product (Raster a) = A.product a
{-# INLINE product #-}
length (Raster a) = (\(Sz2 r c) -> r * c) $ A.size a
{-# INLINE length #-}
lazy :: Source u Ix2 a => Raster u p r c a -> Raster D p r c a
lazy (Raster a) = Raster $ delay a
{-# INLINE lazy #-}
strict :: (Load v Ix2 a, Mutable u Ix2 a) => u -> Raster v p r c a -> Raster u p r c a
strict u (Raster a) = Raster $ computeAs u a
{-# INLINE strict #-}
constant :: (KnownNat r, KnownNat c, Construct u Ix2 a) => u -> Comp -> a -> Raster u p r c a
constant u c a = fromFunction u c (const a)
{-# INLINE constant #-}
fromFunction :: forall u p r c a. (KnownNat r, KnownNat c, Construct u Ix2 a) =>
u -> Comp -> (Ix2 -> a) -> Raster u p r c a
fromFunction u c f = Raster $ makeArrayR u c (Sz sh) f
where sh = fromInteger (natVal (Proxy :: Proxy r)) :. fromInteger (natVal (Proxy :: Proxy c))
{-# INLINE fromFunction #-}
fromVector :: forall v p r c a. (KnownNat r, KnownNat c, GV.Vector v a, Mutable (A.ARepr v) Ix2 a, Typeable v) =>
Comp -> v a -> Either String (Raster (A.ARepr v) p r c a)
fromVector comp v | (r * c) == GV.length v = Right . Raster $ A.fromVector' comp (Sz2 r c) v
| otherwise = Left $ printf "Expected Pixel Count: %d - Actual: %d" (r * c) (GV.length v)
where r = fromInteger $ natVal (Proxy :: Proxy r)
c = fromInteger $ natVal (Proxy :: Proxy c)
{-# INLINE fromVector #-}
data RGBARaster p r c a = RGBARaster { _red :: !(Raster S p r c a)
, _green :: !(Raster S p r c a)
, _blue :: !(Raster S p r c a)
, _alpha :: !(Raster S p r c a) }
fromRGBA :: forall p r c a. (Elevator a, KnownNat r, KnownNat c) => FilePath -> IO (Either String (RGBARaster p r c a))
fromRGBA fp = do
cap <- getNumCapabilities
img <- setComp (bool Par Seq $ cap == 1) <$> readImageAuto fp
let rows = fromInteger $ natVal (Proxy :: Proxy r)
cols = fromInteger $ natVal (Proxy :: Proxy c)
Sz (r :. c) = size img
if r == rows && c == cols
then do
(ar, ag, ab, aa) <- spreadRGBA img
pure . Right $ RGBARaster (Raster ar) (Raster ag) (Raster ab) (Raster aa)
else pure . Left $ printf "Expected Size: %d x %d - Actual Size: %d x %d" rows cols r c
{-# INLINE fromRGBA #-}
spreadRGBA :: (Index ix, Elevator e)
=> A.Array S ix (Pixel RGBA e)
-> IO (A.Array S ix e, A.Array S ix e, A.Array S ix e, A.Array S ix e)
spreadRGBA arr = do
let sz = A.size arr
mr <- A.unsafeNew sz
mb <- A.unsafeNew sz
mg <- A.unsafeNew sz
ma <- A.unsafeNew sz
flip A.imapM_ arr $ \i (PixelRGBA r g b a) -> do
A.unsafeWrite mr i r
A.unsafeWrite mg i g
A.unsafeWrite mb i b
A.unsafeWrite ma i a
ar <- A.unsafeFreeze (getComp arr) mr
ag <- A.unsafeFreeze (getComp arr) mg
ab <- A.unsafeFreeze (getComp arr) mb
aa <- A.unsafeFreeze (getComp arr) ma
return (ar, ag, ab, aa)
{-# INLINE spreadRGBA #-}
fromGray :: forall p r c a. (Elevator a, KnownNat r, KnownNat c) => FilePath -> IO (Either String (Raster S p r c a))
fromGray fp = do
cap <- getNumCapabilities
img <- setComp (bool Par Seq $ cap == 1) <$> readImageAuto fp
let rows = fromInteger $ natVal (Proxy :: Proxy r)
cols = fromInteger $ natVal (Proxy :: Proxy c)
Sz (r :. c) = size img
pure . bool (Left $ printf "Expected Size: %d x %d - Actual Size: %d x %d" rows cols r c) (Right $ f img) $ r == rows && c == cols
where f :: Image S Y a -> Raster S p r c a
f img = Raster . A.fromVector' (getComp img) (size img) . VS.unsafeCast $ A.toVector img
{-# INLINE fromGray #-}
invisible :: Pixel RGBA Word8
invisible = PixelRGBA 0 0 0 0
ramp :: Ord k => [(Word8, Word8, Word8)] -> [k] -> M.Map k (Pixel RGBA Word8)
ramp colours bs = M.fromList . P.zip bs $ P.map (\(r,g,b) -> PixelRGBA r g b maxBound) colours
{-# INLINE ramp #-}
greenRed :: Ord k => [k] -> M.Map k (Pixel RGBA Word8)
greenRed = ramp colours
where colours = [ (0, 48, 0), (31, 79, 20), (100, 135, 68), (148, 193, 28), (193, 242, 3)
, (241, 255, 159), (249, 228, 227), (202, 145, 150), (153, 101, 97), (142, 38 ,18) ]
spectrum :: Ord k => [k] -> M.Map k (Pixel RGBA Word8)
spectrum = ramp colours
where colours = [ (0, 22, 51), (51, 18, 135), (150, 0, 204), (242, 13, 177), (255, 61, 61)
, (240, 152, 56), (248, 230, 99), (166, 249, 159), (184, 249, 212), (216, 230, 253) ]
blueGreen :: Ord k => [k] -> M.Map k (Pixel RGBA Word8)
blueGreen = ramp colours
where colours = [ (29, 43, 53), (37, 44, 95), (63, 70, 134), (89, 112, 147), (87, 124, 143)
, (117, 160, 125), (188, 219, 173), (239, 253, 163), (222, 214, 67), (189, 138, 55) ]
purpleYellow :: Ord k => [k] -> M.Map k (Pixel RGBA Word8)
purpleYellow = ramp colours
where colours = [ (90, 89, 78), (73, 65, 132), (107, 86, 225), (225, 67, 94), (247, 55, 55)
, (251, 105, 46), (248, 174, 66), (249, 219, 25), (255, 255, 0), (242, 242, 242) ]
brownBlue :: Ord k => [k] -> M.Map k (Pixel RGBA Word8)
brownBlue = ramp colours
where colours = [ (27, 36, 43), (86, 52, 42), (152, 107, 65), (182, 176, 152), (215, 206, 191)
, (198, 247, 0), (53, 227, 0), (30, 158, 184), (22, 109, 138), (12, 47, 122) ]
grayBrown :: Ord k => [k] -> M.Map k (Pixel RGBA Word8)
grayBrown = ramp colours
where colours = [ (64, 57, 88), (95, 96, 116), (158, 158, 166), (206, 208, 197), (215, 206, 191)
, (186, 164, 150), (160, 124, 98), (117, 85, 72), (90, 70, 63), (39, 21, 17) ]
greenPurple :: Ord k => [k] -> M.Map k (Pixel RGBA Word8)
greenPurple = ramp colours
where colours = [ (89, 168, 15), (158, 213, 76), (196, 237, 104), (226, 255, 158), (240, 242, 221)
, (248, 202, 140), (233, 161, 137), (212, 115, 132), (172, 67, 123), (140, 40, 110) ]
brownYellow :: Ord k => [k] -> M.Map k (Pixel RGBA Word8)
brownYellow = ramp colours
where colours = [ (96, 72, 96), (120, 72, 96), (168, 96, 96), (192, 120, 96), (240, 168, 72)
, (248, 202, 140), (254, 236, 174), (255, 244, 194), (255, 247, 219), (255, 252, 246) ]
purpleGreen :: Ord k => [k] -> M.Map k (Pixel RGBA Word8)
purpleGreen = ramp colours
where colours = [ (80, 73, 113), (117, 64, 152), (148, 116, 180), (199, 178, 214), (223, 204, 228)
, (218, 234, 193), (171, 214, 155), (109, 192, 103), (13, 177, 75), (57, 99, 83) ]
purpleRed :: Ord k => [k] -> M.Map k (Pixel RGBA Word8)
purpleRed = ramp colours
where colours = [ (51, 60, 255), (76, 60, 233), (99, 60, 211), (121, 60, 188), (155, 60, 155)
, (166, 60, 143), (188, 60, 121), (206, 60, 94), (217, 60, 83), (255, 60, 76) ]
grayscale :: Functor (Raster u p r c) => Raster u p r c a -> Raster u p r c (Pixel Y a)
grayscale = fmap PixelY
{-# INLINE grayscale #-}
png :: (Source u Ix2 (Pixel cs a), ColorSpace cs a) => Raster u p r c (Pixel cs a) -> BL.ByteString
png (Raster a) = encode PNG def a
{-# INLINE png #-}
classify :: (Ord a, Functor f) => b -> M.Map a b -> f a -> f b
classify d m r = fmap f r
where f a = maybe d snd $ M.lookupLE a m
{-# INLINE classify #-}
lmin :: (Ord a, Source u Ix2 a) => Raster u p r c a -> Raster u p r c a -> Raster D p r c a
lmin = zipWith P.min
{-# INLINE lmin #-}
lmax :: (Ord a, Source u Ix2 a) => Raster u p r c a -> Raster u p r c a -> Raster D p r c a
lmax = zipWith P.max
{-# INLINE lmax #-}
lmean :: (Real a, Fractional b, KnownNat r, KnownNat c) => NonEmpty (Raster D p r c a) -> Raster D p r c b
lmean (a :| [b]) = Raster $ A.zipWith (\n m -> realToFrac (n + m) / 2) (_array a) (_array b)
lmean (a :| [b,c]) = Raster $ A.zipWith3 (\n m o -> realToFrac (n + m + o) / 3) (_array a) (_array b) (_array c)
lmean (a :| as) = (\n -> realToFrac n / len) <$> foldl' (+) a as
where len = 1 + fromIntegral (length as)
{-# INLINE lmean #-}
lvariety :: (KnownNat r, KnownNat c, Eq a) => NonEmpty (Raster D p r c a) -> Raster D p r c Word
lvariety = fmap (fromIntegral . length . NE.nub) . P.sequenceA
{-# INLINE lvariety #-}
lmajority :: (KnownNat r, KnownNat c, Ord a) => NonEmpty (Raster D p r c a) -> Raster D p r c a
lmajority = fmap majo . P.sequenceA
{-# INLINE lmajority #-}
majo :: forall t a. (Foldable t, Ord a) => t a -> a
majo = fst . g . f
where
f :: t a -> M.Map a Int
f = foldl' (\m a -> M.insertWith (+) a 1 m) M.empty
g :: M.Map a Int -> (a, Int)
g = L.foldl1' (\(a,c) (k,v) -> if c < v then (k,v) else (a,c)) . M.toList
{-# INLINE majo #-}
lminority :: (KnownNat r, KnownNat c, Ord a) => NonEmpty (Raster D p r c a) -> Raster D p r c a
lminority = fmap mino . P.sequenceA
{-# INLINE lminority #-}
mino :: forall t a. (Foldable t, Ord a) => t a -> a
mino = fst . g . f
where
f :: t a -> M.Map a Int
f = foldl' (\m a -> M.insertWith (+) a 1 m) M.empty
g :: M.Map a Int -> (a, Int)
g = L.foldl1' (\(a,c) (k,v) -> if c > v then (k,v) else (a,c)) . M.toList
{-# INLINE mino #-}
lvariance
:: forall p a r c. (Real a, KnownNat r, KnownNat c)
=> NonEmpty (Raster D p r c a) -> Maybe (Raster D p r c Double)
lvariance (_ :| []) = Nothing
lvariance rs = Just (f <$> P.sequenceA rs)
where
len :: Double
len = realToFrac $ length rs
avg :: NonEmpty a -> Double
avg ns = (\z -> realToFrac z / len) $ foldl' (+) 0 ns
f :: NonEmpty a -> Double
f os@(n :| ns) = num / (len - 1)
where
num = foldl' (\acc m -> acc + ((realToFrac m - av) ^ (2 :: Int))) ((realToFrac n - av) ^ (2 :: Int)) ns
av = avg os
{-# INLINE lvariance #-}
zipWith :: (Source u Ix2 a, Source u Ix2 b) =>
(a -> b -> d) -> Raster u p r c a -> Raster u p r c b -> Raster D p r c d
zipWith f (Raster a) (Raster b) = Raster $ A.zipWith f a b
{-# INLINE zipWith #-}
fsum :: (Num a, Manifest u Ix2 a) => Raster u p r c a -> Raster DW p r c a
fsum (Raster a) = Raster $ mapStencil (Fill 0) (sumStencil $ Sz2 3 3) a
{-# INLINE fsum #-}
fproduct :: (Num a, Manifest u Ix2 a) => Raster u p r c a -> Raster DW p r c a
fproduct (Raster a) = Raster $ mapStencil (Fill 1) (productStencil $ Sz2 3 3) a
{-# INLINE fproduct #-}
fmonoid :: (Monoid a, Manifest u Ix2 a) => Raster u p r c a -> Raster DW p r c a
fmonoid (Raster a) = Raster $ mapStencil (Fill mempty) (foldStencil $ Sz2 3 3) a
{-# INLINE fmonoid #-}
fmean :: (Fractional a, Manifest u Ix2 a) => Raster u p r c a -> Raster DW p r c a
fmean (Raster a) = Raster $ mapStencil (Fill 0) (avgStencil $ Sz2 3 3) a
{-# INLINE fmean #-}
fmax :: (Ord a, Bounded a, Manifest u Ix2 a) => Raster u p r c a -> Raster DW p r c a
fmax (Raster a) = Raster $ mapStencil Edge (maxStencil $ Sz2 3 3) a
{-# INLINE fmax #-}
fmin :: (Ord a, Bounded a, Manifest u Ix2 a) => Raster u p r c a -> Raster DW p r c a
fmin (Raster a) = Raster $ mapStencil Edge (minStencil $ Sz2 3 3) a
{-# INLINE fmin #-}
fvariety :: (Ord a, Default a, Manifest u Ix2 a) => Raster u p r c a -> Raster DW p r c Word
fvariety (Raster a) = Raster $ mapStencil Edge (neighbourhoodStencil f) a
where f nw no ne we fo ea sw so se = fromIntegral . length $ L.nub [ nw, no, ne, we, fo, ea, sw, so, se ]
{-# INLINE fvariety #-}
fmajority :: (Ord a, Default a, Manifest u Ix2 a) => Raster u p r c a -> Raster DW p r c a
fmajority (Raster a) = Raster $ mapStencil Continue (neighbourhoodStencil f) a
where f nw no ne we fo ea sw so se = majo [ nw, no, ne, we, fo, ea, sw, so, se ]
{-# INLINE fmajority #-}
fminority :: (Ord a, Default a, Manifest u Ix2 a) => Raster u p r c a -> Raster DW p r c a
fminority (Raster a) = Raster $ mapStencil Continue (neighbourhoodStencil f) a
where f nw no ne we fo ea sw so se = mino [ nw, no, ne, we, fo, ea, sw, so, se ]
{-# INLINE fminority #-}
fpercentage :: (Eq a, Default a, Manifest u Ix2 a) => Raster u p r c a -> Raster DW p r c Double
fpercentage (Raster a) = Raster $ mapStencil Continue (neighbourhoodStencil f) a
where f nw no ne we fo ea sw so se = ( bool 0 1 (nw == fo)
+ bool 0 1 (no == fo)
+ bool 0 1 (ne == fo)
+ bool 0 1 (we == fo)
+ bool 0 1 (ea == fo)
+ bool 0 1 (sw == fo)
+ bool 0 1 (so == fo)
+ bool 0 1 (se == fo) ) / 8
{-# INLINE fpercentage #-}
fpercentile :: (Ord a, Default a, Manifest u Ix2 a) => Raster u p r c a -> Raster DW p r c Double
fpercentile (Raster a) = Raster $ mapStencil Continue (neighbourhoodStencil f) a
where f nw no ne we fo ea sw so se = ( bool 0 1 (nw < fo)
+ bool 0 1 (no < fo)
+ bool 0 1 (ne < fo)
+ bool 0 1 (we < fo)
+ bool 0 1 (ea < fo)
+ bool 0 1 (sw < fo)
+ bool 0 1 (so < fo)
+ bool 0 1 (se < fo) ) / 8
{-# INLINE fpercentile #-}
flinkage :: (Default a, Eq a, Manifest u Ix2 a) => Raster u p r c a -> Raster DW p r c Line
flinkage (Raster a) = Raster $ mapStencil (Fill def) (neighbourhoodStencil linkage) a
{-# INLINE flinkage #-}
newtype Line = Line { _line :: Word8 }
deriving stock (Eq, Ord, Show)
deriving newtype (Storable, Prim, Default)
instance NFData Line where
rnf (Line a) = deepseq a ()
linkage :: Eq a => a -> a -> a -> a -> a -> a -> a -> a -> a -> Line
linkage nw no ne we fo ea sw so se = Line $ axes + diags
where axes = bool 0 2 (no == fo) + bool 0 8 (we == fo) + bool 0 16 (ea == fo) + bool 0 64 (so == fo)
diags = bool 0 1 (nw == fo && not (testBit axes 1 || testBit axes 3))
+ bool 0 4 (ne == fo && not (testBit axes 1 || testBit axes 4))
+ bool 0 32 (sw == fo && not (testBit axes 3 || testBit axes 6))
+ bool 0 128 (se == fo && not (testBit axes 4 || testBit axes 6))
{-# INLINE linkage #-}
data Direction = East | NorthEast | North | NorthWest | West | SouthWest | South | SouthEast
deriving (Eq, Ord, Enum, Show)
flength :: Functor (Raster u p r c) => Raster u p r c Line -> Raster u p r c Double
flength = fmap f
where half = 1 / 2
root = 1 / sqrt 2
f (Line a) = bool 0 half (testBit a 1)
+ bool 0 half (testBit a 3)
+ bool 0 half (testBit a 4)
+ bool 0 half (testBit a 6)
+ bool 0 root (testBit a 0)
+ bool 0 root (testBit a 2)
+ bool 0 root (testBit a 5)
+ bool 0 root (testBit a 7)
{-# INLINE flength #-}
data Corners = Corners { _topLeft :: !Surround
, _bottomLeft :: !Surround
, _bottomRight :: !Surround
, _topRight :: !Surround } deriving (Eq, Show)
instance NFData Corners where
rnf (Corners tl bl br tr) = tl `deepseq` bl `deepseq` br `deepseq` tr `deepseq` ()
data Surround = Complete
| OneSide
| Open
| RightAngle
| OutFlow
deriving (Eq, Ord, Show)
instance NFData Surround where
rnf s = case s of
Complete -> ()
OneSide -> ()
Open -> ()
RightAngle -> ()
OutFlow -> ()
surround :: Eq a => a -> a -> a -> a -> Surround
surround fo tl tr br
| up && tl == tr && tr == br = Complete
| up && right = RightAngle
| (up && diag) || (diag && right) = OneSide
| diag && fo == tl && fo == br = OutFlow
| otherwise = Open
where up = fo /= tl
diag = fo /= tr
right = fo /= br
{-# INLINE surround #-}
frontage :: Corners -> Double
frontage (Corners tl bl br tr) = f tl + f bl + f br + f tr
where f Complete = 1 / sqrt 2
f OneSide = 1 / 2
f Open = 0
f RightAngle = 1
f OutFlow = 1 / sqrt 2
fpartition :: (Default a, Eq a, Manifest u Ix2 a) => Raster u p r c a -> Raster DW p r c Corners
fpartition (Raster a) = Raster $ mapStencil Reflect partStencil a
{-# INLINE fpartition #-}
partStencil :: (Eq a, Default a) => Stencil Ix2 a Corners
partStencil = makeStencil (Sz2 2 2) (1 :. 0) $ \f -> do
tl <- f (-1 :. 0)
tr <- f (-1 :. 1)
br <- f (0 :. 1)
fo <- f (0 :. 0)
pure $ Corners (surround fo tl tl fo) Open (surround fo fo br br) (surround fo tl tr br)
{-# INLINE partStencil #-}
fshape :: (Default a, Eq a, Manifest u Ix2 a) => Raster u p r c a -> Raster DW p r c Corners
fshape (Raster a) = Raster $ mapStencil Reflect (neighbourhoodStencil f) a
where f nw no ne we fo ea sw so se = Corners (surround fo no nw we)
(surround fo so sw we)
(surround fo so se ea)
(surround fo no ne ea)
{-# INLINE fshape #-}
ffrontage :: Functor (Raster u p r c) => Raster u p r c Corners -> Raster u p r c Double
ffrontage = fmap frontage
{-# INLINE ffrontage #-}
area :: Corners -> Double
area (Corners tl bl br tr) = 1 - f tl - f bl - f br - f tr
where f Complete = 1/8
f OutFlow = -1/8
f _ = 0
{-# INLINE area #-}
farea :: Functor (Raster u p r c) => Raster u p r c Corners -> Raster u p r c Double
farea = fmap area
{-# INLINE farea #-}
fvolume :: (Fractional a, Default a, Manifest u Ix2 a) => Raster u p r c a -> Raster DW p r c a
fvolume (Raster a) = Raster $ mapStencil Reflect (neighbourhoodStencil volume) a
{-# INLINE fvolume #-}
volume :: Fractional a => a -> a -> a -> a -> a -> a -> a -> a -> a -> a
volume tl up tr le fo ri bl bo br =
((fo * 8)
+ nw + no
+ no + ne
+ ne + ea
+ ea + se
+ se + so
+ so + sw
+ sw + we
+ we + nw) / 24
where nw = (tl + up + le + fo) / 4
no = (up + fo) / 2
ne = (up + tr + fo + ri) / 4
we = (le + fo) / 2
ea = (fo + ri) / 2
sw = (le + fo + bl + bo) / 4
so = (fo + bo) / 2
se = (fo + ri + bo + br) / 4
{-# INLINE volume #-}
neighbourhood :: Applicative f => (a -> a -> a -> a -> a -> a -> a -> a -> a -> b) -> (Ix2 -> f a) -> f b
neighbourhood g f = g <$> f (-1 :. -1) <*> f (-1 :. 0) <*> f (-1 :. 1)
<*> f (0 :. -1) <*> f (0 :. 0) <*> f (0 :. 1)
<*> f (1 :. -1) <*> f (1 :. 0) <*> f (1 :. 1)
{-# INLINE neighbourhood #-}
neighbourhoodStencil :: Default a => (a -> a -> a -> a -> a -> a -> a -> a -> a -> b) -> Stencil Ix2 a b
neighbourhoodStencil f = makeStencil (Sz2 3 3) (1 :. 1) (neighbourhood f)
{-# INLINE neighbourhoodStencil #-}
facetStencil :: (Fractional a, Default a) => (a -> a -> a -> a -> a -> a -> a -> a -> a -> b) -> Stencil Ix2 a b
facetStencil f = makeStencil (Sz2 3 3) (1 :. 1) (neighbourhood g)
where g nw no ne we fo ea sw so se = f ((nw + no + we + fo) / 4)
((no + fo) / 2)
((no + ne + fo + ea) / 4)
((we + fo) / 2)
fo
((fo + ea) / 2)
((we + fo + sw + so) / 4)
((fo + so) / 2)
((fo + ea + so + se) / 4)
{-# INLINE facetStencil #-}
leftPseudo :: LA.Matrix Double
leftPseudo = LA.inv (aT <> a) <> aT
where aT = LA.tr' a
a = LA.matrix 3 [ -0.5, -0.5, 1
, -0.5, 0, 1
, -0.5, 0.5, 1
, 0, -0.5, 1
, 0, 0, 1
, 0, 0.5, 1
, 0.5, -0.5, 1
, 0.5, 0, 1
, 0.5, 0.5, 1 ]
fgradient :: (Manifest u Ix2 Double) => Raster u p r c Double -> Raster DW p r c Double
fgradient (Raster a) = Raster $ mapStencil Reflect (facetStencil gradient) a
{-# INLINE fgradient #-}
tau :: Double
tau = 6.283185307179586
gradient :: Double -> Double -> Double -> Double -> Double -> Double -> Double -> Double -> Double -> Double
gradient nw no ne we fo ea sw so se = (tau / 2) - acos (normal vs LA.! 2)
where vs = [ nw, no, ne, we, fo, ea, sw, so, se ]
normal :: [Double] -> LA.Vector Double
normal = LA.normalize . zcoord (-1) . normal'
normal' :: [Double] -> LA.Vector Double
normal' vs = leftPseudo LA.#> LA.vector vs
zcoord :: Double -> LA.Vector Double -> LA.Vector Double
zcoord n v = LA.vector [ v LA.! 0, v LA.! 1, n ]
faspect :: Manifest u Ix2 Double => Raster u p r c Double -> Raster DW p r c (Maybe Double)
faspect (Raster a) = Raster $ mapStencil Reflect (facetStencil f) a
where f nw no ne we fo ea sw so se = case normal' [ nw, no, ne, we, fo, ea, sw, so, se ] of
n | ((n LA.! 0) =~ 0) && ((n LA.! 1) =~ 0) -> Nothing
| otherwise -> Just $ angle (LA.normalize $ zcoord 0 n) axis
axis = LA.vector [1, 0, 0]
{-# INLINE faspect #-}
faspect' :: Manifest u Ix2 Double => Raster u p r c Double -> Raster DW p r c Double
faspect' (Raster a) = Raster $ mapStencil Reflect (facetStencil f) a
where f nw no ne we fo ea sw so se = angle (LA.normalize $ zcoord 0 $ normal' [ nw, no, ne, we, fo, ea, sw, so , se ]) axis
axis = LA.vector [1, 0, 0]
{-# INLINE faspect' #-}
(=~) :: Double -> Double -> Bool
a =~ b = abs (a - b) < 0.0061359
angle :: LA.Vector Double -> LA.Vector Double -> Double
angle u v = acos $ LA.dot u v
newtype Drain = Drain { _drain :: Word8 }
deriving stock (Eq, Ord, Show)
deriving newtype (Storable, Prim, Default)
instance NFData Drain where
rnf (Drain a) = deepseq a ()
fdownstream :: Manifest u Ix2 Double => Raster u p r c Double -> Raster DW p r c Drain
fdownstream (Raster a) = Raster $ mapStencil Reflect (facetStencil downstream) a
{-# INLINE fdownstream #-}
downstream :: Double -> Double -> Double -> Double -> Double -> Double -> Double -> Double -> Double -> Drain
downstream nw no ne we fo ea sw so se = Drain . snd $ foldl' f (0, 0) angles
where f (!curr, !s) (!a, !d) | a =~ curr = (curr, s + d)
| a > curr = (a, d)
| otherwise = (curr, s)
angles = [ (fo - nw, 1)
, (fo - no, 2)
, (fo - ne, 4)
, (fo - we, 8)
, (fo - ea, 16)
, (fo - sw, 32)
, (fo - so, 64)
, (fo - se, 128) ]
fupstream :: Manifest u Ix2 Drain => Raster u p r c Drain -> Raster DW p r c Drain
fupstream (Raster a) = Raster $ mapStencil (Fill $ Drain 0) (neighbourhoodStencil f) a
where f nw no ne we _ ea sw so se = Drain $ bool 0 1 (testBit (_drain nw) 7)
+ bool 0 2 (testBit (_drain no) 6)
+ bool 0 4 (testBit (_drain ne) 5)
+ bool 0 8 (testBit (_drain we) 4)
+ bool 0 16 (testBit (_drain ea) 3)
+ bool 0 32 (testBit (_drain sw) 2)
+ bool 0 64 (testBit (_drain so) 1)
+ bool 0 128 (testBit (_drain se) 0)
{-# INLINE fupstream #-}
direction :: Direction -> Drain -> Bool
direction dir (Drain d) = case dir of
NorthWest -> testBit d 0
North -> testBit d 1
NorthEast -> testBit d 2
West -> testBit d 3
East -> testBit d 4
SouthWest -> testBit d 5
South -> testBit d 6
SouthEast -> testBit d 7
directions :: Drain -> S.Set Direction
directions d = S.fromList $ foldl' (\acc dir -> bool acc (dir : acc) $ direction dir d) [] dirs
where dirs = [NorthWest, North, NorthEast, West, East, SouthWest, South, SouthEast]
drainage :: S.Set Direction -> Drain
drainage = Drain . S.foldl' f 0
where f acc d = case d of
NorthWest -> acc + 1
North -> acc + 2
NorthEast -> acc + 4
West -> acc + 8
East -> acc + 16
SouthWest -> acc + 32
South -> acc + 64
SouthEast -> acc + 128
newtype Histogram = Histogram { _histogram :: VS.Vector Word } deriving (Eq, Show)
histogram :: Source u Ix2 Word8 => Raster u p r c Word8 -> Histogram
histogram (Raster a) = runST $ do
acc <- VSM.replicate 256 0
A.mapM_ (VSM.unsafeModify acc succ . fromIntegral) a
Histogram <$> VS.unsafeFreeze acc
{-# INLINE histogram #-}
breaks :: Histogram -> [Word8]
breaks (Histogram h) = take 10 . (1 :) . P.reverse . snd . VS.ifoldl' f (binWidth, []) . VS.postscanl' (+) 0 $ VS.tail h
where binWidth = VS.sum (VS.tail h) `div` 11
f a@(goal, acc) i n | n > goal = (next n goal, (fromIntegral i + 1) : acc)
| otherwise = a
next n goal | (n - goal) > binWidth = goal + (binWidth * (((n - goal) `div` binWidth) + 1))
| otherwise = goal + binWidth