Harris算子的运用 用于图像配准

Harris算子介绍：

该算子是C．Harris和M．J．Stephens在1988年提出的一种点特征提取算子。这种算子受信号处理中自相关函数的启发，可以给出图像中某一像素点的自相关矩阵肘，其特征值是自相关函数的一阶曲率，如果算，Y两个方向上的曲率值都高，那么就认为该点是角点。Harris角点检测算子

Harris[2]角点检测算子是Moravec角点检测算子的改进.

（1）算子用高斯函数代替二值窗口函数,对离中心点越近的像素赋于越大的权重，以减少噪声影响。

                                                   图1-3高斯函数

(2)算子只考虑了每隔45度方向，Harris算子用Taylor展开去近似任意方向。

写成矩阵形式: :

 

                            式子（1-2）

                           式子（1-3）

式中，Ix为x方向的差分，Iy为y方向的差分，w(x,y)为高斯函数。

(3)Harris采用了一种新的角点判定方法。矩阵M的两个特征向量l1和l2与矩阵M的主曲率成正比。Harris利用l1, l2来表征变化最快和最慢的两个方向.若两个都很大就是角点,一个大一个小就是边缘,两个都小就是在变化缓慢的图像区域.

                                                                                来自文献[11]

图1- 4用矩阵Ｍ的特征向量分类图像像素点

但是解特征向量需要比较多的计算量,且两个特征值的和等于矩阵M的迹，两个特征值的积等于矩阵M的行列式。所以用(1-4）式来判定角点质量。(k常取0.04-0.06)

   （1-4）

(4) Harris算法总结

Step 1:对每一像素点计算相关矩阵M。

Step 2:计算每像素点的Harris 角点响应。

Step 3.在w*w范围内寻找极大值点，若Harris 角点响应大于阀值，则视为角点。

Harris算子对灰度的平移是不变的,因为只有差分，对旋转也有不变性，但是对尺度很敏感,在一个尺度下是角点, 在在另一个尺度下可能就不是了.

                  图1- 5　harris算子对尺度的敏感性

                    图1- 6　harris算子对简单图像的响应

               Harris 算子是一种有效的点特征提取算子,其优点总结起来有:

              ①计算简单:Harris 算子中只用到灰度的一阶差分以及滤波,操作简单。

       ②提取的点特征均匀而且合理:Harris 算子对图像中的每个点都计算其兴趣值,然后在邻域中选择最优点。实验表明,在纹理信息丰富的区域,Harris 算子可以提取出大量有用的特征点,而在纹理信息少的区域,提取的特征点则较少。

       ③稳定:Harris算子的计算公式中只涉及到一阶导数，因此对图像旋转、灰度变化、噪声影响和视点变换不敏感,它也是比较稳定的一种点特征提取算子。

             Harris 算子的局限性有：

      ①它对尺度很敏感，不具有尺度不变性。

      ②提取的角点是像素级的。

代码：

<span style="font-size:18px;">im=imread('lena.jpg');tau=100;im=double(im);keyXs=[];keyYs=[];win=3;[height,width] = size(im);result = zeros(height,width);%Then I will get the gradients of the image along the x and y axises.sobel_x=1/4*[-1 0 1;-2 0 2;-1 0 1];sobel_y=1/4*[-1 0 1;-2 0 2;-1 0 1]';diffx=imfilter(im,sobel_x);         %对图像x方向进行梯度diffy=imfilter(im,sobel_y);       %对图像y方向的梯度进行计算%For smoothing the differentiation of the image along the x and y%direction, the gauss filter of the diffx and diffy is must.gauss_win=win;sigma=1;[x,y]=meshgrid(-gauss_win:gauss_win,-gauss_win:gauss_win);gauss2D=exp(-(x.^2+y.^2)/(2*sigma.^2));  %产生高斯算子gauss2D=gauss2D/(sum(sum(gauss2D)));  %对高斯算子进行归一化%Then calculate the M matrix.A=imfilter(diffx.*diffx,gauss2D);      %二阶x方向梯度进行高斯滤波B=imfilter(diffy.*diffy,gauss2D);      %二阶y方向梯度进行高斯滤波C=imfilter(diffx.*diffy,gauss2D);      %对图像x y方向的梯度进行高斯滤波%Harris mehtods.if(strcmp(tau,'Harris'))k=0.2;threshold=200;%Harris criteria.Harris=A.*B-C.^2-k*(A+B).^2;         %Then I will do the non-maximum supression.supress_win=2;points_count=0;   %%对图像的每个像素进行阈值判断是否是角点for x=supress_win+1:width-supress_win    for y=supress_win+1:height-supress_win        %Then you need to judge if the pixel has the biggest Harris        %response in the (2*supress_win+1)*(2*supress_win+1) neighbour.        temp=Harris(y,x);        %得到图像x  y位置的harri值        if(temp>threshold)          %该点的haari值大于周围像素的阈值时            flag=0;            for i=-supress_win:supress_win                for j=-supress_win:supress_win                    if(temp>=Harris(y+j,x+i))                        flag=flag+1;              %像素的个数加1                    end                end            end            if(flag==((2*supress_win+1)*(2*supress_win+1)))                result(y,x)=1;                points_count=points_count+1;                keyXs(points_count)=x;                keyYs(points_count)=y;    %存储haari角点的坐标            end        end    endend</span>

  另一种代码：

<span style="font-size:18px;">%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   Harris角点提取算法                                          %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%clear;%filename='487_r.png';%X= imread('Circle.bmp');     % 读取图像X= imread('001.jpg');     % 读取图像%imshow(X);%Info=imfinfo(filename);% if Info.BitDepth>8%     f=rgb2gray(X);f=X;%end%% fx = [5 0 -5;8 0 -8;5 0 -5];          % 高斯函数一阶微分，x方向(用于改进的Harris角点提取算法)ori_im=double(f)/255;                   %unit8转化为64为双精度double64fx = [-2 -1 0 1 2];                     % x方向梯度算子(用于Harris角点提取算法)Ix = filter2(fx,ori_im);                % x方向滤波% fy = [5 8 5;0 0 0;-5 -8 -5];          % 高斯函数一阶微分，y方向(用于改进的Harris角点提取算法)fy = [-2;-1;0;1;2];                     % y方向梯度算子(用于Harris角点提取算法)Iy = filter2(fy,ori_im);                % y方向滤波Ix2 = Ix.^2;Iy2 = Iy.^2;Ixy = Ix.*Iy;clear Ix;clear Iy;h= fspecial('gaussian',[7 7],2);        % 产生7*7的高斯窗函数，sigma=2Ix2 = filter2(h,Ix2);Iy2 = filter2(h,Iy2);Ixy = filter2(h,Ixy);height = size(ori_im,1);width = size(ori_im,2);result = zeros(height,width);           % 纪录角点位置，角点处值为1R = zeros(height,width);Rmax = 0;                              % 图像中最大的R值for i = 1:height    for j = 1:width        M = [Ix2(i,j) Ixy(i,j);Ixy(i,j) Iy2(i,j)];             % auto correlation matrix        R(i,j) = det(M)-0.06*(trace(M))^2;                     % 计算R        if R(i,j) > Rmax            Rmax = R(i,j);        end;    end;end;cnt = 0;for i = 2:height-1    for j = 2:width-1        % 进行非极大抑制，窗口大小3*3        if R(i,j) > 0.01*Rmax && R(i,j) > R(i-1,j-1) && R(i,j) > R(i-1,j) && R(i,j) > R(i-1,j+1) && R(i,j) > R(i,j-1) && R(i,j) > R(i,j+1) && R(i,j) > R(i+1,j-1) && R(i,j) > R(i+1,j) && R(i,j) > R(i+1,j+1)            result(i,j) = 1;            cnt = cnt+1;        end;    end;end;i=1;    for j=1:height        for k=1:width            if result(j,k)==1;                corners1(i,1)=j;                corners1(i,2)=k;                i=i+1;            end;        end;    end;[posc, posr] = find(result == 1);cnt                                      % 角点个数imshow(ori_im)hold on;plot(posr,posc,'r+');a=ginput(1);b=ginput(1);j=1;for i=1:cnt    if corners1(i,1)>a(1,2) && corners1(i,1)<b(1,2)        if corners1(i,2)>a(1,1) && corners1(i,2)<b(1,1)            B(j,1)=corners1(i,1);            B(j,2)=corners1(i,2);            j=j+1;        end    endendxlswrite('C:\Documents and Settings\ipsuser\桌面\Harris\ceshidata.xls',B,'Sheet1','C1');</span>