Basics of Paleomagnetism and Rock Magnetism

Circle fitting: Improved method of Taubin (1991)

Taubin method is one of those proposed to improve the Kasa method (see, Chernov and Lesort 2005). Hereafter, a circle equation is denoted as, \[ A(x^2 + y^2) + Bx + Cy + D = 0 \] Transforming this to the ordinary form, \[ \left(x + \frac{B}{2A}\right)^2 + \left(y + \frac{C}{2A}\right)^2 = \frac{B^2 + C^2 - 4AD}{4A^2} \] Hence, \begin{eqnarray*} a & = & -\frac{B}{2A}, \\ b & = & -\frac{C}{2A}, \\ R^2 & = & \frac{B^2 + C^2 - 4AD}{4A^2}. \end{eqnarray*} As \(B^2+C^2-4AD>0\) and the parameters can be scaled by a scalar, the next constraint can be imposed, \begin{equation} B^2 + C^2 - 4AD = 1. \label{eq01} \end{equation} As stated in the previous page, \({\cal F} = \frac{1}{4R^2}{\cal F_K}\) should improve the Kasa method, and Taubin (1991) adopted a mean of \(r_i^2\) for \(R^2\). Hence, Taubin method minimizes, \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} where \(z_i=x_i^2 + y_i^2\) and \(\bar x = \frac{1}{n}\sum_{i=1}^n x_i\), etc. Taking \eqref{eq01} into consideration, minimizing \eqref{eq02} is equivalent to minimizing \begin{equation} {\cal F_1} = \sum_{i=1}^n (Az_i + Bx_i + Cy_i + D)^2 \label{eq03} \end{equation} on the constraint of \begin{equation} 4A^2M_z + 4ABM_x + 4ACM_y + B^2n + C^2n =1, \label{eq04} \end{equation} where \(M_x = \sum_{i=1}^n x_i\), etc. Using matrix notation, \eqref{eq03} and \eqref{eq04} are \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} and \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} Minimizing \eqref{eq05} on the constraint \eqref{eq06} is a typical inverse problem that minimizes \[ {\cal F_*} = {^t{\bf \alpha}}{\bf M}{\bf \alpha} - \eta({^t{\bf \alpha}}{\bf N}{\bf \alpha} - 1), \] where \(\eta\) is a Lagrange multiplier. \(\partial{\cal F_*}/\partial{\bf \alpha}=0\) gives, \[ {\bf M}{\bf \alpha} - \eta{\bf N}{\bf \alpha} = 0. \] This is an eigenvalue problem and \(\eta\) is obtained by \[ \det({\bf M} - \eta{\bf N}) = 0. \] Corresponding \({^t{\bf \alpha}}=(A,B,C,D)\) is determined by solving \[ ({\bf M} - \eta{\bf N}){\bf \alpha} = 0. \]

As stated in the previous page, the algebraic methods are not the most accurate. Nevertheless, they are practical due to their fast calculation speed. If the accuracy is the most important, you should use the geometric methods, in which the simpler algebraic methods are still useful to obtain an initial guess. Excellent programs of both algebraic and geometric methods by Prof. Nikolai Chernov are posted on the web site (Chernov, 2012).

Program of circle fitting by Taubin method

Circle fitting programs in C language which use Taubin' algebraic method are available. The programming code is adapted from the C++ routines of Chernov (2012). Program "pcirc" works on Linux system in which Generic Mapping Tools 6 (or GMT 5) of Wessel et al. (2019) is installed. Program "pcalc" is for paleomagnetic analysis and includes circle fitting. "pcalc" does not need GMT 6 (GMT 5), and hence should work on terminals of most Linux systems.

→ See this page.

References:

↑ Top of this page   ← Previous page