【python】最小二乗法で円のフィッテングを行う方法を解説

円の最小二乗近似もやってみる
- 式変形
- Pythonでやってみる
  - サンプルデータ
  - 最小二乗法を使って近似
終わりに

円の最小二乗近似もやってみる

前回は1次関数、2次関数の最小二乗近似を行いました。
mori-memo.hateblo.jp
今回はそれを円に適応してみます。
1次関数や2次関数、指数関数の場合は、Excel標準の機能でグラフを書いたら簡単につけてくれますが、
円で近似したいときってありませんか？
Excelの標準機能ではそこまで対応してませんが、pythonでは簡単に実装できますのでやってみます。
意外と円の最小二乗近似はよく使われていて、
製品の内径にセンサを当てて、そのデータから製品の内径寸法を求めたり、径の中心座標を求めたりするのに使われたりします。
簡単に原理から見ていきましょう。

式変形

$x, y$ が測定データだとして、データを以下の円の式に近似すればいいわけです。

$(x - c_{x})^2 + (y - c_{y})^2 = r^2$

前回の記事を参考にすると、
$u$ を近似式の係数の行列とした場合に以下のような式に変形できれば係数を求められそうです。

$Au = v$

円の式を展開してみると、以下の形に書き直せます。

$x^2 + y^2 + ax + by + c = 0$

ここで、 $a, b, c$ は下記のとおりです。

$a = -2c_{x}$
$b = -2c_{y}$
$c = c_{x}^2 + x_{y}^2 - r^2$

この $a, b, c$ が求まればいいので、 $n$ 個の実験データ( $x_{1}, x_{2},...,x_{n}, y_{1}, y_{2},...,y_{n}$ )を用いて以下のように書き直すことができます。

$\begin{pmatrix} x_{1} & y_{1} & 1 \\ x_{2} & y_{2} & 1 \\ \vdots & \vdots & \vdots \\ x_{n} & y_{n} & 1 \end{pmatrix} \begin{pmatrix} a \\ b \\ c \end{pmatrix} = -\begin{pmatrix} x_{1}^2 + y_{1}^2 \\ x_{2}^2 + y_{2}^2 \\ \vdots \\ x_{n}^2 + y_{n}^2 \end{pmatrix}$

これで $Au = v$ の形になりました。つまり、以下のように代入すればいいわけです。

$A = \begin{pmatrix} x_{1} & y_{1} & 1 \\ x_{2} & y_{2} & 1 \\ \vdots & \vdots & \vdots \\ x_{n} & y_{n} & 1 \end{pmatrix}, u = \begin{pmatrix} a \\ b \\ c \end{pmatrix}, v = -\begin{pmatrix} x_{1}^2 + y_{1}^2 \\ x_{2}^2 + y_{2}^2 \\ \vdots \\ x_{n}^2 + y_{n}^2 \end{pmatrix}$

ここまでくればあとは前回と同じです。

$u = (A^T A)^{-1}A^T v$

上記で $u$ が求められましたね！
それではプログラムで確認してみましょう。

Pythonでやってみる

サンプルデータ

テキトーに作ります。
今回は $c_{x} = 10, c_{y} = 3, r = 5$ の円にノイズを加えてサンプルデータとしました。

import numpy as np
import matplotlib.pyplot as plt

# 測定データ
theta = np.linspace(-np.pi, np.pi, 50)
r = 5 + np.random.normal(0, 0.2, (50))
cx = 10
cy = 3
x = a + r * np.cos(theta)
y = b + r * np.sin(theta)

# 測定データのプロット
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.scatter(x, y, label='measure data')
ax.legend(loc='best')
ax.set_aspect('equal')
plt.show()

いい感じに近似できそうなデータができました。
ラベルが中央過ぎて気持ち悪いですが……

これを最小二乗法で円に近似しましょう。

最小二乗法を使って近似

A = np.vstack((x, y, np.ones((len(x))))).T
v = -(x ** 2 + y ** 2)
u, residuals, rank, s = np.linalg.lstsq(A, v, rcond=None)

甘ちゃんなので、np.linalg.lstsq使って解きます。
この返り値としてuに $a, b, c$ がリスト形式で入ってます。

ここで、 $a, b, c$ は下記の式だったのでそこから $c_{x}, c_{y}, r$ を求めます。

$a = -2c_{x}$
$b = -2c_{y}$
$c = c_{x}^2 + x_{y}^2 - r^2$

# 係数を求める
cx_pred = u[0] / (-2)
cy_pred = u[1] / (-2)
r_pred = np.sqrt(cx_pred ** 2 + cy_pred ** 2 - u[2])
print(cx_pred, cy_pred, r_pred)
# 9.995265931036581 2.9946774110049157 4.9752197197960495

サンプルデータは $c_{x} = 10, c_{y} = 3, r = 5$ だったので、なかなか近い値にはなってる気がします。
プロットします

# 測定データのプロット
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.scatter(x, y, label='measure data')
ax.plot(cx_pred + r_pred * np.cos(theta), cy_pred + r_pred * np.sin(theta), 'r', label='least square')
ax.legend(loc='best')
ax.set_aspect('equal')
plt.show()