Contents
- Week 3 : The Hough Transform
- Creating the accumulator matrix
- Hough transform example 1
- Hough transform example 2
- Hough transform example 3
- Hough transform example 4 : Image of corridors
- Using the built in function Hough
- Detecting circles with built in Huff transform in MATLAB
- A naive Huff implementation
- Some work for you
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/;
Out of curiosity - has anyone gotten equivalent results to Matlab on the circles task using the randomized Hough transform?