作业一 数码相机定位.doc

对数码相机双目定位的研究林海10009317 夏丹07109125 向仟飞10009310摘要数码相机定位在交通监管（电子警察）等方面有广泛的应用。本文运用针孔成像模型初步研究了数码相机双目定位的解决方法。首先我们将所得图像转化为二值图，并提取出轮廓。对于问题(1)(2)，我们建立了针孔成像的线性模型，假设像图为靶标的线性映射，因此用椭圆曲线拟合像图轮廓。然后利用线性映射中公切线交点为圆心的性质将公切点选取为特征点，求解圆心的像素坐标。最后将圆心的像素坐标转化到像坐标。最后解得圆心像坐标：A(-49.2130,-84.7389)，B(-23.0004,-83.0871)，C(34.2860,-78.7256)，D(19.3701, -2.2521)，E(-59.5371,-2.4689)。对于问题(3)，我们用所求的公切点与其相对应的靶标上的点，将像素坐标到靶标坐标的映射矩阵用最小二乘法求出来，再将像图椭圆上的“圆心”反映射到靶标上，通过比较其与真实坐标的差距分析精度。对于问题四，我们通过求解两个坐标系的单位正交转换矩阵看出两他们之间的相对位置关系。分别有两部相机得到两张图片，可以分别解得两部相机关于世界坐标系的转换矩阵，联立即得相机之间的转换矩阵，从而得到两部相机之间的相对位置关系。关键词：针孔成像；公切线；坐标转换；系统标定；最小二乘法；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个像素单位)，相机分辨率为1024786；（3） 设计一种方法检验你们的模型，并对方法的精度和稳定性进行讨论；（4） 建立用此靶标给出两部固定相机相对位置的数学模型和方法。2.符号说明u,Op,v 像素坐标系横轴，原点，纵轴x,O,y 像坐标系横轴，原点（相片中心），纵轴Oc,xc,yc,zc 相机坐标系原点（焦点），横轴，纵轴，光轴Ow,xw,yw,zw 世界坐标系原点（靶标中心），横轴，纵轴，光轴f 凸透镜焦距dx 表示每一个像素在x轴方向上的物理尺寸dy 表示每一个像素在x轴方向上的物理尺寸M 像平面到物平面的投影矩阵R 旋转矩阵T 三位平移向量3问题分析3.1相机成像原理本问题研究的是数码相机的定位问题，因此我们首先要弄清数码相机的成像原理。由数码相机的结构可以看出，数码相机的成像原理涉及凹凸透镜的光学原理及相机本身的结构，满足凸透镜成像公式1/f=1/u+1/v。但由于一般拍照时物距都远大于像距（1577px）和焦距，且光线穿过光学中心成像，所以数码相机的凸透镜成像模型可以简化为一个针孔成像模型，其原理图如图4。由图可看出，景物所在的面为物面，底片所在的面为像面，BO为物距，OD为像距， 与相似，已知物面任意一点的坐标即可根据相似三角形的原理求得对应像面点的坐标。图4 针孔成像模型3.2求圆心根据针孔成像模型，物体与像的对应关系应该是线性唯一的，这就使得点线面的一些几何性质在此映射下保持不变。直线的投影还是直线，圆的投影是椭圆，切线的投影还是切线，交线的投影还是交线且焦点不变。根据这些性质，我们可以联想到两圆之间外公切点的连线交点即是圆心，也就是说，投影过去后，两椭圆外公切线的切点连线交点即对应着圆心。而已知方程的曲线之间的公切线是很容易联立方程得到的，因此我们可以想办法拟合出轮廓的椭圆曲线方程，进而求得圆心。得到圆心的像素坐标之后，还需经过一系列转换矩阵转换到题目所要求的像坐标系。3.3双目标定几何法可以通过两者所拍的照片，设法求出两者分别关于靶标的距离与角度，解三角形知道两个相机的相对位置。而最数学的方式就是求出两步相机的坐标系之间的转换关系，得到其旋转矩阵和平移矩阵，从而两者之间的距离和角度。通过照片设法求出两者分别关于某一坐标系（如靶标）的转换关系，联立解得两者之间的转换关系。无论如何，问题都被归一到了通过照片求靶标的相对位置或像素坐标系与世界坐标系的转换矩阵。因为呈线性关系，所以可以选取几个合适的特征点（如公切点）代入解得转换矩阵即可。4基本假设1. 假设各个方向上光强度较好，不影响图像成像质量。2. 假设数码相机中图像平面与光轴垂直。3. 假设靶标像的中心恰好在光轴上。4. 以相机焦点为原点，平面平行于像平面建立坐标系。 5. 忽略透镜尺寸及其产生的图像畸变，即物体与图像呈唯一的线性映射关系。6. 部分假设在题中给出。 5模型建立5.1轮廓的提取首先要对图像做一些数值化处理。直接将图3用imread导入MATLAB软件，用edge指令提取轮廓。运行程序见附录一。运行结果见图5。图5 轮廓提取图5.2轮廓的曲线拟合由于我们用针孔成像模型模拟真实的数码相机，并且假设在针孔成像的过程中，物体与像之间呈线性映射关系，不产生畸变。而圆在空间某平面上的投影只能是椭圆，因此我们准备用椭圆曲线来拟合像图轮廓。假设椭圆曲线方程为x2+bxy+cy2+dx+ey+f=0，我们在曲线上选取几个点(x1,y1) ，(x2,y2)， (xn,yn)，那么可以写成矩阵形式(5.2.1)= (5.2.1)对此超定方程，我们选取了足够多的点，用最小二乘法解出方程中的系数，以尽可能减小误差。5.3公切线法求圆心在线性投影过程中，点线面之间的几何关系保持不变，相切依旧相切，相交依旧相交且交点不变。因此在图6中，圆与圆之间的公切线的投影成为椭圆与椭圆之间的公切线，其交点即对应着圆心。首先求公切线。已知两椭圆方程分别为x12+b1x1y1+c1y12+d1x1+e1y1+f1=0，其中一条公切线为，切点分别为(x1,y1)， (x2,y2)。将y=Ax+B代入椭圆方程，并令=0。图6 公切点示意图即有方程组(5.3.1)1=(b1B+2c1AB+d1+e1A)2-4(c1A2+1+b1A)(c1B2+e1B+f1)=0 2=(b2B+2c2AB+d2+e2A)2-4(c2A2+1+b2A)(c2B2+e2B+f2)=0 (5.3.1)解出A、B后代入方程组(5.3.2)后解出切点坐标。 (5.3.2)y1=Ax1+B，y2=A x2+B对椭圆A，设“圆心”A(x,y)，四个公切点A1(x1,y1)，A2(x2,y2)，A3(x3,y3)，A4(x4,y4)，其中A1与A2相对，A3与A4相对，则联立方程组(5.3.3) (5.3.3)解出“圆心”的像素坐标。5.4像素坐标到像坐标的转换像素坐标，就是指即像素的行、列所标识的坐标，如图7所示，用uOpv坐标表示，单位为像素。像坐标，就是坐标系原点取在照片中心的物理坐标。如图7所示，用xOy坐标表示，单位为毫米。设O点在uOpv坐标系中坐标为(u0,v0)。Opuyx(u0,v0)vO图7 像素坐标与像坐标示意图假设靶标像的中心恰好在光轴上且图片大小为1024*768像素，所以O(u0,v0)的像素坐标为(512,384)。如图所示，若用dx和dy分别表示每一个像素在x轴与y轴方向上的物理尺寸，那么像素坐标(u,v)与像坐标(x,y,z)有如下比例关系式(5.4.1) (5.4.1)即 (5.4.2)逆矩阵为 (5.4.3)这就是像素坐标到像坐标的转换。5.5相机坐标到世界坐标的转换相机坐标，就是坐标原点取在相机光学中心，x-y平面平行于像素平面。用Oc-xcyczc表示，单位为毫米。世界坐标，就是依赖于被摄物体（本题为靶标）而不依赖相机的坐标系。本题不妨以靶标中心为原点，AC向为横轴，AE向为纵轴，光轴为z轴，用Ow-xwywzw表示，单位为毫米。由于xc-yc平面平行于像素平面且该相机的像距（即焦点到像平面的距离）是1577个像素，所以在相机坐标系下，像的zc=1577px=1577/3.78mm=417.2mm。设空间中某一点P在世界坐标系下的坐标为(xw,yw,zw)，在数码相机坐标系下的坐标为(xc,yc,zc)，则它们存在关系式(5.5.1)： (5.5.1)其中R为三阶单位正交矩阵，其标准形式为T为三位平移向量(t1,t2,t3)T。也即是 (5.5.2)其中O=(0,0,0) T。这就是任意两个正交坐标系之间的转换。5.6针孔成像模型像坐标到世界坐标的转换假设数码相机中图像平面与光轴垂直，且靶标像的中心恰好在光轴上，那么用针孔模型可以近似表示空间任何一点P在图像上的成像位置，即任何点P在图像上的投影位置p，为光心O与P点的连线OP与图像平面的交点，这种关系也称为中心射影或透视投影。那么如图8所示，像坐标与相机坐标有几何关系(5.6.1)： (5.6.1)也即是 (5.6.2)图8 像坐标与相机坐标的几何关系示意图代入式(5.4.2)和(5.5.2)得(5.6.3)M3*4即像素坐标到世界坐标的转换矩阵。解出这个矩阵，即可得到物体与像的对应关系。5.7 系统标定(1)几何法如图9所示，在不同的相机坐标系下，某点P具有不同的坐标，相对于O1与O2的距离和角度也不同。通过解三角形，即可得到O1与O2的位置关系。图9 双目定标几何示意图(2)代数法首先解出相机1坐标系与相机2坐标系到世界坐标系的转换关系分别为 ， 则有 (5.7.1)由于R1，R2都是单位正交矩阵，所以R1R2-1也是单位正交矩阵，也能写成代入多组数据，用最小二乘法即可解得两个相机的在三个方向上的相互交角以及两个相机之间的距离|。6模型求解6.1像轮廓的椭圆曲线拟合在椭圆A所得的轮廓上分别取三十个尽可能离散的点Pi(xi,yi)(i=1,230)，建立方程xi2+bxiyi+cyi2+dxi+eyi+h=0,即有 (6.1.1)用最小二乘法解得（MATLAB程序见附录二）椭圆的拟合方程为椭圆A:x2+0.0892*x*y+0.9752*y2-244.6511*x-147.2674*y+19518.2075=0; 椭圆B: f=x2+0.0812*x*y+1.0041*y2-316.8527*x-157.9675*y+30135.1418=0;椭圆C: f=x2+0.2717*x*y+0.9323*y2-491.0356*x-209.9722*y+65720.9098=0;椭圆D: f=x2+0.3645*x*y+1.0463*y2-495.6962*x-463.6696*y+95673.1015=0;椭圆E: f=x2+0.2309*x*y+1.2147*y2-252.4949*x-471.2171*y+56409.1379=0;将拟合曲线画出如图10，MATLAB程序见附录三。图10 椭圆曲线拟合图（蓝色为拟合线）由图中可以看出，椭圆曲线拟合程度较高。6.2椭圆之间公切点的求解假设椭圆1 ：x12+b1x1y1+c1y12+d1x1+e1y1+h1=0与椭圆2 ：x22+b2x2y2+c2y22+d2x2+e2y2+h2=0的其中一条外公切线为y=Ax+B，公切点分别为(x1,y1) ，(x2,y2)，则对联立方程组(5.3.1)可用最小二乘法（MATLAB程序见附录四）解得：椭圆A,B的两条外公切线：y= 0.0513x+ 79.0824y= 0.0766x+ 45.7406同理可得：椭圆A,E的两条外公切线：y= -8.5866x+ 965.4573y= -7.5254x+ 1079.5395椭圆B,D的两条外公切线：y= 1.8436x+ -183.0276y= 1.9669x+ -266.3146椭圆B,C的两条外公切线：y= 0.0658x+ 76.8534y= 0.0854x+ 44.3537椭圆C,D的两条外公切线：y= -5.3903x+ 1273.3998y= -4.9070x+ 1295.5962椭圆D,E的两条外公切线：y= -0.0034x+ 197.1047y= 0.0082x+ 170.3364同时求出切点坐标：椭圆公切线切点坐标AAB(117.7878,85.1225),(120.9900,55.0103)AE(104.4058,68.9586),( 133.9654,71.3995)BAB(154.1450, 86.9875),(157.1974,57.7843)BD(142.5287,79.7374),(168.6194, 65.3388)CBC(232.1853,92.1204),(237.7911,64.6693)CD(221.8287, 77.6718),( 247.8585,79.3614)DBD(203.0654,191.3426),(225.5880, 177.3889)DE(212.0984,196.3862),(216.5801, 172.1213)EAE(91.0113,183.9723),(118.9778,184.1872)DE(103.5842,196.7537),(106.6070, 171.2150)表一 公切点坐标6.3圆心的像素坐标的求解通过求公切点连线交点求得椭圆的“圆心”。对椭圆A，设圆心A的投影“圆心”坐标为A(x,y)，四个公切点A1(x1,y1)，A2(x2,y2)，A3(x3,y3)，A4(x4,y4)，其中A1与A2相对，A3与A4相对，则联立方程组(5.3.3)可解得椭圆A的“圆心”像素坐标, 同理可得其余椭圆的像素坐标圆心像素坐标A(325.9747, 191.8576)B(425.0586, 198.1068)C(641.6010,214.6079)D(585.2190,503.9348)E(286.9496, 503.1146)表二 圆心的像素坐标对应到图中也能看出大致位置是对的。6.4像素坐标到像坐标的转换用求解得：圆心像坐标A(-49.2130,-50.8313)B(-23.0004,-49.1781 )C(34.2860, -44.8127)D(19.3701, 31.7288 )E(-59.5371,31.5118)表三 圆心的像坐标6.5 M变换矩阵M变换矩阵，就是像素坐标到世界坐标的转换矩阵。 (6.5.1)因为本题不妨以靶标中心为原点，AC向为横轴，AE向为纵轴，光轴为z轴，所以zw=0，因此化为 (6.5.2)取多个点(ui,vi)，对应(xi,yi),则有 (6.5.3)我们选取了上文中所求得的公切点代入，用最小二乘法（MATLAB程序见附录五）求得即7精度分析及检验7.1误差原因分析(1)畸变产生的误差我们采用的是针孔成像的线性模型，并没有考虑透镜导致的畸变。实际的相机拍的照都是有畸变的，在广角镜头等更加明显。此外，空气折射等因素也会导致非线性，是的结果更加复杂。畸变对精度会产生一定影响。(2)连续图形离散化产生的误差现实中物体的图像一般都是连续的，但在数码相机照射后，相机将连续的图形离散化，转化成一组点阵来存储。在连续图形离散化的过程中，图形边缘的信息难免会有损失，这就带来了一定的误差。(3)过渡色不确定带来的误差实际拍摄过程中，往往会因为对焦没对准，图片上出现一些过渡色，数码相机处理颜色的过程是数字化的、不连续的，以数据形式存储。在确定颜色的过程中，颜色本身的信息会丢失，这也带来了一定误差。(4)模型误差以上原因都是客观条件上产生的，而模型误差则存在于模型本身。比如我们做了一些假设：各个方向上光强度较好，不影响图像成像质量；数码相机中图像平面与光轴垂直；靶标像的中心恰好在光轴上这些假设在实际情况中都是不可能出现的。另外，轮廓提取过程中图像边缘信息的损失，选取特征点对原始轮廓的进一步损失，都造成了一定的误差。(5)图片误差由于1024*768的图片保存格式不对，无法导入MATLAB，所以我们无奈采用了375*281的图片，这也有可能导致较大误差。7.2精度检验用公切点求出像素坐标到世界坐标的转换矩阵M后，再将圆心的像素坐标反映射到世界坐标系，与真实坐标相对比，这就是我们的检验思路。解得A(-51.3252,50.6869)B(-19.0716, 50.0633)C(51.5604,47.7670)D(47.9971,-48.0026)E(-48.1864, -52.0302)距离实际圆心的误差均值为2.14mm相对于圆的半径12mm，这个数值还是可以的。毕竟375*281图片像素相对于1024*768的图片清晰度有点低，对1024*768的图片误差可能会小很多。8稳定性分析及检验稳定性分析就是分析模型对参数、初值的依赖敏感性。为了检验方法的稳定性，我们采用给圆心像素坐标增减一个微小量结果再与原始结果比较的微扰法探究其对结果的影响。为了简便起见，我们只考虑增减一个像素的情况。MATLAB程序见附录六，运行结果如下表四：原始结果像素坐标均增加一个像素像素坐标均减少一个像素A(-51.3252,50.6869)(-50.3070,49.8292)(-52.3434，51.5446)B(-19.0716, 50.0633)(-18.0534,49.2056)(-20,50.9209)C(51.5604,47.7670)(52.8872,46.9094)(50.5422,48.6247)D(47.9971,-48.0026)(49.0153,-48.8603)(46.9789,-47.1449)E(-48.1864,-52.0302)(-47.1682,-52.8878)(-49.2046,-51.1725)表四 稳定性分析五个圆在像素横纵坐标增加一个像素时对应的圆心与开始时所求圆心距离的均值为 0.2302，标准差为0.5251。五个圆在像素横纵坐标减少一个像素时对应的圆心与开始时所求圆心距离的均值为0.9017，标准差为0.0302由此可见，像素坐标改变一点点，不会显著改变结果，偏差也不超过1mm，所以对该模型，稳定性验证合理。9模型评价模型优点：1.用靶标像的四条公切线的交点对应于靶标上的圆的公切线的交点，由此确定靶标平面在相机坐标中的位置，用模拟退火算法和最小二乘法求解，误差较小；2.精度检验的模型中充分考虑了靶标像边界的信息，从整体上把握了圆的作用；3.稳定性检验中，将靶标像上的离散点分别取不同数目，其异概率性和随机性增强了检验的可靠性和说服力。模型缺点：1.公切线法只适用有多个圆的情况，方法的通用性不好。2.本文模型没有考虑相机拍摄时的畸变情况，使得结果存在一定的误差。3. 模型的验证过程无实际实验数据的支持，不能更精确的对模型优化改进，以及在实际使用中给予详细指导性建议。 参考文献：【1】Frank R.Goprdano, William P.Fox，数学建模，机械工业出版社，第四版，2010年。【2】杨淑莹 模式识别与智能计算Matlab技术实现，北京：电子工业出版社，2008年。【3】刘恩惠，数码相机入门，成像原理关键词，/view/2b7ffd1f650e52ea55189845.html，2011.06.26。附录附录一:轮廓提取I = imread(tuxiang.jpg);BW1 = edge(I,prewitt);imshow(BW1)flag=find(BW1=1);xx=floor(flag./1024);% yy=mod(flag,1024);pixval(on); axis on;附录二：求解椭圆拟合方程format long;A=104.481,69.4491;118.01,54.8787;131.54,61.1231;126.336,82.9787;132.581,75.6935;104.481,68.4083;109.684,82.4583;109.164,80.8972;127.377,82.4583;120.092,84.5398;133.621,65.2861;129.458,58.0009;124.255,55.3991;117.491,55.3991;123.734,84.0194;133.621,72.0509;133.621,67.3676;111.766,59.9602;116.969,85.0602;106.562,78.8157;103.96,70.4898;134.142,67.888;130.499,79.8565;126.336,83.4991; 110.205,82.4583;105.001,74.6528;104.481,72.0509;103.96,66.3269;108.123,59.562;112.286,56.9602B=zeros(30,5);for i=1:30 B(i,1)=A(i,2)*A(i,2); B(i,2)=A(i,1)*A(i,2); B(i,3)=A(i,1); B(i,4)=A(i,2); B(i,5)=1;endBC=zeros(30,1);for i=1:30 C(i,1)=-A(i,1)*A(i,1);endC;D=B*B;E=B*C;X=DE 附录三：作出椭圆拟合图像I = imread(tuxiang.jpg);BW1 = edge(I,prewitt);imshow(BW1)flag=find(BW1=1);xx=floor(flag./1024); yy=mod(flag,1024);pixval(on); axis on;hold on;syms x yf=x2+0.0892*x*y+0.9752*y2-244.6511*x-147.2674*y+19518.2075;ezplot(f,0,400,0,300);hold on;f=x2+0.0812*x*y+1.0041*y2-316.8527*x-157.9675*y+30135.1418;ezplot(f,0,400,0,300);hold on;f=x2+0.2717*x*y+0.9323*y2-491.0356*x-209.9722*y+65720.9098;ezplot(f,0,400,0,300);hold on;f=x2+0.2309*x*y+1.2147*y2-252.4949*x-471.2171*y+56409.1379;ezplot(f,0,400,0,300);hold on;f=x2+0.3645*x*y+1.0463*y2-495.6962*x-463.6696*y+95673.1015;ezplot(f,0,400,0,300);附录四：求解公切线方程syms x y;%A,Bx,y=solve(0.0892*y+2*0.9752*x*y-244.6511-147.2674*x)2-4*(0.9752*x*x+1+0.0892*x)*(0.9752*y2-147.2674*y+19518.2075)=0,(0.0812*y+2*1.0041*x*y-316.8527-157.9675*x)2-4*(1.0041*x*x+1+0.0812*x)*(1.0041*y2-157.9675*y+30135.1418)=0)syms x y;%A,Ex,y=solve(0.0892*y+2*0.9752*x*y-244.6511-147.2674*x)2-4*(0.9752*x*x+1+0.0892*x)*(0.9752*y2-147.2674*y+19518.2075)=0,(0.2309*y+2*1.2147*x*y-252.4949-471.2171*x)2-4*(1.2147*x*x+1+0.2309*x)*(1.2147*y2-471.2171*y+56409.1379)=0)syms x y;%B,Dx,y=solve(0.3645*y+2*1.0463*x*y-495.696-463.6696*x)2-4*(1.0463*x*x+1+0.3645*x)*(1.0463*y2-463.6696*y+95673.1015)=0,(0.0812*y+2*1.0041*x*y-316.8527-157.9675*x)2-4*(1.0041*x*x+1+0.0812*x)*(1.0041*y2-157.9675*y+30135.1418)=0)syms x y;%B,Cx,y=solve(0.2717*y+2*0.9323*x*y-491.0356-209.9722*x)2-4*(0.9323*x*x+1+0.2717*x)*(0.9323