2次元の点群にフィットする円を求める - 最小二乗円

はじめに

点群にフィットする円のことを「最小二乗円」というらしいです。
その最小二乗円を以下の3つのパターンで求めたいと思います。
・普通の最小二乗円(任意の点を通らない)
・ある1点を通る最小二乗円
・ある2点を通る最小二乗円
点群は $\{(x_i,y_i)\in{\mathbb R}^2\}_{i=1}^N$ とします。

イメージしやすいように、以下の図に点群の点の数が4の場合の最小二乗円を描画しています。
緑色が「普通の最小二乗円」、赤色が「ある1点を通る最小二乗円」、青色が「ある2点を通る最小二乗円」です。

普通の最小二乗円

重心 $(\bar{x},\bar{y})$ を求めます。

$\begin{eqnarray} \bar{x}=\dfrac{1}{N}\sum_{i=1}^Nx_i,\ \bar{y}=\dfrac{1}{N}\sum_{i=1}^Ny_i.\tag{1} \end{eqnarray}$

$(-\bar{x},-\bar{y})$ だけ平行移動した座標系の点群を $\{(X_i,Y_i)\}_{i=1}^N$ とします。

$\begin{eqnarray} X_i=x_i-\bar{x},\ Y_i=y_i-\bar{y}.\tag{2} \end{eqnarray}$

$X_i,Y_i$ の総和は以下の関係があります。

$\begin{eqnarray} &&\sum_{i=1}^NX_i=\sum_{i=1}^N(x_i-\bar{x})=N\bar{x}-N\bar{x}=0,\\ &&\sum_{i=1}^NY_i=\sum_{i=1}^N(y_i-\bar{y})=N\bar{y}-N\bar{y}=0.\tag{3} \end{eqnarray}$

変換前、変換後の座標系における円の中心座標をそれぞれ $(c_x,c_y),\ (C_x,C_y)$ とすると、以下の式が成り立ちます。

$\begin{eqnarray} C_x=c_x-\bar{x},\ C_y=c_y-\bar{y}.\tag{4} \end{eqnarray}$

円の半径を $r$ とし、 $R=r^2$ とおきます。
目的関数 $E$ を以下のように定めます。

$\begin{eqnarray} E=\sum_{i=1}^N\left((X_i-C_x)^2+(Y_i-C_y)^2-R\right)^2\tag{5} \end{eqnarray}$

式を簡単にするために $T_{mn}=\displaystyle\sum_{i=1}^NX_i^mY_i^n\ (m,n=0,1,2,3)$ とおきます。

$E$ を $R$ で微分して、 $=0$ とおきます。

$\begin{eqnarray} &&\frac{\partial E}{\partial R}=0\\ &&\Rightarrow\frac{\partial}{\partial R}\sum_{i=1}^N\left((X_i-C_x)^2+(Y_i-C_y)^2-R\right)^2=0\\ &&\Rightarrow-2\sum_{i=1}^N\left((X_i-C_x)^2+(Y_i-C_y)^2-R\right)=0\\ &&\Rightarrow\sum_{i=1}^N\left((X_i-C_x)^2+(Y_i-C_y)^2-R\right)=0\tag{6}\\ &&\Rightarrow\sum_{i=1}^NX_i^2-2C_x\underbrace{\sum_{i=1}^NX_i}_{=0\ \because eq.(3)}+NC_x^2+\sum_{i=1}^NY_i^2-2C_y\underbrace{\sum_{i=1}^NY_i}_{=0\ \because eq.(3)}+NC_y^2-NR=0\\ &&\Rightarrow\sum_{i=1}^NX_i^2+NC_x^2+\sum_{i=1}^NY_i^2+NC_y^2-NR=0\\ &&\Rightarrow R=C_x^2+C_y^2+\dfrac{1}{N}\left(\sum_{i=1}^NX_i^2+\sum_{i=1}^NY_i^2\right)\\ &&\Rightarrow R=C_x^2+C_y^2+\dfrac{1}{N}(T_{20}+T_{02})\tag{7} \end{eqnarray}$

$E$ を $C_x$ で微分して、 $=0$ とおきます。

$\begin{eqnarray} &&\frac{\partial E}{\partial C_x}=0\\ &&\Rightarrow\frac{\partial}{\partial C_x}\sum_{i=1}^N\left((X_i-C_x)^2+(Y_i-C_y)^2-R\right)^2=0\\ &&\Rightarrow2\sum_{i=1}^N\left(-2(X_i-C_x)\right)\left((X_i-C_x)^2+(Y_i-C_y)^2-R\right)=0\\ &&\Rightarrow\sum_{i=1}^NX_i\left((X_i-C_x)^2+(Y_i-C_y)^2-R\right)-C_x\underbrace{\sum_{i=1}^N\left((X_i-C_x)^2+(Y_i-C_y)^2-R\right)}_{=0\ \because\ eq.(6)}=0\\ &&\Rightarrow\sum_{i=1}^NX_i\left((X_i-C_x)^2+(Y_i-C_y)^2-R\right)=0\\ &&\Rightarrow\sum_{i=1}^NX_i^3-2C_x\sum_{i=1}^NX_i^2+\sum_{i=1}^NX_iY_i^2-2C_y\sum_{i=1}^NX_iY_i+\left(C_x^2+C_y^2-R\right)\underbrace{\sum_{i=1}^NX_i}_{=0\ \because\ eq.(3)}=0\\ &&\Rightarrow\sum_{i=1}^NX_i^3-2C_x\sum_{i=1}^NX_i^2+\sum_{i=1}^NX_iY_i^2-2C_y\sum_{i=1}^NX_iY_i=0\\ &&\Rightarrow T_{30}-2C_xT_{20}+T_{12}-2C_yT_{11}=0\tag{8} \end{eqnarray}$

$E$ を $C_y$ で微分して、 $=0$ とおき、 $C_x$ で微分した結果 $(8)$ を用いると以下の式が成り立ちます。

$\begin{eqnarray} T_{03}-2C_yT_{02}+T_{21}-2C_xT_{11}=0\tag{9} \end{eqnarray}$

式 $(8),(9)$ を変形して、 $C_x,C_y$ を求めます。

$\begin{eqnarray} &&\left\{ \begin{array}{l} T_{30}-2C_xT_{20}+T_{12}-2C_yT_{11}=0\\ T_{03}-2C_yT_{02}+T_{21}-2C_xT_{11}=0 \end{array} \right.\\ &&\Rightarrow\begin{pmatrix} T_{20} & T_{11}\\ T_{11} & T_{02} \end{pmatrix} \begin{pmatrix} C_x\\ C_y \end{pmatrix} =\dfrac{1}{2} \begin{pmatrix} T_{30} + T_{12}\\ T_{03} + T_{21} \end{pmatrix}\\ &&\Rightarrow \begin{pmatrix} C_x\\ C_y \end{pmatrix} = \begin{pmatrix} T_{20} & T_{11}\\ T_{11} & T_{02} \end{pmatrix}^{-1} \dfrac{1}{2} \begin{pmatrix} T_{30} + T_{12}\\ T_{03} + T_{21} \end{pmatrix}\\ &&\Rightarrow \begin{pmatrix} C_x\\ C_y \end{pmatrix} = \dfrac{1}{2(T_{20}T_{02}-T_{11}^2)} \begin{pmatrix} T_{02} & -T_{11}\\ -T_{11} & T_{20} \end{pmatrix} \begin{pmatrix} T_{30} + T_{12}\\ T_{03} + T_{21} \end{pmatrix}\\ &&\Rightarrow \begin{pmatrix} C_x\\ C_y \end{pmatrix} = \dfrac{1}{2(T_{20}T_{02}-T_{11}^2)} \begin{pmatrix} T_{02}(T_{30} + T_{12})-T_{11}(T_{30} + T_{12})\\ -T_{11}(T_{03} + T_{21})+T_{20}(T_{03} + T_{21}) \end{pmatrix}\tag{10} \end{eqnarray}$

式 $(10)$ より、 $C_x,C_y$ が求まったので、これを式 $(7)$ に代入すると $R$ が求まります。
$c_x,c_y,r$ は式 $(4)$ と $R=r^2$ により以下のように求まります。

$\begin{eqnarray} \left\{ \begin{array}{l} c_x=C_x+\bar{x}\\ c_y=C_y+\bar{y}\\ r=\sqrt{R} \end{array} \right.\tag{11} \end{eqnarray}$

以上より、点群にフィットする円が定まりました。

円が定まらないときについての考察 - 普通の最小二乗円

式 $(10)$ の分母が $0$ の時は計算不能です。
$T_{20}T_{02}-T_{11}^2 = 0 \Rightarrow\left(\displaystyle\sum_{i=1}^NX_i^2\right)\left(\displaystyle\sum_{i=1}^NY_i^2\right)-\left(\displaystyle\sum_{i=1}^NX_iY_i\right)\left(\displaystyle\sum_{i=1}^NX_iY_i\right)=0$ となるのは、シュワルツの不等式の等号が成り立つときです。
それは、全ての点 $(X_i,Y_i)$ と原点を通る直線が存在する時です。この時、元の座標系の全ての点 $(x_i,y_i)$ を通る直線が存在しています。
よって、全ての点 $(x_i,y_i)$ を通る直線が存在しない時に、フィットする円は計算可能であることが分かります。

計算機で「解なし（点群にフィットする円を求められない）」を判定するときは許容誤差用定数 $\epsilon$ を用いて $|T_{20}T_{02}-T_{11}^2|<\epsilon$ とするのではなく、
点群の覆う矩形の対角線の長さ $l$ をかけて、 $|T_{20}T_{02}-T_{11}^2|< l \epsilon$ とするなどの工夫が必要かと思います。

ある1点を通る最小二乗円

「普通の最小二乗円」のアルゴリズムを使いまわすため、円が通る点を $(\bar{x},\bar{y})$ とおきます。
そうすると、式 $(3)$ が成り立たないので、 $(7)$ が成り立ちません。 $(6)$ は成り立ちます。
変換後の座標系では円が原点を通るので、以下が成り立ちます。

$\begin{eqnarray} R=C_x^2+C_y^2\tag{12} \end{eqnarray}$

$E$ を $R$ で微分して、 $=0$ とおきます。式 $(8)$ の計算を途中まで使って計算します。

$\begin{eqnarray} &&\frac{\partial E}{\partial C_x}=0\\ &&\Rightarrow\sum_{i=1}^NX_i^3-2C_x\sum_{i=1}^NX_i^2+\sum_{i=1}^NX_iY_i^2-2C_y\sum_{i=1}^NX_iY_i+\underbrace{\left(C_x^2+C_y^2-R\right)}_{=0\ \because\ eq.(12)}\sum_{i=1}^NX_i=0\\ &&\Rightarrow\sum_{i=1}^NX_i^3-2C_x\sum_{i=1}^NX_i^2+\sum_{i=1}^NX_iY_i^2-2C_y\sum_{i=1}^NX_iY_i=0\\ &&\Rightarrow T_{30}-2C_xT_{20}+T_{12}-2C_yT_{11}=0\tag{13} \end{eqnarray}$

式 $(13)$ と式 $(8)$ は全く同じ式なので、 $C_x,C_y$ は式 $(10)$ で求まり、これを式 $(12)$ に代入すると $R$ が求まります。
あとは式 $(11)$ で $c_x,c_y,r$ が求まります。

以上より、点群にフィットする「ある1点を通る円」が定まりました。

円が定まらないときについての考察 - ある1点を通る最小二乗円

前の議論同様、円が定らないのは全ての点 $(X_i,Y_i)$ と原点を通る直線が存在する時です。この時、元の座標系の全ての点 $(x_i,y_i)$ と $(\bar{x},\bar{y})$ を通る直線が存在しています。
よって、全ての点 $(x_i,y_i)$ と $(\bar{x},\bar{y})$ を通る直線が存在しない時に、フィットする円は計算可能であることが分かります。

ある2点を通る最小二乗円

「普通の最小二乗円」のアルゴリズムを使いまわしません。
円が通る2点を $(x_a,y_a),\ (x_b,y_b)$ とおき、2点の中間点を $(x_m,y_m)$ とおきます。
ただし、円が通る2点は等しくないものとします。 $(x_a,y_a)\not=(x_b,y_b)$

$\begin{eqnarray} x_m=\dfrac{x_a+x_b}{2},\ y_m=\dfrac{y_a+y_b}{2}.\tag{14} \end{eqnarray}$

$(x_a-x_b,y_a-y_b)$ と $(1,0)$ のなす角を $\theta$ とおくと、 $\sin\theta,\cos\theta$ は以下のようになります。

$\begin{eqnarray} \cos\theta=\dfrac{x_a-x_b}{\sqrt{(x_a-x_b)^2+(y_a-y_b)^2}},\ \sin\theta=\dfrac{y_a-y_b}{\sqrt{(x_a-x_b)^2+(ya_-y_b)^2}}\tag{15} \end{eqnarray}$

$(-x_m,-y_m)$ だけ平行移動した後、 $-\theta$ 回転した座標系の点群を $\{(X_i,Y_i)\}_{i=1}^N$ とします。

$\begin{eqnarray} \begin{pmatrix} X_i\\ Y_i \end{pmatrix} &=& \begin{pmatrix} \cos(-\theta) & -\sin(-\theta)\\ \sin(-\theta) & \cos(-\theta) \end{pmatrix} \begin{pmatrix} x_i-x_m\\ y_i-y_m \end{pmatrix}\\ &=& \begin{pmatrix} \cos\theta & \sin\theta\\ -\sin\theta & \cos\theta \end{pmatrix} \begin{pmatrix} x_i-x_m\\ y_i-y_m\\ \end{pmatrix}\\ &=& \begin{pmatrix} (x_i-x_m)\cos\theta+(y_i-y_m)\sin\theta\\ -(x_i-x_m)\sin\theta+(y_i-y_m)\cos\theta\tag{16} \end{pmatrix} \end{eqnarray}$

式 $(16)$ の逆変換は $\theta$ 回転した後、 $(x_m,y_m)$ だけ平行移動したものです。

$\begin{eqnarray} \begin{pmatrix} x_i\\ y_i \end{pmatrix} &=& \begin{pmatrix} \cos\theta & -\sin\theta\\ \sin\theta & \cos\theta \end{pmatrix} \begin{pmatrix} X_i\\ Y_i \end{pmatrix} + \begin{pmatrix} x_m\\ y_m \end{pmatrix}\\ &=& \begin{pmatrix} X_i\cos\theta -Y_i\sin\theta\\ X_i\sin\theta + Y_i\cos\theta \end{pmatrix} + \begin{pmatrix} x_m\\ y_m\\ \end{pmatrix}\\ &=& \begin{pmatrix} X_i\cos\theta -Y_i\sin\theta+x_m\\ X_i\sin\theta + Y_i\cos\theta+y_m\tag{17} \end{pmatrix} \end{eqnarray}$

式 $(16),(17)$ の写像を $f,f^{-1}$ とおきます。

$\begin{eqnarray} &&f:{\mathbb R}^2\rightarrow{\mathbb R}^2,\ \begin{pmatrix}x\\y\end{pmatrix} \mapsto \begin{pmatrix} (x-x_m)\cos\theta+(y-y_m)\sin\theta\\ -(x-x_m)\sin\theta+(y-y_m)\cos\theta \end{pmatrix}\tag{18}\\ &&f^{-1}:{\mathbb R}^2\rightarrow{\mathbb R}^2,\ \begin{pmatrix}x\\y\end{pmatrix} \mapsto \begin{pmatrix} x\cos\theta - y\sin\theta+x_m\\ x\sin\theta + y\cos\theta+y_m \end{pmatrix}\tag{19} \end{eqnarray}$

$(x_a,y_a)$ は変換後、 $x$ 軸上に存在するので、変換後の座標を $(X_a,0)$ とおくと、以下のようになります。

$\begin{eqnarray} \begin{pmatrix}X_a\\0\end{pmatrix}=f \begin{pmatrix}x_a\\y_a\end{pmatrix} = \begin{pmatrix} (x_a-x_m)\cos\theta+(x_a-x_m)\sin\theta\\ 0 \end{pmatrix}\tag{20} \end{eqnarray}$

変換前、変換後の座標系における円の中心座標をそれぞれ $(c_x,c_y),\ (C_x,C_y)$ とします。
変換後の円の中心は、 $y$ 軸上に存在するので以下の式が成り立ちます。

$\begin{eqnarray} C_x=0\tag{21} \end{eqnarray}$

円の半径を $r$ とし、 $R=r^2$ とおくと、 $R,X_a,C_y$ に以下の関係が成り立ちます。

$\begin{eqnarray} R=X_a^2+C_y^2\tag{22} \end{eqnarray}$

目的関数 $E$ を以下のように定めます。

$\begin{eqnarray} E&=&\sum_{i=1}^N\big((X_i-\underbrace{C_x}_{eq.(21)})^2+(Y_i-C_y)^2-\underbrace{R}_{eq.(22)}\big)^2\\ &=&\sum_{i=1}^N\left(X_i^2+Y_i^2-2Y_iC_y+C_y^2-X_a^2-C_y^2\right)^2\\ &=&\sum_{i=1}^N\left(-2Y_iC_y+X_i^2+Y_i^2-X_a^2\right)^2\tag{23} \end{eqnarray}$

式を簡単にするために $T_{mn}=\displaystyle\sum_{i=1}^NX_i^mY_i^n\ (m,n=0,1,2,3)$ とおきます。

$E$ を $C_y$ で微分して、 $=0$ とおきます。

$\begin{eqnarray} &&\frac{\partial E}{\partial C_y}=0\\ &&\Rightarrow\sum_{i=1}^N-4Y_i\left(-2Y_iC_y+X_i^2+Y_i^2-X_a^2\right)=0\\ &&\Rightarrow-2C_y\sum_{i=1}^NY_i^2+\sum_{i=1}^NX_i^2Y_i+\sum_{i=1}^NY_i^3-X_a^2\sum_{i=1}^NY_i=0\\ &&\Rightarrow-2T_{02}C_y+T_{21}+T_{03}-X_a^2T_{01}=0\\ &&\Rightarrow C_y=\dfrac{T_{21}+T_{03}-X_a^2T_{01}}{2T_{02}}\tag{24} \end{eqnarray}$

式 $(24)$ より、 $C_y$ が求まったので、これを式 $(22)$ に代入すると $R$ が求まります。
$r$ は $R=r^2$ より、 $R$ の平方根を取ることで求まります。

$\begin{eqnarray} r=\sqrt{R}\tag{25} \end{eqnarray}$

$c_x,c_y$ は $(0,C_y)$ を逆変換することで求まります。

$\begin{eqnarray} \begin{pmatrix}c_x\\ c_y\end{pmatrix} =f^{-1}\begin{pmatrix}0\\ C_y\end{pmatrix} =\begin{pmatrix} - C_y\sin\theta+x_m\\ C_y\cos\theta+y_m \end{pmatrix}\tag{26} \end{eqnarray}$

以上より、点群にフィットする「ある2点を通る円」が定まりました。

円が定まらないときについての考察 - ある2点を通る最小二乗円

式 $(24)$ の分母が $0$ の時は計算不能です。
$T_{02} = 0 \Rightarrow\displaystyle\sum_{i=1}^NY_i^2=0$ となるのは、
全ての点 $(X_i,Y_i)$ がX軸上に存在する時です。この時、元の座標系の全ての点 $(x_i,y_i)$ と $(x_a,y_a)$ と $(x_b,y_b)$ を通る直線が存在しています。
よって、全ての点 $(x_i,y_i)$ と $(x_a,y_a)$ と $(x_b,y_b)$ を通る直線が存在しない時に、フィットする円は計算可能であることが分かります。

TODO

図を挿入する
普通の最小二乗円：変換の図、E用の円の図、円が計算できない場合の図(変換後の座標系⇒変換前の座標系)
1点を通る最小二乗円：同様
2点を通る最小二乗円：変換の図を2段階で書く

参考リンク

円の最小二乗法を解説
 二点を固定した円の最小二乗法による近似
 Cauchy–Schwarz 不等式の等号成立条件

機械学習基礎理論独習

誤りがあればご指摘いただけると幸いです。数式が整うまで少し時間かかります。リンクフリーです。