轉自：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('最終得到的匹配的角點');