harris角点检测与ncc匹配

转自：http://zixuanjinan.blog.163.com/blog/static/11543032620097510122831/

file1:--------------------------------------------------------------------------------------

function [y1,y2,r,c]=harris(X)
% 角点的检测，利用harris 算法
% 输出的是一幅图像
% [result,cnt,r,c]=harris(X)

% clc,clear all;
% filename='qiao1.bmp';
% X= imread('filename.bmp'); % 读取图像
% Info=imfinfo(filename);

% f=rgb2gray(X);
f=X;
%
ori_im=double(f)/255;%unit8转化为64为双精度double64
fx = [-2 -1 0 1 2]; % x方向梯度算子(用于Harris角点提取算法)
Ix = filter2(fx,ori_im); % x方向滤波 善于使用filter
% 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',[10 10 ],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)]; %2*2的矩阵
R(i,j) = det(M)-0.06*(trace(M))^2;% 计算R 求得RMAX，看来是整体求得的角点响应函数
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;

[posr, posc] = find(result == 1); % 角点x、y坐标
% imshow(ori_im) %和 X的效果是一样的
% hold on;
% plot(posr,posc,'r.');
y1=result;
y2=cnt;
r=posr;c=posc;
return;

file2----------------------------------------------------------------------------------------------------

function res=match(a1,cnt1,r1,c1,a2,cnt2,r2,c2)
% res=match(a1,a2)
% 将从a1寻找a2中的最佳匹配点，得到从a2中抽取的res，也就是单向搜索
% [result1,cnt1,r11,c11]=harris(a1);
% [result2,cnt2,r22,c22]=harris(a2);%可以保证想匹配哪些点就匹配哪些
% figure;
% imshow(result1);title('result1角点位置');
% figure;title('result2角点位置');
% imshow(result2);

win=[1/9 1/9 1/9;1/9 1/9 1/9;1/9 1/9 1/9];
u1=filter2(win,a1);
u2=filter2(win,a2); %求均值

%基于点特征的图像配准算法研究 山大

a1=double(a1);
a2=double(a2);
A=filter2(win,(a1-u1).^2);%求方差
B=filter2(win,(a2-u2).^2);
[m1,n1]=size(a1);
[m2,n2]=size(a2);
res1=zeros(m1,n1);
res2=zeros(m2,n2); %寻找的匹配的点

for s=1:cnt1
max=0; p=0;q=0;i=r1(s,1);j=c1(s,1); %p.q存放坐标
for v=1:cnt2
m=r2(v,1);n=c2(v,1);
k1=(a1(i-1,j-1)-u1(i,j))*(a2(m-1,n-1)-u2(m,n)); %用result是找不到什么信息，3*3对于如此稀疏点没有用
k2=(a1(i-1,j)-u1(i,j))*(a2(m-1,n)-u2(m,n));
k3=(a1(i-1,j+1)-u1(i,j))*(a2(m-1,n+1)-u2(m,n)); %用循环要好一些
k4=(a1(i,j-1)-u1(i,j))*(a2(m,n-1)-u2(m,n));
k5=(a1(i,j)-u1(i,j))*(a2(m,n)-u2(m,n));
k6=(a1(i,j+1)-u1(i,j))*(a2(m,n+1)-u2(m,n));
k7=(a1(i+1,j-1)-u1(i,j))*(a2(m+1,n-1)-u2(m,n));
k8=(a1(i+1,j)-u1(i,j))*(a2(m+1,n)-u2(m,n));
k9=(a1(i+1,j+1)-u1(i,j))*(a2(m+1,n+1)-u2(m,n));
num=k1+k2+k3+k4+k5+k6+k7+k8+k9;

den=sqrt(A(i,j)*B(m,n));
ncc=num/den;
if ncc>max
max=ncc;p=m;q=n;
end

end

res2(p,q)=1;

end
%做不出来可以先搜索一个点
res=res2;
return

file3----------------------------------------------------------------------------------------------------
%特征点的匹配，主要应用 harris 角点的检测,match单向匹配 函数
%适合于有白边的图像，因为，在加窗滤波的时候，没有限定范围，尽量保证角点不在边上
clc,clear all;
a1=imread('qiao1.bmp');
a2=imread('qiao2.bmp'); %double的作用很大
subplot(1,2,1); imshow(a1);subplot(1,2,2);imshow(a2); title('原来的图像');

[result1,cnt1,r1,c1]=harris(a1);%角点检测，得到原始的焦点位置图result
[result2,cnt2,r2,c2]=harris(a2);
figure;subplot(1,2,1);imshow(a1);hold on;plot(c1,r1,'g.');
subplot(1,2,2);imshow(a2);hold on;plot(c2,r2,'g.'); title('图1，2的角点图');

res2=match(a1,cnt1,r1,c1,a2,cnt2,r2,c2);%从result1中开始搜索在result2中可能达到的
[r22,c22]=find(res2==1);
[m22,n22]=size(r22);
cnt22=m22;

res1=match(a2,cnt22,r22,c22,a1,cnt1,r1,c1);%反向搜索res2--result1
res1=and(res1,result1); %保证反向匹配不会出现不可能的点
[r11,c11]=find(res1==1);
[m11,n11]=size(r11);
cnt11=m11;

res22=match(a1,cnt11,r11,c11,a2,cnt22,r22,c22);%从res1中开始搜索在res2中可能达到的
res22=and(res22,res2);
[r222,c222]=find(res22==1);
[m222,n222]=size(r222);
cnt222=m222;
%当然，匹配的次数越多，点数会越来越少,这个匹配永远都只是相对的
%注意最终结果的点数不一定一样的哈，也没有必要一一对应，几个点就可以决定参数了的
figure;subplot(1,2,1);imshow(a1);hold on;plot(c11,r11,'g.');
subplot(1,2,2);imshow(a2);hold on;plot(c222,r222,'g.'); %单向搜索和双向搜索到的结果
title('最终得到的匹配的角点');