遥感原理实习报告openCV

1、遥感原理与应用课程设计实习报告 学 院:遥感信息工程学院 班 级: 学 号: 姓 名: 实习地点: 5-325机房 指导教师: 石文轩 2016年 6 月 16 日基于多项式的遥感图像几何纠正摘 要： 本次实习大致分为两个部分，第一部分是通过调用OpenCV函数实现遥感图像的基本处理，我们小组选择的题目是专题二，对彩色图像进行缩放、旋转等几何操作。第二部分是自主编程实现遥感图像处理，我们小组抽到的是基于多项式的遥感图像几何纠正，包括纠正多项式系数的解算，纠正变换函数的建立，坐标转换和灰度重采样的算法原理与编程实现。本文主要介绍实现第二个任务的主要原理，具体步骤，试验结果展示以及遇到的问题分析和

2、心得总结。关键词：多项式，几何纠正，重采样1 原理介绍1.1 利用地面已知控制点求解多项式系数一般多项式纠正的变换公式如下：上式中，x，y为像素的原始图像坐标，X，Y为同名像素点的地面坐标。利用地面已知控制点求解多项式系数的具体步骤如下：1.1.1 列出误差方程式1.1.2 构成法方程1.1.3计算多项式系数1.1.4 精度判定其中，n为控制点的个数，N为纠正系数的个数，n-N 为多余观测数。评定精度的步骤是设定一个限差，作为精度评定的标准，如果通过上式计算出的>,则说明计算结果存在一定程度的粗差，精度超限结果不可取。并对每一个控制点上的平差残余误差Vx，Vy进行比较检查，将最大者视为粗

3、差并将其剔除或重新选点进行平差，直至满足<的精度评定标准为止。1.2 遥感图像的几何纠正变换在进行遥感图像的几何纠正变换之前，需要先在计算机中为输出图像开辟一定的存储空间，即确定输出图像的空间边界以及地图地面坐标值。1.2.1 把原始图像的四个角点坐标按照纠正变换函数计算出投影到地图坐标系统中的坐标值：(X1,Y1), (X2,Y2), (X3,Y3), (X4,Y4)。1.2.2 对上述四个坐标按照横、纵坐标分别分组并求出最大最小值，所得结果即为纠正后图像的边界的地图坐标系的值。Xmin = min (min (X1,X2),min (X3,X4)Xmax= max (max (X1,

4、X2), max (X3,X4)Ymin = min (min (Y1,Y2),min (Y3,Y4)Ymax = max (max (y1,Xy2), max (y3,y4)1.2.3 在利用计算机进行纠正后图像的存储之前，需要把上述边界范围转换为计算机中的储存数组空间，将边界范围划分为格网，每一个格网点代表一个像素值，因此需要根据精度要求定义输出图像的地面尺寸。这里的地面点尺寸即图像分辨率为30米。随后利用给出的左上角边界点为输出图像的坐标原点，并利用其他的坐标信息定义图像的行数和列数。最后计算得到的图像总的行列数如下所示：Nrows = (Ymax - Ymin)/30 + 1;Ncol

5、s = (Xmax - Xmin)/30 + 1。1.2.4 图像纠正变换函数的目的是确定原始图像的坐标值与地面坐标值之间的变换关系，为了进一步表达原始图像与输出图像之间的对应关系，我们需要把地面坐标进行转换，首先转换为输出图像坐标，转换公式如下：在输出图像的边界坐标确定之后，就可以建立坐标系统并按照纠正变换函数把输入图像的每个像素逐个变换到输出图像空间中，这一步有两种可供选择的纠正方案，直接法和间接法。直接法的原理是对应输入图像的每一个像素按照变换函数找到其在输出图像上的对应位置，间接法的原理是对于输出图像的每一个像素点，反算得到其对应的输入图像的像素。1.3 数字图像亮度值重采样1.3.1

6、 最邻近像元法最邻近像元采样发的实质是采取距离采样点最近的像素的亮度值作为重采样的亮度值。这种方法的采样方法最简单，辐射保真度比较好，但是这种方法将会造成像素点在一个像素范围之内的位移，与双线性内插法、三次卷积法相比，几何精度较差。1.3.2 双线性内插法这种方法的原理可以用下图所示的一个三角形线性函数表达：图1 双线性内插法原理图在使用双线性内插采样法时，需要综合计算被采样点周围四个已知像素点的亮度值。这种方法的计算较为简单，与最邻近像元法相比具有更好的亮度采样精度，但图像略变模糊。因此双线性内插采样法也是实践中常用的方法。1.4 纠正结果评价即通过查阅相关资料和相应的规范，对于进行了几何纠

7、正后的影像结果进行精度评价，如果符合精度的要求，则完成了遥感图像的几何纠正操作。2 实现方法与过程1. 利用地面已知控制点求解多项式系数总体思路是利用最小二乘法求解多项式系数，通过建立误差方程，建立法方程，求解多项式的系数，其中矩阵相关的运算封装在函数里。/求取变换系数/初始化数组double POINT = ;double point = ;/建立矩阵并开辟空间_Matrix POINTMat(7,6),pointMat(7,2),solveMat(6,2),POINTTMat(6,7);POINTTMat.init_matrix();POINTMat.init_matrix();point

8、Mat.init_matrix();solveMat.init_matrix(); _Matrix_Calc cal;int i,j,k;/初始化矩阵k=0;for (i = 0 ; i < 7;i+)for (j = 0 ;j< 6;j+)POINTMat.write(i,j,POINTk);k+;k = 0;for (i = 0;i<7;i+)for (j = 0;j<2;j+)pointMat.write(i,j,pointk);k+;cal.transpos(&POINTMat,&POINTTMat);/矩阵求逆/最小二乘法求矩阵系数，结果存在了

9、solveMat中，为6行二列的矩阵_Matrix mul1Mat(6,6);mul1Mat.init_matrix();cal.multiply(&POINTTMat,&POINTMat,&mul1Mat);_Matrix mul2Mat(6,2);mul2Mat.init_matrix();cal.multiply(&POINTTMat,&pointMat,&mul2Mat);_Matrix invMat(6,6);invMat.init_matrix();cal.inverse(&mul1Mat,&invMat);cal.m

10、ultiply(&invMat,&mul2Mat,&solveMat);2. 遥感图像的几何纠正变换首先确定输出地图坐标的边界坐标，通过给出的输入图像左上角端点的坐标值求解出原始图像的各个端点坐标，通过多项式进行转换得到输入地图的端点坐标。Mat image;image = imread("./left-3.bmp",0);if(image.empty()return;getArea(X1,Y1,X2,Y2);/得到输出图像的地理坐标范围int row,col;row = image.rows;col = image.cols; /得到输出图像的行列

11、号，已知该图像分辨率为30mint M = (Y2 - Y1)/ 30 + 1;int N = (X2 - X1) / 30 + 1;unsigned char * im = new unsigned charM*N;/存储输出图像像素for(i = 0 ; i < M;i+)for(j = 0 ; j<N;j+)imi*N + j = 0;/先赋值为03. 数字图像亮度值重采样数字图像亮度重采样的方法有很多，这里采用了两种方法。（1）最邻近像元法，采取距离被采样点最近的像素值作为采样亮度。for( i = 1 ; i <= M ;i+)/重采样 for( j = 1 ; j

12、<= N;j+) double tempX = X1 + (j - 1)*30; double tempY = Y2 - (i - 1)*30;double geo6 = 1,tempX,tempY,tempX*tempY,tempX*tempX,tempY*tempY; _Matrix geoMat(1,6);geoMat.init_matrix();int what;int whatk = 0;for (what = 0 ; what < 6;what+) geoMat.write(0,what,geowhatk); whatk +;_Matrix imgMat(1,2);im

13、gMat.init_matrix();cal.multiply(&geoMat,&solveMat,&imgMat);double imgX = imgMat.arr0;/原始图像坐标double imgY = imgMat.arr1;int t1 = (int(imgX+0.5) - 800);int t2 = (-(int(imgY+0.5 ) + 3600);/最邻近插值if(t1 >= 0 &&t1 < col && t2 >= 0 && t2 < row) im(i - 1) * N +

14、(j - 1) = image.datat1 + t2* col ;geoMat.free_matrix();imgMat.free_matrix(); （2）双线性内插法。对被采样点周围的四个像素的亮度值进行双线性内插，具体原理可以用图1所示的三角形线性函数来表示。for(i=0;i<M;i+)for(j=0;j<N;j+) Xp = X1 + j*30; Yp = Y2 - i*30;calculate(&x1,&y1,Xp,Yp,a,b);if(x1>=xa&&x1<=xc && y1>=ya&&

15、;y1<=yb)double r,c;c = x1 - xa;r = yb - y1;double deltax,deltay;deltax = int(r)+1-r;deltay = int(c)+1-c;double Wx1,Wy1,Wx2,Wy2;Wx1 = 1-deltax; Wy1 = 1-deltay;Wx2 = deltax; Wy2 = deltay;int I11,I12,I21,I22;I11 = lpDatanWidth*int(r) + int(c);I12 = lpDatanWidth*int(r) + int(c)+1;I21 = lpDatanWidth*(

16、int(r)+1) + int(c);I22 = lpDatanWidth*(int(r)+1) + int(c)+1;lpRDataN*i + j = Wx1*Wy1*I11 + Wx2*Wy1*I21 + Wx1*Wy2*I12 + Wx2*Wy2*I22;else lpRDataN*i + j = 0;3 软件使用方法1. 打开Microsoft Visual Studio 2010 软件，新建一个Win32控制台应用程序。2. 配置相应的OpenCV环境。首先将OpenCV2.2 文件夹放在与工程根目录平行的目录之下。回到工程，在“解决方案资源管理器”中，右键单击工程名，弹出菜单中点选

17、“属性”。在弹出的属性页对话框中，选择“VC+目录”-点击“包含目录”右边的倒三角-点击“编辑”。在弹出的包含目录对话框中，分别加入“.OpenCV2.2include”，“.OpenCV2.2includeopencv”,“.OpenCV2.2includeopencv2”三个文件夹。最后点击“确定”。点击“库录”右边的倒三角-点击“编辑”。在弹出的库目录对话框中，加入“.OpenCV2.2lib”文件夹。最后点击“确定” 在属性页对话框中，选择“链接器”-“输入”。点击“附加依赖项”右边的倒三角-点击“编辑”。在弹出的“附加依赖项”对话框中，分别填入如下文件名： opencv_core22

18、0d.libopencv_highgui220d.libopencv_video220d.libopencv_ml220d.libopencv_legacy220d.libopencv_imgproc220d.lib点击“确定”。最后点击“应用”-“确定”。通过上述方法的操作，即完成了OpenCV环境的配置。3. 程序的编写。基于多项式的遥感图像几何纠正处理可以大致分为三个步骤：通过地面已知控制点解算多项式系数、遥感图像的纠正变换、数字图像亮度值的重采样。4实验结果与分析基于多项式的遥感图像几何纠正处理前后的效果图如下：图2：几何纠正效果图从这张结果图像可以看出，原先具有几何变形的图像得到了一

19、定程度的纠正，图像信息得到了良好的还原。在亮度重采样的方法选择上，我们采用了最邻近像元法和双线性插值法两种方法进行对比。效果图如下： 图3：最邻近像元法 图4：双线性插值法从对比上可以看出，图像的几何变形都得到了较为完善的纠正，目视判读的纠正结果基本一致，但是仔细对比可以发现后者相对于前者在图像细节上要模糊一些，但是根据理论后者的几何精度更高。因此纠正方法的选择还是应该综合考虑实际因素。5 讨论和总结 5.1 小组分工本实习小组共有贾丹、刘爽、石东博和吴京航四人。贾丹主要负责答辩PPT 的主要制作。刘爽主要负责答辩PPT的参与制作，实习报告的参与编写。石东博主要负责实习报告的主要编写，答辩PP

20、T的参与制作。吴京航主要负责用于答辩展示的代码编写。此外组内四人都已各自独立完成彩色图像的几何变形、基于多项式的遥感图像几何纠正的代码的编写。5.2 实习总结本次遥感原理课程设计的编程部分为期五天，共六个课时，我们的收获有很多。首先对于第一个任务，也就是彩色图像的缩放、旋转和平移等操作，这个任务的要求是需要调用OpenCV库自带的函数来完成代码的编写，由于我们在之前的数字图象处理的实习中有接触过OpenCV对图像进行处理，因此对于OpenCV具有一定的了解，而这次实习通过给同学们任务让同学们自主查阅调用函数的方法来完成任务，则是让我们对于遥感图像的处理有了更深层次的了解。举例来说，对于彩色图像

21、缩放函数，在查阅相关资料的时候我就发现不只有一种函数具有图像缩放功能，比如resize（），Resize（），虽然看上去只有一个字母的差异，不过在使用的过程中，可以发现在参数调用，函数使用上还是具有一定程度的区别的，这就需要我们通过比对和学习选择更加适合的函数来完成彩色图像几何处理的任务。这次实习我们用了绝大多数时间来完成第二个任务，即基于多项式的遥感图像几何校正。这一部分的知识我们在遥感原理与应用的理论课程上已经进行了比较详细的学习，因此在算法架构的理解上的难度并不是很大。在重温了课本上的基本原理并进行系统的理解后，我们开始了代码的编写。第一步操作是利用地面已知控制点来求解多项式的系数。这部

22、分几乎完全是数学公式的代码实现。多项式次数的选择有一次，二次和三次，而相对应的系数的个数分别是三个，六个和十个，对应的公式是k = (n+1)(n+2)/2。这次实习给出的控制点一共有七个，因此可以用一次或二次纠正。在这里我们小组的成员使用的方法不尽相同。在求解多项式系数矩阵的过程中，主要涉及的操作即矩阵的转置、求逆等等运算，所以这一部分的重点难点主要在于这些矩阵运算的实现。第二步操作是对于图像进行几何纠正。首先需要确定输出地图的边界坐标值。这里用到的方法即利用上一步结算处的多项式系数，把给出的原图端点坐标进行转换得到输出图像地图的坐标。这里的变换步骤比较琐碎，一定要十分细心和认真才能得到正确

23、的结果，否则很可能会出现坐标反转、区域大小不正确等问题。这里稍不留神就可能出现上述失误。因此我们感悟到，在代码的编写过程中，算法流程的掌握，编码经验都很重要，但是一颗精益求精的心，一个谨小慎微的态度，也是同样重要的。第三步是在确定了输出图像边界后进行纠正变换。这里有两种方法可以采用，直接法和间接法。这两种方法的原理大致相同，只是函数的参数选择略有不同。这里我们采用的是间接法，从输出图像的像素点出发，寻找对应的原始图像坐标进行纠正处理。第四步是数字图像的灰度值重采样。这里我们学习了三种方法，分别是最邻近像元法、双线性内插法和双三次卷积法。这三种方法的原理相似，都是通过对邻近像元亮度的综合取舍进行输出。最邻近像元法是采取距离最接近的像素值，双线性内插法则综合考虑了周围四个像素的像素值进行加权处理，双三次卷积法则是综合考虑周围4x4范围内的像素值。这三种方法的几何精度越来越高，计算精度也在逐渐增大，而由于采用的更大的模版进行像素取值则会使得图像模糊度变高。这里我们采用了前两种方法，在输出结果上可以看出比较明显的区别，几何精度上双线性内插法应该比最邻近像元法更高，但是目视难以判读出区别，模糊程度上双线性内插法与最邻近像元法相比具有明显的模糊。因此在实际应用上还是要综合权衡各方面因素，选取最适合的方法。5.3 遇到的问题在本次实习中我们遇到