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