問(wèn)題

給定離散點(diǎn)集$X_i=(x_i,y_i)$，我們希望找到最好的橢圓去擬合這些離散點(diǎn)。

方法

通常我們使用最小二乘法求解如下的最優(yōu)化問(wèn)題：

$$Min\sum _{i=1}^{N}f(x_i,E{)}^2$$
這里$f(x_i,E)$表示點(diǎn)$x_i$到E（指待擬合的橢圓）的最小距離。

通常我們有兩種方法來(lái)表達(dá)$f(x_i,E)$，分別是：幾何距離和代數(shù)距離。前面我們已經(jīng)描述了基于代數(shù)距離的橢圓擬合，下面我們將重點(diǎn)介紹基于幾何距離的橢圓擬合方法，這種方法也是橢圓擬合方法中最魯棒、精度最高的擬合方法。

我們主要參考論文《Least-squares orthogonal distances fitting of circle,
sphere, ellipse, hyperbola, and parabola》。

論文的核心思路

顯然，由于離散點(diǎn)到橢圓的幾何距離是非線性的，因此橢圓擬合是一種非線性優(yōu)化的問(wèn)題，需要通過(guò)多次迭代求解。

因此下面我們闡述論文的思路時(shí)，將分成兩個(gè)片段來(lái)講解。第一段，主要介紹基于高斯-牛頓法的非線性優(yōu)化方法。第二段，具體介紹這種方法在橢圓擬合中的應(yīng)用。

非線性最小二乘擬合

問(wèn)題

考慮一般的非線性最小二乘問(wèn)題如下：
$$\underset{a}{\min }||X?F(a)|{|}^2（1）$$
這里$a\in R^q$為優(yōu)化參數(shù)，F為非線性函數(shù)，$X\in R^p$是已知向量。如下我們給出基于梯度的迭代公式：
$$a_{k+1}=a_k+\lambda A{J}_F^T(X?F(a_k))（2）$$
這里$\lambda$為步長(zhǎng)，A為縮放因子，J_F為F在當(dāng)前參數(shù)a下的Jacobian矩陣。

各種優(yōu)化方法不同，主要是A的選擇不同。具體如下：

• 當(dāng)$A=H^{?1}$時(shí)，為牛頓迭代法。
• 當(dāng)$A=(J_F^T{J}_F{)}^{?1}$時(shí)，為高斯-牛頓迭代法。
• 當(dāng)$A=I$時(shí)，為梯度下降法。
  這里我們選擇使用高斯-牛頓迭代法。

我們將A帶入（2）,即可以改寫為：
$$\{\begin{array}{c}{a}_{k+1}=a_k+\lambda \Delta a(3)\\ J_F\Delta a=X?F(a_k)(4)\end{array}$$
這里$J_F=(\frac{?F_i}{?a_j}),i=1,2,...,p,j=1,2,...,q$
并且，我們根據(jù)（1）可以寫出迭代算法的性能表現(xiàn)指數(shù)，其實(shí)就是參差的表達(dá)式：
$${\sigma }_0^2=||X?F(a)|{|}^2=[X?F(a){]}^T[X?F(a)](5)$$

求解

接下來(lái)，我們需要求解出$\Delta a$,才能進(jìn)行迭代。

事實(shí)上求解$\Delta a$只需要求解方程組（4)即可.我們可以使用奇異值分解求解方程組。

基于幾何距離的橢圓擬合

平面上的橢圓可以使用5個(gè)參數(shù)唯一地表達(dá)，即圓心$(X_c,Y_c)$、主軸$a,b$,角度$\alpha (?\pi /2<\alpha \le \pi /2)$.

對(duì)于橢圓的最小二乘正交距離擬合，我們將使用第一段所講的方法。此時(shí)，我們定義
$$\hat{a}=(X_c,Y_c,a,b,\alpha {)}^T$$
$X$可以看成給定的離散點(diǎn)$X_i$，而自然$F(\hat{a})$就是定點(diǎn)$X_i$在橢圓上的最近點(diǎn)$X_i^{\prime }$(也稱為正交關(guān)聯(lián)點(diǎn))。迭代公式就變成了如下：
$$\{\begin{array}{ccc}{\hat{a}}_{k+1}& =& {\hat{a}}_k+\lambda \Delta \hat{a}(6)\\ J_{X_i^{\prime },\hat{a}}\Delta \hat{a}& =& X_i?X_i^{\prime }，i=1,2....m(7)\end{array}$$
這里m指給定的離散點(diǎn)的個(gè)數(shù)。各個(gè)離散點(diǎn)均滿足（7）,因此關(guān)于（7）可以聯(lián)立。實(shí)際上關(guān)于（7）的方程是2m個(gè)，而未知數(shù)a的維數(shù)是5，顯然2m>>5,因此解是不唯一的，我們需要求出最小二乘解即可。

顯然我們必須計(jì)算$X_i^{\prime }$以及在$X_i^{\prime }$點(diǎn)上的Jacobian矩陣$J_{X_i^{\prime },\hat{a}}$。下面我們將闡述如何求解這兩個(gè)量。

坐標(biāo)系轉(zhuǎn)換

為了便于后面的求解，我們需要考慮,將原來(lái)的坐標(biāo)系XY，旋轉(zhuǎn)$\alpha$角，建立一個(gè)位于$(0,0)$的臨時(shí)坐標(biāo)系xy。見(jiàn)Fig.3

則有
$$x=R(X?X_c)(8)$$
or
$$X=R^{?1}x+X_c(9)$$
這里
$$R=\left(\begin{array}{cc}C& S\\ ?S& C\end{array}\right)和R^{?1}=R^T$$
$C=\cos (\alpha ),S=\sin (\alpha )$

橢圓上的正交關(guān)聯(lián)點(diǎn)

經(jīng)過(guò)坐標(biāo)軸轉(zhuǎn)換，5個(gè)橢圓參數(shù)變成了2個(gè)，僅僅包含了主軸$a,b$。橢圓上的點(diǎn)可以用標(biāo)準(zhǔn)方程描述如下：
$$\frac{x^2}{a^2}+\frac{y^2}{b^2}=1(10)$$
從Fig.3上，我們可以看到位于xy坐標(biāo)系上的正交關(guān)聯(lián)點(diǎn)$（x_i^{\prime },y_i^{\prime }）$的切向量與兩點(diǎn)（即$（x_i,y_i）、（x_i^{\prime },y_i^{\prime }）$）的連線是正交的。因此，有：
$$\frac{dy}{dx}?\frac{y_i?y}{x_i?x}=\frac{?b^{2}x}{a^{2}y}?\frac{y_i?y}{x_i?x}=?1(11)$$
重寫（10,11）有：

$$\begin{gathered} f_1(x,y)=\frac{1}{2}(a^2{y}^2+b^2{x}^2?a^2{b}^2)=0(12) \\ {f}_2(x,y)=b^{2}x(y_i?y)?a^{2}y(x_i?x)=0(13) \end{gathered}$$
上面兩式稱為正交關(guān)聯(lián)條件,我們將使用廣義牛頓法來(lái)求解上面方程組的解。
即使用如下的迭代公式求解：
$$\{\begin{array}{c}{x}_{k+1}=x_k+\Delta x\\ Q_k\Delta x=?f(x_k)\end{array}(14)$$
$$Q=\left(\begin{array}{cc}\frac{?f_1}{?x}& \frac{?f_1}{?y}\\ \frac{?f_2}{?x}& \frac{?f_2}{?y}\end{array}\right)=\left(\begin{array}{cc}{b}^{2}x& a^{2}y\\ (a^2?b^2)y+b^2{y}_i& (a^2?b^2)x?a^2{x}_i\end{array}\right)(15)$$
對(duì)于迭代公式（14）,我們提供初始解$x_0$,如下：(見(jiàn)Fig.3)
$$x_0=\frac{1}{2}({\hat{x}}_{k1}+{\hat{x}}_{k2})$$
這里
$${\hat{x}}_{k1}=\left(\begin{array}{c}{x}_{k1}\\ y_{k1}\end{array}\right)=\left(\begin{array}{c}{x}_i\\ y_i\end{array}\right)ab/\sqrt{b^2{x}_i^2+a^2{y}_i}$$

一般經(jīng)過(guò)3-4輪迭代就可以提供了一個(gè)很高精度的解。

我們先把$X_i$通過(guò)轉(zhuǎn)換公式(8)轉(zhuǎn)換成xy坐標(biāo)系的$x_i$,然后通過(guò)求解廣義的牛頓法求解得到正交關(guān)聯(lián)點(diǎn)$x_i^{\prime }$,最后再通過(guò)轉(zhuǎn)換公式（9）轉(zhuǎn)換成XY坐標(biāo)系的$X_i^{\prime }$,最后給出正交誤差距離向量為：
$$X_i^{\prime \prime }=X_i?X_i^{\prime }(16)$$

橢圓上的正交關(guān)聯(lián)點(diǎn)的Jacobian矩陣

下面我們直接給出橢圓上的正交關(guān)聯(lián)點(diǎn)的Jacobian矩陣，（注：本部分推導(dǎo)比較復(fù)雜，貼出本人的詳細(xì)推導(dǎo)過(guò)程，需要的可以參考推導(dǎo)1,2,3,4,5）
$$J_{X_i^{\prime },\hat{a}}=(R^{?1}{Q}^{?1}B){|}_{x=x_i^{\prime }}(17)$$
這里$Q$見(jiàn)（15），$B=(B_1,B_2,B_3,B_4,B_5)$
此時(shí)

橢圓的正交距離擬合

給定m個(gè)二維點(diǎn)，我們可以利用(16)、（17）給定的誤差距離向量$X_i^{\prime \prime }$和Jacobian矩陣$J_{X_i^{\prime },\hat{a}}$,構(gòu)造p(=2m)個(gè)線性方程組。如下：

,當(dāng)我們使用奇異值分解求解出上述方程組的解時(shí)，就可以進(jìn)行高斯-牛頓迭代。
此時(shí)我們還需要一個(gè)初值，一般可以選擇基于代數(shù)擬合的橢圓或者使用圓作為初始值來(lái)進(jìn)行迭代。通常我們推薦使用圓來(lái)作為初始值參數(shù)。

顯然，從圓擬合中獲得的初始參數(shù)為：
$(X_c,Y_c{)}_{ellipse}=(X_c,Y_c{)}_{circle},a=b=R,\alpha =0$
并且在迭代過(guò)程中，如果$a<b$，則需要令$\alpha \leftarrow \alpha ?sign(\alpha )\frac{\pi }{2}$(保持a為主軸長(zhǎng))

標(biāo)準(zhǔn)坐標(biāo)軸下的橢圓擬合

有時(shí)候，我們也需要為橢圓擬合增加一些約束，經(jīng)典的就是要求在標(biāo)準(zhǔn)坐標(biāo)軸下進(jìn)行橢圓擬合（$\alpha =0或者\pi /2$），或者要求增加橢圓面積約束.此時(shí)，我們只需要在原來(lái)的第（7）式增加如下約束即可。

此時(shí)$w_3,w_4$為權(quán)重參數(shù)，一般可以取$1.0\times 10^6$.

實(shí)驗(yàn)結(jié)果

例1.

我們對(duì)Table 7的離散點(diǎn)進(jìn)行橢圓擬合，其中初始參數(shù)$a_0$源自基于幾何距離的圓擬合。
取$\lambda =1.2$，在經(jīng)過(guò)19次迭代后，最終的校正向量的范數(shù)為$||\Delta a||=4.2\times 10^{?6}$,并且得到誤差指數(shù)${\sigma }_0=1.1719$.見(jiàn)如下：

為了更好地比較，我們也考慮了初始參數(shù)$a_0$源自基于代數(shù)距離的橢圓擬合，而達(dá)到與上述相同的估計(jì)，需要進(jìn)行21次迭代（$||\Delta a||=1.1\times 10^{?6}$）.
我們也比較了我們的結(jié)果與Gander的結(jié)果，而后者達(dá)到相同的估計(jì)，需要進(jìn)行71次迭代（$||\Delta a||=1.0\times 10^{?6}$）.具體可參考Table 8 的III,IV.

Gander的算法有一個(gè)缺點(diǎn)就是不能使用圓的結(jié)果作為初始參數(shù)，為了克服這個(gè)困難，一般取$ a=R, b=R/2$.為了更加直接驗(yàn)證我們算法的魯棒性，我們也考慮了使用兩種設(shè)置作為初始值，分別見(jiàn)Table 8的II，III。即：
II：
$$a=R,b=R/2,\alpha =0$$
III:
$$a=R,b=R/2,\alpha =?1.211$$.
我們分別使用17次和23次迭代達(dá)到了相同的估計(jì)，其中校正向量的范數(shù)分別是$1.2\times 10^{?6},5.2\times 10^{?6}$.

從以上結(jié)果來(lái)看，我們的算法是魯棒的，即使初始值不夠好，最終也能迭代到較好的效果。

創(chuàng)作挑戰(zhàn)賽新人創(chuàng)作獎(jiǎng)勵(lì)來(lái)咯，堅(jiān)持創(chuàng)作打卡瓜分現(xiàn)金大獎(jiǎng)

總結(jié)

以上是生活随笔為你收集整理的基于几何距离的椭圆拟合的全部?jī)?nèi)容，希望文章能夠幫你解決所遇到的問(wèn)題。

如果覺(jué)得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。