Gabor filter

.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,17 +101,19 @@ 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"
   #   os: osx
 
   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
 

hip.cabal

@@ -44,6 +44,7 @@ Library
                , repa            >= 3.2.1.1 && < 4
                , temporary       >= 1.1.1
                , vector          >= 0.10
+               , array
   Other-Extensions: BangPatterns
                   , ConstraintKinds
                   , CPP
@@ -67,6 +68,8 @@ Library
                  , Graphics.Image.Processing.Binary
                  , Graphics.Image.Processing.Complex
                  , Graphics.Image.Processing.Filter
+                 , Graphics.Image.Processing.Hough
+                 , Graphics.Image.Processing.Ahe
                  , Graphics.Image.Types
   Other-Modules:   Graphics.Image.ColorSpace.Binary
                  , Graphics.Image.ColorSpace.CMYK

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).
+--
+-- <<images/yield.jpg>>   <<images/yield_ahe.png>>
+--
+-- 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   
+

src/Graphics/Image/Processing/Filter.hs

@@ -14,7 +14,7 @@
 --
 module Graphics.Image.Processing.Filter
   ( -- * Filter
-    Filter
+    Filter (..)
   , applyFilter
   , Direction(..)
     -- * Gaussian
@@ -26,6 +26,14 @@ module Graphics.Image.Processing.Filter
     -- * Prewitt
   , prewittFilter
   , prewittOperator
+    -- * Laplacian
+  , laplacianFilter
+    -- * Laplacian of Gaussian
+  , logFilter
+    -- * Gaussian Smoothing
+  , gaussianSmoothingFilter
+    -- * Mean 
+  , meanFilter  
   ) where
 
 
@@ -46,8 +54,6 @@ data Direction
   = Vertical
   | Horizontal
 
-
-
 -- | Create a Gaussian Filter.
 --
 -- @since 1.5.3
@@ -95,7 +101,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 +124,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 +136,98 @@ 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 <https://homepages.inf.ed.ac.uk/rbf/HIPR2/log.htm>
+--
+-- <<images/yield.jpg>>   <<images/yield_laplacian.png>>
+--
+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 <https://homepages.inf.ed.ac.uk/rbf/HIPR2/log.htm>
+--
+-- <<images/yield.jpg>>   <<images/yield_log.png>>
+--
+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 <https://homepages.inf.ed.ac.uk/rbf/HIPR2/gsmooth.htm>
+-- 
+-- <<images/GSM_gsn_yield_IP.jpg>> <<images/GSM_gsn_yield_OP.png>> 
+--
+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 <http://homepages.inf.ed.ac.uk/rbf/HIPR2/mean.htm>
+-- 
+-- <<images/yield.jpg>>   <<images/yield_mean.png>>
+-- 
+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 #-}
+

src/Graphics/Image/Processing/Gabor.hs

@@ -0,0 +1,44 @@
+{-# LANGUAGE RankNTypes #-}
+{-# LANGUAGE FlexibleContexts #-}
+{-# LANGUAGE ScopedTypeVariables #-}
+{-# LANGUAGE BangPatterns #-} 
+
+module Graphics.Image.Processing.Gabor where
+
+import Control.Monad (forM_, when)
+import Control.Monad.ST
+import Data.STRef 
+
+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)
+import Data.Complex (Complex((:+)))
+
+gaborfn
+  :: (RealFloat p, Fractional p) => p -> p -> p -> p -> p -> p -> p -> Complex p  
+gaborfn λ θ ψ σ γ x y = exp ( (-0.5) * ((x'^2 + γ^2*y'^2) / (σ^2)) :+ 0) * exp ( 0 :+ (2*pi*(x'/λ+ψ)) )
+    where x' =  x * cos θ + y * sin θ
+          y' = -x * sin θ + y * cos θ
+          λ = 10.0
+          θ = pi
+          ψ = 0 
+          σ = 4.0
+          γ = 0.5
+
+gaborFilter :: (Array arr cs e, Array arr X e, Floating e, Fractional e) =>
+                   Border (Pixel cs e) -- ^ Border resolution technique.
+                -> Filter arr cs e
+gaborFilter !border =
+  Filter (correlate border gV' . correlate border gV)
+  where
+    !gV = compute $ (gabor / scalar weight)
+    !gV' = compute $ transpose gV
+    !gabor = makeImage (1, n) (promote gaborfn)
+    !weight = I.fold (+) 0 gabor
+    !n = 5
+
+
+