数码相机数学建模

1、2008高教社杯全国大学生数学建模竞赛承诺书我们仔细阅读了中国大学生数学建模竞赛的竞赛规则.我们完全明白，在竞赛开始后参赛队员不能以任何方式（包括电话、电子邮件、网上咨询等）与队外的任何人（包括指导教师）研究、讨论与赛题有关的问题。我们知道，抄袭别人的成果是违反竞赛规则的，如果引用别人的成果或其他公开的资料（包括网上查到的资料），必须按照规定的参考文献的表述方式在正文引用处和参考文献中明确列出。我们郑重承诺，严格遵守竞赛规则，以保证竞赛的公正、公平性。如有违反竞赛规则的行为，我们将受到严肃处理。我们参赛选择的题号是（从A/B/C/D中选择一项填写）：A我们的参赛报名号为（如果赛区设置报名号的话

2、）：1324615所属学校（请填写完整的全名）：中国地质大学（武汉）参赛队员（打印并签名）：1.H2.3. 代永力指导教师或指导教师组负责人（打印并签名）：付丽华日期：2008年9月22日赛区评阅编号（由赛区组委会评阅前进行编号）：2008高教社杯全国大学生数学建模竞赛编号专用页赛区评阅编号（由赛区组委会评阅前进行编号）：赛区评阅记录（可供赛区评阅时使用）：评阅人评分备注全国统一编号（由赛区组委会送交全国前编号）：全国评阅编号（由全国组委会评阅前进行编号）：数码相机定位摘要相机定位是计算机视觉领域里从二维图像获取三维信息的基本要求，是完成许多视觉工作必不可少的步骤。题目利用由五个圆组成的靶标模

3、板（图2）对相机进行标定。对于给出的相机拍摄的靶标的像（图3）,本文对其进行了预处理：首先应用MATLAB软件对图3进行去噪，然后利用自行设计的算法（扫描法）提取去噪后图形的边缘（图8）。事实上，标定就是找到空间坐标系和像平面坐标系之间的相互关系。而这是由相机成像的几何模型决定的，具成像参数包括相机内部几何和光学特性（内部参数），以及相机世界坐标系的三维位置和方向（外部参数）。求解相机内外参数是解决问题的关键。本文考虑了理想的情况，即直线的投影仍是直线，以及交点的投影仍然是像的交点。利用图2中A、B、C、D、E两两圆的外公切线的切点与图3中像的外公切线的切点相对应，通过最小二乘法求出该相机的内

4、、外参数。然后，利用相机成像原理，将图2中圆心的坐标Ac、Bc、Cc、Dc、Ec和相机内外参数代入，即可得到靶标上圆的圆心在该相机像平面的像坐标Ac、Bc、Cc、Dc、Ec（见表2）。为了验证该模型，本文利用留一法，用四个圆（共有三组：A、B、D、E;B、C、D、E;A、C、D、E）的切点与对应像点求相机内外参数，剩下的一个圆用于做检验。以第一组为例，具体过程是：根据A、B、D、E算出内外参数，然后求出C圆圆心的像的坐标CJ,并与Ac、Bc所确定的直线的距离来检验方法的精度和稳定性。结果显示（见本文13-14页），本文所用方法精度较高，稳定性较好。双目标定时，将双目系统的坐标系建立在左相机上，

5、把右相机的坐标系相对于其的旋转和平移参数利用最小二乘法求出，从而可以得到两部相机的相对位置。关键词：标定；内、外参数；深度优先搜索；最小二乘法一、问题的提出1 .背景目前，交通管理中使用人工监控的方式已经不适合日益壮大的交通系统，而有效的电子监管（例如数码相机定位）能更有效的管理交通系统运行，于是电子监管逐渐被人们所认识。随着数码相机的普及，数码相机定位在交通监管（电子警察）等方面应用越来越广泛。所谓数码相机定位是指用数码相机摄制物体的相片确定物体表面某些特征点的位置。当前，最常用的定位方法是双目定位，即用两部相机来进行定位。对物体上一个特征点，用两部固定于不同位置的相机摄得物体的像，分别获得

6、该点在两部相机像平面上的坐标。如果知道两部相机精确的相对位置，就可使用几何的方法得到该特征点在固定一部相机的坐标系中的坐标，即确定了特征点的位置。由此可见，对于双目定位，精确地确定两部相机的相对位置就是关键，这一过程称为系统标定。系统标定方法很多，常用的有张正友平面标定方法、孟胡平面圆标定方法、平行圆标定方法1。目前有人提出这样一种方法：在一块平板上画若干个点，同时用这两部相机照相，分别得到这些点在它们像平面上的像点，利用这两组像点的几何关系就可以得到这两部相机的相对位置。然而，无论在物平面或像平面上我们都无法直接得到没有几何尺寸的“点”。实际的做法是在物平面上画若干个圆（称为靶标），它们的圆

7、心就是几何的点。而它们的像一般会变形，如图1所示。问题的关键是，能否从靶标上的这些圆的像中把圆心的像精确地找到，如果能找到，标定就可实现。图1像的变形示意图2 .问题有人设计靶标如下，取1个边长为100mm的正方形，分别以四个顶点（对应为A、C、D、E）为圆心，12mm为半径作圆。以AC边上距离A点30mm处的B为圆心，12mm为半径作圆，如图2所示。图2靶标示意图用一个位置固定的数码相机摄得其像，如图3所示。图3靶标的像研究的问题如下：（1）建立数学模型和算法以确定靶标上圆的圆心在该相机像平面的像坐标，这里坐标系原点取在该相机的光学中心，x-y平面平行于像平面；（2）对由图2、图3分别给出的

8、靶标及其像，计算靶标上圆的圆心在像平面上的像坐标，该相机的像距（即光学中心到像平面的距离）是1577个像素单位（1毫米约为3.78个像素单位），相机分辨率为1024X768;（3）设计一种方法检验模型，并对方法的精度和稳定性进行讨论；（4）建立用此靶标给出两部固定相机相对位置的数学模型和方法。二、符号定义OwXwywZw世界坐标系Oc数码相机透视中心OcXcycZc数码相机坐标系Ou主光轴与像平面的交点OuXuyu坐标轴分别平行于像素坐标系的坐标轴建立的物理坐标系fx0U0F数码相机内参矩阵，F0fyV0,fx，fy分别表示在X001方向和y方向上像点的物理坐标到图像坐标的比例系数，U0,V0

9、是主光轴与像平面交点的图像像素坐标RT数码相机的外参矩阵，R和T分别为模板相对数码相机的r1,1r1,2r1,3t1旋转分量和平移分量，R3收22,3,Tt2r3,1r3,2r3,3t3Ac、Bc、Cc、Dc、Ec圆A、B、C、D、E在靶标平面上的圆心A、Bc、Cc、Dc、Ec圆心Ac、Bc、Cc、Dc、Ec在像平面上对应的点三、问题的分析根据题意，要找到靶标圆的圆心在像平面中的对应点，就应该确定数码相机的内、外参数，然后利用所得到的内、外参数来确定靶标平面上圆心在像平面上对应点的坐标。只要寻找靶标上的点在像平面上的对应点，就可以利用成像原理求出相机的内、外参数。利用圆及其切线对投影面的特点，

10、即两圆及其公切线在投影面上的图形仍保持直线的投影公切于两圆的投影，且它们在靶标平面上的切点对应于投影面上相应的切点，利用这一性质就可以找到靶标平面与投影面上部分相对应的点，然后就可以用所建立的模型来求出数码相机的内、外参数，从而就找到了靶标平面上点与投影面上点的对应关系，再代入数据就可以求出靶标上圆的圆心在像平面上的像坐标。模型稳定性和精度的验证，主要是考查该模型对数据变化这一影响下所得结果的考量。这里，我们将利用留一法，同时用3组数据来进行验证。问题4中两部相机相对位置的确定，主要是两相机像平面坐标系之间的转化，我们将以左相机为基准，将右相机坐标系统一到左相机上，从而得到两者相对位置。四、模

11、型假设1 .假设数码相机成像的原理为针孔相机模型；2 .假设相机本身的配置没有误差；3 .假设数码相机成像产生较小的畸变，对结果影响不大，忽略不计；4 .假设像平面中心与光学中心连线垂直于像平面。五、模型的建立与问题的解（一）模型的建立建立物点（世界坐标）和像点（像坐标）位置对应关系方程。如图4所示：在模板平面上，分别以两条正交的直线为xy轴，以靶标上正方形的中心为原点，并按右手准则建立世界坐标系OwXwyw4,然后以数码相机透视中心Oc为坐标原点，并使其x轴和y轴分别平行于像面坐标系的x轴和y轴,且方向相同，z轴为数码相机镜头的主光轴，按右手准则建立数码相机坐标系OcXcycZc0以像平面上

12、像素坐标值0,0的点为坐标原点O,以水平像元阵列方向为X方向，以垂直像元阵列方向为y方向，建立像素平面坐标系OUV。最后以主光轴与像平面的交点Ou为坐标原点，并使坐标轴分别平行于像素坐标系的坐标轴来建立物理坐标系Ou凡九，模板平面上Pw的理想成像点为R。图4透视投影模型示意图假设模型为理想透视模型根据针孔成像原理，由世界坐标点到理想像素点的齐次变换2,3如下：u s v1F R T ywzw 1其中，s为照相机坐标向像平面坐标转化的比例系数，F为数码相机内参矩阵,RT为数码相机的外参矩阵，R和T分别为模板相对数码相机的旋转分量和平1,11,21,32,12,22,3,Tr3,13,23,3t1

13、t2t3移分量，分别定义如下3,4fx0UoF0fyVo,R001其中，fx,fy分别表示在x方向和y方向上像点的相机坐标到物理坐标的比例系数，U0,V0表示主光轴与像平面交点的图像坐标。由于世界坐标系的xy轴在模板平面上，因此对于模板上的点zw00这样式（1）可重写为:fx000fy0Uori,iV01,1,1消去s即得到以下理想透视数学模型：fxri,iXwri,2ywr1,2r2,2,2tit2t3xwyw13,1xwr3,2ywt3fy2,1xw2,2Ywt23,1Xw3,2ywt3t1U0V0（二）模型的求解1.图像处理题目中给定一幅使用待标定数码相机拍摄的图片（图为1024*768

14、大小的位图。图3中有5个黑色非正规圆斑，分别标号为ABC、黑色非正规圆斑周围为暗灰色，与圆斑颜色相差很大。使用3靶标的像），规格D、E,如图5所示。MATLAB数学软件读入图片后，将图片转化为1024*768的0-1矩阵，图中黑色部分标记为0,非黑色部分标记为1。这样，将图片转化为一个离散化的矩阵，并用Morlet小波进行去噪处理，去噪结果如图5所示。图5黑斑标号图2,边缘提取边缘提取就是提取图片中特征区域的边界，该处指提取图5中黑斑的边界。边缘提取算法很多，MATLAB软件中有边缘提取功能，但效果并不理想。由于黑斑的特点，本文设计新的算法进行边缘提取，称为扫描法（算法实现C程序码见附录一）。

15、扫描法的主要原理是：分两次扫描，第一次按行扫描，第二次按列扫描。按行扫描时，对每一行从左往右遍历，如果第I个点到第J个点全为黑色点，且第I-1、J+1点不为黑色，则标记第I点和第J点，其余点不标记。按列扫描与按行扫描方法一样。按行按列扫描完后，被标记的点即为边缘点。图6单行像素扫描前示意图图6表小一行像素点，黑球表小黑色点，白球表小不为黑色的点。用按行扫描扫描图6,从左往右遍历，可得到两段黑色点段，分别为36、911,标记点3、6、9、11,其余点不标记。保留下来的边缘点如图7所示：图7单行像素扫描后示意图利用扫描法对图5进行处理，得到效果图如图8。图8对图5进行边缘提取的效果图通过上面算法，

16、得到黑斑边缘矩阵后，再对该矩阵使用深度优先搜索将5个黑斑边缘分离出来，存储各自的边缘点，并将5个黑斑边缘标号（见图8）。3.计算机模拟求切点对图8中任意两边缘曲线的外公切线，其他点都在该直线的同一侧，利用该特点，通过计算机模拟来求两边缘曲线的公切线及其切点。具体的模拟算法：对两个边缘曲线，以A、B为例，在曲线A上任取一点a,曲线B上任取一点b,作直线L,然后穷举A和B上的其他点，看是否在该直线的同侧，如果在同侧，既可认为该直线为A、B的公切线，这两点即为切点。例如，根据A边缘点集和B边缘点集，求两边缘曲线的两条外公切线，并求得公切点1,3,5,6。根据上述算法，编写部分C+程序（见附录二），求

17、得的A和B外公切线效果图如图9:图9切线效果图根据已建立的坐标系，将此算法应用于AB、C、D、E圆，可求得像上18个点坐标，这18个点分布如图10所示。图10像坐标系及切点对应点靶标中各圆之间外公切线与圆切点标号，如图11所示。图11靶标坐标系及切点编号圆上切点及其像图10和图11上切点在相应坐标系内的坐标见表1。表1圆上切点及对应像的在相应坐标系中的坐标、坐标切）圆切点在世界坐标系中的坐标（mm）圆上切点的像在像坐标系中的坐标（mm）xyxy1-506248.14815-62.433862-6250-61.11111-53.439153-5038-52.38095-41.0052914-38

18、50-39.41799-50.264555-2062-22.48677-60.317466-20P38-26.98413-39.4179917385024.07407-46.031758506236.77249-55.291019625042.85714-43.9153410503832.01058-35.4497411-50-38-57.6719621.6931212-62-50-70.6349230.6878313-50P-62-61.1111140.21164114-38-50-50.2645531.746031538-509.5238130.952381650-3822.2222222

19、.751321762-5027.2486832.010581850-6218.7830739.682544.数码相机的内、外参数在实际求解的过程中最终的像坐标采用物理坐标形式（x,y）,单位为mm,且fxfxd,fyfyd,其中d为单位像素的物理尺寸（1/3.78）。这样靶标上点与像平面上的点对应关系如下：Xsy1fX0X0fy000011,1r2,13,1r1,2t1r2,2t23,2t3Xwyw1(5)由上式可得：Xsy1一fX1,1fy2,1r3,1一fX1,2fyr2,2r3,2一fXt1fyt：t3IXw2yw1(6)令单应矩阵45A3,4a1a4a2a3a5a6fX1,1fy2,1

20、fX1,2fy2,2fXt1fyt2(7)a7a813,13,2t3则可把参数求解问题转化为先求单应矩阵A,然后从A矩阵分离出待求参数的问题。为求取A阵，变换式(6),消去s/3,4,并分离待求参数，则可得aiXwyw1000XwXywXa2x(8)000Xwyw1XwyywyMya8本文获得了18个切点，将其带入该方程组，共36个方程。已知点的个数远远超过解该线性方程组的需要，这是一个超定方程组。在此使用最小二乘法，解这个超定方程组，求出最优的参数(程序见附录三)。求出单应矩阵A后，即可提取出fx，fy以及外部参数矩阵：,42222a1a2a4a5a1a2a4a51/1/f21 Xf2y22

21、a8a7122a2a5fX2fy22a82a1fX212a42一,2a7fy1,12,13,11,22,23,2t1t2t3fX 0 00 fy 0100 A1(10)(11)旋转矩阵R的第三列可由其正交性求得：这里为向量叉积。根据表1中点的坐标求出相机的内、外参数:5.求靶标圆心像在像坐标系中的坐标将靶标圆心在世界坐标系中的坐标代入方程，使用物理坐标公式，求靶标圆心像在像坐标系中的坐标:fxrm例L(12)3,lXwr3,2ywt3.fyr2,lXwr2,2ywt2r3,lXwr3,2ywt3求出的圆心的像坐标晨bC、cC、dC、eC见表2(根据题目，若以相机光学中6.两相机相对位置的确定假

22、设两部相机的内部参数相同，两部相机的坐标系建立同已提出的模型。为 了确定两部相机的相对位置，将左、右两部相机中左相机固定，即以左相机坐标 系为基准，则右相机相对于左相机的坐标有如下转换关系：XxrV R yr tzzrxryrzr(13)心为原点，则像点的z坐标均为像距417.19576mm):表2靶标圆心及其像点坐标坐标圆编号、靶标圆心坐标(mm)靶标圆心的像坐标(mm)XYXYA-5050-50.460651.7391B-2050-23.801549.6625C505033.897445.1681D-50-50-59.9735-30.9809E50-5019.3974-31.5215其中x

23、l,y,zT,xr,yr,zrT分别是左相机和右相机的像点坐标，r1,1r1,2r1,3Rr2,1r2,2r2,3为变换的旋转分量,r3,1r3,2r3,3t1Tt2为变换的位移分量t3由于在建立坐标系时,相机的像平面上的点位于xy平面上，故z坐标为零，因此转换关系简化为：xlr1,1r1,2Lxrylr2,1r2,2t2yr(14)1r3,1r3,2t31通过模板上圆切点在两相机上的像，可用最小二乘法把九个参数标定。旋转矩阵的第三列可根据其正交性由前两列算出。两坐标系的转换关系求出，即可确定它们的相对位置关系。六、模型的精度和稳定性分析在上一节模型的求解过程中，我们利用A、B、C、D、E5个

24、圆上18个切点与像坐标求出相机内、外参数，为了验证该模型，我们利用留一法，只利用4个圆(A、B、D、E;B、C、D、E;A、C、D、E)来计算另一个圆心像坐标，以A、B、D、E这一组为例，代入圆A、B、D、E上切点坐标及其相对应的投影点坐标，求出数码相机的一组内、外参数，利用这组参数求出圆C的圆心在投影面上对应点的坐标C；,该点应该在圆A和圆B的圆心在投影面上像点，即aC,bC所成直线上面，然而由于不同的参数导致了误差，通过求C；到该直线的距离dc应反映所求参数误差的大小。同理，用B、C、D、E和A、C、D、E求出另外两组参数后，进行类似的处理。三组数据如下：1 .用AB、D、E上切点求相关参

25、数得：Cc的坐标为(32.5312,45.2120),与AB。所成直线距离为0.062355mm2 .用RC、D、E上切点求相关参数得：Ac的坐标为(-51.5430,51.8816),与BC。所成直线距离为0.058033mm3 .用AC、D、E上切点求相关参数：B；的坐标为(-25.3578,52.7422),与AC；所成直线距离为2.949541mm分析数据可知本文提出的模型精确度比较高，前两组误差均不到一个像素，即基本重合，完全可以满足要求。第三组没有考虑B圆相应的切点，误差较大，可见B圆的设置对计算精度有较大影响。这可能与B圆离中心较近，产生畸变较小所致5,6。七、模型的优缺点优点：

26、本模型在提取图像边缘时，采用扫描法，对于噪音极小的简单的图案，提取完全精确，比一般的通用边缘检测算法效果要好；该模型采用最小二乘法求取变换参数，容易实施，且不失其精确性和稳定性；本文采用留一法进行检验，检验结果精度高稳定性好。缺点：模型没有考虑相机物理因素及制作工艺缺陷而导致的成像畸变，所以和实际的成像结果会有一定的偏差；在使用最小二乘法计算相机内、外参数时，由线性方程组所得的线性方程没有经过规格化处理，所求结果会有一定的影响；计算旋转矩阵第三列是直接使用了其正交性，使用叉积计算，而由于计算误差，已得参数可能并不严格满足正交性，这样计算得出的结果并不十分精确(但对以后计算圆心的像坐标没有影响)

27、。参考文献1李伟、吕晓旭、钱晓凡、钟丽云。基于平面模板的摄像机标定方法比较J,激光杂志，第27卷，第2期，54-55,2006年。2毛剑飞、邹细勇、诸静。改进的平面模板两步法标定摄像机J,中国图象图形学报，第9卷，第7期，846852,2004年。3雷彦章、赵慧洁、姜宏志。一种单双目视觉系统结合的三维测量方法J,光学学报，第28卷，第7期，1338-1342,2008年。4吴谨，王俊。基于投影矩阵的摄像机参数标定J,福建电脑，第8期，1527,2008年。5孟晓桥，胡占义。一种新的基于圆环点的摄像机自标定方法J,软件学报，第13卷，第5期，957965,2002年。6J.Heikkila.Ge

28、ometriccameracalibrationusingcircularcontrolpointsJ.IEEETransactionsonPatternAnalysisandMachineIntelligence,October,22(10),10661077,2000.附录一：扫描法算法C代码由图像处理已求出图片的0-1矩阵，设矩阵为a口口，aIJ表示位置为第I行第J列的值，0为黑色点，1为非黑色点。处理完后，图像矩阵保存在s口口中，值为0表示为边缘点，值为-1的点位非边缘点。include#inclueinta1025769;ints1025769;intmain#()inti,j；fr

29、eopen(data.txt,r,stdin);freopen(point.txt,w,stdout);memset(s,-1,sizeof(s);for(i=0;i768;i+)for(j=0;j1024;j+)intt;scanf(%d,&t);a皿i=t;flagji=false;for(intii=0;ii1024;ii+)for(intjj=0;jj768;jj+)if(aiijj=0)sii皿=0;intkk=jj;while(aiikk=0&kk768)kk+;kk-;jj=kk+1;siikk=0;for(intjss=0;jss768;jss+)for(intiss=0;is

30、s1024;iss+)if(aissjss=0)sissjss=0;intkss=iss;while(akssjss=0&kss1024)kss+;kss-;iss=kss+1;skssjss=0;for(inti=0;i1024;i+)for(intj=0;j768;j+)printf(%d,sij);printf(n);return0;附录二：求公切线C+弋码structnodeintx,y;r2500;/r0保存A的所有边缘点，r1保存B所有边缘点intt2;/t0表示A边缘点个数，t1表示B边缘点个数boolCheck(doublea,doubleb,doublec,intt1,int

31、t2,doublew)inti;for(i=0;itt1;i+)doublew1=a*rt1i.y+b*rt1i.x+c;if(w1*w0)returnfalse;for(i=0;itt2;i+)doublew1=a*rt2i.y+b*rt2i.x+c;if(w1*w0)returnfalse;returntrue;voidGetQieDian()boolok=false;for(inti=0;it0;i+)for(intj=0;jt1;j+)if(r0i.x=r1j.x|r0i.y=r1j.y)continue;doublek=(double)(r0i.y-r1j.y)/(r0i.x-r1j.x);doublea=1;doubleb=-k;doublec=k*r0i.x-r0i.y;doublew=c;ok=Check(a,b,c,0,1,-w);/abc为直线方程参数if(ok)point1=r0i;/保存切点坐标point5=r1j;break;if(ok)break;for(inti=0;it0;i+)for(intj=0;jt1;j+)if(r0i.x=r1j.x|r0i.y=r1j.y)contin