21 April, 2014

Understanding the Hough transform #2

In a previous post the basics of the Hough transform were explained. Although we have the built-in hough function in MATLAB to do this operation, it is definitely worth to write our optimized version for study purposes. In addition this lets us to have a function meeting our needs better if needed.

First have a look at the next figure illustrating the θ, r space equation of lines:

We want to get the corresponding θ, r pairs of all lines going through point P at x0, y0. In the previous post it was shown, that the equation of the line can be rearranged to:

$$ r_\theta = x_0\cdot cos(\theta)+y_0\cdot sin(\theta) $$

To get all lines going through P we have to loop through all θ values from 0 to 179 degrees, and calculate the appropriate r value for each θ. It is worth to mention, that r can be negative, too.

In the application in the previous post negative r occurs, when the direction of the normal vector of the line and the vector pointing from the origin to the line has opposite direction.

The input of the Hough transform is usually a binary image containing the edge pixels. During the transform we iterate through the edge points and calculate all θ, r pairs for them: a sinusoid curve will belong to each point. The crossing of two sinusoids means, that the two corresponding points are on the same line. The more sinusoids cross a given point, the more edge points are on the same line. Our task is to find the crossings. How shall we do it? There are infinite number of possible θ, r pairs...

The answer is quantization: we calculate the r values at given discrete θ angles and then quantize r, too. This way we reduce the number of the possible outcomes. After processing all the edge points, there will be lots of quantized θ, r pairs. The only remaining step is to create a two-dimensional histogram from them, and find the peaks to identify the lines on the image.

To understand it better, have a look at the following picture:

Two sinusoids are visualized, a red and a blue one, while the small rectangles are the bins of the two-dimensional histogram. The empty bins have zero value, while the blue and red ones contain one. The intersection of the sinusoids are marked by green having value two.

Turning this to a MATLAB function

The first task is to calculate the quantized angles:

% quantized angles with step dtheta
theta = 0 : dtheta : pi;

% number of angles
N = length(theta);

Next we have to get the coordinates of all the edge points:

% x, y are column vectors containing the coordinates of
% the non-zero edge points
[y, x] = find(image);

% the number of edge points on the image
P = length(x);

Before the next step, have a look at the following equation again:

$$ r_\theta = x_0\cdot cos(\theta)+y_0\cdot sin(\theta) $$

This is the way to calculate the r value for a given x0, y0 point at angle θ. Remember, we have lots of points, so optimization is very important.

The first thing we realize, that it is unnecessary to calculate the sine and cosine values again and again, since they are the same for all points. Another important consideration is, that we can get all the r values by a single matrix multiplication.

Suppose, that matrix P contains the coordinates of the edge points, while in matrix A we have an angle step of 30 degrees:

$$ P = \begin{bmatrix} x_0 & y_0 \\ x_1 & y_1 \\ \vdots & \vdots \\ x_p & y_p \\ \end{bmatrix} \\ A = \begin{bmatrix} sin(0^o) & sin(30^o) & \cdots & sin(150^o)\\ cos(0^o) & cos(30^o) & \cdots & cos(150^o)\\ \end{bmatrix} $$

Do the following matrix multiplication:

$$ R = PA $$

Some elements of the result:

$$ R_{0,0} = x_0 \cdot sin(0^o) + y_0 \cdot cos(0^o) \\ R_{0,1} = x_0 \cdot sin(30^o) + y_0 \cdot cos(30^o) \\ R_{1,0} = x_1 \cdot sin(0^o) + y_1 \cdot cos(0^o) \\ R_{1,1} = x_1 \cdot sin(30^o) + y_1 \cdot cos(30^o)  $$

This way we get all the r values by a single matrix multiplication. Turn it to MATLAB code:

% calculate sine and cosine maps
s = sin(theta);
c = cos(theta);

% do the matrix multiplication
rho = [x, y] * [s; c];
% d is the length of the diagonal of the image
[h, w] = size(image);
d = norm([h, w]) + 1;

% rho contains all the r values
% shift it by d to get only positive distances
% floor to get integer values rho = floor(rho + d);

In rho each row belongs to an edge point, describing the r values of its sinusoidal through θ. As written in the comments, all these values are positive integers, so they can be used directly as row indices for creating the two dimensional histogram. The column indices are simply from 1 to N – the number of the discrete angles – for each row.

For this step please have a look at the following post: Create two dimensional histograms in MATLAB.

The final step is creating the two dimensional histogram:

% create the two dimensional histogram
% repmat(1 : N, [P, 1]) creates the following matrix
% [1 2 3 4 ... N
%  1 2 3 4 ... N
%  . . . . ... .
%  1 2 3 4 ... N]
% this matrix is used for x coordinates
% the appropriate y coordinates are in matrix rho map = full(sparse(rho, repmat(1 : N, [P, 1]), 1));

Putting it all together:

function [theta, map] = hough(image, dtheta)
% angle map
theta = 0 : dtheta : pi;
N = length(theta);

% coordinates of edge points
[y, x] = find(image);
P = length(x);

% sine, cosine map
s = sin(theta);
c = cos(theta);

% the shifting value to get positive r values
[h, w] = size(image);
d = norm([h, w]) + 1;

% get all r values
rho = floor([x, y] * [s; c] + d);

% create the two dimensional histogram
map = full(sparse(rho, repmat(1 : N, [P, 1]), 1));
end

In theta we get back the quantized angle values, while map contains the two-dimensional histogram in which we have to find the peaks to identify the lines on the image.

         

New comment

comments powered by Disqus