Week 3 : The Hough Transform, Updated

Contents

Week 3 : The Hough Transform

Update: The Naive Hough transform, see towards the end, has been updated. I have uploaded a version provided by Cameron Palmer which seems to be correct when compared to MATLAB's implementation. Thank you, Cameron!

This week we are looking at the Hough transform. We will start of by looking at some simple examples on how the Hough Transform works, and continue with some examples on how to use the Hough Transform.

This code can be downloaded from https://github.com/olemarius90/INF4300

clear all
close all

% Start by making a simple image
im=zeros(16);
im(1,1)=1;
im(5,5)=1;
im(8,8)=1;
im(1,16)=1;
im(16,1) = 1;
figure(1)
imagesc(im)
title('A simple image');
xlabel('x');
ylabel('y');

% The question is, are these point on a line?
% That is, is there a line y=ax+b that passes
% through these points for some values of
% a and b?

% Take every x,y coordinate from the given image and plug it into
% the equation b=-xa+y AS PARAMETERS. Then vary a over some
% predefined interval and take a look at the lines produced by this:
% y = ax + b

figure(2);clf
% Point 1, x=1, y=1, b=-a+1

a=-5:5;
b=-a+1;
plot(a,b,'r');
hold on

% Point 2, x=5, y=5, b=-5a+5

a=-5:5;
b=-5*a+5;
plot(a,b,'g');

% Point 3, x=9, y=9, b=-9a+9

a=-5:5;
b=-8*a+8;
plot(a,b,'b');

% Point 4, x = 1, y = 16; b = -a + 16
a=-5:5;
b=-a+16;
plot(a,b,'k');

% Point 5 x = 16, y = 1; b = -16a + 1
a=-5:5;
b=-16*a+1;
plot(a,b,'y');

grid on
axis([-5 5 -50 50])
xlabel('a')
ylabel('b')
title('ab space');

Creating the accumulator matrix

Here we are creating a matrix to represent the "ab-space". If we find the maximum of this matrix we will find the most dominant line in the original image.

% The accumulator matrix
acc=zeros(11,161);

% We need to introduce some offsets since the values for a and b will be
% negative as well.
b_offset = 80;
a_offset = 6;

%For point 1, x=1, y=1, b=-a+1
a=-5:5;
b=-a+1;
for x=1:length(a)
    acc(a(x)+a_offset,b(x)+b_offset)=acc(a(x)+a_offset,b(x)+b_offset)+1;
end

% Point 2, x=5, y=5, b=-5a+5
a=-5:5;
b=-5*a+5;
for x=1:length(a)
    acc(a(x)+a_offset,b(x)+b_offset)=acc(a(x)+a_offset,b(x)+b_offset)+1;
end

% Point 3, x=9, y=9, b=-9a+9
a=-5:5;
b=-9*a+9;
for x=1:length(a)
    acc(a(x)+a_offset,b(x)+b_offset)=acc(a(x)+a_offset,b(x)+b_offset)+1;
end

% Point 4, x = 1, y = 16; b = -a + 16
a=-5:5;
b=-a+16;
for x=1:length(a)
    acc(a(x)+a_offset,b(x)+b_offset)=acc(a(x)+a_offset,b(x)+b_offset)+1;
end

% Point 5 x = 16, y = 1; b = -16a + 1
a=-5:5;
b=-16*a+1;
for x=1:length(a)
    acc(a(x)+a_offset,b(x)+b_offset)=acc(a(x)+a_offset,b(x)+b_offset)+1;
end

figure(3);clf
imagesc((acc))
xlabel('b'),ylabel('a')
title('ba-space');

figure(4);clf
imagesc(rot90(acc))
xlabel('a'),ylabel('b')
title('ab-space (rotated)');

% Lets find the maximum of the matrix
[v,idx] = max(acc(:));
% And we want the index of the maximum
[a_max,b_max] = ind2sub(size(acc),idx)
%But we need to remember to account for the offsets we introduced
a=a_max-a_offset
b=b_max-b_offset

disp(sprintf('We have found the line y=%dx+%d.',a,b))
disp('Does this fit our original image?');

a_max =

     7


b_max =

    80


a =

     1


b =

     0

We have found the line y=1x+0.
Does this fit our original image?

Hough transform example 1

Creating a simple image with one line

img=zeros(11);
for j=3:9
    img(j,j)=1;
end

figure(5)
colormap(gray(2))
imagesc(img);
grid on
axis image
title('Another simple image (image 1)');

% Using the built in function Hough
[H,theta,rho] = hough(img);
figure(6)
imagesc(theta,rho,H);
colormap jet
xlabel('\theta (degrees)')
ylabel('\rho (pixels from center)')
title('Hough transform of image 1')

Hough transform example 2

Creating a simple image with one line

img=zeros(11);
for j=3:9
    img(j,j)=1;
end

% Putting some holes in the line
img(5,5)=0;
img(7,7)=0;

figure(7)
colormap(gray(2))
imagesc(img);
grid on
axis image
title('Another simple image (image 2)');

% Using the built in function Hough
[H,theta,rho] = hough(img);
figure(8)
imagesc(theta,rho,H);
colormap jet
xlabel('\theta (degrees)')
ylabel('\rho (pixels from center)')
title('Hough transform of image 2')

Hough transform example 3

Creating a simple image with two lines

img=zeros(12);
for j=3:9
    img(j,j)=1;
    img(j+3,j) = 1;
end

figure(9)
colormap(gray(2))
imagesc(img);
grid on
axis image
title('Another simple image (image 3)');

% Using the built in function Hough
[H,theta,rho] = hough(img);
figure(10)
imagesc(theta,rho,H);
colormap jet
xlabel('\theta (degrees)')
ylabel('\rho (pixels from center)')
title('Hough transform of image 3')

figure(11)
imagesc(theta,rho,H);
colormap jet
xlim([-60 -20]);
ylim([-10 10]);
xlabel('\theta (degrees)')
ylabel('\rho (pixels from center)')
title('Hough transform of image 3, zoomed in')

Hough transform example 4 : Image of corridors

Now lets look at the image of corridors

img=imread('corridor.png');
img=double(rgb2gray(img));

figure(12)
imshow(img,[])
title('Original image');

pause(1)

% Lets filter the original image with a Sobel filter to find the Sobel
% magnitude
h1=fspecial('sobel');
h2=h1';
igh=imfilter(img,h1);
igv=imfilter(img,h2);
igs=abs(igh)+abs(igv);
figure(13)
imshow(igs,[])
title('Sobel magnitude');

pause(1)

%Lets treshold the Sobel image
%And lets skip the border
igsT=igs(5:end-5,5:end-5)>170;
figure(14)
imshow(igsT)
title('Sobel image tresholded');

pause(1)
%resize the image as well
img = img(5:end-5,5:end-5);

Using the built in function Hough

[H,theta,rho] = hough(igsT);
figure(15);clf
imagesc(theta,rho,H);
colormap hot
colorbar
caxis([0 200 ])
xlabel('\theta (degrees)')
ylabel('\rho (pixels from center)')
title('Hough transform of image 1')

pause(1)

% The rest of this code is taken from a MATLAB example from the "Help" of
% houghlines
% First we find and indicate the 5 most dominant lines
P = houghpeaks(H,5);

figure(16);clf
imshow(H,[],'XData',theta,'YData',rho,'InitialMagnification','fit');
xlabel('\theta'), ylabel('\rho');
axis on, axis normal, hold on;
plot(theta(P(:,2)),rho(P(:,1)),'s','color','white');
title('Houghtransform with 5 most dominant lines indicated');
xlabel('\theta (degrees)')
ylabel('\rho (pixels from center)')
pause(1)
% Find lines and plot them
lines = houghlines(img,theta,rho,P,'FillGap',5,'MinLength',7);
figure(17), imshow(img,[]), hold on
title('Image with lines found indicated');
max_len = 0;
for k = 1:length(lines)
    xy = [lines(k).point1; lines(k).point2];
    plot(xy(:,1),xy(:,2),'LineWidth',2,'Color','green');

    % Plot beginnings and ends of lines
    plot(xy(1,1),xy(1,2),'x','LineWidth',2,'Color','yellow');
    plot(xy(2,1),xy(2,2),'x','LineWidth',2,'Color','red');

    % Determine the endpoints of the longest line segment
    len = norm(lines(k).point1 - lines(k).point2);
    if ( len > max_len)
        max_len = len;
        xy_long = xy;
    end
end

% highlight the longest line segment
plot(xy_long(:,1),xy_long(:,2),'LineWidth',2,'Color','blue');
pause(1)

Detecting circles with built in Huff transform in MATLAB

Get image and display

i=imread('coins2.jpg');
ig=double(rgb2gray(i));

figure(18)
imshow(ig,[])
title('Image of coins');

pause(0.5);

% Make Sobelgradient image and display
h1=fspecial('sobel');
h2=h1';
igh=imfilter(ig,h1);
igv=imfilter(ig,h2);
igs=abs(igh)+abs(igv);

figure(19);
imshow(igs,[]);
title('Sobel gradient magnitude');

pause(1)

%Treshold
igsT=igs>170;

figure(20);
imshow(igsT,[]);
title('Sobel gradient magnitude tresholded');

pause(1)

figure(21);
imshow(ig,[])
[centers, radii, metric] = imfindcircles(igsT,[15 30]);
number_of_circles = 6;
centersStrong5 = centers(1:number_of_circles,:);
radiiStrong5 = radii(1:number_of_circles);
metricStrong5 = metric(1:number_of_circles);
viscircles(centersStrong5, radiiStrong5,'EdgeColor','b');
title('Image with circles indicated');

pause(1)

A naive Huff implementation

Fixed and provided by Cameron Palmer. This seems to be correct when compared to MATLAB's implementation. Thank you Cameron!

%Reading image
img=imread('corridor.png');
img=double(rgb2gray(img));

% Lets filter the original image with a Sobel filter to find the Sobel
% magnitude
h1=fspecial('sobel');
h2=h1';
igh=imfilter(img,h1);
igv=imfilter(img,h2);
igs=abs(igh)+abs(igv);

%Lets treshold the Sobel image
%And lets skip the border
igsT=igs(5:end-5,5:end-5)>170;

I = igsT;

%Using the tresholded Sobel magnitude of the corridor from earlier
[rows, cols] = size(I);

theta_maximum = 90;
rho_maximum = floor(sqrt(rows^2 + cols^2)) - 1;
theta_range = -theta_maximum:theta_maximum - 1;
rho_range = -rho_maximum:rho_maximum;

Hough = zeros(length(rho_range), length(theta_range));

wb = waitbar(0, 'Naive Hough Transform');

for row = 1:rows
    waitbar(row/rows, wb);
    for col = 1:cols
        if I(row, col) > 0
            x = col - 1;
            y = row - 1;
            for theta_ = theta_range
                rho_ = round((x * cosd(theta_)) + (y * sind(theta_)));
                rho_index = rho_ + rho_maximum + 1;
                theta_index = theta_ + theta_maximum + 1;
                Hough(rho_index, theta_index) = Hough(rho_index, theta_index) + 1;
            end
        end
    end
end

close(wb);

figure(22)
imagesc(theta,rho,H);
title('Matlab Hough');

figure(23)
imagesc(theta_range,rho_range,Hough)
title('Naive Hough: Cameron Palmer');

Some work for you

I'll leave it to you to convert this code into the randomised Huff transform. Try to look at the provided PDF and the provided solution and try to make your own implementation. If you get a solution you are happy with feel free to send it to me, and I'll post it here and give credits to you! post_id = 346; %delete this line to force new post; permaLink = http://inf4300.olemarius.net/2015/09/17/omhr_solution-m/;

Published with MATLAB® R2015a