Harris 角点检测

一 、Motivation

对于做图像处理的人来说，Harris角点检测肯定听过，1988年发表的文章“A combined corner and edge detector”描述了这种角点检测方法，这篇论文朴实无华，对于图像处理入门来说，非常值得读一读。

Harris角点检测的提出是图像匹配问题的需求，在立体视觉（stereo vision)和运动估计（motion estimation）中，常常需要在两个view（立体视觉）或者同一视频的两帧（运动估计）中找到对应的特征（correspondence feature），如下图1.1所示。

                           图 1.1

以patch matching 为例，若在两个view中提取出来的patch 如下图1.2，那么匹配两幅图中相似的patch是比较容易的，

                         图 1.2

而如果两个view中提取出来的patch如下图1.3，那么匹配就不那么容易了，

                        图 1.3

为什么呢？因为图1.2中的patch很独特，信息丰富，图1.3 中的patch单独看来毫无特点，极易混淆。我们称图1.2中的特征为“好特征”，图1.3中的特征是”坏特征“。

那么什么是好特征，什么是坏特征？我认为有两个要考虑的：1 稳定，对缩放，视角变换，光线变化等稳定 2.易区分 。

角点就具有这样的特征，角点如何描述，请看图1.4，

                   图1.4

上图具体解释是这样的，给定一个窗口，如果包含角点，那么这个窗口平移（u，v）个单位，不管这个平移是往哪个方向，窗口中像素对应位置的变化都比较大，而如果包含的是一条边缘，在沿着边缘平移窗口时，窗口中像素强度变化基本没有，而垂直于边缘移动时，变化强烈，对于平坦区域，怎么移动都没有多大变化，当然，这里的平移都是小范围平移。

二、Mathematics representation

数学描述这种强度变化如下图2.1.

                              图2.1

可以看到，这个公式表示往各个方向移动时强度变化的累加和，控制w就可以控制平移后强度累加的方式。

然后用一级泰勒展开近似I(x+u,y+v)-I(x,y)，并将上图的公式用矩阵的形式表达出来，有：

                             图2.2

最后 E(u,v) 可以表示为:

                        图 2.3

注意到此时M是对称矩阵，可以表示为M = Q A Q^T 的形式 ，A 为对角矩阵，因此，对角中即为M的特征值，因此，一定要M的2个特征值都比较大才能保证E总是很大。

实际计算过程中，用高斯核来表示w(x,y).

Harris 角点检测的过程如下：

                                图2.4

需要注意的是，求微分图像和第三步的W矩阵都是可以调节或者换成其他形式的，W换成高斯核主要是利用了它各项同性的性质。

三、implementation

最终matlab代码（转自网上）实现如下：

% Harris detector
% The code calculates
% the Harris Feature Points(FP)
%
% When u execute the code, the test image file opened
% and u have to select by the mouse the region where u
% want to find the Harris points,
% then the code will print out and display the feature
% points in the selected region.
% You can select the number of FPs by changing the variables
% max_N & min_N
% A. Ganoun

load Imag

I =double(frame);
%****************************
imshow(frame);
k = waitforbuttonpress;
point1 = get(gca,'CurrentPoint');  %button down detected
rectregion = rbbox;  %%%return figure units
point2 = get(gca,'CurrentPoint');%%%%button up detected
point1 = point1(1,1:2); %%% extract col/row min and maxs
point2 = point2(1,1:2);
lowerleft = min(point1, point2);
upperright = max(point1, point2);
ymin = round(lowerleft(1)); %%% arrondissement aux nombrs les plus proches
ymax = round(upperright(1));
xmin = round(lowerleft(2));
xmax = round(upperright(2));
%***********************************
Aj=6;
cmin=xmin-Aj; cmax=xmax+Aj; rmin=ymin-Aj; rmax=ymax+Aj;
min_N=12;max_N=16;
%%%%%%%%%%%%%%Intrest Points %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
sigma=2; Thrshold=20; r=6; disp=1;
dx = [-1 0 1; -1 0 1; -1 0 1]; % The Mask
    dy = dx';
    %%%%%%
    Ix = conv2(I(cmin:cmax,rmin:rmax), dx, 'same');
    Iy = conv2(I(cmin:cmax,rmin:rmax), dy, 'same');
    g = fspecial('gaussian',max(1,fix(6*sigma)), sigma); %%%%%% Gaussien Filter

    %%%%%
    Ix2 = conv2(Ix.^2, g, 'same');
    Iy2 = conv2(Iy.^2, g, 'same');
    Ixy = conv2(Ix.*Iy, g,'same');
    %%%%%%%%%%%%%%
    k = 0.04;
    R11 = (Ix2.*Iy2 - Ixy.^2) - k*(Ix2 + Iy2).^2;
    R11=(1000/max(max(R11)))*R11;
    R=R11;
    ma=max(max(R));
    sze = 2*r+1;
    MX = ordfilt2(R,sze^2,ones(sze));
    R11 = (R==MX)&(R>Thrshold);
    count=sum(sum(R11(5:size(R11,1)-5,5:size(R11,2)-5)));

    figure;plot(R11);
    loop=0;
    while (((count<min_N)|(count>max_N))&(loop<30))
        if count>max_N
            Thrshold=Thrshold*1.5;
        elseif count < min_N
            Thrshold=Thrshold*0.5;
        end

        R11 = (R==MX)&(R>Thrshold);
        count=sum(sum(R11(5:size(R11,1)-5,5:size(R11,2)-5)));
        loop=loop+1;
    end


    R=R*0;
    R(5:size(R11,1)-5,5:size(R11,2)-5)=R11(5:size(R11,1)-5,5:size(R11,2)-5);
    [r1,c1] = find(R);
    PIP=[r1+cmin,c1+rmin];%% IP


   %%%%%%%%%%%%%%%%%%%% Display

   Size_PI=size(PIP,1);
   for r=1: Size_PI
   I(PIP(r,1)-2:PIP(r,1)+2,PIP(r,2)-2)=255;
   I(PIP(r,1)-2:PIP(r,1)+2,PIP(r,2)+2)=255;
   I(PIP(r,1)-2,PIP(r,2)-2:PIP(r,2)+2)=255;
   I(PIP(r,1)+2,PIP(r,2)-2:PIP(r,2)+2)=255;

   end

   imshow(uint8(I))