Harris角點演算法實現

基本原理

人眼對角點的識別通常是在一個區域性的小區域或小視窗完成的。如果在各個方向上移動這個特徵的小視窗，視窗內區域的灰度發生了較大的變化，那麼就認為在視窗內遇到了角點。如果這個特定的視窗在影象各個方向上移動時，視窗內影象的灰度沒有發生變化，那麼視窗內就不存在角點；如果視窗在某一個方向移動時，視窗內影象的灰度發生了較大的變化，而在另一些方向上沒有發生變化，那麼，視窗內的影象可能就是一條直線的線段。

Harris角點演算法實現

1. 計算影象I(x,y)在X和Y兩個方向的梯度Ix、Iy
2. 計算影象兩個方向梯度的乘積
3. 使用高斯函式對Ix^2、Iy^2和Ixy進行高斯加權（取σ=1），生成矩陣M的元素A、B和C
4. 計算每個畫素的Harris響應值R，並對小於某一閾值t的R置為零
5. 在3×3或5×5的鄰域內進行非最大值抑制，區域性最大值點即為影象中的角點

Harris角點的性質

引數α對角點檢測的影響:增大α的值，將減小角點響應值R，降低角點檢測的靈性，減少被檢測角點的數量；減小α值，將增大角點響應值R，增加角點檢測的靈敏性，增加被檢測角點的數量。
Harris角點檢測運算元對亮度和對比度的變化不敏感:這是因為在進行Harris角點檢測時，使用了微分運算元對影象進行微分運算，而微分運算對影象密度的拉昇或收縮和對亮度的擡高或下降不敏感。換言之，對亮度和對比度的仿射變換並不改變Harris響應的極值點出現的位置，但是，由於閾值的選擇，可能會影響角點檢測的數量。
**Harris角點檢測運算元具有旋轉不變性**Harris角點檢測運算元使用的是角點附近的區域灰度二階矩矩陣。而二階矩矩陣可以表示成一個橢圓，橢圓的長短軸正是二階矩矩陣特徵值平方根的倒數。當特徵橢圓轉動時，特徵值並不發生變化，所以判斷角點響應值R也不發生變化，由此說明Harris角點檢測運算元具有旋轉不變性。
Harris角點檢測運算元不具有尺度不變性

function [posr,posc]=Harris1(in_image,a)
% 功能：檢測影象harris角點
% in_image-待檢測的rgb影象陣列
% a--角點引數響應，取值範圍：0.04~0.06
% [posr，posc]-角點座標
in_image=rgb2gray(in_image);
I=double(in_image);
%%%
 
%計算xy方向梯度%%%%%

fx=[-1,0,1];%x方向梯度模板
Ix=filter2(fx,I);%x方向濾波
fy=[-1;0;1];%y方向梯度模板(注意是分號)
Iy=filter2(fy,I);
%%%%計算兩個方向梯度的乘積%%%%%
Ix2=Ix.^2;
Iy2=Iy.^2;
Ixy=Ix.*Iy;
%%%%使用高斯加權函式對梯度乘積進行加權%%%%
%產生一個7*7的高斯窗函式，sigma值為2
h=fspecial('gaussian',[7,7],2);
IX2=filter2(h,Ix2);
IY2=filter2(h,Iy2);
IXY=filter2(h,Ixy);
%%%%%計算每個像元的Harris響應值%%%%%
[height,width]=size(I);
R=zeros(height,width);
%畫素(i,j)處的Harris響應值
for i=1:height
    for j=1:width
        M=[IX2(i,j) IXY(i,j);IXY(i,j) IY2(i,j)];
        R(i,j)=det(M)-a*(trace(M))^2;
    end
end
%%%%%去掉小閾值的Harris值%%%%%
Rmax=max(max(R));
%閾值
t=0.01*Rmax;
for i=1:height
    for j=1:width
        if R(i,j)<t
            R(i,j)=0;
        end
    end
end
%%%%%進行3*3領域非極大值抑制%%%%%%%%%
corner_peaks=imregionalmax(R);
%imregionalmax對二維圖片，採用8領域（預設，也可指定）查詢極值，三維圖片採用26領域
%極值置為1，其餘置為0
num=sum(sum(corner_peaks));
%%%%%%顯示所提取的Harris角點%%%%
[posr,posc]=find(corner_peaks==1);
figure
imshow(in_image);
hold on
for i=1:length(posr)
    plot(posc(i),posr(i),'r+');
end
end