\documentclass[a4paper,11pt,leqno]{article}
\usepackage{color}
\usepackage{graphics}
\usepackage{amsmath}
\usepackage{times}
\usepackage{mathptm}
\usepackage{url}
\usepackage{float}
\floatstyle{boxed}
\newfloat{algorithm}{p}{loa}\floatname{algorithm}{Algorithm}
\begin{document}
\title{Finding the circle that best fits a set of points}
\author{L. \textsc{Maisonobe}\thanks{luc@spaceroots.org}}
\date{October 25${}^\text{th}$ 2007}\maketitle\tableofcontents
\clearpage

\section{Introduction}
We consider the following 2D problem: given a set of $n$ points
$P_i(x_i, y_i)$ ($i=1\ldots n$) find the center $C(x_c, y_c)$ and the
radius $r$ of the circle that pass closest to all the points.

\begin{figure}[htbp]
\begin{center}\begin{minipage}{93.75mm}
\caption{\label{fig:problem-statement}closest circle problem}
\setlength{\unitlength}{0.075mm}\begin{picture}(1250,1400)
  \textcolor[gray]{0.5}{
    \put(0,300){\vector(1,0){1250}}\put(1240,320){\mbox{$x$}}
    \put(60,0){\vector(0,1){1200}}\put(80,1190){\mbox{$y$}}
  }%

  \textcolor[rgb]{0.13,0.55,0.13}{
    \put(360,980){\circle{20}}\put(380,950){\mbox{$P_1$}}
    \put(560,240){\circle{20}}\put(490,190){\mbox{$P_2$}}
    \put(1160,100){\circle{20}}\put(1090,130){\mbox{$P_3$}}
    \put(410,450){\circle{20}}\put(430,450){\mbox{$P_4$}}
    \put(510,1270){\circle{20}}\put(420,1270){\mbox{$P_5$}}}%

  \textcolor[rgb]{0.15,0.13,0.73}{
    \put(1022.5,782.2){\circle*{15}}\put(1050,770){\mbox{$C$}}
    \put(709.1,938.9){\vector(-2,1){313.4}}
    \put(709.1,938.9){\vector(2,-1){313.4}}
    \put(712,945){\mbox{$r$}}
    \qbezier[0](600.7,1342)(551.4,1304.8)(509.4,1259.7)
    \qbezier[0](509.4,1259.7)(467.4,1214.5)(433.9,1162.7)
    \qbezier[0](433.9,1162.7)(400.4,1110.9)(376.4,1054)
    \qbezier[0](376.4,1054)(352.5,997.1)(338.9,936.9)
    \qbezier[0](338.9,936.9)(325.3,876.7)(322.4,815.1)
    \qbezier[0](322.4,815.1)(319.5,753.5)(327.4,692.3)
    \qbezier[0](327.4,692.3)(335.3,631.1)(353.8,572.2)
    \qbezier[0](353.8,572.2)(372.3,513.4)(400.8,458.6)
    \qbezier[0](400.8,458.6)(429.3,403.9)(466.9,355)
    \qbezier[0](466.9,355)(504.5,306.1)(550.1,264.5)
    \qbezier[0](550.1,264.5)(595.7,222.9)(647.8,189.9)
    \qbezier[0](647.8,189.9)(700,156.9)(757.1,133.6)
    \qbezier[0](757.1,133.6)(814.2,110.2)(874.5,97.2)
    \qbezier[0](874.5,97.2)(934.8,84.2)(996.5,81.9)
    \qbezier[0](996.5,81.9)(1058.2,79.6)(1119.3,88.1)
    \qbezier[0](1119.3,88.1)(1180.4,96.6)(1239.1,115.7)}%

\end{picture}\end{minipage}\end{center}
\end{figure}

For the latest version of this document, please check the downloads
page of the spaceroots site at
\url{http://www.spaceroots.org/downloads.html}. The paper can be
browsed on-line or retrieved as a PDF, compressed PostScript or LaTeX
source file. An implementation in the Java language of the algorithms
presented here is also available in source form as a stand-alone file
at \url{http://www.spaceroots.org/documents/circle/CircleFitter.java}.


\section{Solving the problem}

\subsection{Principles}
The underlying principles of the proposed fitting method are to first
compute an initial guess by averaging all the circles that can be
built using all triplets of non-aligned points, and then to
iteratively reduce the distance between the circle and the complete
set of points using a minimization method.

\subsection{Initialization}
For any given triplet of non-aligned points, there is a single circle
passing through all three points: the triangle circumcircle. We will
use this property to build an initial circle.

Let $P_i(x_i, y_i)$, $P_j(x_j, y_j)$ and $P_k(x_k, y_k)$ be three
points. The triangle circumcenter is at the intersection of the
perpendicular bisectors of the triangle sides. Equation
(\ref{eqn:bisector-1-2}) holds as the circumcenter belongs to the
perpendicular bisector of side $(P_i,P_j)$ and equation
(\ref{eqn:bisector-2-3}) holds it also belongs to the perpendicular
bisector of side $(P_j,P_k)$:

\begin{gather}
\left\{\begin{aligned}
x_{C_{i,j,k}} &= \frac{(x_i+x_j)+\alpha_{i,j}(y_j-y_i)}{2}\\
y_{C_{i,j,k}} &= \frac{(y_i+y_j)-\alpha_{i,j}(x_j-x_i)}{2}
\end{aligned}\right.\label{eqn:bisector-1-2}\\
\left\{\begin{aligned}
x_{C_{i,j,k}} &= \frac{(x_j+x_k)+\alpha_{j,k}(y_k-y_j)}{2}\\
y_{C_{i,j,k}} &= \frac{(y_j+y_k)-\alpha_{j,k}(x_k-x_j)}{2}
\end{aligned}\right.\label{eqn:bisector-2-3}
\end{gather}

Solving this set of linear equations is straightforward:
\begin{gather}
\left\{\begin{aligned}
\alpha_{i,j} &= \frac{(x_k-x_i)(x_k-x_j)+(y_k-y_i)(y_k-y_j)}{\Delta}\\
\alpha_{j,k} &= \frac{(x_k-x_i)(x_j-x_i)+(y_k-y_i)(y_j-y_i)}{\Delta}
\end{aligned}\right.\nonumber\\
\text{where}\quad
\Delta=(x_k-x_j)(y_j-y_i)-(x_j-x_i)(y_k-y_j)\label{eqn:determinant}
\end{gather}

Hence the circumcenter coordinates are:
\begin{equation}\label{eqn:circumcenter}
\left\{\begin{aligned}
x_{C_{i,j,k}} &= \phantom{-}\frac{(y_k-y_j)(x_i^2+y_i^2)
                                 +(y_i-y_k)(x_j^2+y_j^2)
                                 +(y_j-y_i)(x_k^2+y_k^2)}{2\Delta}\\
y_{C_{i,j,k}} &=          - \frac{(x_k-x_j)(x_i^2+y_i^2)
                                 +(x_i-x_k)(x_j^2+y_j^2)
                                 +(x_j-x_i)(x_k^2+y_k^2)}{2\Delta}
\end{aligned}\right.\\
\end{equation}

These coordinates can be computed as long as the determinant $\Delta$
is non-null, which means the three points are not aligned.

Since the method is intended to be general, we do not make any
assumptions on the points. Some of them may have small (or even large)
errors, some may be aligned with other ones, some may appear several
times in the set. We cannot reliably select one triplet among the
other ones. We build the initial guess of the fitting circle by
averaging the coordinates of the circumcenters of all triplets of non
aligned points.

This entire initialization process is implemented in algorithm
\ref{alg:initialization} (page \pageref{alg:initialization}).

\subsection{Issues}
The initialization process described above may be inadequate in some
cases.

If the number of points is very large, the number of triplets will be
overwhelmingly large as it grows in $n^3$. In this case, only a subset
of the complete triplets set must be used. Selecting the subset
implies knowing the distribution of the points, in order to avoid
using only points in the same area for initialization.

If the points sample contains really far outliers in addition to
regular points, the initial center may be offset by a large
amount. This is due to the fact the mean is not a robust
statistic. This can be circumvented by either using a median rather
than a mean or by dropping outliers (either before or after
circumcenter determination).

If the points distribution is known to cover exactly one circle with a
consistent step between the points and no outliers, a much simpler
initialization process is simply to use the mean of the x and y
coordinates for all points.

\subsection{Iterative improvement}
\subsubsection{Problem statement}
Once we have an initial guess for the circle center, we want to
improve it for some definition of \emph{best fit} against the points
set. Using a least squares estimator based on the euclidean distance
between the points and the circle is a common choice.

We try to minimize the cost function $J$:
\begin{equation*}
J=\sum_{i=1}^{n}(d_i-r)^2
\quad\text{where}\quad
d_i=\sqrt{(x_i-x)^2+(y_i-y)^2}
\end{equation*}
$d_i$ is the euclidean distance between the point $P_i(x_i,y_i)$ and
the circle center $C(x, y)$ and $r$ is the circle radius.

For a given circle center, we compute the optimum circle radius $\hat{r}$ by
solving:
\begin{equation}\label{eqn:circle-radius}
\left.\frac{\partial J}{\partial r}\right|_{r=\hat{r}} = 0
\Rightarrow \hat{r} = \frac{\sum_{i=1}^{n}d_i}{n}
\end{equation}

This means that both $d_i$ and $\hat{r}$ can be considered as
functions of the circle center coordinates $x$ and $y$ which from now
on will be the free parameters of our model.

Using the choice
$r=\hat{r}$, the cost function can be rewritten:
\begin{equation}
J=\sum_{i=1}^{n}(d_i-\hat{r})^2
 =\sum_{i=1}^{n}d_i^2 - \frac{\left(\sum_{i=1}^{n}d_i\right)^2}{n}
 \label{eqn:cost-value}
\end{equation}
from this expression we can compute the gradient of the cost function
with respect to the free parameters:
\begin{equation}
\nabla J\left\{\begin{aligned}
  \frac{\partial J}{\partial x}
  &= 2\sum_{i=0}^{n}\frac{(x-x_i)(d_i-\hat{r})}{d_i}\\
  \frac{\partial J}{\partial y}
  &= 2\sum_{i=0}^{n}\frac{(y-y_i)(d_i-\hat{r})}{d_i}
\end{aligned}\right.\label{eqn:cost-gradient}
\end{equation}

\subsubsection{Direct solving}
Equation (\ref{eqn:cost-gradient}) can easily be solved analytically. However
experience shows that even for small errors with the initial guess,
this leads to a solution $C(x, y)$ that is very far from the validity
domain of the partial derivatives values.

This sensitivity to the initial guess is a common issue with
\textsc{Gauss-Newton} based least squares solvers. Iterating as in a
fixed point method does not help either, iterations diverge once we
are too far. So we have to refrain from using this easy but wrong
solution and use a more robust approach.

\subsubsection{Levenberg-Marquardt}
The \textsc{Levenberg-Marquardt} method is a combination of the
\textsc{Gauss-Newton} and steepest descent methods. It benefits from
the strength of both methods and is both robust even for starting
points far from the solution and efficient near the solution. It is a
sure bet and is highly recommended.

This method is however quite complex to implement, so it should be
used only if a validated implementation of the method is easily
available. The MINPACK
implementation\footnote{\url{http://www.netlib.org/minpack/lmder.f}}
in fortran is widely used. Our own
Mantissa\footnote{\url{http://www.spaceroots.org/software/mantissa/index.html}}
Java library also provides a Java implementation based on the MINPACK
one as of version~6.0. Many other implementations are available in
various languages.

For cost function $J=\sum_{i=1}^{n}(d_i-\hat{r})^2$, the
\textsc{Levenberg-Marquardt} method needs the jacobian matrix of
the deviations $d_i-\hat{r}$. The two elements of row $i$ of this matrix is given by
equation~ (\ref{eqn:deviation-derivative}):
\begin{align}
\nabla (d_i-\hat{r})
&= \left\{\begin{gathered}
          \frac{\partial\left(d_i-\frac{1}{n}\sum_{k=0}^{n}d_k\right)}{\partial x}\\
         \frac{\partial\left(d_i-\frac{1}{n}\sum_{k=0}^{n}d_k\right)}{\partial y}
         \end{gathered}\right.\nonumber\\
&= \left\{\begin{gathered}
          \frac{x-x_i}{d_i}-\frac{1}{n}\sum_{k=0}^{n}\frac{x-x_k}{d_k}\\
          \frac{y-y_i}{d_i}-\frac{1}{n}\sum_{k=0}^{n}\frac{y-y_k}{d_k}
          \end{gathered}\right.\label{eqn:deviation-derivative}
\end{align}

This solver is far more efficient than the simple one we will present
in the following section (by several orders of magnitude). The drawbacks
are the implementation complexity.

\subsubsection{Conjugate gradient}
If no \textsc{Levenberg-Marquardt} based solver are available, we can
still find the values of the free parameters that minimize the cost
function $J$ using the simpler conjugate gradient method with
\textsc{Polak} and \textsc{Ribi\`ere} factor. This will be less
efficient than using \textsc{Levenberg-Marquardt} especially when the
points are far from a real circle, but is very simple to do.

The conjugate gradient method is also an iterative one, using two
embedded loops. The outer loop computes a new search direction by
combining the cost function gradient at the current step and the
previous search direction. The inner loop roughly minimizes the cost
function in the search direction. This method is very easy to
implement and allows to have a self-contained algorithm.

The rough minimization of the cost function along the search direction
can be done by performing a few steps of a \textsc{Newton} solver on the
derivative of the cost function $J(x + \lambda u, y + \lambda v)$
with respect to the step $\lambda$.

The following equations are the exact derivatives used for this rough
inner minimization (we use $J_{u,v}(\lambda)=J(x + \lambda u, y +
\lambda v)$, ${d_i}_{u,v}(\lambda)={d_i}(x + \lambda u, y +
\lambda v)$ and $\hat{r}_{u,v}(\lambda)=\hat{r}(x + \lambda u, y +
\lambda v)$ as a way to simplify the equation):

\begin{equation*}
\left\{\begin{aligned}
J_{u,v}(\lambda)
&= \sum_{i=0}^{n}{d_i}_{u,v}^2(\lambda)
 - \frac{1}{n}\left(\sum_{i=0}^{n}{d_i}_{u,v}(\lambda)\right)^2\\
\frac{dJ_{u,v}(\lambda)}{d\lambda}
&= 2\sum_{i=0}^{n}
      \frac{\big[(x+\lambda u-x_i)u+(y+\lambda v-y_i)v\big]
            \big[{d_i}_{u,v}(\lambda)-\hat{r}_{u,v}(\lambda)\big]}
           {{d_i}_{u,v}(\lambda)}\\
\frac{d^2J_{u,v}(\lambda)}{d\lambda^2}
&= \begin{aligned}[t]
  &2 (u^2+v^2)\sum_{i=0}^{n}\frac{{d_i}_{u,v}(\lambda)-\hat{r}_{u,v}(\lambda)}
                                 {{d_i}_{u,v}(\lambda)}\\
  \hspace{-1em}-&\frac{2}{n}\left(
    \sum_{i=0}^{n}\frac{(x+\lambda u-x_i)u+(y+\lambda v-y_i)v}
                       {{d_i}_{u,v}(\lambda)}\right)^2\\
  \hspace{-1em}+&2\hat{r}_{u,v}(\lambda)
    \sum_{i=0}^{n}\frac{\big((x+\lambda u-x_i)u+(y+\lambda v-y_i)v\big)^2}
                        {{d_i}_{u,v}^3(\lambda)}
\end{aligned}
\end{aligned}\right.
\end{equation*}
hence
\begin{equation}\label{eqn:first-derivative}
\left.\frac{dJ_{u,v}(\lambda)}{d\lambda}\right|_{\lambda=0}
= 2\sum_{i=0}^{n}
   \frac{\big[(x-x_i)u+(y-y_i)v\big]\big[d_i-\hat{r}\big]}{d_i}
\end{equation}
and
\begin{equation}\label{eqn:second-derivative}\begin{split}
\left.\frac{d^2J_{u,v}(\lambda)}{d\lambda^2}\right|_{\lambda=0}
  =&\phantom{-}2 (u^2+v^2)\sum_{i=0}^{n}\frac{d_i-\hat{r}}{d_i}\\
  \hspace{-1em}&-\frac{2}{n}\left(
    \sum_{i=0}^{n}\frac{(x-x_i)u+(y-y_i)v}{d_i}\right)^2\\
  \hspace{-1em}&+2\hat{r}
    \sum_{i=0}^{n}\frac{\big((x-x_i)u+(y-y_i)v\big)^2}{d_i^3}
\end{split}\end{equation}

The preceding equations are used in a traditional conjugate gradient
method, assuming $x$ and $y$ have been initialized, as shown in
algorithm \ref{alg:conjugate-gradient} (page
\pageref{alg:conjugate-gradient}).

The parameters of the algorithms are the maximal number of iterations
$i_\text{max}$, and the two convergence thresholds
$\varepsilon_{inner}$ and $\varepsilon_{outer}$. Their values depend
on the problem at hand (number of points, expected error at each
points, required accuracy). Finding the values for a given problem is
a matter of trial and errors. As always with the conjugate gradient
method, there is no real need to have an accurate minimum along the
search direction, hence $\varepsilon_{inner}$ may be quite large (I
used values between 0.02 and 0.1). It is even possible to remove
completely the convergence check and replace it with an arbitrary
small number of \textsc{Newton} steps.

\section{Numerical example}
This section shows an example application of the algorithms described
in this note. It is based on a simple case with only five points
and an order of magnitude of about 1\% of the circle radius for the
errors in the points locations.

We start with the following set of five points:

\begin{center}\begin{tabular}{|c|r|r|}
\hline
point & \multicolumn{1}{c|}{$x$} & \multicolumn{1}{c|}{$y$}\\
\hline
$P_1$ & 30  &  68\\
$P_2$ & 50  &  -6\\
$P_3$ & 110 & -20\\
$P_4$ & 35  &  15\\
$P_5$ & 45  &  97\\
\hline
\end{tabular}\end{center}

Since there are five points in the set, there are ten different
triplets, none of them have aligned points:
\begin{center}\begin{tabular}{|c|r|r||c|r|r|}
\hline
circumcenter & \multicolumn{1}{c|}{$x$} & \multicolumn{1}{c||}{$y$}&
circumcenter & \multicolumn{1}{c|}{$x$} & \multicolumn{1}{c|}{$y$}\\
\hline
$C_{1,2,3}$ &  93.650 & 45.500 & $C_{1,4,5}$ & 103.768 & 48.223 \\
$C_{1,2,4}$ & 103.704 & 48.217 & $C_{2,3,4}$ &  92.400 & 40.142 \\
$C_{1,2,5}$ & 103.753 & 48.230 & $C_{2,3,5}$ &  94.178 & 47.765 \\
$C_{1,3,4}$ &  95.821 & 47.473 & $C_{2,4,5}$ & 103.720 & 48.229 \\
$C_{1,3,5}$ &  99.228 & 50.571 & $C_{3,4,5}$ &  96.580 & 49.100 \\
\hline
\end{tabular}\end{center}

we observe that with this data set, four circumcenters ($C_{1,2,4}$,
$C_{1,2,5}$, $C_{1,4,5}$ and $C_{2,4,5}$) are very close to each
other. This is only a coincidence. This explains why few points show
up in figure \ref{fig:initialization}: the four rightmost points are
almost superimposed.

\begin{figure}[htbp]
\begin{center}\begin{minipage}{62.5mm}
\caption{\label{fig:initialization}initialization}
\setlength{\unitlength}{0.05mm}\begin{picture}(1250,1400)
  \textcolor[gray]{0.5}{
    \put(0,300){\vector(1,0){1250}}\put(1240,320){\mbox{$x$}}
    \put(60,0){\vector(0,1){1200}}\put(80,1190){\mbox{$y$}}
  }%
  \textcolor[rgb]{0.13,0.55,0.13}{
    \put(360,980){\circle{20}}\put(380,950){\mbox{$P_1$}}
    \put(560,240){\circle{20}}\put(490,190){\mbox{$P_2$}}
    \put(1160,100){\circle{20}}\put(1090,130){\mbox{$P_3$}}
    \put(410,450){\circle{20}}\put(430,450){\mbox{$P_4$}}
    \put(510,1270){\circle{20}}\put(420,1270){\mbox{$P_5$}}
    \put( 996.5,755.0){\circle{15}}
    \put(1097.0,782.1){\circle{15}}
    \put(1097.5,782.3){\circle{15}}
    \put(1018.2,774.7){\circle{15}}
    \put(1052.2,805.7){\circle{15}}
    \put(1097.6,782.2){\circle{15}}
    \put( 984.0,701.4){\circle{15}}
    \put(1001.7,777.6){\circle{15}}
    \put(1097.2,782.2){\circle{15}}
    \put(1025.8,791.0){\circle{15}}
    \put(1046.8,773.4){\circle{30}}
    \put(1022.4,782.2){\circle*{20}}
    \put(750,1070){\circle{15}}\put(800,1070){\mbox{circumcenters}}
    \put(750,1000){\circle{30}}\put(800,1000){\mbox{$C_{\textrm{init}}$}}
    \put(750,930){\circle*{20}}\put(800,930){\mbox{$C$}}}%
  \textcolor[rgb]{0.15,0.13,0.73}{
    \qbezier[0](600.7,1342)(551.4,1304.8)(509.4,1259.7)
    \qbezier[0](509.4,1259.7)(467.4,1214.5)(433.9,1162.7)
    \qbezier[0](433.9,1162.7)(400.4,1110.9)(376.4,1054)
    \qbezier[0](376.4,1054)(352.5,997.1)(338.9,936.9)
    \qbezier[0](338.9,936.9)(325.3,876.7)(322.4,815.1)
    \qbezier[0](322.4,815.1)(319.5,753.5)(327.4,692.3)
    \qbezier[0](327.4,692.3)(335.3,631.1)(353.8,572.2)
    \qbezier[0](353.8,572.2)(372.3,513.4)(400.8,458.6)
    \qbezier[0](400.8,458.6)(429.3,403.9)(466.9,355)
    \qbezier[0](466.9,355)(504.5,306.1)(550.1,264.5)
    \qbezier[0](550.1,264.5)(595.7,222.9)(647.8,189.9)
    \qbezier[0](647.8,189.9)(700,156.9)(757.1,133.6)
    \qbezier[0](757.1,133.6)(814.2,110.2)(874.5,97.2)
    \qbezier[0](874.5,97.2)(934.8,84.2)(996.5,81.9)
    \qbezier[0](996.5,81.9)(1058.2,79.6)(1119.3,88.1)
    \qbezier[0](1119.3,88.1)(1180.4,96.6)(1239.1,115.7)}%
\end{picture}\end{minipage}\end{center}
\end{figure}

Averaging these circumcenters gives the initial value for the circle
center:
\begin{equation*}
C_\text{init} (98.680, 47.345)
\end{equation*}

The initial value of the cost function is 13.407.

Using the \textsc{Levenberg-Marquardt} method with very stringent
convergence thresholds ($10^{-10}$ for all tolerance settings:
relative cost, relative parameters and orthogonality), convergence is
reached after only 4 cost function evaluations and 4 jacobians
evaluations. Using the conjugate gradient method with a loose
convergence setting for the inner loop and a tight one for the outer
loop ($\varepsilon_\text{inner}=0.1$ and
$\varepsilon_\text{outer}=10^{-12}$), convergence is reached after 7
iterations.

Both methods give the same results\footnote{a previous version of this
paper found slight differences between the results and a large
difference in the number of iterations, this was due to an
implementation bug in the conjugate gradient method}. The minimal
value found for the cost function is 3.127. The parameters for the
best fitting circle and the corresponding residuals are:
\begin{equation*}
\left\{\begin{gathered}
C (96.076, 48.135)\\
r =69.960
\end{gathered}\right.
\Rightarrow\left\{\begin{aligned}
d_1-r &= -0.963\\
d_2-r &=  1.129\\
d_3-r &= -0.417\\
d_4-r &= -0.475\\
d_5-r &=  0.725
\end{aligned}\right.
\end{equation*}

This example is only one example, no conclusions should be drawn from
the similar results and iteration numbers. Some key features that
depend on the problems and are not addressed in this paper are the
following ones:

\begin{description}
\item[number of points:] this example uses 5 points, other may have
  hundreds of points,
\item[error level:] errors range to about 1\% of
  the circle radius, I have encountered real life applications were
  errors were around 0.001\%, other were errors were about 30\%,
\item[error distribution is homogeneous:] but in some application
  errors may be highly correlated (dented or squizzed circles),
\item[circle coverage:] the example points all belong to the same left
  half of the circle, in some cases points will be regularly
  distributed over the circle, and in some cases there will be some
  strange distribution like 3.5 turns.
\end{description}

\clearpage\appendix
\section{Algorithms}


\begin{algorithm}[hp]\caption{initialization}\label{alg:initialization}
\begin{center}\begin{minipage}{0pt}\begin{tabbing}
xx\=xx\=xx\=xx\=xx\kill
$\sigma_x \leftarrow 0$\\
$\sigma_y \leftarrow 0$\\
$q \leftarrow 0$\\
\textbf{loop for} $i=1$ \textbf{to} $i=n-2$\+\\
  \textbf{loop for} $j=i+1$ \textbf{to} $j=n-1$\+\\
    \textbf{loop for} $k=j+1$ \textbf{to} $k=n$\+\\
      compute $\Delta$ using equation (\ref{eqn:determinant})\\
      \textbf{if} $|\Delta| > \varepsilon$ \textbf{then}\+\\
        \textsl{// we know the points are not aligned}\\
        compute $C_{i,j,k}$ using equation (\ref{eqn:circumcenter})\\
        $\sigma_x \leftarrow \sigma_x + x_{C_{i,j,k}}$\\
        $\sigma_y \leftarrow \sigma_y + y_{C_{i,j,k}}$\\
        $q \leftarrow q + 1$\-\\
      \textbf{end if}\-\\
    \textbf{end loop}\-\\
  \textbf{end loop}\-\\
\textbf{end loop}\\
\textbf{if} $q = 0$ \textbf{then}\+\\
  \textbf{error} all points are aligned\-\\
\textbf{end if}\\
\textbf{return} circle center $C(\sigma_x/q, \sigma_y/q)$\\
\end{tabbing}\end{minipage}\end{center}
\end{algorithm}

\begin{algorithm}[hp]\caption{conjugate gradient}\label{alg:conjugate-gradient}
\begin{center}\begin{minipage}{0pt}\begin{tabbing}
xx\=xx\=xx\=xx\=xx\kill
compute $\hat{r}$ using equation (\ref{eqn:circle-radius})\\
compute $J$ using equation (\ref{eqn:cost-value})\\
compute $\nabla J$ using equation (\ref{eqn:cost-gradient})\\
\textbf{if}
  $|J|<\varepsilon_1$ \textbf{or} $|\nabla J|<\varepsilon_2$
  \textbf{then}\+\\
  \textbf{return}\textsl{ // we consider we already have the minimum}\-\\
\textbf{end if}\\
$J_{\text{prev}} \leftarrow J$\\
$\nabla J_{\text{prev}} \leftarrow \nabla J$\\
$\vec{u}_{\text{prev}} \leftarrow \vec{0}$\\
\textbf{loop from} $i=1$ \textbf{to} $i=i_\text{max}$\+\\
  \\
  \textsl{// search direction}\\
  $\vec{u} \leftarrow -\nabla J$\\
  \textbf{if} $i > 1$ \textbf{then}\+\\
    $\beta \leftarrow
      \nabla J^T (\nabla J-\nabla J_\text{prev})/\nabla J_\text{prev}^2$
     \textsl{// Polak-Ribiere coefficient}\\
    $\vec{u} \leftarrow \vec{u} + \beta \vec{u}_\text{prev}$\-\\
  \textbf{end if}\\
  $\nabla J_{\text{prev}} \leftarrow \nabla J$\\
  $\vec{u}_{\text{prev}} \leftarrow \vec{u}$\\
  \\
  \textsl{// rough minimization along the search direction
          (a few Newton steps)}\\
  \textbf{loop}\+\\
    $J_\text{inner} \leftarrow J$\\
    compute $dJ/d\lambda$ using equation (\ref{eqn:first-derivative})\\
    compute $d^2J/d\lambda^2$ using equation (\ref{eqn:second-derivative})\\
    $\lambda \leftarrow - (dJ/d\lambda)/(d^2J/d\lambda^2)$\\
    $C \leftarrow C + \lambda \vec{u}$\\
    compute $\hat{r}$ using equation (\ref{eqn:circle-radius})\\
    compute $J$ using equation (\ref{eqn:cost-value})\\
    compute $\nabla J$ using equation (\ref{eqn:cost-gradient})\-\\
  \textbf{while} $i\le i_\text{max}$
    \textbf{and} $|J-J_\text{inner}|/J\ge\varepsilon_{inner}$\\
  \\
  \textbf{if} $|J-J_\text{prev}|/J<\varepsilon_{outer}$\+\\
    \textbf{return}\textsl{ // convergence reached}\-\\
  \textbf{end if}\-\\
  \\
\textbf{end loop}\\
\textbf{error} unable to converge in $i_\text{max}$ iterations
\end{tabbing}\end{minipage}\end{center}
\end{algorithm}

\end{document}