Harris Corner（Harris角檢測）

在做圖像匹配時，常需要對兩幅圖像中的特征點進行匹配。為了保證匹配的准確性，所選擇的特征必須有其獨特性，角點可以作為一種不錯的特征。

那么為什么角點有其獨特性呢？角點往往是兩條邊緣的交點，它是兩條邊緣方向變換的一種表示，因此其兩個方向的梯度變換通常都比較大並且容易檢測到。

這里我們理解一下Harris Corner 一種角點檢測的算法

角點檢測基本原理：

人們通常通過在一個小的窗口區域內觀察點的灰度值大小來識別角點，如果往任何方向移動窗口都會引起比較大的灰度變換那么往往這就是我們要找的角點。如下圖右

下面我們看一下Harris的數學公式，對於[x,y]平移[u,v]個單位后強度的變換有下式,I(x+u,y+v)是平移后的強度，I(x,y)是原圖像像素。對於括號里面的值，如果是強度恆定的區域，那么它就接近於零，反之如果強度變化劇烈那么其值將非常大，所以我們期望E(u,v)很大。

 

 其中w是窗函數，它可以是加權函數，也可以是高斯函數

利用二維泰勒展開式我們有

 

 所以其中一階可以近似為

於是我們可以給出Harris Corner的如下推導，其中Ix，Iy是x，y方向的Gradient模，乘以位移得到位移后的量

對於小的位移，我們可以用雙線性插值方法近似：

 

其中M為2*2矩陣如下

 在本質上我們可以把二次項看成一個橢圓函數，我們對M進行特征值分析有λ1，λ2

 

根據λ1，λ2的值我們可以把其分為三類：

1.λ1，λ2都很小且近似，E在所以方向接近於常數；

2.λ1>>λ2,或者λ2>>λ1, E將在某一方向上很大；

3.λ1，λ2都很大且近似，E將在所以方向上很大；

如圖所示：

最后我們通過計算角點響應值R來判斷其屬於哪個區間

其中k一般為常數取在0.04-0.06間。

 

算法步驟：

1.計算圖像x,y方向的梯度I_x，I_y

2.計算每個像素點的梯度平方

3.計算梯度在每個像素點的和

4.定義在每個像素點的矩陣H，也就是前面的M

5.計算每個像素的角點響應

6.設置閾值找出可能點並進行非極大值抑制

 

代碼：

close all
clear all

I = imread('empire.jpg');
I = rgb2gray(I);
I = imresize(I,[500,300]);
imshow(I);

sigma = 1;
halfwid = sigma * 3;

[xx, yy] = meshgrid(-halfwid:halfwid, -halfwid:halfwid);

Gxy = exp(-(xx .^ 2 + yy .^ 2) / (2 * sigma ^ 2));
Gx = xx .* exp(-(xx .^ 2 + yy .^ 2) / (2 * sigma ^ 2));
Gy = yy .* exp(-(xx .^ 2 + yy .^ 2) / (2 * sigma ^ 2));

%%apply sobel in herizontal direction and vertical direction compute the
%%gradient
%fx = [-1 0 1;-1 0 1;-1 0 1];
%fy = [1 1 1;0 0 0;-1 -1 -1];
Ix = conv2(I,Gx,'same');
Iy = conv2(I,Gy,'same');
%%compute Ix2, Iy2,Ixy
Ix2 = Ix.*Ix;
Iy2 = Iy.*Iy;
Ixy = Ix.*Iy;

%%apply gaussian filter
h = fspecial('gaussian',[6,6],1);
Ix2 = conv2(Ix2,h,'same');
Iy2 = conv2(Iy2,h,'same');
Ixy = conv2(Ixy,h,'same');
height = size(I,1);
width = size(I,2);
result = zeros(height,width);
R = zeros(height,width);
Rmax = 0;
%% compute M matrix and corner response
for i = 1:height
    for j =1:width
        M = [Ix2(i,j) Ixy(i,j);Ixy(i,j) Iy(i,j)];
        R(i,j) = det(M) - 0.04*(trace(M)^2);
        if R(i,j)> Rmax
            Rmax = R(i,j);
        end
    end
end
%% compare whith threshold
count = 0;
for i = 2:height-1
    for j = 2:width-1
        if R(i,j) > 0.01*Rmax
            result(i,j) = 1;
            count = count +1;
        end
    end
end

%non-maxima suppression 
result = imdilate(result, [1 1 1; 1 0 1; 1 1 1]);

[posc,posr] = find(result == 1);
imshow(I);
hold on;
plot(posr,posc,'r.');

 

 

本文原創，轉載請注明出處