function features = harris_detector(input_image, varargin)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%input parameters:
%input_image grayscale image
%segment_length input the block size for weighting the image for filter
%k determine the parameter of value k formula H = det(G)-k*tr(G)^2 G present the harris matrix
%tau Threshold for H the point is a feature if H > tau
% do_plot true indicate plot the corners on the image
%min_dist the minimal distance of two feature points
%tile_size the image will be divided into blocks with size tile_size, in one block int can only exist N features
%N the maximal features in a block
%output parameter:
% features the locations of the features
% the names and values of the parameters should be given in pairs
%default value
%'segment_length',15,'k',0,05,'tau',1000000,'do_plot',faluse,'min_dist',20,'tile_size',[200,200],'N',5
% by 南市吃主爷
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% set constraints of input values, defined as optioinal
p = inputParser;
addParameter(p,'segment_length',15,@(x)validateattributes(x,{'numeric'},{'odd','>',1}))
addParameter(p,'k',0.05,@(x)validateattributes(x,{'numeric'},{'>',0,'<',1}))
addParameter(p,'tau',1000000,@(x)validateattributes(x,{'numeric'},{'>',0}))
addParameter(p,'do_plot',false,@(x)validateattributes(x,{'logical'},{}))
addParameter(p,'min_dist',20,@(x)validateattributes(x,{'numeric'},{'>',1}));
addParameter(p,'tile_size',[200,200],@(x)validateattributes(x,{'numeric'},{}));
addParameter(p,'N',5,@(x)validateattributes(x,{'numeric'},{'>=',1}));
parse(p,varargin{:})
% if tile_size has only one value, the height is equal to width
if size(p.Results.tile_size) == [1,1]
tile_size = [p.Results.tile_size, p.Results.tile_size];
else
tile_size = p.Results.tile_size;
end
min_dist = p.Results.min_dist;
N = p.Results.N;
segment_length = p.Results.segment_length;
k = p.Results.k;
tau = p.Results.tau;
do_plot = p.Results.do_plot;
[len,col,de] = size(input_image);
% test if the image is grayscale
if de ~= 1
error('Image format has to be NxMx1');
end
% Approximation of the image gradient
input_image = double(input_image);
sobel_x = [1,0,-1;2,0,-2;1,0,-1];
Ix = conv2(input_image,sobel_x,'same');
sobel_y = [1,2,1;0,0,0;-1,-2,-1];
Iy = conv2(input_image,sobel_y,'same');
% Weighting
w = fspecial('gaussian',[1,segment_length],segment_length/3);
w_2 = w'*w;
% Harris Matrix G
G11 = conv2(Ix.*Ix,w_2,'same');
G12 = conv2(Ix.*Iy,w_2,'same');
G22 = conv2(Iy.*Iy,w_2,'same');
det_G = G11.*G22-G12.^2;
trace_G = G11+G22;
% calculate the matrix
H = det_G-k*trace_G.^2;
[len_H,col_H] = size(H);
% add a border of size segment_length/2
border_segment = ceil(segment_length/2);
H = padarray(H,[border_segment,border_segment]);
% extract the features which H > tau
sub_H_lo = gt(H,tau);
sub_H_lo1 = sub_H_lo(border_segment+1:border_segment+len_H,border_segment+1:border_segment+col_H);
%[x,y] = find(sub_H_lo1);
%features = [y';x'];
corners = H.*sub_H_lo;
corners = corners(border_segment+1:border_segment+len_H,border_segment+1:border_segment+col_H);
[len_corners,col_corners] = size(corners);
%add a border of size min_dist
%corners_b = zeros(2*min_dist+len_corners,2*min_dist+col_corners);
%corners_b(min_dist+1:min_dist+len_corners,min_dist+1:min_dist+col_corners) = corners;
%corners = corners_b;
corners = padarray(corners,[min_dist,min_dist]);
%sort the corners by the intensities of features
[sorted_corners,sorted_index] = sort(corners(:),'descend');
greater = find(sorted_corners);
sorted_index = sorted_index(greater);
% calculate the block size and divide the image into blocks
num_x = ceil(len/tile_size(1));
num_y = ceil(col/tile_size(2));
acc_array = zeros(num_x,num_y);
% linear index to 2-D index
[ix,iy] = ind2sub(size(corners),sorted_index);
ct = 1;
%location of two features should be larger than min_dist
Cakes = ones(min_dist*2+1,min_dist*2+1);
for i = 1:min_dist*2+1
for j = 1:min_dist*2+1
dis = sqrt((min_dist+1-i)^2+(min_dist+1-j)^2);
if dis <= min_dist
Cakes(i,j) = 0;
end
end
end
Cakes = logical(Cakes);
% the center feature can be kept.
Cakes(min_dist+1,min_dist+1) = 1;
for i = 1:numel(sorted_index)
temp = corners(ix(i)-min_dist:ix(i)+min_dist,iy(i)-min_dist:iy(i)+min_dist);
if corners(ix(i),iy(i) ) ~= 0
temp = temp.*Cakes;
corners(ix(i)-min_dist:ix(i)+min_dist,iy(i)-min_dist:iy(i)+min_dist) = temp;
end
end
corners = corners(min_dist+1:end-min_dist,min_dist+1:end-min_dist);
%every block only contains N features
[sorted_corners,sorted_index1] = sort(corners(:),'descend');
greater1 = find(sorted_corners);
sorted_index1 = sorted_index1(greater1);
[x1,y1] = ind2sub(size(corners),sorted_index1);
for po = 1:numel(sorted_index1)
block_x = ceil((x1(po))/tile_size(1));
block_y = ceil((y1(po))/tile_size(1));
if (acc_array(block_x,block_y) < N)
if (corners(x1(po),y1(po)) ~= 0)
acc_array(block_x,block_y) = acc_array(block_x,block_y)+1;
features(1:2,ct) = [x1(po);y1(po)];
ct = ct+1;
end
else
corners(x1(po),y1(po)) = 0;
end
end
[x1,y1] = find(corners);
features = [y1';x1'];
if do_plot
input_image = uint8(input_image);
imshow(input_image);
hold on
plot(features(1,:),features(2,:),'go')
end
end
