资源描述
课程编号：20162024123 课程性质：专业必修数字摄影测量实习实习报告学院： 测绘学院 专业： 测绘工程 地点： 校 内 班级： 组号： 姓名： 学号： 教师： 2017年 6月 19日 至 2017年 7 月 8日第 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 内定向：确定与物点相对应的像点坐标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编辑的理解82.2.5 Dem拼接裁剪82.3 DOM的生成编辑82.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.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 函数文件25第 1 章 空三报告1.1 NAT空中三角成果数据流程1. 设置测区基本参数2. 建立相机文件，导入老师给的文件3. 导入外业控制点4. 影像格式转换5. 自动转点6. 进行PATB平差计算1.2 空三加密需要注意1. 像素大小、焦距、像主点坐标、畸变参数，像幅尺寸 (数码相机)2. 模型连接点平面及高程限差、地面定向点残差限差、地面检查点残差限差3. 控制点坐标、刺点片及点之记1.3 空三流程图1.4 空三评查步骤1.4.1 需要平差的原因1 同一物体在不同视角下的图像有极大的不同，物体 遮挡而丢失信息2 光照、场景中物体的几何形状和物理特性、物体和摄像机之间的空间关系3 从带畸变和噪声的、二维的、单一的灰度值中提取出尽可能不带畸变和噪声的、三维场景中的诸多因素1.4.2 粗差的影响1观测值中存在粗差，则某一观测的误差对所有残差都有影响2 一个观测的残差包含着所有粗差的影响1.4.3 平差方法-循环粗差剔除1 先外围加4个点进行平差解算像方的粗差点排除到一定程度后，开始加预测的控制点，不要一次性刺入所有预测控制点，建议从外向内回字形一圈圈刺预测控制点，留区域中的3到4个控制点做检核点。每刺入一轮都先平差，编辑得差不多再进行下一轮。2 继续调整控制点关掉粗差探测，勾选验后方差，再平差一次，平差结束。3 上述为精度报告检查平差结果，最后需要到立体测图软件进行检查。1.4.4 自己的方案在本次测区的15个点中，我分了三轮进行平差，每次选取4个点最后保留3个点作为检查点，难处在于有的点有灰尘遮盖造成精度不准，我的精度指标如下所示。第 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 绝对定向：借助物空间坐标为已知的控制点来确定空间辅助坐标系与实际物空间坐标的变换关系称为立体模型的绝对定向。2.1.6 后方交会后方交会：如果知道每张影像的6个外方位元素就可以知道确定被摄物体与航射影像的关系，利用影像范围内的控制点的空间坐标与影像坐标根据共线方程式，反求影像的外方位元素。2.1.7 后方交会绝对定向刺6个点能成功：绝对定向有七个参数，空间相似变换的3个角元素，3个平移量，1个比例尺的缩放。2 个平高+1个高程控制点。独立模型法区域空中三角测量的基本思想，把一个单元模型视为刚体，一个单元模型由两个或更多个立体像对构成，因为一个立体像对需要两张影像，所以如图所示，两张影像的拼接必须要有三个同名像点作为公共点，满足X,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第一条航带的DEM编辑，数字高程模型（Digital Elevation Model)，简称DEM，是通过有限的地形高程数据实现对地面地形的数字化模拟（即地形表面形态的数字化表达），它是用一组有序数值阵列形式表示地面高程的一种实体地面模型。2.2.4 问老师得到的对Dem编辑的理解第一条航带编辑必须要把浮动到眼前，通过3D眼睛看到的明显具有高程的房屋建筑改化到地面高程中去，生成的是地物地貌模型，而不是具有人工建筑的地物模型。将房屋等人工建筑周围的地形高程进行内插，求出该店的实际高程。方法：匹配的点，没贴准地面需要编辑内插的点没贴准地面需要编辑匹配点内插方式：不量测地面的点，用两边的点对中间点内插量测点内插方式：量测地面的点2.2.5 Dem拼接裁剪（Dem成图）2.3 DOM的生成编辑2.3.1 DOM概念DOM-数字正射影像图是对航空（或航天）相片进行数字微分纠正和镶嵌，按一定图幅范围裁剪生成的数字正射影像集。它是同时具有地图几何精度和影像特征的图像。2.3.2 选取新建的DOM产品自动生成调整扫描分辨率与参数 沿影像边缘生成，且原始影像单独生成 DOM2.3.3 规格要求 数据能正确衔接 建筑等无错开、扭曲变形等现象 色调一致、合理1 YLD 裁切后的影像（0.1米）（用DSM生成的拼接后的DOM不用单独提交，主要用于DOM建筑等人工地物的变形分析）2 hammer数据拼接后DOM 0.5米（不用编辑，主要用于DOM色彩分析） 2.3.4 生成72张影像的DOM导入易拼图进行正射影像编辑（影像镶嵌）（待编辑的正射影像）2.3.5 EPT编辑工具使用总结将72张影像生成的DOM拼接导入EPT进行编辑编辑注意要不过水体，不过房屋，可以借助工具进行镶嵌线编辑，添加点、插入点、移动点、删除点、删除临时线和修补选区。注意事项：对于某些四点公用点的镶嵌线编辑1.在与关键点相连的任意一条镶嵌线上单击左键添加点，开始绘制。2.到第四条镶嵌线上，单击左键添加结束点。3.最后单击右键，结束镶嵌线编辑。4.其中绘制镶嵌线时，开始点和结束点必须捕捉到镶嵌线上，所绘制的镶嵌线必须跨越其他两条镶嵌线，否则无法结束此条镶嵌线的编辑。2.3.6 DOM建筑等人工地物变形对比分析纠正前后对比图样纠正前 纠正后进行镶嵌线编辑后要进行正射影像修补，用修补工具围起来变形的房子，右键即可完成修补工作。2.3.7 DOM色彩对比分析DOM编辑对于某些地域进行了拼接线调整扭曲、变形、错位调整局部调色。如下面的hammer测区最后生成的DOM图所示，航带衔接的航空图片有着明显的区分镶嵌线，如右图所示，明显看出本条航带有三张影像。（hammer生成的DOM影像）2.3.8 最后生成的产品 可见未经房屋过滤前生成的DEM图像，房屋屋顶以黄褐色标示，月亮岛小区平地多为黄色或者是浅黄色，地势较为低洼的湖泊池塘多用绿色来渲染，象征了高程的不同。房屋过滤后的DEM 影像较为单一，没有了大的凸起，进行了一层过滤。（经过房屋过滤后的DEM）（最终的Dom成果图）第 3 章 DLG的立测生成3.1 DLG的概念DLG（Digital Line Graphic，数字线划地图）是与现有线划基本一致的各地图要素的矢量数据集，且保存了各要素间的空间关系信息和相关的属性信息，是4D产品的一种。3.2 本次试验生成的DLG成果，导出为DXF格式3.3 立测的精度检核3.4 DLG生成1 要求地物丰富2 培养立体感3“符号化配置”成“MapMatrix2007 国标旧版”3.4.1 经验总结在绘制矢量线划图时，注意要绘制一个高程平面的地物，如在控制点附近可以参考控制点的高程坐标。实验的每一步都是建立在上一步的基础上， DEM-DOM-DLG这样的顺序。至于4D的另一个产品DRG数字栅格地图则是DOM,DEM派生出新的可视化信息。绘制房子与绘制道路平面时要注意在不同的高程，否则立测看起来道路仿佛具有了高程，“浮在”了地面之上，不注意调整高程还会导致不同地物之间的覆盖，不能正确的表达地貌。在这次实验认识到了摄影测量空三成果的基础性以及对于后面的成果有着重大影响，要求我们每一步都要有合适的精度作为生产标准，理解了平差精度的重要。 最大的收获是培养了立体感，当真真切切的带着立测眼睛看自己的工作时，感觉自己仿佛置身于正在测绘的无人机上一样，以一种高度来俯视测区，山丘的断崖，房屋的起伏都是真真切切的，航空摄影测量以后将会有着很好的前途，逐步替代低精度要求却外业工作很繁重的测绘任务。实习又认识到了3R之间的异同处，无论是通过全站仪，GPS，无人机，都是为了生产出合格的dlg地图，都可以导入到cass里进行编辑，这都是殊途同归的。第 4 章 编程实习4.1 实验目的在掌握共线方程原理的基础上，自己编写后方交会程序，在计算机上调试，输出计算结果：像片的外方位元素。4.2 具体要求已知值像点坐标为扫描坐标，原点在像片的左上角，单位为pixel。应将坐标原点平移至像片中心，作为像平面坐标系（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 = 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,，所以有三个不在一条直线上的已知地面点坐标就可以求出像片的这六个外方位元素。由于共线条件方程是非线性方程，为了便于迭代计算，需要把方程用泰勒级数展开，取一次项得到线型表达式，有用新的符号表示各偏导数后为其中（x）、（y）是函数近似值，,是外方位元素近似值的改正数，它们的系数为函数的偏导数。为了便于推导，令=对于竖直摄影而言，像片的角方位元素都是小值，因而各系数的近似值为为了提高精度和可靠性，通常需要测量四个或更多的地面控制点和对应的像点坐标，采用最小二乘平差方法解算。此时像点坐标（x,y）作为观测值，加入相应的偶然误差改正数,，可列出每个点的误差方程式用矩阵表示为那么由 4.4 算法解释每一个物方坐标就有一个未知的x,y就可以列出两个误差方程，因为有6个外方位元素所以只有3个点就可以了，本次试验多余的数据作为平差多余计算参与平差解算。4.5 算法流程：获取已知数据 m, x0 , y0 , f , Xtp, Ytp, Ztp 量测控制点像点坐标 x,y 确定未知数初值 Xs0, Ys0, Zs0, j0, w0, k0 组成误差方程式并法化解求外方位元素改正数检查迭代是否收敛具体的流程程序框图如下：输入原始数据像点坐标计算，系统误差改正确定外方位元素初始值组成旋转矩阵R结束并显示错误信息是逐点组成误差方程式并法化迭代次数小于n所以像点完否否否解法方程，求外方位元素改正数是计算改正后的外方位元素外方位元素改正数是否小于限差否是输出计算成果，计算结束4.6 MATLAB源代码%aa = textread(物方控制点坐标.txt); fid=fopen(CONTROL.txt); %打开数据总文件B=textscan(fid,%f %f %f %f);%把每一列的数据读入到读入到单元数组B中C=B2 B3 B4; %从单元数组B中提取每列数据赋值给矩阵Cn=max(size(C); %N为矩阵的行总数也就是控制点位个数Xa fid=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,3)/n; %要求待定参数Xs=3373.40082 mm,Ys=-141.55657 mm,Zs=92.77483 mm,Phi,Omega,Kappa%像方坐标，控制点，Xa,Ya,Za,delta %设定外方位元素初始值 x0 = 47.4857 ;y0 = 12.0276 ;f = 4547.9352 ;Xs=3370;Ys=-140;Zs=92;Phi=0;Omega=0;Kappa=0; %RPhi=cos(Phi),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 1 for j=n:-1:1i=n-j+1;a(1)=cos(q)*cos(k)-sin(q)*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)-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赋予初值%误差方程的系数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,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;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)0.1&PD(2)0.1&PD(3)2) Baseimg=rgb2gray(Baseimg);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(Bpoint,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,Warpimg); %最小二乘匹配前点对 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; end if(Bpr-r-1Bimagesizer|Bpc-r-1Bimagesizec)%改成到四条边直线距离0825 %error(输入错误Wpc出界 ); BW(num_Bpoint)=1; 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; bb2=zeros(Num_it+1,1); bb2(1)=b2; %获取初始块，以匹配点为中心，以r为窗口半径 for i=1:2*r+3 for j=1:2*r+3 wx(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); end end for i=2:2*r+2 for j=2:2*r+2 WarpBxBy(i-1,j-1)=bilinc(Warpimg,wx(i),wy(j); end end % BaseBxBy1=BaseBxBy; % Lo(1)=Loxiangguan(Warp1,BaseBxBy);%计算初始块之间相关系数 Lo(1)=0.6; %进入迭代环节% for k=1:Num_it for i=1:2*r+3 for j=1:2*r+3 Bx=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;% end BaseBxBy1(i,j)=bilinc(Baseimg,Bx,By); end end for i=2:2*r+2 for j=2:2*r+2 BaseBxBy(i-1,j-1)=hh0(k)+hh1(k)*BaseBxBy1(i,j); end end Lo(k+1)=Loxiangguan(WarpBxBy,BaseBxBy);%计算变形参数更改后相关系数 if (Lo(k+1)=0.95) break; end if (Lo(k+1)=Lo(k)&Lo(k)0.95)%修改一下，判断Lo/0825 %BW(num_Bpoint)=1; break; end c=zeros(2*r+1)2,8);%设置初始C矩阵大小，详情请参考武大摄影测量最小二乘影像匹配教材等 L=zeros(2*r+1)2,1);%设置初始L矩阵大小，l=g1-g2 t=1; for i=2:2*r+2 for j=2:2*r+2 gy=(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); %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; end end LX=inv(c*c)*(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+1)=bb1(k)+aa1(k)*LX(7)+bb1(k)*LX(8); bb2(k+1)=bb2(k)+aa2(k)*LX(7)+bb2(k)*LX(8); end if max(Lo)0.95 continue; end XT=0; YT=0; GX=0; GY=0; for i=2:2*r+2 for j=2:2*r+2 gy=(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; end end XXT=XT/GX; YYT=YT/GY; WpointR(num_Bpoint,1)=XXT;WpointR(num_Bpoint,2)=YYT; BXT=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(:,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(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 & dety=0) Phuidu=Baseimg(Bx,By);else Phuidu=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
展开阅读全文