生如蚁,美如神

众里寻她千百度,蓦然回首,那人却在灯火阑珊处

基于几何距离的椭圆拟合

问题

给定离散点集Xi=(xi,yi),我们希望找到最好的椭圆去拟合这些离散点。

方法

通常我们使用最小二乘法求解如下的最优化问题:

Mini=1Nf(xi,E)2

这里f(xi,E)表示点xi到E(指待拟合的椭圆)的最小距离。

通常我们有两种方法来表达f(xi,E),分别是:几何距离代数距离。前面我们已经描述了基于代数距离的椭圆拟合,下面我们将重点介绍基于几何距离的椭圆拟合方法,这种方法也是椭圆拟合方法中最鲁棒、精度最高的拟合方法。

我们主要参考论文《Least-squares orthogonal distances “tting of circle,
sphere, ellipse, hyperbola, and parabola》

论文的核心思路

显然,由于离散点到椭圆的几何距离是非线性的,因此椭圆拟合是一种非线性优化的问题,需要通过多次迭代求解。

因此下面我们阐述论文的思路时,将分成两个片段来讲解。第一段主要介绍基于高斯-牛顿法的非线性优化方法第二段,具体介绍这种方法在椭圆拟合中的应用

非线性最小二乘拟合

问题

考虑一般的非线性最小二乘问题如下:

mina  ||XF(a)||2                         1

这里aRq为优化参数,F为非线性函数,XRp是已知向量。如下我们给出基于梯度的迭代公式:
ak+1=ak+λAJTF(XF(ak))                         2

这里λ为步长,A为缩放因子,J_F为F在当前参数a下的Jacobian矩阵。

各种优化方法不同,主要是A的选择不同。具体如下:

  • A=H1时,为牛顿迭代法。
  • A=(JTFJF)1时,为高斯-牛顿迭代法。
  • A=I时,为梯度下降法。
    这里我们选择使用高斯-牛顿迭代法。

我们将A带入(2),即可以改写为:

{ak+1=ak+λΔa                 (3)JFΔa=XF(ak)              (4)

这里JF=(Fiaj),i=1,2,...,p,j=1,2,...,q
并且,我们根据(1)可以写出迭代算法的性能表现指数,其实就是参差的表达式:
σ20=||XF(a)||2=[XF(a)]T[XF(a)]              (5)

求解

接下来,我们需要求解出Δa,才能进行迭代。

事实上求解Δa只需要求解方程组(4)即可.我们可以使用奇异值分解求解方程组。

基于几何距离的椭圆拟合

平面上的椭圆可以使用5个参数唯一地表达,即圆心(Xc,Yc)、主轴a,b,角度α(π/2<απ/2).

对于椭圆的最小二乘正交距离拟合,我们将使用第一段所讲的方法。此时,我们定义

a^=(Xc,Yc,a,b,α)T

X可以看成给定的离散点Xi,而自然F(a^)就是定点Xi在椭圆上的最近点Xi(也称为正交关联点)。迭代公式就变成了如下:
{a^k+1JXi,a^Δa^==a^k+λΔa^                      (6)XiXii=1,2....m    (7)

这里m指给定的离散点的个数。各个离散点均满足(7),因此关于(7)可以联立。实际上关于(7)的方程是2m个,而未知数a的维数是5,显然2m>>5,因此解是不唯一的,我们需要求出最小二乘解即可。

显然我们必须计算Xi以及在Xi点上的Jacobian矩阵JXi,a^。下面我们将阐述如何求解这两个量。

坐标系转换

为了便于后面的求解,我们需要考虑,将原来的坐标系XY,旋转α角,建立一个位于(0,0)的临时坐标系xy。见Fig.3


椭圆坐标轴

则有

x=R(XXc)              (8)

or
X=R1x+Xc              (9)

这里
R=(CSSC)R1=RT

C=cos(α),S=sin(α)

椭圆上的正交关联点

经过坐标轴转换,5个椭圆参数变成了2个,仅仅包含了主轴a,b。椭圆上的点可以用标准方程描述如下:

x2a2+y2b2=1              (10)

从Fig.3上,我们可以看到位于xy坐标系上的正交关联点xi,yi的切向量与两点(即xi,yixi,yi)的连线是正交的。因此,有:
dydxyiyxix=b2xa2yyiyxix=1              (11)

重写(10,11)有:

f1(x,y)=12(a2y2+b2x2a2b2)=0              (12)f2(x,y)=b2x(yiy)a2y(xix)=0           (13)

上面两式称为正交关联条件,我们将使用广义牛顿法来求解上面方程组的解。
即使用如下的迭代公式求解:
{xk+1=xk+ΔxQkΔx=f(xk)                      (14)

Q=f1xf2xf1yf2y=(b2x(a2b2)y+b2yia2y(a2b2)xa2xi)               (15)

对于迭代公式(14),我们提供初始解x0,如下:(见Fig.3)
x0=12(x^k1+x^k2)

这里
x^k1=(xk1yk1)=(xiyi)ab/b2x2i+a2yi


ak

一般经过3-4轮迭代就可以提供了一个很高精度的解。

我们先把Xi通过转换公式(8)转换成xy坐标系的xi,然后通过求解广义的牛顿法求解得到正交关联点xi,最后再通过转换公式(9)转换成XY坐标系的Xi,最后给出正交误差距离向量为:

X′′i=XiXi                      (16)

椭圆上的正交关联点的Jacobian矩阵

下面我们直接给出椭圆上的正交关联点的Jacobian矩阵,(注:本部分推导比较复杂,贴出本人的详细推导过程,需要的可以参考推导1,2,3,4,5

JXi,a^=(R1Q1B)|x=xi                      (17)

这里Q见(15),B=(B1,B2,B3,B4,B5)
此时



椭圆的正交距离拟合

给定m个二维点,我们可以利用(16)、(17)给定的误差距离向量X′′i和Jacobian矩阵JXi,a^,构造p(=2m)个线性方程组。如下:



,当我们使用奇异值分解求解出上述方程组的解时,就可以进行高斯-牛顿迭代。
此时我们还需要一个初值,一般可以选择基于代数拟合的椭圆或者使用圆作为初始值来进行迭代。通常我们推荐使用圆来作为初始值参数。

显然,从圆拟合中获得的初始参数为:
(Xc,Yc)ellipse=(Xc,Yc)circle,a=b=R,α=0
并且在迭代过程中,如果a<b,则需要令ααsign(α)π2(保持a为主轴长)

标准坐标轴下的椭圆拟合

有时候,我们也需要为椭圆拟合增加一些约束,经典的就是要求在标准坐标轴下进行椭圆拟合(α=0π/2),或者要求增加椭圆面积约束.此时,我们只需要在原来的第(7)式增加如下约束即可。



此时w3,w4为权重参数,一般可以取1.0×106.

实验结果

例1.



我们对Table 7的离散点进行椭圆拟合,其中初始参数a0源自基于几何距离的圆拟合。
λ=1.2,在经过19次迭代后,最终的校正向量的范数为||Δa||=4.2×106,并且得到误差指数σ0=1.1719.见如下:



为了更好地比较,我们也考虑了初始参数a0源自基于代数距离的椭圆拟合,而达到与上述相同的估计,需要进行21次迭代(||Δa||=1.1×106).
我们也比较了我们的结果与Gander的结果,而后者达到相同的估计,需要进行71次迭代(||Δa||=1.0×106).具体可参考Table 8 的III,IV.

Gander的算法有一个缺点就是不能使用圆的结果作为初始参数,为了克服这个困难,一般取a=R,b=R/2.为了更加直接验证我们算法的鲁棒性,我们也考虑了使用两种设置作为初始值,分别见Table 8的II,III。即:
II:

a=R,b=R/2,α=0

III:
a=R,b=R/2,α=1.211
.
我们分别使用17次和23次迭代达到了相同的估计,其中校正向量的范数分别是1.2×106,5.2×106.

最后我们给出拟合图:



从以上结果来看,我们的算法是鲁棒的,即使初始值不够好,最终也能迭代到较好的效果。
例2.



我们对Table 9的离散点进行椭圆拟合,其中初始参数a0源自基于几何距离的圆拟合。
λ=1.2.

我们本例的目的是验证本算法在标准坐标系以及增加面积约束条件下的拟合情况。
为了做了两组实验:
I:

α=π2,w3=1.0×106

II:
α=0,ab=20,w3=w4=1.0×106

我们将本算法与Spa¨thLeast-squares fitting of ellipses and hyperbolas做了比较,具体结果可参考:




注意:Table11(b)提供了各参数之间的相关系数,发现a与b是强的负相关,显然这是由于面积约束造成的。

从最终的结果上看,我们的算法在标准坐标系以及增加面积约束条件下,依然是鲁棒的。

阅读更多
版权声明:本文为博主原创文章,转载请注明原地址。 https://blog.csdn.net/xiamentingtao/article/details/54934467
文章标签: 椭圆拟合
个人分类: 图像处理
想对作者说点什么? 我来说一句

matlab 拟合椭圆

椭圆拟合

Phoebexiang Phoebexiang

2017-08-02 16:25:10

阅读数:2254

没有更多推荐了,返回首页

加入CSDN,享受更精准的内容推荐,与500万程序员共同成长!
关闭
关闭