1、1课程编号：20162024123 课程性质：专业必修数字摄影测量实习实习报告学院： 测绘学院 专业： 测绘工程 地点： 校 内 班级： 组号： 姓名： 学号： 教师： 2017年 6 月 19 日 至 2017 年 7 月 8 日2第 1 章 空三报告 .41.1 NAT 空中三角成果数据流程 41.2 空三加密需要注意 41.3 空三流程图 41.4 空三评查步骤 51.4.1 需要平差的原因 .51.4.2 粗差的影响 .51.4.3 平差方法-循环粗差剔除 51.4.4 自己的方案 .5第 2 章 Mapmatrix 进行 DOM,DEM 生成 62.1 外方位元素的确定 62.1.1
2、 内定向：确定与物点相对应的像点坐标 .62.1.2 扫描分辨率为什么不能错： .62.1.3 创建立体像对： .62.1.4 相对定向： .62.1.5 绝对定向： .6借助物空间坐标为已知的控制点来确定空间辅助坐标系与实际物空间坐标的变换关系称为立体模型的绝对定向。 .62.1.6 后方交会 .62.1.7 后方交会绝对定向刺 6 个点能成功： .62.1.8 核线、匹配采样 .72.1.9 影像匹配 .72.2 DEM 编辑生成 .72.2.1 选取新建的 DEM 成品自动生成 72.2.2 DEM 编辑 72.2.3 DEM 介绍 72.2.4 问老师得到的对 Dem 编辑的理解 .8
3、2.2.5 Dem 拼接裁剪 .82.3 DOM 的生成编辑 .832.3.1 DOM 概念 82.3.2 选取新建的 DOM 产品自动生成 8调整扫描分辨率与参数 .82.3.3 规格要求 .82.3.4 生成 72 张影像的 DOM 导入易拼图进行正射影像编辑 92.3.5 EPT 编辑工具使用总结 92.3.6 DOM 建筑等人工地物变形对比分析 102.3.7 DOM 色彩对比分析 102.3.8 最后生成的产品 .10（经过房屋过滤后的 DEM） .10第 3 章 DLG 的立测生成 .113.1 DLG 的概念 113.2 本次试验生成的 DLG 成果，导出为 DXF 格式 113
4、.3 立测的精度检核 123.4 DLG 生成 123.4.1 经验总结 .12第 4 章 编程实习 .134.1 实验目的 134.2 具体要求 134.3 计算原理 134.4 算法解释 154.5 算法流程： 154.6 MATLAB 源代码 .164.7 影像匹配原理 19常用的影像匹配方法可分为两种：基于影像灰度匹配的影像匹配算法和基于影像特征的匹配算法。 .194.8 主程序 204.9 函数文件 254第 1 章 空三报告1.1 NAT 空中三角成果数据流程1. 设置测区基本参数2. 建立相机文件，导入老师给的文件3. 导入外业控制点4. 影像格式转换5. 自动转点6. 进行 P
5、ATB 平差计算1.2 空三加密需要注意1. 像素大小、焦距、像主点坐标、畸变参数，像幅尺寸 (数码相机)2. 模型连接点平面及高程限差、地面定向点残差限差、地面检查点残差限差3. 控制点坐标、刺点片及点之记1.3 空三流程图51.4 空三评查步骤1.4.1 需要平差的原因1 同一物体在不同视角下的图像有极大的不同，物体 遮挡而丢失信息2 光照、场景中物体的几何形状和物理特性、物体和摄像机之间的空间关系3 从带畸变和噪声的、二维的、单一的灰度值中提取出尽可能不带畸变和噪声的、三维场景中的诸多因素1.4.2 粗差的影响1 观测值中存在粗差，则某一观测的误差对所有残差都有影响2 一个观测的残差包含
6、着所有粗差的影响1.4.3 平差方法-循环粗差剔除1 先外围加 4 个点进行平差解算像方的粗差点排除到一定程度后，开始加预测的控制点，不要一次性刺入所有预测控制点，建议从外向内回字形一圈圈刺预测控制点，留区域中的 3 到 4 个控制点做检核点。每刺入一轮都先平差，编辑得差不多再进行下一轮。2 继续调整控制点关掉粗差探测，勾选验后方差，再平差一次，平差结束。3 上述为精度报告检查平差结果，最后需要到立体测图软件进行检查。1.4.4 自己的方案在本次测区的 15 个点中，我分了三轮进行平差，每次选取 4 个点最后保留3 个点作为检查点，难处在于有的点有灰尘遮盖造成精度不准，我的精度指标如下所示。6
7、第 2 章 Mapmatrix 进行 DOM,DEM 生成2.1 外方位元素的确定2.1.1 内定向：确定与物点相对应的像点坐标本次实习测量了四个框标，故采用仿射变换形式X=a0+a1*x+a2*y;Y=b0+b1*y+b2*y;2.1.2 扫描分辨率为什么不能错：量测的像点坐标存在着从扫描坐标岛像点坐标的转换，如果没有调整合适的分辨率就会对影像内定向产生影响。2.1.3 创建立体像对：2.1.4 相对定向：1.单独像对相对定向2.连续像对相对定向目的是为了恢复摄影时两张影像光束的相互关系，从而使同名光线对对相交。5 个元素2.1.5 绝对定向：借助物空间坐标为已知的控制点来确定空间辅助坐标系
8、与实际物空间坐标的变换关系称为立体模型的绝对定向。2.1.6 后方交会后方交会：如果知道每张影像的 6 个外方位元素就可以知道确定被摄物体与航射影像的关系，利用影像范围内的控制点的空间坐标与影像坐标根据共线方程式，反求影像的外方位元素。2.1.7 后方交会绝对定向刺 6 个点能成功：绝对定向有七个参数，空间相似变换的 3 个角元素，3 个平移量，1 个比例尺的缩放。2 个平高+1 个高程控制点。独立模型法区域空中三角测量的基本思想，把一个单元模型视为刚体，一个7单元模型由两个或更多个立体像对构成，因为一个立体像对需要两张影像，所以如图所示，两张影像的拼接必须要有三个同名像点作为公共点，满足 X
9、,Y,Z 的传递。连接过程中每个模型可以视为刚体，只能做平移，旋转，缩放。两张相片外方位元素 12 个，绝对定向的 3 个角元素，3 个平移量，1 个比例尺的缩放，相对定向 5 个元素。通过相对定向绝对定向能够恢复两张影像的外方位元素。2.1.8 核线、匹配采样核线的性质：同名像点必位于同名核线上。通过摄影基线与地面所作的平面称为核面核面与影像面交线称为核线同名像点 必定在同名核线上。点击定义核线范围，重采样即可完成。2.1.9 影像匹配2.2 DEM 编辑生成2.2.1 选取新建的 DEM 成品自动生成2.2.2 DEM 编辑（对 hammer 编辑）2.2.3 DEM 介绍Hammer 第
10、一条航带的 DEM 编辑，数字高程模型（Digital Elevation Model)，简称 DEM，是通过有限的地形高程数据实现对地面地形的数字化模拟（即地形表面形态的数字化表达） ，它是用一组有序数值阵列形式表示地面高程的一种实体地面模型。82.2.4 问老师得到的对 Dem 编辑的理解第一条航带编辑必须要把浮动到眼前，通过 3D 眼睛看到的明显具有高程的房屋建筑改化到地面高程中去，生成的是地物地貌模型，而不是具有人工建筑的地物模型。将房屋等人工建筑周围的地形高程进行内插，求出该店的实际高程。方法：匹配的点，没贴准地面需要编辑内插的点没贴准地面需要编辑匹配点内插方式：不量测地面的点，用两
11、边的点对中间点内插量测点内插方式：量测地面的点2.2.5 Dem 拼接裁剪（Dem 成图）2.3 DOM 的生成编辑2.3.1 DOM 概念DOM-数字正射影像图是对航空（或航天）相片进行数 字微分纠正和镶嵌，按一定图幅范围裁剪生成的数字正射影 像集。它是同时具有地图几何精度和影像特征的图像。2.3.2 选取新建的 DOM 产品自动生成调整扫描分辨率与参数 沿影像边缘生成，且原始影像单独生成 DOM2.3.3 规格要求 数据能正确衔接9 建筑等无错开、扭曲变形等现象 色调一致、合理1 YLD 裁切后的影像（0.1 米） （用 DSM 生成的拼接后的 DOM 不用单独提交，主要用于 DOM 建筑
12、等人工地物的变形分析）2 hammer 数据拼接后 DOM 0.5 米（不用编辑，主要用于 DOM 色彩分析） 2.3.4 生成 72 张影像的 DOM 导入易拼图进行正射影像编辑（影像镶嵌）（待编辑的正射影像）2.3.5 EPT 编辑工具使用总结将 72 张影像生成的 DOM 拼接导入 EPT 进行编辑编辑注意要不过水体，不过房屋，可以借助工具进行镶嵌线编辑，添加点、插入点、移动点、删除点、删10除临时线和修补选区。注意事项：对于某些四点公用点的镶嵌线编辑1.在与关键点相连的任意一条镶嵌线上单击左键添加点，开始绘制。2.到第四条镶嵌线上，单击左键添加结束点。3.最后单击右键，结束镶嵌线编辑。
13、4.其中绘制镶嵌线时，开始点和结束点必须捕捉到镶嵌线上，所绘制的镶嵌线必须跨越其他两条镶嵌线，否则无法结束此条镶嵌线的编辑。2.3.6 DOM 建筑等人工地物变形对比分析纠 正 前 后 对 比 图 样纠正前 纠正后进行镶嵌线编辑后要进行正射影像修补，用修补工具围起来变形的房子，右键即可完成修补工作。2.3.7 DOM 色彩对比分析DOM 编辑对于某些地域进行了拼接线调整扭11曲、变形、错位调整局部调色。如下面的 hammer 测区最后生成的 DOM 图所示，航带衔接的航空图片有着明显的区分镶嵌线，如右图所示，明显看出本条航带有三张影像。（hammer 生成的 DOM 影像）2.3.8 最后生成
14、的产品 可见未经房屋过滤前生成的 DEM 图像，房屋屋顶以黄褐色标示，月亮岛小区平地多为黄色或者是浅黄色，地势较为低洼的湖泊池塘多用绿色来渲染，象征了高程的不同。房屋过滤后的 DEM 影像较为单一，没有了大的凸起，进行了一层过滤。（经过房屋过滤后的 DEM）（最终的 Dom 成果图）12第 3 章 DLG 的立测生成3.1 DLG 的概念DLG（ Digital Line Graphic，数字线划地图）是与现有线划基本一致的各地图要素的矢量数据集，且保存了各要素间的空间关系信息和相关的属性信息，是4D 产品的一种。3.2 本次试验生成的 DLG 成果，导出为 DXF 格式3.3 立测的精度检核
15、3.4 DLG 生成1 要求地物丰富2 培养立体感3“符号化配置”成“MapMatrix2007 国标旧版”133.4.1 经验总结在绘制矢量线划图时，注意要绘制一个高程平面的地物，如在控制点附近可以参考控制点的高程坐标。实验的每一步都是建立在上一步的基础上， DEM-DOM-DLG 这样的顺序。至于 4D 的另一个产品 DRG 数字栅格地图则是 DOM,DEM 派生出新的可视化信息。绘制房子与绘制道路平面时要注意在不同的高程，否则立测看起来道路仿佛具有了高程， “浮在”了地面之上，不注意调整高程还会导致不同地物之间的覆盖，不能正确的表达地貌。在这次实验认识到了摄影测量空三成果的基础性以及对于
16、后面的成果有着重大影响，要求我们每一步都要有合适的精度作为生产标准，理解了平差精度的重要。最大的收获是培养了立体感，当真真切切的带着立测眼睛看自己的工作时，感觉自己仿佛置身于正在测绘的无人机上一样，以一种高度来俯视测区，山丘的断崖，房屋的起伏都是真真切切的，航空摄影测量以后将会有着很好的前途，逐步替代低精度要求却外业工作很繁重的测绘任务。实习又认识到了 3R 之间的异同处，无论是通过全站仪，GPS，无人机，都是为了生产出合格的 dlg 地图，都可以导入到 cass 里进行编辑，这都是殊途同归的。第 4 章 编程实习4.1 实验目的在掌握共线方程原理的基础上，自己编写后方交会程序，在计算机上调试
17、，输出计算结果：像片的外方位元素。4.2 具体要求已知值像点坐标为扫描坐标，原点在像片的左上角，单位为 pixel。应将坐标原点平14移至像片中心，作为像平面坐标系（x 轴向右、y 轴向上） ，代入共线方程。像片高 h = 4008 （pixel ）y 方向像宽片 w = 5344 （pixel ）x 方向内方位元素x0 = 47.4857 pixely0 = 12.0276 pixelf = 4547.9352 pixel对共线方程进行线性化，列立误差方程式，系数应采用严密式，按照最小二乘原则进行解。参考结果Xs = 3373.40082 mm； Ys = -141.55657 mm；Zs
18、= 92.77483 mm；Phi = -0.01205弧度；Omega = 0.09997 弧度；Kappa = -0.04101 弧度。4.3 计算原理如图所示，物点 A 和摄影中心 S 在地面摄影测量坐标系中的坐标依次是（X,Y,Z ） 、 （ XS,YS,ZS） ；像点 a 在像空间坐标系中的坐标是（x,y,-f） 。那么由共线条件方程知，其中 ai，bi， ci 是只含三个独立参数 ， 的九个方向余弦。在方程中共有六个未知参数 XS,YS,ZS,， 所以有三个不在一条直线上的已知地面点坐标就可以求出像片的这六个外方位元素。由于共线条件方程是非线性方程，为了便于迭代计算，需要把方程用泰
19、勒级数展开，取一次项得到线型表达式，有)()()( )()()( 333 2220 333 1110 sss sss ZcYbafyXZcYbfx sss sss dZyYdXydyy xxxx )()(15用新的符号表示各偏导数后为 dadaZdYaXyxsss 26252423221 11111)( 其中（x） 、 （y）是函数近似值， , 是外方位元素近似ssZYX,值的改正数，它们的系数为函数的偏导数。为了便于推导，令=X)()()(111 sss ZcYba=Y222 sss=Z)()()(333 sssc对于竖直摄影而言，像片的角方位元素都是小值，因而各系数的近似值为 xafyaf
20、xya yxHHfxaaHfa 262524 111 2322 111 ),(, ,),(,0,0为了提高精度和可靠性，通常需要测量四个或更多的地面控制点和对应的像点坐标，采用最小二乘平差方法解算。此时像点坐标（x,y）作为观测值，加入相应的偶然误差改正数 , ，可列出每个点的误差方程式xvy ysssy xx ldadaZdYaXvd 26252423221 11111用矩阵表示为1626524321 110, aayxZYXvssyxAlV那么由 lxV)(T1TlA4.4 算法解释每一个物方坐标就有一个未知的 x,y 就可以列出两个误差方程，因为有 6 个外方位元素所以只有 3 个点就可
21、以了，本次试验多余的数据作为平差多余计算参与平差解算。4.5 算法流程：获取已知数据 m, x0 , y0 , f , Xtp, Ytp, Ztp 量测控制点像点坐标 x,y 确定未知数初值 Xs0, Ys0, Zs0, 0, 0, 0 组成误差方程式并法化解求外方位元素改正数检查迭代是否收敛具体的流程程序框图如下：输入原始数据像点坐标计算，系统误差改正确定外方位元素初始值组成旋转矩阵 R17是否 否是否是4.6 MATLAB 源代码%aa = textread(物方控制点坐标.txt);fid=fopen(CONTROL.txt); %打开数据总文件B=textscan(fid,%f %f
22、%f %f);%把每一列的数据读入到读入到单元数组B中C=B2 B3 B4; %从单元数组B中提取每列数据赋值给矩阵Cn=max(size(C); %N为矩阵的行总数也就是控制点位个数Xafid=fopen(IAMGE.txt); %打开数据总文件J=textscan(fid,%f %f %f);%把每一列的数据读入到读入到单元数组B中I=J2 J3; %从单元数组B中提取每列数据赋值给矩阵C %确定读入数据的坐标数目 %Cp为平均值Cpx1=sum(C);%控制点的X,Y,Z总和Xap=Cpx1(1,1)/n;%控制点平均X averageYap=Cpx1(1,2)/n;Zap=Cpx1(1
23、,3)/n;%要求待定参数Xs=3373.40082 mm,Ys=-141.55657 mm,Zs=92.77483 mm,Phi,Omega,Kappa%像方坐标，控制点，Xa,Ya,Za,delta%设定外方位元素初始值 x0 = 47.4857 ;逐点组成误差方程式并法化所以像点完否解法方程，求外方位元素改正数计算改正后的外方位元素外方位元素改正数是否小于限差输出计算成果，计算结束迭代次数小于n结束并显示错误信息18y0 = 12.0276 ;f = 4547.9352 ;Xs=3370;Ys=-140;Zs=92;Phi=0;Omega=0;Kappa=0;%RPhi=cos(Phi)
24、,0,-sin(Phi);0,1,0;sin(Phi),0,cos(Phi);%ROmega=1,0,0;0,cos(Omega),-sin(Omega);0,sin(Omega),cos(Omega);%RKappa=cos(Kappa),-sin(Kappa),0;sin(Kappa),cos(Kappa),0;0,0,1;%R=RPhi*ROmega*RKappa;%a1=R(1,1);%a2=R(1,2);%a3=R(1,3);q=Phi;w=Omega;k=Kappa;zy=0;while 1for j=n:-1:1i=n-j+1;a(1)=cos(q)*cos(k)-sin(q)*
25、sin(w)*sin(k);a(2)=-cos(q)*sin(k)-sin(q)*sin(w)*cos(k);a(3)=-sin(q)*cos(w);b(1)=cos(w)*sin(k);b(2)=cos(w)*cos(k);b(3)=-sin(w);c(1)=sin(q)*cos(k)+cos(q)*sin(w)*sin(k);c(2)=-sin(q)*sin(k)+cos(q)*sin(w)*cos(k);c(3)=cos(q)*cos(w);%以i开始迭代Control坐标%开始迭代Xbar(i)=a(1)*(C(i,1)-Xs)+b(1)*(C(i,2)-Ys)+c(1)*(C(i,3
26、)-Zs);Ybar(i)=a(2)*(C(i,1)-Xs)+b(2)*(C(i,2)-Ys)+c(2)*(C(i,3)-Zs);Zbar(i)=a(3)*(C(i,1)-Xs)+b(3)*(C(i,2)-Ys)+c(3)*(C(i,3)-Zs);H(i)=-Zbar(i);%H赋予初值19%误差方程的系数aij改为xsijxs(2*i-1,1)=-f/H(i);xs(2*i-1,2)=0;xs(2*i-1,3)=-(I(i,1)-x0)/H(i);xs(2*i-1,4)=-(f+(I(i,1)-x0)*(I(i,1)-x0)/f);xs(2*i-1,5)=-(I(i,1)-x0)*(I(i,
27、2)-y0)/f;xs(2*i-1,6)=I(i,2)-y0;xs(2*i,1)=0;xs(2*i,2)=-f/H(i);xs(2*i,3)=-(I(i,2)-y0)/H(i);xs(2*i,4)=-(I(i,1)-x0)*(I(i,2)-y0)/f;xs(2*i,5)=-(f+(I(i,2)-y0)*(I(i,2)-y0)/f);xs(2*i,6)=-(I(i,1)-x0);%L=x0-I(i,1),y0-I(i,2);L(2*i-1,1)=x0-I(i,1);L(2*i,1)=y0-I(i,2);%此时L为矩阵2n*1end%=deltaXs;deltaYs;deltaZs;deltaq;
28、deltaw;deltak;N=xs*xs;X=pinv(N)*xs*L;deltaXs=X(1,1);deltaYs=X(2,1);deltaZs=X(3,1);deltaq=X(4,1);deltaw=X(5,1);deltak=X(6,1);Xs=deltaXs+Xs;Ys=deltaYs+Ys;Zs=deltaZs+Zs;q=q+deltaq;w=w+deltaw;k=k+deltak;zy=zy+1;PD(1)=deltaq/pi*60;PD(2)=deltaw/pi*60;PD(3)=deltak/pi*60;if PD(1)2)22Baseimg=rgb2gray(Baseimg
29、);endif (size(size(Warpimg),2)2)Warpimg=rgb2gray(Warpimg);endBimagesizer,Bimagesizec=size(Baseimg);Wimagesizer,Wimagesizec=size(Warpimg);Baseimg=double(Baseimg);Warpimg=double(Warpimg);%高斯滤波BA=fspecial(gaussian,11,11);Baseimg=filter2(BA,Baseimg);Warpimg=filter2(BA,Warpimg);for num_Bpoint=1:size(Bpoi
30、nt,1)Lo=zeros(Num_it+1,1);%Lo(1)=0;% if (size(size(Baseimg),2)2)% Baseimg=rgb2gray(Baseimg);% end% if (size(size(Warpimg),2)2)% Warpimg=rgb2gray(Warpimg);% end% Baseimg=double(Baseimg);% Warpimg=double(Warpimg);% %高斯滤波% BA=fspecial(gaussian,11,11);% Baseimg=filter2(BA,Baseimg);% Warpimg=filter2(BA,W
31、arpimg);%最小二乘匹配前点对Wpr=Wpoint(num_Bpoint,1);Wpc=Wpoint(num_Bpoint,2);Bpr=Bpoint(num_Bpoint,1);Bpc=Bpoint(num_Bpoint,2);if(Wpr-r-1Wimagesizer|Wpc-r-1Wimagesizec)%error(输入错误Wpc出界 );BW(num_Bpoint)=1;continue;endif(Bpr-r-1Bimagesizer|Bpc-r-1Bimagesizec)%改成到四条边直线距离082523%error(输入错误Wpc出界 );BW(num_Bpoint)=1
32、;continue;end%设置初值h0=0;h1=1;a0=Bpr-Wpr;a1=1;a2=0;b0=Bpc-Wpc;b1=0;b2=1;%储存变量hh0=zeros(Num_it+1,1);hh0(1)=h0;hh1=zeros(Num_it+1,1);hh1(1)=h1;aa0=zeros(Num_it+1,1);aa0(1)=a0;aa1=zeros(Num_it+1,1);aa1(1)=a1;aa2=zeros(Num_it+1,1);aa2(1)=a2;bb0=zeros(Num_it+1,1);bb0(1)=b0;bb1=zeros(Num_it+1,1);bb1(1)=b1;b
33、b2=zeros(Num_it+1,1);bb2(1)=b2;%获取初始块，以匹配点为中心，以r为窗口半径for i=1:2*r+3for j=1:2*r+3wx(i)=Wpr-r-2+i;wy(j)=Wpc-r-2+j;Warp1(i,j)=bilinc(Warpimg,wx(i),wy(j);% Bx=Bpr-r-2+i;% By=Bpc-r-2+i;% BaseBxBy(i,j)=bilinc(Baseimg,Bx,By);endendfor i=2:2*r+2for j=2:2*r+2WarpBxBy(i-1,j-1)=bilinc(Warpimg,wx(i),wy(j);24ende
34、nd% BaseBxBy1=BaseBxBy;% Lo(1)=Loxiangguan(Warp1,BaseBxBy);%计算初始块之间相关系数Lo(1)=0.6;%进入迭代环节%for k=1:Num_itfor i=1:2*r+3for j=1:2*r+3Bx=aa0(k)+aa1(k)*wx(i)+aa2(k)*wy(j);By=bb0(k)+bb1(k)*wx(i)+bb2(k)*wy(j);% if( Bxsize(Baseimg,1)|Bx=NaN|Bysize(Baseimg,2)|By=NaN)% break;% endBaseBxBy1(i,j)=bilinc(Baseimg,
35、Bx,By);endendfor i=2:2*r+2for j=2:2*r+2BaseBxBy(i-1,j-1)=hh0(k)+hh1(k)*BaseBxBy1(i,j);endendLo(k+1)=Loxiangguan(WarpBxBy,BaseBxBy);%计算变形参数更改后相关系数if (Lo(k+1)=0.95)break;endif (Lo(k+1)=Lo(k)break;endc=zeros(2*r+1)2,8);%设置初始C矩阵大小，详情请参考武大摄影测量最小二乘影像匹配教材等L=zeros(2*r+1)2,1);%设置初始L矩阵大小，l=g1-g2t=1;for i=2:2*
36、r+2for j=2:2*r+2gy=(BaseBxBy1(i,j+1)-BaseBxBy1(i,j-1)/2;gx=(BaseBxBy1(i+1,j)-BaseBxBy1(i-1,j)/2;g2=BaseBxBy1(i,j);g1=Warp1(i,j);25%c(t,:)=1;g2;hh1(k)*gx;hh1(k)*i*gx;hh1(k)*j*gx;hh1(k)*gy;hh1(k)*i*gy;hh1(k)*j*gy;c(t,:)=1;g2;gx;wx(i)*gx;wy(j)*gx;gy;wx(i)*gy;wy(j)*gy;L(t)=g1-g2;t=t+1;endendLX=inv(c*c)*
37、(c*L);%计算X=dh0,dh1,da0,da1,da2,db0,db1,db2%计算变形参数%hh0(k+1)=hh0(k)+LX(1)+hh0(k)*LX(2);hh1(k+1)=hh1(k)+hh1(k)*LX(2);aa0(k+1)=aa0(k)+LX(3)+aa0(k)*LX(4)+bb0(k)*LX(5);aa1(k+1)=aa1(k)+aa1(k)*LX(4)+bb1(k)*LX(5);aa2(k+1)=aa2(k)+aa2(k)*LX(4)+bb2(k)*LX(5);bb0(k+1)=bb0(k)+LX(6)+aa0(k)*LX(7)+bb0(k)*LX(8);bb1(k+
38、1)=bb1(k)+aa1(k)*LX(7)+bb1(k)*LX(8);bb2(k+1)=bb2(k)+aa2(k)*LX(7)+bb2(k)*LX(8);endif max(Lo)0.95continue;endXT=0;YT=0;GX=0;GY=0;for i=2:2*r+2for j=2:2*r+2gy=(Warp1(i,j+1)-Warp1(i,j-1)/2;gx=(Warp1(i+1,j)-Warp1(i-1,j)/2;XT=XT+wx(i)*gx2;YT=YT+wy(j)*gy2;GX=GX+gx2;GY=GY+gy2;endendXXT=XT/GX;YYT=YT/GY;Wpoin
39、tR(num_Bpoint,1)=XXT;WpointR(num_Bpoint,2)=YYT;26BXT=aa0(k)+aa1(k)*XXT+aa2(k)*YYT;BYT=bb0(k)+bb1(k)*XXT+bb2(k)*YYT;BpointR(num_Bpoint,1)=BXT;BpointR(num_Bpoint,2)=BYT; end%ind=find(BW);index=find(WpointR(:,1)=0);WpointR(index,:)=;BpointR(index,:)=;% Bblkpos(:,1)=BpointR(:,2);% Bblkpos(:,2)=BpointR(:
40、,1);% Wblkpos(:,1)=WpointR(:,2);% Wblkpos(:,2)=WpointR(:,1); % % BW1=Bblkpos Wblkpos;% dlmwrite(E:02cBeijing72457myfilename0824lms.txt, BW1, t);4.9 函数文件function Lo=Loxiangguan(A,B)if size(A)=size(B)error(AB 大小不一样);endr=size(A,2);sum_A=sum(sum(A);sum_B=sum(sum(B);sum_A_B=sum(sum(A.*B);sum_A_A=sum(sum
41、(A.*A);sum_B_B=sum(sum(B.*B);Lo=(sum_A_B-sum_A*sum_B/r/r)/sqrt(sum_A_A-sum_A*sum_A/r/r)*(sum_B_B-sum_B*sum_B/r/r);%Lo=sum_A_B/sqrt(sum_A_A*sum_B_B);function Phuidu=bilinc(Baseimg,Bx,By)Bxz=fix(Bx);Byz=fix(By);detx=Bx-Bxz;dety=By-Byz;if (detx=0 27elsePhuidu=Baseimg(Bxz,Byz)*(1-detx)*(1-dety)+Baseimg(Bxz,Byz+1)*(1-detx)*dety+detx*(1-dety)*Baseimg(Bxz+1,Byz)+detx*dety*Baseimg(Bxz+1,Byz+1);end
