上 海 海 事 大 学 学 报
Vo l . 27 No . 3
第 27卷 第 3期
2006年 9 月 JOURNAL O F S HAN GHA IMAR IT I M E UN I V ER S I TY
Sep. 2006
文章编号 : 1672 29498 ( 2006) 0320090 205
基于子空间投影技术的椭圆拟合算法
任 蕾 , 杨忠根
(上海海事大学 信息工程学院 , 上海 200135)
摘 要 ：通过两次子空间投影 ，把拟合误差矢量正交分解为 3个分量 ，并使它们的范数同时最小化 ,
导致一个有约束总体最小化过程 , 该过程在同时约束其中的 2 个正交分量 (即误差均值和误差矢 量在中
心化一阶子空间的投影分量 )为 0 的条件下使另 1 个分量 (即误差矢量在中心化一阶子空 间的正交投影分量 )的范数最小化 . 通过结合子空间投影与奇异值分解 ，采用渐进计算技术来完成 椭圆参数矢量的 3个分量的最优估计. 理论分析和实验表明 ，该算法具有快速 、精确 、抗噪能力强 、 拟合成功率高等优点 .
关键词 : 计算机视觉 ; 椭圆拟合 ; 子空间投影 ; 奇异值分解 ; 正则化 中图分类号 : TP391. 41
文献标志码 : A
Elli pse f i tti n g a l gor i th m b a s ed on subsp ace projecti o n techn i que
R E N L e i , Y AN G Zhonggen
( I nf o r m a t ion Eng . Co l lag e, S han g ha i M a r iti m e U n i v . , S h ang ha i 200135 , Ch i na )
A b s tra c t :
B y m e an s of t w i ce sub s p a ce p r o j ec t i o n s , fitti ng e r r o r vec t o r is o rthogona l l y decompo s ed a s th ree componen ts . The op ti m a l e sti m a ti o n of the e lli p se p a ram e te r
vec t o r is gi ven by si m u ltaneou sl y m i n i 2 m izi ng the no r m s of th ree componen ts . It re su lts i n a t o ta l op ti m iz a ti o n p r oce ss i n wh i ch the no r m of one componen t, i . e . o rthogona l p r o j ec ti o n of fitti ng e rr o r vec t o r i n cen tra liz ed o rde r 2one sub sp ace, i s m i n i 2 m iz ed sub j ec t t o the con stra i n ts tha t the o the r t w o componen ts, i . e . the e rr o r ave rage and the p r o j e c t i o n of fitti ng e rr o r vec t o r i n cen tra liz ed o rde r 2one sub sp ace, a re equa l t o z e r o s . The op ti m a l p r oce s s can be com b i ned w ith si ngu l a r va l ue decompo siti o n and p r ogre ssi ve l y comp l e ted. The theo re ti ca l ana l ysis and ex 2 p e ri m en ta l demon stra ti o n p r ove tha t the new a l go rithm is fa st and exac t, its an ti 2no ise ab ility is str o ng and its ra t e of succe s sfu l l y fitti ng is h i gh.
Key word s : comp u t e r visi o n; e lli p se fitti ng; sub s p a ce p r o j ec t i o n; si ngu l a r va l ue decompo s iti o n; regu l a r 2 iz a t i o n
常见的封闭曲线基元 . 椭圆提取的核心技术是椭圆
拟合 . 长久以来 ，椭圆拟合一直受到学术界的高度重
视 ，许多对椭圆提取的研究 [ 1～9 ]
都涉及到它 . 传统技 术是标准奇异值分解 ( SVD ) 技术 ，即椭圆参数矢量
0 引 言
椭圆提取在图像识别与计算机视觉中有着特别
重要的意义 ，因为包括圆在内的椭圆是最简单的 、最
收稿日期 : 2005 210 208；修回日期 : 2006 207 203
基金项目 ：上海市高等学校科学技术发展基金 ( 01 G 02 )
作者简介 ：任 蕾 ( 1979 2) ，女 ，山东淄博人 ，助教 ，在读博士 ，研究方向为通信与信息系统 , ( E 2m ail ) le ir en@ c ie . shm tu . edu . cn
是与数据阵的最小奇异值对应的右奇异矢量 . 理论 分析和实验结果都表明 ，数据阵的条件数过大会导 致标准 S VD 算 法的 实 际性 能较 差 , 存在 估 计 偏 差 大 、均方拟合误差大和噪声灵敏度高的缺点. HA R T 2
L EY [ 10 ]的正则化技术是解决数字椭圆拟合在数 值
计算上病态性的有效手段 . 实践表明 ，其效果仍然受
限 .
现代信号处理认为 ，观测空间由信号子空间和 噪声子空间构成 ，把观测数据阵投影到噪声子空间 或信号子空间可开发出行之有效的噪声子空间法或 信号子空间法 . [ 11 ]
子空间投影技术已广泛应用于现
代谱估计和自适应滤波等领域 . 能否利用它来提高 椭圆拟合的性能并简化计算仍是个有待回答的 、有 益的研究课题 . 通过分析可以看出 ，椭圆拟合中的数 据空间由零阶子空间 、一阶子空间和二阶子空间构 成 ，正是这 3个能量差异相当大的子空间之间的紧 耦合导致传统椭圆拟合算法的数值特性不佳 . 利用 两次子空间投影技术把拟合误差矢量正交分解为 3 个分量 ，并使它们的范数同时最小化 ，导致在同时约 束其中的两个正交分量 (即误差均值和误差矢量在 中心化一阶子空间的投影分量 )为 0 的条件下使另 一个分量 (即误差矢量在中心化一阶子空间的正交 投影分量 ) 的范数最小化 . 这是一个有约束总体最 小化过程 ，并且可结合奇异值分解 ，采用渐进技术完 成椭圆参数矢量的 3 个分量的最优估计. 这就是本 文提出的椭圆拟合新算法 .
的右奇异矢量 θ^, 此时 ,σ2
m in 给出拟合误差 ; 步 3. 检验 所拟 合 的二 次曲 线 是否 为椭 圆 , 若 是 , 则椭圆拟合成功 , 并计算椭
圆的特征参数. 2 结合 HAR T L E Y 正则化变换的 S VD 椭圆拟合 理论与实践表明 ，数字椭圆参数矢量的拟合误 差正比于 数 据 阵 D 的 条 件 数 K ( D ) = ‖D +
‖ ·
σ1 (D )
‖D ‖ = (D ) , 本文规定矩阵的奇异值按降序 排
σ5 列 .
事实上 , 可把数据阵 D 分解为 3 个有不同基底
的子空间 D = [ D ( 1 ) D ( 2 ) D ( 3 ) ], 其中 , D ( 1 ) = 1N 为 N 个 1 (坐标的零次项 ) 组成的全 1 矢量 , 被称为
x 1
x 2 y 1
y 2
零阶子空间 ; D ( 2 ) =
由坐标的一次项构成 ,
x y N N
2
2
x 1
x 1 y 1
x 2 y 2
y 1 2 2 x 2 y 2 被称为一阶子空间 ; 而 D =
由坐
( 3 ) 2
正则化坐标x 2
y x y N
N N
N
标的二次项构成 , 被称为二阶子空间. 这时 , 椭圆拟
合方程可改为
D ( 1 ) θ1 + D ( 2 ) θ2 + D ( 3 ) θ3 ( 4 )
a 1
= 0
T
其中 ,θ1 = [ a 5 ], θ2 = [ a 3
a 4 ] 和 θ3 = [ a 0
1 基于 S V D 技术的椭圆拟合
椭圆拟合是用方程 :
T
a 2 ] .
正是这 3个能量相差极大的子空间的紧耦合使
得数据阵 D 的条件数极大 [ 10 ] , 所以 SS VD 算法的误
差往往较 大 , 须 增 加 正 则 化 技 术 才 能 实 用 . HA R T 2 L EY 的坐 标 正 则 化 技 术 是 简 单 实 用 的 正 则 化 技 术 [ 10 ] ：首先采用预白化变换坐标 ，使得变换后的坐 标中心位于坐标原点 ，散布半径为 1; 然后在变域中 提取椭圆 ；最后把它映射回原域中 ，得到待估的椭圆 参数 . 这给出了算法 2: 结合 HA R TL EY 正则化变换 的 SVD 椭圆 拟合 算 法 ———H S VD 算法 . H S VD 的 条 件数远远低于 S VD 的条件数 ，所以 H SVD 算法的性 能显著地优于 SSVD 算法.
2
2
a 0 x + a 1 xy + a 2 y + a 3 x + a 4 y + a 5 = 0 ( 1 )
, N } , 其
T
拟合二维点集 { p i = [ x i y i 1 ] | i = 1, 2, 等价的线性矩阵 - 矢量方程为
T
D i θ = 0
其中 , 数据矢量 D i = [ 1 x i y i ( 2 )
2
2 T
y i ] , 参
x i x i y i T
数矢量 θ = [ a 5 a 3 a 4 a 0 a 1 a 2 ]
, 其最优解 可用标准 S VD 技术给出 , 即 θ的最优估计 θ^为数据
T
D N ] 的最小奇异值 σm in 所 阵 D = [ D 1 D 2 对应的右奇异矢量 νD , m in ,
即 T
i f D = U D ΣD V D , then θ^ =νD , m in 这样 , 有如下所述的算法 :
( 3 )
3 基于子空间投影的椭圆拟合
拟合误差矢量为
算法 1. 基 于 标 准 SVD 技 术 的 椭 圆 拟 合 算 法 ———SSVD 算法 : 步 1. 把二维点集转换成数据矢量集合 , 并构成 数据阵 D ;
步 2. 计算数据阵 D 的最小奇异值 σm in 所对应
= D ( 1 ) θ1 + D ( 2 ) θ2 + D ( 3 ) θ3
( 5 ) e 1 T ⊥
使用投影算子 = 1N 1 和正交投影算子 P D ( 1) P N D ( 1)
N
上海海事大学学报第27卷92
]
1
1T
θ
3
的最优估计θ^3 为数据阵D (3 ) = U D , 2 D o(3 ) 的最小= I - 1T 可把拟合误差矢量正交分解成
N N 2 N
]
奇异值σm in (D (3 ) )所对应的右奇异矢量时, 拟合误
差矢量在中心化一阶子空间的正交投影分量的长度⊥
e = e
1
+ e
1
, 其中,
T T
e
1
= P
D ( 1)
e = 1
N
(θ1 + D¯(2 )θ2 + D¯(3 )θ3 ) ,
(6)]
⊥
‖e1 ⊥2 ‖达到最小值σm in (D ( 3 ) ).
总之，有如下所述的算法:
算法3. 基于子空间投影SVD 技术的椭圆拟合
算法———SSPSVD 算法:
步1. 把二维点集转换成数据矢量集合, 并构成e⊥⊥
= P
D (1)
e = D o(2 )θ2 + D o(3 )θ3
1
这儿, D¯T = [ x¯y¯], D¯T = [ x22y ],x¯=
xy
( 2 )( 3 )
1
6 x , y¯= 1 6 y , x2 = 1 6 x2 , xy = 1 6
N N N N
x y , y2 =
i i i i i
N N i = 1 N i = 1 N i = 1
i = 1
1 6N2y , D o
= P⊥- 1 D¯T 和D o=
= D
D
i ( 2 ) D ( 1)( 2 ) ( 2 ) N ( 2 )( 3 )子数据阵D 和D . 通过中心化运算，把它们变为N i = 1 ( 2 ) ( 3 )
P⊥T
- 1
N
D¯( 3 ) 分别为中心化一阶子空间D o( 2 ) 和D o(3 ) .
= D
( 3 )
D ( 1)
D (
3 )
2. 对D o( 2 ) 进行瘦型S VD 分解: D o( 2 ) = 和中心化二阶子空间.
同理, 使用投影算子
步
T ]
= D o( 2 ) ( D o( 2 )
P
D o
V T T
D 2, 1
Σ
D 2 D 2
, 计算D( 3 ) = U D
2, 2
D o(3 ) , 取θ^3 为它的最小
U
( 2)
- 1 T ⊥T
D o( 2 ) ) D o(2 ) 和正交投影算子= I - P 可把]
奇异值σ
m in
(D (3 ) )所对应的右奇异矢量.
P
D o D o(2)
( 2)
e⊥正交分解成- 1 T
步3. 计算θ^2 = - V D(ΣD)U D , 1 D o(3 )θ^3 和θ^1 1
2 2 2
e⊥⊥
= e
1 ⊥2
+ e
1 ⊥2
,
其中,= - (D¯T + D¯θ^ ) .
( 3 ) 3
T
( 2 )
θ^
2
1
⊥
= P
D o(2)
e1 = D o( 2 ) θ2 + P D o
(2)
D o( 3 ) θ3
e1 ⊥2步4
. 检验所拟合的二次曲线是否为椭圆, 若
是，则椭圆拟合成功，并计算椭圆的特征参数.
上述分析表明, 基于子空间投影理论的椭圆拟
合等价于一个在约束拟合误差均值和拟合误差矢量
在中心化一阶子空间的投影分量为0的条件下使拟
合误差矢量在中心化一阶子空间的正交投影分量的
范数最小化. 这是一个渐进完成的有约束总体最优
化过程, 而且参数矢量的渐进式求解可结合SVD 技
术进行，以达到最佳的拟合效果. 该过程利用子空间
投影技术把 3 个能量相差相当大的子空间完全解
(7)
e⊥⊥⊥⊥
= P
D o(2)
D o(3 )θ3
= P
D o(2)
e
1
1 ⊥2
通过两次子空间投影把拟合误差矢量正交分解
成 3 个分量, 使拟合误差矢量 e 的范数(即矢量长
度) ‖e ‖最小化等价于使其 3 个正交分量e1 , e1 ⊥2
和e⊥的范数‖e ‖, ‖e‖和‖e⊥‖同时最小
1 ⊥
2 1 ⊥2 1 ⊥2
1
化.这就把一个较复杂的最优化问题分解为
对简单的、互相独立的子最优化问题.
显然，当且仅当
3 个相
T T
θ
1
= - (D¯( 2)θ2 + D¯(3 ) θ3 )(8)
]
时, e1 = P D
(1)
e = 1N e¯= 0, 使得‖e1 ‖达到最小值0,⊥
耦, 使得最后的数据阵D ( 1) = P D D o(1 ) 仅属于与零
( 2)
阶和一阶子空间都正交的中心化二阶子空间，达到
大大降低拟合条件数的目的. 如果把它与HAR T2
L EY的预正则化技术相结合(把点坐标集合变换成
具有零均值和么方差阵，在此变域中用SSPSVD 算
法计算变域椭圆参数，把它映射为原域椭圆参数. 这
就是算法 4 ———结合HA R T L EY正则化变换的子空
间投影SV D 椭圆拟合算法H S SPSVD ) ，就可进一步
降低条件数，从而进一步改善算法性能.
其中, 帄均拟合误差e¯=
1 6N
e . e¯为0 意味着拟合
i
N i = 1
误差具有零均值, 这是被线性最优滤波理论预期的
结果. 在拟合误差具有零均值的条件下, e⊥= P⊥e
1 D ( 1)
⊥⊥
= e o= e，并且, e1 ⊥2 = P D o e和e
1 ⊥2
= P
D o
e分别为拟
( 2)( 2)
合误差矢量在中心化一阶子空间的投影分量和正交
投影分量.
同样地, 容易证明, 当且仅当
- 1 T
θ
2U D , 1 D
o(
3 )
θ
3
= - V
D
(Σ)(9)
2 D 2 2
4 实验
时, 拟合误差矢量在中心化一阶子空间的投影分量
e1 ⊥2为0, 使得‖e1 ⊥2 ‖也达到最小值0. 式中的矩阵为了分析所推荐技术的工作性能，以椭心( x ,
c U D
2, 1
,ΣD
2
和V D
2
由D o(2 ) 的瘦型SVD 分解式D o(2 ) =
= 15, 10 ，长短半轴a, b = 25, 12 )和倾角α
y
c
)()()(
T ⊥
U D
2, 1
Σ
D 2
V D
2
给出.约束e1 = 0 和= 0 时 e = e1 ⊥2 ,
e1 ⊥2
= 30 °的椭圆为例, 比较 4 种椭圆拟合算法( SSV D
算法、H S VD 算法、SSPSVD 算法和H S SPSVD 算法)
的性能.
图1和图2分别示出4个算法对一个噪声强度使得拟合误差矢量长度为‖e‖= ‖e⊥‖.
1 ⊥2
e⊥P⊥T U
考虑到= e = U e =
1 ⊥2D o(2) D 2,
2 D 2, 2
T ⊥T
U
D 2, 2
U
D 2, 2
D o(3 )θ3使得‖e1 ⊥2 ‖= ‖U D
2, 2
D o(3 )θ3 ‖，其
中，矩阵U D , 2是矩阵U D , 1 的正交补矩阵.因此, 当 2 2
σ= 2. 00的加噪的全椭圆和对一个噪声强度σ=
2 2
0. 05的加噪的 1 /4 椭圆的拟合情况 . 其中 , 实线为 理想椭圆 ，十字点为 加 噪的 椭圆 样本 点 , 细 虚线 为
SSVD 拟合的椭圆 ，点划线为 H S VD 拟合的椭圆 ，长
虚线为 SSPSVD 拟合的椭圆 ，加粗的虚线为 H S SPS 2 VD 拟合的椭圆 .
‖
θ‖) )的数学期望 , 相应的性能比较见表 1. 图 3 4种算法对加噪 1 /4椭圆的拟合误差
图 1 4种算法对加噪全椭圆的拟合
图 4 4种算法对加噪全椭圆的拟合误差
注意 ，这里所列的加噪情况实际上属于中等噪
声或强噪声情况 .
实验中取 500 次实验结果计算帄均参数矢量误 差和拟合成功率 .
图 2 4种算法对加噪 1 /4椭圆的拟合
图 3和图 4 示出对不同噪声强度而言的拟合全 椭圆时和拟合 1 /4椭圆时 4 种算法的参数矢量帄均
误 差 , 即 参 数 矢 量 误 差 角 a r cco s (θ^T
θ / ( ‖θ^ ‖ ·
表 1 拟合全椭圆时和拟合 1 /4椭圆时 4种算法的性能比较
拟合全椭圆时的噪声强度 σ2
拟合 1 / 4 椭圆时的噪声强度 σ2
性能 指标 使用算法
0. 00 0. 05 0. 50 1. 00 2. 00 4. 00 0. 00 0. 05 0. 50 1. 00 2. 00 4. 00 SS VD H S V D SSPS V D H SSPS VD 0. 000
0. 000 0. 000
0. 000 0. 100 0. 022 0. 020
0. 012 0. 128 0. 035 0. 031
0. 020 0. 140 0. 052 0. 048
0. 030 0. 154 0. 076 0. 070
0. 049 0. 182
0. 097 0. 092
0. 063 0. 000 0. 000 0. 000
0. 000 0. 196 0. 106 0. 105
0. 064 0. 202 0. 118 0. 115
0. 072 0. 209 0. 123 0. 120
0. 081 0. 217 0. 134 0. 132
0. 091 0. 231 0. 142 0. 139
0. 100 参数 矢量 帄均 误差
SS VD
100 100 100 99. 7 98. 6 95. 7
100 94. 3 91. 6 88. 1 85. 5 80. 4 H S V D 100 100 100 100 100 100 100 100 100 100 99. 8 99. 3 拟合
成功率
SSPS V D 100 100 100 100 100 100 100 100 100 100 99. 9 99. 6 H SSPS VD
0. 00
0. 05
0. 50
1. 00
2. 00
4. 00
100
100
100
100
100
100
实验表明 , SSVD 在 所有 加 噪情 况下 都 给出 最 差的性能 ，不但参数矢量误差较大 ，而且有拟合不出 椭圆的情 况 , H S VD 的 性 能 较 之 有 了 明 显 的 改 进 ,
SSPSVD 在 性 能 上 与 H S V D 相 同 或 略 优 于 H S VD , H SSPSVD 在所有情况下给出最好的拟合性能.
上 海 海 事 大 学 学 报 第 27卷
94
实验表明 , H S VD 的性能显著优于 SSVD , SSPS 2 VD 的性能与 H SVD 相近 ，有时还略优 ; H SSPSVD 的
性能 较 H S VD 有 可 观 测 到 的 改 进 , SSPSVD 常 与
H SSPSVD 有相近的性能 . 这说明 , 子空间投影技术
有强自正则化作用 ；而且 ，它通过把一个较复杂的最 优化问题分解为 3 个相对简单的 、互相独立的子最 优化问题 ，渐进地计算椭圆参数 ，达到简化算法 、加
速计算和改善性能的目的 .
事实上 ，子空间投影技术可推广应用于其他领 域 ，如计算机视觉中的摄像机自标定和三维重建 ，我 们正在开发和实现有关算法 .
5 结论与展望
通过两次子空间投影 ，把拟合误差矢量正交分 解为 3个分量 ，并使它们的范数同时最小化 ，导致一 个有约束总体最小化过程 , 该过程在同时约束其中 的两个正交分量 (即误差均值和误差矢量在中心化 一阶子空间的投影分量 )为 0 的条件下使另一个分 量 (即误差矢量在中心化一阶子空间的正交投影分 量 )的
范数最小化 . 通过结合子空间投影与奇异值 分解结合 ，采用渐进计算技术来完成椭圆参数矢量 的 3 个分量的最优估计 .
参考文献 :
[ 1 ] 杨忠根 , 栾晓明. 利用椭圆性质提取椭圆 [ J ]. 中国图像图形学报 , 1997 , 2 ( 7 ) : 517 - 519.
[ 2 ] 杨忠根 , 栾晓明 , 赵昶冰. 直线和二次曲线的实时鲁棒精确提取 [ J ]. 哈尔滨工程大学学报 , 1997 , 18 ( 4 ) : 69 - 77.
[ 3 ] Y A M G Z G , J I A N G G X, R EN L. E l lip se de tec t ion ba sed on i mp r oved 2GEV D techn ique [ C ] / /W C I CA ’2004 F if th W o rld Congr ess on I n telligence Con tr o l and A u to m a tion Conf e r ence P r oceed ings, 2004: 4 181 - 4 185. [ 4 ] 杨忠根 ，马彦. 使用广义正交概念的 K 2RA N S AC 椭圆提取 [ J ]. 自动化学报 , 2002 , 28 ( 4 ) : 520 - 526.
[ 5 ] CH EN G Y C, L EE S C . A ne w m e thod f o r quad ra tic cu rve detec tion u sing K 2RAN S AC w ith acce le r ation techn ique s [ J ]. Pa tte r n R ec ogn ition,
1995 , 28 ( 5 ) : 663 - 682.
[ 6 ] OL S ON C F . Con s tr a i ned Hough t r an sf o r m f o r cu rve de tec tion [ J ]. Comp u t e r V ision & I m age U nde rstand ing, 1999 , 73 ( 3 ) : 329 - 345. [ 7 ] CH U TA TA P E O , G UO L. A mod if ied Hough tr an s f o r m f o r line detec tion and its p e r f o r m ance [ J ]. Pa tte r n R ec ogn ition, 1999 , 32 ( 2 ) : 181 - 192. [ 8 ] B ENN ETT N , BURR I D G E R , S A I TO N. A m e t hod to de tect and charac te riz e e llip se s u sing the Hough tr an sf o r m [ J ]. I EEE Tr an s Patte r n A na ly 2
sis and M a ch ine I n te lligence, 1999 , 21 ( 7 ) : 652 - 657.
[ 9 ] F I TZG I B BAN A , P I L U M , F I SH E R R B. D ir ec t lea st squa r e f itting of e l lip se s[ J ]. I EEE Tr an s Patte r n A na lysis and M ach ine I n t e lligence, 1999 ,
21 ( 5 ) : 476 - 480.
[ 10 ] H AR TL EY R I . I n def en s e of e i gh t 2po i n t a lgo rithm [ J ]. I EEE Tr an s Pa t te r n A n alysis and M ach ine I n te lligence, 1997 , 19 ( 6 ) : 580 - 593. [ 11 ] 张贤达. 现代信号处理 [M ]. 北京 : 科学出版社 , 1997.
(编辑 李佩芬 )
版权声明：本站内容均来自互联网，仅供演示用，请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198，我们将在24小时内删除。
发表评论