最小二乗法による円の近似

icon 項目のみ表示/展開表示の切り替え

最小二乗法による円の近似方法

n個の円周付近の点xi,yiから中心及び半径を最小二乗法により近似してみます。
各点にばらつきがある場合は誤差がある場合に有効な方法となります。
最小二乗法は読んで字のごとく差の二乗の和が最小となるように係数を計算する方法です。
ただしあくまでも近似ですので、さらに精度を求めるときは最小二乗法の値をもとにさらに別の方法で近似する必要がるかもしれません。
中心(a,b)で半径r、円周の任意の点を(xi,yi)とすると円は以下の式で表せます。
\sqrt{(x_i-a)^2+(y_i-b)^2}=r
rと中心とxi,yi間の距離の差を出す式をとりあえず差を0として作成します。
(x_i-a)^2+(y_i-b)^2=r^2
(x_i-a)^2+(y_i-b)^2-r^2=0

上記の式を二乗すると
\bigl\{ (x_i-a)^2+(y_i-b)^2-r^2 \bigr\}^2=0
(c-d)^2=c^2-2cd+d^2  左式を使って上式からカッコを取り除きシンプルにします。
例えば(c-d)^2  をc^2+Dc+E  に置き換えるためには、
D=-2d,E=d^2  とすれば可能です。従って下式の様に変形できます。
\bigl\{ {x_i}^2+{y_i}^2+Ax_i+By_i+C \bigr\}^2=0
A=-2a
B=-2b
C=a^2+b^2-r^2

上式を二乗してi=1～nについて合計しA,B,Cについてそれぞれ偏微分します。
合成関数の微分の公式と指数の微分の式を使用します。
y=f(u),u=g(x),y=f(g(x))
\displaystyle \frac{dy}{dx}=\frac{dy}{du} \cdot \frac{du}{dx}
(x^a)'=ax^{a-1}
Aについて偏微分
\displaystyle \frac{\partial}{\partial A}=\sum_{i=1}^n{2x_i({x_i}^2+{y_i}^2+Ax_i+By_i+C)}
=\sum_{i=1}^n{2({x_i}^3+x_i {y_i}^2+A{x_i}^2+Bx_i y_i+Cx_i)}
Bについて偏微分
\displaystyle \frac{\partial}{\partial B}=\sum_{i=1}^n{2y_i({x_i}^2+{y_i}^2+Ax_i+By_i+C)}
=\sum_{i=1}^n{2(y_i{x_i}^2+{y_i}^3+Ax_i y_i+B{y_i}^2+Cy_i)}
Cについて偏微分
\displaystyle \frac{\partial}{\partial C}=\sum_{i=1}^n{2({x_i}^2+{y_i}^2+Ax_i+By_i+C)}
A,B,CについてA,B,C=の形式に式を並び替える。
\displaystyle \frac{\partial}{\partial A}=\sum_{i=1}^n{A{x_i}^2}+\sum_{i=1}^n{Bx_iy_i}+\sum_{i=1}^n{Cx_i}=-\sum_{i=1}^n{{x_i}^3}-\sum_{i=1}^n{x_i{y_i}^2}
\displaystyle \frac{\partial}{\partial B}=\sum_{i=1}^n{Ax_i y_i +\sum_{i=1}^n{B{y_i}^2}}+\sum_{i=1}^n{Cy_i}=-\sum_{i=1}^n{y_i {x_i}^2}-\sum_{i=1}^n{{y_i}^3}
\displaystyle \frac{\partial}{\partial C}=\sum_{i=1}^n{Ax_i}+\sum_{i=1}^n{By_i}+\sum_{i=1}^n{C}=-\sum_{i=1}^n{{x_i}^2}-\sum_{i=1}^n{{y_i}^2}
行列にしてA,B,Cを解く。
\left( \begin{array}{lll} \sum_{i=1}^n{{x_i}^2} & \sum_{i=1}^n{x_i y_i} & \sum_{i=1}^n{x_i} \\ \sum_{i=1}^n{x_i y_i} & \sum_{i=1}^n{{y_i}^2} & \sum_{i=1}^n{y_i} \\ \sum_{i=1}^n{x_i} & \sum_{i=1}^n{y_i} & \sum_{i=1}^n{1} \end{array} \right) \left( \begin{array}{lll} A \\ B \\ C \end{array} \right)=\left( \begin{array}{lll} -\sum_{i=1}^n{({x_i}^3-x_i {y_i}^2)} \\ -\sum_{i=1}^n{({x_i}^2y_i-{y_i}^3)} \\ -\sum_{i=1}^n{({x_i}^2-{y_i}^2)} \end{array} \right)
A,B,C=に変形すると
\left( \begin{array}{lll} A \\ B \\ C \end{array} \right)={\left( \begin{array}{lll} \sum_{i=1}^n{{x_i}^2} & \sum_{i=1}^n{x_i y_i} & \sum_{i=1}^n{x_i} \\ \sum_{i=1}^n{x_i y_i} & \sum_{i=1}^n{{y_i}^2} & \sum_{i=1}^n{y_i} \\ \sum_{i=1}^n{x_i} & \sum_{i=1}^n{y_i} & \sum_{i=1}^n{1} \end{array} \right) }^{-1} \left( \begin{array}{lll} -\sum_{i=1}^n{({x_i}^3-x_i {y_i}^2)} \\ -\sum_{i=1}^n{({x_i}^2y_i-{y_i}^3)} \\ -\sum_{i=1}^n{({x_i}^2-{y_i}^2)} \end{array} \right)
行列に-1乗と表現しているものは指数ではなく逆行列というものです。イメージは逆数です。
A \cdot A^{-1}=1

A=\begin{pmatrix} a & b & c \\ d & e & f \\ g & h & i \end{pmatrix}
\displaystyle A^{-1}=\frac{1}{a e i+b f g + c d h - c e g - b d i - a f h} \begin{pmatrix} e i-f h & -(b i-c h) & b f-c e \\ -(d i-f g) & a i-c g & -(a f-c d) \\ d h-e g & -(a h-b g) & a e -b d \\ \end{pmatrix}

行列の積
\left( \begin{array}{lll} a&b&c\\ d&e&f\\ g&h&i \end{array} \right) \left( \begin{array}{lll} j\\ k\\ l \end{array} \right) =\left( \begin{array}{lll} cl+bk+aj \\ fl+ek+dj \\ il+hk+gj \end{array} \right)

行列の計算をするとA,B,Cが算出できます。
これよりa,b,rが算出できます。

最小二乗法による円の近似の計算例

以下のフォームに座標の一覧を入力して計算ボタンをクリックすると円の近似の結果が表示されます。
入力された各点は黒色の小さい円、近似により算出された円は黒い大きい円、赤色の円が近似された円の中心です。
下表は計算の途中結果を表示しています。一番下の行が合計となります。

座標リスト



SVGの代替画像
点番号を表示
i Xi Yi xi2 yi2 xi3 yi3 xiyi xi2y xiyi2 ri

\left( \begin{array}{lll} A \\ B \\ C \end{array} \right)={\left( \begin{array}{lll} \sum_{i=1}^n{{x_i}^2} & \sum_{i=1}^n{x_i y_i} & \sum_{i=1}^n{x_i} \\ \sum_{i=1}^n{x_i y_i} & \sum_{i=1}^n{{y_i}^2} & \sum_{i=1}^n{y_i} \\ \sum_{i=1}^n{x_i} & \sum_{i=1}^n{y_i} & \sum_{i=1}^n{1} \end{array} \right) }^{-1} \left( \begin{array}{lll} -\sum_{i=1}^n{({x_i}^3-x_i {y_i}^2)} \\ -\sum_{i=1}^n{({x_i}^2y_i-{y_i}^3)} \\ -\sum_{i=1}^n{({x_i}^2-{y_i}^2)} \end{array} \right)