Showing posts with label Harris Corner detector. Show all posts
Showing posts with label Harris Corner detector. Show all posts

Sunday, January 17, 2010

Rose Mary, she likes berries, believes in Fairies and is in love with Harry

We will explain here, how to carry out the Harris Corner Detector Algorithm in Matlab.
We will divide the task in 2 parts, one part will calculate the corner points of the image, and the other part, will draw small squares around those corner points with the purpose of having a way of displying them to the user.
The part of detecting the corner points of the image is the following:



%we create a function named harris that receives an rgb image
%and the desired k value. We will return a matrix named A that
%has information about whether or not a certain pixel represents
%a corner.
function [A] = harris(rgb, k)


%If we received an RGB image we convert it to Gray Scale.
%We do this, because if we were to work with an RGB image it
%would be necessary to work with 3 channels, the red, green and blue channel, %instead of just one, which is possible with the gray scale image .

if( size(rgb,3) >= 3 )
img = double( rgb2gray(rgb) );
end

%Before , we continue , it is important to recall that in the Harris Corner Detector, %we obtained a weighted sum of squared differences between an area (uv) and %another area that was obtained by shifting the original area by (x,y), and thus we %had the following form:

S(x,y) = \sum_u \sum_v w(u,v) \, \left( I(u,v) - I(u+x,v+y) \right)^2
%w(u,v) represented a weighted sum. It is important to note that I(u + x,v + y) %can be approximated by a Taylor expansion.
I(u+x,v+y) \approx I(u,v) + I_x(u,v)x+I_y(u,v)y
%Where Ix and Iy are the partial derivatives of I.
%Resulting in:
S(x,y) \approx \sum_u \sum_v w(u,v) \, \left( I_x(u,v)x + I_y(u,v)y \right)^2,
%Since we will carry out an operation that involves partial derivatives we need to %carry out a smoothing, because computation of derivatives generally involves a %stage of scale-space smoothing. For this,we will use the convolution of the gray %image with a Gaussian kernel.
So the next steps that we will carry out are:
a) Calculate the partial derivatives with respect to X and to respect to Y of the image. This will give us the gradients with respect to X and with respect to Y.
b) convolve these gradients with a Gaussian Kernel.
The following is our Gaussian Kernel, which will give us blurring on both directions%

g = 1/16 * [1 2 1; 2 4 2; 1 2 1];

%Here we define the Gradients Operators, we will use the Prewitt Gradient Kernel to obtain the gradients. If we pay more attention to this matrix, we notice that this kernel considers that the orthogonal and diagonal pixel differentials equally

dx = [-1 0 1; -1 0 1; -1 0 1]; %Prewitt Gradient Kernel in X
dy = dx';

% We obtain all the partial derivatives in x and y of the image. These are the %Gradient values.
Ix = conv2(img, dx, 'same');
Iy = conv2(img, dy, 'same');
% We obtain a matrix, that will have the product of Ix*Iy for all Ix and all Iy
Ixy= Ix .* Iy;

% We will now obtain the square value of Ix and of Iy and we will obtain a blurring of Ix square,Iy square and of Ixy.The blurring will be carried out by using the gaussian kernel.We need those square values , since let's recall we had the following:
I(u+x,v+y) \approx I(u,v) + I_x(u,v)x+I_y(u,v)y

Which produces the approximation


S(x,y) \approx \sum_u \sum_v w(u,v) \, \left( I_x(u,v)x + I_y(u,v)y \right)^2,


which can be written in matrix form:

S(x,y) \approx  \begin{pmatrix} x & y \end{pmatrix} A \begin{pmatrix} x \\ y \end{pmatrix},
where A is:

A = \sum_u \sum_v w(u,v) \begin{bmatrix} I_x^2 & I_x I_y \\ I_x I_y & I_y^2  \end{bmatrix} = \begin{bmatrix} \langle I_x^2 \rangle & \langle I_x I_y \rangle\\ \langle I_x I_y \rangle & \langle I_y^2 \rangle \end{bmatrix}
%So we calculate the values that are inside this matrix.

Ix2 = conv2(Ix .* Ix, g, 'same');
Iy2 = conv2(Iy .* Iy, g, 'same');
Ixy = conv2(Ixy, g, 'same');


%We now have 3 matrices , Ix2 that hold the values of all the xs in the image, but they are squared and have been convolved with a Gaussian, so all the xto the square values are blurred a bit. It is important to note, that this operation helps to reduce noise. It smooths the image. We have another matrix Iy2 with the ys to the square value and also blurred as well as third matrix Ixy that holds the values of all the x's times their corresponding y. All 3 of these matrices will help us to calculate the Harris corner response that each pixel of the image has.
The Harris Corner response for each pixel will come out of a matrix A, that as we had said before handwas conformed of:
A = \sum_u \sum_v w(u,v) \begin{bmatrix} I_x^2 & I_x I_y \\ I_x I_y & I_y^2  \end{bmatrix} = \begin{bmatrix} \langle I_x^2 \rangle & \langle I_x I_y \rangle\\ \langle I_x I_y \rangle & \langle I_y^2 \rangle \end{bmatrix}
%Where Ix represents the x value of the pixel, and Iy the y value of the pixel. The Harris Corner response %of a pixel will be Mc:
M_c = \lambda_1 \lambda_2 - \kappa \, (\lambda_1 + \lambda_2)^2 = \operatorname{det}(A) - \kappa \, \operatorname{trace}^2(A)
%What our Harris Function will return will be a matrix that will hold the harris corner response for each %one of the pixels that conform the image:



A = (Ix2.*Iy2 - Ixy.^2) - k * (Ix2 + Iy2).^2;

end


We will now write the function that draws a small red square around the detected corner points.

%This function receives the image we wish to detect and draw the corner points of, %it also receives the desired k value to use, and the desired threshold. It is %important to remember that in the Harris corner detector, we consider a corner %to be a corner when the measure of corner response surpasses a certain %threshold. This measure is computed through the determinant and the trace of the % matrix.

function img_h = project1(img, k, threshold)
%We store the original image in img_h, we need to store it, since we will draw the %squares denoting the corners above it.
img_h = img;
%We use the function we had defined above and with it obtain a matrix that holds %all the corner response measures for all the pixels of the image%
M = my_harris(img, k);

%We iterate through the whole matrix
for x = 2 : size( M, 1 )
for y = 2 : size(M, 2)
%If we find a point of the matrix, that has a value above the threshold, then that %point is a corner and we will draw a rectangle on that pixel.
if M(x,y) > threshold
for xpos = x - 1:x+1
img_h(xpos,y-1,1) = 255;
img_h(xpos,y+1,1) = 255;
end

for ypos = y - 1:y+1
img_h(x-1,ypos,1) = 255;
img_h(x+1,ypos,1) = 255;
end
end
end
end

end



We now present a image that was used for this purpose.
The original image is:

And the image with the corners detected is:

Saturday, January 16, 2010

The Harris Marris whose not Ferris post...

In many robotic problems it is necessary for the machine to be able to detect the depth and distance of certain objects sometimes to avoid obstacles and other times to retrieve with precision a certain object from the scene.
To accomplish this, the robot is usually equipped with 2 cameras that take pictures of the environment they are in. These 2 cameras commonly hold a distance from each other, a distance similar to the one we present with our eyes, due to it, the pictures from one camera are slightly shifted with respect to the ones taken by the other camera. This shift is usually denominated disparity and is what the computer uses to know whether an object is close by or far away.
One major problem that the machine encounters while carrying out this task, is how does it detect that the red cup in picture 1 is moved , let's say, 4 cm with respect to where the red cup is in picture 2?
These vision tasks require finding corresponding features across 2 or more views.


Therefore the first necessary step is to find the features of a scene. But how do we do this?

What we can start doing is making image patches.Elements to be matched are image patches of a fixed size..The task is therefore to find the best (most similar) patch in the second image, it is clear that the chosen patch should be very distinctive (there should only be one patch in the second picture that looks similar). One good patch is one that presents large variation in the neighborhood of a point in all directions.
For example Take the following 2 images:








A good patch image patch could be:


while a bad one, because it has many matchings is:


We are looking for stable features over changes of view points. One type of features that maintain this type of characteristic are Corners.
The Harris Corner detector provides a mathematical tool for finding them.
With an image patch, we can have the following cases:
a) The patch represents a 'flat' zone.
b)The patch represents an edge .
c) The patch represents a corner .



A Flat region as we can see from the above image, presents no change in all directions, an edge presents no change along the edge direction, and a corner presents significant change in all directions. This means that if we shift the window of where we are gathering the patch image, we should perceive a large change in appearance.
The Harris Corner Detector gives us a way to determine which of the above cases hold.
But how does it do it exactly???!

The Harris Corner Detector utilizes the following expression:
E (u,v)=∑ W(Xi,Yi)[I1(Xi+U,Yi+V)-I0(Xi,Yi)]^2
W(Xi,Yi) is a window function. Which sets:


I0(Xi,Yi) is the intensity that is present in the pixel located in (Xi,Yi) and I1(Xi+U,Yi+V) is the intensity located in the pixel Xi+U,Yi+V it is called the shifted or displaced Intensity.
It is easy to see that to detect corners, we want points where E(u,v) is very large.
Using Taylor's first order approximation and matrix algebra the above expression can be rewritten as:


Where M is a 2x2 matrix computed from image derivatives.


The classification can be carried out by analyzing the eigenvalues of the M matrix.


The measure of the corner response is actually set by:
R=(determinant of M)+k(trace of M)^2
Determinant of M=λ1λ2
Trace of M=λ1+λ2
K is an empirical constant that varies from .04-.06

For corners R tends to be very large.
For edges, R tends to be a very large negative number.
For flat areas R tends to zero.