哈工大遥感图像配准报告

HarbinInstituteofTechnology多时相遥感图像配准实验报告课程名称：院系：姓名：学号：授课教师：哈尔滨工业大学1.实验目的（1）了解图像配准原理、导致几何失真的原因及矫正方法（2）掌握两幅图像配准的基本方法（3）学会评价配准结果2.实验原理图像配准就是将不同时间、不同传感器或不同条件下（天候、照度、摄像位置和角度等）获取的两幅或多幅图像进行匹配、叠加的过程。图像配准的基本过程是寻找把待配准图像中的点映射到参考图像中对应点的最佳变换，从而去除待配准图像和参考图像之间在几何上的不一致，如平移、旋转和形变，最终使得同一目标在不同的图像上具有相同的坐标位置。图像配准中的主要问题就是确定这种映射关系，一般来讲我们将需要进行变换的图像成为输入图像或待配准图像，将用于矫正输入输入图像的图像成为参考图像。图像配准的基本步骤包括：（1）特征提取：在进行图像配准之前根据实际需要选择图像中待匹配特征。（2）特征匹配：在确定了来自不同图像的特征集后，建立特征的对应关系。对应关系建立可以通过使用各种特征描述、相似性测量以及空间关联的方法实现。（3）变换模型估计：通过估计或者最优搜索等方法得到映射函数的类型和参数，以便能够对齐待配准图像和基准图像。（4）图像重采样和变换：依靠映射函数对输入图像进行图像变换，使用适当的插值算法计算出非整数点坐标的图像值。常用的重采样方法包括：（1）最近邻重采样法：简单地选择其中心距图像上的点最近的实际像素，然后将这个像素变换到相应的显示网格上。这种技术适用于新图像用于分类的情况，原因是该图像是由原始像素的亮度构成的，只是按照正确的几何关系将像素进行重新排列。（2）双线性插值法：对于图像中与给定的显示网格位置对应点的周边四个像素使用三次线性插值。（3）三次卷积插值法：利用周边16个像素，在图像中一点四周各取4个像素，利用三次多项式对这些像素确定的四条线对进行拟合以形成4个插值点，然后利用第5个三次多项式对这四个插值点进行拟合，最终合成在显示网格对应位置的亮度值。利用抽样定理以及由一个离散样本集构建连续函数的方法，可以推导出用于插值的多项式的实用形式。3.实验步骤（1）实验图像为：参考图像待配准图像（2）选取控制点，一共手动选取了20个控制点对：（3）双线性映射模型为：12345678xcvcwcvwcycvcwcvwc利用最小二乘法计算模型参数得：123456780.9995,0.0034,0.0000,58.90260.0008,0.9982,0.0000,53.5062cccccccc2阶多项式映射模型为:2212345622789101112xcvcwcvwcvcwcycvcwcvwcvcwc利用最小二乘法计算模型参数得：1234567891011120.9950,0.0093,1.6155,7.4176,1.36025,57.96760.0045,1.0062,1.98196,8.98226,1.77005,52.3716cccecececcccececec（4）获得配准结果：双线性映射模型：配准后图像参考图像二阶多项式映射模型：配准后图像参考图像4.实验结果与分析选取10对同名点进行定量分析：双线性映射模型：x方向标准差y方向标准差综合标准差0.60200.14000.6180二阶多项式映射模型：x方向标准差y方向标准差综合标准差0.62700.33520.7110由分析结果可以看出两种模型的处理效果并无明显区别，但在一些边缘部位如图中黄色方框所示，二阶多项式映射模型先比喻会产生一些杂散点。从计算得的标准差来看，双线性模型优于二阶多项式模型，这是因为本次实验中两幅图片的差异主要由于几何平移与旋转造成，不存在扭曲等变化，使用双线性模型更合适，可以保证全局达到最优，而二阶多项式模型只能保证局部参考点最优。5.实验程序clearall;closeall;image1=imread('E:\遥感导论\实验一\subcut1.tif');%参考图像image2=imread('E:\遥感导论\实验一\subcut2.tif');%输入图像load('controlpoints.mat')load('testpoints.mat')%双线性模型n=size(fixedPoints_control,1);emin=1;sume=0;RMS=0;jianyan_point=[];jbase_point=[];C1=[];C2=[];fori=1:njianyan_point=fixedPoints_control(i,:);fixedPoints_control(i,:)=[];jbase_point=movingPoints_control(i,:);movingPoints_control(i,:)=[];d=ones(n-1,4);d(:,1)=fixedPoints_control(:,1);d(:,2)=fixedPoints_control(:,2);d(:,3)=fixedPoints_control(:,1).*fixedPoints_control(:,2);X=movingPoints_control(:,1);Y=movingPoints_control(:,2);C1(:,i)=(inv(d'*d))*d'*X;C2(:,i)=(inv(d'*d))*d'*Y;I=[jianyan_point(1,1),jianyan_point(1,2),jianyan_point(1,1)*jianyan_point(1,2),1];X_r=C1(:,i)'*I';Y_r=C2(:,i)'*I';e(i)=(jbase_point(1,1)-X_r)^2+(jbase_point(1,2)-Y_r)^2;ife(i)eminC1_best=C1(:,i);C2_best=C2(:,i);emin=e(i);endifi1;fixedPoints_control=[fixedPoints_control(1:i-1,:);jianyan_point;fixedPoints_control(i:n-1,:)];movingPoints_control=[movingPoints_control(1:i-1,:);jbase_point;movingPoints_control(i:n-1,:)];elsefixedPoints_control=[jianyan_point;fixedPoints_control(i:n-1,:)];movingPoints_control=[jbase_point;movingPoints_control(i:n-1,:)];endsume=sume+e(i);endRMS=sqrt(sume/n)eminC1_bestC2_bestimage1_1=image1(:,:,1);image2_1=image2(:,:,1);image3=zeros(512,512);fori0=1:512;forj0=1:512;I=[i0,j0,i0*j0,1];i=round(C1_best'*I');j=round(C2_best'*I');if(i=0||j=0||i512||j512)image3(j0,i0)=image1_1(j0,i0);elseimage3(j0,i0)=image2_1(j,i);endendendfigure(1)imshow(image3,[])title('配准后图像')figure(2)imshow(image1(:,:,1),[])title('参考图像')figure(3)imshow(image2_1,[])title('待配准图像')ex=[];ey=[];sumx=0;sumy=0;sum=0;fori=1:10;S=[fixedPoints_test(i,1),fixedPoints_test(i,2),fixedPoints_test(i,1)*fixedPoints_test(i,2),1];M(i)=C1_best'*S';N(i)=C2_best'*S';ex(i)=movingPoints_test(i,1)-M(i);ey(i)=movingPoints_test(i,2)-N(i);sumx=sumx+ex(i)^2;sumy=sumy+ey(i)^2;sum=sum+ex(i)^2+ey(i)^2;endmsex=sqrt(sumx/10)msey=sqrt(sumy/10)mse=sqrt(sum/10)clearall;closeall;image1=imread('E:\遥感导论\实验一\subcut1.tif');%参考图像image2=imread('E:\遥感导论\实验一\subcut2.tif');%输入图像%figure,imshow(image1,[])%title('参考图像'),axison%figure,imshow(image2,[])%title('待配准图像'),axisoncpselect(image2,image1)load('controlpoints.mat')load('testpoints.mat')%二阶多项式模型n=size(fixedPoints_control,1);emin=1;sume=0;RMS=0;jianyan_point=[];jbase_point=[];C1=[];C2=[];fori=1:njianyan_point=fixedPoints_control(i,:);fixedPoints_control(i,:)=[];jbase_point=movingPoints_control(i,:);movingPoints_control(i,:)=[];d=ones(n-1,6);d(:,1)=fixedPoints_control(:,1);d(:,2)=fixedPoints_control(:,2);d(:,3)=fixedPoints_control(:,1).*fixedPoints_control(:,2);d(:,4)=fixedPoints_control(:,1).*fixedPoints_control(:,1);d(:,5)=fixedPoints_control(:,2).*fixedPoints_control(:,2);X=movingPoints_control(:,1);Y=movingPoints_control(:,2);C1(:,i)=(inv(d'*d))*d'*X;C2(:,i)=(inv(d'*d))*d'*Y;I=[jianyan_point(1,1),jianyan_point(1,2),jianyan_point(1,1)*jianyan_point(1,2),jianyan_point(1,1)*jianyan_point(1,1),...jianyan_point(1,2)*jianyan_point(1,2),1];X_r=C1(:,i)'*I';Y_r=C2(:,i)'*I';e(i)=(jbase_point(1,1)-X_r)^2+(jbase_point(1,2)-Y_r)^2;ife(i)eminC1_best=C1(:,i);C2_best=C2(:,i);emin=e(i);endifi1;fixedPoints_control=[fixedPoints_control(1:i-1,:);jianyan_point;fixedPoints_control(i:n-1,:)];movingPoints_control=[movingPoints_control(1:i-1,:);jbase_point;movingPoints_control