{-# LANGUAGE Rank2Types, DataKinds, KindSignatures, ScopedTypeVariables #-}
{-# LANGUAGE FlexibleInstances, FlexibleContexts #-}
{-# LANGUAGE ApplicativeDo, BangPatterns, UnboxedTuples, TypeInType #-}
{-# LANGUAGE DerivingStrategies, GeneralizedNewtypeDeriving, DeriveAnyClass #-}
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.Default (Default, def)
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 qualified Data.Massiv.Array as A
import Data.Massiv.Array hiding (zipWith)
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 (Elevator, RGBA, Y, Pixel(..), ColorSpace)
import qualified Numeric.LinearAlgebra as LA
import qualified Prelude as P
import Prelude hiding (zipWith)
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 (EltRepr u Ix2) Ix2 a, Size u Ix2 a) => Show (Raster u p r c a) where
show (Raster a) = show . computeAs B $ extract' zeroIndex (r :. c) a
where (r :. c) = liftIndex (P.min 10) $ 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) = (\(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 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 (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)
(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.imapP_ 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)
(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) . 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 . sequenceA
{-# INLINE lmajority #-}
majo :: (Foldable t, Ord a) => t a -> a
majo = fst . g . f
where f = foldl' (\m a -> M.insertWith (+) a 1 m) M.empty
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 . sequenceA
{-# INLINE lminority #-}
mino :: (Foldable t, Ord a) => t a -> a
mino = fst . g . f
where f = foldl' (\m a -> M.insertWith (+) a 1 m) M.empty
g = L.foldl1' (\(a,c) (k,v) -> if c > v then (k,v) else (a,c)) . M.toList
{-# INLINE mino #-}
lvariance :: (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 <$> sequenceA rs)
where len = realToFrac $ length rs
avg ns = (\z -> realToFrac z / len) $ foldl' (+) 0 ns
f os@(n :| ns) = foldl' (\acc m -> acc + ((realToFrac m - av) ^ 2)) ((realToFrac n - av) ^ 2) ns / (len - 1)
where 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, Default a, Manifest u Ix2 a) => Raster u p r c a -> Raster DW p r c a
fsum (Raster a) = Raster $ mapStencil (neighbourhoodStencil f (Fill 0)) a
where f nw no ne we fo ea sw so se = nw + no + ne + we + fo + ea + sw + so + se
{-# INLINE fsum #-}
fproduct :: (Num a, Default a, Manifest u Ix2 a) => Raster u p r c a -> Raster DW p r c a
fproduct (Raster a) = Raster $ mapStencil (neighbourhoodStencil f (Fill 1)) a
where f nw no ne we fo ea sw so se = nw * no * ne * we * fo * ea * sw * so * se
fmonoid :: (Monoid a, Default a, Manifest u Ix2 a) => Raster u p r c a -> Raster DW p r c a
fmonoid (Raster a) = Raster $ mapStencil (neighbourhoodStencil f (Fill mempty)) a
where f nw no ne we fo ea sw so se = fo `mappend` nw `mappend` no `mappend` ne `mappend` we `mappend` ea `mappend` sw `mappend` so `mappend` se
fmean :: (Real a, Fractional b, Default a, Manifest u Ix2 a) => Raster u p r c a -> Raster DW p r c b
fmean = fmap (\n -> realToFrac n / 9) . fsum
{-# INLINE fmean #-}
fmax :: (Ord a, Default a, Manifest u Ix2 a) => Raster u p r c a -> Raster DW p r c a
fmax (Raster a) = Raster $ mapStencil (neighbourhoodStencil f Edge) a
where f nw no ne we fo ea sw so se = P.max nw . P.max no . P.max ne . P.max we . P.max fo . P.max ea . P.max sw $ P.max so se
{-# INLINE fmax #-}
fmin :: (Ord a, Default a, Manifest u Ix2 a) => Raster u p r c a -> Raster DW p r c a
fmin (Raster a) = Raster $ mapStencil (neighbourhoodStencil f Edge) a
where f nw no ne we fo ea sw so se = P.min nw . P.min no . P.min ne . P.min we . P.min fo . P.min ea . P.min sw $ P.min so se
{-# 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 (neighbourhoodStencil f Edge) 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 (neighbourhoodStencil f Continue) 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 (neighbourhoodStencil f Continue) 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 (neighbourhoodStencil f Continue) 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 (neighbourhoodStencil f Continue) 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 (neighbourhoodStencil linkage (Fill def)) 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 partStencil a
{-# INLINE fpartition #-}
partStencil :: (Eq a, Default a) => Stencil Ix2 a Corners
partStencil = makeStencil Reflect (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 (neighbourhoodStencil f Reflect) 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 (neighbourhoodStencil volume Reflect) 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) -> Border a -> Stencil Ix2 a b
neighbourhoodStencil f b = makeStencil b (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 Reflect (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 (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 (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 (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 (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 (neighbourhoodStencil f $ Fill (Drain 0)) 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 :) . 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