高教社杯全国大学生数学建模竞赛数码相机定位

2008高教社杯全国大学生数学建模竞赛——数码相机定位承诺书我们仔细阅读了中国大学生数学建模竞赛的竞赛规则.我们完全明白，在竞赛开始后参赛队员不能以任何方式（包括电话、电子邮件、网上咨询等）与队外的任何人（包括指导教师）研究、讨论与赛题有关的问题。我们知道，抄袭别人的成果是违反竞赛规则的,如果引用别人的成果或其他公开的资料（包括网上查到的资料），必须按照规定的参考文献的表述方式在正文引用处和参考文献中明确列出。我们郑重承诺，严格遵守竞赛规则，以保证竞赛的公正、公平性。如有违反竞赛规则的行为，我们将受到严肃处理。我们参赛选择的题号是（从A/B/C/D中选择一项填写）：A我们的参赛报名号为（如果赛区设置报名号的话）：所属学校（请填写完整的全名）：福建工程学院参赛队员(打印并签名)：1.李林芳指导教师或指导教师组负责人(打印并签名)：林志雄日期：2008年9月21日赛区评阅编号（由赛区组委会评阅前进行编号）：2008高教社杯全国大学生数学建模竞赛编号专用页赛区评阅编号（由赛区组委会评阅前进行编号）：赛区评阅记录（可供赛区评阅时使用）：全国统一编号（由赛区组委会送交全国前编号）：全国评阅编号（由全国组委会评阅前进行编号）：数码相机定位摘要随着数码技术的迅速发展，数码相机定位在交通监管（电子警察）等方面得到了泛的应用，针对本题的数码相机的定位问题，我们建立了一个多目标优化的数学模型I用来确定像的坐标与靶标的坐标之间对应关系，即通过像坐标可以直接求出靶标的坐标。这种方法可以较好的解决由于像坐标存在误差，而引起靶标坐标能否精确计算的问题。我们用此模型，比较准确的还原出靶标上的点。给定靶标上的点，我们可以对应的求出像面上的点，即得到了一个像面上的点与靶标上的点的一一对应的较准确的关系。我们首先要确定出像面上椭圆的中心坐标，因此我们采用了几何方法，根据椭圆任意两点的切线的交点与两切点的中点的连线必经过椭圆中心，所以椭圆中心即可确定下来，由以上的多目标优化模型即可求得靶标的中心，但由于多目标模型运行时间过长，所以我们改进模型为单目标优化模型，经改进以后的模型就可用较短的时间求出靶标的中心。由于本文采用的是一个优化模型，求出的是其近似解，与实际的原坐标位置有一定的偏移，所以我们需检验其精度，采用两种方法检验：1、求出5个点的共面度，结果为mm；2、通过用B点坐标计算出像的坐标，结果与实际相差1点像素，精度相当不错的。对于由两部相机摄的像确定两部相机的相对位置及方向，我们给出一个模型II：利用模型1所求的两组坐标，及坐标系的平移和旋转关系确定。并指出该模型的最小二乘求法，该方法可以较好的处理误差所引起的方程不相容问题。关键词几何变换多目标优化数码成像一、问题的重述数码相机定位数码相机定位在交通警察（电子警察）等方面有广泛的应用。所谓数码相机定位是指用数码相机摄制物体的相片确定物体表面某些特征点的位置。最常用的定位方法是双目定位，即用两部相机来定位。对物体上一个特征点，用两部固定于不同位置的相机摄得物体的像，分别获得该点在两部相机像平面上的坐标。只要知道两部相机精确的相对位置，就可用几何的方法得到该特征点在固定一部相机的坐标系中的坐标，即确定了特征点的位置。于是对双目定位，精确地确定两部相机的相对位置就是关键，这一过程称为系统标定。标定的一种做法是：在一块平板上画若干个点，同时用这两部相机照相，分别得到这些点在它们像平面上的像点，利用这两组像点的几何关系就可以得到这两部相机的相对位置。然而，无论在物平面或像平面上我们都无法直接得到没有几何尺寸的“点”。实际的做法是在物平面上画若干个圆（称为靶标），它们的圆心就是几何的点了。而它们的像一般会变形，如图1所示，所以必须从靶标上的这些圆的像中把圆心的像精确地找到，标定就可实现。图1靶标上圆的像有人设计靶标如下，取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×786；（3）设计一种方法检验你们的模型，并对方法的精度和稳定性进行讨论；（4）建立用此靶标给出两部固定相机相对位置的数学模型和方法。二、问题的分析针孔成像：假设物体表面的反射光都经过一个针孔中心而投影到像平面上，即满足光的直线传播条件。针孔成像主要有光心（投影中心）、成像面和光轴组成。本文中的针孔是指光学中心。本问题研究的是相机的定位问题，数码相机成像原理图如下图（1），又光线穿过光学中心成像，所以凸透镜成像模型可以简化为一个针孔成像模型，其原理图如图（2），由图可看出，景物所在的面为物面，底片所在的面为像面，BO为物距，OD为像距，?OAB与?OCD相似，已知物面任意一点的坐标即可根据相似三角形的原理求得对应像面点的坐标。照相机成像原理图（1）图（2）针孔成像的原理图形三、模型的建立与求解但是由于像坐标的计算上存在误差，这种误差有可能是相机系统造成，也可能是计算标靶圆的像中心时算法所产生的。这种误差可定是存在的，由数码相机摄制物体的相片确定物体表面某些特征点的位置，则必将导致原物象的变形，所以我们就要求的一种误差最小的方法，该问题可转化一个多目标化的最优化问题。根据多目标求误差最小化，该问题可分四步来进行求解：1.求出物面中的点空间坐标(x,y,z)在像面中对应的点(x',y',z')的对应关系。2.由已知物面中圆的圆心求出对应像面中五个椭圆的中心点(x',y',z')。3.该问题相当是第一问的反函数，我们已知了像平面上各点的坐标，反推出靶标的位置。先通过任意3个点的像坐标确定对应靶标的位置，确定一张平面，再验算其他两点是否也在同一张平面上。检验其精度。4.用两部相机拍同一实物，得到两张像，在第三问的基础上，找出两个两部相机的夹角，从而确定相机的相对位置。要确定靶标上圆的圆心在该相机像平面的像坐标，由于我们不能确定相机与实物的相对位置，现以光学中心为原点建立三维空间坐标o-xyz如下图，对于坐标系中的任意一点设其坐标为(x,y,z)，设像距为v，根据成像原理（针孔模型）可以确定像的坐标(x?,y?,z?)，图33．1靶标点与像点坐标的优化模型思路与分析：首先，让我们分析数码相机的成像原理。数码相机是通过以图片文件的格式保存图片的，也就是说，数码相机所产生的图片，不是我们透镜照相机经过折射后所成的照片，这就是我们通常所说的针孔成像，它的照片可以理解为在某一像距处的一种想象的图片。比如：在本文中的照片是357*281（像素），可以理解为在像距为1577/2.73（像素）处的照片，如果用1024*768分辩率计算，则可以理解为像距是1577(像素)的照片。由此可见，数码像机只能用针孔（中心）成像原理的方法来讨论问题，而物理学中的111??是不适合数码像机的成像原理。uvf由于像坐标的计算上存在误差，这种误差有可能是相机系统造成，也可能是计算标靶圆的像中心时算法所产生的。而且是肯定存在的。因此简单的通过几何方法,是无法计算出标靶各圆的圆心坐标。五个标靶圆的圆心是共面的，四条边是垂直的，每条边的长度都是精确的，但是，用像点坐标和几何方法所列出的等式方程由于像坐标的误差而产生互相矛盾，从而使得无解。这种简单几何比例的方法在理论上是精确的，但像坐标存在误差，必将导致所求出的标靶圆的圆心的坐标是不共面的，它们的边也会不垂直，甚至所求出的各边的长度也不是原来的数值。所以，要还原出标靶圆的坐标，我们只能寻求一个误差最小的优化方法。设光学中心与像平面的中心的连线为X轴,建立XYZ直角坐标系,五个标靶圆的圆心坐标设为A(xa,ya,za),B(xb,yb,zb),C(xc,yc,zc),D(xd,yd,zd),E(xe,ye,ze)由于像平面上的横坐标一样，均为v=1577/3.78,五个像的中心坐标设为A'(v,ya,za),B'(v,yb,zb),C'(v,yc,zc),D'(v,yd,zd),E'(v,ye,ze)（1）根据几何学中相似三角形的比例关系，我们很容易的建立下列几个等式关系：222ya(y2'?v2)?y2'(ya?xa)aa222za(y2'?v2)?z2'(ya?xa)aa222xb(y2'?z2')?v2(yb?zb)bb222zb(y2'?v2)?z2'(yb?xb)bbxc2(y2'?z2')?v2(yc2?zc2)cczc2(y2'?v2)?z2'(yc2?xc2)cc222zd(y2'?v2)?z2'(yd?xd)dd222xd(y2'?z2')?v2(yd?zd)dd222222xe(y'?z')?v(ye?ze)eeze2(y2'?v2)?z2'(ye2?xe2)ee（2）根据边的长度，也可以列出下列几个方程：(xa?xb)2?(ya?yb)2?(za?zb)2?900(xa?xc)2?(ya?yc)2?(za?zc)2?10000(xa?xe)2?(ya?ye)2?(za?ze)2?10000(xd?xc)2?(yd?yc)2?(zd?zc)2?10000(xe?xd)2?(ye?yd)2?(ze?zd)2?10000(3)A,B,C三点共线,方向向量比例关系3(xc?xa)?10(xb?xa)3(yc?ya)?10(yb?ya)3(zc?za)?10(zb?za)(4)四条边互相垂直，可列出下列等式方程：(xc?xa)(xe?xa)?(yc?ya)(ye?yc)?(zc?za)(ze?za)?0(xc?xa)(xc?xd)?(yc?ya)(yc?yd)?(zc?za)(zc?zd)?0(xe?xa)(xe?xd)?(ye?ya)(ye?yd)?(ze?za)(ze?zd)?0当然，我们还可以列出其它的方程，但是，这些方程都是等式方程，由于像坐标的误差，使得这些方程之间产生矛盾，从而无法还原出标靶圆的圆心坐标，也就无法确定标靶所在的平面，相机与标靶平面的夹角难以确定。根据上述的分析，我们采用多目标的方法，来确定标靶的圆心坐标。既然误差会产生标靶圆和边的变形，我们不强求上的各等式方程都能同时成立，在上述的方程中选取若干个柔性条件的方程转变为目标函数，其他的方程作为约束条件。这样就可以建立一个多目标的优化数学模型。比如：由于标靶的边在还原时可能不会等于100，我们将他们改为与100的差的平方最小，从而，这些方程就转变为目标函数，其他方程仍为约束条件。目标函数：(xc?xd)(xe?xd)?(yc?yd)(ye?yd)?(zc?zd)(ze?zd)?0((xa?xb)2?(ya?yb)2?(za?zb)2?30)2?min((xa?xc)2?(ya?yc)2?(za?zc)2?100)2?min((xa?xe)2?(ya?ye)2?(za?ze)2?100)2?min((xd?xc)2?(yd?yc)2?(zd?zc)2?100)2?min((xe?xd)2?(ye?yd)2?(ze?zd)2?100)2?min约束条件：上面的其余方程均为约束条件方程。通过求解上述的优化问题，我们可以求出标靶圆的圆心坐标。进一步的具体做法在后面讨论。对于上面的模型，即物面上的坐标与像面上的坐标得近似对应关系，也就是用给定的物面上的坐标，就可对应的求得像面上的坐标，已知像面上的坐标，也可对应的求得物面坐标。3．2确定像平面上各椭圆的中心坐标这是一张分辨率为1024*768的375*281的真彩24的图片，如果我们把像片看作二维矩阵，由于他们的点的相对是已知的，所以不难求出各点的坐标。1、我们在椭圆边缘上随机选取三个点，要求有一定距离并且不共线，在每个点处拟合切线。选取其中两个点，设他们切线的交点为A，两点连线的中点为B，则直线AB过椭圆中心。利用另外一点再求出这样一条直线，两条直线的交点就是椭圆的中心。如下图椭圆O的中心为O点，直线AD、DE、CE与椭圆O相切，切点分别为A、B、C，AB中点为F，BC的中点为G，则DF与CG相交于O点。2、由于相机拍摄的RGB图形，通过用C语言对图形的读取、处理，得到0-1矩阵（0表示的是黑色区域，即是椭圆形的像，程序见附录1），分别用三条两两互不平行的直线在0-1矩阵内部移动，记录每条直线第一次接触0的点（切点）的位置，选取其中任意两个切点，取出两个切点的中点，则两个切点对应直线的交点和中点的连线经过椭圆的中心，利用另外一点再求出这样一条直线，两条直线的交点就是椭圆的中心。图4但是在本题中给出的像的图片的像素是375×281，而实际题目中给出的像素是1024×768，所以Y，Z轴的每一个元素都应乘以2.73倍，经过上述方法的3．3优化模型的求解（靶标圆心坐标的还原）有了五点像的中心坐标，从理论上，用上述的优化模型求出五个标靶圆的圆心坐标是可行的，但是，约束条件方程是等式方程，像坐标存在误差，使得该模型的优化变成不可行。另一方面，我们所建立的优化模型是一个二次的多目标优化模型，需要较大的求解时间。基于这二点的原因，我们必须先对上述的多目标优化进行改进，才能解决求解问题。1、把多目标转成单目标处理原多目标函数都是求最小化，所以，可以将他们求和处理，即将几个多目标的函数取和求最小化。这样多目标问题就转变为单目标问题。2、等式方程太多，容易使得方程不相容，仍然无法求解。因此，必须对约束条件的若干个方程进行适当处理，才能顺利的进行求解。（1）先用优化方法还原出A,C,D,E四点，即先不考虑B点。这样等式约束条件方程将减少几个，可以降低由于误差所引起的不相容性，使得求解能继续下去。（2）从像上观察可看出：A点的像失真较其他的小，我们认为AC?AE，是可以成立的，而其他三个(AC?CD,CD?DE,DE?AE)正交性的约束条件方程很可能不会成立，我们把这三个正交性的约束方程转成目标函数，只要求它们尽可能的正交即可。（3）经过（1）、（2）步骤处理后的优化模型如下：目标函数:MIN(S)?((xa?xc)2?(ya?yc)2?(za?zc)2?100)2?((xa?xe)2?(ya?ye)2?(za?ze)2?100)2?((xd?xc)2?(yd?yc)2?(zd?zc)2?100)2?((xe?xd)2?(ye?yd)2?(ze?zd)2?100)2?((xa?xd)2?(ya?yd)2?(za?zd)2?2?100)2?((xe?xc)2?(ye?yc)2?(ze?zc)2?2?100)2?((xc?xa)(xc?xd)?(yc?ya)(yc?yd)?(zc?za)(zc?zd))2?((xc?xd)(xe?xd)?(yc?yd)(ye?yd)?(zc?zd)(ze?zd))2?((xe?xa)(xe?xd)?(ye?ya)(ye?yd)?(ze?za)(ze?zd))2objectto?(x?x)(x?x)?(y?y)(y?y)?(z?z)(z?z)?0aeacaeccaea?c222?ya(y2'?v2)?y2'(ya?xa)aa?222?za(y2'?v2)?z2'(ya?xa)aa??yc2(y2'?v2)?y2'(yc2?xc2)cc??222222?zc(yc'?v)?za'(yc?xc)?222222?yd(yd'?v)?yd'(yd?xd)?222222?zd(yd'?v)?zd'(yd?xd)?222222?ye(ye'?v)?ye'(ye?xe)?22ze(y'?v2)?z2'(ye2?xe2)?ee?（4）经过（3）的优化模型和LINGO8编程后，可求出A,C,D,E的坐标。B点坐标可以利用A,B,C共线，AB/AC?10/3求出。3、标靶上A,B,C,D,E的坐标如下表：（程序见附录2）目标优化值为：0.40334673．4模型的精度分析有了标靶的五点坐标，就可以通过分析这五点是否共面，来确定模型的计算精度。一般的，这五点是不会共面的，因此，我们先取A,C,D三点做一个平面?1，定义余下二点到平面?1的距离和d,用d的大小作为五点共面度，d越接近与零，它们的共面程度越高，反之，共面程度越差。1、下面是计算共面度的过程：取A,C,D三点，作一个平面?1，法向量?i60.793j8.72419k93.55391n?AC?AD??33.17222.18352125.6168?(?979.45,9854.383,?1266.24)法向量的单位向量e?nn，则E到该平面的距离为d?AE?e?8.0944（），d的值是比较小的，所以共面程度比较高。2、还原B点的像坐标拟合222222??xb(yb'?zb')?v(yb?zb);返回去计算其像的像素坐标，?22，得到2222??zb(yb'?v)?zb'(yb?xb)B??(1577,69.3488,34.2036),而已知第二问中的B点像素坐标为（-1577，69，32），两个坐标非常接近，模型较好。3.5两个相机距离间隔的标定的数学模型每个相机对标靶都有一张照片，每张照片应用上述优化方法可以确定一组标靶圆的圆心坐标。形成两组五对坐标组，不妨记为：{(xa,ya,za),(xb,yb,zb),(xc,yc,zc),(xd,yd,zd),(xe,ye,ze)}{(Xa,Ya,Za),(Xb,Yb,Zb),(Xc,Yc,Zc),(Xd,Yd,Zd),(Xe,Ye,Ze)}这两组表示同样的对应点，只是它们相对于坐标系不同。要标定两个相机的相对位置（间距，方向），可以设变换?X????:?Y???Z????x??f?????A?y???g?，?z??h??????f???其中：A表示3*3坐标旋转矩阵，?g?表示一个相机在另一个相机的坐标系中?h???的坐标。将两组坐标代入变换?中，可以得到15个方程的线性方程组，其中只有12个未知变量（A中有9个，加上f,g,h）。这个方程组是线性的，方程个数比未知个数多，由于坐标在计算上存在误差，可能导致这些方程存在不相容的情况。为了能够计算未知变量，可以利用MATHLAB的最小二乘法求解。可以求出?a11?A??a21?a?31a12a22a32a13??f????a23?中的9个变量，以及?g?的值。?h?a33????令，d?f2?g2?h2，即表示两个相机的间距。?cos?cos?sin??cos?sin?cos?cos?cos??sin?sin?cos?sin?sin?sin????cos?sin???cos???cos?cos??cos?sin?sin?令：A??cos?sin??cos?cos?sin???sin?sin??则?cos?cos??cos?sin?sin??A??cos?sin??cos?cos?sin??sin?sin???cos?cos?sin??cos?sin?cos?cos?cos??sin?sin?cos?sin?sin?sin????cos?sin???cos??=?a11??a21?a?31a12a22a32a13??a23?a33??同样用最小二乘法从矩阵A中可以求出两个相机的在三个方向上的相互交角?,?,?。四、模型评价与改进本题我们通过建立一个多目标优化模型来确定靶标上圆的中心与其像的中心的对应关系，其中利用简单的几何关系组成了优化模型的约束条件，使原本复杂的物理模型更易于理解，推广。文中在求解一个实际问题时，将多目标优化模型改进为单目标优化模型，使程序运行时间缩短，增强了模型的可行性。为了分析结果的精确性，我们先用三点确定一个平面，再用第四个点与平面的距离来定义其精确度。这样模型的优劣就有了一个好的评判标准。文中求得的结果的精确度较高，说明此模型较好。模型为了简化，没有考虑数码相机自身参数的影响，如镜头畸变，这样结果的准确性欠佳，在物体识别应用系统中和视觉精密测量中（物体特征的相对位置必须要精确计算），此模型的结果就不十分适用。在此模型中，我们没有得到靶标上圆的中心与其像的中心的对应关系的确定式子，如果题目给出一些数码相机的自身参数，我们可以通过建立非线性模型来求解，这样可以使数码相机的定位精度更高，增大它的适用范围。五、参考文献《摄像机的标定》/blog/userl/98/index.html六、附录附录一：#include<stdio.h>#include<stdlib.h>#include<math.h>#include<string.h>main(){FILE*fp,*fp1;chars[1],s1[1],b[1024];intx,a[8],kk;longinti,j,k,len,bmp_offset;longwidth,height,llq=-1;fp=fopen("tu2.bmp","rb");fp1=fopen("tu1.txt","w");fseek(fp,10L,SEEK_SET);fread(b,4,1,fp);bmp_offset=b[3]*256*256*256+b[2]*256*256+b[1]*256+b[0];fseek(fp,18L,SEEK_SET);fread(b,4,1,fp);width=b[3]*256*256*256+b[2]*256*256+b[1]*256+b[0];fread(b,4,1,fp);height=b[3]*256*256*256+b[2]*256*256+b[1]*256+b[0];fseek(fp,bmp_offset,SEEK_SET);len=(width+31)/8;len=len/4*4;s1[0]=13;for(i=height-1;i>=0;i--){fseek(fp,bmp_offset+i*len,SEEK_SET);fread(b,len,1,fp);llq=-1;for(j=0;j<len;j++){for(k=0;k<8;k++){a[7-k]=b[j]%2;b[j]=b[j]/2;}for(kk=0;kk<8;kk++){llq++;if(llq>=width){fwrite(s1,1,1,fp1);break;}s[0]=a[kk]+48;fwrite(s,1,1,fp1);}}}fclose(fp);}附录二：!求A,C,D,E的坐标的LINGO程序;model:!五个像点的坐标;A1y=72/3.78;A1z=68/3.78;B1y=69/3.78;B1z=32/3.78;C1y=63/3.78;C1z=-47/3.78;D1y=-43/3.78;D1z=81/3.78;E1y=-43/3.78;E1z=-29/3.78;v=-1577/3.78;init:ax=1200;cx=1000;endinit!以四条边的长度尽量接近100为目标,含义是尽量保持标靶不边形;min=(((ax-Cx)^2+(ay-Cy)^2+(az-Cz)^2)^0.5-100)^2+(((cx-dx)^2+(cy-dy)^2+(cz-dz)^2)^0.5-100)^2+(((ex-dx)^2+(ey-dy)^2+(ez-dz)^2)^0.5-100)^2+(((ex-ax)^2+(ey-ay)^2+(ez-az)^