数学建模优秀论文-数码相机定位的数学模型

数码相机定位的数学模型 摘 要 随着数码相机定位在各领域的广泛应用，对相关问题机器视觉的研究也成为热点。因此建立一个精度较高，稳定性好的数码相机定位的数学模型，具有很好的现实意义。 问题 1 要 求 给出确定靶标上圆的圆心在给定相机像平面的像坐标的算法 ，问题 2利用问题 1 的模型对给定数据求解。 为此 ， 首先建立了四个空间 直角 坐标系，在把图 3 的数字信息提取出来，主要是五个椭圆的边缘点的信息；同时为了便于运算，通过坐标变换 将 计算机图像坐标变换为图像坐标；并用提取的图像边界坐标拟合出 5 个椭圆的方程，利用 “ 曲线 切线的投影仍与曲线的投影相切，而且切点的投影仍为投影的切点 ” 这一引理，提取出靶标上圆及其像 上 的公切点的坐标作为特征点，利用 步法标定过程和最小二乘法建立了计算世界坐标系到相机坐标系的旋转变换矩阵 R 和平移向量 T 及径向畸变系数 k 的算法 。 利用 16 个公切点作为特征点，通过 程求得靶标上圆的圆心在文中给定相机像平面的五个坐标 (单位： A( B( C( D( E( 问题 3 的解决 分为两步： 一是通过对模型计算出的 焦距及畸变系数 及上面五个坐标值的分析得出模型的精度较高的结论；二是采用改变特征点数的方法 或 利用“ A ,B ,C 三个标靶的中心的像应在一 条 直线上 ” 验证模型的稳定性。问题 4 采用二目立体视觉模型确定了给出两部固定相机相对位置的数学模型和方法。 本文建立的算法可操作性强，精度较高， 稳定性好，对解决类似问题的计算有一定的推广价值。 关键词： 拟合椭圆 特征点提取 步法 坐标旋转矩阵 公切点 数码相机定位的数学模型 一 问题的 提出 数码相机定位 的数学模型来源于 2008 年全国大学生数学建模竞赛的 A 题。 一般地 在物平面上画若干个圆（称为靶标），它们的圆心就是几何的点了 ， 它 们的像一般会变形 为椭圆， 从靶标上的这些圆的像中把圆心的像精确地找到，标定就可实现。 现设计靶标如下，取 1 个边长为 100正方形，分别以四个顶点（对应为 A、 C、 D、E）为圆心， 12半径作圆。以 上距离 A 点 30的 B 为圆心， 12如图 1 所示，用一位置固定的数码相机摄得其像，如图 2 所示 。 (1) 建立数学模型和算法以确定靶标上圆的圆心在该相机像平面的像坐标 , 这里坐标系原点取在该 相机的光学中心 ， 面平行于像平面； (2) 对由图 2、图 3 分别给出的靶标及其 像，计算靶标上圆的圆心在像平面上的像坐标 , 该相机的像距（即 光学中心 到像平面的距离）是 1577 个像素单位 (1 毫米约为 像素单位 )，相机分辨率为 1024 768； (3) 设计一种方法检验你们的模型，并对方法的精度和稳定性进行讨论； (4) 建立用此靶标给出两部固定相机相对位置的数学模型和方法。图 1 靶标示意图 图 2 靶标的像 二 模型的假设及符号说明 1 模型的假设 (1) 假设题目中给出的图 3的尺寸是实际物理尺寸； (2) 图像坐标系的坐标原点是图 3的中心，称为主心坐标； (3) 相机 不需做任何的运动； (4) 标定物为共面点， 将 靶标所在的面作为世界坐标系的坐标平面； (5) 相机 的有效焦距不变 。 2. 符号说明 (1) R ：表示世界坐标系到 相机 坐标系旋转变换矩阵 ，333231232221131211 (2) T ：表示世界坐标系到 相机 坐标系的平移向量 ， (3) f ：表示 相机 的有效焦距 ， 57 7 ； (4) ),(示世界坐标系下物点 P 的坐标； (5) ),(X：表示图像坐标系下物点 P 有径向畸变 的实际像坐标； (6) ),(X：表示 图像 坐标系下物点 P 的针孔成像的理想坐标； (7) ),( 表示 相机 坐标系下物点 P 的坐标； (8) ),(00 点 坐标，即图像坐标系下坐标原点的像素坐标； (9) ),( 表示计算机坐标 系 下 P 的 像的像素坐标。 (10) k ： 表示径向畸变系数。 (11) 每个 像素在 X 轴与 Y 轴方向上的物理尺寸 。 三 问题的分析 为了确定靶标上圆的圆心在文中给定相机像平面的像坐标， 要把图 3 中的数 据 信息提取出来。 由于图 3 中的数 据 信息是以像素为单位的，为了与图 2 的毫米单位一致，同时便于运算，通过坐标变换 ， 将 计算机图像坐标变换为图像坐标 ， 并用图像边界坐标的信息拟合出 5 个椭圆的方程，从而可以求出椭圆的几何中心。 1 坐标系的 建立 计算过程中 需要建立如下的 四个坐标系 ，如图 4 所示。 (1) 相机坐标系 ：原点 O 定义在相机的光 学中心 ， z 轴与光轴重合； (2) 图像坐标系 1 ：原点 1O （主点）定义为相机光轴与图像平面的交点， 轴与 x ， y 轴平行， 1相机的 有效焦距 57 7 ； (3) 计算机图像坐标系 2 ：原点 2O 位于 u 轴和 v 轴分别与 X 轴和 Y 轴平行。 u 、 v 分别表示该像素在数组中的列数和行数，且以像素为单位 。在 计算机图像坐标系 2 中 ， 主点坐标记为 1O ),(00 假设主点 坐标是图像坐标的中点，本文中可表示为 3842/768,5122/1 0 2 400 (4) 世界坐标系：原点在正方形的 中心 ，C 、 轴的正向为向量的正方向），因此物点 P 的坐标可表示为)0,( X 。 设 ),( 在小孔成像的理想图像坐标中的像点坐标， ),( 实际图像点，考虑畸变的径向摄像机模型图，如图 4所示2 。一般情况下， 头畸变主要为一阶径向畸变，在这里我们只考虑一阶径向畸变。用一个二阶多项式近似： )(1 22)(1 22式中 k 为畸变系数，X ，Y 为畸变坐标。 2 提取图像数据 (1) 提取图 3并以 (2) 在 用 (3) 使用 (4) 使用 令对二值图像矩阵进行分块标识（ 每个椭圆 域标记 为一块，共 5块 ） 。 (5) 对每一块椭圆域用 令提取边缘点，并求出这些边缘点在计算 机图像坐标系下的坐标。 3 计算机图像坐标与图像坐标的转换 如图 4所示，计算机图像坐标系以像素为单位，为便于计算，需将其转换为图像坐标系。设 1 坐标系原点在 坐标系中坐标为 ),(00 个像素在 X 轴与Y 轴方向上的物理尺寸为 则图像中任意一个像素在两个坐标系下的坐标有如下关系： ,( , 1O ),( ,( 0,0( z 图 4 考虑畸变的径向摄像机模型图 00) 用齐次坐标与矩阵形成，将式 (1)表示为 11001001100110000100(2) 已知相机分辨率为 7681024 ，若题中给定图 3尺寸为实际物理尺寸，则 7/10 241 7/7681 方便计算，取 1 像素单位。 设 ),(00 据式 (2)， 可 将像平面五个椭圆边缘点坐标变换到图像坐标系下的坐标。 4 椭圆拟合 1 椭圆拟合 是对 已提取的 图 3的五个椭圆 边缘点 的数据 .,2,1),( ，进行曲线 拟合 。 设 椭圆 曲线的一般表达式为 0),( 22 其椭圆的中心点 ),(式 (3)求得： 2 2 ， )04(42 22 c (3) 为了用数据 .,2,1),( 作曲线拟合，首先设 ，1 22 ，利用 线性 最小二乘 法 拟合 ， 用 程系数 如 表 1所示，椭圆 所示。 表 1 拟合椭圆曲线 方程系数 a b c d e g A 000 000 000D 000 000 椭圆方程标准化 ， 由式 (3)得到 A、 B、 C、 D、 何 中心点 ),( 2所示。 表 2 五个靶标的像的几何中心点的坐标 A B C D E =1577/ 单位： 用椭圆的公切线提取特征点 (1) 引理： 曲线切线的投影仍与曲线的投影相切，而且切点的投影仍为投影的切点。 3 (2)两椭圆公切线的计算 设像平面上任意两个椭圆（基于图像坐标系）的曲线方程为： 065423221 065423221 (4) (5) 其公切线方程为： (6) 将式 (6)代入式 (4)并整理得： 02 (7) 其中， 2321 ， ，6523 图 5 椭圆 根据切线的判别法则： 042 得 ： 065423221 (8) 其中，63251 4 ，52432 24 ，31223 4 ，62544 42 ， 51425 42 ，61246 4 同理，将式 (6)代入式 (5)，可得与式 (8)形式相同的结果： 065423221 (9) 其中，63251 4 ，52432 24 ，31223 4 ，62544 42 ， 51425 42 ，61246 4 求解 (8)、 (9)联立的方程组，可得出四组 ),( 即确立了任意两个椭圆的四条公切线。 对于题中给出的靶标，我们仅提取出如图 6所示的四条外围公切线的投影（如图7所示），并计算出八个 公 切点的坐标。 图 7 对应的 所示。 图 8 图 6 对应的 解 A B w E 6 (1) (2) (9) (10) (11) (3) (12) (4) (14) (13) (6) (5) (8) (16) (7) (15) A B C E D 图 7 四条外围公切线的投影 方程分别为： : y= : y= : y= : y=条 内 围公切线的投影 方程分别为： : y= : y= : y= : y=6 个公切点及其对应像的坐标如表 3 所示。 表 3 16 个公切点及其对应像的坐标 序号 (1) (2) (3) (4) 原坐标 (50) (62) 像坐标 ( 序号 (5) (6) (7) (8) 原坐标 (50， 62) (62， 50) (62， (50， 像坐标 ( 序号 (9) (10) (11) (12) 原坐标 (38) (50) 像坐标 ( 序号 (13) (14) (15) (16) 原坐标 (38， 50) (50， 38) (50， (38， 像坐标 (四、模型的建立与求解 1 问题 (1)的求解算法 为了数学模型和算法以确定靶标上圆的圆心在该相机像平面 像坐标 ， 采用 出的 步法 标定过程， 由于 此方法 考虑 了 镜头畸变， 使精度有了一定的提高。 此 方法 采用径向排列约束关系求出世界的坐标系到相机坐标系得旋 转矩阵R 和 平移向量 T ， 然后求得利用 从而求出世界的坐标系与 相机 坐标系之间的对应关系，同时又求得 镜头畸变 系数 。 设 ),(建立的世界坐标系 下物点 P 的 坐标 ， ),( 示物点 机 坐标系下的 像 坐标。则世界 坐标系到 相机 坐标系之间的变换关系为 ： (10) 其中，333231232221131211示旋转矩阵，示平移向量。 下面采用径向 约束的两步法（ 行计算 ： 第一步 ：计算 R 和yx、渡参数的求解 利用径向排列约束 原理及 公式 (10)得到 ： 232221131211(11) 由于 0 (11)可表示为 22211211 )(12) 即 : 2211211 (13) 若对每一个特征点来说，知道了其空间坐标和相应 于 图像坐标 系 下坐标 ，就 有一个方程 (13)与之对应， 取 5个以上的点 ，式 (12)即 成为超定方程组，利用最小二乘法求解得到 R 和 T 。 利用最小二乘法求解 式 (13)的 具体 方法：取 N 个物点 ),(11 X， ),(,X和N 个对应的图像 坐标 点 ),( 11 X ， , ),( X ，得超定方程组为 ： 其中 222222222111111111 22211211 ， 21 则 X 的最小二乘估计为 T 1 )( 。 记：252143122111 , 。 2) 求出上述过度参数后，就可以根据 当参数 1l ， 2l ， 4l ， 5时： )4( 2/122 25242221 B= 2451 )( 2 否则 2 3） 判断在运用径向平行约束时，包括两平行向量同向和反向 的两种情况，故采用下列方法判断： 先设过渡参数求出 1r ， 2r ， 4r ，5r，此将标定点再投影到图像平面上。计算出对应得图像坐标cx，较两者的符号，如果么则上面求出的其余外部参数都作相应的改变。 4） 根据 算 21221113 1 22222123 1 （如果 22122111 的符号为正，则在23 2213231231 231121131332 2112221133 如果下面计算的 负号，则 32312313 , 符号与上面相反。 4） 过渡 焦距 过渡 径向畸变系数 过渡 平移分量 求解 由 于 323112112 1 323122212 1经 整理可得如下含未知数 K 得线 性 方程组： 22 (14) 式中 1211， 1121，231 ， ，222 。 用最小二乘法联合解 (14)这个线形方程组，即可 得到 解。进一步判断 f 的符号，如果 负，则将 号。进一步解得 ，其余参数符号按上面3)所述作相应得改变。 至此， R 与 第二步： 计算平移 分量 径向 畸变系数 k 由于 R ， f 已知，根据 (14)式可以得到： 22 (15) 式中 1211， 1121，231 ， 222。 用最小二乘法联合解 (15)这个线形方程组，即可 得到 k 的解。 以上即为 步法的完整的求解过程。 这样 通过式 (10)得到 物 点 ),(,( 再 利用关系式 ,/ (16) 得到物点 ),(,( 问题 (1)得到解决。 2 问题（ 2）的求解 (1) 利 用上述方法计算靶标上圆的圆心在像平面上的像坐标 在前面 利用椭圆的公切线 ，已经 找出的 8 个 公切点的 世界坐标和图像坐标 如下： 序号 (1) (2) (3) (4) 原坐标 (50) (62) 像坐标 ( 序号 (5) (6) (7) (8) 原坐标 (50， 62) (62， 50) (62， (50， 像坐标 (把它们 作为特征点 代入 公式 (13)、 (15)得到超定方程组 ，利用 程 求解超定方程组得到 和k 71003 k f 最后利用 公式 和,/ ，得到 五个 靶标的中心点的像在相机坐标系下坐标如 表 4所示。 表 4 利用 8 个特征点求得 五个 靶标的中心点的像在相机坐标系下坐标 A B C D E X =1577/ 单位： 2) 计算结果的分析 f 的值： 通过 模型 (14)计算 f 的 过渡 值为 f ， 而实际的 ， f 值的相对误差： % f 说明了计算方法的可靠性 。 k 的值：通过模型 (14)计算径向畸变系数 k 的值为 71003 k ,说明在此问题的求解过程中可以忽略径向畸变带来的误差。 靶标 中心点的像的坐标：通过对表 2 和表 4 的比较可以看出，靶标的像的几何中心和靶标的像的中心坐标并 不 一致，但误差并不 是 很大，说明如果要求精度不高时，也可以用靶标的像的几何中心代替靶标的像的中心。 3 问题 (3)： 模型的检验 通过上面的计算结果的分析可以看出 算 法的精度较高。 下面用以下两种方法检验上述 模型的稳定性 (1) 利用改变特征点数的方法验证模型的稳定性 在计算靶标上圆的圆心在像平面上的像坐标时，前面利用 8 个公切点作为特征点计算了 结果。现在加大和减少特征点的个数计算模型，如果加大和减少特征点对计算结果的影响不大说明，我们建立的模型是稳定的。 通过 对靶标上的 , 的 16 个外公切点作为特征点，计算靶标上圆的圆心在像平面上的像坐标及相应的参数 k ， 8 7 0 1 6 4 2 0 8 0 7 和 k 71053 k f 表 5 利用 16 个特征点求得 五个 靶标的中心点的像在相机坐标系下坐标 A B C D E X =1577/ 单位： 算结果与 利用 8 个特征点求得的结果比较更精确，但差别并不大，说明我们的算法的稳定性。 (2) 利用成像原理检验 引理： 一 条空间直线的投影到任意 一个不与 之 垂直的 平面上像仍是一条直线 。 已知 A ,B ,C 三个标靶的中心在一 条 直线上，从而其像的三个中心坐标 , 应该在一 条 直线上。若利用我们的模型求得的结果， , 中心在一 条 直线或近似一 条 直线，说明模型较稳定。 检验方法：取 , 心的三个坐标 ),(ii 3,2,1i ，拟合直线 ， 31 2)(i 很小，说明算法稳定。 算例： 利用 16 个 特征点 ，根据 模型 (13)求得 的 , ,E 中心的像坐标如表 5所示， 即 有： A ( B ( C (, 拟合方程为： 4 05 4 残差 2810553 。 从图 9 中可以看出 , 点几乎构成一条直线。 同理对 特征点数为 8 的情形， 拟合直线如图 10 所示。 对应的 , 拟合方程为： 5 42 5 残差 2810068 。 4 问题（ 4）的模型 (1) 二目立体 视觉模型 二目立体测试系统利用两台不同角度和位置安放的相机来模拟人眼的功能，获取图 9 对 16个特征点的情形三点拟合的直线 图 10 对 8个特征点的情形三点拟合的直线 像点与相机的三维空间中物体的几何信息。未来利用二目立体测量系统获得三维几何信息，需要是两个相机在一个统一的坐标系下工作，因此不仅要得到单个相机的相机模型，还需要精确确定两个相机的相互位置关系，即建立精确的立体视觉模型。 立体视觉模型以单相机模型为基础来表示两个相机的位置几何关系。如果左相机的外部参数为lR，相机的外部参数为 则两个相机之间的几何关系表示如下（以左相机为基准）： 1 (2) 双目视觉特征点三维重建算法 利用基于视差原理的双目视觉测量原理。对一般情况进行考虑， 如图 11 所示 ，假定左相机 位于世界坐标系的原点处且无旋转，图像坐标系为，焦距为相机坐标系为 ，图像坐标系为 ，焦距为 根据相机透视变换模型有 与 之间可通过空间转换矩阵 M 表示为： 11 333231232221131211 与 之间的旋转矩阵和原点之间的平移变量矢量。 对于 坐标系中的空间点，两相机像平面点的对应关系为： P 图 11 1O y O 11/1 333231232221131211由此，空间点三维坐标可表示为： ll ll )()3()(131211931 )()( )( 232322333231 当已知 两 相机焦距空间低昂在左右相机中的图像坐标时，通过标定出的旋转矩阵 R 和平移矢量 T 就可以得到被测物体点的三维空间坐标。 五、模型的评价及推广 优点：算法简单，适合用计算机编写程序计算，误 差积累少，精度较高，稳定性好，算法可操作性强。 缺点：对算法精度的分析只是利用数值分析得直观方法，缺少理论依据，有待进一步分析。 鉴于以上特点，本文建立的模型对解决类似问题的计算有一定的推广价值。 参考文献 1 张广军， 机器视觉，北京 ： 科学出版社， 2005； 2 谭晓军、余志、李军， 一 种 改 进 的立体 摄 像机 标 定方法 ，测绘学报，第 35 卷第2期， 3 ； 4 苏金明，王永利， 京 ： 电子工业出版社， 5 黄宣国，空间解析几何与微分几何，上海 ： 复旦大学出版社， 6 胡培成等，一种改进的基于圆环点的摄像机自标定方法，光电工程，第 34卷 12期； 7 邹凤娇，苏显渝，李美菊，基于共面点的摄像机线性标定法，光电工程，第 32卷第 4期， 8 何文章等，数学建模与实验 ，哈尔滨：哈尔滨工程大学出版社， 2002 附录 主程序，文件名： m %启动计时器 %= I=; J=); K=,100/255); K=K; I1,); J %* : %依次对块进行处理 %= r,c = 1= %寻找块内点的坐标 n1,r); n,m=); M=n,m); i=1:(r(i),c(i)=1; %构建单个椭圆的图像 = M1=,; %寻找边界线 %1) y,x=1=1); %= y=(n/2)/ %坐标变换 x=()/%= bi,r,s=x,y); %椭圆拟合 )=s; %r,%= %存储拟合方程的系数 %系数顺序： +dx+ey+f(=1)=0 XS()=%= xc, %求椭圆中心点 %= if 1 %存储拟合椭圆的中心点坐标 )=)= rp=x)x); %绘制拟合椭圆 x):x)+); %确定横轴的范围和步长 xn, %= if 1 %存储各椭圆边的拟合点坐标 g=():(g)=xn( g)= g)= = x,y,+,xc,+,xn,o,xn,o); %= rp xc yc r c n1 m2 n m M M1 xn J K I1 s %清理内存 %= 拟合椭圆方程的系数 ) XS 拟合椭圆的中心点坐标 ) %= %求椭圆 A,B,C 几何中心的拟合直线 x=) ) ); y=) ) ); n=y); X=x,n,1); Y=y; bi,r,X); 椭圆 A,B,