Hough transform
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.
![Task](http://static.miraheze.org/rosettacodewiki/thumb/b/ba/Rcode-button-task-crushed.png/64px-Rcode-button-task-crushed.png)
You are encouraged to solve this task according to the task description, using any language you may know.
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 -space, where the line corresponds to points of the form ). 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.
![](http://static.miraheze.org/rosettacodewiki/thumb/c/c6/Pentagon.png/300px-Pentagon.png)
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.
C
(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>
Tcl
<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 < 720} {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,
or
<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>