harries算子特征提取与图像匹配

概述：图像处理过程中经常要遇到图像匹配的问题，通过灰度相关法可以解决比较简单的平移、旋转、亮度差异匹配的问题，当平移角度过大，灰度差异过大、旋转角度较大时，灰度相关法会变得比较麻烦，采用角点匹配的方法能取得不错的成果，高级的有sift和suft，这里我们介绍一种比较简单的角点检测方法harris
，对harris原理上的东西不多讲，主要谈一下图像匹配注意点，以及我自己的处理方式，方式是以前学习时自己想的，没有查阅相关资料，仅供参考。
一、Harris角点检测
注意点1，非极大值抑制，采用阈值分割的方法得到的角点往往都扎堆存在，进行非极大值抑制，得到最合适的角点，如图。

注意点2，边界处理，由于原始图像可能是经过拼接后的图像，形状不规则，其边缘可能不是矩形，在边界上非常容易产生角点，这些角点都要除去，如下图绿色圈内的点

注意点3，但是由于先验知识，我们可以知道图像的相对位置关系，比如一张在左一张在右 ，在选取角点 时，对于左图我们可以只保留靠近右边的角点。右下图为删除靠左角点后的，这两张图是在第12展图匹配拼接后求的Harris角点用来和第三张图匹配

二、Harris角点粗匹配
得到两张图的harris角点特征后，我们便可以开始粗匹配了。
粗匹配时其实还是采用的灰度相关加上双向最大原则的匹配方法，即取每个角点周围的一个邻域，然后求图一中角点和图二中角点的互相关函数。互为最大值则视为两点匹配成功。
这里要特别关心在边界周围的点，由于matlab中第一次拼接后的图像中含有nan（表示没有数据，但是运算时又会看成max）因此在角点检测的时候建议把邻域内含nan的点全部去掉。

三、基于randsac方法的Harris角点精确匹配
在精匹配前我们需要选定变换的模型，这里我们采用仿射变换，关于仿射变换，不多介绍。仿射变换可以用下面这个式子表示
P1=F*P2 F是仿射矩阵，P2是原始坐标，P1是变换后坐标
randsac方法也不多介绍，原理很简单，随机选取满足模型的最少匹配点（仿射模型时为3对），进行逆变换得到仿射矩阵F，然后将所有P2变换得到P1’，比较P1与P1’，记录误差小于阈值的点的个数，点数最多的那组点解出来的F便是最优模型。
得到的精匹配图如下
四、反变换更改图像坐标系
这里其实就是知道模型，进行仿射变换的过程。
注意，要将两张图调整到统一坐标系下，因为仿射变换的过程中可能会产生负的坐标值，因此在一开始处理的时候，我就将基准图，向上平移了，消除坐标中的幅值。
变换后图片如下，可以看出其实前两张图片主要是经过了平移变换

五、图像拼接
图像拼接的方式有很多，可以采用距离重叠中心的方式进行加权。
第一步，将两张图调整为一样的大小，注意调整大小时，要保持坐标不变。
第二步，计算每一行中重叠元素个数，
第三步，计算每一行中第一个重叠的位置
第四步，生成权值矩阵
第五步，加权求和。
拼接结果如图

六、附件：
原始图片三张

结果
程序代码

clc
clear
close all

%已知条件一，由于在统一位置拍摄，因此图像变换后面积变化不大
path1=imgetfile;
path2=imgetfile;
path3=imgetfile;
imrgb1=im2double(imread(path1));
imrgb1(size(imrgb1,1)+50,size(imrgb1,2))=nan;%上拉50个点，在仿射变换中，防止先到原图坐标系有负数
harrispipei(imrgb1,im2double(imread(path2)),5,1);%打包好的函数实现,不足的是点是自己找的
harrispipei(ans,im2double(imread(path3)),4,2);%打包好的函数实现,不足的是点是自己找的


function imo=harrispipei(imrgb1,imrgb2,r,cnt)
    [imrgb1,imrgb2]=im2samesizeup(imrgb1,imrgb2,nan);
    img1=(rgb2gray(imrgb1));
    img2=(rgb2gray(imrgb2));
    img1(find(img1==0))=nan;
    img2(find(img2==0))=nan;
    xyh1=HarrisPoint(img1,0.8,r,1,cnt-1);%  
    xyh2=HarrisPoint(img2,0.8,r,2,2*cnt-2);
    [xya1,xya2]=HarrisPointMosaic(img1,img2,xyh1,xyh2,r);%得到粗匹配点
    [A,xy1,xy2]=randsac([xya1,xya2],img1,img2)
    
%     if cnt==1
%         xy1=[147,158;139 193;216 296;187 207];
%         xy2=[253,163;246 199;320 296;289 212]-[224 0]+[0 50];
%     end
%     if cnt==2
%         xy1=[273 133;320 190;273 235];
%         xy2=[399 183;434 239;392 285]-[374 0];
%     end
%     if cnt==3
%         xy1=[469 245;356 422;370 490];
%         xy2=[683 220;576 397;589 471]-[500,0];
%     end
    A=rfangshe(xy2,xy1);
    [im2to1 dx1 dy1]=fangshe(img2,A);
    HarrisPoint(im2to1,0.8,r,3,0);
    plotcpl(img1,img2,xya1,xya2);
    plotcpl(img1,img2,xy1,xy2);
    %拼接
    imrgb21=fangshe(imrgb2,A);
    [imp1,imp2]=im2samesizeup(imrgb1,imrgb21,nan);%调整到同一个大小，这个是在图像上方补nan
    imo=pj12(imp1,imp2);%拼接在一起
    figure; imshow(imo);
end
%辅助函数
function plot_rxm_wh(img,pos,fig)%画红叉,pos为[w,h],支持多个图画到一个里面
    if fig==1
        figure
    end
    subplot(1,3,fig);
    imshow(img);hold on;
    for k=1:size(pos,1)
        plot(pos(k,2),pos(k,1),'r+');
    end
    hold off;
end
function [x y]=exchange(x,y)
    temp=x;x=y;y=temp;
end
function hwo=dropoutedge(hwi,mat2,r);%去掉靠近边界的点
    [hm,wm]=size(mat2);
    h=find(hwi(:,1)>r);    hwi=hwi(h,:);%删掉不满足的点
    h=find(hwi(:,1)<hm-r);    hwi=hwi(h,:);%删掉不满足的点
    h=find(hwi(:,2)>r);    hwi=hwi(h,:);%删掉不满足的点
    h=find(hwi(:,2)<wm-r); 
    hwo=hwi(h,:);%删掉不满足的点
end
function [imo1,imo2]=im2samesizeup(im1,im2,num)%将两张图调整到一样大小
    [h1,w1,d1]=size(im1),[h2,w2,d2]=size(im2);
    hm=max(h1,h2);
    wm=max(w1,w2);
    if h1~=hm||w1~=wm
        temp=ones(hm,wm,d1)*num;
        im1=flipud(im1);
        temp(1:h1,1:w1,1:d1)=im1;
        im1=flipud(temp);
    end
    if h2~=hm||w2~=wm
        temp=ones(hm,wm,d2)*num;
        im2=flipud(im2);
        temp(1:h2,1:w2,1:d2)=im2;
        im2=flipud(temp);
    end
    imo1=im1;
    imo2=im2;
end
function [imo1,imo2]=im2samesizedown(im1,im2,num)%将两张图调整到一样大小
    [h1,w1,~]