MathGroup Archive 2003

[Date Index] [Thread Index] [Author Index]

Search the Archive

Re: Circle Fit

  • To: mathgroup at smc.vnet.net
  • Subject: [mg38633] Re: [mg38555] Circle Fit
  • From: Daniel Lichtblau <danl at wolfram.com>
  • Date: Wed, 1 Jan 2003 03:40:32 -0500 (EST)
  • References: <200212260934.EAA05128@smc.vnet.net>
  • Sender: owner-wri-mathgroup at wolfram.com

Dieter Palme wrote:
> 
> Hi experts,
> I have a set of points (xi, yi), i>30. I search for a solution to fit a
> circle (x-x0)^2 + (y-y0)^2 = r^2 to these points. I need r and sigma_r.
> It could be helpful to get x0,y0 also, but it is not nessecary.
> I found a algorythm for another system (LeastSquareFit) but not for
> Mathematica 4.
> Who can help?
> Thanks in advance
> Dieter
> 
> --
> Dieter Palme
> dl7udp
> 
> mailto:dl7udp at darc.de
> mailto:dieter_palme at t-online.de

I don't know what sigma_r is but presumably what you want to estimate
are radius and coordinates of the center. The best approach I know
starts with the observation that it can be cast as a linear least
squares problem. To do this you write the circle equation as

x^2+y^2 = a*x+b*y+c where you want to find the parameters {a,b,c}. So we
set up a matrix equation using the given data as follows.

{{x1,y1,1},{x2,y2,1},...,{xn,yn,1}} . {a,b,c} =
  {x1^2+y1^2,x2^2+y2^2,...,xn^2+yn^2}

This can be solved in least squares by a number of techniques using e.g.
singular values or Cholesky or Q-R decomposition. For some code details
see the Further Examples in the Mathematica Help Browser under each of
those various functions. In the code below I use the QRDecomposition
approach but the others are substantially similar.

The next step is to convert to the form (x-x0)^2 + (y-y0)^2 == r^2. This
is readily accomplished once we have the parameters {a,b,c} from the
original form of the equation. Note that this readily extends to higher
dimension so the code does not insist that the data be in two
dimensions. The code below returns a pair of the form {radius, center
coordinates}. We could make the code a bit more terse but it is already
reasonably short.

hypersphereFit[data_List] := Module[
  {mat, rhs, qq, rr, params, cen},
  mat = Map[Append[#,1]&, data];
  rhs = Map[#.#&, data];
  {qq, rr} = QRDecomposition[mat];
  params = LinearSolve[rr, qq.rhs];
  cen = Drop[params,-1]/2;
  {Sqrt[Last[params] + cen.cen], cen}
  ]

We generate some data to check this. It is intentionally made to lie
near just a part of the circle; we keep the angle between 0 and 1 in
radian measure and we use random noise so that it does not lie exactly
on the circle. The "base" circle has radius 2 and center {1,-3}.

points = Table[theta=Random[];
  {2*Cos[theta]+1+.05*Random[],2*Sin[theta]-3+.05*Random[]}, {30}]

You can readily chec that the points lie "almost" on a circle.

ListPlot[points]

Now we recover our approximate circle.

In[55]:= hypersphereFit[points]
Out[55]= {1.95099, {1.06397, -2.9453}}


There was a prior discussion on this topic in MathGroup around April or
May of 1995. I'm not sure which parts of that thread can be found in a
web search but I do recall that some of the techniques presented in
responses made use of linear least squares.


Daniel Lichtblau
Wolfram Research


  • Prev by Date: Re: Numbering of figure captions
  • Next by Date: Re: OOP in Mathematica
  • Previous by thread: Re: Circle Fit
  • Next by thread: Re: Circle Fit