对于数码相机定位问题的分析.doc

高教社杯全国大学生数学建模竞赛承 诺 书我们仔细阅读了中国大学生数学建模竞赛的竞赛规则.我们完全明白，在竞赛开始后参赛队员不能以任何方式（包括电话、电子邮件、网上咨询等）与队外的任何人（包括指导教师）研究、讨论与赛题有关的问题。我们知道，抄袭别人的成果是违反竞赛规则的, 如果引用别人的成果或其他公开的资料（包括网上查到的资料），必须按照规定的参考文献的表述方式在正文引用处和参考文献中明确列出。我们郑重承诺，严格遵守竞赛规则，以保证竞赛的公正、公平性。如有违反竞赛规则的行为，我们将受到严肃处理。我们参赛选择的题号是（从A/B/C/D中选择一项填写）： 我们的参赛报名号为（如果赛区设置报名号的话）： 所属学校（请填写完整的全名）： 参赛队员 (打印并签名) ：1. 2. 3. 指导教师或指导教师组负责人 (打印并签名)： 日期： 年 月 日赛区评阅编号（由赛区组委会评阅前进行编号）：高教社杯全国大学生数学建模竞赛编 号 专 用 页赛区评阅编号（由赛区组委会评阅前进行编号）：赛区评阅记录（可供赛区评阅时使用）：评阅人评分备注全国统一编号（由赛区组委会送交全国前编号）：全国评阅编号（由全国组委会评阅前进行编号）：对于数码相机定位问题的分析摘要：本文针对双目定位的数码相机相对位置的标定问题进行了深入研究。对于问题一，建立了世界坐标系和相机成像坐标系之间的转换矩阵，可以根据相机与靶面任意的相对位置关系求得靶标圆心的像坐标；对于问题二，针对所给靶标的特征提出了质心模型、切线模型、解析法模型三种求靶标圆心在像坐标投影点的模型，并就所给的靶标相片求出了圆心投影的像坐标。对于问题三，应用了问题一中转换矩阵，建立了判定问题二中模型优劣的仿真算法，对问题二中提出的三种模型进行了优劣分析。对于问题四，应用Roger Y. Tsai的单部相机内部和外部参数的标定算法，不用求解多元非线性方程组，直接用最小二乘法求解超定线性方程组，即可求得相机坐标系的变换矩阵和像距，从而能够确定两相机的相对方位和位置关系。本文得到的主要结论如下：问题一：求解靶标圆心像坐标的算法为世界坐标系和相机成像坐标系之间的坐标变换；问题二：利用三种模型得出问题所给靶标相片中靶标圆心投影位置的像坐标，如下表所示；质心模型切线模型解析法模型XYXYXYA322.8948189.4935323.2771189.9167323.0474189.6938B422.9960196.9423423.2943197.3401423.0209197.0611C639.8994213.1522640.1433213.3946640.1127213.2289D582.7329502.9824583.0062503.2098582.8821503.1650E284.6657501.7731284.9940502.0426284.8034502.0957问题三：在像平面和物平面夹角不是很大的情况下，三种模型的精度相差不大；像平面和物平面夹角比较大时，切线模型的精度大于解析法模型，质心模型的精度最差。问题四：运用相机内外参数的标定算法和问题二得到的靶标圆心和切点坐标，得题图三对应相机的相机坐标系的旋转矩阵与平移向量分别为，。关键词：机器视觉、坐标变换、透视投影、相机标定对于数码相机定位问题的分析一、 问题重述双目视觉定位系统是机器视觉学科的主要研究内容，系统通过处理放在两个不同位置的摄像机捕捉到的图像，对空间中的某点或某物体进行定位。这种技术在交通、医疗、工业、军事等领域都有广泛的应用。要使定位系统进行准确的定位，必须知道两台相机精确的相对位置关系，可以使用如下方法对两台相机进行标定以获得相机的相对位置：在一块平板上画若干个点， 同时用这两部相机照相，分别得到这些点在它们像平面上的像点，利用这两组像点的几何关系就可以得到这两部相机的相对位置。然而，无论在物平面或像平面上我们都无法直接得到没有几何尺寸的“点”。实际的做法是在物平面上画若干个圆（称为靶标），它们的圆心就是几何的点了。而它们的像一般会变形，所以必须从靶标上的这些圆的像中把圆心的像精确地找到，标定就可实现。本课题要求完成下列任务：（1） 建立数学模型和算法以确定靶标上圆的圆心在该相机像平面的像坐标, 这里坐标系原点取在该相机的焦点，x-y平面平行于像平面；（2） 由靶标及其像，计算靶标上圆的圆心在像平面上的像坐标（3） 设计一种方法检验你们的模型，并对方法的精度和稳定性进行讨论；（4） 建立用此靶标给出两部固定相机相对位置的数学模型和方法。二、 问题分析该问题实际上是求解不同坐标系中的点的对应关系的问题。对于第一个问题，首先需要解决的是空间圆成像的映射问题。由于实际情况中维度的差距，可以考虑将成像的模型设置为小孔成像的模型，从而省略了成像的畸变对定位造成的影响对问题进行简化，对于由此造成的误差在后面的模型分析中进行讨论。空间圆小孔成像中有一些基本的原则和理论需要进行分析，然后以此进行建模，考虑空间中各个参考系之间的坐标变换关系，通过旋转、平移等方式达到点的映射的目的。对于问题二，即考虑上述过程的反向过程：已知像的信息和原象的信息，要求得出原象的圆心在像上的位置。一个朴素的想法就是考虑质心的对应关系，如果能找到像和原象的质心的对应关系，就求出了圆心在像平面上的坐标。但是此种做法缺乏理论依据，只是根据生活经验得出，但不失为一种参考的方法，可以将其与其他方法进行对比，选出更优的方法。再者，还可以通过切线的垂线的交点来确定圆心。首先对给出的像的信息进行处理，拟合出像中的曲线方程，接着在物的由物平面中的切点求出它们在相平面中的坐标（用问题一的方法），之后就可以由切点求出切线的方程从而得出圆心的坐标。另一个想法是解析的方法。通过空间解析几何和平面上线、角的关系计算出圆心相对于空间坐标系的坐标，并由第一问中各坐标系的坐标变换关系就可以得出圆心在像平面上的坐标。以上方法都在理论推导的过程中进行了一些近似从而会导致系统误差，所以在由上述一些方法对问题进行分析过后还需要对它们的精确程度和稳定性进行分析。此即问题三的要求。而对于不同的方法，应采用不同的误差分析方法以得出他们的精度和稳定性。对于质心的方法，由于其缺乏理论根据，所以只能将其结果与其他的方法的结果进行比较以考虑其精度；对于切线的方法，由于计算机拟合的过程很精确（点足够多的情况下），且在之后的处理过程中没有进行近似处理，其精度的讨论应该也着重与其他方法的比较；对于解析的方法，由于其通过公式的推导，经过适当的近似得到结果，因此应着重讨论其近似时产生的误差。而对于稳定性的讨论，应朝着考虑靶标的移动、标定形状、噪声对方法精度的影响。最后一问要求通过靶标的世界坐标系坐标和靶标在两部相机的成像坐标系中的像点坐标来确定两部相机的相对位置。具体来说，即已知靶标上若干物点的世界坐标和成像坐标，寻找标定两个相机坐标系相对世界坐标系的旋转矩阵 和平移向量 以及相机的内部参数的算法。由这四个矩阵能够自然地求得两相机的相对位置。三、 模型假设1本题中不考虑镜头畸变；2本题中相机模型为小孔模型；3各空间坐标轴均为右手系。四、 符号与术语说明透射投影：坐标变换中的旋转矩阵。：坐标变换中的平移向量。：像距，光心到像平面的距离。世界坐标系：，被拍摄物体所在坐标系。相机坐标系：，以相机光轴为，光轴平行于相机像平面。成像坐标系：，像平面坐标系。五、 模型的建立与求解问题一该问要求靶标上圆的圆心在相机像平面上的坐标，即世界坐标系（）、相机坐标系（）和成像坐标系（）三个坐标系之间的点坐标的转换。1 世界坐标系和相机坐标系之间的转换该转换用旋转矩阵和平移向量来描述，世界坐标和相机坐标间满足下式：其中（具体确定见下文），即世界坐标先旋转矩阵，再将坐标原点平移向量成为相机坐标。为方便表达，上式亦可表示为，其中，而旋转矩阵的确定有多种方法，如绕坐标轴的旋转变换，欧拉角旋转变换等，考虑到直观易用性，本文选取绕坐标轴的旋转变换来确定矩阵，分别为坐标系旋转变换时绕固定的世界坐标系坐标轴的旋转夹角，此时为，图1 坐标系的旋转又应用z-x-y旋转顺序，即先绕世界坐标系的z轴旋转，再绕其x轴旋转，最后绕其y轴旋转，但每次旋转都是以固定的世界坐标系坐标轴为转轴的。图2 相机坐标系和成像坐标系2 相机坐标系和成像坐标系之间的转换对于该问题，用小孔模型代表相机工作原理，相机坐标系以小孔（光心）为坐标原点，以光轴为轴正方向，而成像平面坐标系则平行于相机坐标系的平面，如图2所示。原理上说，成像平面应当在物点P 和光心O连线的延长线上，完成小孔成像，但由于成倒像，为处理和讨论方便，将成像平面放置在光心和物点之间，距光心的位置上，此处为像距。由小孔成像原理和相似三角形知识，得到相机坐标系和成像坐标系之间的关系为：，其中，为P点的成像坐标，为P点在相机坐标系中的空间坐标。矩阵形式表示为：。综上所述，可得世界坐标系和成像坐标系之间的变换关系为： (1)在相机内部参数都确定的情况下，M为一常矩阵，为随着世界坐标系中点坐标变化的量，为列向量的第三项元素。因而得到了世界坐标系和成像坐标系之间的变换关系。对于本小问，输入为世界坐标系中靶标的5个圆心位置以及数码相机的6个位置参数（旋转矩阵中的3个旋转角度和平移向量的3个变量），要求输出靶标圆心在像平面上的坐标。对于靶标，定义如图3的世界坐标系，5个圆心放置在xoy平面上，其坐标如表1所示；相机的位置参数则由外部设定6个参数。用Matlab编写程序（代码见附录），可得到应用上述坐标变换式所得圆心对应的像坐标的和。由于题中要求坐标系的原点定在光心，因而最终像坐标为。图3 世界坐标系表1世界坐标系中靶标圆心的坐标圆编号圆心x坐标圆心y坐标圆心z坐标A-50500B-20500C50500D50-500E-50-500问题二问题二要求根据靶标图像和靶标的像，靶标上圆的圆心在像平面上的像坐标。1 图像的预处理由于题中靶标图像的照片给出的是灰度图像，并且实际摄像头拍摄的照片也是灰度图像，所以提取圆心工作以前，必须对图像进行预处理工作。图像预处理流程如图4所示，读入图像利用软件提取图形文件中像素的颜色信息，由于圆心投影只与边缘有关，所以用软件算法实现边缘的提取，本文利用了Sobel边界探测器，由于原图像灰度及阈值设定的原因，完成边界探测后会有许多与边界不连续的点，我们称之为噪声，这里我们手工去除这些噪声点，然后将边界点坐标存入数组。同时我们将图像做二值化处理，作为求质心的原始数据。图5显示了图像预处理的效果，A为初始读入图像数据，是灰度图像，B为经过Sobel边界探测后的结果，有噪点存在，C为手工去除噪声点后的图像，D为原始图像做二值化处理后的图像。图4 图像预处理流程 A B C D图5 图像预处理阶段示意2 圆空间透视投影的分析本小节介绍推导模型必须证实的几个定理定理1：空间内圆的投射投影为椭圆证明：如图6所示以为圆心的圆经过小孔投影到像平面上，以透镜的光心为坐标原点建立空间三位坐标系，Z轴垂直于像平面向下，O到像平面的距离为，O到圆所在的平面的距离为。图6 空间圆的透视投影设圆以为轴旋转角度，再像平面形成一个封闭曲线，下面证明该封闭曲线为椭圆。圆的参数方程为：根据相似三角形比例原则，圆上一点与象平面中该点的坐标有如下关系整理得：此方程为一椭圆的标准方程，由此可知圆通过照相机在像平面上成的像为椭圆（在时为圆）。定理2：透射投影中，物圆的圆心的像不是所成椭圆像的中心。证明：图7 XOZ平面AB剖面图过图6的的面作圆锥的剖面如图7所示： 其中，为圆在像平面上所成的像，过作的平行线分别交的反向延长线为，交于。假设圆心的像为椭圆的中心，则有由相似三角形边的关系有：将上式代入下式有。分别从、向作垂线垂足分别为、，由于为的中点，所以有，由图易知：。与上述矛盾，故假设不成立，即：物圆的圆心的像不是所成椭圆像的中心定理3：空间圆圆心的投影在投影椭圆的短轴上证明：首先说明圆的圆心的像、椭圆中心都在。由于是直径，易知圆的圆心的像上。图8平行于转轴的剖面图如图8所示，由于圆以为轴旋转角度，而轴过圆心，所以椭圆的中心在圆上的原象一定在垂直于轴的直径上，在图8中的表述即为椭圆的中心在的中点，其原象在上，所以为过椭圆中心的一条轴。下面证明该轴是椭圆的最短轴：设、的坐标分别为、由相似三角形边的关系及1中坐标的关系有：由上两式可得，在上为单调递减，又由于过椭圆中心，所以转角最大的直径所对应的像即为椭圆的最短轴。即物圆的圆心的像在其像椭圆的最短轴上。3 模型一像素质心模型31模型建立一般情况下，平面上一薄片的质心计算公式如下，为薄片的范围，为薄片在(x,y)的面密度，为质心坐标。本模型因为是基于像素的模型，像素是一个个离散的点，又由于图形做了二值化处理，这里为二值量。定义：则在像素意义下的中心坐标定义式如下所示：其中(x1,y1)，(x2,y2)为包含此椭圆而不包含其他椭圆的一个矩形的两个对角顶点的坐标。经过2的分析，虽然空间圆圆心的投射投影一般不在投影图形的几何中心，而且利用像素所作的离散计算存在着一定的误差，但是作为一种近似的方法，也是可行的，后文将对这种模型造成的误差做定量的分析。32模型求解按上述模型求得的圆心坐标如表2所示，表2 切线模型求点坐标圆XYA322.8948189.4935B422.9960196.9423C639.8994213.1522D582.7329502.9824E284.6657501.77314 模型二切线模型在上述的质心模型中，仅利用了一个椭圆自己的性质就给出了该椭圆内圆心的投影，而本问题给出了五个相等大小并且相对位置已知的圆，所以如果利用圆之间相互关系求解各圆圆心的投影坐标将更为精确。41 切点位置不变的性质空间的投射投影有如下的性质，如图9所示，如果空间存在圆A，B，这两个圆有公切线，切点a,b,c,d。如果将这两个圆做某透视投影如图所示，此时在做投影椭圆A，B的公切线，则得到的四个切点a,b,c,d分别是原来a,b,c,d的投影。图9切点位置不变的性质示意图投射投影变换属于仿射变换，点在线上的位置关系不会改变，在本问题的靶标图形中，有四个相等大小的圆(A,C,D,E)正好在正方形的四个定点上，两两做这四个圆的外公切钱，如图，可得相对的切点的连线经过圆心，两条连线的交点就是圆的圆心。所以在投射的象的图形中，只需找到这些切线与切点，连接各个切点形成交线就可以求得圆心。对于问题中的圆B，由于不存在相互垂直的切线，所以只能够确定一条直径的位置，另外根据2的叙述，圆心的投影肯定在投影椭圆的短轴上，这样用已知直径的投影和椭圆的短轴相交，就可以求得B圆圆心投影的坐标。42 椭圆方程的拟合本模型要求作椭圆的曲线，由于所给图像分辨率不高，仅为1024*768像素，所以如果直接用图像中图形边界做切线，精度将会变得非常低，会造成很大的误差，所以在本模型中，先要利用所给图像中图形的边界(在1中提取)拟合出椭圆的方程。平面中，二次曲线的一般方程是：其中将空间中任意一点带入F，被称为空间点到圆锥曲线的“代数距离”。拟合一条圆锥曲线可以看作是求所有被拟合点代数距离平方和的最小值。如下式所示：设用N个点拟合一个椭圆，又根据文献2中描述如果限制可行解的向量肯定是椭圆的参数向量。由此可以把椭圆拟合的问题转换为一个带二次限制条件的最小二乘问题，如下：其中图中五个椭圆方程拟合的结果如表3所示。表3 椭圆方程参数圆abcdefA0.49083680.045627320.5103946-325.6590-208.190271450.98B0.50054520.071117930.5019815-437.4560-227.8848114165.5C0.52528670.13234670.4842668-700.4724-291.0644254440.3D0.47898800.18026210.5388937-648.9240-647.1483351269.2E0.45211530.10417250.5589569-309.6662-590.5596191540.0图10为实际过程中椭圆拟合的效果图，从图中可以看出，方程对于椭圆做了很好的拟合，图10 椭圆拟合效果图表4显示了拟合过程中代数距离平方和的最小值，可看出大约在数量级，由于每个椭圆用大约200个点进行拟合，所以应该在20左右，x,y坐标平均的误差应该在几个像素以内。表4 代数距离平方和圆ABCDE平方和52242.5745703.2339525.2030956.1140892.5243 公切线的求取图11靶标的公切线图设两椭圆方程分别为A椭圆在处的切线斜率为：B椭圆在处的切线斜率为：如果，连线是切线，那么联立上述五式可求得，和的解，其中，关于，是四次方程，既可以求得四组解，这和切线的物理意义（两组内切，两组外切）是一致的。实际求解过程当中，由于计算规模过大，不宜求解，所以这里我们采用了将最后一个条件改为的优化条件，然后限定，的范围，用LINGO软件解优化(程序见附录)问题，经试验，效果良好，不仅求解速度很快，而且每次优化目标都等于0，说明求得了正确解。44 模型求解按上述模型求得的圆心坐标如表5所示，表5 切线模型求点坐标圆XyA323.2771189.9167B423.2943197.3401C640.1433213.3946D583.0062503.2098E284.9940502.04265 模型三 解析法模型51 模型建立在图12所示在照相机坐标系中圆的半径是，设该圆绕 y轴旋转角，图中坐标系原点设在焦距中心，圆心，距坐标系原点的距离为，摄像机焦距为，圆的参数方程为根据相似三角形比例原则，圆上一点与图像 中该点的坐标有如下关系整理化成椭圆的标准方程（在像平面上）为：图12空间圆的透视投影图过图12中的面作圆锥的剖面如图7所示：，（设圆的半径为），易知AB在像平面上的像为所形成椭圆的短轴。而过圆心且与AB垂直的半径CD在像平面上的像为一条平行于圆在像平面上的像（椭圆）的长轴的直线，设其长度为，设其像椭圆的长轴长为，短轴长为。如图8所示，设圆心到原点的距离为，圆心到像平面的距离为。则由相似三角形的关系有： （2）又由小孔成像的原理，对应的直线在成像后大小变化的比例是相等的，反映在上左图中即： （3）由上面的椭圆标准方程可知圆所成的像的中心并不是圆的圆心在像平面上的对应点，圆心的坐标为与上面（*）（*）式联立可以得出如下关系： （4）为了确定圆心的位置，由于的纸无法直接得出，所以取椭圆的半长轴代替，可以得出圆心的近似位置。52 参数求解对于式（*），像距f已知，需要解得椭圆的半长轴和半短轴，对于任意椭圆的一般式，总可以将其化成如下形式：根据系数相等可得：解得由此可得椭圆的圆心为，椭圆的半长轴为，半短轴为椭圆旋转角度为53 模型求解通过5.2求得拟合的各椭圆的长轴和短轴，转角，中心，然后通过5.1中的(2)-(4)式求得空间圆圆心在象平面投影的坐标。表6 空间圆圆心在象平面投影的坐标圆XYA323.0474189.6938B423.0209197.0611C640.1127213.2289D582.8821502.0957E284.8034502.0426问题三1模型的检验方法该方法的核心是基于问题一的结论，在问题一中，已经获得靶标上的点和像平面上像坐标的对应关系，因此可以按以下步骤完成基于仿真的检验：步骤1 固定靶标，确定相机相对于靶标的姿态，包括旋转矩阵和位移矢量。步骤2通过问题一结论中的变换公式，求出靶标上所有点映射到相机像平面的坐标，这里必须使用离散化算法。（离散化算法解释）这一步实际上就仿真求出了在此相机姿态下，靶标在相机像平面所成的像，同时，可以准确的知道靶标上圆的圆心在相机像平面投影的坐标。步骤3 用问题二中的多种方法对步骤2中求出的靶标的像求圆心步骤4 用所求得圆心和映射在相机像平面的准确的圆心的像做比较，两者相差越小说明方法越准确。2检验结论(a) (b)图13 成像检验图表7 图13(a)检验结果标准圆心质心模型切线模型解析法模型XyXYxyxy519.9263309.5831519.9351310.3686519.9333310.4732520.1446309.9762589.1007311.0271589.0412311.8446589.1136312.0014588.9056311.4932747.6679314.3308747.5765315.1844747.5763315.4256747.7089315.1042722.0924518.4924505.9950518.2396723.0138518.3788721.8761518.3584506.1380517.4491722.1016519.3891504.1380517.7691503.8780517.4360误差平方和3.4754167.9130146.229838表8 图13(b)检验结果 标准圆心质心模型切线模型解析法模型Xyxyxyxy484.6105361.5820486.8028361.5050485.7325361.5549485.6051361.5320616.7665348.9821616.7102349.1967616.7449349.0115616.7233349.1549914.1542320.6321914.3273320.8787914.2246320.6122914.2781320.6257813.8254506.6421813.8603507.2887813.8434506.6281813.8398506.6377466.4812544.9864466.3966545.6107466.4572545.1974466.4355545.4902误差平方和5.7683241.3119191.294977表7-8显示了用上述方法检验问题二中方法的结果，生成的虚拟照片如图13所示。其中图13(a)为像平面与物平面近似平行的情况，从结果中可以看出，在这种情况下，利用质心模型求解投影模型的误差是十分小的，在某些情况下，质心模型求解的误差小于用理论上更精确的切线模型和解析法模型。经过分析，我们认为，这是因为由于有限分辨率造成的量化误差其实可以等效为均值在质心的二维高斯噪声，而求质心相当于对噪声作了平均，一次可以抵消一部分噪声的影响，而切线模型和解析法模型是基于二次曲线最小二乘拟合的，这不能消除高斯噪声的影响。但图13(b)中，像平面与物平面有较大的角度，这时，质心模型有较大的误差，而切线模型和解析法模型的误差较小。问题四此小问即要求通过靶标的世界坐标系坐标和靶标在两部相机的成像坐标系中的像点坐标来确定两部相机的相对位置。具体来说，即已知靶标上若干物点的世界坐标和成像坐标，标定两个相机坐标系相对世界坐标系的旋转矩阵和平移向量以及相机的内部参数。由这四个矩阵能够自然地求得两相机的相对位置。由问题一可以看到，本小问是问题一的反向求解，即已知两种坐标系下点的坐标，求转移矩阵的过程。但直接通过式(1)来求解是非常困难的。原因有二：第一，是3*4的矩阵，有12个未知数， 需要解12个方程，即使借助计算机软件也很难完成任务，Matlab的Solve函数求解能力不足以解决12个方程，Lingo也应变量过多无法在有限时间内给出结果；第二，通过消去，每个靶标上的物点及其对应成像点可提供2个方程，因而共需6个点才能解出未知数，但本题每个靶标上只有5个圆，直接提供5个圆心，提供信息不足。因此，需要更直接的标定相机位置的方法，在参考了文献14后，我们发现了不用解多元或非线形方程的标定方法。该算法依然基于问题一的式(1)，同时充分利用了靶标的自由设定来减少未知量。现描述标定一部相机的算法。、1标定单部相机内部和外部参数的算法步骤1：令靶标上个物点的世界坐标为，物点在成像坐标上的坐标为，矩阵和的定义同问题一，则当时即可用最小二乘法由上式解得；步骤2：令当矩阵没有任何一行或一列为全0时（有一行或一列为全0的情况见文献1），其中，。由此可得的大小；步骤3：假设的符号为正，在个物点中选择一个成像坐标远离坐标原点的物点，计算如果与同号并且与同号，则的符号为正，否则为负；步骤4：由于旋转矩阵为单位正交阵，因此由四个量即可确定出， 而可由与叉乘得到，的符号的具体确定请请参阅文献1；步骤5：构造矩阵，其行为，再构造向量，然后解线性方程，即得相机像距和转移矩阵中的。至此成功标定了相机的内部参数和全部外部参数。算法MATLAB代码见附录。2标定两台照相机相对位置的方法使用此算法可以得到两部相机各自的和，即得到如下方程组：消去后得到：即得两相机坐标系之间的相对位置关系：3针对本题的算法实例使用问题二得到的成像坐标系中五个圆的圆心坐标，再结合问题一中的世界坐标系的假设，使用上述算法，编制Matlab程序得到了拍摄题目中图三的相机的旋转矩阵和平移向量： 再用问题一的算法，使用世界坐标系中的靶标圆心坐标和矩阵、，运用Matlab反推出其在成像坐标系中的像（重塑值），及其与问题二中得到的圆心像坐标（理论值）对比如图14。可见两者基本重合，1中所述算法可认为有效。图14靶标圆心点的理论值和重塑值的像坐标对比尽管五个圆心可以标定出相机的位置，但上述算法的特点是参与标定的物点越多，得到的变换矩阵越精确。因而，我们将第二问中得到的四个圆的十六个切点的两套坐标也代入算法中，此时得到的和为： 可以发现此时的变换矩阵和只用5个圆心得到的矩阵近似相等，但显然和更为精确，因此将其作为题中拍摄出图3的相机的标定位置。图15是圆心和切点的理论值和重塑值对比图。由图中两者的近似重合亦可验证算法的有效性。图15靶标圆心点和切点的理论值和重塑值的像坐标对比六、 模型评价和模型改进1关于问题一，完全按照坐标变换和几何关系进行推导，若不考虑相机镜头的畸变，则模型是精确的，因而改进方法为尽量研究抗量化噪声的算法。但对于问题一的坐标变换模型，世界坐标系到相机坐标系的变换为先旋转后平移，不如先平移后旋转直观。因为旋转后的平移相对于新坐标轴，这个平移量比基于原始坐标轴，即未旋转的坐标轴，要来得复杂和抽象。因此，先平移后旋转的坐标变换方法能使得用户或研究人员更容易地确定相机坐标系方位。2关于问题二提出的三种计算圆心投影的方法中，都遇到了相片分辨率有限存在量化误差的问题，这会造成圆心投影计算的误差。其中质心法和解析法均已做了近似处理，而切线法的精度仅取决于二次曲线拟合的精确与否。一般情况下，精度切线模型 解析法模型 质心模型，运算复杂性也是切线模型 解析法模型 质心模型，所以在精度要求不高，并且物像平面基本平行的情况下，质心模型也不失为一种好的选择。3对于问题四的相机参数标定算法，参量符号的确定方法不够严谨而比较随意，特别是转移矩阵中的和的符号有两种组合方式，要根据后面算出的的符号来确定，不具有优秀算法所必须的确定性。参考文献：1Roger Y. Tsai，A Versatile Camera Calibration Technique forHigh-Accuracy 3D MachineVision MetrologyUsing Off-the-shelf TV Cameras and Lenses，IEEE JOURNAL OF ROBOTICS AND AUTOMATION, VOL. RA-3, NO. 4, AUGUST 19872 Andrew Fitzgibbon, Maurizio Pilu, and Robert B. Fisher. Direct Least Square Fitting of Ellipses, PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. 21, NO. 5, MAY 19993刘红宁，基于双目立体视觉拟人机器人定位关键技术的研究，河北工业大学硕士学位论文，p7-28，2006.124方超，双目CCO成像目标器识别算法研究，南京理工大学硕士学位论文，p4-28，2007.65Zhang, Z，A Flexible New Technique for Camera Calibration，Pattern Analysis and Machine Intelligence，Volume 22, Issue 11，p2-8，Nov. 20006韩庆瑶 万书亭，机器视觉中空间圆透视投影特性研究，华北电力大学学报，第25卷第1期，p79-83，1998.1附录：1.问题一Matlab程序代码clcclear%ques1 using coordinates trans.f=1577/3.78; %mmtrans1=f 0 0 0;0 f 0 0;0 0 1 0;a=0; %rotate angle inputb=pi/2;c=-pi/2;%rotation matrixR= cos(b)*cos(c)+sin(a)*sin(b)*sin(c), cos(a)*sin(c), -cos(c)*sin(b)+sin(a)*sin(c)*cos(b); -cos(b)*sin(c)+sin(a)*sin(b)*cos(c), cos(a)*cos(c), sin(c)*sin(b)+sin(a)*cos(c)*cos(b); cos(a)*sin(b), -sin(a), cos(a)*cos(b); t=0 0 -400; %displacementtrans2=R t;0 0 0 1;targ1=0 -50 50 1;targ2=0 -20 50 1;targ3=0 50 50 1;targ4=0 -50 -50 1;targ5=0 50 -50 1;targ=targ1 targ2 targ3 targ4 targ5;img=zeros(5,2); %image pointss coordi.for i=1:5 mytemp=trans2*targ(:,i); zc=mytemp(3); w=trans1*trans2*targ(:,i)./zc; img(i,1)=w(1); img(i,2)=w(2);endplot(img(:,1),img(:,2),.);grid on;axis(-520 520 -520 520);% px=-50 -20 50 -50 50;% py=50 50 50 -50 -50;% plot(px,py,.);2LINGO拟合椭圆方程的代码MODEL: SETS:POINT/R1.R184/:X,Y; ENDSETS DATA:X=OLE(a5,X);Y=OLE(a5,Y); ENDDATA MIN=SUM(POINT:(A*X*X+B*X*Y+C*Y*Y+D*X+E*Y+F)2); B2-4*A*C=-1; FREE(A); FREE(B); FREE(C); FREE(D); FREE(E); FREE(F);END3问题四算法的MATLAB代码clcclear%ques1 using coordinates trans.f=1577/3.78; %mmtrans1=f 0 0 0;0 f 0 0;0 0 1 0;targ1=-50 50 0 1;targ2=-20 50 0 1;targ3=50 50 0 1;targ4=50 -50 0 1;targ5=-50 -50 0 1;targ=targ1 targ2 targ3 targ4 targ5;img=zeros(10,5,2); %image pointss coordi.%for k=1:10 %rotate angle input a=rand*pi; b=rand*pi; c=rand*pi; TT(:,k)=rand(3,1)*300; %displacement %rotation matrix% R= cos(b)*cos(c)+sin(a)*sin(b)*sin(c), cos(a)*sin(c), -cos(c)*sin(b)+sin(a)*sin(c)*cos(b);% -cos(b)*sin(c)+sin(a)*sin(b)*cos(c), cos(a)*cos(c), sin(c)*sin(b)+sin(a)*cos(c)*cos(b);% cos(a)*sin(b), -sin(a), cos(a)*cos(b); R = cos(b)*cos(c) cos(b)*sin(c) -sin(b); sin(a)*sin(b)*cos(c)-cos(a)*sin(c) sin(a)*sin(b)*s