Task

Implement the Hough transform, which is used as part of feature extraction with digital images.

It is a tool that makes it far easier to identify straight lines in the source image, whatever their orientation.

The transform maps each point in the target image, 



(
ρ
,
θ
)


{\displaystyle (\rho ,\theta )}

, to the average color of the pixels on the corresponding line of the source image (in 



(
x
,
y
)


{\displaystyle (x,y)}

-space, where the line corresponds to points of the form 



x
cos
⁡
θ
+
y
sin
⁡
θ
=
ρ


{\displaystyle x\cos \theta +y\sin \theta =\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.

 Sample PNG image to use for the Hough transform.
The target space actually uses polar coordinates, but is conventionally plotted on rectangular coordinates for display. There's no specification of exactly how to map polar coordinates to a flat surface for display, but a convenient method is to use one axis for 



θ


{\displaystyle \theta }

 and the other for 



ρ


{\displaystyle \rho }

, with the center of the source image being the origin.

There is also a spherical Hough transform, which is more suited to identifying planes in 3D data.




import Control.Monad (forM_, when)
import Data.Array ((!))
import Data.Array.ST (newArray, writeArray, readArray, runSTArray)
import qualified Data.Foldable as F (maximum)
import System.Environment (getArgs, getProgName)
 
-- Library JuicyPixels:
import Codec.Picture
       (DynamicImage(ImageRGB8, ImageRGBA8), Image, PixelRGB8(PixelRGB8),
        PixelRGBA8(PixelRGBA8), imageWidth, imageHeight, pixelAt,
        generateImage, readImage, pixelMap, savePngImage)
import Codec.Picture.Types (extractLumaPlane, dropTransparency)
 
dot
  :: Num a
  => (a, a) -> (a, a) -> a
dot (x1, y1) (x2, y2) = x1 * x2 + y1 * y2
 
mag
  :: Floating a
  => (a, a) -> a
mag a = sqrt $ dot a a
 
sub
  :: Num a
  => (a, a) -> (a, a) -> (a, a)
sub (x1, y1) (x2, y2) = (x1 - x2, y1 - y2)
 
fromIntegralP
  :: (Integral a, Num b)
  => (a, a) -> (b, b)
fromIntegralP (x, y) = (fromIntegral x, fromIntegral y)
 
{-
  Create a Hough space image with y+ measuring the distance from
  the center of the input image on the range of 0 to half the hypotenuse
  and x+ measuring from [0, 2 * pi].
  The origin is in the upper left, so y is increasing down.
  The image is scaled according to thetaSize and distSize.
-}
hough :: Image PixelRGB8 -> Int -> Int -> Image PixelRGB8
hough image thetaSize distSize = hImage
  where
    width = imageWidth image
    height = imageHeight image
    wMax = width - 1
    hMax = height - 1
    xCenter = wMax `div` 2
    yCenter = hMax `div` 2
    lumaMap = extractLumaPlane image
    gradient x y =
      let orig = pixelAt lumaMap x y
          x_ = pixelAt lumaMap (min (x + 1) wMax) y
          y_ = pixelAt lumaMap x (min (y + 1) hMax)
      in fromIntegralP (orig - x_, orig - y_)
    gradMap =
      [ ((x, y), gradient x y)
      | x <- [0 .. wMax] 
      , y <- [0 .. hMax] ]
    -- The longest distance from the center, half the hypotenuse of the image.
    distMax :: Double
    distMax = (sqrt . fromIntegral $ height ^ 2 + width ^ 2) / 2
    {-
      The accumulation bins of the polar values.
      For each value in the gradient image, if the gradient length exceeds
      some threshold, consider it evidence of a line and plot all of the
      lines that go through that point in Hough space.
    -}
    accBin =
      runSTArray $
      do arr <- newArray ((0, 0), (thetaSize, distSize)) 0
         forM_ gradMap $
           \((x, y), grad) -> do
             let (x_, y_) = fromIntegralP $ (xCenter, yCenter) `sub` (x, y)
             when (mag grad > 127) $
               forM_ [0 .. thetaSize] $
               \theta -> do
                 let theta_ =
                       fromIntegral theta * 360 / fromIntegral thetaSize / 180 *
                       pi :: Double
                     dist = cos theta_ * x_ + sin theta_ * y_
                     dist_ = truncate $ dist * fromIntegral distSize / distMax
                     idx = (theta, dist_)
                 when (dist_ >= 0 && dist_ < distSize) $
                   do old <- readArray arr idx
                      writeArray arr idx $ old + 1
         return arr
    maxAcc = F.maximum accBin
    -- The image representation of the accumulation bins.
    hTransform x y =
      let l = 255 - truncate ((accBin ! (x, y)) / maxAcc * 255)
      in PixelRGB8 l l l
    hImage = generateImage hTransform thetaSize distSize
 
houghIO :: FilePath -> FilePath -> Int -> Int -> IO ()
houghIO path outpath thetaSize distSize = do
  image <- readImage path
  case image of
    Left err -> putStrLn err
    Right (ImageRGB8 image_) -> doImage image_
    Right (ImageRGBA8 image_) -> doImage $ pixelMap dropTransparency image_
    _ -> putStrLn "Expecting RGB8 or RGBA8 image"
  where
    doImage image = do
      let houghImage = hough image thetaSize distSize
      savePngImage outpath $ ImageRGB8 houghImage
 
main :: IO ()
main = do
  args <- getArgs
  prog <- getProgName
  case args of
    [path, outpath, thetaSize, distSize] ->
      houghIO path outpath (read thetaSize) (read distSize)
    _ ->
      putStrLn $
      "Usage: " ++ prog ++ " <image-file> <out-file.png> <width> <height>"