Hough transform: Difference between revisions
Pointed to sub page |
Transfer C example to Example:Hough transform/C |
||
Line 62: | Line 62: | ||
=={{header|C}}== |
=={{header|C}}== |
||
[[file:penta-hugh.png|thumb]][[file:hugh-lines-in.png|thumb]][[file:hugh-lines-out.png|thumb]] |
|||
This code is a little to long to my liking, because I had to put some ad hoc stuff that should be better served by libraries. But you don't want to see libpng code here, trust me. |
|||
<lang C>#include <stdio.h> |
|||
#include <stdlib.h> |
|||
#include <unistd.h> |
|||
#include <string.h> |
|||
#include <ctype.h> |
|||
#include <sys/types.h> |
|||
#include <sys/stat.h> |
|||
#include <fcntl.h> |
|||
#include <err.h> |
|||
#include <math.h> |
|||
* see [[Example:Hough transform/C]] |
|||
/* start of utility functions: not interesting */ |
|||
typedef unsigned char uchar; |
|||
typedef unsigned long ulong; |
|||
typedef struct intensity_t { |
|||
double **pix; |
|||
long width, height; |
|||
} *intensity; |
|||
double PI; |
|||
#define decl_array_alloc(type) \ |
|||
type ** type##_array(long w, long h) { \ |
|||
int i; \ |
|||
type ** row = malloc(sizeof(type*) * h); \ |
|||
type * pix = malloc(sizeof(type) * h * w); \ |
|||
for (i = 0; i < h; i++) \ |
|||
row[i] = pix + w * i; \ |
|||
memset(pix, 0, sizeof(type) * h * w); \ |
|||
return row; \ |
|||
} |
|||
decl_array_alloc(double); |
|||
decl_array_alloc(ulong); |
|||
intensity intensity_alloc(long w, long h) |
|||
{ |
|||
intensity x = malloc(sizeof(struct intensity_t)); |
|||
x->width = w; |
|||
x->height = h; |
|||
x->pix = double_array(w, h); |
|||
return x; |
|||
} |
|||
long get_num(uchar **p, uchar *buf_end) |
|||
{ |
|||
uchar *ptr = *p, *tok_end; |
|||
long tok; |
|||
while (1) { |
|||
while (ptr < buf_end && isspace(*ptr)) ptr++; |
|||
if (ptr >= buf_end) return 0; |
|||
if (*ptr == '#') { /* ignore comment */ |
|||
while (ptr++ < buf_end) { |
|||
if (*ptr == '\n' || *ptr == '\r') break; |
|||
} |
|||
continue; |
|||
} |
|||
tok = strtol((char*)ptr, (char**)&tok_end, 10); |
|||
if (tok_end == ptr) return 0; |
|||
*p = tok_end; |
|||
return tok; |
|||
} |
|||
return 0; |
|||
} |
|||
/* Note: not robust. A robust version would be to long for example code */ |
|||
intensity read_pnm(char *name) |
|||
{ |
|||
struct stat st; |
|||
uchar *fbuf, *ptr, *end; |
|||
long width, height, max_val; |
|||
int i, j; |
|||
intensity ret; |
|||
int fd = open(name, O_RDONLY); |
|||
if (fd == -1) err(1, "Can't open %s", name); |
|||
/* from now on assume all operations succeed */ |
|||
fstat(fd, &st); |
|||
fbuf = malloc(st.st_size + 1); |
|||
read(fd, fbuf, st.st_size); |
|||
*(end = fbuf + st.st_size) = '\0'; |
|||
close(fd); |
|||
if (fbuf[0] != 'P' || (fbuf[1] != '5' && fbuf[1] != '6') || !isspace(fbuf[2])) |
|||
err(1, "%s: bad format: can only do P5 or P6 pnm", name); |
|||
ptr = fbuf + 3; |
|||
width = get_num(&ptr, end); |
|||
height = get_num(&ptr, end); |
|||
max_val = get_num(&ptr, end); |
|||
if (max_val <= 0 || max_val >= 256) |
|||
err(1, "Can't handle pixel value %ld\n", max_val); |
|||
fprintf(stderr, "[Info] format: P%c w: %ld h: %ld value: %ld\n", |
|||
fbuf[1], width, height, max_val); |
|||
ret = intensity_alloc(width, height); |
|||
ptr ++; /* ptr should be pointint at the first pixel byte now */ |
|||
if (fbuf[1] == '5') { /* graymap, 1 byte per pixel */ |
|||
for (i = 0; i < height; i++) { |
|||
for (j = 0; j < width; j++) { |
|||
ret->pix[i][j] = (double)*(ptr++) / max_val; |
|||
} |
|||
} |
|||
} else { /* pnm, 1 byte each for RGB */ |
|||
/* hocus pocus way of getting lightness from RGB for us */ |
|||
for (i = 0; i < height; i++) { |
|||
for (j = 0; j < width; j++) { |
|||
ret->pix[i][j] = (ptr[0] * 0.2126 + |
|||
ptr[1] * 0.7152 + |
|||
ptr[2] * 0.0722) / max_val; |
|||
ptr += 3; |
|||
} |
|||
} |
|||
} |
|||
free(fbuf); |
|||
return ret; |
|||
} |
|||
void write_pgm(double **pix, long w, long h) |
|||
{ |
|||
long i, j; |
|||
unsigned char *ptr, *buf = malloc(sizeof(double) * w * h); |
|||
char header[1024]; |
|||
sprintf(header, "P5\n%ld %ld\n255\n", w, h); |
|||
ptr = buf; |
|||
for (i = 0; i < h; i++) |
|||
for (j = 0; j < w; j++) |
|||
*(ptr++) = 256 * pix[i][j]; |
|||
write(fileno(stdout), header, strlen(header)); |
|||
write(fileno(stdout), buf, w * h); |
|||
free(buf); |
|||
} |
|||
/* Finally, end of util functions. All that for this function. */ |
|||
intensity hugh_transform(intensity in, double gamma) |
|||
{ |
|||
long i, j, k, l, m, w, h; |
|||
double bg, r_res, t_res, rho, r, theta, x, y, v, max_val, min_val, *pp; |
|||
intensity graph; |
|||
/* before anything else, legalize Pi = 3 */ |
|||
PI = atan2(1, 1) * 4; |
|||
/* first, run through all pixels and see what the average is, |
|||
* so we can take a guess if the background is black or white. |
|||
* a real application wouldn't do silly things like this */ |
|||
for (i = 0, bg = 0; i < in->height; i++) |
|||
for (j = 0; j < in->width; j++) |
|||
bg += in->pix[i][j]; |
|||
fprintf(stderr, "[info] background is %f\n", bg); |
|||
bg = (bg /= (in->height * in->width) > 0.5) ? 1 : 0; |
|||
/* if white, invert it */ |
|||
if (bg) { |
|||
for (i = 0; i < in->height; i++) |
|||
for (j = 0; j < in->width; j++) |
|||
in->pix[i][j] = 1 - in->pix[i][j]; |
|||
} |
|||
/* second, decide what resolution of rho and theta should be. |
|||
* here we just make the rho/theta graph a fixed ratio |
|||
* of input, which is dumb. It should depend on the application. |
|||
* finer bins allow better resolution between lines, but will |
|||
* lose contrast if the input is noisy. Also, lower resolution, faster. |
|||
*/ |
|||
# define RRATIO 1.5 |
|||
# define TRATIO 1.5 |
|||
x = in->width - .5; |
|||
y = in->height - .5; |
|||
r = sqrt(x * x + y * y) / 2; |
|||
w = in->width / TRATIO; |
|||
h = in->height / RRATIO; |
|||
r_res = r / h; |
|||
t_res = PI * 2 / w; |
|||
graph = intensity_alloc(w, h); |
|||
for (i = 0; i < in->height; i++) { |
|||
y = i - in->height / 2. + .5; |
|||
for (j = 0; j < in->width; j++) { |
|||
x = j - in->width / 2 + .5; |
|||
r = sqrt(x * x + y * y); |
|||
v = in->pix[i][j]; |
|||
/* hackery: sample image is mostly blank, this saves a great |
|||
* deal of time. Doesn't help a lot with noisy images */ |
|||
if (!v) continue; |
|||
/* at each pixel, check what lines it could be on */ |
|||
for (k = 0; k < w; k++) { |
|||
theta = k * t_res - PI; |
|||
rho = x * cos(theta) + y * sin(theta); |
|||
if (rho >= 0) { |
|||
m = rho / r_res; |
|||
l = k; |
|||
} else { |
|||
m = -rho / r_res; |
|||
l = (k + w/2.); |
|||
l %= w; |
|||
} |
|||
graph->pix[m][l] += v * r; |
|||
} |
|||
} |
|||
/* show which row we are precessing lest user gets bored */ |
|||
fprintf(stderr, "\r%ld", i); |
|||
} |
|||
fprintf(stderr, "\n"); |
|||
max_val = 0; |
|||
min_val = 1e100; |
|||
pp = &(graph->pix[graph->height - 1][graph->width - 1]); |
|||
for (i = graph->height * graph->width - 1; i >= 0; i--, pp--) { |
|||
if (max_val < *pp) max_val = *pp; |
|||
if (min_val > *pp) min_val = *pp; |
|||
} |
|||
/* gamma correction. if gamma > 1, output contrast is better, noise |
|||
is suppressed, but spots for thin lines may be lost; if gamma < 1, |
|||
everything is brighter, both lines and noises */ |
|||
pp = &(graph->pix[graph->height - 1][graph->width - 1]); |
|||
for (i = graph->height * graph->width - 1; i >= 0; i--, pp--) { |
|||
*pp = pow((*pp - min_val)/ (max_val - min_val), gamma); |
|||
} |
|||
return graph; |
|||
} |
|||
int main() |
|||
{ |
|||
//intensity in = read_pnm("pent.pnm"); |
|||
intensity in = read_pnm("lines.pnm"); |
|||
intensity out = hugh_transform(in, 1.5); |
|||
/* binary output goes straight to stdout, get ready to see garbage on your |
|||
* screen if you are not careful! |
|||
*/ |
|||
write_pgm(out->pix, out->width, out->height); |
|||
/* not going to free memory we used: OS can deal with it */ |
|||
return 0; |
|||
}</lang> |
|||
This program takes a pnm file (binary, either P5 or P6) and does the transformation, then dump output onto stdout. Sample images below are output from the pentagram; sample lines with added noise; output of processing that. Both output were with 1.5 gamma. |
|||
=={{header|D}}== |
=={{header|D}}== |
Revision as of 15:29, 30 March 2013
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 -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.
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.
BBC BASIC
BBC BASIC uses Cartesian coordinates so the image is 'upside down' compared with some other solutions.
<lang bbcbasic> 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</lang>
C
D
This uses the module from the Grayscale image Task. The output image is the same as in the Go solution. <lang d>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.max) // Not pure. if (im[x, y] == Gray(255)) 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() {
auto im = new Image!RGB; loadPPM6(im, "Pentagon.ppm") .rgb2grayImage() .houghTransform() .savePGM("Pentagon_hough.pgm");
}</lang>
Go
<lang go>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) }
}</lang>
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>
Example use:
<lang j> require 'viewmat media/platimg'
Img=: readimg jpath '~temp/pentagon.png' viewmat 460 360 houghTransform _1 > Img</lang>
Java
Code: <lang Java>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; }
}</lang>
Example use:
java HoughTransform pentagon.png JavaHoughTransform.png 640 480 100
Mathematica
<lang Mathematica> Radon[image, Method -> "Hough"] </lang>
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>
Python
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. <lang python> 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()
</lang>
Ruby
<lang 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 </lang>