Find the approximate center of a circle passing through more than three points

Question

Consider n point $(x_1,y_1), (x_2,y_2),\ldots, (x_n,y_n)$. For $n = 3$ it is easy to find the center of the circle passing through the three points. I wanted find the approximate center of the circle passing through more than three points.

The application is as follows. I have circles on a sheet of metal. I can poll different points on the circle to get the their co-ordinates. By polling three points on every circle I can predict the co-ordinates of the center. But I wanted to improve the accuracy of the prediction. A Google search reveals only the center of the passing through three points.

Edit 1: In response to Hoda's comment, I am adding that the measurements related to the points on the circle have errors. The objective is to minimize the error in the position of the center by polling more points on the circle.

Answer 1

One can find on the web several methods of cercle fitting. Most of them are itterative. A straightforward method, without trial and error process, is described pp.12-13 in the paper "Régressions coniques, quadriques, circulaire, ..." : http://www.scribd.com/JJacquelin/documents No need to read the theoretical part written in French. Just apply the formulas (1), (2) and (3) which allows to very easily compute the coordinates of the center and the radius of the fitted circle. A numerical example is provided page 15.

Answer 2

Introducing a bias, consider the equations $$F_i=(X-x_i)^2+(Y-y_i)^2-r^2=0$$ You can write $$F_j-F_i=2(x_i-x_j) X+2(y_i-y_j)Y=(x_i^2+y_i^2)-(x_j^2+y_j^2)$$ This is identical to what Steven Taschuk wrote.

Generate all equations for $(i=1,2,\cdots,n-1)$, $(j=i+1,i+2,\cdots,n)$. This will create $\frac{n(n-1)}2$ equations and defining the regressors you basically face the problem of a least square fit (with no intercept) and the calculation is immediate (using matrix will make the problem very easy). So you will have very good estimates of the coordinates of the center and the full nonlinear regression will converge in a couple of iterations since averaging all the $F_i$'s, you have a good estimates of $r^2$ then $r$.

Applied to the data points I used in my previous answer, this procedure leads to $X=5.40887$, $Y=8.42447$. The average radius is $r=3.41043$.

Now, any procedure will solve in a couple of iterations to $$\begin{array}{clclclclc} \text{} & \text{Estimate} & \text{Standard Error} & \text{Confidence Interval} \\ X & 5.39989 & 0.0928415 & \{5.10443,5.69535\} \\ Y & 8.42763 & 0.128379 & \{8.01907,8.83619\} \\ r & 3.41031 & 0.0760165 & \{3.16839, 3.65223\} \\ \end{array}$$

Answer 3

Recalling the circle equation with center $(h,k)$ and radius $r$

$$ (x-h)^2+(y-k)^2 = r^2 .$$

The above equation has three unknowns which can be determined uniquely by knowing three points lie on the circle. However, if you have more than three points, then you may need to appeal to least squares method. See related technique.

Answer 4

Assuming normally distributed errors, a maximum-likelihood estimate for the circle parameters $(x_0, y_0, r)$ is one that minimizes the sum of the squared distances from the circle: $$ SS(x_0,y_0,r)=\sum_{i=1}^{n}\left(\sqrt{(x_i-x_0)^2 + (y_i-y_0)^2}-r\right)^2. $$ Taking the derivative with respect to $r$ gives $$ 0=\frac{\partial}{\partial r}SS(x_0,y_0,r)=-2\left(\sum_{i=1}^{n}\sqrt{(x_i-x_0)^2+(y_i-y_0)^2}-nr\right), $$ so $$ r^{*}=\frac{1}{n}\sum_{i=1}^{n}\sqrt{(x_i-x_0)^2+(y_i-y_0)^2}. $$ A numerical search for coordinates $(x_0,y_0)$ that minimize $SS(x_0,y_0,r^*)$ should be relatively straightforward.

Answer 5

The answer given by mjqxxxx is very interesting. What he shows is that the problem can be reduced to $x_0$ and $y_0$ since $r$ is eliminated. Keeping the same objective function as the one he proposed, just use Newton-Raphson method to compute the values of $x_0$ and $y_0$ for which the derivatives cancel. $$\frac{\partial}{\partial x_0}SS(x_0,y_0,r)=0$$ $$\frac{\partial}{\partial y_0}SS(x_0,y_0,r)=0$$

Since I suppose that you have "reasonable" estimates of the variables, the convergence will be quite easy. You could even be lazy and use numerical derivatives to compute the Jacobian of the system.

I applied the method using six data points : $(2.5,6.5), (2.0,8.5), (4.0,11.5), (7.5,11.0), (9.0,9.0), (8.0,6.5)$.
The calculations solved in three iterations and the circle is such that $r=3.41031$, $x_0=5.39989$, $y_0=8.42763$.

If you apply to the same data the method proposed by JJacquelin, you directly arrive to $r=3.41287$, $x_0=5.40887$, $y_0=8.42447$. This method is from far superior to what I proposed since it is a direct and explicit method which just requires the computation of a few sums over the coordinates of the data points.

Answer 6

The perpendicular bisector of the line segment joining point $(x_i,y_i)$ and point $(x_j,y_j)$ has equation $$ (x_i-x_j)x+(y_i-y_j)y = \tfrac12(x_i^2+y_i^2-x_j^2-y_j^2) $$ If all the given points fell on one circle, all these lines would intersect, and you could just solve the system. If they do not intersect, you have an overdetermined system of linear equations, for which you can compute an approximate solution by least squares.

This method has the advantage of only using very standard tools.

Answer 7

You have $n$ equations that look like this, with $x_i$ and $y_i$ being the known points, $a$ and $b$ being the unknown center coordinates, and $r$ the unknown radius. $$ (x_i - a)^2 + (y_i - b)^2 - r^2 = 0$$

Multiply out all the squared differences, and consolidate all the constant terms on the right. Now you have $n$ equations that look like this:

$$e_i a + a^2 + f_i b + b^2 - r^2 = g_i$$

Subtract the first equation from all the rest. Now you have the first equation same as before, but the rest for $2 \le i \le n$ look like: $$(e_i - e_1)a + (f_i - f_1)b = g_i - g_1$$

Since $n \gt 3$, you have at least $3$ of these linear equations with 2 unknowns, which you can solve for $a$ and $b$, giving you an overdetermined linear system. Rewrite this system as $\mathrm{A}x = g$, then you can solve for $x$ in a least squares sense with $x = (\mathrm{A}^T\mathrm{A})^{-1}\mathrm{A}^Tg$.

Once you have estimated $a$ and $b$, you can easily estimate $r$ with root-sum-square on all the calculated distances between the $(x_i,y_i)$ and $(a,b)$.

Answer 8

First you need to define what your metric for error is, and then you can talk about minimizing it. For example, in simple linear regression, the error is measured by summing the squares of the vertical displacement from each known datum $(x_i,y_i)$ to the line. But in exponential regression, the error is measured by summing the squares of the logarithms of the ratios between the known data $(x_i,y_i)$ and the exponential curve.

I can imagine several candidates here for quantifying error. Simply summing the distances from the data to the circle (where by distance, I mean shortest distance). Or summing the squares of these distances. The somehow doesn't feel right though, since it seems reasonable to make a point at the circle's center contribute a lot more error than one which is one radius removed externally from the circle. Some metric that put the center of the circle infinitely far away from the circle itself might be nice.

Answer 9

Since you did not specify in which way the center should be uniquely defined, I'll just present some thoughts based on the observation that the usual method of finding the center of a circle through three points is really nothing else but finding the intersection of two or three perpendicular bisectors of the lines connecting those points.

One method often used in triangulation on a map (using a magnetic compass etc) is to take more than the strictly necessary two lines, and take all the intersections as approximations to the true position. You could do some statistical analysis to discard outliers, you could take the barycenter of their convex hull, or simply their average.

Answer 10

Just including a python (v2.7) code below (with the least squares method) for anyone who wants to try it out using a computer:

# Original source: http://www.had2know.com/academics/best-fit-circle-least-squares.html
# Notes: I have slightly modified the original code to suit my ends. There are 
# quite a few errors on the source link. Try out their example and it should be clear.

import numpy as np

# Initialize the lists for the abscissae and the ordinates
x_l=[]
y_l=[]
# Get the x and y coordinates of the points
while True:
    x_temp = raw_input("Enter the x-cordinate (press n to initiate fitting):\n")                       
    if x_temp.strip() == "n":
        break
    y_temp = raw_input("Enter the y-cordinate:\n")
    print "------------------------------------------------------"
    x_l.append(float(x_temp))
    y_l.append(float(y_temp))
nop = len(x_l)                          # number of points

# Primary raw material
x = np.array(x_l)                       #numpy array
y = np.array(y_l)                       #numpy array

x_y = np.multiply(x,y)                  #elementwise multiplication
x_2 = np.square(x)                      #elementwise square
y_2 = np.square(y)                      #elementwise square

# Secondary raw material
x_2_plus_y_2 = np.add(x_2,y_2)          #elementwise addition
x__x_2_plus_y_2 = np.multiply(x,x_2_plus_y_2)   #elementwise multiplication
y__x_2_plus_y_2 = np.multiply(y,x_2_plus_y_2)   #elementwise multiplication

# Summations (by default the result of summation is integer)
sum_x = x.sum(dtype=float)
sum_y = y.sum(dtype=float)
sum_x_2 = x_2.sum(dtype=float)
sum_y_2 = y_2.sum(dtype=float)
sum_x_y = x_y.sum(dtype=float)
sum_x_2_plus_y_2 = x_2_plus_y_2.sum(dtype=float)
sum_x__x_2_plus_y_2 = x__x_2_plus_y_2.sum(dtype=float)
sum_y__x_2_plus_y_2 = y__x_2_plus_y_2.sum(dtype=float)

# Comment the line below; don't delete. Use for debugging.
# print sum_x,sum_y,sum_x_2,sum_y_2,sum_x_y,sum_x_2_plus_y_2,sum_x__x_2_plus_y_2,sum_y__x_2_plus_y_2

m3b3 = np.array([[sum_x_2,sum_x_y,sum_x],
        [sum_x_y,sum_y_2,sum_y],
        [sum_x,sum_y,nop]])
invm3b3 = np.linalg.inv(m3b3)
m3b1 = np.array([sum_x__x_2_plus_y_2,sum_y__x_2_plus_y_2,sum_x_2_plus_y_2])
A=np.dot(invm3b3,m3b1)[0]
B=np.dot(invm3b3,m3b1)[1]
C=np.dot(invm3b3,m3b1)[2]
CenterX = A/2
CenterY = B/2
Radius = np.sqrt(4*C+A**2+B**2)/2
print "Center's x-cordinate is :",CenterX,"\nCenter's y-cordinate is : ",CenterY,"\nCircle's radius is : ",Radius