遥感图像配准matlab

首先进行Harris角点特征提取，而后利用NCC算法进行粗匹配，剔除误匹配和不匹配向量，基于灰度相关系数，计算配准误差， 从而得到配准后叠加图像。适用于输入图像存在平移变换和旋转变换；能够实现可见光区图像配准，可实现时间间隔不是很久的多时相遥感影像配准

function [v1_gcpu1_gcpv2_gcpu2_gcp]=GCP_Select2（v1u1v2u2I1I2）
[height1width1]=size（I1）;
[height2width2]=size（I2）;
%采用NCC（归一化互相关）算法对两幅图像的角点进行粗匹配%
k=15;                                                                      %邻域窗口大小（2k+1）*（2k+1））
p=1;
R=zeros（1length（u2）） * （1e-10）;
[v1_new u1_newv2_newu2_new] = deal（zeros（1length（v1）））;

for i=1:length（v1）
    if（（（u1（i）+k）0）&&（（v1（i）-k）>0））
        G1_window = I1（ v1（i）-k:v1（i）+k u1（i）-k:u1（i）+k ）;
        G1_avg = mean（G1_window（:））;
        for j=1:length（v2）
            G_sum=0;G1_sub=0;G2_sub=0;
            if（（（u2（j）+k）0）&&（（v2（j）-k）>0））&&（abs（v1（i）-v2（j））<80）&&（abs（u1（i）-u2（j））<120）
                G2_window = I2（ v2（j）-k:v2（j）+k u2（j）-k:u2（j）+k ）;     
                G2_avg = mean（G2_window（:））;
                
                G_sum = sum（（G1_window（:） - G1_avg）.*（G2_window（:） - G2_avg））;
                G1_sub = sum（（G1_window（:） - G1_avg）.^2）;
                G2_sub = sum（（G2_window（:） - G2_avg）.^2）;
            end
            R（j） = G_sum / （sqrt（G1_sub*G2_sub） + 1e-10）;
        end
        R_max=max（R）;
        if （R_max>0.6）                                                   %阈值0.5，可调，视需要粗配准点数而定
            v1_new（p）=v1（i）;
            u1_new（p）=u1（i）;
            [~y]=find（R==R_max）;                                        %构建新的角点位置向量，该向量中I1I2角点一一对应
            v2_new（p）=v2（y）;
            u2_new（p）=u2（y）;
            p=p+1;
        end
    end
end

p = p-1;
[v1_new u1_newv2_newu2_new] = deal（v1_new（1:p） u1_new（1:p）v2_new（1:p）u2_new（1:p））;

%通过简单方法对匹配点对进行筛选，剔除错误匹配点对（剔除距离出现次数最多之外的点），统计两图像坐标差值最多的坐标为新的匹配点集%
x_sub=zeros（1length（u1_new））;
y_sub=zeros（1length（u1_new））;
for i=1:length（u1_new）
    x_sub（i）=u1_new（i）-u2_new（i）;
    y_sub（i）=v1_new（i）-v2_new（i）;
end
x=unique（x_sub）;
[mn]=hist（x_subx）;
x_sub_most=n（find（m==max（m）））;
[mn]=hist（y_suby）;
y_sub_most=n（find（m==max（m）））;
p=1;
[v1_new2 u1_new2v2_new2u2_new2] = deal（zeros（1length（u1_new）））;%剔除错误匹配点对之后的新点对起始坐标存储位置
for i=1:length（u1_new）
    if（ （abs（x_sub（i）-x_sub_most（1））<1）&&（abs（y_sub（i）-y_sub_most（1））<1））
        v1_new2（p）=v1_new（i）;u1_new2（p）=u1_new（i）;
        v2_new2（p）=v2_new（i）;u2_new2（p）=u2_new（i）;
        p=p+1;
    end
end
p = p-1;
[v1_new2 u1_new2v2_new2u2_new2] = deal（v1_new2（1:p） u1_new2（1:p）v2_new2（1:p）u2_new2（1:p））;
[v1_gcp u1_gcp v2_gcp u2_gcp] = deal（v1_new2 u1_new2v2_new2u2_new2）;

% % 将I1作为参考图像，在I2中的对应角点附近寻找I1角点的同名点 %
[ghz] = deal（1）;  %搜索窗口大小为2*k+1  %灰度窗口大小为2*h+1    %重新筛选之后的地面控制点起始存储位置
R_G1G2=zeros（33）;
for i=1:length（u1_new2）
    G1_gary = I1（v1_new2（i）-h:v1_new2（i）+h u1_new2（i）-h:u1_new2（i）+h）;
    G1_avg = mean（G1_gary（:））;                                                                %生成I1中灰度窗口中的灰度向量
    for k=-g:1:g
        for l=-g:1:g
            G2_gary = I2（v2_new2（i）+k-h:v2_new2（i）+k+h u2_new2（i）+l-h:

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----

     文件       4035  2019-05-29 20:39  registration\GCP_Select2.m

     文件       2107  2018-01-12 11:01  registration\Harris.m

     文件     786624  2018-01-12 11:01  registration\I1.tif

     文件     786624  2018-01-12 11:01  registration\I2.tif

     文件     184572  2018-01-12 11:01  registration\I3.tif

     文件     186756  2018-01-12 11:01  registration\I4.tif

     文件      89584  2018-01-12 11:01  registration\I5.jpg

     文件      86935  2018-01-12 11:01  registration\I6.jpg

     文件       2452  2019-05-29 19:59  registration\ImgRegistration.m

     文件       3296  2018-01-12 11:01  registration\Registration.m

     文件       1758  2019-05-29 15:23  registration\registration_manual.m

     文件        673  2018-09-07 14:17  registration\数据\1\1.txt

     文件     254517  2018-09-24 21:07  registration\数据\1\Sar（从）.png

     文件     233601  2018-09-24 19:46  registration\数据\1\ZY（主）.png

     文件        679  2018-09-07 14:25  registration\数据\2\2.txt

     文件     190726  2018-05-18 16:12  registration\数据\2\主.png

     文件     258692  2018-09-26 10:13  registration\数据\2\从.png

     文件       1281  2018-05-18 20:20  registration\数据\3\lt5tosar.txt

     文件      40831  2018-09-25 20:12  registration\数据\3\主.png

     文件      60592  2018-05-18 18:14  registration\数据\3\从.png

     文件      63481  2018-05-18 20:38  registration\数据\4\dataenvisat.png

     文件      47437  2018-09-25 20:09  registration\数据\4\datalt5.png

     文件       1298  2018-05-18 21:29  registration\数据\4\lt5tosar.txt

     文件        390  2018-09-07 14:29  registration\数据\5\5.txt

     文件      36428  2018-09-07 09:07  registration\数据\5\主.png

     文件      83199  2018-09-07 09:07  registration\数据\5\从.png

     目录          0  2019-09-05 08:24  registration\数据\1

     目录          0  2019-09-05 08:24  registration\数据\2

     目录          0  2019-09-05 08:24  registration\数据\3

     目录          0  2019-09-05 08:24  registration\数据\4

............此处省略6个文件信息