python 擬合圓 最小二乘法擬合圓 轉

2021-10-11 16:55:23 字數 4552 閱讀 2404

有一系列的資料點 ,我們知道這些資料點近似的落在乙個圓上,根據這些資料估計這個圓的引數就是乙個很有意義的問題。今天就來講講如何來做圓的擬合。圓擬合的方法有很多種,最小二乘法屬於比較簡單的一種。今天就先將這種。

我們知道圓方程可以寫為:

(x−xc)2+(y−yc)2=r2(x−xc)2+(y−yc)2=r2

通常的最小二乘擬合要求距離的平方和最小。也就是

f=∑((xi−xc)2+(yi−yc)2−−−−−−−−−−−−−−−−−−√−r)2f=∑((xi−xc)2+(yi−yc)2−r)2

最小。這個算起來會很麻煩。也得不到解析解。所以我們退而求其次。

f=∑((xi−xc)2+(yi−yc)2−r2)2f=∑((xi−xc)2+(yi−yc)2−r2)2

這個式子要簡單的多。我們定義乙個輔助函式:

g(x,y)=(x−xc)2+(y−yc)2−r2g(x,y)=(x−xc)2+(y−yc)2−r2

那麼上面的式子可以表示為:

f=∑g(xi,yi)2f=∑g(xi,yi)2

按照最小二乘法的通常的步驟,可知ff 取極值時對應下面的條件。

∂f∂xc=0∂f∂yc=0∂f∂r=0∂f∂xc=0∂f∂yc=0∂f∂r=0

先來化簡 ∂f∂r=0∂f∂r=0

∂f∂r=−2r×∑((xi−xc)2+(yi−yc)2−r2)=−2r×∑g(xi,yi)=0∂f∂r=−2r×∑((xi−xc)2+(yi−yc)2−r2)=−2r×∑g(xi,yi)=0

我們知道半徑 rr 是不能為 0 的。所以必然有:

∑g(xi,yi)=0∑g(xi,yi)=0

這是個非常有用的結論。剩下的兩個式子:

∂f∂xc=−4∑((xi−xc)2+(yi−yc)2−r2)(xi−xc)=−4∑g(xi,yi)(xi−xc)=−4∑xig(xi,yi)=0∂f∂yc=−4∑((xi−xc)2+(yi−yc)2−r2)(yi−yc)=−4∑g(xi,yi)(yi−yc)=−4∑yig(xi,yi)=0∂f∂xc=−4∑((xi−xc)2+(yi−yc)2−r2)(xi−xc)=−4∑g(xi,yi)(xi−xc)=−4∑xig(xi,yi)=0∂f∂yc=−4∑((xi−xc)2+(yi−yc)2−r2)(yi−yc)=−4∑g(xi,yi)(yi−yc)=−4∑yig(xi,yi)=0

也就是下面兩個等式:

∑xig(xi,yi)=0∑yig(xi,yi)=0∑xig(xi,yi)=0∑yig(xi,yi)=0

這裡設:

ui=xi−x¯uc=xc−x¯vi=yi−y¯vc=yc−y¯ui=xi−x¯uc=xc−x¯vi=yi−y¯vc=yc−y¯

其中:x¯=∑xi/ny¯=∑yi/nx¯=∑xi/ny¯=∑yi/n

那麼簡單的推導一下,就會發現:

∑uig(ui,vi)=0∑vig(ui,vi)=0∑uig(ui,vi)=0∑vig(ui,vi)=0

其實都不需要推導,這個變數替換只不過是修改了座標原點的位置。而我們的等式根本就與座標原點的具體位置無關。所以必然成立。

這兩個式子展開寫是這樣的:

∑((ui−uc)2+(vi−vc)2−r2)ui=0∑((ui−uc)2+(vi−vc)2−r2)vi=0∑((ui−uc)2+(vi−vc)2−r2)ui=0∑((ui−uc)2+(vi−vc)2−r2)vi=0

進一步展開:

∑(u3i−2u2iuc+uiu2c+uiv2i−2uivivc+uiv2c−uir2)=0∑(u2ivi−2uiviuc+viu2c+v3i−2v2ivc+viv2c−vir2)=0∑(ui3−2ui2uc+uiuc2+uivi2−2uivivc+uivc2−uir2)=0∑(ui2vi−2uiviuc+viuc2+vi3−2vi2vc+vivc2−vir2)=0

我們知道 ∑ui=0∑ui=0, ∑vi=0∑vi=0。所以上面兩個式子是可以化簡的。

∑(u3i−2u2iuc+uiv2i−2uivivc)=0∑(u2ivi−2uiviuc+v3i−2v2ivc)=0∑(ui3−2ui2uc+uivi2−2uivivc)=0∑(ui2vi−2uiviuc+vi3−2vi2vc)=0

為了簡化式子,我們定義幾個引數:

suuu=∑u3isvvv=∑v3isuu=∑u2isvv=∑v2isuv=∑uivisuuv=∑u2ivisuvv=∑uiv2isuuu=∑ui3svvv=∑vi3suu=∑ui2svv=∑vi2suv=∑uivisuuv=∑ui2visuvv=∑uivi2

那麼上面的式子可以寫為:

suuuc+suvvc=suuu+suvv2suvuc+svvvc=suuv+svvv2suuuc+suvvc=suuu+suvv2suvuc+svvvc=suuv+svvv2

至此,就可以解出ucuc 和vcvc 了。

uc=suuvsuv−suuusvv−suvvsvv+suvsvvv2(s2uv−suusvv)vc=−suusuuv+suuusuv+suvsuvv−suusvvv2(s2uv−suusvv)uc=suuvsuv−suuusvv−suvvsvv+suvsvvv2(suv2−suusvv)vc=−suusuuv+suuusuv+suvsuvv−suusvvv2(suv2−suusvv)

那麼:xc=uc+x¯yc=vc+y¯xc=uc+x¯yc=vc+y¯

還剩下個 rr 沒求出來, 也很簡單。

∑g(xi,yi)=0∑((xi−xc)2+(yi−yc)2−r2)=0∑g(xi,yi)=0∑((xi−xc)2+(yi−yc)2−r2)=0

所以:r2=∑((xi−xc)2+(yi−yc)2)r2=∑((xi−xc)2+(yi−yc)2)

好了。下面給出個**,這個**的具體公式和我這裡給出的有一點小差異,但是原理是相同的。

* 最小二乘法擬合圓

* 擬合出的圓以圓心座標和半徑的形式表示

* 此**改編自 newsmth.net 的 jingxing 在 graphics 版貼出的**。

* 版權歸 jingxing, 我只是搬運工外加一些簡單的修改工作。

typedef complex point;

bool circleleastfit(const std::vector &points, double ¢er_x, double ¢er_y, double &radius)

center_x = 0.0f;

center_y = 0.0f;

radius = 0.0f;

if (points.size() < 3)

return false;

double sum_x = 0.0f, sum_y = 0.0f;

double sum_x2 = 0.0f, sum_y2 = 0.0f;

double sum_x3 = 0.0f, sum_y3 = 0.0f;

double sum_xy = 0.0f, sum_x1y2 = 0.0f, sum_x2y1 = 0.0f;

int n = points.size();

for (int i = 0; i < n; i++)

double x = points[i].real();

double y = points[i].imag();

double x2 = x * x;

double y2 = y * y;

sum_x += x;

sum_y += y;

sum_x2 += x2;

sum_y2 += y2;

sum_x3 += x2 * x;

sum_y3 += y2 * y;

sum_xy += x * y;

sum_x1y2 += x * y2;

sum_x2y1 += x2 * y;

double c, d, e, g, h;

double a, b, c;

c = n * sum_x2 - sum_x * sum_x;

d = n * sum_xy - sum_x * sum_y;

e = n * sum_x3 + n * sum_x1y2 - (sum_x2 + sum_y2) * sum_x;

g = n * sum_y2 - sum_y * sum_y;

h = n * sum_x2y1 + n * sum_y3 - (sum_x2 + sum_y2) * sum_y;

a = (h * d - e * g) / (c * g - d * d);

b = (h * c - e * d) / (d * d - g * c);

c = -(a * sum_x + b * sum_y + sum_x2 + sum_y2) / n;

center_x = a / (-2);

center_y = b / (-2);

radius = sqrt(a * a + b * b - 4 * c) / 2;

return true;

下圖是個實際測試的結果。對這種均勻分布的雜訊資料計算的結果還是很準確的。

但是當資料中有部分偏向某乙個方向的干擾資料時。結果就不那麼樂觀了。下圖就很說明問題。

資料點中有 20% 是干擾資料。擬合出的圓就明顯被拽偏了。

之所以出現這個問題就是因為最小二乘擬合的平方項對離群點非常敏感。解決這個問題就要用其他的擬合演算法,比如用距離之和作為擬合判據。等我有空了再寫一篇部落格介紹其他幾種方法。

python最小二乘法擬合圓 最小二乘法擬合圓

有一系列的資料點 我們知道這些資料點近似的落在乙個圓上。依據這些資料預計這個圓的引數就是乙個非常有意義的問題。今天就來講講怎樣來做圓的擬合。圓擬合的方法有非常多種,最小二乘法屬於比較簡單的一種。今天就先將這樣的。我們知道圓方程能夠寫為 x?xc 2 y?yc 2 r2 通常的最小二乘擬合要求距離的平...

python最小二乘法擬合圓 最小二乘法擬合圓

有一系列的資料點 我們知道這些資料點近似的落在乙個圓上。依據這些資料預計這個圓的引數就是乙個非常有意義的問題。今天就來講講怎樣來做圓的擬合。圓擬合的方法有非常多種,最小二乘法屬於比較簡單的一種。今天就先將這樣的。我們知道圓方程能夠寫為 x?xc 2 y?yc 2 r2 通常的最小二乘擬合要求距離的平...

改進最小二乘法缺陷圓輪廓擬合

機器視覺影象處理中常用到圓形檢測,準確地擬合圓輪廓是測量圓直徑尺寸的前提和保證。霍夫圓檢測容易受到變形圓或缺損圓的影響,改進的最小二乘法可以提高擬合和尺寸測量精度。1 公式推導 最小二乘法通過最小化誤差的平方和找到一組資料的最佳函式匹配。最小二乘法是用最簡的方法求得一些絕對不可知的真值,而令誤差平方...