最近工作中遇到一個擬合圓的問題，通過找到的輪廓點（存在缺失的情況），需要找出圓心及半徑。這裡採用最小二乘法進行擬合，並記錄一下具體的推導過程。最小二乘法是解決曲線擬合問題最常用的方法，其基本思路是：令

\[f(x) = \alpha_1 \phi_1(x) + \alpha_2 \phi_2(x) + ... + \alpha_m \phi_m(x) \]

其中\(\phi_k(x)\)是事先選擇的一組線性無關函式，\(\alpha_k\)是待定係數 ，擬合準則是使\(y_i(i = 1, 2, ..., n)\)與\(f(x_i)\)的距離\(\delta_i\)的平方和最小。注意：最小二乘擬合對於離群點很敏感，如果存在離群點可以採用RANSAC(Random Sample Consensus)演算法進行擬合

。

最小二乘法擬合圓推導流程

1、確定待定係數

設圓的圓心為\((A, B)\)，半徑為\(R\)，則圓的方程可以寫成

\[ \begin{align} R^2 &= (x-A)^2 + (y - B)^2 \tag{1}\\ &= x^2 - 2Ax + A^2 + y^2 - 2By + B^2 \tag{2} \end{align} \]

令\(a = -2A, b = -2B, c = A^2 + B^2 - R^2\)，則圓的曲線方程可寫成(3)式：

\[ \begin{align} x^2 + y^2 +ax + by + c = 0 \tag{3} \end{align} \]

只需要求出引數\(a, b, c\)

即可得到圓心及半徑：

\[ \begin{align} A &= -\frac{a}{2} \tag{4}\\ B &= -\frac{b}{2} \tag{5}\\ R &= \frac{1}{2}\sqrt{a^2 + b^2 - 4c} \tag{6} \end{align} \]

所以這裡待定係數為\(a, b, c\)。

2、確定距離和

樣本集\((X_i, Y_i), i = {1, 2, 3, ..., N}\)中第\(i\)個樣本到圓心的距離記為\(d_i\),

\[{d_i}^2 = (X_i - A)^2 + (Y_i - B)^2 \tag{7} \]

則其與半徑之間的誤差為：

\[\delta_i = {d_i}^2 - R^2 = {X_i}^2 + {Y_i}^2 + aX_i + bY_i + c \tag{8} \]

所以所有樣本的誤差和\(Q(a, b, c)\)為

\[ \begin{align} Q(a, b, c) &= \sum_{i = 1}^{N}{\delta_i}^2 \tag{9} \\ &= \sum_{i = 1}^{N}{({X_i}^2 + {Y_i}^2 + aX_i + bY_i + c)}^2 \tag{10} \end{align} \]

現在問題變成求引數\(a,b,c\)使得誤差和\(Q(a, b, c)\)最小。

3、求解

\(Q(a, b, c)\)的每一項都大於等於0，所以函式存在大於或等於零的極小值，極大值為正無窮大。
\(Q(a, b, c)\)分別對引數\(a,b,c\)求偏導，然後令偏導數為0，即可求出極值點。

\[ \begin{align} \frac{{\rm d}Q(a, b, c)}{{\rm d}a} &= \sum_{i = 1}^{N} 2({X_i}^2 + {Y_i}^2 + aX_i + bY_i + c)X_i = 0 \tag{11} \\ \frac{{\rm d}Q(a, b, c)}{{\rm d}b} &= \sum_{i = 1}^{N}2({X_i}^2 + {Y_i}^2 + aX_i + bY_i + c)Y_i = 0 \tag{12} \\ \frac{{\rm d}Q(a, b, c)}{{\rm d}c} &= \sum_{i = 1}^{N}2({X_i}^2 + {Y_i}^2 + aX_i + bY_i + c) = 0 \tag{13} \end{align} \]

$式(11) * N - 式(13) * \sum_{i = 1}^{N} X_i $消去c

\[ \begin{align} 0 &= N\sum_{i = 1}^{N} 2({X_i}^2 + {Y_i}^2 + aX_i + bY_i + c)X_i - \sum_{i = 1}^{N}2({X_i}^2 + {Y_i}^2 + aX_i + bY_i + c) * \sum_{i = 1}^{N}X_i \tag{14} \\ &= N\sum_{i = 1}^{N} 2({X_i}^2 + {Y_i}^2 + aX_i + bY_i)X_i - \sum_{i = 1}^{N}2({X_i}^2 + {Y_i}^2 + aX_i + bY_i) * \sum_{i = 1}^{N}X_i \tag{15}\\ &= (N\sum_{i = 1}^{N}{X_i}^2 - \sum_{i = 1}^{N}X_i \sum_{i = 1}^{N}X_i)a + (N\sum_{i = 1}^{N}X_i Y_i - \sum_{i = 1}^{N}X_i \sum_{i = 1}^{N}Y_i)b + N\sum_{i = 1}^{N}{X_i}^3 + N\sum_{i = 1}^{N}X_i{Y_i}^2 - N\sum_{i = 1}^{N}({X_i}^2 + {Y_i}^2)\sum_{i = 1}^{N}X_i \tag{16} \end{align} \]

$式(12) * N - 式(13) * \sum_{i = 1}^{N} Y_i $得

\[ \begin{align} 0 &= N\sum_{i = 1}^{N} 2({X_i}^2 + {Y_i}^2 + aX_i + bY_i + c)Y_i - \sum_{i = 1}^{N}2({X_i}^2 + {Y_i}^2 + aX_i + bY_i + c) * \sum_{i = 1}^{N}Y_i \tag{17} \\ &= N\sum_{i = 1}^{N} 2({X_i}^2 + {Y_i}^2 + aX_i + bY_i)Y_i - \sum_{i = 1}^{N}2({X_i}^2 + {Y_i}^2 + aX_i + bY_i) * \sum_{i = 1}^{N}Y_i \tag{18}\\ &= (N\sum_{i = 1}^{N}{X_i}{Y_i} - \sum_{i = 1}^{N}X_i \sum_{i = 1}^{N}Y_i)a + (N\sum_{i = 1}^{N}{Y_i}^2 - \sum_{i = 1}^{N}Y_i \sum_{i = 1}^{N}Y_i)b + N\sum_{i = 1}^{N}{Y_i}^3 + N\sum_{i = 1}^{N}{X_i}^2{Y_i} - N\sum_{i = 1}^{N}({X_i}^2 + {Y_i}^2)\sum_{i = 1}^{N}Y_i \tag{19} \end{align} \]

令

\[ \begin{align} C &= N\sum_{i = 1}^{N} {X_i}^2 - \sum_{i = 1}^{N}X_i \sum_{i = 1}^{N}X_i \tag{20} \\ D &= N\sum_{i = 1}^{N} {X_i}{Y_i} - \sum_{i = 1}^{N}X_i \sum_{i = 1}^{N}Y_i \tag{21} \\ E &= N\sum_{i = 1}^{N} {X_i}^3 + N\sum_{i = 1}^{N} {X_i}{Y_i}^2 - \sum_{i = 1}^{N}({X_i}^2 + {Y_i}^2) \sum_{i = 1}^{N}X_i \tag{22} \\ G &= N\sum_{i = 1}^{N} {Y_i}^2 - \sum_{i = 1}^{N}Y_i \sum_{i = 1}^{N}Y_i \tag{23} \\ H &= N\sum_{i = 1}^{N} {X_i}^2Y_i + N\sum_{i = 1}^{N}{Y_i}^3 - \sum_{i = 1}^{N}({X_i}^2 + {Y_i}^2) \sum_{i = 1}^{N}Y_i \tag{24} \end{align} \]

可得：

\[ \begin{align} 0 &= Ca + Db + E \tag{25} \\ 0 &= Da + Gb + H \tag{26} \\ a &= \frac{HD - EG} {CG - D^2} \tag{27} \\ b &= \frac{HC - ED} {D^2 - GC} \tag{28} \\ c &= -\frac{\sum_{i = 1}^{N}({X_i}^2 + {Y_i}^2) + a\sum_{i = 1}^{N}X_i + b\sum_{i = 1}^{N}Y_i}{N} \tag{29} \end{align} \]

再結合式(4)(5)(6)即可求得圓心和半徑。

程式碼

參考文件

圓擬合算法
最小二乘法擬合圓公式推導及其實現
最小二乘法擬合圓