harris.m bản nháp

function [cim, r, c] = harris(im, sigma, thresh, radius)
% Arguments:
%            im     - image to be processed.
%            sigma  - standard deviation of smoothing Gaussian. Typical
%                     values to use might be 1-3.
%            thresh - threshold (optional). Try a value ~50
%            radius - radius of region considered in non-maximal
%                     suppression (optional). Typical values to use might
%                     be 1-3.
%
% Returns:
%            cim    - binary image marking corners.
%            r      - row coordinates of corner points.
%            c      - column coordinates of corner points.
error(nargchk(2,5,nargin));
%dx = [-1 0 1; -1 0 1; -1 0 1];
dx = fspecial('sobel');
dy = dx';
Ix = conv2(im, dx, 'same');
Iy = conv2(im, dy, 'same');
g = fspecial('gaussian',max(1,fix(6*sigma)), sigma);
Ix2 = conv2(Ix.^2, g, 'same');
Iy2 = conv2(Iy.^2, g, 'same');
Ixy = conv2(Ix.*Iy, g, 'same');
%-----------
cim = (Ix2.*Iy2 - Ixy.^2)./(Ix2 + Iy2 + eps);
if nargin > 2
sze = 2*radius+1;
mx = ordfilt2(cim,sze^2,ones(sze));
cim = (cim==mx)&(cim>thresh);       % Find maxima
[r,c] = find(cim);
end
%-----------
max = 0;
[sy, sx] = size(Ix2);
K = zeros(sy, sx);
for i = 1:sx,
    for j = 1:sy,
        M = [Ix2(j,i), Ixy(j,i)
            Ixy(j,i), Iy2(j,i)];
        R(j,i) = det(M) - k*(trace(M))^2;
        if R(j,i) > max,
            max = R(j,i);
        end
    end
end