Hough transform

From Rosetta Code
Jump to: navigation, search
Task
Hough transform
You are encouraged to solve this task according to the task description, using any language you may know.

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, (ρ,θ), to the average color of the pixels on the corresponding line of the source image (in (x,y)-space, where the line corresponds to points of the form xcosθ + ysinθ = ρ). 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 θ and the other for ρ, 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.

Contents

[edit] BBC BASIC

BBC BASIC uses Cartesian coordinates so the image is 'upside down' compared with some other solutions.

Hough bbc.gif
      Width% = 320
Height% = 240
 
VDU 23,22,Width%;Height%;8,16,16,128
*DISPLAY Pentagon.bmp
OFF
 
DIM hist%(Width%-1, Height%-1)
 
rs = 2 * SQR(Width%^2 + Height%^2) / Height% : REM Radial step
ts = PI / Width% : REM Angular step
h% = Height% / 2
 
REM Hough transform:
FOR y% = 0 TO Height%-1
FOR x% = 0 TO Width%-1
IF TINT(x%*2, y%*2) = 0 THEN
FOR t% = 0 TO Width%-1
th = t% * ts
r% = (x%*COS(th) + y%*SIN(th)) / rs + h% + 0.5
hist%(t%,r%) += 1
NEXT
ENDIF
NEXT
NEXT y%
 
REM Find max:
max% = 0
FOR y% = 0 TO Height%-1
FOR x% = 0 TO Width%-1
IF hist%(x%,y%) > max% max% = hist%(x%,y%)
NEXT
NEXT y%
 
REM Plot:
GCOL 1
FOR y% = 0 TO Height%-1
FOR x% = 0 TO Width%-1
c% = 255 * hist%(x%,y%) / max%
COLOUR 1, c%, c%, c%
LINE x%*2,y%*2,x%*2,y%*2
NEXT
NEXT y%
 
REPEAT
WAIT 1
UNTIL FALSE

[edit] C

[edit] D

Translation of: Go

This uses the module from the Grayscale image Task. The output image is the same as in the Go solution.

import std.math, grayscale_image;
 
Image!Gray houghTransform(in Image!Gray im,
in size_t hx=460, in size_t hy=360)
pure nothrow in {
assert(im !is null);
assert(hx > 0 && hy > 0);
assert((hy & 1) == 0, "hy argument must be even.");
} body {
auto result = new Image!Gray(hx, hy);
result.clear(Gray.white);
 
immutable double rMax = hypot(im.nx, im.ny);
immutable double dr = rMax / (hy / 2.0);
immutable double dTh = PI / hx;
 
foreach (immutable y; 0 .. im.ny) {
foreach (immutable x; 0 .. im.nx) {
if (im[x, y] == Gray.white)
continue;
foreach (immutable iTh; 0 .. hx) {
immutable double th = dTh * iTh;
immutable double r = x * cos(th) + y * sin(th);
immutable iry = hy / 2 - cast(int)floor(r / dr + 0.5);
if (result[iTh, iry] > Gray(0))
result[iTh, iry]--;
}
}
}
return result;
}
 
void main() {
(new Image!RGB)
.loadPPM6("Pentagon.ppm")
.rgb2grayImage()
.houghTransform()
.savePGM("Pentagon_hough.pgm");
}

[edit] Go

Output png
Translation of: Python
package main
 
import (
"fmt"
"image"
"image/color"
"image/draw"
"image/png"
"math"
"os"
)
 
func hough(im image.Image, ntx, mry int) draw.Image {
nimx := im.Bounds().Max.X
mimy := im.Bounds().Max.Y
mry = int(mry/2) * 2
him := image.NewGray(image.Rect(0, 0, ntx, mry))
draw.Draw(him, him.Bounds(), image.NewUniform(color.White),
image.ZP, draw.Src)
 
rmax := math.Hypot(float64(nimx), float64(mimy))
dr := rmax / float64(mry/2)
dth := math.Pi / float64(ntx)
 
for jx := 0; jx < nimx; jx++ {
for iy := 0; iy < mimy; iy++ {
col := color.GrayModel.Convert(im.At(jx, iy)).(color.Gray)
if col.Y == 255 {
continue
}
for jtx := 0; jtx < ntx; jtx++ {
th := dth * float64(jtx)
r := float64(jx)*math.Cos(th) + float64(iy)*math.Sin(th)
iry := mry/2 - int(math.Floor(r/dr+.5))
col = him.At(jtx, iry).(color.Gray)
if col.Y > 0 {
col.Y--
him.SetGray(jtx, iry, col)
}
}
}
}
return him
}
 
func main() {
f, err := os.Open("Pentagon.png")
if err != nil {
fmt.Println(err)
return
}
pent, err := png.Decode(f)
if err != nil {
fmt.Println(err)
return
}
if err = f.Close(); err != nil {
fmt.Println(err)
}
h := hough(pent, 460, 360)
if f, err = os.Create("hough.png"); err != nil {
fmt.Println(err)
return
}
if err = png.Encode(f, h); err != nil {
fmt.Println(err)
}
if cErr := f.Close(); cErr != nil && err == nil {
fmt.Println(err)
}
}

[edit] Haskell

Library: JuicyPixels
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>"

Example use:

HoughTransform Pentagon.png hough.png 360 360

[edit] J

Solution:

NB.*houghTransform v Produces a density plot of image y in hough space
NB. y is picture as an array with 1 at non-white points,
NB. x is resolution (width,height) of resulting image
houghTransform=: dyad define
'w h'=. x NB. width and height of target image
theta=. o. (%~ 0.5+i.) w NB. theta in radians from 0 to π
rho=. (4$.$. |.y) +/ .* 2 1 o./theta NB. rho for each pixel at each theta
'min max'=. (,~-) +/&.:*: $y NB. min/max possible rho
rho=. <. 0.5+ h * (rho-min) % max-min NB. Rescale rho from 0 to h and round to int
|.([: <:@(#/.~) (i.h)&,)"1&.|: rho NB. consolidate into picture
)
Resulting viewmat image from J implementation of Hough Transform on sample pentagon image
Example use:
   require 'viewmat media/platimg'
Img=: readimg jpath '~temp/pentagon.png'
viewmat 460 360 houghTransform _1 > Img


[edit] Java

Code:

import java.awt.image.*;
import java.io.File;
import java.io.IOException;
import javax.imageio.*;
 
public class HoughTransform
{
public static ArrayData houghTransform(ArrayData inputData, int thetaAxisSize, int rAxisSize, int minContrast)
{
int width = inputData.width;
int height = inputData.height;
int maxRadius = (int)Math.ceil(Math.hypot(width, height));
int halfRAxisSize = rAxisSize >>> 1;
ArrayData outputData = new ArrayData(thetaAxisSize, rAxisSize);
// x output ranges from 0 to pi
// y output ranges from -maxRadius to maxRadius
double[] sinTable = new double[thetaAxisSize];
double[] cosTable = new double[thetaAxisSize];
for (int theta = thetaAxisSize - 1; theta >= 0; theta--)
{
double thetaRadians = theta * Math.PI / thetaAxisSize;
sinTable[theta] = Math.sin(thetaRadians);
cosTable[theta] = Math.cos(thetaRadians);
}
 
for (int y = height - 1; y >= 0; y--)
{
for (int x = width - 1; x >= 0; x--)
{
if (inputData.contrast(x, y, minContrast))
{
for (int theta = thetaAxisSize - 1; theta >= 0; theta--)
{
double r = cosTable[theta] * x + sinTable[theta] * y;
int rScaled = (int)Math.round(r * halfRAxisSize / maxRadius) + halfRAxisSize;
outputData.accumulate(theta, rScaled, 1);
}
}
}
}
return outputData;
}
 
public static class ArrayData
{
public final int[] dataArray;
public final int width;
public final int height;
 
public ArrayData(int width, int height)
{
this(new int[width * height], width, height);
}
 
public ArrayData(int[] dataArray, int width, int height)
{
this.dataArray = dataArray;
this.width = width;
this.height = height;
}
 
public int get(int x, int y)
{ return dataArray[y * width + x]; }
 
public void set(int x, int y, int value)
{ dataArray[y * width + x] = value; }
 
public void accumulate(int x, int y, int delta)
{ set(x, y, get(x, y) + delta); }
 
public boolean contrast(int x, int y, int minContrast)
{
int centerValue = get(x, y);
for (int i = 8; i >= 0; i--)
{
if (i == 4)
continue;
int newx = x + (i % 3) - 1;
int newy = y + (i / 3) - 1;
if ((newx < 0) || (newx >= width) || (newy < 0) || (newy >= height))
continue;
if (Math.abs(get(newx, newy) - centerValue) >= minContrast)
return true;
}
return false;
}
 
public int getMax()
{
int max = dataArray[0];
for (int i = width * height - 1; i > 0; i--)
if (dataArray[i] > max)
max = dataArray[i];
return max;
}
}
 
public static ArrayData getArrayDataFromImage(String filename) throws IOException
{
BufferedImage inputImage = ImageIO.read(new File(filename));
int width = inputImage.getWidth();
int height = inputImage.getHeight();
int[] rgbData = inputImage.getRGB(0, 0, width, height, null, 0, width);
ArrayData arrayData = new ArrayData(width, height);
// Flip y axis when reading image
for (int y = 0; y < height; y++)
{
for (int x = 0; x < width; x++)
{
int rgbValue = rgbData[y * width + x];
rgbValue = (int)(((rgbValue & 0xFF0000) >>> 16) * 0.30 + ((rgbValue & 0xFF00) >>> 8) * 0.59 + (rgbValue & 0xFF) * 0.11);
arrayData.set(x, height - 1 - y, rgbValue);
}
}
return arrayData;
}
 
public static void writeOutputImage(String filename, ArrayData arrayData) throws IOException
{
int max = arrayData.getMax();
BufferedImage outputImage = new BufferedImage(arrayData.width, arrayData.height, BufferedImage.TYPE_INT_ARGB);
for (int y = 0; y < arrayData.height; y++)
{
for (int x = 0; x < arrayData.width; x++)
{
int n = Math.min((int)Math.round(arrayData.get(x, y) * 255.0 / max), 255);
outputImage.setRGB(x, arrayData.height - 1 - y, (n << 16) | (n << 8) | 0x90 | -0x01000000);
}
}
ImageIO.write(outputImage, "PNG", new File(filename));
return;
}
 
public static void main(String[] args) throws IOException
{
ArrayData inputData = getArrayDataFromImage(args[0]);
int minContrast = (args.length >= 4) ? 64 : Integer.parseInt(args[4]);
ArrayData outputData = houghTransform(inputData, Integer.parseInt(args[2]), Integer.parseInt(args[3]), minContrast);
writeOutputImage(args[1], outputData);
return;
}
}
Output from example pentagon image
Example use:
java HoughTransform pentagon.png JavaHoughTransform.png 640 480 100


[edit] Mathematica

 
Radon[image, Method -> "Hough"]
 

[edit] MATLAB

[edit] Python

Library: PIL

This is the classical Hough transform as described in wikipedia. The code does not compute averages; it merely makes a point on the transformed image darker if a lot of points on the original image lie on the corresponding line. The output is almost identical to that of the Tcl code. The code works only with gray-scale images, but it is easy to extend to RGB.

 
from math import hypot, pi, cos, sin
import Image
 
 
def hough(im, ntx=460, mry=360):
"Calculate Hough transform."
pim = im.load()
nimx, mimy = im.size
mry = int(mry/2)*2 #Make sure that this is even
him = Image.new("L", (ntx, mry), 255)
phim = him.load()
 
rmax = hypot(nimx, mimy)
dr = rmax / (mry/2)
dth = pi / ntx
 
for jx in xrange(nimx):
for iy in xrange(mimy):
col = pim[jx, iy]
if col == 255: continue
for jtx in xrange(ntx):
th = dth * jtx
r = jx*cos(th) + iy*sin(th)
iry = mry/2 + int(r/dr+0.5)
phim[jtx, iry] -= 1
return him
 
 
def test():
"Test Hough transform with pentagon."
im = Image.open("pentagon.png").convert("L")
him = hough(im)
him.save("ho5.bmp")
 
 
if __name__ == "__main__": test()
 
 

[edit] Racket

[edit] Ruby

 
require 'mathn'
require 'rubygems'
require 'gd2'
include GD2
 
def hough_transform(img)
mx, my = img.w*0.5, img.h*0.5
max_d = Math.sqrt(mx**2 + my**2)
min_d = max_d * -1
hough = Hash.new(0)
(0..img.w).each do |x|
puts "#{x} of #{img.w}"
(0..img.h).each do |y|
if img.pixel2color(img.get_pixel(x,y)).g > 32
(0...180).each do |a|
rad = a * (Math::PI / 180.0)
d = (x-mx) * Math.cos(rad) + (y-my) * Math.sin(rad)
hough["#{a.to_i}_#{d.to_i}"] = hough["#{a.to_i}_#{d.to_i}"] + 1
end
end
end
end
heat = GD2::Image.import 'heatmap.png'
out = GD2::Image::TrueColor.new(180,max_d*2)
max = hough.values.max
p max
hough.each_pair do |k,v|
a,d = k.split('_').map(&:to_i)
c = (v / max) * 255
c = heat.get_pixel(c,0)
out.set_pixel(a, max_d + d, c)
end
out
end
 

[edit] Tcl

Library: Tk
Personal tools
Namespaces

Variants
Actions
Community
Explore
Misc
Toolbox