diff --git a/.travis.yml b/.travis.yml index 34744e2..494a4b8 100644 --- a/.travis.yml +++ b/.travis.yml @@ -36,12 +36,6 @@ matrix: include: # We grab the appropriate GHC and cabal-install versions from hvr's PPA. See: # https://github.com/hvr/multi-ghc-travis - - env: BUILD=cabal GHCVER=7.8.4 CABALVER=1.24 HAPPYVER=1.19.5 ALEXVER=3.1.7 - compiler: ": #GHC 7.8.4" - addons: - apt: - packages: [cabal-install-1.24,ghc-7.8.4,happy-1.19.5,alex-3.1.7] - sources: [hvr-ghc] - env: BUILD=cabal GHCVER=7.10.3 CABALVER=1.24 HAPPYVER=1.19.5 ALEXVER=3.1.7 compiler: ": #GHC 7.10.3" @@ -76,22 +70,18 @@ matrix: # The Stack builds. We can pass in arbitrary Stack arguments via the ARGS # variable, such as using --stack-yaml to point to a different file. - - env: BUILD=stack ARGS="--stack-yaml stack.lts-3.yaml" - compiler: ": #stack GHC-7.10.2 (lts-3.22)" + - env: BUILD=stack ARGS="--stack-yaml stack.lts-6.yaml" + compiler: ": #stack GHC-7.10.3 (lts-6.30)" addons: apt: packages: [libgmp-dev] - - env: BUILD=stack ARGS="--stack-yaml stack.lts-6.yaml" - compiler: ": #stack GHC-7.10.3 (lts-6.30)" + - env: BUILD=stack ARGS="--resolver lts-7" + compiler: ": #stack GHC-8.0.1 (lts-7.24)" addons: apt: packages: [libgmp-dev] - # - env: BUILD=stack ARGS="--resolver lts-7" - # compiler: ": #stack GHC-8.0.1 (lts-7)" - # addons: {apt: {packages: [libgmp-dev]}} - - env: BUILD=stack ARGS="--resolver lts-8.24" compiler: ": #stack GHC-8.0.2 (lts-8.24)" addons: @@ -111,9 +101,11 @@ matrix: packages: [libgmp-dev] # Nightly builds are allowed to fail - # - env: BUILD=stack ARGS="--resolver nightly" - # compiler: ": #stack nightly" - # addons: {apt: {packages: [libgmp-dev]}} + - env: BUILD=stack ARGS="--resolver nightly" + compiler: ": #stack nightly" + addons: + apt: + packages: [libgmp-dev] # - env: BUILD=stack ARGS="--resolver lts-8" # compiler: ": #stack 8.0.2 osx" @@ -121,7 +113,7 @@ matrix: allow_failures: - env: BUILD=cabal GHCVER=head CABALVER=head HAPPYVER=1.19.5 ALEXVER=3.1.7 - #- env: BUILD=stack ARGS="--resolver nightly" + - env: BUILD=stack ARGS="--resolver nightly" branches: except: @@ -152,11 +144,6 @@ before_install: mkdir -p $HOME/.cabal echo 'remote-repo: hackage.haskell.org:http://hackage.fpcomplete.com/' > $HOME/.cabal/config echo 'remote-repo-cache: $HOME/.cabal/packages' >> $HOME/.cabal/config - - # if [ "$CABALVER" != "1.16" ] - # then - # echo 'jobs: $ncpus' >> $HOME/.cabal/config - # fi ;; esac @@ -177,7 +164,7 @@ install: cabal --version travis_retry cabal update - # Get the list of packages from the stack.yaml file + # Set the list of packages we will build PACKAGES="$TRAVIS_BUILD_DIR" echo $PACKAGES diff --git a/hip.cabal b/hip.cabal index 015728d..9db522e 100644 --- a/hip.cabal +++ b/hip.cabal @@ -44,6 +44,8 @@ Library , repa >= 3.2.1.1 && < 4 , temporary >= 1.1.1 , vector >= 0.10 + , array + , random Other-Extensions: BangPatterns , ConstraintKinds , CPP @@ -67,6 +69,9 @@ Library , Graphics.Image.Processing.Binary , Graphics.Image.Processing.Complex , Graphics.Image.Processing.Filter + , Graphics.Image.Processing.Hough + , Graphics.Image.Processing.Ahe + , Graphics.Image.Processing.Noise , Graphics.Image.Types Other-Modules: Graphics.Image.ColorSpace.Binary , Graphics.Image.ColorSpace.CMYK diff --git a/images/GSM_gsn_yield_IP.jpg b/images/GSM_gsn_yield_IP.jpg new file mode 100644 index 0000000..0d68c12 Binary files /dev/null and b/images/GSM_gsn_yield_IP.jpg differ diff --git a/images/GSM_gsn_yield_OP.png b/images/GSM_gsn_yield_OP.png new file mode 100644 index 0000000..477cebc Binary files /dev/null and b/images/GSM_gsn_yield_OP.png differ diff --git a/images/yield_ahe.png b/images/yield_ahe.png new file mode 100644 index 0000000..0b4a8e4 Binary files /dev/null and b/images/yield_ahe.png differ diff --git a/images/yield_gray.png b/images/yield_gray.png new file mode 100644 index 0000000..0cf094d Binary files /dev/null and b/images/yield_gray.png differ diff --git a/images/yield_hough.png b/images/yield_hough.png new file mode 100644 index 0000000..4fb378c Binary files /dev/null and b/images/yield_hough.png differ diff --git a/images/yield_laplacian.png b/images/yield_laplacian.png new file mode 100644 index 0000000..18c87c8 Binary files /dev/null and b/images/yield_laplacian.png differ diff --git a/images/yield_log.png b/images/yield_log.png new file mode 100644 index 0000000..be60aab Binary files /dev/null and b/images/yield_log.png differ diff --git a/images/yield_mean.png b/images/yield_mean.png new file mode 100644 index 0000000..cbe5086 Binary files /dev/null and b/images/yield_mean.png differ diff --git a/images/yield_snp.png b/images/yield_snp.png new file mode 100644 index 0000000..ab0e856 Binary files /dev/null and b/images/yield_snp.png differ diff --git a/images/yield_unsharpMasking.png b/images/yield_unsharpMasking.png new file mode 100644 index 0000000..c865897 Binary files /dev/null and b/images/yield_unsharpMasking.png differ diff --git a/src/Graphics/Image/Processing/Ahe.hs b/src/Graphics/Image/Processing/Ahe.hs new file mode 100644 index 0000000..b50bc41 --- /dev/null +++ b/src/Graphics/Image/Processing/Ahe.hs @@ -0,0 +1,85 @@ +{-# LANGUAGE RankNTypes #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE BangPatterns #-} + +-- | Adaptive histogram equalization is used to improve contrast in images. +-- It adjusts image intensity in small regions (neighborhood) in the image. +module Graphics.Image.Processing.Ahe where + +import Control.Monad (forM_, when) +import Control.Monad.ST +import Data.STRef +import Debug.Trace (trace) + +import Prelude as P hiding (subtract) +import Graphics.Image.Processing.Filter +import Graphics.Image.Interface as I +import Graphics.Image +import Graphics.Image.Types as IP +import Graphics.Image.ColorSpace (X) + +-- | Supplementary function for applying border resolution and a general mask. +simpleFilter :: (Array arr cs e, Array arr X e) => Direction -> Border (Pixel cs e) -> Filter arr cs e +simpleFilter dir !border = + Filter (correlate border kernel) + where + !kernel = + case dir of + Vertical -> fromLists $ [ [ 0, -1, 0 ], [ -1, 4, -1 ], [ 0, -1, 0 ] ] + Horizontal -> fromLists $ [ [ 0, -1, 0 ], [ -1, 4, -1 ], [ 0, -1, 0 ] ] + +-- | 'ahe' operates on small 'contextual' regions of the image. It enhances the contrast of each +-- region and this technique works well when the distribution of pixel values is similar throughout +-- the image. +-- +-- The idea is to perform contrast enhancement in 'neighborhood region' of each pixel and the size +-- of the region is a parameter of the method. It constitutes a characteristic length scale: contrast +-- at smaller scales is enhanced, while contrast at larger scales is reduced (For general purposes, a size +-- factor of 5 tends to give pretty good results). +-- +-- <> <> +-- +-- Usage : +-- +-- >>> img <- readImageY VU "images/yield.jpg" +-- >>> input1 <- getLine +-- >>> input2 <- getLine +-- >>> let thetaSz = (P.read input1 :: Int) +-- >>> let distSz = (P.read input2 :: Int) +-- >>> let neighborhoodFactor = (P.read input2 :: Int) +-- >>> let aheImage :: Image VU RGB Double +-- >>> aheImage = ahe img thetaSz distSz neighborhoodFactor +-- >>> writeImage "images/yield_ahe.png" (toImageRGB aheImage) +-- +ahe + :: forall arr e cs . ( MArray arr Y Double, IP.Array arr Y Double, IP.Array arr Y Word16, MArray arr Y Word16, Array arr X Double) + => Image arr Y Double + -> Int -- ^ width of output image + -> Int -- ^ height of output image + -> Int -- ^ neighborhood size factor + -> Image arr Y Word16 +ahe image thetaSz distSz neighborhoodFactor = I.map (fmap toWord16) accBin + where + ip = applyFilter (simpleFilter Horizontal Edge) image -- Pre-processing (Border resolution) + widthMax, var1, heightMax, var2 :: Int + var1 = ((rows ip) - 1) + widthMax = ((rows ip) - 1) + var2 = ((cols ip) - 1) + heightMax = ((cols ip) - 1) + + accBin :: Image arr Y Word16 + accBin = runST $ -- Core part of the Algo begins here. + do arr <- I.new (thetaSz, distSz) -- Create a mutable image with the given dimensions. + forM_ [0 .. var1] $ \x -> do + forM_ [0 .. var2] $ \y -> do + rankRef <- newSTRef (0 :: Int) + let neighborhood a maxValue = filter (\a -> a >= 0 && a < maxValue) [a-5 .. a+5] + forM_ (neighborhood x var1) $ \i -> do + forM_ (neighborhood y var2) $ \j -> do + when (I.index ip (x, y) > I.index ip (i, j)) $ modifySTRef' rankRef (+1) + rank <- readSTRef rankRef + let px = ((rank * 255)) + I.write arr (x, y) (PixelY (fromIntegral px)) + freeze arr + diff --git a/src/Graphics/Image/Processing/Filter.hs b/src/Graphics/Image/Processing/Filter.hs index 5643b47..83eeb4e 100644 --- a/src/Graphics/Image/Processing/Filter.hs +++ b/src/Graphics/Image/Processing/Filter.hs @@ -14,7 +14,7 @@ -- module Graphics.Image.Processing.Filter ( -- * Filter - Filter + Filter (..) , applyFilter , Direction(..) -- * Gaussian @@ -26,6 +26,16 @@ module Graphics.Image.Processing.Filter -- * Prewitt , prewittFilter , prewittOperator + -- * Laplacian + , laplacianFilter + -- * Laplacian of Gaussian + , logFilter + -- * Gaussian Smoothing + , gaussianSmoothingFilter + -- * Mean + , meanFilter + -- * Unsharp Masking + , unsharpMaskingFilter ) where @@ -46,8 +56,6 @@ data Direction = Vertical | Horizontal - - -- | Create a Gaussian Filter. -- -- @since 1.5.3 @@ -95,7 +103,7 @@ sobelFilter dir !border = , [ 0, 0, 0 ] , [ 1, 2, 1 ] ] Horizontal -> fromLists $ [ [ -1, 0, 1 ] - , [ -2, 0, 1 ] + , [ -2, 0, 2 ] , [ -1, 0, 1 ] ] {-# INLINE sobelFilter #-} @@ -118,8 +126,6 @@ sobelOperator !img = sqrt (sobelX ^ (2 :: Int) + sobelY ^ (2 :: Int)) {-# INLINE sobelOperator #-} - - prewittFilter :: (Array arr cs e, Array arr X e) => Direction -> Border (Pixel cs e) -> Filter arr cs e prewittFilter dir !border = @@ -132,10 +138,121 @@ prewittFilter dir !border = {-# INLINE prewittFilter #-} - prewittOperator :: (Array arr cs e, Array arr X e, Floating e) => Image arr cs e -> Image arr cs e prewittOperator !img = sqrt (prewittX ^ (2 :: Int) + prewittY ^ (2 :: Int)) where !prewittX = applyFilter (prewittFilter Horizontal Edge) img !prewittY = applyFilter (prewittFilter Vertical Edge) img {-# INLINE prewittOperator #-} + +-- |The Laplacian of an image highlights regions of rapid intensity change +-- and is therefore often used for edge detection. It is often applied to an +-- image that has first been smoothed with something approximating a +-- Gaussian smoothing filter in order to reduce its sensitivity to noise. +-- More info about the algo at +-- +-- <> <> +-- +laplacianFilter :: (Array arr cs e, Array arr X e) => + Border (Pixel cs e) -> Filter arr cs e +laplacianFilter !border = + Filter (correlate border kernel) + where + !kernel = fromLists $ [ [ -1, -1, -1 ] -- Unlike the Sobel edge detector, the Laplacian edge detector uses only one kernel. + , [ -1, 8, -1 ] -- It calculates second order derivatives in a single pass. + , [ -1, -1, -1 ]] -- This is the approximated kernel used for it. (Includes diagonals) +{-# INLINE laplacianFilter #-} + +-- | 'Laplacian of Gaussian' (LOG) filter is a two step process of smoothing +-- an image before applying some derivative filter on it. This comes in +-- need for reducing the noise sensitivity while working with noisy +-- datasets or in case of approximating second derivative measurements. +-- +-- The LoG operator takes the second derivative of the image. Where the image +-- is basically uniform, the LoG will give zero. Wherever a change occurs, the LoG will +-- give a positive response on the darker side and a negative response on the lighter side. +-- More info about the algo at +-- +-- <> <> +-- +logFilter :: (Array arr cs e, Array arr X e) => + Border (Pixel cs e) -> Filter arr cs e +logFilter !border = + Filter (correlate border kernel) + where + !kernel = fromLists $ [ [ 0, 1, 1, 2, 2, 2, 1, 1, 0 ] + , [ 1, 2, 4, 5, 5, 5, 4, 2, 1 ] + , [ 1, 4, 5, 3, 0, 3, 5, 4, 1 ] + , [ 2, 5, 3, -12, -24, -12, 3, 5, 2 ] + , [ 2, 5, 0, -24, -40, -24, 0, 5, 2 ] + , [ 2, 5, 3, -12, -24, -12, 3, 5, 2 ] + , [ 1, 4, 5, 3, 0, 3, 5, 4, 1 ] + , [ 1, 2, 4, 5, 5, 5, 4, 2, 1 ] + , [ 0, 1, 1, 2, 2, 2, 1, 1, 0 ] ] +{-# INLINE logFilter #-} + +-- | The Gaussian smoothing operator is a 2-D convolution operator that is used to +-- `blur' images and remove detail and noise. The idea of Gaussian smoothing is to use +-- this 2-D distribution as a `point-spread' function, and this is achieved by convolution. +-- Since the image is stored as a collection of discrete pixels we need to produce a +-- discrete approximation to the Gaussian function before we can perform the convolution. +-- More info about the algo at +-- +-- <> <> +-- +gaussianSmoothingFilter :: (Fractional e, Array arr cs e, Array arr X e) => + Border (Pixel cs e) -> Filter arr cs e +gaussianSmoothingFilter !border = + Filter (I.map (/ 273) . correlate border kernel) + where + !kernel = fromLists $ [[ 1, 4, 7, 4, 1 ] -- Discrete approximation to the Gaussian function. + ,[ 4, 16, 26, 16, 4 ] -- 273 is the sum of all values in the mask and hence used in rescaling. + ,[ 7, 26, 41, 26, 7 ] + ,[ 4, 16, 26, 16, 4 ] + ,[ 1, 4, 7, 4, 1 ]] + +{-# INLINE gaussianSmoothingFilter #-} + + +-- | The mean filter is a simple sliding-window spatial filter that replaces the +-- center value in the window with the average (mean) of all the pixel values in +-- the window. The window, or kernel, can be any shape, but this one uses the most +-- common 3x3 square kernel. +-- More info about the algo at +-- +-- <> <> +-- +meanFilter :: (Fractional e, Array arr cs e, Array arr X e) => + Border (Pixel cs e) -> Filter arr cs e +meanFilter !border = + Filter (I.map (/ 9) . correlate border kernel) + where + !kernel = fromLists $[ [ 1, 1, 1 ] -- Replace each pixel with the mean value of its neighbors, including itself. + , [ 1, 1, 1 ] + , [ 1, 1, 1 ]] + +{-# INLINE meanFilter #-} + +-- | The unsharp-masking filter is a sharpening operator which derives its name from +-- the fact that it enhances edges (and other high frequency components in an image) +-- via a procedure which subtracts an unsharp, or smoothed, version of an image from +-- the original image. It is commonly used in the photographic and printing industries +-- for crispening edges. +-- More info about the algo at +-- +-- <> <> +-- +unsharpMaskingFilter :: (Fractional e, Array arr cs e, Array arr X e) => + Border (Pixel cs e) -> Filter arr cs e +unsharpMaskingFilter !border = + Filter (I.map (/256) . correlate border kernel) + where + !kernel = fromLists $ [[ -1, -4, -6, -4, -1 ] + ,[ -4, -16, -24, -16, -4 ] -- Uses negative image to create a mask of the original image. + ,[ -6, -24, 476, -24, -6 ] -- The unsharped mask is then combined with the positive (original) image. + ,[ -4, -16, -24, -16, -4 ] -- So, the resulting image is less blurry, i.e clearer. + ,[ -1, -4, -6, -4, -1 ]] + +{-# INLINE unsharpMaskingFilter #-} + + diff --git a/src/Graphics/Image/Processing/Hough.hs b/src/Graphics/Image/Processing/Hough.hs new file mode 100644 index 0000000..089c5b9 --- /dev/null +++ b/src/Graphics/Image/Processing/Hough.hs @@ -0,0 +1,109 @@ +{-# LANGUAGE RankNTypes #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE ScopedTypeVariables #-} + +-- | Hough Transform is used as a part of feature extraction in images. +-- It is a tool that makes it far easier to identify straight lines in +-- the source image, whatever their orientation. +module Graphics.Image.Processing.Hough where + +import Control.Monad (forM_, when) +import qualified Data.Foldable as F (maximum) +import Data.Array +import Data.Array.ST (newArray, writeArray, readArray, runSTArray) + +import Prelude as P hiding (subtract) +import Graphics.Image +import Graphics.Image.Interface as I +import Graphics.Image.Types as IP + +-- | Some helper functions : +-- | Trivial function for subtracting co-ordinate pairs +sub :: Num x => (x, x) -> (x, x) -> (x, x) +sub (x1, y1) (x2, y2) = (x1 - x2, y1 - y2) + +-- | Compute the sum of squares or dot product of a given pair of co-ordinates +dotProduct :: Num x => (x, x) -> (x, x) -> x +dotProduct (x1, y1) (x2, y2) = (x1 * x2) + (y1 * y2) + +-- | Conversion of pair fromIntegral +fromIntegralP :: (Integral x, Num y) => (x, x) -> (y, y) +fromIntegralP (x1, y1) = (fromIntegral x1, fromIntegral y1) + +-- | Compute magnitude +mag :: Floating x => (x, x) -> x +mag x = sqrt (dotProduct x x) + +-- | 'hough' computes the Linear Hough Transform and maps each point in the target image, ​ (ρ, θ) ​ +-- to the average color of the pixels on the corresponding line of the source image ​(x,y) ​- space, +-- where the line corresponds to points of the form ​(xcosθ + ysinθ = ρ(rho)). +-- +-- The idea is that where there is a straight line in the original image, it corresponds to a +-- bright (or dark, depending on the color of the background field) spot; by applying a suitable +-- filter to the results of the transform, it is possible to extract the locations of the lines in the original image. +-- +-- <> <> +-- +-- Usage : +-- +-- >>> frog <- readImageRGB VU "yield.jpg" +-- >>> input1 <- getLine +-- >>> input2 <- getLine +-- >>> let thetaSz = (P.read input1 :: Int) +-- >>> let distSz = (P.read input2 :: Int) +-- >>> let houghImage :: Image VU RGB Double +-- >>> houghImage = hough frog thetaSz distSz +-- >>> writeImage "test.png" houghImage +-- +hough + :: forall arr . ( MArray arr Y Double, IP.Array arr Y Double, IP.Array arr Y Word8) + => Image arr Y Double + -> Int + -> Int + -> Image arr Y Word8 +hough image thetaSz distSz = I.map (fmap toWord8) hImage + where + widthMax, xCtr, heightMax, yCtr :: Int + widthMax = ((rows image) - 1) + xCtr = (widthMax `div` 2) + heightMax = ((cols image) - 1) + yCtr = (heightMax `div` 2) + + slope :: Int -> Int -> (Double, Double) + slope x y = + let PixelY orig = I.index image (x, y) + PixelY x' = I.index image (min (x+1) widthMax, y) + PixelY y' = I.index image (x, min (y+1) heightMax) + in (orig - x', orig - y') + + slopeMap :: [ ((Int, Int), (Double, Double)) ] + slopeMap = [ ((x, y), slope x y) | x <- [0 .. widthMax], y <- [0 .. heightMax] ] + + distMax :: Double -- Compute Maximum distance + distMax = (sqrt . fromIntegral $ (heightMax + 1) ^ (2 :: Int) + (widthMax + 1) ^ (2 :: Int)) / 2 + + accBin = runSTArray $ -- Core part of Algo begins here. Working in a safe way with a mutable array. + do arr <- newArray ((0, 0), (thetaSz, distSz)) (0 :: Double) -- Build a new array, with every element initialised to the supplied value. + forM_ slopeMap $ \((x, y), gradient) -> do + let (x', y') = fromIntegralP $ (xCtr, yCtr) `sub` (x, y) + when (mag gradient > 0) $ + forM_ [0 .. thetaSz] $ \theta -> do + let theta_ = + fromIntegral theta * 360 / fromIntegral thetaSz / 180 * + pi :: Double + distance = cos theta_ * x' + sin theta_ * y' -- (ρ(rho) = xcosθ + ysinθ) + distance_ = truncate $ distance * fromIntegral distSz / distMax -- returns the nearest integer + idx = (theta, distance_) + when (distance_ >= 0 && distance_ < distSz) $ + do old <- readArray arr idx -- read an element at 'idx' from mutable array 'arr' + writeArray arr idx (old + 1) + return arr + + maxAcc = F.maximum accBin + hTransform (x, y) = + let l = 255 - truncate ((accBin ! (x, y)) / maxAcc * 255) -- pixel generating function + in PixelY l + + hImage :: Image arr Y Word8 + hImage = makeImage (thetaSz, distSz) hTransform + diff --git a/src/Graphics/Image/Processing/Noise.hs b/src/Graphics/Image/Processing/Noise.hs new file mode 100644 index 0000000..0d9d1bd --- /dev/null +++ b/src/Graphics/Image/Processing/Noise.hs @@ -0,0 +1,100 @@ +{-# LANGUAGE RankNTypes #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE BangPatterns #-} + +module Graphics.Image.Processing.Noise where + +import Control.Monad (forM_) +import Control.Monad.ST +import System.Random +import Data.Random.Normal +import Data.Array.MArray + +import Prelude as P hiding (subtract) +import Graphics.Image.Interface as I +import Graphics.Image +import Graphics.Image.Types as IP + +-- | Helper function for generating a list of random co-ordinates. +randomCoords :: StdGen -> Int -> Int -> [(Int,Int)] +randomCoords a width height = (rnx1, rny1) : randomCoords g2 width height + where + (rnx1, g1) = randomR (0, width) a + (rny1, g2) = randomR (0, height) g1 + +-- | Salt and pepper noise or impulse noise is a form of noise seen on images. +-- It is mainly caused by sharp and sudden disturbances in the image signal. +-- +-- 'saltAndPepper' generates this particular type of noise by introducing a sparse +-- distribution of white and black pixels in the input image. The level or intensity +-- of the noise to be introduced is a parameter of this method and is scaled +-- between 0 and 1, that is the input Noise Intensity has a domain : (0, 1). +-- +-- <> <> +-- +-- Usage : +-- +-- >>> img <- readImageY VU "images/yield.jpg" +-- >>> input1 <- getLine +-- >>> g <- newStdGen +-- >>> let noiseLevel = (P.read input1 :: Float) +-- >>> let snpImage :: Image VU Y Double +-- >>> snpImage = saltAndPepper img noiseLevel g +-- >>> writeImage "images/yield_snp.png" snpImage +-- +saltAndPepper + :: forall arr e cs . (IP.MArray arr Y Double, IP.Array arr Y Double) + :: forall arr e cs . (MArray arr Y Double, IP.Array arr Y Double) + => Image arr Y Double + -> Float -- ^ Noise Intensity -> Domain : (0, 1) + -> StdGen -- ^ Instance of RandomGen + -> Image arr Y Double +saltAndPepper image noiseLevel = accBin + where + widthMax, heightMax, noiseIntensity :: Int + widthMax = ((rows image) - 1) + heightMax = ((cols image) - 1) + noiseIntensity = round (noiseLevel * (fromIntegral widthMax) * (fromIntegral heightMax)) + + accBin :: StdGen -> Image arr Y Double + accBin g = runST $ + do arr <- I.thaw image + let coords = take (noiseIntensity + 1) (randomCoords g widthMax heightMax) + forM_ coords $ \i -> do + let a :: Int + a = uncurry (+) i + if (a `mod` 2 == 0) + then do let px = 0 + I.write arr i px + else do let px = 1.0 + I.write arr i px + I.freeze arr + +gaussianNoise + :: forall arr e cs . (IP.MArray arr Y Double, IP.Array arr Y Double) + => Image arr Y Double + -> Float -- ^ Mean + -> Float -- ^ Sigma + -> StdGen -- ^ Instance of RandomGen + -> Image arr Y Double +gaussianNoise image mean sigma g = accBin + where + widthMax, heightMax :: Int + widthMax = ((rows image) - 1) + heightMax = ((cols image) - 1) + samples = normals' (mean, sigma) g + new = newListArray (widthMax, heightMax) samples + + accBin :: [(Int,Int)] -> Image arr Y Double + accBin new = runST $ + do arr <- I.thaw image + forM_ [0 .. widthMax] $ \x -> do + forM_ [0 .. heightMax] $ \y -> do + let value = (fromIntegral (readArray new (x, y) )) + let px = ( (I.index image (x, y)) + (PixelY value) ) + I.write arr (x, y) px + I.freeze arr + + + freeze arr