{-# LANGUAGE RankNTypes #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE ScopedTypeVariables #-}
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
sub :: Num x => (x, x) -> (x, x) -> (x, x)
sub (x1, y1) (x2, y2) = (x1 - x2, y1 - y2)
dotProduct :: Num x => (x, x) -> (x, x) -> x
dotProduct (x1, y1) (x2, y2) = (x1 * x2) + (y1 * y2)
fromIntegralP :: (Integral x, Num y) => (x, x) -> (y, y)
fromIntegralP (x1, y1) = (fromIntegral x1, fromIntegral y1)
mag :: Floating x => (x, x) -> x
mag x = sqrt (dotProduct x x)
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
distMax = (sqrt . fromIntegral $ (heightMax + 1) ^ (2 :: Int) + (widthMax + 1) ^ (2 :: Int)) / 2
accBin = runSTArray $
do arr <- newArray ((0, 0), (thetaSz, distSz)) (0 :: Double)
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'
distance_ = truncate $ distance * fromIntegral distSz / distMax
idx = (theta, distance_)
when (distance_ >= 0 && distance_ < distSz) $
do old <- readArray arr idx
writeArray arr idx (old + 1)
return arr
maxAcc = F.maximum accBin
hTransform (x, y) =
let l = 255 - truncate ((accBin ! (x, y)) / maxAcc * 255)
in PixelY l
hImage :: Image arr Y Word8
hImage = makeImage (thetaSz, distSz) hTransform