数学建模：数码相机定位.doc
【精品文档】如有侵权，请联系网站删除，仅供学习与交流数学建模：数码相机定位.精品文档.高教社杯全国大学生数学建模竞赛承 诺 书我们仔细阅读了中国大学生数学建模竞赛的竞赛规则.我们完全明白，在竞赛开始后参赛队员不能以任何方式（包括电话、电子邮件、网上咨询等）与队外的任何人（包括指导教师）研究、讨论与赛题有关的问题。我们知道，抄袭别人的成果是违反竞赛规则的, 如果引用别人的成果或其他公开的资料（包括网上查到的资料），必须按照规定的参考文献的表述方式在正文引用处和参考文献中明确列出。我们郑重承诺，严格遵守竞赛规则，以保证竞赛的公正、公平性。如有违反竞赛规则的行为，我们将受到严肃处理。我们参赛选择的题号是（从A/B/C/D中选择一项填写）： 我们的参赛报名号为（如果赛区设置报名号的话）： 所属学校（请填写完整的全名）： 参赛队员 (打印并签名) ：1. 2. 3. 指导教师或指导教师组负责人 (打印并签名)： 日期： 年 月 日赛区评阅编号（由赛区组委会评阅前进行编号）：高教社杯全国大学生数学建模竞赛编 号 专 用 页赛区评阅编号（由赛区组委会评阅前进行编号）：赛区评阅记录（可供赛区评阅时使用）：评阅人评分备注全国统一编号（由赛区组委会送交全国前编号）：全国评阅编号（由全国组委会评阅前进行编号）：数码相机定位摘要柯达于1975年开发世界第一部数码相机。由此，数码照相机便家喻户晓起来。数码相机定位在交通监管（电子警察）等方面有广泛的应用。所谓数码相机定位是指用数码相机摄制物体的相片确定物体表面某些特征点的位置。最常用的定位方法是双目定位，即用两部相机来定位。对物体上一个特征点，用两部固定于不同位置的相机摄得物体的像，分别获得该点在两部相机像平面上的坐标。只要知道两部相机精确的相对位置，就可用几何的方法得到该特征点在固定一部相机的坐标系中的坐标，即确定了特征点的位置。于是对双目定位，精确地确定两部相机的相对位置就是关键，这一过程称为系统标定。 标定的一种做法是：在一块平板上画若干个点， 同时用这两部相机照相，分别得到这些点在它们像平面上的像点，利用这两组像点的几何关系就可以得到这两部相机的相对位置。然而，无论在物平面或像平面上我们都无法直接得到没有几何尺寸的“点”。实际的做法是在物平面上画若干个圆（称为靶标），它们的圆心就是几何的点了。而它们的像一般会变形，所以必须从靶标上的这些圆的像中把圆心的像精确地找到，标定就可实现。关键词：针孔成像，坐标变换，图像处理，相机镜头畸变，双目定位一、 问题的提出与重述1.数码相机监视系统是由：景点（scene）方位、相机参数以及方位、成像参数部分组成的系统，系统的标定就是要确定系统各部分的相互关系（主要是几何、数量关系），系统的参数主要有两部分：1)相机的内参数：用于描述相机本身特定属性的参数以及在空间中定位的参数，2)相机的外参数：是指相机相对与世界坐标系（用于公共参考的坐标系）的位置关系，主要由相机坐标系到世界坐标系的变换（平移、旋转）来描述。确定了相机的内参数和外参数，系统就标定成功。 2.图示 u -景点在像平面中的像；uOC -主光轴（相机坐标系Z轴）上的点在像平面的像点（殴氏坐标）；u0a -主光轴（相机坐标系Z轴）上的点在像平面的像点（仿射坐标）； 3. 坐标系统 OwXwYwZw - 世界坐标系 OcXcYcZc - 相机坐标系 OiXiYiZi - 像欧氏坐标系 OaXaYaZa - 像仿射坐标系 注：世界坐标系是系统的一个客观的参考系； 相机坐标系原点在相机光心（焦点）；像的欧氏坐标系与相机坐标系的关系是：Z轴平行且同向，X-Y面平行；像仿射坐标系与像欧氏坐标系关系密切，Z轴，Y轴平行，X轴有个倾斜，主要考虑是，像素的方快长和宽可以不等，而且，视觉效果上可能会出现倾斜情况。 4.射影几何简介主要介绍如何通过2D图像信息实现3D世界的自动测量，这里的测量主要指，3D中点的空间坐标，以及通过2D图像两点位置关系测量三维距离信息，这里限于针孔模型（thin lens）也称中心投影（central projection）针孔模型的图像信息中，3D中的平行线不在保持平行。射影空间的概念考虑不包括坐标原点的n+1维空间，Rn+1-0,0,.0，定义一个等价关系，x1,x2,.,xnT等价于x1',x2',.,xn'T当且仅当存在非零数值t,x1,x2,.,xn,xn+1T=t*x1',x2',.,xn',xn+1'T,射影空间Pn等于Rn+1-0,0,.0 关于此等价关系的商空间,射影空间中的点称为齐性类,射影空间中的点的坐标通常用齐性坐标表示为，x*=x1,x2,.,xn,1T,最后一个坐标为1,事实上,通过原点的任意直线上的点(原点除外)属于同一个等价类. 于是,对相机坐标系的过原点的任意射线上的点,是等价类, 因为他们的像点相同. 于是, 射影空间Pn可以和Rn建立起一一对应, ，x1,x2,.,xn,1T-x1,x2,.,xnT 注：在这样的表示下的好处是，坐标变换中的平移、旋转的表达形式达到一致，后面会看到这一点。 一个射影变换是一个(n+1)*(n+1)矩阵A 使y*=Ax*, 与A相差一个数值因子的变换也是射影变换 光学中心、像平面、场景示意图 二、求解的思想 1. 建立系统的坐标变换描述, 坐标间的位置关系(主要有: 1) 世界坐标系到相机坐标系的平移和旋转变换,2) 相机坐标系到像坐标系的仿射变换, 确定需描述的系统参数. 2. 根据已知靶标上的景点坐标与像平面对应的像素坐标,建立方程组,求解方程组确定系统参数.完成系统的初步标定. 3. 根据标定的系统,系统参数已知, 计算相应景点的像素坐标用于检验偏差,或采用特殊的坐标系或特殊的位置关系检验系统标定的效果. 4. 上述是单相机监视系统, 根据单相机系统, 确定两个相机的摄象机矩阵P，P' 由此, 确定两部相机的相对位置. 注: 把标是平面图形, 因此, 两幅图片存在单应关系, 两幅图片间存在单应矩阵 三、单相机定位系统的描述 1. 世界坐标系中的坐标 到 相机坐标系的坐标的转换： Xw是景点X在世界坐标系中的表示，点X的坐标由世界坐标系转换至相机坐标系的变换为：Xc = R（Xw - t），t 为平移列向量，R为三阶正交矩阵，而Xc 是点X在相机坐标系下的坐标； 2. 景点以及对应像点的三维坐标： X的三维成像坐标，的推导：Xc的三维像点Uc 设 Xc=xc，yc，zcT ，Uc = -fxc/zc，-fyc/zc，-f T Uc的坐标推导，见下述示意图（应用成像原理针孔模型）相机坐标系的原点是相机的焦点，则由原点出发的射线上的任何点具有相同的像素坐标，这样的点在相机坐标系下的坐标具有性质：xc，yc，zc，xc', yc', zc'对应坐标成比例。 既存在a非零常数，：xc，yc，zc = axc', yc', zc'，这种关系是定义在相机坐标系下三维空间中点坐标的一个等价关系，利用等价关系将R3等价类的商空间，其中的类中每个景点成像具有相同的二维像素坐标。 3. 三维像点坐标到像素坐标的转换 1）OiXiYiZi (像平面欧几里得坐标系)： OXY面作为成像像素坐标平面，为此根据像素坐标特点，坐标系原点取在左上角，Z轴与相机坐标系Z轴有相同的指向，相机坐标系的Z轴与像素平面交点是像素平面的主点，像素坐标记为u0，v0； 2）OUVW坐标系(像素平面仿射坐标系)： W与Zi相同，V轴是Yi的伸缩，UV面是XiYi面的仿射 3）由Xc=xc，yc，zcT ，Uc = -fxc/zc，-fyc/zc，-f T 到U，V，W的变换上述的 变换矩阵记为 K 如果已知u,v坐标与xc,yc,zc 的数据, 有一对点的数据就可以得到两个方程, 对足够的信息量,就可以将参数确定下来 4. 景点坐标到像素坐标的变换 设世界坐标系转换至相机坐标系的变换为：Xc = R（Xw - t），t 为平移列向量，R为三阶正交矩阵，而Xc 是点X在相机坐标系下的坐标；则有5. 建立求解M的方程组上述可建立2k个12元的齐次方程组，可以通过Matlab中的 null 函数求解 6. K，R，t 矩阵的计算 M 是 3 行 4 列矩阵，前3列构成的方阵是A = KR，后1列为b = -KRt 1）利用Matlab中的函数 k，r = qr（A）可以实现从矩阵A得到，上三角矩阵K，和正交矩阵R； 2）t = A-1*b 这样就完成了系统的标定。 实现由景点的三维坐标计算像点的像素坐标，但其逆不是一一对应的，原因在于景点与像点是多对一的。 四、标定示例的计算 1.取世界坐标系OWXWYWZW原点为靶标中心，靶标平面为XOY面，Z轴指向相机方向。 于是可确定其上5个圆的圆心在世界坐标系中的坐标，在加上靶标中心的世界坐标系的原点，六个点的坐标列表。 xyz= 2. 靶标上的六点对应的像坐标的确定 1) 将靶标的像(题目本身给出的靶标的像就是1024*768分辨率）按1024*768分辨率,建立图像文件（复制到画板中，保存为24色真彩图或256色或16色索引图像，保存到Matlab中work文件夹。 24色真彩图的图像文件没有颜色表，图像矩阵为 1024*768*3 大小，256色或16色位图包含图像矩阵1024*768大小，和颜色表分表为 256*3，16*3。 2）转换为灰度图像 对24位真彩图 w = imread('babiaoxiang.bmp') %读出图像矩阵w wgray = rgb2gray(w) % 将w矩阵转换为灰度矩阵(注, w 为1024*768*3, 而wgray 是1024*768) imshow(wgray) % 显示灰度图像 对于索引图像(256色或16色色拉图） w, map = imread('babiaoxiang.bmp') % 读出图像矩阵和颜色表 wgray = ind2gray(w,map) % 转换为灰度矩阵 imshow(wgray) % 显示灰度图像 3) 转换为二值图像 选取适当阈值将灰度图像转换为二值图像 n=size(wgray) for i=1:n(1) for j=1:n(2) if (wgray(i,j)<50) wgray(i,j)=0; else wgray(i,j)=255; end end end 4) 提取边缘 （1）利用edge函数直接提取边缘（或定义其它提取边缘的算子） wgray=edge（wgray） （2）将边缘点集合定义出来 逐个像素判断是否是边缘，是，将像素坐标累计追加到边缘点坐标集合中。 5）将边缘点分类 （1）最小距离聚类 两个边缘点距离小于5个像素就被分成一类，按这个原则，就可以把所有边缘点分成若干类，每一类计算类中心像素坐标。 i）计算两两边缘点的距离矩阵 ii）初始时各点成一类，将距离最小的两点合并成一类，类数减少 iii）计算各类距离，再将距离最小的两类合并，直至最后合并为一类。 本题可以实施到类间距离大于5，不在合并。 （2） 已知分为五类， i）每一类中确定一个两个坐标平均值猜测的初始向量。 ii）然后，将其它边缘点按距离初始坐标平均向量距离最短的原则，进行第一次分类。 iii）将第一次分类结果得到的各类边缘点，求坐标平均值向量作为下一次分类的各类中心，继续分类。iv）迭代一定次数为止。 （3）有逐个边缘点按距离小于一定阈值进行分类6) 计算各类像素平均坐标，写出u(i) ，v(i) i=1,2,3,4,5,6 按均值聚类并计算各类边界点平均坐标的程序：clearw=imread('camrea.bmp');w=rgb2gray(w)n=size(w);for i=1:n(1) for j=1:n(2) if(w(i,j)<50) w(i,j)=0; else w(i,j)=255; end endendww=edge(w);n=size(w)ed=;imshow(w);for i=1:n(1) for j=1:n(2) if(w(i,j)=1) ed=ed,i;j; end endend%center=190,195,210,501,500;324,426,642,583,285;n=size(ed);m=104;for k=1:5 % 按K-均值聚类算法进行五次迭代 l=1; ed1=; for i=1:n(2) for j=1:5 mm=max(abs(ed(:,i)-center(:,j); if(mm<m) m=mm; l=j; end end m=104; ed1=ed1,ed(:,i);l; l=1; end % nn=size(ed1); center1=zeros(3,5); for i=1:nn(2) for j=1:5 if(ed1(3,i)=j) center1(:,j)=center1(:,j)+ed1(1:2),i);1; end end end for i=1:2 for j=1:5 center1(i,j)=center1(i,j)/center1(3,j); end end center=center1(1:2),:)endcentercenter = 190.9779 196.8821 213.4818 503.4054 501.2403 323.4081 425.2890 639.8097 581.7838 285.2747center = 190.9779 196.8821 213.4818 503.4054 501.2403 323.4081 425.2890 639.8097 581.7838 285.2747center = 190.9779 196.8821 213.4818 503.4054 501.2403 323.4081 425.2890 639.8097 581.7838 285.2747center = 190.9779 196.8821 213.4818 503.4054 501.2403 323.4081 425.2890 639.8097 581.7838 285.2747center = 190.9779 196.8821 213.4818 503.4054 501.2403 323.4081 425.2890 639.8097 581.7838 285.2747center = 190.9779 196.8821 213.4818 503.4054 501.2403 323.4081 425.2890 639.8097 581.7838 285.2747注：程序中所使用的图像为靶标的像保存为位真彩图像，文件名称为camreadingwei3. 单相机系统参数的标定 1) center =233.9779 239.8821 256.4818 544.2403 546.4054; 344.4081 446.2890 660.8097 306.2747 602.7838/3.78 单位换算为毫米center = 61.8989 91.1133 63.4609 118.0659 67.8523 174.8174 143.9789 81.0251 144.5517 159.4666 2) 靶标上圆心三维坐标: scenecoord=-50,50,0;-30,50,0;0,50,0;50,50,0;-50,-50,0;50,-50,0; 3) 建立单相机系统参数方程组 Mx=0;还应有一个条件，考虑到M = KR| -KRt，且K的第三行为 0 0 1，R的第三行为单位向量，故M矩阵的第三行前三列三个元素m31 m32 m33 平方和为 1 即m312+m322+m332=1系数矩阵的确定: center =233.9779 239.8821 256.4818 544.2403 546.4054; 344.4081 446.2890 660.8097 306.2747 602.7838/3.78'scenecoord=-50,50,0;-30,50,0;50,50,0;-50,-50,0;50,-50,0mm11=scenecoord,ones(5,1);mm12=zeros(5,4)mm13=(-center(:,1)*ones(1,4).*mm11mm21=zeros(5,4);mm22=mm11;mm23=(-center(:,2)*ones(1,4).*mm11;mm=mm11,mm12,mm13;mm21,mm22,mm23, 为参数方程的系数矩阵方程的求解: 在只知道五组点的对应坐标情况下,无法确定12个参数, 但由于是齐次方程, 因此只有11个参数独立, 即便这样,也无法确定所有参数,为此利用下述方法求近似解: 为求 mm*x=0;的解 改为寻求满足条件 m312+m322+m332=1 的条件下的 |mm*x|2的最小值 而|mm*x|2 = (mm*x)'(mm*x) = x'*mm'mm*x 问题转化为: 求解 min x'*mm'mm*x s.t. m312+m322+m332=1 约束条件也可以改写为 x'ax=1 其中a为12阶方阵, 其中,除a(9,9)=a(10,10)=a(11,11)=1外,其余元素皆为0 利用优化函数 fmincon 求解约束极值问题 编写目标函数文件 myfundingwei.mfunction y=myfundingwei(x)center =233.9779 239.8821 256.4818 544.2403 546.4054; 344.4081 446.2890 660.8097 306.2747 602.7838'/3.78scenecoord=-50,50,0;-30,50,0;50,50,0;-50,-50,0;50,-50,0mm11=scenecoord,ones(5,1);mm12=zeros(5,4)mm13=(-center(:,1)*ones(1,4).*mm11mm21=zeros(5,4);mm22=mm11;mm23=(-center(:,2)*ones(1,4).*mm11;mm=mm11,mm12,mm13;mm21,mm22,mm23;y=x'*mm'*mm*x; 编写约束函数文件 conobjdingwei.mfunction c,ceq=conobjdingwei(x)a=zeros(12,12);a(9,9)=1;a(10,10)=1;a(11,11)=1;c=;ceq=x'*a*x-1; 编写主程序文件 xiangjidingwei.mx0=(1/sqrt(3)*ones(12,1);options=optimset('largescale','off','tolfun',10(-10),'maxiter',20000,'tolx',10(-10),'maxfunevals',1000000,'tolcon',10(-10);x,fval=fmincon(myfundingwei,x0,conobjdingwei,options)mmm=x(1:4),x(5:8),x(9:12)'q,r=qr(mmm(1:3),(1:3)t=-inv(mmm(1:3),(1:3)*mmm(:,4) 4)参数的具体标定可由前期推导与上面的程序运行结果给出mmm = 0.0027 -0.0150 0.5774 1.8664 0.0166 0.0014 0.5774 2.2909 0.0000 -0.0000 1.0000 0.0176q = -0.1598 0.9872 -0.0008 -0.9872 -0.1598 -0.0011 -0.0012 0.0006 1.0000r = -0.0168 0.0010 -0.6634 0 -0.0150 0.4783 0 0 0.9989t = -145.3987 97.8651 -0.0136 5) 系统标定效果的验证将靶标圆心坐标输入，计算输出点的像素坐标：center =233.9779 239.8821 256.4818 544.2403 546.4054; 344.4081 446.2890 660.8097 306.2747 602.7838'/3.78scenecoord=-50,50,0;-30,50,0;50,50,0;-50,-50,0;50,-50,0mm11=scenecoord,ones(5,1);mm12=zeros(5,4)mm13=(-center(:,1)*ones(1,4).*mm11mm21=zeros(5,4);mm22=mm11;mm23=(-center(:,2)*ones(1,4).*mm11;mm=mm11,mm12,mm13;mm21,mm22,mm23;u=(mmm(1,:)*mm11')./(mmm(3,:)*mm11')v=(mmm(2,:)*mm11')./(mmm(3,:)*mm11')输出 u，v 坐标u = 61.1530 62.8877 69.0360 145.1199 143.5296v = 95.1520 112.9529 176.0459 81.3373 159.3397 6) 系统标定的矫正 将上述得到的标定后的系统输出的（u，v）像素坐标作为靶标上圆心的像点的像素坐标的估计，按同样方法重新标定系统，并计算靶标上圆心的重新映像坐标得到： x = -0.0017 0.0093 0.5774 -1.1562 -0.0103 -0.0009 0.5774 -1.4195 -0.0000 0.0000 1.0000 -0.0109fval = 8.6984e-007mmm = -0.0017 0.0093 0.5774 -1.1562 -0.0103 -0.0009 0.5774 -1.4195 -0.0000 0.0000 1.0000 -0.0109q = -0.1586 0.9873 -0.0008 -0.9873 -0.1586 -0.0011 -0.0012 0.0006 1.0000r = 0.0104 -0.0006 -0.6628 0 0.0093 0.4791 0 0 0.9989t = -145.9233 98.1776 0.0085 u = 61.1793 62.8987 69.0027 145.0867 143.5592v = 95.1861 112.9509 176.0193 81.3089 159.3670 而前一次标定输出的（u，v）为center = 61.1530 95.1520 62.8877 112.9529 69.0360 176.0459 145.1199 81.3373 143.5296 159.3397 4. 成像畸变的矫正 五、系统标定的效果检验 将标定好的系统，计算各景点的像素坐标，并与前期计算的像素坐标比较其偏差。给出效果分析。 六、非理想透镜成像畸变问题的考虑 上述的相机系统的描述是在透镜很细微时建立起来的，相机成像时不发生畸变，但对于一般情况，并不如此，通常发生几个像素大小的偏差，肉眼一般不容易观察到，但作为相机系统的测量与矫正，这种畸变给予一定的修正、补偿还是应该考虑的。 一般的透镜成像会发生两种情况的畸变，径项畸变（radical distortion）主要是线的弯曲，矩形的像边缘向内或外弯曲。 另一方面是主点（相机坐标系Z轴对应的景点成像点）成像的偏离，前面K矩阵中相机的五个内参数中，有一个焦距长度 f ，应用上一般取代 f 而引进相机常数（仍记为 f ）的参数，理想情况下，相机常数与焦距相同，而实际中，相机常数会略小于焦距。于是主点的成像也会轻微地发生偏离。 应用上通常设计适当的标定模式图像估计参数的变化。 两种畸变多数情况下采用多项式修正模型进行修正，若在前述相机系统的标定过程中，以上述公式进行修正建立求解标定矩阵，会增加一个待定参数。 七、 利用现有把标以及在两部相机的图像确定两部相机相对位置的数学模型。 设两摄象机坐标变换分别为: Xa = K1(R1XW-t1) Xa = K1(R1XW-t1) 三、 两视点几何两视点几何双相机系统 简介： 立体（stereo）视觉可以根据不同位置的视觉信息给出空间场景的深度信息，标定（Calibration）一部相机的系统，需要弄清一条光线的走向，而标定定两部相机的系统，确定三维场景点的三维坐标，需要计算两条光线的交点，这就是立体视觉的基本原理，他主要有三个问题： 1）相机标定（确定相机的内外参数）； 2）创建两部相机拍摄的图像上点的对应关系(图像配准）； 3）重新构建景点的三维坐标。一、极几何：1. 极平面: 通过两个摄象机光心的平面, 两个摄象机光心的连线称为 基线, 极平面的全体构成共基线的平面束;2. 极线: 极平面与摄象机象平面的交线 称为极线, 同一张极平面与两张像平面交线称为 一对极线对应.3. 极点: 基线与像平面的交点称为极点, 两个极点分别为摄象机光心在另一个摄象机像平面上的投影.二、极几何约束 设 m，m是两幅图像上的一个点对应，极几何约束给出了一个必要条件，这在作两幅图像的配准时，很有用。 m，m是两幅图像上的一个点对应，则 m 位于m' 对应的极线上, m' 位于m 所对应的极线上, 即 m lm' m'lm' 极几何约束不涉及几何结构, 只反映射影空间的属性。三、极几何的代数表示- 基本矩阵 假定两个摄象机矩阵分别为 P，P' ，两个摄象机像平面分别为：I，I ' 则对任意的m I , 反投影线 的参数方程：X（s）= P+m+sC 注：设 X 是空间中一点，其在 世界坐标系、相机坐标系中的坐标分别为 XC = (xc,yc,zc,1)，X = (x,y,z,1) 则 XC=PX 设C为摄象机光心在世界坐标系中的齐次坐标, C 为非齐次坐标, C=C ', 1 ' P = KR(I-C ), R(I-C ) 称为摄象机的外参数矩阵, K 称为内参数矩阵. P+是P的广义逆矩阵, 即PP+=I , C是第一个摄象机的光心, 即PC=0. m点出发的反投影线 的 方程 X(s)= P+m+sC , s为参数. 从而 P(X(s) = P(P+m+sC) = PP+m+sPC = m , X(s) 的像是 m 基本矩阵 F 描述了像点与其对应的极线之间的关系： lm= Fm 又由于m 在第二个像平面上的对应点m' 在 lm上 所以必有 m'Fm=0 , 且 Fe=0, FTe'=0 四、基本矩阵的几个例子 例1 假定第一个摄象机内参数矩阵 K 第二个摄象机内参数矩阵K' ，第二个摄象机相对于第一个摄象机的方位为 （R，t） 即：两部摄象机之间的坐标转换矩阵为 X' =RX+t , R为旋转, t为平移向量. 于是两摄象机矩阵可表示为：P=K(I，0)，P'= K'(R,t) 例2, 纯平移运动下的基本矩阵 当摄象机沿x轴做纯平行移动时, 得到的两幅图像, e' =e = (1, 0, 0)T
