円近似: 改良された Taubin (1991) の方法
Taubin 法は Kasa 法を改良するべく提案された多くの方法の一つです (Chernov and Lesort 2005)．ここからは円の方程式を次式で表します． \[ A(x^2 + y^2) + Bx + Cy + D = 0 \] 通常の形に直すと， \[ \left(x + \frac{B}{2A}\right)^2 + \left(y + \frac{C}{2A}\right)^2 = \frac{B^2 + C^2 - 4AD}{4A^2} \] 故に， \begin{eqnarray*} a & = & -\frac{B}{2A}, \\ b & = & -\frac{C}{2A}, \\ R^2 & = & \frac{B^2 + C^2 - 4AD}{4A^2}. \end{eqnarray*} \(B^2+C^2-4AD>0\) であり，パラメーターはスカラー倍が許されるので，次の条件を課します． \begin{equation} B^2 + C^2 - 4AD = 1. \label{eq01} \end{equation} 前ページで述べたように， \({\cal F} = \frac{1}{4R^2}{\cal F_K}\) は Kasa 法を改良するはずで， Taubin (1991) は \(R^2\) として \(r_i^2\) の平均を採用しました．故に， Taubin 法は次式を最小にします． \begin{eqnarray} {\cal F_T} & = & \frac{\sum_{i=1}^n [(x_i-a)^2 + (y_i-b)^2 -R^2]^2} {4n^{-1} \sum_{i=1}^n [(x_i-a)^2 + (y_i-b)^2]} \nonumber \\ & = & \frac{\sum_{i=1}^n (Az_i + Bx_i + Cy_i + D)^2} {4A^2\bar z + 4AB\bar x + 4AC\bar y + B^2 + C^2} \label{eq02} \end{eqnarray} 但し， \(z_i=x_i^2 + y_i^2\) で \(\bar x = \frac{1}{n}\sum_{i=1}^n x_i\) 等です． \eqref{eq01} を考慮すると， \eqref{eq02} を最小にすることは次式， \begin{equation} {\cal F_1} = \sum_{i=1}^n (Az_i + Bx_i + Cy_i + D)^2 \label{eq03} \end{equation} を次の条件のもとで最小にすることと同じです． \begin{equation} 4A^2M_z + 4ABM_x + 4ACM_y + B^2n + C^2n =1, \label{eq04} \end{equation} 但し， \(M_x = \sum_{i=1}^n x_i\) 等です. \eqref{eq03} と \eqref{eq04} は行列表記を用いると，最小にすべき次式と， \begin{eqnarray} {\cal F_1} & = & {^t{\bf \alpha}} {\bf M} {\bf \alpha} \nonumber \\ & = & \left(\begin{array}{cccc} A & B & C & D \end{array}\right) \left(\begin{array}{cccc} M_{zz} & M_{xz} & M_{yz} & M_z \\ M_{xz} & M_{xx} & M_{xy} & M_x \\ M_{yz} & M_{xy} & M_{yy} & M_y \\ M_z & M_x & M_y & n \end{array}\right) \left(\begin{array}{c} A \\ B \\ C \\ D \end{array}\right) \label{eq05} \end{eqnarray} 次の条件式， \begin{equation} {^t{\bf \alpha}} {\bf N} {\bf \alpha} = \left(\begin{array}{cccc} A & B & C & D \end{array}\right) \left(\begin{array}{cccc} 4M_z & 2M_x & 2M_y & 0 \\ 2M_x & n & 0 & 0 \\ 2M_y & 0 & n & 0 \\ 0 & 0 & 0 & 0 \end{array}\right) \left(\begin{array}{c} A \\ B \\ C \\ D \end{array}\right) = 1 \label{eq06} \end{equation} となります．条件式 \eqref{eq06} のもとに式 \eqref{eq05} を最小にすることは典型的な逆問題で， \(\eta\) をラグランジュ乗数として次式を最小にします． \[ {\cal F_*} = {^t{\bf \alpha}}{\bf M}{\bf \alpha} - \eta({^t{\bf \alpha}}{\bf N}{\bf \alpha} - 1). \] \(\partial{\cal F_*}/\partial{\bf \alpha}=0\) より， \[ {\bf M}{\bf \alpha} - \eta{\bf N}{\bf \alpha} = 0. \] を得ます．これは固有値問題で \(\eta\) は次式から得られます． \[ \det({\bf M} - \eta{\bf N}) = 0. \] 対応する \({^t{\bf \alpha}}=(A,B,C,D)\) は次式を解くことで得られます． \[ ({\bf M} - \eta{\bf N}){\bf \alpha} = 0. \]
前ページで述べたように，代数学的方法は最も正確というわけではありません．しかしこの方法は速い計算速度のために実践的です．もし正確さが最も大切な場合は，幾何学的方法を使用するべきですが，その場合にも，より単純な代数学的方法は初期値を得るために役に立ちます． Nikolai Chernov 教授は両者の方法の優れたプログラムをウェブ・サイトに載せています (Chernov, 2012)．
Taubin 法による円近似のプログラム
Taubin の代数学的方法を用いた C 言語によるプログラムが利用可能です．プログラム・コードは Chernov (2012) の C++ ルーチンを C 言語用に変更しました．プログラム "pcirc" は Wessel et al. (2019) の Generic Mapping Tools 6 (または GMT 5) がインストールされた Linux システム上で動きます．プログラム "pcalc" は古地磁気データの計算を目的としていますが，円近似も含まれています． "pcalc" は GMT 6 (GMT 5) を必要としないため，多くの Linux システムの端末上で動くと思われます．
参考文献:
- Chernov, N., and C. Lesort, Least squares fitting of circles, J. Math. Imaging Vision, 23, 239-251, 2005.
- Chernov, N, Fitting ellipses, circles, and lines by least squares, 2012. (URL: →http://people.cas.uab.edu/~mosya/cl/ index.html)
- Taubin, G., Estimation of planar curves, surfaces and nonplanar space curves defined by implicit equations, with applications to edge and range image segmentation, IEEE Trans. Pattern Analysis Machine Intelligence, 13, 1115-1138, 1991.
- Wessel, P., Luis, J.F., Uieda, L., Scharroo, R., Wobbe, F., Smith, W.H.F., and Tian, D., The Generic Mapping Tools version 6, Geochemistry, Geophysics, Geosystems, 20, 5556–5564, 2019. (doi: 10.1029/2019GC008515; Manuals: https://docs.generic-mapping-tools.org/6.0/)