数学建模：数码相机定位模型.doc

高教社杯全国大学生数学建模竞赛承 诺 书我们仔细阅读了中国大学生数学建模竞赛的竞赛规则.我们完全明白，在竞赛开始后参赛队员不能以任何方式（包括电话、电子邮件、网上咨询等）与队外的任何人（包括指导教师）研究、讨论与赛题有关的问题。我们知道，抄袭别人的成果是违反竞赛规则的, 如果引用别人的成果或其他公开的资料（包括网上查到的资料），必须按照规定的参考文献的表述方式在正文引用处和参考文献中明确列出。我们郑重承诺，严格遵守竞赛规则，以保证竞赛的公正、公平性。如有违反竞赛规则的行为，我们将受到严肃处理。我们参赛选择的题号是（从A/B/C/D中选择一项填写）： 我们的参赛报名号为（如果赛区设置报名号的话）： 所属学校（请填写完整的全名）： 参赛队员 (打印并签名) ：1. 2. 3. 指导教师或指导教师组负责人 (打印并签名)： 日期： 年 月 日赛区评阅编号（由赛区组委会评阅前进行编号）：高教社杯全国大学生数学建模竞赛编 号 专 用 页赛区评阅编号（由赛区组委会评阅前进行编号）：赛区评阅记录（可供赛区评阅时使用）：评阅人评分备注全国统一编号（由赛区组委会送交全国前编号）：全国评阅编号（由全国组委会评阅前进行编号）：数码相机定位模型摘要 本文利用光学原理建立物体成像位置关系的数学模型，解决了相机系统标定中的如何确定靶标圆心对应像坐标和固定相机的相对位置等问题。针对如何确定靶标圆心对应像坐标问题，我们分析讨论了各种影像分析方法，基于针孔成像原理和相机坐标系变换矩阵的思想，建立了线性相机模型，利用透视变换矩阵确定摄像机的位置标定思想给出算法，通过找靶标和像平面上对应公切线切点的方法取点，利用取定的点和算法求得相机成像内外参数，利用MATLAB编程求解，求得靶标圆心对应像坐标。结果如下：在问题1要求的三维坐标系下A,B,C,D,E圆心对应的像坐标是（单位：mm）：（-50.78，50.93，417.20），（-24.42，48.88，417.20），（32.77，44.44，417.20），（18.33，-32.20，417.20），（-60.37，-31.53，417.20）.针对设计一种方法检验我们模型的问题，我们采用最小二乘法拟合椭圆，确定中心点，改变所选取的标定点，求得多组中心的值，比较每组中心值与拟合的中心点，差别不大。并对方法的精度和稳定性进行了讨论。针对如何确定两部固定相机相对位置的问题，以物坐标为中间变量，通过坐标矩阵的变换建立模型，给出算法，最终确定两部相机的相对位置。本论文通过两种方法对比分析，结果相对吻合，模型的可信度较高。考虑到畸变因素的影响，本论文进行了模型改进，使模型具有更高的精度，有较强的推广性。关键词：针孔成像 线性相机模型 坐标系变换矩阵 公切线切点 最小二乘法拟合椭圆1、问题重述数码相机定位在交通监管（电子警察）等方面有广泛的应用。所谓数码相机定位是指用数码相机摄制物体的相片确定物体表面某些特征点的位置。最常用的定位方法是双目定位，即用两部相机来定位。对物体上一个特征点，用两部固定于不同位置的相机摄得物体的像，分别获得该点在两部相机像平面上的坐标。只要知道两部相机精确的相对位置，就可用几何的方法得到该特征点在固定一部相机的坐标系中的坐标，即确定了特征点的位置。于是对双目定位，精确地确定两部相机的相对位置就是关键，这一过程称为系统标定。标定的一种做法是：在一块平板上画若干个点， 同时用这两部相机照相，分别得到这些点在它们像平面上的像点，利用这两组像点的几何关系就可以得到这两部相机的相对位置。然而，无论在物平面或像平面上我们都无法直接得到没有几何尺寸的“点”。实际的做法是在物平面上画若干个圆（称为靶标），它们的圆心就是几何的点了。而它们的像一般会变形，如图1所示，所以必须从靶标上的这些圆的像中把圆心的像精确地找到，标定就可实现。 图 1 靶标上圆的像有人设计靶标如下，取1个边长为100mm的正方形，分别以四个顶点（对应为A、C、D、E）为圆心，12mm为半径作圆。以AC边上距离A点30mm处的B为圆心，12mm为半径作圆，如图2所示（略）。用一位置固定的数码相机摄得其像，如图3所示（略）。问题：（1） 建立数学模型和算法以确定靶标上圆的圆心在该相机像平面的像坐标, 这里坐标系原点取在该相机的焦点，x-y平面平行于像平面；（2） 对由图2、图3分别给出的靶标及其像，计算靶标上圆的圆心在像平面上的像坐标, 该相机的像距（即焦点到像平面的距离）是1577个像素单位(1毫米约为3.78个像素单位)，相机分辨率为1024786；（3） 设计一种方法检验你们的模型，并对方法的精度和稳定性进行讨论；（4） 建立用此靶标给出两部固定相机相对位置的数学模型和方法。2、模型假设2.1假设相机镜头（凸透镜）前后两侧弧度、介质均相同，即镜头前后两侧折射率相同。2.2 假设相机镜头是很完善的，没有缺损或倾斜等外部形变。2.3 假设物平面和像平面的中心与镜头的光学中心在同一条直线（光轴）上。2.4 假设物平面的点经透镜折射后都可以被像平面接到。2.5 假设问题1中靶标和像平面的图形是已知的。2.6 像平面上点的像坐标的原点为该相机的光学中心，x-y平面平行与像平面。2.7 假设不考虑照相机畸变。2.8 有效焦距为相机镜头光学中心到像平面的距离，即像距。3、符号说明：图像像素坐标系 ：像点像素坐标：图像物理坐标系：像点物理坐标：物点相机坐标系：世界坐标系：物点世界坐标：旋转正交变换矩阵：平移变换矩阵：有效焦距：为相机成像内部参数阵：为相机成像外部参数阵：称为投影矩阵4、模型建立及求解4.1问题14.1.1 问题分析本问题的整体思路是：在仅已知物平面和像平面，相机的其他属性参数未知的前提下，建立照相机模型，确定靶标（物平面）上的点到该照相机像平面的像坐标的映射关系，其中包括照相机成像系统内外几何及光学参数的标定，从而达到只要给定标靶上的点，就可直接通过模型求得其像坐标的目的。 在具体建立照相机模型时，由于照相机与摄相机在成像原理上没有本质区别，因此沿用现有比较成熟的摄机标定技术的理论思想。对于问题1，根据假设，靶标和像平面的图形是已知的，利用针孔成像原理建立模型，引入常用的相机模型坐标体系，确定各坐标体系的转换关系，建立线性相机模型，利用透视变换矩阵的相机标定技术的思想给出算法，求取相机模型的内部参数和外部参数。考虑到上述模型是基于靶标点与像点已确定的前提下才能实现的，而实际上我们看到的只是靶标上的若干个圆和其拍摄的像（圆的变形图），不是模型中要求的“点”。因此我们结合光学成像原理和两图形几何位置关系确定若干组特殊靶标点和对应像点。将靶标点与对应像点带入模型，求得相机成像系统的内外参数。此时只要给出靶标上的点，即可求得其像坐标。建模流程图如下： 建立线性相机模型，给出求解相机成像内外参数的算法靶标点、像点坐标带入模型，求得相机成像内外参数在靶标和像图上确定特殊标靶点及其对应像点给定靶标上圆心，得其像坐标4.1.2 问题1模型建立：4.1.2.1坐标系为了定量的描述光学成像过程，计算机一般采用右手定则定义坐标体系，即图像像素坐标系，图像物理坐标系，相机坐标系和世界坐标系，如图1。【1】 图1 图像像素坐标系：在图像上定义直角坐标系，本模型中原点在图像平面的左下角，是以像素为单位的图像坐标系的坐标。 图像物理坐标系：用物理单位表示图像像素位置，定义坐标系，原点定义在相机与图像平面交点，是以毫米为单位的图像坐标系的坐标。 相机坐标系：相机坐标系是固定在相机上的直角坐标系，其原点定义在相机的光心，轴分别平行于图像物理坐标系的轴，与光轴重合，是物点相机坐标。 世界坐标系：坐标系描述环境中任何物体的位置，根据具体情况而定。4.1.2.2常用坐标系变换关系【2】（1）从世界坐标系 到相机坐标系的变换，世界坐标系中的点到相机坐标系的变换可由一个正交变换矩阵和一个平移变换矩阵表示： （1）用齐次坐标和矩阵形式表示为： （2） 其中：是世界坐标系原点在相机坐标系中的坐标，正交旋转矩阵 且矩阵元素满足：正交旋转矩阵中的三个独立变量，加上平移矩阵中的总共有6 个参数决定了相机光轴在世界坐标系中空间位置，因此这六个参数称为相机外部参数。 (2) 从相机坐标系到图像物理坐标系的理想投影变换(为有效焦距)即针孔模型下的理想透视投影变换，有下式成立 （3）同样用齐次坐标和矩阵表示上式为: (4) 将上式的图像坐标系进一步转化为图像坐标系： (5) 齐次坐标表示为： (6)其中,是图像中心（光轴与图像平面的交点）坐标，分别为一个像素在X 与Y 方向上的物理尺寸，分别为X 与Y 方向上的采样频率，即单位长度的像素个数。由此可得物点P与图像像素坐标系中像点P的变换关系： （7）其中， 分别定义为X 和Y 方向的等效焦距。 等4 个参数只与相机内部结构有关，因此称为相机内部参数。（3） 世界坐标系与图像物理坐标系变换关系（共线方程） （8）齐次坐标表示为： (9) 上式就是摄影测量学中最基本的共线方程,说明物点、光心和像点这三点必须在同一条直线上。这是针孔模型或者中心投影的数学表达式。根据共线方程，在相机内部参数确定的条件下，利用若干个已知的物点和相应的像点坐标，就可以求解出相机的六个外部参数，即相机的光心坐标和光轴方位的信息。4.1.3 算法设计利用透视变换矩阵相机标定思想【3】，写成具体算法如下： 将（9）式写为： （10）其中是空间三维点的世界坐标，为相应的图像坐标，为变换矩阵的元素。它包含三个方程，整理消去后，可得到如下两个关于关于的线性方程。 (11)方程（11）描述了三维世界坐标点与相应图像点(u,v,1) 之间的关系。如果已知三维世界坐标和相应的图像坐标，将变换矩阵看作未知数，则共有12 个未知数。对于每一个物体点，都有如上的两个方程，因此，理论上，取6 个物体点，就可以得到12 个方程，从而求得变换矩阵M 的系数。但特殊的是，由于该方程没有常数项，即方程形式为,因此无法直接求解。（1）求解一般可做如下变换令,即求解时令=1，则共有11个未知数，取个目标点可得个方程，是一个超定方程，表示成矩阵形式如下： (12)其中利用最小二乘法求出上述线性方程组的解为：实际上由等式矩阵即可确定物体世界坐标到图像点的映射关系。（2）求解，和将和矩阵分块，其中代表矩阵中第行，代表矩阵中第行的前三个元素组成的行向量,即 = （13） 由矩阵变换可得： （14）化简结果如下： 综上，带入已知参数，可得到投影矩阵M（），同时可以求得相机成像内外参数。（3）在靶标图和像图上确定特殊标点及其对应像点线性相机模型成像原理本质是相同的，因此基于小孔成像的原理，靶标中两圆间公切线所在的切点在经过镜头折射后，必然与像中对应的两个形变后的圆的公切线所在切点是对应的，如下图所示: 由左图（物图）到右图（像图），圆的形状虽然有改变，但是公切线所在切点还是对应不变的。对于物平面上的若干个圆，任取两圆可得四个切点，这样可得若干个切点及其对像平面上的像点，实际取多少切点可根据图形和解题需要具体确定。由此可确定靶标中切点的世界坐标和对应像中切点像素坐标。4.1.4 求解通过上述模型和算法（程序见附录pro1）得到投影矩阵，也就是确定了原模型中相机成像的内外参数。此时，只要给出物平面的圆心，就可直接求得靶标上圆心在像平面上像坐标。4.2问题24.2.1 问题分析根据题目给出的已知条件，首先利用问题1取点方法得到多组靶标点与其对应像点，从而确定靶标点的世界坐标及其像点的像素坐标，再将其带入问题1模型，可以得到相机成像内外参数。最后将图2的圆心的世界坐标带入模型，得到对应像平面上的像坐标。4.2.2 模型的建立（1）选取特殊物点及其对应像点本题中选取的世界坐标系是以物平面的左下角E为原点，EA方向为EB方向为。利用问题1中公切线取点的方法，结合问题本身，可选取的切点及其像平面对应的像点如下（图中仅给出圆两两外切情况，公切线内切选点也可用）：靶标上的取点：像平面上靶标点相对应的像点具体取点的算法和流程图如下（编程见附录program2）：算法：Step1：根据图3靶标的像的每个像素颜色的不同将bmp图转化为0/1编码Step2：运用广度搜索算法扩展节点，建立在一个区域图形的相互关系Step3：运用并查集算法压缩路径，将每个区域图形独立Step4：用一个像点周围四个点来确定该点是否为区域图形的边缘点，如果为边缘点则存入v数组中Step5：对于任意两个区域，各取一个边缘点，确定一条直线y=k*x-bStep6：将这两个区域上所有的边缘点的横坐标带入直线，检测该直线是否与该区域相切，如果该直线与两个区域都相切，则这对边缘点为两个区域的一对切点Step7：对求出的切点进行处理，确定任意两个区域的所有切点注： 任意两个区域应有两对切点，但由于精度原因，两个区域我们求出多对解。将横纵坐标不超过2个像素的点归为一类点，将这一类点求平均值，抛弃距离平均值点10个像素点，再求剩余点横纵坐标的平均值，取整，确定两个区域的两对切点。N将图片转换为0/1编码，0为像的点，1为非像的点读入靶标的像的图片广度搜索，将距离不超过1.5的两个像的点建立关系，认为其在一个区域图象内运用并查集压缩路径，确定每个像点所属区域图像对于任意两个区域图像，分别任取一个边缘点，确定一条直线,将两个区域边缘的点横坐标代入直线，求得纵坐标，与图像纵坐标的比值在一个精度范围内，则该点在直线上一个区域的边缘点只有一个点在直线上Y确定这对点为这两个区域的切点N检测完所有的点对处理数据，得到所有切点，结束Y得到如下特殊切点的像点坐标（像平面左下角为原点，单位是像素）：uvuvuvP1326619P7288301P13364580P2425610P8591296P14323266P3645591P9278232P15604554P4312539P10577233P16549253P5416532P11282579P17676559P6636517P12246262P18617270（2）确定相机内外参数确定特殊点以后，带入问题1中线性相机模型，利用问题1中具体算法可以求得相机成像内外参数：注: 线性方程组中未知数的数目大于实际自由度的数目，因此，线性方程中的未知数不是相互独立的，在有噪声的情况下，解的线性方程中的未知数也许能很好地符合这一组线性方程，但由此分解得到的参数值未必能与实际情况相符，使得精度受到一定限制。（3） 靶标上圆心在像平面上对应的像坐标利用MATLAB编程求解，（程序见pro2）对于靶标上的圆，通过矩阵变换求得其在像平面上的对应点的像素坐标为(以左下角为坐标原点,单位：像素)：ABCDEu320.07419.69635.88581.32283.79v576.51568.78552262.3264.82问题1要求的坐标系下，A,B,C,D,E的相机坐标为（单位：mm）：ABCDEx-50.7751-24.420632.7724918.33862-60.373y50.9285748.883644.44444-32.1958-31.5291z417.1958417.1958417.1958417.1958417.1958利用MATLAB做图，在像平面的圆上标出求得的像点，如下：从图形可以直观的看出，求得的像点在图形（近似看成椭圆）的中心附近，偏差不大，这是符合实际中照相机拍照原理的，由此也体现了模型的合理性和有效性。4.3 问题34.3.1 问题分析 在前两问中，我们建立的模型主要是从光学成像原理方面，利用坐标系的变换关系来找到靶标上圆的像坐标。实际中，由于物距比相距大得多，一般可以近似的认为，圆在经过透镜成像以后变为椭圆，若不考虑畸变因素，可以认为椭圆的中心就是原物平面圆的中心。基于这种思路，在本问题中我们首先提取椭圆边缘点，然后采用最小二乘法提取像平面中椭圆，进而求得椭圆中心坐标，即为靶标圆心的像坐标，并基于此解对问题2进行检测。4.3.2 模型建立（1） 建立椭圆最小二乘法拟合模型：先对采集的像平面的bmp图像转化为0/1编码，通广度搜索和并查集确定椭圆区域边缘点，跟据边缘点数据进行最小二乘拟合。椭圆的平面方程为【4】：=0根据该方程所确定的椭圆的中心坐标为：假设一般形式的椭圆方程如所列，为了避免零解，并将解的任何整数倍都视为对同一椭圆的表述，对参数做一些限制 ，约束条件设为 A=1显然，直接应用上述方程对边缘检测后的离散点进行最小二乘处理，就可以得到方程中的各系数 也即求目标函数( B,C,D,E,F)= 的最小值来确定各系数再由极值原理，欲使f(A ，B，C，D，E)值为最小，必有,由此可得一个线性方程组【5】：然后应用求解线性方程组的算法(如全主元高斯消去法)，结合约束条件，用MATLAB编程求解（程序见pro3），就可以求得方程系数 A，B，C，D，E，F的值，同时求得椭圆圆心（2） 结果检验及精度和稳定性分析Step1:判断区域图像内的每个像点是否在拟合的椭圆内，像素点坐标为，代入，若0,则点在拟合的椭圆外，否则则在椭圆内。Step2:确定像点在拟合椭圆内的个数M，像点个数为N,则该区域的拟合精度为Step3:求出所有区域图像的拟合精度，求平均值，即为此最小二次拟合椭圆的精度利用MATLAB编程求解（程序见pro4），拟合图像为：拟合度为：椭圆12345精度0.97250.90370.98160.97670.9903则平均椭圆拟合度为：=0.96496，由此可以看出拟合度较高，拟合效果很好。（3）利用拟合椭圆对问题2结果的检验由于对相机来说，物距远远大于焦距，圆经小孔成像的投影接近椭圆，我们对成像进行椭圆的最小二乘拟合，以椭圆的中心作为圆心的参考投影点，评价模型结果的精确性以及稳定性。 模型2结果的精确性检验根据椭圆拟合的结果，以椭圆左下角为原点，自左向右为轴，自下向上为轴建立图像素坐标系，得到椭圆中心的坐标ABCDE325.4300424.9500641.0200583.7100287.5200586.9800576.0500557.0700266.4000271.7900定义平均偏差则求得 像素单位若以分辨率，=768作为权重则误差率两种方法得到的圆心坐标基本一致，证明模型结果比较精确。 模型2结果的稳定性分析在问题二的模型中，参数是由切点的世界坐标和图像像素坐标决定的。考察模型结果的稳定性应考察取不同的切点组求得的圆心坐标的差异。在取100组()不同的切点,每组10个切点：第一组：P1,P2,P3,P4,P5,P6,P7,P8,P9,P10第二组：P1,P2,P3,P4,P5,P6,P7,P8,P9,P11第三组：P1,P2,P3,P4,P5,P6,P7,P8,P9,P12第一组ABCDE323.07415.69639.88584.32283.79578.51568.78552.10261.30264.82第二组ABCDE317.09417.52635.88580.32286.79576.52566.30551.98264.30262.81第三组ABCDE319.09419.52628.88586.32289.79582.52585.30558.98269.30266.81 求得误差率分别为： 的取值在0.48,0.78区间内，变化范围不大，证明模型结果有较强的稳定性。4.3.3对误差原因的探讨误差主要是未考虑畸变因素造成的，引入畸变影响因素，可建立非线性模型对其进行修正，进而得到精度更高的解。主要的畸变误差分为三类：径向畸变、偏心畸变和薄棱镜畸变。第一类只产生径向位置的偏差，后两类则既产生径向偏差，又产生切向偏差。（1）径向畸变：光学镜头径向曲率的变化是引起径向变形的主要原因。这种变形会引起图像点沿径向移动，离中心点越远，其变形量越大。正的径向变形量会引起点向远离图像中心的方向移动，其比例系数增大；负的径向变形量会引起点向靠近图像中心的方向移动，其比例系数减小。（2）偏心畸变：由于装配误差，组成光学系统的多个光学镜头的光轴不可能完全共线，从而引起偏心变形，这种变形是由径向变形分量和切向变形分量共同构成。（3）薄棱镜变形：薄棱镜变形是指由光学镜头制造误差和成像敏感阵列制造误差引起的图像变形，这种变形是由径向变形分量和切向变形分量共同构成图像平面理想图像点坐标与实际图像点坐标与畸变误差之和，即：综上，根据实际情况引入不同畸变因素，用非线性相机模型求解，具体可见模型改进。4.4 问题44.4.1 问题分析当两部相机位置固定时确定两部相机的相对位置，就是要确定两部相机间的旋转角度和平移距离。根据其拍摄的图像，利用问题1模型可以分别得到在同一世界坐标系下的物与两部相机的位置关系，然后以物的世界坐标为中间变量，利用旋转矩阵和平移矩阵，对物与相机坐标系进行转换，得到两相机间的相对位置关系。4.4.2 模型建立立体相机模型如图双目立体相机系统不同于单相机的一个主要方面是,除了每个相机的内部参数需要标定之外,还要精确确定两个相机之间的相对位置和姿态, 称为外部参数标定【6】。如图，设P 点在世界坐标系、左镜头坐标系、右镜头坐标系下的坐标分别为，故由问题1模型中从世界坐标系到相机坐标系的坐标变换有 （15） （16）1其中分别为左右摄像机在世界坐标系下的方向矩阵，为左右摄像机在世界坐标系下的平移矩阵则由左摄像机坐标变换可得 将该式带入右坐标系，可得 令 （17） （18）则用表示从相机左坐标系到右相机坐标系的旋转矩阵,表示从左相机坐标系到右相机坐标系的平移向量,则有如下变换关系和精确确定了两个相机的相对位置。4.4.3 模型求解 综合以上模型的分析，可按如下步骤进行求解：第一步：在两个相机照下来的图像上提取标定点，可按问题1模型中选取公切线上切点的方法取点。第二步：将标定点对应的图像像素坐标和世界坐标代入问题1中坐标系变换模型，利用透视变变换矩阵的算法求解，得到各自相机的内外部参数，包括。第三步：将代入到（17）式，求得，确定两部相机（镜头）的相对旋转角度。 将 代入到（18）式，求得 ，确定两部相机的平移向量。4.4.4 模型扩展与应用设空间点P的对应点已经检出。但在实际应用中，需要对两幅摄影图上的对应点进行匹配。对应点的配准(也称对应点的匹配)就是对空间P点在第一张照片中的图像点p1，在第二张照片中找出与p1点对应的P在第二张照片中的图像点p2。对应点的建立还需利用其