mapalgebra 0.1.0 → 0.1.1

raw patch · 10 files changed

+147/−54 lines, 10 filesdep ~basebinary-added

Dependency ranges changed: base

Files

CHANGELOG.md +6/−0
README.md +20/−3
bench/Bench.hs +19/−11
data/1024x1024.tif too large to diff
data/256x256.tif binary
data/gray.png binary
data/spectrum.png binary
lib/Geography/MapAlgebra.hs +60/−22
mapalgebra.cabal +9/−10
test/Test.hs +33/−8

CHANGELOG.md view

@@ -1,3 +1,9 @@+## 0.1.1++- Added `histogram` and `breaks` functions for automatically producing+  "break point" lists for functions like `greenRed`. Currently works+  only for `Word8` Rasters.+ ## 0.1.0  - Initial release.

README.md view

@@ -21,7 +21,7 @@ can hold any kind of data, and are aware of their projection and dimensions at the type level. This means that imagery of different size or projection are considered completely different types, which prevents an entire class-of potential bugs.+of bugs.  `Raster`s have types signatures like this: @@ -39,6 +39,8 @@  -- | A "windowed" Raster of an ADT, the result of some Focal Operation. -- Waiting to be evaluated by the `strict` function.+--+-- A generic `p` means we don't care about Projection here. Raster DW p 1024 1024 (Maybe Double) ``` @@ -93,9 +95,24 @@  -- | Given a list of "breaks", forms a colour ramp to be passed -- to `classify`.-greenRed :: Ord k => [k] -> Map k (Pixel RGBA Word8)+spectrum :: Ord k => [k] -> Map k (Pixel RGBA Word8) ``` +We can generate the breaks via `histogram` and `breaks` (currently supports `Word8` only):++```haskell+-- | Input here is `Raster S`, meaning this image was freshly read,+-- or has already been fully computed with `strict`.+colourIt :: Raster S p r c Word8 -> Raster S p r c (Pixel RGBA Word8)+colourIt r = strict S . classify invisible cm $ lazy r+  where cm = spectrum . breaks $ histogram r+```++Before-and-after:++![original image](data/gray.png)+![colour via spectrum ramp](data/spectrum.png)+ ### Local Operations  All Local Operations defined in *GIS and Cartographic Modeling* are available.@@ -142,7 +159,7 @@ ## Future Work  - Projection handling at IO time-- Histograms for colour ramp generation+- Histogram generation for more data types - Reprojections - Extended neighbourhoods for Focal Ops - Upsampling and Downsampling

bench/Bench.hs view

@@ -5,7 +5,7 @@ import           Criterion.Main import           Data.List.NonEmpty (NonEmpty(..)) import           Data.Massiv.Array as A hiding (zipWith)-import           Data.Monoid ((<>), Sum(..))+import           Data.Monoid ((<>)) import qualified Data.Vector as V import qualified Data.Vector.Unboxed as U import           GHC.TypeLits@@ -38,6 +38,7 @@     , conversions img     , focalOps img imgF     , compositeOps a512 -- a1024 a2048 a4096+    , histograms img     ]  creation :: Benchmark@@ -74,9 +75,11 @@ localOps (RGBARaster r g b _) (RGBARaster rF gF bF _) img = bgroup "Local Operations"   [ bench "fmap (+ 17) . lazy" $ nf (_array . strict S . fmap (+ 17) . lazy) img   , bench "zipWith (+)"    $ nf (_array . strict S . zipWith (+) r) g-  , bench "zipWith (/)"    $ nf (_array . strict S . zipWith (/) doubles) doubles+  , bench "zipWith (/)"    $ nf (_array . strict S . zipWith (/) rF) gF   , bench "(+)"            $ nf (_array . strict S . (+ lazy r)) (lazy g)-  , bench "(/)"            $ nf (_array . strict S . (/ lazy doubles)) (lazy doubles)+  , bench "(/)"            $ nf (_array . strict S . (/ lazy rF)) (lazy gF)+  , bench "(.+)"           $ nf (\g' -> computeAs S $ (_array r) .+ g') (_array g)+  , bench "(./)"           $ nf (\g' -> computeAs S $ (_array rF) ./ g') (_array gF)   , bench "lmax"           $ nf (_array . strict S . lmax img) img   , bench "lmin"           $ nf (_array . strict S . lmin img) img   , bench "lmean (Word8)"  $ nf (_array . strict S . lmean @Word8 @Double) rs@@ -108,7 +111,7 @@                  [ bench "512"    $ nf (_array . strict S . fsum) img                  -- , bench "46500"  $ nf (_array . strict S . fsum) huge                  -- , bench "maybing" $ nf (_array . strict B . fmap Just . lazy) img-                 , bench "nodata" $ nf (_array . strict S . nodatafsum) img+                 -- , bench "nodata" $ nf (_array . strict S . nodatafsum) img                  ]                , bgroup "fmean"                  [ bench "Word8"  $ nf (_array . strict S . fmean @Word8 @Double) img@@ -193,8 +196,8 @@ vectorB' :: V.Vector Int -> Raster B p 256 256 Int vectorB' = fromRight . fromVector Par -doubles :: Raster U p 512 512 Double-doubles = fromRight . fromVector Par $ U.fromList ([1..262144] :: [Double])+-- doubles :: Raster U p 512 512 Double+-- doubles = fromRight . fromVector Par $ U.fromList ([1..262144] :: [Double])  zing :: LA.Matrix Double zing = LA.matrix 3 [ -0.5, -0.5, 1@@ -243,8 +246,13 @@         denom = nir + (6 * lazy r) - (7.5 * lazy b) + 1 {-# INLINE evi #-} -nodatafsum :: Raster S p 512 512 Word8 -> Raster DW p 512 512 Word8-nodatafsum = fmap (maybe 0 getSum) . fmonoid . strict B . fmap check . lazy-  where check 0 = Nothing-        check n = Just $ Sum n-{-# INLINE nodatafsum #-}+-- nodatafsum :: Raster S p 1753 1760 Double -> Raster DW p 1753 1760 Double+-- nodatafsum = fmap (maybe 0 ((/9) . getSum)) . fmonoid . strict B . fmap check . lazy+--   where check 0 = Nothing+--         check n = Just $ Sum n+-- {-# INLINE nodatafsum #-}++histograms :: Raster S p 512 512 Word8 -> Benchmark+histograms r = bgroup "Histogram Algorithms"+  [ bench "Storable Super Mutable" $ nf (_histogram . histogram) r+  ]

− data/1024x1024.tif

file too large to diff

− data/256x256.tif

binary file changed (263242 → absent bytes)

+ data/gray.png view

binary file changed (absent → 47430 bytes)

+ data/spectrum.png view

binary file changed (absent → 122990 bytes)

lib/Geography/MapAlgebra.hs view

@@ -83,19 +83,27 @@   -- *** Creation   , constant, fromFunction, fromVector, fromRGBA, fromGray   -- *** Colouring-  -- | The `M.Map`s here can be used with `classify` to-  --   transform /ranges/ of values into certain colours in \(\mathcal{O}(n\log(n))\).-  --   Each Map-generating function (like `greenRed`) creates a "colour ramp" of 10 colours. So, it expects-  --   to be given a list of 10 "break points" which become the Map's keys. Any less than 10 will result-  --   in the later colours not being used. Any more than 10 will be ignored. The list of break points is-  --   assumed to be sorted.-  --   `invisible` can be used as the default value to `classify`, to make invisible any value that falls outside-  --   the range of the Maps.+  -- | To colour a `Raster`:   --+  -- 1. Fully evaluate it with `strict` `S`.+  -- 2. Use `histogram` to analyse the spread of its values. Currently only `Word8` is supported.+  -- 3. Use `breaks` on the `Histogram` to generate a list of "break points".+  -- 4. Pass the break points to a function like `greenRed` to generate a "colour ramp" (a `M.Map`).+  -- 5. Use this ramp with `classify` to colour your `Raster` in \(\mathcal{O}(n\log(n))\).+  --    `invisible` can be used as the default colour for values that fall outside the range+  --    of the Map.+  --+  -- @+  -- colourIt :: Raster D p r c Word8 -> Raster S p r c (Pixel RGBA Word8)+  -- colourIt (strict S -> r) = strict S . classify invisible cm $ lazy r+  --   where cm = blueGreen . breaks $ histogram r+  -- @+  --   -- If you aren't interested in colour but still want to render your `Raster`,   -- consider `grayscale`. Coloured `Raster`s can be unwrapped with `_array` and then   -- output with functions like `writeImage`.   , grayscale+  , Histogram(..), histogram, breaks   , invisible   , greenRed, spectrum, blueGreen, purpleYellow, brownBlue   , grayBrown, greenPurple, brownYellow, purpleGreen, purpleRed@@ -104,9 +112,16 @@   -- @Raster u p r c (Pixel RGBA Word8)@. Then unwrap it with `_array` and output   -- with something like `writeImage`.   ---  -- For quick debugging, you can visualize a `Raster` with `display`.+  -- For quick debugging, you can visualize a `Raster` with `displayImage`:+  --+  -- @+  -- raster :: Raster D p 512 512 Word8+  --+  -- >>> displayImage . _array . strict S . grayscale+  -- @   , writeImage, writeImageAuto-  , png, display+  , displayImage+  , png   -- ** Projections   , Projection(..)   , reproject@@ -244,10 +259,16 @@   , fvolume, fgradient, faspect, faspect', fdownstream, fupstream   -- * Utilities   , leftPseudo, tau+  -- * Massiv Re-exports+  -- | For your convenience.+  , 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@@ -268,6 +289,7 @@ 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)@@ -531,11 +553,11 @@   mb <- A.unsafeNew sz   mg <- A.unsafeNew sz   ma <- A.unsafeNew sz-  flip A.imapP_ arr $ \ix (PixelRGBA r g b a) -> do-    A.unsafeWrite mr ix r-    A.unsafeWrite mg ix g-    A.unsafeWrite mb ix b-    A.unsafeWrite ma ix a+  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@@ -564,7 +586,7 @@ -- | Construct a colour ramp. -- ramp :: Ord k => [(Word8, Word8, Word8)] -> [k] -> M.Map k PixelRGBA8 ramp :: Ord k => [(Word8, Word8, Word8)] -> [k] -> M.Map k (Pixel RGBA Word8)-ramp colours breaks = M.fromList . P.zip breaks $ P.map (\(r,g,b) -> PixelRGBA r g b maxBound) colours+ramp colours bs = M.fromList . P.zip bs $ P.map (\(r,g,b) -> PixelRGBA r g b maxBound) colours {-# INLINE ramp #-}  -- | From page 32 of /Cartographer's Toolkit/.@@ -633,12 +655,6 @@ grayscale = fmap PixelY {-# INLINE grayscale #-} --- | View a `Raster` as grayscale with the default image viewer of your OS.------ For more direct control, consider `displayImage` from 'massiv-io'.-display :: (Functor (Raster u p r c), Load u Ix2 (Pixel Y a), Elevator a) => Raster u p r c a -> IO ()-display = displayImage . computeAs S . _array . grayscale- -- | Render a PNG-encoded `BL.ByteString` from a coloured `Raster`. -- Useful for returning a `Raster` from a webserver endpoint. png :: (Source u Ix2 (Pixel cs a), ColorSpace cs a) => Raster u p r c (Pixel cs a) -> BL.ByteString@@ -1245,3 +1261,25 @@           SouthWest -> acc + 32           South     -> acc + 64           SouthEast -> acc + 128++-- | A count of `Word8` values across some `Raster`.+newtype Histogram = Histogram { _histogram :: VS.Vector Word } deriving (Eq, Show)++-- | Given a `Raster` of byte data, efficiently produce a `Histogram` that+-- describes value counts across the image. To be passed to `breaks`.+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 #-}++-- | Given a `Histogram`, produce a list of "colour breaks" as needed by+-- functions like `greenRed`.+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  -- Yes, 11 and not 10.+        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

mapalgebra.cabal view

@@ -1,15 +1,15 @@--- This file has been generated from package.yaml by hpack version 0.20.0.+-- This file has been generated from package.yaml by hpack version 0.28.2. -- -- see: https://github.com/sol/hpack ----- hash: 8c5d8cd118654ac1b12534eabd521ebaabddb34bc070fb3780ac606cf3d2be9a+-- hash: 5424873aec53550165fe1a8e3fe8c7b2105ac12764efc49bc9608080d70c79df  name:           mapalgebra-version:        0.1.0+version:        0.1.1 synopsis:       Efficient, polymorphic Map Algebra. description:    Efficient, polymorphic Map Algebra.                 .-                This library is an implementation of /Map Algebra/ as described in the book /GIS and Cartographic Modeling/ by Dana Tomlin. The fundamental type, the `Raster`, is typesafe. Rasters of different size and projection are considered different types, and so cannot be combined in any way.+                This library is an implementation of /Map Algebra/ as described in the book /GIS and Cartographic Modeling/ by Dana Tomlin. The fundamental type, the `Raster`, is typesafe. Rasters of different size and projection are considered different types, and so cannot be combined mistakenly in any way.                 .                 Also featured are op fusion (i.e. "lazy Rasters"), extremely fast Focal Operations, and typesafe NoData handling. Please see the main module for a more detailed introduction. category:       Geography@@ -21,13 +21,12 @@ license-file:   LICENSE build-type:     Simple cabal-version:  >= 1.10- extra-source-files:     CHANGELOG.md-    data/1024x1024.tif-    data/256x256.tif     data/512x512.tif+    data/gray.png     data/gray512.tif+    data/spectrum.png     README.md  library@@ -35,7 +34,7 @@       lib   ghc-options: -fwarn-unused-imports -fwarn-unused-binds -fwarn-name-shadowing -fwarn-unused-matches -fwarn-incomplete-patterns -Wincomplete-uni-patterns   build-depends:-      base >=4.9 && <4.12+      base >=4.10 && <4.12     , bytestring     , containers     , data-default@@ -57,7 +56,7 @@   build-depends:       HUnit-approx >=1.1 && <1.2     , QuickCheck-    , base >=4.9 && <4.12+    , base >=4.10 && <4.12     , bytestring     , containers     , data-default@@ -79,7 +78,7 @@       bench   ghc-options: -fwarn-unused-imports -fwarn-unused-binds -fwarn-name-shadowing -fwarn-unused-matches -fwarn-incomplete-patterns -Wincomplete-uni-patterns -threaded -O2 -with-rtsopts=-N   build-depends:-      base >=4.9 && <4.12+      base >=4.10 && <4.12     , bytestring     , containers     , criterion >=1.1 && <1.4

test/Test.hs view

@@ -6,8 +6,9 @@ module Main ( main ) where  import           Data.Int+import           Data.List (sort) import           Data.List.NonEmpty (NonEmpty(..))-import           Data.Massiv.Array as A+import           Data.Massiv.Array (getComp, index') import qualified Data.Set as S import qualified Data.Vector as V import qualified Data.Vector.Unboxed as U@@ -25,10 +26,12 @@ ---  main :: IO ()-main = defaultMain suite+main = do+  img <- fromRight <$> fromGray "data/gray512.tif"+  defaultMain $ suite img -suite :: TestTree-suite = testGroup "Unit Tests"+suite :: Raster S p 512 512 Word8 -> TestTree+suite r = testGroup "Unit Tests"   [ testGroup "Raster Creation"     [ testCase "constant (256x256)"     $ length (lazy small) @?= 65536     , testCase "constant (2^16 x 2^16)" $ length lazybig @?= 4294967296@@ -96,8 +99,19 @@       , testCase "3x3 Flat" fupstreamFlat       ]     ]+  , testGroup "Histograms"+    [ testCase "Total Sum"     $ VS.sum (_histogram $ histogram r) @?= 262144+    , testCase "10 Breaks"     $ length (breaks $ histogram r) @?= 10+    , testCase "Sorted Breaks" $ do+        let bs = breaks $ histogram r+        sort bs @?= bs+    ]   ] +fromRight :: Either a b -> b+fromRight (Right b) = b+fromRight _ = error "Was Left"+ one :: Raster P p 7 7 Word one = constant P Seq 1 @@ -116,6 +130,9 @@ -- bog :: Raster P p 65536 65536 Word8 -- bog = constant P Par 10 +-- indices :: Raster S p 512 512 Word8+-- indices = fromFunction S Par (\(r :. c) -> fromIntegral $ r + c)+ fileRGBA :: IO (Either String (RGBARaster p 512 512 Word8)) fileRGBA = fromRGBA "data/512x512.tif" @@ -171,10 +188,6 @@         expected :: Raster U p 3 3 Double         expected = fromRight . fromVector Seq $ U.fromList [ 0, 0.5, 0, 0.5, 2, 0.5, 0, 0.5, 0 ] -fromRight :: Either a b -> b-fromRight (Right b) = b-fromRight _ = error "Was Left"- fpartitionTest :: Assertion fpartitionTest = actual @?= expected   where expected :: Raster B p 2 2 Corners@@ -336,3 +349,15 @@ fupstreamPeak = (flip index' (1 :. 1) $ _array actual) @?= Drain 0   where actual :: Raster S p 3 3 Drain         actual = strict S . fupstream . strict S . fdownstream . fromRight . fromVector Seq $ U.fromList ([1,1,1,1,3,1,1,1,1] :: [Double])++{-+hists :: Raster S p 512 512 Word8 -> Assertion+hists r = do+  u  <- histU r+  um <- histMut r+  histU' r @?= u+  histMut' r @?= um+  histMut'' r @?= um+  histMut'' r @?= u+  um @?= u+-}