Implementation of the Harris Corner Detection Algorithm

harris_corner_detector.m

function [ response, corners ] = harris_corner_detector( image, varargin )

%   MATLAB FUNCTION DEFINITIONS
%       function [ response, corners ] = harris_corner_detector( image, varargin )
%    MATLAB CALLING SEQUENCE
%       response = harris_corner_detector( image )
%       response = harris_corner_detector( image, sigma )
%       [ response, corners ] = harris_corner_detector( image )
%       [ response, corners ] = harris_corner_detector( image, sigma)  





% Sort out the inputs. Parse inputs based on the syntax, this way sigma can
% be defined. 
if nargin         == 1
    sigma         =  1
elseif nargin     == 2
    sigma         = varargin{1'};
end

% Set the inputs.
% Make sure that the iamge is a double.
image_in          = image;
image             = double(image_in);

% Minimum points within the image to start running the filter over
min_coord_point   = 10;

% Define the first set of points for the iamge to go over. These images
% will be referenced again and reset to a new variable value once they are
% evaluated by the variable Sigma.
y_origin          = min_coord_point;
y_edge            = 255;
x_origin          = min_coord_point;       
x_edge            = 255;

% Find the minimum and maximum values as they pertain to sigma. Sigma is
% pre-defined in your function call. This is the step in which the values
% are evaluated by sigma.
c_origin          =x_origin - sigma; 
c_edge            =x_edge   + sigma; 
r_origin          =y_origin - sigma; 
r_edge            =y_edge   + sigma;

% Reset the variable for organization sake. It variable strings started 
% to get unreasonably long if left the way they were.   
r                 =min_coord_point; 


% Bring in the filter here that Carl references on the slides. This is a
% basic edge detection filter, generally used and widely used for edge
% detection.
mask_filter       = [-1 0 1;
                     -1 0 1; 
                     -1 0 1]; 

                 
% Take a line of the mask filter from above to be used as the X vector.                 
mask_x_vector     = [-1 0 1];


% The same filter applies to the Y dimensions of the image but reshape
% using the tic operator to serve as the new Y vector.
mask_y_vector     = mask_x_vector';


% Use the convolve filter function to move across the image from the
% pre-defined origin point. Make sure to keep the 'same' var
IDXx              = conv2(image(x_origin:x_edge,y_origin:y_edge),...
                    mask_x_vector, 'same');   
IDXy              = conv2(image(x_origin:x_edge,y_origin:y_edge),...
                    mask_y_vector, 'same');
 

% Assemble a basic filter to run over the image. This set of equations takes the 
% equations that Carl and the forums provided to find the general
% equation and turn it into a matlab function to run here.
g1                = (-mask_filter^2);
g2                = (2*sigma^2);
g3                = (sigma*sqrt(2*pi));
% Preliminary variables defined, now run the math
g4                = (g1./g2);
g5                = (g4./g3);
    
gaussian_filter   = exp(g5);


% Use the colvolve function again to perform the 2D convolution of the indexed 
% R values and indexed C values 
A              = conv2(IDXx.^2,  gaussian_filter);  
B              = conv2(IDXy.^2,  gaussian_filter);
C              = conv2(IDXx.*IDXy, gaussian_filter);
    
% Displacement for the points on the sides. Upon advice from Carl, use
% either 1 or 0.01 for the K value.
k                 = 0.01;

% Run all of these values through the noise reduction formula provided in
% the forum. I wrote it as a separate function to make the code run faster,
% also to make it easier to debug.
[ R ] = noise_reduction_sigma(A, B, C, k);

    
% Take the maximum value of the value R that was found using Carl's equation
% and assign it to the new variable.
maximum_Range_Values = max(max(R));
max_size_values      = r; 


% Important step, take the eigenvalues of each of the determined values for
% A, B and C. This is what allows the points to be mapped by finding their
% relative distances across the filter. The eigen matlab function was used
% to save time, it was also ok'd by Carl to be used.

eigenA = eig(A);
eigenB = eig(B);
eigenC = eig(C);

% Re-classify the variable r back to its original call name, for better
% organization during . I wrote two sub functions here, one using the
% ordfilt function to insert each of my values
overall_size      = r; 
[ filtered_area ] = order_operation(R,overall_size);    
final_R_values    = (R == filtered_area)&(R>0); 
[ R ]             = harris_computer(R, final_R_values);


% Note that the variable 'lowercase r is recycled here'
% This was a huge boo boo on my part, never going to do that again but it
% was easier to leave it since it made more sense in my head than changing
% it to something more obscure, after all of this time referring to rows
% and columns as [r, c]


[r,c]             = find(R);
corner_points     =[r+x_origin,...
                    c+y_origin]; 
% Begin indexing the loop to run the algorithm over the entire image from 
% origin to end.  
idx_total         =size(corner_points,1);   
    
    for r         =1: idx_total
       % Run through the rows   
       image(corner_points(r,1):corner_points(r,1)+2,...
           corner_points(r,2))  = 256;
   
       % Run through the rows        
       image(corner_points(r,1):corner_points(r,1)+2,...
           corner_points(r,2))  = 256;
   
       % Run through the cols
       image(corner_points(r,1),...
           corner_points(r,2)-2:corner_points(r,2))  = 256;
   
       % Run through the cols
       image(corner_points(r,1),...
           corner_points(r,2)-2:corner_points(r,2))  = 256;
    end
   
    
   % Find all of the date for representing
   % I chose the variable U and O because every other row x col letter set
   % has already been used.
   
   [U, O]            = size(corner_points);
   number_of_corners = U
   corners           = corner_points; 
   
   % Display the response, which now is just a map of the points over the
   % corners.
   response          = imshow(R, []);
   
   
   
end