Rosetta Code

Sandbox/Realazthat/Hough transform: Difference between revisions

Page Discussion

Sandbox/Realazthat/Hough transform (view source)

Revision as of 17:42, 22 October 2010

9,961 bytes removed ,  13 years ago
no edit summary
(Created page with '{{task|Image processing}}Category:Graphics algorithms Implement the Hough transform, which is used as part of feature extraction with digital images. I…')
 
No edit summary
 
(6 intermediate revisions by the same user not shown)
Line 1:
{{task/realazthat
{{task|Image processing}}[[Category:Graphics algorithms]]
|name=Hough transform
|type=
|description=
 
Implement the [[wp:Hough transform|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.
 
Line 8 ⟶ 12:
 
There is also a spherical Hough transform, which is more suited to identifying planes in 3D data.
|algorithms=
 
=={{header|C}}==
{{trans|Tcl}}
 
{{libheader|cairo}}
 
(Tested only with the pentagon image given)
<lang c>#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <string.h>
#include <math.h>
 
#include <cairo.h>
 
#ifndef M_PI
#define M_PI 3.1415927
#endif
 
#define GR(X,Y) (d[(*s)*(Y)+bpp*(X)+((2)%bpp)])
#define GG(X,Y) (d[(*s)*(Y)+bpp*(X)+((1)%bpp)])
#define GB(X,Y) (d[(*s)*(Y)+bpp*(X)+((0)%bpp)])
#define SR(X,Y) (ht[4*tw*((Y)%th)+4*((X)%tw)+2])
#define SG(X,Y) (ht[4*tw*((Y)%th)+4*((X)%tw)+1])
#define SB(X,Y) (ht[4*tw*((Y)%th)+4*((X)%tw)+0])
#define RAD(A) (M_PI*((double)(A))/180.0)
uint8_t *houghtransform(uint8_t *d, int *w, int *h, int *s, int bpp)
{
int rho, theta, y, x, W = *w, H = *h;
int th = sqrt(W*W + H*H)/2.0;
int tw = 360;
uint8_t *ht = malloc(th*tw*4);
memset(ht, 0, 4*th*tw); // black bg
for(rho = 0; rho < th; rho++)
{
for(theta = 0; theta < tw/*720*/; theta++)
{
double C = cos(RAD(theta));
double S = sin(RAD(theta));
uint32_t totalred = 0;
uint32_t totalgreen = 0;
uint32_t totalblue = 0;
uint32_t totalpix = 0;
if ( theta < 45 || (theta > 135 && theta < 225) || theta > 315) {
for(y = 0; y < H; y++) {
double dx = W/2.0 + (rho - (H/2.0-y)*S)/C;
if ( dx < 0 || dx >= W ) continue;
x = floor(dx+.5);
if (x == W) continue;
totalpix++;
totalred += GR(x, y);
totalgreen += GG(x, y);
totalblue += GB(x, y);
}
} else {
for(x = 0; x < W; x++) {
double dy = H/2.0 - (rho - (x - W/2.0)*C)/S;
if ( dy < 0 || dy >= H ) continue;
y = floor(dy+.5);
if (y == H) continue;
totalpix++;
totalred += GR(x, y);
totalgreen += GG(x, y);
totalblue += GB(x, y);
}
}
if ( totalpix > 0 ) {
double dp = totalpix;
SR(theta, rho) = (int)(totalred/dp) &0xff;
SG(theta, rho) = (int)(totalgreen/dp) &0xff;
SB(theta, rho) = (int)(totalblue/dp) &0xff;
}
}
}
 
*h = th; // sqrt(W*W+H*H)/2
*w = tw; // 360
*s = 4*tw;
return ht;
}
 
int main(int argc, char **argv)
{
cairo_surface_t *inputimg = NULL;
cairo_surface_t *houghimg = NULL;
uint8_t *houghdata = NULL, *inputdata = NULL;
int w, h, s, bpp;
 
if ( argc < 3 ) return EXIT_FAILURE;
inputimg = cairo_image_surface_create_from_png(argv[1]);
w = cairo_image_surface_get_width(inputimg);
h = cairo_image_surface_get_height(inputimg);
s = cairo_image_surface_get_stride(inputimg);
bpp = cairo_image_surface_get_format(inputimg);
switch(bpp)
{
case CAIRO_FORMAT_ARGB32: bpp = 4; break;
case CAIRO_FORMAT_RGB24: bpp = 3; break;
case CAIRO_FORMAT_A8: bpp = 1; break;
default:
fprintf(stderr, "unsupported\n");
goto destroy;
}
 
inputdata = cairo_image_surface_get_data(inputimg);
houghdata = houghtransform(inputdata, &w, &h, &s, bpp);
printf("w=%d, h=%d\n", w, h);
houghimg = cairo_image_surface_create_for_data(houghdata,
CAIRO_FORMAT_RGB24,
w, h, s);
cairo_surface_write_to_png(houghimg, argv[2]);
destroy:
if (inputimg != NULL) cairo_surface_destroy(inputimg);
if (houghimg != NULL) cairo_surface_destroy(houghimg);
 
return EXIT_SUCCESS;
}</lang>
 
Output image (but with white background):
 
[[Image:Houghtrasf-c.png|thumb|left|360x200px|Output image when input is the given Pentagon.]]
<br style="clear:both" />
 
=={{header|J}}==
'''Solution:'''
<lang j>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
)</lang>
[[Image:JHoughTransform.png|320px200px|thumb|right|Resulting viewmat image from J implementation of Hough Transform on sample pentagon image]]'''Example use:'''
<lang j> require 'viewmat media/platimg'
Img=: readimg jpath '~temp/pentagon.png'
viewmat 460 360 houghTransform _1 > Img</lang>
<br style="clear:both" />
 
=={{header|MATLAB}}==
This solution takes an image and the theta resolution as inputs. The image itself must be a 2-D boolean array. This array is constructed such that all of the pixels on an edge have the value "true." This can be done for a normal image using an "edge finding" algorithm to preprocess the image. In the case of the example image the pentagon "edges" are black pixels. So when the image is imported into MATLAB simply say any pixel colored black is true. The syntax is usually, cdata < 255. Where the vale 255 represents white and 0 represents black.
 
<lang MATLAB>function [rho,theta,houghSpace] = houghTransform(theImage,thetaSampleFrequency)
 
%Define the hough space
theImage = flipud(theImage);
[width,height] = size(theImage);
rhoLimit = norm([width height]);
rho = (-rhoLimit:1:rhoLimit);
theta = (0:thetaSampleFrequency:pi);
numThetas = numel(theta);
houghSpace = zeros(numel(rho),numThetas);
%Find the "edge" pixels
[xIndicies,yIndicies] = find(theImage);
%Preallocate space for the accumulator array
numEdgePixels = numel(xIndicies);
accumulator = zeros(numEdgePixels,numThetas);
%Preallocate cosine and sine calculations to increase speed. In
%addition to precallculating sine and cosine we are also multiplying
%them by the proper pixel weights such that the rows will be indexed by
%the pixel number and the columns will be indexed by the thetas.
%Example: cosine(3,:) is 2*cosine(0 to pi)
% cosine(:,1) is (0 to width of image)*cosine(0)
cosine = (0:width-1)'*cos(theta); %Matrix Outerproduct
sine = (0:height-1)'*sin(theta); %Matrix Outerproduct
accumulator((1:numEdgePixels),:) = cosine(xIndicies,:) + sine(yIndicies,:);
 
%Scan over the thetas and bin the rhos
for i = (1:numThetas)
houghSpace(:,i) = hist(accumulator(:,i),rho);
end
 
pcolor(theta,rho,houghSpace);
shading flat;
title('Hough Transform');
xlabel('Theta (radians)');
ylabel('Rho (pixels)');
colormap('gray');
 
end</lang>
 
Sample Usage:
<lang MATLAB>>> uiopen('C:\Documents and Settings\owner\Desktop\Chris\MATLAB\RosettaCode\180px-Pentagon.png',1)
>> houghTransform(cdata(:,:,1)<255,1/200); %The image from uiopen is stored in cdata. The reason why the image is cdata<255 is because the "edge" pixels are black.</lang>
[[Image:HoughTransformHex.png|thumb|left|360x200px|Image produced by MATLAB implementation of the Hough transform when applied to the sample pentagon image.]]
<br style="clear:both" />
 
=={{header|Tcl}}==
{{libheader|Tk}}
 
<lang tcl>package require Tk
 
set PI 3.1415927
proc HoughTransform {src trg {fieldColor "#000000"}} {
global PI
 
set w [image width $src]
set h [image height $src]
set targetH [expr {int(hypot($w, $h)/2)}]
 
# Configure the target buffer
$trg configure -width 360 -height $targetH
$trg put $fieldColor -to 0 0 359 [expr {$targetH-1}]
 
# Iterate over the target's space of pixels.
for {set rho 0} {$rho < $targetH} {incr rho} {
set row {}
for {set theta 0} {$theta < 360} {incr theta} {
set cos [expr {cos($theta/180.0*$PI)}]
set sin [expr {sin($theta/180.0*$PI)}]
set totalRed 0
set totalGreen 0
set totalBlue 0
set totalPix 0
 
# Sum the colors of the line with equation x*cos(θ) + y*sin(θ) = ρ
if {$theta<45 || ($theta>135 && $theta<225) || $theta>315} {
# For these half-quadrants, it's better to iterate by 'y'
for {set y 0} {$y<$h} {incr y} {
set x [expr {
$w/2 + ($rho - ($h/2-$y)*$sin)/$cos
}]
if {$x < 0 || $x >= $w} continue
set x [expr {round($x)}]
if {$x == $w} continue
incr totalPix
lassign [$src get $x $y] r g b
incr totalRed $r
incr totalGreen $g
incr totalBlue $b
}
} else {
# For the other half-quadrants, it's better to iterate by 'x'
for {set x 0} {$x<$w} {incr x} {
set y [expr {
$h/2 - ($rho - ($x-$w/2)*$cos)/$sin
}]
if {$y < 0 || $y >= $h} continue
set y [expr {round($y)}]
if {$y == $h} continue
incr totalPix
lassign [$src get $x $y] r g b
incr totalRed $r
incr totalGreen $g
incr totalBlue $b
}
}
 
# Convert the summed colors back into a pixel for insertion into
# the target buffer.
if {$totalPix > 0} {
set totalPix [expr {double($totalPix)}]
set col [format "#%02x%02x%02x" \
[expr {round($totalRed/$totalPix)}] \
[expr {round($totalGreen/$totalPix)}] \
[expr {round($totalBlue/$totalPix)}]]
} else {
set col $fieldColor
}
lappend row $col
}
$trg put [list $row] -to 0 $rho
}
}</lang>
 
===Demonstration Code===
Takes the name of the image to apply the transform to as an argument. If using PNG images,
 
{{works with|Tk|8.6}} or {{libheader|TkImg}}
<lang tcl># Demonstration code
if {[catch {
package require Tk 8.6; # Just for PNG format handler
}] == 1} then {catch {
package require Img
}}
# If neither Tk8.6 nor Img, then only GIF and PPM images can be loaded
 
set f [lindex $argv 0]
image create photo srcImg -file $f
image create photo targetImg
pack [labelframe .l1 -text Source] [labelframe .l2 -text Target]
pack [label .l1.i -image srcImg]
pack [label .l2.i -image targetImg]
# Postpone until after we've drawn ourselves
after idle HoughTransform srcImg targetImg [lrange $argv 1 end]</lang>
 
[[Image:Hough-Pentagon-Tcl-Results.gif|thumb|left|360x200px|Image produced by Tcl implementation of the Hough transform when applied to the sample pentagon image.]]
Retrieved from "https://rosettacode.org/wiki/Special:MobileDiff/160363...160370"