数码相机数学建模

2008高教社杯全国大学生数学建模比赛许诺书我们认真阅读了中国大学生数学建模比赛的比赛规则.我们完好理解，在比赛开始后参赛队员不可以以任何方式（包含电话、电子邮件、网上咨询等）与队外的任何人（包含指导教师）研究、议论与赛题有关的问题。我们知道，剽窃他人的成就是违犯比赛规则的,假如引用他人的成就或其余公然的资料（包含网上查到的资料），一定依据规定的参照文件的表述方式在正文引用途和参照文件中明确列出。我们郑重许诺，严格恪守比赛规则，以保证比赛的公正、公正性。若有违犯比赛规则的行为，我们将遇到严肃办理。我们参赛选择的题号是（从A/B/C/D中选择一项填写）：A我们的参赛报名号为（假如赛区设置报名号的话）：1324615所属学校（请填写完好的全名）：中国地质大学（武汉）参赛队员(打印并署名)：1.王飞2.李丁3.代永力指导教师或指导教师组负责人(打印并署名)：付丽华日期：2008年9月22日赛区评阅编号（由赛区组委会评阅行进行编号）：2008高教社杯全国大学生数学建模比赛编号专用页赛区评阅编号（由赛区组委会评阅行进行编号）：赛区评阅记录（可供赛区评阅时使用）：评阅人评分备注全国一致编号（由赛区组委会送交全国前编号）：全国评阅编号（由全国组委会评阅行进行编号）：数码相机定位纲要相机定位是计算机视觉领域里从二维图像获取三维信息的基本要求，是达成很多视觉工作必不行少的步骤。题目利用由五个圆构成的靶标模板(图2)对相机进行标定。关于给出的相机拍摄的靶标的像(图3)，本文对其进行了预办理：第一应用MATLAB软件对图3进行去噪，而后利用自行设计的算法(扫描法)提取去噪后图形的边沿(图8)。事实上，标定就是找到空间坐标系和像平面坐标系之间的互相关系。而这是由相机成像的几何模型决定的，其成像参数包含相机内部几何和光学特征(内部参数)，以及相机世界坐标系的三维地点和方向(外面参数)。求解相机内外参数是解决问题的重点。本文考虑了理想的状况，即直线的投影还是直线，以及交点的投影仍旧是像的交点。利用图2中A、B、C、D、E两两圆的外公切线的切点与图3中像的外公切线的切点相对应，经过最小二乘法求出该相机的内、外参数。而后，利用相机成像原理，将图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圆圆心的像的坐标Cc''，并与Ac'、Bc'所确立的直线的距离来查验方法的精度和稳固性。结果显示（见本文13-14页），本文所用方法精度较高，稳固性较好。双目标准时，将双目系统的坐标系成立在左相机上，把右相机的坐标系有关于其的旋转和平移参数利用最小二乘法求出，进而能够获取两部相机的相对地点。重点词：标定；内、外参数；深度优先搜寻；最小二乘法一、问题的提出1．背景目前，交通管理中使用人工监控的方式已经不合适日趋壮大的交通系统，而有效的电子看管（比如数码相机定位）能更有效的管理交通系统运转，于是电子看管渐渐被人们所认识。跟着数码相机的普及，数码相机定位在交通看管（电子警察）等方面应用愈来愈宽泛。所谓数码相机定位是指用数码相机摄制物体的相片确立物体表面某些特色点的地点。目前，最常用的定位方法是双目定位，即用两部相机来进行定位。对物体上一个特色点，用两部固定于不一样地点的相机摄得物体的像，分别获取该点在两部相机像平面上的坐标。假如知道两部相机精准的相对地点，便可使用几何的方法获取该特色点在固定一部相机的坐标系中的坐标，即确立了特色点的地点。因而可知，关于双目定位，精准地确立两部相机的相对地点就是重点，这一过程称为系统标定。系统标定方法好多，常用的有张正友平面标定方法、孟胡平面圆标定方法、平行圆标定方法[1]。目前有人提出这样一种方法：在一块平板上画若干个点，同时用这两部相机照相，分别获取这些点在它们像平面上的像点，利用这两组像点的几何关系就能够获取这两部相机的相对地点。但是，不论在物平面或像平面上我们都没法直接获取没有几何尺寸的“点”。实质的做法是在物平面上画若干个圆（称为靶标），它们的圆心就是几何的点。而它们的像一般会变形，如图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分别给出的靶标及其像，计算靶标上圆的圆心在像平面上的像坐标,该相机的像距（即光学中心到像平面的距离）是1577个像素单位(1毫米约为3.78个像素单位)，相机分辨率为1024×768；3）设计一种方法查验模型，并对方法的精度和稳固性进行议论；4）成立用此靶标给出两部固定相机相对地点的数学模型和方法。二、符号定义Owxwywzw世界坐标系Oc数码相机透视中心Ocxcyczc数码相机坐标系Ou主光轴与像平面的交点Ouxuyu坐标轴分别平行于像素坐标系的坐标轴成立的物理坐标系fx0u0F数码相机内参矩阵，F0fyv0，fx,fy分别表示在001x方向和y方向上像点的物理坐标到图像坐标的比率系数，u0,v0是主光轴与像平面交点的图像像素坐标RT数码相机的外参矩阵,R和T分别为模板相对数码相机r1,1r1,2r1,3t1的旋转重量和平移重量，Rr2,1r2,2r2,3,Tt2r3,1r3,2r3,3t3Ac、Bc、Cc、Dc、Ec圆A、B、C、D、E在靶标平面上的圆心Ac'、Bc'、Cc'、Dc'、Ec'圆心Ac、Bc、Cc、Dc、Ec在像平面上对应的点三、问题的剖析依据题意，要找到靶标圆的圆心在像平面中的对应点，就应当确立数码相机的内、外参数，而后利用所获取的内、外参数来确立靶标平面上圆心在像平面上对应点的坐标。只需找寻靶标上的点在像平面上的对应点，就能够利用成像原理求出相机的内、外参数。利用圆及其切线对投影面的特色，即两圆及其公切线在投影面上的图形仍保持直线的投影公切于两圆的投影，且它们在靶标平面上的切点对应于投影面上相应的切点，利用这一性质就能够找到靶标平面与投影面上部分相对应的点，而后就能够用所成立的模型来求出数码相机的内、外参数，进而就找到了靶标平面上点与投影面上点的对应关系，再代入数据就能够求出靶标上圆的圆心在像平面上的像坐标。模型稳固性和精度的考证，主假如考察该模型对数据变化这一影响下所得结果的考量。这里，我们将利用留一法，同时用3组数据来进行考证。问题4中两部相机相对地点确实定，主假如两相机像平面坐标系之间的转变，我们将以左相机为基准，将右相机坐标系一致到左相机上，进而获取二者相对地点。四、模型假定1．假定数码相机成像的原理为针孔相机模型；2．假定相机自己的配置没有偏差；3．假定数码相机成像产生较小的畸变，对结果影响不大，忽视不计；4．假定像平面中心与光学中心连线垂直于像平面。五、模型的成立与问题的解（一）模型的成立成立物点（世界坐标）和像点（像坐标）地点对应关系方程。如图4所示：在模板平面上，分别以两条正交的直线为xy轴，以靶标上正方形的中心为原点，并按右手准则成立世界坐标系Owxwywzw，而后以数码相机透视中心Oc为坐标原点，并使其x轴和y轴分别平行于像面坐标系的x轴和y轴，且方向同样，z轴为数码相机镜头的主光轴，按右手准则成立数码相机坐标系Ocxcyczc。以像平面上像素坐标值0,0的点为坐标原点O，以水平像元阵列方向为x方向，以垂直像元阵列方向为y方向，成立像素平面坐标系OUV。最后以主光轴与像平面的交点Ou为坐标原点，并使坐标轴分别平行于像素坐标系的坐标轴来成立物理坐标系Ouxuyu，模板平面上Pw的理想成像点为Pu。图4透视投影模型表示图假定模型为理想透视模型依据针孔成像原理，由世界坐标点到理想像素点的齐次变换2,3以下：uxwywsvFRT(1)1

zw1此中，s为照相机坐标向像平面坐标转变的比率系数，F为数码相机内参矩阵，T为数码相机的外参矩阵，R和T分别为模板相对数码相机的旋转重量和平移重量，分别定义以下3,4：fx0u0r1,1r1,2r1,3t1F0fyv0,Rr2,1r2,2r2,3,Tt2(2)001r3,1r3,2r3,3t3此中，fx,fy分别表示在x方向和y方向上像点的相机坐标到物理坐标的比率系数，u0,v0表示主光轴与像平面交点的图像坐标。因为世界坐标系的xy轴在模板平面上，所以关于模板上的点zw0。这款式(1)可重写为：ufx0u0r1,1r1,2t1xwsv0fyv0r2,1r2,2t2yw(3)1001r3,1r3,2t31消去s即获取以下理想透视数学模型:fxr1,1xwr1,2ywt1u0ur3,2ywt3r3,1xw(4)fyr2,1xwr2,2ywt2v0vr3,1xwr3,2ywt3（二）模型的求解1．图像办理题目中给定一幅使用待标定数码相机拍摄的图片（图3靶标的像），规格为1024*768大小的位图。图3中有5个黑色非正规圆斑，分别标号为A、B、C、D、E，如图5所示。黑色非正规圆斑四周为暗灰色，与圆斑颜色相差很大。使用MATLAB数学软件读入图片后，将图片转变为1024*768的0-1矩阵，图中黑色部分标志为0，非黑色部分标志为1。这样，将图片转变为一个失散化的矩阵，并用Morlet小波进行去噪办理，去噪结果如图5所示。图5黑斑标号图2．边沿提取边沿提取就是提取图片中特色地区的界限，该处指提取图5中黑斑的界限。边沿提取算法好多，MATLAB软件中有边沿提取功能，但成效其实不理想。因为黑斑的特色，本文设计新的算法进行边沿提取，称为扫描法（算法实现C程序码见附录一）。扫描法的主要原理是：分两次扫描，第一次按行扫描，第二次按列扫描。按行扫描时，对每一行从左往右遍历，假如第I个点到第J个点全为黑色点，且第I-1、J+1点不为黑色，则标志第I点和第J点，其余点不标志。按列扫描与按行扫描方法同样。按行按列扫描完后，被标志的点即为边沿点。图6单行像素扫描前表示图图6表示一行像素点，黑球表示黑色点，白球表示不为黑色的点。用按行扫描扫描图6，从左往右遍历，可获取两段黑色点段，分别为36、911，标志点3、6、9、11，其余点不标志。保存下来的边沿点如图7所示：图7单行像素扫描后表示图利用扫描法对图5进行办理，获取成效图如图8。图8对图5进行边沿提取的成效图经过上边算法，获取黑斑边沿矩阵后，再对该矩阵使用深度优先搜寻将5个黑斑边沿分别出来，储存各自的边沿点，并将5个黑斑边沿标号（见图8）。3．计算机模拟求切点对图8中随意两边沿曲线的外公切线，其余点都在该直线的同一侧，利用该特色，经过计算机模拟来求两边沿曲线的公切线及其切点。详细的模拟算法：对两个边沿曲线，以A、B为例，在曲线A上任取一点a，曲线B上任取一点b，作直线L，而后穷举A和B上的其余点，看能否在该直线的同侧，假如在同侧，既可以为该直线为A、B的公切线，这两点即为切点。比如，依据A边沿点集和B边沿点集，求两边沿曲线的两条外公切线，并求得公切点1，3，5，6。依据上述算法，编写部分C++程序（见附录二），求得的A和B外公切线成效图如图9：图9切线成效图依据已成立的坐标系，将此算法应用于A、B、C、D、E圆，可求得像上18个点坐标，这18个点散布如图10所示。图10像坐标系及切点对应点靶标中各圆之间外公切线与圆切点标号，如图11所示。图11靶标坐标系及切点编号圆上切点及其像图10和图11上切点在相应坐标系内的坐标见表1。表1圆上切点及对应像的在相应坐标系中的坐标坐圆切点在世界坐标系中的坐圆上切点的像在像坐标系中的坐标(mm)标(mm)标xyxy切点1-506248.14815-62.433862-6250-61.11111-53.439153-5038-52.38095-41.005294-3850-39.41799-50.264555-2062-22.48677-60.317466-2038-26.98413-39.417997385024.07407-46.031758506236.77249-55.291019625042.85714-43.9153410503832.01058-35.4497411-50-38-57.6719621.6931212-62-50-70.6349230.6878313-50-62-61.1111140.2116414-38-50-50.2645531.746031538-509.5238130.952381650-3822.2222222.751321762-5027.2486832.010581850-6218.7830739.682544．数码相机的内、外参数在实质求解的过程中最后的像坐标采纳物理坐标形式(x,y)，单位为mm，且fx'fxd,fy'fyd，此中d为单位像素的物理尺寸（1/3.78）。这样靶标上点与像平面上的点对应关系以下：xfx'00r1,1r1,2t1xwsy0fy'0r2,1r2,2t2yw(5)1001r3,1r3,2t31由上式可得：xfx'r1,1fx'r1,2fx't1xwsyfy'r2,1fy'r2,2fy't2yw(6)1r3,1r3,2t31令单应矩阵45a1a2a3fx'r1,1fx'r1,2fx't1Ar3,4a4a5a6fy'r2,1fy'r2,2fy't2(7)a7a81r3,1r3,2t3则可把参数求解问题转变为先求单应矩阵A,而后从A矩阵分别出待求参数的问题。为求取A阵,变换式(6),消去s,r3,4，并分别待求参数，则可得a1xwyw1000xwxywxa2x000xwyw1xwyywyM(8)ya8本文获取了18个切点，将其带入该方程组，共36个方程。已知点的个数远远超出解该线性方程组的需要，这是一个超定方程组。在此使用最小二乘法，解这个超定方程组，求出最优的参数（程序见附录三）。求出单应矩阵A后，即可提拿出fx',fy'以及外面参数矩阵：a12a22a42a521/fx'2a82a72(9)a1a2a4a51/fy'2a7a8r3,411(10)a22a522a12a422fx'2fy'2a8fx'2fy'2a7r1,1r1,2t1fx'001r2,1r2,2t20fy'0A(11)r3,1r3,2t3001旋转矩阵R的第三列可由其正交性求得：这里为向量叉积。依据表1中点的坐标求出相机的内、外参数：5．求靶标圆心像在像坐标系中的坐标将靶标圆心在世界坐标系中的坐标代入方程，使用物理坐标公式，求靶标圆心像在像坐标系中的坐标：fx'r1,1xwr1,2ywt1xr3,2ywt3r3,1xw(12)fy'r2,1xwr2,2ywt2yr3,1xwr3,2ywt3求出的圆心的像坐标Ac'、Bc'、Cc'、Dc'、Ec'见表2（依据题目，若以相机光学中心为原点，则像点的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.52156．两相机相对地点确实定假定两部相机的内部参数同样，两部相机的坐标系成立同已提出的模型。为了确立两部相机的相对地点，将左、右两部相机中左相机固定，即以左相机坐标系为基准，则右相机有关于左相机的坐标有以下变换关系：xlxrxrRTyrylRyr(13)T1zrzlzr01此中xl，，zlT，，yr，zrT分别是左相机和右相机的像点坐标，ylxrr1,1r1,2r1,3t1Rr2,1r2,2r2,3为变换的旋转重量，Tt2为变换的位移重量。r3,1r3,2r3,3t3因为在成立坐标系时，相机的像平面上的点位于xy平面上，故z坐标为零，因此变换关系简化为：xlr1,1r1,2t1xrylr2,1r2,2t2yr(14)1r3,1r3,2t31经过模板上圆切点在两相机上的像，可用最小二乘法把九个参数标定。旋转矩阵的第三列可依据其正交性由前两列算出。两坐标系的变换关系求出，即可确立它们的相对地点关系。六、模型的精度和稳固性剖析在上一节模型的求解过程中，我们利用A、B、C、D、E5个圆上18个切点与像坐标求出相机内、外参数，为了考证该模型，我们利用留一法，只利用4个圆（A、B、D、E；B、C、D、E；A、C、D、E）来计算另一个圆心像坐标，以A、B、D、E这一组为例，代入圆A、B、D、E上切点坐标及其相对应的投影点坐标，求出数码相机的一组内、外参数，利用这组参数求出圆C的圆心在投影面上对应点的坐标Cc''，该点应当在圆A和圆B的圆心在投影面上像点，即Ac',Bc'所成直线上边，但是因为不一样的参数致使了偏差，经过求Cc''到该直线的距离dc应反应所求参数偏差的大小。同理，用B、C、D、E和A、C、D、E求出此外两组参数后，进行近似的办理。三组数据以下：1．用A、B、D、E上切点求有关参数得：Cc''的坐标为(32.5312，45.2120)，与Ac'Bc'所成直线距离为0.062355mm。2．用B、C、D、E上切点求有关参数得：Ac''的坐标为(-51.5430，51.8816)，与Bc'Cc'所成直线距离为0.058033mm。3．用A、C、D、E上切点求有关参数：Bc'的坐标为(-25.3578，52.7422)，与Ac'Cc'所成直线距离为2.949541mm。剖析数据可知本文提出的模型精准度比较高，前两组偏差均不到一个像素，即基本重合，完好能够知足要求。第三组没有考虑B圆相应的切点，偏差较大，可见B圆的设置对计算精度有较大影响。这可能与B圆离中心较近，产生畸变较小所致5,6。七、模型的优弊端长处：本模型在提取图像边沿时，采纳扫描法，关于噪音极小的简单的图案，提取完好精准，比一般的通用边沿检测算法成效要好；该模型采纳最小二乘法求取变换参数，简单实行，且不失其精准性和稳固性；本文采纳留一法进行查验，查验结果精度高稳固性好。弊端：模型没有考虑相机物理要素及制作工艺缺点而致使的成像畸变，所以和实质的成像结果会有必定的偏差；在使用最小二乘法计算相机内、外参数时，由线性方程组所得的线性方程没有经过规格化办理，所求结果会有必定的影响；计算旋转矩阵第三列是直接使用了其正交性，使用叉积计算，而因为计算偏差，已得参数可能其实不严格知足正交性，这样计算得出的结果其实不十分精准（但对此后计算圆心的像坐标没有影响）。参照文件李伟、吕晓旭、钱晓凡、钟丽云。鉴于平面模板的摄像机标定方法比较[J]，激光杂志，第27卷，第2期，54－55，2006年。毛剑飞、邹细勇、诸静。改良的平面模板两步法标定摄像机[J]，中国图象图形学报，第9卷，第7期，846－852，2004年。雷彦章、赵慧洁、姜宏志。一种单双目视觉系统联合的三维丈量方法[J]，光学学报，第28卷，第7期，1338-1342,2008年。吴谨，王俊。鉴于投影矩阵的摄像机参数标定[J]，福建电脑，第8期，1527，2008年。孟晓桥,胡占义。一种新的鉴于圆环点的摄像机自标定方法[J]，软件学报，第13卷，第5期，957－965，2002年。[6]J.Heikkila.Geometriccameracalibrationusingcircularcontrolpoints[J].IEEETransactionsonPatternAnalysisandMachineIntelligence,October,22(10)，－，2000.10661077附录附录一：扫描法算法C代码由图像办理已求出图片的0-1矩阵，设矩阵为a[][],a[I][J]表示地点为第I行第J列的值，0为黑色点，1为非黑色点。办理完后，图像矩阵保存在s[][]中，值为0表示为边沿点，值为-1的点位非边沿点。include<stdio.h>#inclue<string.h>inta[1025][769];ints[1025][769];intmain#( ){inti,j;freopen("data.txt","r",stdin);freopen("point.txt","w",stdout);memset(s,-1,sizeof(s));for(i=0;i<768;i++){for(j=0;j<1024;j++){intt;scanf("%d",&t);a[j][i]=t;flag[j][i]=false;}}for(intii=0;ii<1024;ii++){for(intjj=0;jj<768;jj++){if(a[ii][jj]==0){s[ii][jj]=0;intkk=jj;while(a[ii][kk]==0&&kk<768)kk++;kk--;=kk+1;s[ii][kk]=0;}}}for(intjss=0;jss<768;jss++){for(intiss=0;iss<1024;iss++){if(a[iss][jss]==0){s[iss][jss]=0;intkss=iss;while(a[kss][jss]==0&&kss<1024)kss++;kss--;iss=kss+1;s[kss][jss]=0;}}}for(inti=0;i<1024;i++){for(intj=0;j<768;j++)printf("%d",s[i][j]);printf("\n");}return0;}附录二：求公切线C++代码structnode{intx,y;}r[2][500];//r[0]保存A的全部边沿点，r[1]保存B全部边沿点intt[2];//t[0]表示A边沿点个数，t[1]表示B边沿点个数boolCheck(doublea,doubleb,doublec,intt1,intt2,doublew){inti;for(i=0;i<t[t1];i++){doublew1=a*r[t1][i].y+b*r[t1][i].x+c;if(w1*w<0)returnfalse;}for(i=0;i<t[t2];i++){doublew1=a*r[t2][i].y+b*r[t2][i].x+c;if(w1*w<0)returnfalse;}returntrue;}voidGetQieDian( ){boolok=false;for(inti=0;i<t[0];i++){for(intj=0;j<t[1];j++){if(r[0][i].x==r[1][j].xr[1][j].y)continue;doublek=r[1][j].y)/(r[0][i].x-r[1][j].x);doublea=1;doubleb=-k;doublec=k*r[0][i].x-r[0][i].y;doublew=c;

||r[0][i].y(double)(r[0][i].y

==-ok=Check(a,b,c,0,1,-w);//abcif(ok){

为直线方程参数point[1]=r[0][i];

//

保存切点坐标point[5]=r[1][j];break;}}if(ok)break;}for(inti=0;i<t[0];i++){for(intj=0;j<t[1];j++){if(r[0][i].x==r[1][j].xr[1][j].y)continue;doublek=(double)(r[0][i].y(r[0][i].x-r[1][j].x);doubl