2008高教社杯全国大学生数学建模竞赛
承 诺
我们仔细阅读了中国大学生数学建模竞赛的竞赛规则.
我们完全明白,在竞赛开始后参赛队员不能以任何方式(包括电话、电子邮件、网上咨
询等)与队外的任何人(包括指导教师)研究、讨论与赛题有关的问题。
我们知道,抄袭别人的成果是违反竞赛规则的, 如果引用别人的成果或其他公开的资料
(包括网上查到的资料),必须按照规定的参考文献的表述方式在正文引用处和参考文献中
明确列出。
我们郑重承诺,严格遵守竞赛规则,以保证竞赛的公正、公平性。如有违反竞赛规则的
行为,我们将受到严肃处理。
我们参赛选择的题号是(从 A/B/C/D中选择一项填写): A
我们的参赛报名号为(如果赛区设置报名号的话): 2108373
所属学校(请填写完整的全名): 南京大学
参赛队员 (打印并签名) :1. 笪庆
2. 周超
3. 俞庆进
指导教师或指导教师组负责人 (打印并签名):
日期: 2008 年 9 月 22 日
赛区评阅编号(由赛区组委会评阅前进行编号):
数学中国
共享
www.madio.cn
2008高教社杯全国大学生数学建模竞赛
编 号 专 用 页
赛区评阅编号(由赛区组委会评阅前进行编号):
赛区评阅
(可供赛区评阅时使用):
评
阅
人
评
分
备
注
全国统一编号(由赛区组委会送交全国前编号):
全国评阅编号(由全国组委会评阅前进行编号):
数学中国论文共享
www.madio.cn
1
数码相机定位
摘要
本文假设数码相机成像原理为小孔成像,在此基础上,通过两种合理的模型
对数码相机定位问题进行了较深入的研究。
针对问题一和二,我们建立了两种不同模型——变换矩阵模型和公切线模
型。在变换矩阵模型中,建立了物、像、相机三个坐标系,分别称为世界坐标系,
像坐标系和光心坐标系。研究世界坐标系向像坐标系的变换矩阵 3 4( )ijM a ´= ,推
导出圆在像坐标系中的像为椭圆。利用灰度检测可以得到像中各椭圆圆周上各点
的坐标,通过多元线性回归拟合出各椭圆方程;对单独一个圆进行研究时,在合
理的近似前提下,以圆心为世界坐标系的原点,可求出该圆心所成像在像坐标系
中的坐标 14 24,u a v a= = 。最后我们求得 5个圆心所成像在光心坐标系中的坐标分
别为(单位:mm):(-50.00,51.32,-417.20)、(-23.54,49.47,-417.20)、(33.86,
45.24,-417.20)、(-60.05,-31.22,-417.20)、(18.52,-31.48,-417.20)。在公
切线模型中,通过简单几何证明,得出在小孔成像时,公切线交点的像就是公切
线像的交点,联系题目中所给标靶的特殊性(所有圆全等),得出像平面中公切
线交点连线的交点就是标靶中对应圆心的像,并
了一种算法得到 5个圆心所
成像在光心坐标系中的坐标分别为(单位:mm):(-49.92,51.36,-417.20)、
(-23.47,49.34,-417.20)、(33.88,45.05,-417.20)、(-60.04,-31.29,-417.20)、
(18.58,-31.56,-417.20)。
在问题三中,我们用计算机模拟的方法,统计和分析了我们模型的在不同的
情况下所得到的结果与理论值之间的误差,并着重研究了相机与标靶的距离和像
平面与圆平面之间的偏角对结果的影响。结果表明在一定的前提下,当相机与标
靶的距离大于 200毫米, 0.5 0.5a- £ £ 以及 1 0.5b- £ £ (单位为弧度)时,我们
的结果与理论值相差不到一个像素,有着较好的稳定性和精度。
问题四中,通过每个相机旋转变换矩阵 R和平移向量 T,可以得到两相机的
变换关系: 1 11 2 1 2 2 1R R R T R R T T
- -= = +、 ,即相对位置关系,并理论推导了从两相
机中像在光心坐标系中的参数得到物在世界坐标中的参数,实现双目定位。另外
在相机的光心和像屏中心的连线垂直于象平面基础上,我们还给出另外一种合理
模型,通过矢量的方法求出物相对于光心坐标系的精确位置,从而可以得到两相
机的相对位置。
关键词:相机定位、小孔成像、变换矩阵、公切线、计算机模拟
数学中国论文共享
www.madio.cn
2
一、问题重述
数码相机定位在交通监管(电子警察)等方面有广泛的应用。所谓数码相机
定位是指用数码相机摄制物体的相片确定物体表面某些特征点的位置。最常用的
定位方法是双目定位,即用两部相机来定位。对物体上一个特征点,用两部固定
于不同位置的相机摄得物体的像,分别获得该点在两部相机像平面上的坐标。只
要知道两部相机精确的相对位置,就可用几何的方法得到该特征点在固定一部相
机的坐标系中的坐标,即确定了特征点的位置。于是对双目定位,精确地确定两
部相机的相对位置就是关键,这一过程称为系统标定。
标定的一种做法是:在一块平板上画若干个点,同时用这两部相机照相,分
别得到这些点在它们像平面上的像点,利用这两组像点的几何关系就可以得到这
两部相机的相对位置。然而,无论在物平面或像平面上我们都无法直接得到没有
几何尺寸的“点”。实际的做法是在物平面上画若干个圆(称为靶标),它们的圆
心就是几何的点了。而像平面上的园一般会变形,所以必须从靶标上的这些圆的
像中把圆心的像精确地找到,标定就可实现。
有人设计靶标如下,取 1个边长为 100mm的正方形,分别以四个顶点(对
应为 A、C、D、E)为圆心,12mm为半径作圆。以 AC边上距离 A点 30mm处
的 B为圆心,12mm为半径作圆。用一位置固定的数码相机摄得其像。利用所得
图像, 具体解决如下几个问题:
(1) 建立数学模型和算法以确定靶标上圆的圆心在该相机像平面的像坐标,
这里坐标系原点取在该相机的焦点,x-y平面平行于像平面;
(2) 对由图 2、图 3分别给出的靶标及其像,计算靶标上圆的圆心在像平面上
的像坐标, 该相机的像距(即焦点到像平面的距离)是 1577个像素单位(1
毫米约为 3.78个像素单位),相机分辨率为 1024×786;
(3) 设计一种方法检验你们的模型,并对方法的精度和稳定性进行讨论;
(4) 建立用此靶标给出两部固定相机相对位置的数学模型和方法。
二、模型假设
a.本题中数码相机成像系统看成是小孔成像;
b.灰度检测前将图像改成黑白图的误差不予考虑。
三、符号说明
WO :世界坐标系
CO :光心坐标系
PO :像坐标系
R:旋转矩阵
T :平移向量
M :空间变换矩阵
数学中国论文共享
www.madio.cn
3
f :焦距(mm)
L:每毫米的像素单位
P:图像矩阵
kB :第 k个圆的边缘点集
四、模型的建立与求解
引理:
我们针对小孔模型提出一些内在性质并给予简单证明。
图 1
性质 1:直线经小孔后所得到的像仍为直线。
证明:如图,a g、 分别为物平面和像平面,b为与像平面平行的平面,点 O为
小孔,AC为a上的线段,B为 AC上任一点,经小孔 O后得 g上直线 2 2A C 。由
AC和 O可以确定一个平面 OAC,所以 2 2A C 为平面 OAC和平面 g的交线;设 B
经小孔成像到g上的 2B 点。因为 2B 同时在平面 OAC 和平面g上,所以 2B 必在
两平面交线即 2 2A C 上。因为 B为 AC上任一点,所以直线 AC经小孔 O成像为 g
平面上的直线 2 2A C ,故直线经小孔后所得到的像仍为直线。
性质 2:线段中点经小孔成像所得点不一定还是像线段(经小孔成像后所得到的
线段)的中点。
证明:如图, 1 1AC 为平面 OAC与 b平面的交线, 1B为 B所对应点。因为 b平行
于g,平面 OAC分别交 b、g与 1 1AC 、 2 2A C ,则 1 1AC 与 2 2A C 相互平行。
当a与 b不重合即q不为 0时,如 B为 AC中点,因为 1 1CC AA与 相交,则 1B必然
数学中国论文共享
www.madio.cn
4
不是 1 1AC 中点,所以 2B 也必然不是 2 2A C 中点,也就是说线段中点经小孔成像所
得点不一定还是像线段的中点,仅当该线段平行于像平面时才仍是中点。
性质 3:两直线交点的像仍是两直线像的交点。
证明:如下图:a 平面上两直线 AB与CD交于点E,经小孔O两直线分别得到 b平
面上的像 1 1 1 1A B C D和 ,点 E在 b平面上的像为 1E ,则平面 1 1 1 1ABB A CDD C与 交于
直线OE,平面 1 1 1 1ABB A CDD C与 分别交平面 b于直线 1 1 1 1A B C D与 ,则 1E 同时在平
面 1 1 1 1ABB A CDD C b、 与 上,则 1E 同时在 1 1 1 1A B C D和 上,即 b平面上, 1 1 1 1A B C D和
交于 1E 点,两直线交点的像仍是两直线像的交点。
我们将此结论推广到曲线的情形,很显然结论也必然成立。只要在曲线相交处对
两曲线取无限小段,那么就可以看成是直线相交的情形。
图二
性质 4:圆经小孔成像为椭圆.
图三
数学中国论文共享
www.madio.cn
5
证明:如图所示,圆O经小孔成像如图。在一个与圆O所在平面平行的平面上所
成像仍然是圆,如图中圆 2O (这点很容易通过以上三个性质得到)。那么对于在
与圆O所在平面不平行的平面上的像,则可以看成是用一个平面去截图中的圆锥
(当然图中圆锥可以无限的延长),根据圆锥曲线的定义可知,所截出来的图形
就是椭圆。且可以看出圆心O在截面上所成的像 1O并不是椭圆中心M (只有截
面与圆平面平行时才是椭圆中心),所以圆经小孔成像为椭圆得证。
性质 5:圆的某一条切线的切点的像,仍然是椭圆的切线,而且切点的像就是椭
圆的切点。
证明:由性质 3的推广可知,任意两条曲线的交点的像就是他们两个像的交点。
因为圆中切线与圆是相交于一点的,那么像中切线的像与圆的像(椭圆)至
少也会有个交点。假设圆的切线的像不再是像当中椭圆的切线,则切线的像与椭
圆必有两个交点。根据光路可逆原理,我们可以把原来的圆看成是像(椭圆)经
小孔所得到的像。那么对于另外一个交点,它的原像也必然是原像平面上两曲线
的交点,则原像平面中处切点外还有另外一交点,产生矛盾,故假设不成立,所
以结论得证。
模型一:变换矩阵模型
问题一
在标靶上,以某个圆的圆心为原点建立空间直角坐标系,由 W WX ,Y , WZ 轴组
成,称其为世界坐标系;在像空间上建立像坐标系,由u v、 轴组成;由于摄像
机可以安放在环境中的任何一个位置,我们也建立一个坐标系来描述,由
C C CX Y Z、 、 轴组成,原点位于光心,称其为光心坐标系。光心坐标系与世界坐
标系间的转换可以同过旋转矩阵 R和平移矩阵 T来实现。如空间一点 P在世界坐
标系和光心坐标系中的坐标分别为 ( )( ) TTW W W C C CX Y Z X Y Z 、 ,于是存在关系:
[ ] [ ]C 1[ 1] 1 1 ...........(0)0 1
T TT
C C W W W W W WT
R T
X Y Z X Y Z M X Y Zé ù = = ê ú
ë û
其中 R为 3×3的正交单位矩阵, ( )0 0 0 0 TT = , 1M 为 4×4矩阵,0T和 1的
加入只是为了方便以后的计算。空间点 p的像在像坐标系的位置与 p在光心坐标
系中的关系如图可得:
,C C
C C
fX fYx y
Z Z
= =
其中 ( , )x y 为点 p的像在像坐标系的坐标,写成矩阵形式就是:
数学中国论文共享
www.madio.cn
6
0 0 0
0 0 0 ....................................(1)
1 0 0 1 0
1
C
C
C
C
X
x f
Y
Z y f
Z
é ù
é ù é ù ê ú
ê ú ê ú ê ú=ê ú ê ú ê ú
ê ú ê ú ê úë û ë û
ë û
所以可以得到:
0 0 0
,
0 0 0 ...............................(2)
0 ,1
1 0 0 1 0
1
W
W
C T
W
X
x f
R T Y
Z y f
Z
é ù
é ù é ù ê úé ùê ú ê ú ê ú= ê úê ú ê ú ê úë ûê ú ê ú ê úë û ë û
ë û
参照附录中对于旋转矩阵的说明,我们可以得到:
sin
sin cos sin sin cos
cos cos
W
C W Z W W Z
W
X
Z Y T X Y T
Z
b
a b b a b
a b
é ù é ù
ê ú ê ú= - - = - -ê ú ê ú
ê ú ê úë û ë û
所以由:
0
0 0 0
,
0 0 0
0 ,1
1 0 0 1 0
1 1
W W
W W
C T
W W
X X
x f
R T Y Y
Z y f M
Z Z
é ù é ù
é ù é ù ê ú ê úé ùê ú ê ú ê ú ê ú= =ê úê ú ê ú ê ú ê úë ûê ú ê ú ê ú ê úë û ë û
ë û ë û
可得:
11 12 14
21 22 24
( sin sin cos )
( sin sin cos )
W W Z W W
W W Z W W
X Y T u a X a Y a
X Y T v a X a Y a
b a b
b a b
- - = + +
- - = + +
可以通过上式解得:
22 14 14 22 12 24 12 24
22 11 11 22 12 21 21 12
21 14 14 21 11 24 11 24
22 11 1
( sin cos ) ( sin cos )
( sin sin cos ) ( sin sin cos )
( sin ) ( sin )
( sin sin cos
Z Z
W
Z Z
W
uT a va a a vT a ua a aX
ua va a a va ua a a
uT a va a a vT a ua a aY
ua va a
a b a b
b a b b a b
b b
b a b
+ + - + +
=
- - - - -
- + - - - + -
=
- - 1 22 12 21 21 12) ( sin sin cos )a va ua a ab a b- - -
因为 W WX Y、 为圆周上的点,所以在世界坐标系满足:
2 2 2W WX Y r+ =
带入可得到二次曲线方程:
2
22 24 12 14 14 22 12 24
2
21 24 14 11 14 21 11 24
2 2
22 21 11 12 11 22 21 12
[( sin cos ) ( sin cos ) ( )]
[ ( sin ) ( sin ) ( )]
[( sin sin cos ) ( sin cos sin ) ( )]
Z Z
Z Z
T a a u T a a v a a a a
T a a u a T a v a a a a
r a a u a a v a a a a
a b a b
b b
b a b a b b
- - - + +
+ - + + + - +
= + - + - -
其代表一椭圆。
注意到此处的 CZ 是标靶上圆上点在光心坐标系中 CZ 方向上的坐标。对于本
数学中国论文共享
www.madio.cn
7
题来说,因为标靶在光心坐标系中 CZ 方向上的坐标应该大于相机的两倍焦距,
即应在 1米附近;而对于标靶上一个圆,其半径仅为 12mm,当标靶平面与光心
坐标系中 C C CX O Y 平面存在一定夹角q时,那么一个圆上所有点在 CZ 上坐标的差
异只是12sinq,与 1m相比较而言,其误差是比较小的,所以我们可以近似认为
对于标靶上同一个圆上的点,其 CZ 是相同的(后面我们将来讨论这种近似所带
来的误差,会发现其误差是非常小的,可见这种近似的合理性),
所以在此情况下,我们认为 CZ 是不变的值,于是有:
0 0 0
,1 0 0 0 .........................(3)
0 ,1
1 0 0 1 0
1
W
W
T
C W
X
x f
R T Y
y f
Z Z
é ù
é ù é ù ê úé ùê ú ê ú ê ú= ê úê ú ê ú ê úë ûê ú ê ú ê úë û ë û
ë û
由于在处理像平面上的点时,我们经常用的单位为像素,对像平面坐标为
( , )x y 的点,改用像素为单位后,其坐标为:
u xL
v yL
=ì
í =î
其中 L为每毫米的像素单位。于是:
0 0 0 0 0 0 0
,10 0 0 0 0 0 0 ...(5)
0 ,1
1 0 0 1 1 0 0 1 0 0 1 0
1 1
W W
W W
T
C W W
X X
u L x L f
R T Y Y
v L y L f M
Z Z Z
é ù é ù
é ù é ù é ù é ù é ù ê ú ê úé ùê ú ê ú ê ú ê ú ê ú ê ú ê ú= = =ê úê ú ê ú ê ú ê ú ê ú ê ú ê úë ûê ú ê ú ê ú ê ú ê ú ê ú ê úë û ë û ë û ë û ë û
ë û ë û
其中M 为一 3×4的坐标变换矩阵,而且在相机本身内部参数和相对标靶的
位置不变时,M 是一个常数矩阵。由上式可以看出对于世界坐标系中的任意一
点 ( , , )W W WX Y Z ,经坐标变换矩阵M 变换,便可以得到其像在像坐标系中的坐标
( , )u v 。
针对题目所给信息,为简化模型,我们可分别随标靶中每个圆分别讨论(因
为对每个圆讨论的方法完全相同,故本文只详细讨论一个圆,其他的可完全类
比)。对某个圆讨论时,取该圆所处坐标系为世界坐标系,坐标原点取在圆心处
(后面会发现这样选取的精妙之处),所得像处于像坐标系,数码相机处于摄像
坐标系。可以很直觉的发现圆周上的点是至关重要的,那么我们先来探讨圆周上
的各点及其所成像的位置。
因为在圆周上,注意到此时的世界坐标系里只需其二维情形,取 ZW=0,故其
曲线方程为:
2 2 2W WX Y r+ =
数学中国论文共享
www.madio.cn
8
因为变换矩阵M 为 3×4矩阵,可设
11 12 13 14
21 22 23 24
31 32 33 34
a a a a
M a a a a
a a a a
é ù
ê ú= ê ú
ê úë û
所以圆周上任意一点 ( , , )W W WX Y Z ,经坐标变换矩阵M 变换后得到像在像坐
标系的坐标:
11 12 13 14
21 22 23 24
31 32 33 34
..............(6)
0 0
1
1 1
W W
W W
X X
u a a a a
Y Y
v M a a a a
a a a a
é ù é ù
é ù é ùê ú ê ú
ê ú ê úê ú ê ú= =ê ú ê úê ú ê ú
ê ú ê úê ú ê úë û ë û
ë û ë û
有:
11 12 14
21 22 24
...................(7)W W
W W
u a X a Y a
v a X a Y a
= + +ì
í = + +î
与圆的曲线方程联立可得像的曲线方程:
2 2
222 14 12 24 12 14 11 24
11 22 21 12 12 21 22 11
( ) ( ) ( ) ( )a u a a v a a u a a v a r
a a a a a a a a
é ù é ù- - - - - -
+ =ê ú ê ú- -ë û ë û
从上式不能很显然的发现什么性质,所以我们将其展开:
2 2 2 2 2 2 2 2
21 22 11 12 11 21 12 22 14 21 22 24 11 21 12 22
2 2 2 2 2 2 2 2
24 11 12 14 11 21 12 22 21 22 14 11 12 24 11 21
2 2
12 22 14 24 12 21 22 11
( ) ( ) 2( ) [2 ( ) ( )]
2[ ( ) ( )] [( ) ( ) 2(
) ( ) ] 0
a a u a a v a a a a uv a a a a a a a a u
a a a a a a a a v a a a a a a a a
a a a a a a a a r
+ + + - + - + + +
- + + + + + + + - +
- - = ......(8)
简化上式为:
2 21 2 3 4 5 6 0..........................(9)k k v k u k uv k v k u+ + + + + =
可见这是一个一般二次曲线方程,故其图像为一椭圆(当然也可能会是一条
直线,暂不考虑这种情况)。其中 1k 至 6k 分别为上式的系数。观察(8)式可以发
现这个重要信息:u v、 都是 14 24a a、 和
2 2u v uv、 、 的系数组成,而 14 24a a、 正是世
界坐标系的圆心经转换矩阵得到的像坐标系的坐标。所以其简化式可变为:
2 21 24 5 14 4 14 6 24 4 4 5 6(2 ) (2 ) 0.......(10)k a k a k v a k a k u k uv k v k u- + - + + + + =
我们所需要求的是圆心在像坐标系的像坐标的坐标,由于我们将世界坐标系
的原点放在圆心,故圆心在世界坐标系中的坐标为 (0,0),根据(7)式,在像坐
标系的像坐标的坐标为:
14
24
o
o
u a
v a
=ì
í =î
数学中国论文共享
www.madio.cn
9
于是我们的问题转化为求解 14 24a a、 ,而由我们(10)式的分析,我们可以
得到:
24 5 14 4 2
14 6 24 4 3
(2 )
.............(11)
(2 )
a k a k k
a k a k k
- + =ì
í- + =î
上式为关于 14 24,a a 的二元一次方程组,于是问题又转化为求出 (2 6)ik i£ £ ,
我们知道(9)式是圆的像在像平面的曲线方程,而 ik 为方程的系数,下面我们
来求解这个方程。
我们的数据来源于数码相机拍摄的照片,所以我们首先需要从照片中获取数
据。我们知道,现代的数码照片的存储原理基于点阵。例如在一张灰度照片中存
储了n m´ 的矩阵 P,矩阵中的每个元素 ( , )P i j 记录的是一个灰度值,表示这张
照片在该点的明暗程度(而在一张彩色照片里每个点存储的是一组 RGB)。所以
我们的直接数据是一个768 1024´ 的矩阵,下面我们需要在这个矩阵中提取圆的
像的边缘以便求解该圆的像的曲线方程。
传统的图像边缘检测的算法有基于灰度直方图的边缘检测、基于梯度的边缘
检测、 Laplacan 边缘算子、Canny 边缘检测算子、模糊推理的边缘检测和
Mallat小波边缘检测算子等,不同的算法有着不一样的应用背景,而在本文中
所需处理的照片经过简单的处理后可以转化为单色的照片,所以用最简单的基于
灰度直方图的边缘检测算法便可得出较理想的效果。
图四 经过黑白处理后的照片
经过黑白处理后,所有点的灰度值只有 0(黑色)和 255(白色)两种,我
们用下面的算法进行边缘检测(对圆 B C E DA、 、 、 、 依次编号为圆 1、2、3、4、
5):
边缘检测算法:
Step 1 设灰度阀值 T,当某点的灰度小于 T时,表示其为黑色,反之则
为白色。在本文中我们取 T=128;
数学中国论文共享
www.madio.cn
10
Step 2 为每个圆的像估算一个矩形邻域,使得该矩形邻域可以完全包含
该圆的像且仅包含该圆的像;
Step 3 对第 k个圆的像,执行第 4-6步(k=1,2,3,4,5);
Step 4 对矩形邻域进行逐行扫面,当出现 ( 1, )P i j T- ³ 且 ( , )P i j T<
或者 ( , )P i j T< 且 ( 1, )P i j T+ ³ 时将 ( , )P i j 加入该圆的像的边缘
点集 k
B
;
Step 5 对矩形邻域进行逐列扫面,当出现 ( , 1)P i j T- ³ 且 ( , )P i j T<
或者 ( , )P i j T< 且 ( , 1)P i j T+ ³ 时将 ( , )P i j 加入该圆的像的边缘
点集 kB ;
Step 6 除去 kB 中重复的点。
我们得到的点集 {( , ) | ( , ) i }k i i i iB u v i Z u v
+= Î 且 位与第个圆的像的边缘 ,这里
,i iu v 的单位为像素,像坐标系以左上角为原点,图片的正下方为u轴,正右方为
v轴。
对于每一个圆,我们将其所有边缘点带入其曲线方程:
2 2
1 2 3 4 5 6 0, 1 | | .........(12)i i i i i i kk v k u k u v k v k u k i B+ + + + + = £ £
我们得到了关于 (1 6)ik i£ £ 的线性齐次方程组,而且在这个方程组中,方程
的个数远远大于未知数的个数,而且我们发现若直接求解往往只能得到零解,于
是我们将其化为多元线性回归模型求解。
由上文讨论,像的曲线为椭圆,故平方项系数不为零,故原方程可化为:
2 23 51 2 4
6 6 6 6 6
0, 1 | | ..........(13)i i i i i i k
k kk k kv u u v v u i B
k k k k k
+ + + + + = £ £
令
6
1 5jj
k
b j
k
= , £ £ ,则有:
2 2
1 2 3 4 5 , 1 | | ...........(14)i i i i i i ku b v b u b u v b v b i B- = + + + + £ £
上式是典型的多元线性回归模型,应用最小二乘法可以求解。当我们解出了
1 2 3 4 5, , , ,b b b b b 后,由(11)式得:
24 5 14 4 2
14 24 4 3
(2 )
(2 )
a b a b b
a a b b
- + =ì
í- + =î
数学中国论文共享
www.madio.cn
11
v
u
于是我们得到了圆心的像在uov平面上的坐标:
3 5 2 4
14 2
4 5
2 3 4
24 2
4 5
2
4
.......................(15)
2
4
o
o
b b b bu a
b b
b b bv a
b b
-ì = =ï -ï
í -ï = =
ï -î
最后,我们需要将坐标转换到题目中给定的坐标系中,(x与 v同向,y
与 u反向)即:
( )
( 512) /
(384 ) / .......................... 16
/
o o
o o
o
x v L
y u L
z f L
= -ì
ï = -í
ï = -î
问题二
首先,我们用 matlab实现边缘检测算法,我们得到了每个圆的像的边缘点
集 kB ,如下图:
图五
然后我们对每个圆的像的边缘点用多元线性回归去求该圆的像的曲线方程
(置信度取 0.05),如下表所示:
表 1 第 1个圆的像的曲线方程参数
参数 参数估计值 参数置信区间
b1 140397.20023 [139748.34930 141046.05115]
b2 -640.44473 [-644.04310 -636.84637]
b3 -407.58657 [ -409.32833 -405.84480]
b4 0.08784 [ 0.08246 0.09323]
b5 0.96592 [ 0.96064 0.97120]
R2=0.99999393 , F= 11279037.31254 , p=0
数学中国论文共享
www.madio.cn
12
表 2 第 2个圆的像的曲线方程参数
参数 参数估计值 参数置信区间
b1 227111.11949 [ 225858.66719 228363.57179]
b2 -869.84886 [-875.31270 -864.38502]
b3 -453.55647 [-456.21530 -450.89765]
b4 0.14049 [ 0.13421 0.14677]
b5 0.99545 [ 0.98926 1.00165]
R
2
=0.999993179 , F=9786431.732757331 , p=0
表 3 第 3个圆的像的曲线方程参数
参数 参数估计值 参数置信区间
b1 521021.19332 [ 517789.67136 524252.71528]
b2 -1431.47741 [ -1441.06500 -1421.88982]
b3 -603.80738 [-608.27782 -599.33694]
b4 0.27703 [ 0.27005 0.28402]
b5 1.07232 [1.06508 1.07957]
R2=0.99999442 , F=11019117.97897741 , p=0
表 4 第 4个圆的像的曲线方程参数
参数 参数估计值 参数置信区间
b1 343623.44295 [ 342627.40480 344619.48111]
b2 -557.42512 [-561.92415 -552.92609]
b3 -1057.92607 [-1059.63303 -1056.21910]
b4 0.19027 [0.18428 0.19625]
b5 0.81131 [0.80593 0.81670]
R
2
=0.99999918 , F=75335174.25731819 , p=0
表 5 第 5个圆的像的曲线方程参数
参数 参数估计值 参数置信区间
b1 660749.59267 [ 657278.55376 6.64220.63157]
b2 -1228.17779 [-1237.71197 -1218.64362]
b3 -120816690 [-1212.41769 -1203.91612]
b4 0.34664 [0.33934 0.35393]
b5 0.90414 [0.89715 0.91114]
R2=0.99999914 , F=65011077.78295716 , p=0
在上表最后一行中,R2是回归方程的决定系数,F 为 F 统计量值,p 为与 F
统计量对应的概率值。从数据看来,R2都达到了 99.999%以上,表明 2u- 的 99.999%
以上的可以由我们求出来的方程确定,F 值远远超过 F检验的临界值,p更是远
小于置信度 0.05,所以我们求出的曲线方程是高度精确的。
将以上数据代入(15)式,我们便得到了五个圆心的像的坐标:
数学中国论文共享
www.madio.cn
13
表 6 圆心像的坐标
圆 ou ov ( )ou 像素 ov(像素)
1 189.61023 322.89682 190 323
2 197.06280 423.00273 197 423
3 213.26380 639.91289 213 640
4 501.87983 284.67960 502 285
5 503.07979 582.75218 503 582
用(16)式将其转换到题目所给的坐标系:
表 7 光心坐标系下圆心像的坐标
圆 ( )ox mm ( )oy mm ( )oz mm
1 -50.00000 51.32275 -417.19577
2 -23.54497 49.47089 -417.19577
3 33.86243 45.23809 -417.19577
4 -60.05291 -31.21693 -417.19577
5 18.51851 -31.48148 -417.19577
以上我们就求解出标靶上各圆圆心的像在像坐标系下的具体坐标。
模型二:公切线模型
图六 图七
如图所示,为标靶平面与像平面的示意图。图中直线为公切线。物平面内,圆 1O
与 2 3O O、 的公切点分别为 D AB C、 和 、 ,对应像平面中为椭圆 1 'O 与椭圆
2 3' 'O O、 的公切点分别为 ' B' C' A'D 、 和 、 ,由本文所给出的性质可以很容易得出
' B' C' A'D 、 、 、 分别为 B C AD、 、 、 的像点。
数学中国论文共享
www.madio.cn
14
图八 图九
在标靶平面,显然四边形 PQRS为正方形,则 PR QS与 连线交于圆心 1O,由引理
中性质 3 可知 ' ' ' 'P Q R S 分别为 PQRS的像,则 ' ' ' 'P R Q S与 也分别为 PR QS与 的
像,则圆心 1O 在像平面所对应的像为 ' ' ' 'P R Q S与 的交点 1 'O (注意这里 1 'O 并不
一定就是该椭圆中心)。
两椭圆公切线求解算法如下:
Step 1:由模型一中灰度检测,我们已经得到像平面上椭圆圆周上各点坐标;
Step 2:在两椭圆上各任取一点,连成直线 l,固定其中一点(静点),扫描另一
点(动点)所在椭圆的圆周上各点,每扫描到一点在直线 l上方(或下方,对应
的是另一条切线),用该点取代动点成为新的动点,此时直线 l也变成静点和此新
的动点的连线,扫描一周后停止。
Step 3:固定此时的动点(这时就称之为静点),扫描另一椭圆圆周上各点,当该
点在直线 l上方时,取代动点成为新的动点,此时直线 l也变成静点和新的动点
的连线,扫描一周后停止;
Step 4:循环 step 2和 step 3,直至两椭圆圆周上都不再存在位于直线 l上方的点,
此时的 l就是两椭圆的公切线。
对该算法而言,收敛速度很快,一次迭代后就几乎接近了切点,所以可以用
来高效地计算切点。
我们按照上述思想编写了程序,所得的切线图如下:
图十
数学中国论文共享
www.madio.cn
15
在我们的程序中,我们还计算了切线之间的交点及其交点连线的交点,即上
文所分析的圆心的像,结果如下:
表 8 圆心像的坐标
圆 ou ov ( )ou 像素 ov(像素)
1 189.87734 323.30371 190 323
2 197.51240 423.29359 198 423
3 213.70159 640.08497 214 640
4 502.26962 285.03296 502 285
5 503.31215 582.23246 503 582
表 9 光心坐标系下圆心像的坐标
圆 ( )ox mm ( )oy mm ( )oz mm
1 -49.91965 51.35520 -417.19577
2 -23.46730 49.33535 -417.19577
3 33.88491 45.052489 -417.19577
4 -60.04419 -31.28799 -417.19577
5 18.58001 -31.56405 -417.19577
问题三
在这个部分中,我们针对我们在问题 1,2 中的方法可能带来的误差及稳定
性提出了一种基于计算机模拟的检验方法。
在我们的前面的推导中,我们知道任意一个点在世界坐标系上的点的坐标和
其像在像坐标系的点的坐标满足(0)式,在模型中我们为了简化计算把 cZ 作为
了一个常量来计算(而事实上当拍摄距离较远,标靶上的圆又较小时,这种假设
是合理的),这带来了我们的第一个误差;我们的模型中的另外两个个误差来自
于边缘检测的精确度和用多元线性回归来计算像的曲线方程所带来的误差,而实
际上边缘检测的误差实际上来自于数字图像本身(由于一张图片只能存储有限的
点,所以将具体的像点映射到像素矩阵时会带来误差,但不会超过一个像素),
而从问题二的结果中我们可以得知多元线性回归求曲线方程有着很高的精度,所
以这两个误差相对于我们的第一个误差影响非常的小。
我们知道,(2)式在我们的针孔模型中是精确成立的,如果我们知道了相机
和标靶的确切的位置关系,那么对于标靶上的任意一点,我们总可以根据(2)
式计算出该点的像在像坐标系中的位置,也就是说在这个前提下,圆心的像在像
坐标系中的位置可以精确计算出来,这就是我们这个模型中的理论解。我们的目
的就是分析我们模型的解和理论解之间的关系。
在提出检验方法之前,我们先提出下面的引理。
引理 1:若世界坐标系上按照旋转矩阵 3 3R ´ 进行旋转,然后按向量
数学中国论文共享
www.madio.cn
16
( , , )T T TT X Y Z= 进行平移得到像坐标系,那么原来在世界坐标系中的一点
( , , )W W Wp X Y Z 在光心坐标系中的坐标的 Z轴分量 CZ 满足:
31 32 33C W W W TZ R X R Y R Z Z= + + +
证明:直接由(0)式化简可得。
在我们下面的方法中,我们的世界坐标系 WO 的原点在半径为 12mm的圆的圆
心上, W WX OY 平面与圆所在平面平行,光心坐标系 CO 的原点及相机的光心,
W WX OY 平面与像平面平行,将光心坐标系 CO 沿主光轴平移到主光轴与像平面的
交点即得像坐标系 PO ,如下图:
图十一
引理 2:圆心的像在像平面的坐标的理论值为 ( , ,0)T T
T T
LfX LfY
Z Z
- -
证明:由:
0 0
0 0 0 0 0 0 0 0
, ,0 01 10 0 0 0 0 0 0 0
0 00 ,1 0 ,1
0 0 1 0 0 1 0 0 0 1 01
1 1
o
o T T
C T
u L f Lf
R T R T
v L f Lf
Z Z
é ù é ù
é ù é ù é ù é ùê ú ê úé ù é ùê ú ê ú ê ú ê úê ú ê ú= =ê ú ê úê ú ê ú ê ú ê úê ú ê ú-ë û ë ûê ú ê ú ê ú ê úê ú ê úë û ë û ë ûë û
ë û ë û
得:
T
o
T
T
o
T
LfXu
Z
LfYv
Z
ì = -ïï
í
ï = -
ïî
数学中国论文共享
www.madio.cn
17
检验方法的提出:
Step 1 选取合理的旋转矩阵 R 和平移向量T,即初始化相机与标靶的相对
位置。
Step 2 在圆周上进行等间隔采样,对于每一个采样点,其在世界坐标系的
坐标可以表示为 (12cos ,12sin ,0)q q ,那么可以得到其像点在像坐标系的坐标为:
0 0
31 32 33
12cos
12sin1 1
0
1
11
W
W
C W W W W T
X
u
Y
v M M
Z Z R X R Y R Z Z
q
q
é ù é ù
é ù ê ú ê ú
ê ú ê ú ê ú= =ê ú ê ú ê ú+ + -
ê ú ê ú ê úë û
ë ûë û
Step 3 经过第二步后我们得到了圆的像的边缘点集,这些点是精确的,也
是我们的模型中可以从图片中准确获取的,所以我们将这个边缘点集作为我们模
型的输入,然后利用该边缘点集按照我们的模型求解出圆心的像在像坐标中的坐
标 ( , ,0)u v 。
Step 4 计算 ( , ,0)u v 与理论值 ( , ,0)o ou v 的误差 D(这里用欧氏距离计算)。
Step 5 重复执行 1到 4步,统计误差数据。
在用计算机进行模拟的时候,我们特别关心两个问题,一是相机与标靶的距
离对结果的影响;二是像平面与圆平面之间的偏角(由旋转矩阵 R决定)对结果
的影响。
接下来我们来研究 R和T的合理选取的问题。
由上文的分析可知,旋转矩阵 R中的三个参数 , ,a b l分别代表 y轴向 z轴旋
转的角度, z轴向 x轴旋转的角度, x轴向 y轴旋转的角度,例如下图:
图十二
我们几乎不用考虑 0l ¹ 的情况,因为这实际上是将相机左右倾斜,这并不
影响拍出的物体的形状和大小,而只会影响拍出的物体在整张照片中的位置。而
数学中国论文共享
www.madio.cn
18
0a ¹ 对应的情景时我们在物体的上(下)前方将镜头对着其进行拍摄, 0b ¹ 对
应的情景时我们在物体左(右)前方将镜头对着其进行拍摄, ,a b 体取值范围为
[ , ]
2 2
p p
- 。
对 ( , , )T T TT X Y Z= ,由于圆心在光心坐标系的坐标即 ( , , )T T TX Y Z- - - ,而圆必
须在相机的前方,故 0TZ- > ,即有 0TZ < 。除此之外,还需考虑T的长度,即
相机与标靶的距离,从我们前一问的数据中,可以估算出这段距离大约在 500
毫米左右。
基于以上分析,下面我们用计算机模拟的方式对我们的模型进行检验。
1.相机与标靶的距离对结果的影响
令 0a g= = , 0.7b = - , ( , 0, ),1 1000T k k k= - - £ £ ,我们对每一个 k对应的
相机和标靶的相对位置计算出我们模型的解与理论解的误差。
实验结果表明,当1 20k£ £ ,即相机与标靶的距离大约小于 30毫米时,存
在较大误差,模型很不稳定,误差图如下:
图十三
特别指出的是,当 k在某个数附近时出现了极不稳定的情况,误差也达到了
极大值。
同时,模拟结果也表明,当 140k ³ 时,即相机与模型的距离大约大于 200
毫米时,误差始终保持在 1个像素以内,而且随着 k的增加误差 D在不断减小。
数学中国论文共享
www.madio.cn
19
图十四
这个结论表明,在相机距离标靶在 200毫米(即 20厘米)以上时有着非常
高的精度以及稳定性,而我们实际拍摄时的物距很少小于 20厘米,而且本文中
的物距大约是 500毫米左右,由图可以知道当距离大于 500毫米时,误差仅仅在
0.1个像素之内,所以可以认为我们的模型是比较理想的。
2.像平面与圆平面之间的偏角对结果的