diff --git a/CHANGELOG.md b/CHANGELOG.md
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -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.
diff --git a/README.md b/README.md
--- a/README.md
+++ b/README.md
@@ -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
diff --git a/bench/Bench.hs b/bench/Bench.hs
--- a/bench/Bench.hs
+++ b/bench/Bench.hs
@@ -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
+  ]
diff --git a/data/1024x1024.tif b/data/1024x1024.tif
deleted file mode 100644
# file too large to diff: data/1024x1024.tif
diff --git a/data/256x256.tif b/data/256x256.tif
deleted file mode 100644
Binary files a/data/256x256.tif and /dev/null differ
diff --git a/data/gray.png b/data/gray.png
new file mode 100644
Binary files /dev/null and b/data/gray.png differ
diff --git a/data/spectrum.png b/data/spectrum.png
new file mode 100644
Binary files /dev/null and b/data/spectrum.png differ
diff --git a/lib/Geography/MapAlgebra.hs b/lib/Geography/MapAlgebra.hs
--- a/lib/Geography/MapAlgebra.hs
+++ b/lib/Geography/MapAlgebra.hs
@@ -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
diff --git a/mapalgebra.cabal b/mapalgebra.cabal
--- a/mapalgebra.cabal
+++ b/mapalgebra.cabal
@@ -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
diff --git a/test/Test.hs b/test/Test.hs
--- a/test/Test.hs
+++ b/test/Test.hs
@@ -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
+-}