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('最终得到的匹配的角点');