From: Andrew on
Hi Roger,
Thank you. I found your short sample code and explanation very useful. With your permission, I would like to use it in an open-source (GPL) analysis tool I am writing. I would acknowledge you and reference this posting. I can't seem to contact you via email directly. Please let me know how best to reach you.

Best,
Andy


"Roger Stafford" <ellieandrogerxyzzy(a)mindspring.com.invalid> wrote in message <h554pn$dmf$1(a)fred.mathworks.com>...
> "Ruchir Srivastava" <ruchir(a)nus.edu.sg> wrote in message <h418sh$p59$1(a)fred.mathworks.com>...
> > Hi Roger,
> > The code is good. Can you give me any reference for the mathematics used. I could not understand the use of variances.
> > Regards,
> > Ruchir
> ------------------------------
> Hello Ruchir Srivastava,
>
> This two-year-old thread predates my Watch List, and I didn't discover your question until today. My apologies for the two-week delay.
>
> As you recall, I defined X and Y (uppercase) by X = x - mean(x) and Y = y - mean(y) where x and y (lowercase) are the original coordinate vectors. We can thus express the quantity whose mean square value is to be minimized as
>
> (X-a0)^2 + (Y-b0)^2 - r^2.
>
> For any given fixed values of a0 and b0, it can readily be seen by taking the partial derivative of its mean square with respect to r^2, the mean square will reach a minimum with respect to r variation at the value
>
> r^2 = mean((X-a0)^2+(Y-b0)^2).
>
> Substituting this for r^2 in the above gives an expression involving only the two parameters a0 and b0:
>
> (X-a0)^2 - mean((X-a0)^2) + (Y-b0)^2 - mean((Y-b0)^2)
> = -2*X*a0 - 2*Y*b0 + X^2 - mean(X^2) + Y^2 - mean(Y^2)
>
> for which we must find the least mean square value (taking advantage of the fact that mean(X) = mean(Y) = 0.) This is equivalent to finding a0 and b0 for the least mean square approximation to the equations
>
> a0*X + b0*Y = (X^2-mean(X^2)+Y^2-mean(Y^2)) / 2
>
> As you can see, this is the least squares problem that I used the backslash operator to solve. The final determinations of a, b, and r then readily follow from this.
>
> Roger Stafford
From: Roger Stafford on
"Andrew " <leifer(a)fas.harvard.edu.removethis> wrote in message <i1n85c$dkn$1(a)fred.mathworks.com>...
> Hi Roger,
> Thank you. I found your short sample code and explanation very useful. With your permission, I would like to use it in an open-source (GPL) analysis tool I am writing. I would acknowledge you and reference this posting. I can't seem to contact you via email directly. Please let me know how best to reach you.
>
> Best,
> Andy
- - - - - - - - -
Andy, I am not very familiar with GPL licensing provisions, but assuming that they do not restrict the further free distribution of such computer code, I would have no objection to your plans to include my piece of matlab code from this thread in your "tool".

It is important to realize both the strengths and weaknesses of this circle-fitting code. It does well for points that lie in an approximate arc which covers a large fraction, say a half or a third, of a complete circle even if there is considerable noise present in the points. It also does well for points that lie in an arc which constitutes only a small fraction of a circle provided there is comparatively little noise present - that is, provided the points do actually lie close to a circular arc.

It does not do so well in cases where there is substantial noise present in the points' locations and only a small fraction of a full circle is occupied. By looking at matlab plots I have made of an ideal circle with associated noisy points on a small arc of it and comparing these with the algorithm results, it is clear that in such difficult cases, a human could very likely draw a better circle to fit the points than is generated by the code. To come up to a human's level of performance would probably require a comparatively sophisticated pattern recognition algorithm rather than the present algorithm's simple least squares approach.

It is also important to point out to potential users that the code assumes the cartesian coordinates of points in two-dimensional space to be fitted to a circle have been previously stored in two column vectors which are referred to as x and y. Row vectors would come to grief here with matlab's backslash operator.

Roger Stafford