哈里斯角点检测

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