Harris 角点检测

时间:2023-10-19 21:39:13

一 、Motivation

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

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

Harris 角点检测

                           图 1.1

 

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

Harris 角点检测

Harris 角点检测

                         图 1.2

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

Harris 角点检测

Harris 角点检测

                        图 1.3

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

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

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

Harris 角点检测

                   图1.4

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

 

二、Mathematics representation

 

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

 

Harris 角点检测

                              图2.1

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

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

Harris 角点检测

                             图2.2

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

Harris 角点检测

                        图 2.3

 

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

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

Harris 角点检测的过程如下:

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))