問題

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

方法

通常我們使用最小二乘法求解如下的最優化問題：

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

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

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

論文的核心思路

顯然，由于離散點到橢圓的幾何距離是非線性的，因此橢圓擬合是一種非線性優化的問題，需要通過多次迭代求解。

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

非線性最小二乘擬合

問題

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

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

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

我們將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$
并且，我們根據（1）可以寫出迭代算法的性能表現指數，其實就是參差的表達式：
$${\sigma }_0^2=||X?F(a)|{|}^2=[X?F(a){]}^T[X?F(a)](5)$$

求解

接下來，我們需要求解出$\Delta a$,才能進行迭代。

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

基于幾何距離的橢圓擬合

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

對于橢圓的最小二乘正交距離擬合，我們將使用第一段所講的方法。此時，我們定義
$$\hat{a}=(X_c,Y_c,a,b,\alpha {)}^T$$
$X$可以看成給定的離散點$X_i$，而自然$F(\hat{a})$就是定點$X_i$在橢圓上的最近點$X_i^{\prime }$(也稱為正交關聯點)。迭代公式就變成了如下：
$$\{\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指給定的離散點的個數。各個離散點均滿足（7）,因此關于（7）可以聯立。實際上關于（7）的方程是2m個，而未知數a的維數是5，顯然2m>>5,因此解是不唯一的，我們需要求出最小二乘解即可。

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

坐標系轉換

為了便于后面的求解，我們需要考慮,將原來的坐標系XY，旋轉$\alpha$角，建立一個位于$(0,0)$的臨時坐標系xy。見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 )$

橢圓上的正交關聯點

經過坐標軸轉換，5個橢圓參數變成了2個，僅僅包含了主軸$a,b$。橢圓上的點可以用標準方程描述如下：
$$\frac{x^2}{a^2}+\frac{y^2}{b^2}=1(10)$$
從Fig.3上，我們可以看到位于xy坐標系上的正交關聯點$（x_i^{\prime },y_i^{\prime }）$的切向量與兩點（即$（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}$$
上面兩式稱為正交關聯條件,我們將使用廣義牛頓法來求解上面方程組的解。
即使用如下的迭代公式求解：
$$\{\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)$$
對于迭代公式（14）,我們提供初始解$x_0$,如下：(見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}$$

一般經過3-4輪迭代就可以提供了一個很高精度的解。

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

橢圓上的正交關聯點的Jacobian矩陣

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

橢圓的正交距離擬合

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

,當我們使用奇異值分解求解出上述方程組的解時，就可以進行高斯-牛頓迭代。
此時我們還需要一個初值，一般可以選擇基于代數擬合的橢圓或者使用圓作為初始值來進行迭代。通常我們推薦使用圓來作為初始值參數。

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

標準坐標軸下的橢圓擬合

有時候，我們也需要為橢圓擬合增加一些約束，經典的就是要求在標準坐標軸下進行橢圓擬合（$\alpha =0或者\pi /2$），或者要求增加橢圓面積約束.此時，我們只需要在原來的第（7）式增加如下約束即可。

此時$w_3,w_4$為權重參數，一般可以取$1.0\times 10^6$.

實驗結果

例1.

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

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

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

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

創作挑戰賽新人創作獎勵來咯，堅持創作打卡瓜分現金大獎

總結

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

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。