《遥感原理与应用课程设计》
实习报告
学院: 遥感信息工程学院
班级:
学号:
姓名:
实习地点: 5-325机房
指导教师: 石文轩
2016年 6 月16 日
基于多项式的遥感图像几何纠正
摘要：本次实习大致分为两个部分，第一部分是通过调用OpenCV 函数实现遥感图像的基本处理，我们小组选择的题目是专题二，对彩色图像进行缩放、旋转等几何操作。第二部分是自主编程实现遥感图像处理，我们小组抽到的是基于多项式的遥感图像几何纠正，包括纠正多项式系数的解算，纠正变换函数的建立，坐标转换和灰度重采样的算法原理与编程实现。本文主要介绍实现第二个任务的主要原理，具体步骤，试验结果展示以及遇到的问题分析和心得总结。
关键词：多项式，几何纠正，重采样
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进行比较检查，将最大者视为粗差并将其剔除或重新选点进行平差，直至满足δ<ε的精度评定标准为止。
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,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;
Ncols = (Xmax - Xmin)/30 + 1。
1.2.4 图像纠正变换函数的目的是确定原始图像的坐标值与地面坐标值之间的变换关系，为了进一步表达原始图像与输出图像之间的对应关系，我们需要把地面坐标进行转换，首先转换为输出图像坐标，转换公式如下：
在输出图像的边界坐标确定之后，就可以建立坐标系统并按照纠正变换函数把输入图像的每个像素逐个变换到输出图像空间中，这一步有两种可供选择的纠正方案，直接法和间接法。直接法的原理是对应输入图像的每一个像素按照变换函数找到其在输出图像上的对应位置，间接法的原理是对于输出图像的每一个像素点，反算得到其对应的输入图像的像素。
1.3 数字图像亮度值重采样
1.3.1 最邻近像元法
最邻近像元采样发的实质是采取距离采样点最近的像素的亮度值作为重采样的亮度值。
这种方法的采样方法最简单，辐射保真度比较好，但是这种方法将会造成像素点在一个像素范围之内的位移，与双线性内插法、三次卷积法相比，几何精度较差。
1.3.2 双线性内插法
这种方法的原理可以用下图所示的一个三角形线性函数表达：
图1 双线性内插法原理图
在使用双线性内插采样法时，需要综合计算被采样点周围四个已知像素点
的亮度值。
这种方法的计算较为简单，与最邻近像元法相比具有更好的亮度采样精
度，但图像略变模糊。因此双线性内插采样法也是实践中常用的方法。
1.4 纠正结果评价
即通过查阅相关资料和相应的规范，对于进行了几何纠正后的影像结果进行精度评价，如果符合精度的要求，则完成了遥感图像的几何纠正操作。
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();pointMat.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,POINT[k]);
k++;
}
}
k = 0;
for (i = 0;i<7;i++)
{
for (j = 0;j<2;j++)
{
pointMat.write(i,j,point[k]);
k++;
}
}
cal.transpos(&POINTMat,&POINTTMat);//矩阵求逆
//最小二乘法求矩阵系数，结果存在了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.multiply(&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;
//得到输出图像的行列号，已知该图像分辨率为30m
int M = (Y2 - Y1)/ 30 + 1;
int N = (X2 - X1) / 30 + 1;
unsigned char * im = new unsigned char[M*N];//存储输出图像像素
for(i = 0 ; i < M;i++)
{
for(j = 0 ; j { im[i*N + j] = 0; } }//先赋值为0 3. 数字图像亮度值重采样 数字图像亮度重采样的方法有很多，这里采用了两种方法。 （1）最邻近像元法，采取距离被采样点最近的像素值作为采样亮度。for( i = 1 ; i <= M ;i++)//重采样 { for( j = 1 ; j<= N;j++) { double tempX = X1 + (j - 1)*30; double tempY = Y2 - (i - 1)*30; double geo[6] = {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,geo[whatk]); whatk ++; } _Matrix imgMat(1,2); imgMat.init_matrix(); cal.multiply(&geoMat,&solveMat,&imgMat); double imgX = imgMat.arr[0];//原始图像坐标 double imgY = imgMat.arr[1]; 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 + (j - 1)] = image.data[t1 + t2* col] ; } geoMat.free_matrix(); imgMat.free_matrix(); } } （2）双线性内插法。对被采样点周围的四个像素的亮度值进行双线性内插，具体原理可以用图1所示的三角形线性函数来表示。 for(i=0;i for(j=0;j calculate(&x1,&y1,Xp,Yp,a,b); if(x1>=xa&&x1<=xc && y1>=ya&&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 = lpData[nWidth*int(r) + int(c)]; I12 = lpData[nWidth*int(r) + int(c)+1]; I21 = lpData[nWidth*(int(r)+1) + int(c)]; I22 = lpData[nWidth*(int(r)+1) + int(c)+1]; lpRData[N*i + j] = Wx1*Wy1*I11 + Wx2*Wy1*I21 + Wx1*Wy2*I12 + Wx2*Wy2*I22;} else lpRData[N*i + j] = 0;} } 3 软件使用方法 1. 打开Microsoft Visual Studio 2010 软件，新建一个Win32控制台应用程序。 2. 配置相应的OpenCV环境。 首先将OpenCV2.2 文件夹放在与工程根目录平行的目录之下。 回到工程，在“解决方案资源管理器”中，右键单击工程名，弹出菜单中点选“属性”。 在弹出的属性页对话框中，选择“VC++目录”-点击“包含目录”右边的倒三角-点击“编辑”。 在弹出的包含目录对话框中，分别加入“...\OpenCV2.2\include”，“...\OpenCV2.2\include\opencv”,“...\OpenCV2.2\include\openc v2”三个文件夹。最后点击“确定”。 点击“库录”右边的倒三角-点击“编辑”。在弹出的库目录对话框中，加入“...\OpenCV2.2\lib”文件夹。最后点击“确定” 在属性页对话框中，选择“链接器”-“输入”。点击“附加依赖项”右边的倒三角-点击“编辑”。在弹出的“附加依赖项”对话框中，分别填入如下文件名：opencv_core220d.lib opencv_highgui220d.lib opencv_video220d.lib opencv_ml220d.lib opencv_legacy220d.lib opencv_imgproc220d.lib 点击“确定”。最后点击“应用”-“确定”。 通过上述方法的操作，即完成了OpenCV环境的配置。 3. 程序的编写。基于多项式的遥感图像几何纠正处理可以大致分为三个步骤：通过地面已知控制点解算多项式系数、遥感图像的纠正变换、数字图像亮度值的重采样。 4．实验结果与分析 基于多项式的遥感图像几何纠正处理前后的效果图如下： 图2：几何纠正效果图 从这张结果图像可以看出，原先具有几何变形的图像得到了一定程度的纠正，图像信息得到了良好的还原。 在亮度重采样的方法选择上，我们采用了最邻近像元法和双线性插值法两种方法进行对比。效果图如下： 图3：最邻近像元法图4：双线性插值法从对比上可以看出，图像的几何变形都得到了较为完善的纠正，目视判读的纠正结果基本一致，但是仔细对比可以发现后者相对于前者在图像细节上要模糊一些，但是根据理论后者的几何精度更高。因此纠正方法的选择还是应该综合考虑实际因素。 5 讨论和总结 5.1 小组分工 本实习小组共有贾丹、刘爽、石东博和吴京航四人。 贾丹主要负责答辩PPT 的主要制作。 刘爽主要负责答辩PPT的参与制作，实习报告的参与编写。 石东博主要负责实习报告的主要编写，答辩PPT的参与制作。 吴京航主要负责用于答辩展示的代码编写。 此外组内四人都已各自独立完成彩色图像的几何变形、基于多项式的遥感图像几何纠正的代码的编写。 5.2 实习总结 本次遥感原理课程设计的编程部分为期五天，共六个课时，我们的收获有很多。 首先对于第一个任务，也就是彩色图像的缩放、旋转和平移等操作，这个任务的要求是需要调用OpenCV库自带的函数来完成代码的编写，由于我们在之前的数字图象处理的实习中有接触过OpenCV对图像进行处理，因此对于OpenCV具有一定的了解，而这次实习通过给同学们任务让同学们自主查阅调用函数的方法来完成任务，则是让我们对于遥感图像的处理有了更深层次的了解。 举例来说，对于彩色图像缩放函数，在查阅相关资料的时候我就发现不只有一种函数具有图像缩放功能，比如resize（），Resize（），虽然看上去只有一个字母的差异，不过在使用的过程中，可以发现在参数调用，函数使用上还是具有一定程度的区别的，这就需要我们通过比对和学习选择更加适合的函数来完成彩色图像几何处理的任务。 这次实习我们用了绝大多数时间来完成第二个任务，即基于多项式的遥感图像几何校正。这一部分的知识我们在遥感原理与应用的理论课程上已经进行了比较详细的学习，因此在算法架构的理解上的难度并不是很大。在重温了课本上的基本原理并进行系统的理解后，我们开始了代码的编写。 第一步操作是利用地面已知控制点来求解多项式的系数。这部分几乎完全是数学公式的代码实现。多项式次数的选择有一次，二次和三次，而相对应的系数的个数分别是三个，六个和十个，对应的公式是k = (n+1)(n+2)/2。这次实习给出的控制点一共有七个，因此可以用一次或二次纠正。在这里我们小组的成员使用的方法不尽相同。在求解多项式系数矩阵的过程中，主要涉及的操作即矩阵的转置、求逆等等运算，所以这一部分的重点难点主要在于这些矩阵运算的实现。 第二步操作是对于图像进行几何纠正。首先需要确定输出地图的边界坐标值。这里用到的方法即利用上一步结算处的多项式系数，把给出的原图端点坐标进行转换得到输出图像地图的坐标。这里的变换步骤比较琐碎，一定要十分细心和认真才能得到正确的结果，否则很可能会出现坐标反转、区域大小不正确等问题。这里稍不留神就可能出现上述失误。因此我们感悟到，在代码的编写过程中，算法流程的掌握，编码经验都很重要，但是一颗精益求精的心，一个谨小慎微的态度，也是同样重要的。 第三步是在确定了输出图像边界后进行纠正变换。这里有两种方法可以采用，直接法和间接法。这两种方法的原理大致相同，只是函数的参数选择略有不同。这里我们采用的是间接法，从输出图像的像素点出发，寻找对应的原始图像坐标进行纠正处理。 第四步是数字图像的灰度值重采样。这里我们学习了三种方法，分别是最邻近像元法、双线性内插法和双三次卷积法。这三种方法的原理相似，都是通过对邻近像元亮度的综合取舍进行输出。最邻近像元法是采取距离最接近的像素值，双线性内插法则综合考虑了周围四个像素的像素值进行加权处理，双三次卷积法则是综合考虑周围4x4范围内的像素值。这三种方法的几何精度越来越高，计算精度也在逐渐增大，而由于采用的更大的模版进行像素取值则会使得图像模糊度变高。这里我 们采用了前两种方法，在输出结果上可以看出比较明显的区别，几何精度上双线性内插法应该比最邻近像元法更高，但是目视难以判读出区别，模糊程度上双线性内插法与最邻近像元法相比具有明显的模糊。因此在实际应用上还是要综合权衡各方面因素，选取最适合的方法。 5.3 遇到的问题 在本次实习中我们遇到了如下一些问题： 首先是方法的选择。在纠正法案上有直接法和间接法两种方法，纠正次数上有一次和二次两种方法，亮度重采样上有最邻近像元法等三种方法，在面临这些选择的时候，我们查阅了大量资料，对于不同的选择进行多方面比对，并最终选择了更为适合本次实习的一些方法，在辩证的学习中我们也对于这些方法有了更深入的理解。更难能可贵的是，编写程序时遇到的方法选择，和我们在人生中遇到的诸多抉择是共通的，只有进行了多方面的调查，并认真考虑自身因素，才能筛选出最适合自己的选择。 其次是在纠正变换的过程中，坐标变换方面是一个困扰我们的问题。虽然这里并不存在逻辑上知识上的漏洞，问题仅仅出在编码人的细心与认真程度上，但这个问题也绝对不能被忽视。在代码的编写过程中，绘制出简要的流程图是十分必要的，因此我们得到的感悟是，编码所需要的并不仅仅是一个大脑，一台电脑，更需要一支笔一个本，做更为精细的记录。 5.4 实习的建议 这次实习的时间安排正好处于考试周，导致学生需要在复习计划和考试安排之间合理的安排时间，而往往会让学生顾此失彼。希望实习时间能够尽可能的避开考试。我的建议是实习的时间能够安排的更加动态，比如在闲暇的时候集中实习，事情多的时候少安排课时，尽管现实因素也有很多比如机房的预定，考试的不确定性，不过让实习和考试能够分开来进行是绝大多数同学的心愿。 本次编程实习共安排了十五个专题任务，大部分都是和遥感原理与相关方面的应用密切相关的，可是同学们真正动手实践过的只有1到2个，那么对于其他的专题并没有得到深入的理解和动手实现的机会，希望能给同学们安排更多的专题任务，对遥感原理算法方面有更全面的掌握。 参考文献 [1] 孙家抦，遥感原理与应用，武汉：武汉大学出版社，2013 [2] 贾永红，张谦，崔卫红，余卉，数字图像处理实习教程，2012 综合评语： 姓名 考 勤 代 码 检 查 PPT 与实习报 告 总 成 绩10% 3 0% 20% 指导教师： 年月日
相关主题
文本预览