一、原理

二、實現

close all;
clear all;
I=imread('test.tif');
[posX,posY]=harris(I);
figure;imshow(I);
hold on; plot(posX, posY, 'g*');

function [posX,posY]=harris(I)
%Harris角點檢測
%I:輸入影象
%posX:角點X座標
%posY:角點Y座標
I=double(I);
[m,n]=size(I);
hx=[-1,0,1;-1,0,1;-1,0,1];
Ix=imfilter(I,hx,'replicate','same');%X方向差分影象
Iy=imfilter(I,hx','replicate','same');%Y方向差分影象
Ix2=Ix.^2;
Iy2=Iy.^2;
Ixy=Ix.*Iy;
h=fspecial('gaussian',3,2);
Ix2=imfilter(Ix2,h,'replicate','same');%高斯濾波
Iy2=imfilter(Iy2,h,'replicate','same');
Ixy=imfilter(Ixy,h,'replicate','same');
R=zeros(m,n);
k=0.04;%建議值(0.04--0.06)
for i=1:m
    for j=1:n
       R(i,j)=(Ix2(i,j)*Iy2(i,j)-Ixy(i,j)*Ixy(i,j))-k*((Ix2(i,j)+Iy2(i,j))^2);%角點響應值
    end
end
T=0.1*max(R(:));% 閾值，可控制返回的角點個數
result=zeros(m,n);
%非極大值抑制（3*3視窗中大於閾值T的區域性極大值點被認為是角點）
for i=2:m-1
    for j=2:n-1
        tmp=R(i-1:i+1,j-1:j+1);
        tmp(2,2)=0;
        if(R(i,j)>T&&R(i,j)>max(tmp(:)))
            result(i,j)=1;
        end
    end
end
[posY,posX]=find(result);