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

（1-4）

<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
    end
end</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为双精度double64
fx = [-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=2

Ix2 = 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);           % 纪录角点位置，角点处值为1

R = 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
    end
end
xlswrite('C:\Documents and Settings\ipsuser\桌面\Harris\ceshidata.xls',B,'Sheet1','C1');</span>