直線と同じように、最小二乗法を用いて点P1,P2, ... ,Pnから円を近似することができます。一般的な手法は他の記事を参照してもらうことにして、ここでは始点P1と終点Pnを必ず通り、かつ誤差が小さくなるような円を求めます。また、そのような関数をpython(numpy)で実装してみます。
方針
直線と同様に、各点に対して理論値と測定値の差（残差）を求め、それらの平方和を小さくするような値を求めます。
一般に円を表すためには中心座標$(x,y)$、半径$r$の3つのパラメータが必要ですが、2点$P_1,P_n$を通る円の中心は必ずこの2点の垂直二等分線上にあります。中心が決まれば半径も決まりますから、今回求めるパラメータは一つにできます。
方法
円の中心と半径をベクトルを用いてパラメータ表示し……などとやっても不可能ではないですが、計算が煩雑になってしまうので、別の方法を取ります。ひとまず$P_1(-a,0),P_n(a,0)$の場合を考えましょう。このとき、中心はy軸上の点$(0,y_0)$、半径$r^2 = y_0^2+a^2$となり、求める円の方程式は
x^2+y^2-2yy_0-a^2=0
となります。この$y_0$を求めます。
各点に対する残差は、各点から中心までの距離と半径の差ですが、計算の簡単のためにあらかじめ自乗してから差を取ったものを考えます。つまり各点の残差は上の式に$(x_k,y_k)$を代入して
x_k^2+y_k^2-2y_ky_0-a^2
です。これの平方和、すなわち
\sum \{ x_k^2+y_k^2-2y_ky_0-a^2 \} ^2
を最小にするような$y_0$を求めます。この式は$y_0$の二次関数ですから、$y_0$の偏微分を0とおいて$y_0$を求めます。
\begin{align}
\sum -4y_k( x_k^2+y_k^2-2y_ky_0-a^2 ) &= 0 \\
\sum x_k^2y_k+\sum y_k^3-2y_0\sum y_k^2-a^2\sum y_k &= 0 \\
\end{align}
よって$y_0$は以下です。
y_0 = \frac{\sum x_k^2y_k+\sum y_k^3-a^2\sum y_k}{2\sum y_k^2}
座標変換
さて、これを任意の点群に対して適用するには、点群に対して点$P_1(x_1,y_1)$が$(-a,0)$へ、点$P_n(x_n,y_n)$が$(a,0)$へ移るような変換をしないといけません。そのためには、
①$P_1$と$P_n$の中点$P_c$が原点に移るような平行移動を行う
②$P_1$と$P_n$がx軸上に移るような回転を行う
の2ステップが必要です。細かい計算は省きますが、次の変換を行えばいいことがわかります。
\left( \begin{array}{c} x' \\ y' \end{array} \right) =
-\frac{1}{a}\left(
\begin{array}{cc}
x_1-x_c & -(y_1-y_c)\\
y_1-y_c & x_1-x_c
\end{array}
\right)
\left( \begin{array}{c} x-x_c \\ y-y_c \end{array} \right)
ただし、$x_c,y_c$は$P_c$の座標、$a$は$|P_1P_c|$です。
この変換をすべての点に行ってから近似を行い、逆変換すれば中心・半径が求まります。
実装
pythonで実装してみました。
def approxByArc(points):
start = points[0]
end = points[-1]
sample = points[1:-1]
midpoint = (start + end) /2
start2 = start - midpoint
a = np.linalg.norm(start2)
x = start2[0]
y = start2[1]
rotateMatrix = np.array([[x,y],[-y,x]]) * -1 / a
def conversion(point):
point = point - midpoint
p = np.array([[point[0]], [point[1]]])
p = rotateMatrix @ p
p = np.array([p[0][0],p[1][0]])
return p
sample = np.apply_along_axis(conversion, 1, sample)
xs = np.array([i[0] for i in sample])
ys = np.array([i[1] for i in sample])
p1 = np.sum(xs*xs*ys)
p2 = np.sum(ys*ys*ys)
p3 = np.sum(ys) * a * a
p4 = np.sum(ys*ys) * 2
y0 = (p1+p2-p3)/p4
center = np.array([[0],[y0]])
center = np.linalg.inv(rotateMatrix) @ center
center = np.array([center[0][0], center[1][0]])
center = center + midpoint
radius = np.linalg.norm(center - start)
return center,radius
ちゃんと動くか試してみましょう。
points = []
points.append([60,50])
n = 40
for i in range(1,n):
angle = math.pi / n * i
radius = 10 + random.uniform(-0.4,0.4)
x = 50 + radius * math.cos(angle)
y = 50 + radius * math.sin(angle)
points.append([x,y])
points.append([40,50])
r,c = approxByArc(np.array(points))
print(r)
print(c)
[50. 49.99319584]
10.00000231483068
うまくいきました。
おわりに
単純な座標に変換してから計算する方法はプログラムがわかりやすくなる反面、手計算には向きません。ですが、このような複雑なステップを手計算する人はいないので大丈夫です。たぶん。