Saturday, August 27, 2011

Best Fit Circle to Measured Points

(For a practical finished solution in Excel see here.  This was my first crack at the problem and lacks mathematical rigour and it isn't very precise.  But, if all you need is "close enough" and don't want to take the time to understand how to use the spreadsheet, this might be of some help.)

Suppose you are given a bunch of points - (x,y) coordinates - that you are told are on a circle, or at least very close.  Perhaps it is a very big circle and it is not realistic to measure it with a tape measure.  (Otherwise, you could take a long tape and go around the outside of the circle to determine its circumference, C.  From this you determine its radius, r = C/(2π).)  Or perhaps you think it is a circle, but you'd like to verify that it is sufficiently close.  At any rate, somebody shoots it with a total station and gives you the points.

Now, there you are in AutoCAD with your points.  Sooo...., what's the radius of this thing?  Your first guess might be to try to find points that are almost directly across from each other and call the line between them the diameter.  Then try it with a bunch more.  See if your best estimates are all close to each other and average the close ones out - leaving out the outliers.  Okay, that's not bad, er, not terrible anyways, but also not very "scientific."  (Actually, you should probably take the longer lengths because the longest chord or a circle is its diameter.)  One big disadvantage to this is that you don't have an estimate of how good your approximation is or how well the shape you just "shot" (got points from a total station) conforms to a circle.  Furthermore, how would you make this work if you only have part of an arc and don't even know how many degrees of the circle the arc covers? 90°, 60°, 132.734°?

Let's throw some geometry into the mix.  In AutoCAD I can easily draw "perpendicular bisectors" to several of the chords (just draw lines between several pairs of the given points).  In theory, these perpendicular bisectors should all intersect at the same point.  Since your data is not perfect and the thing being measured is not a 100% perfect circle, there will be many different intersection points.  But, if you have a bonafide circle, you can average out all of these intersection points and call it the center of the circle.  If you're of only average persistence you end up with something like this:
You can enter the coordinates of the different intersections into a spreadsheet and average the x and y coordinates individually and that will give you an average center.  On the other hand, you could put all of the points into a spreadsheet and calculate all of the intersection points you desire - from every combination of 3 points and every combination of 4 points, and then average them all out.  Either way, you have an estimate of the "true", "best", center of the circle.  (The last method - taking all of the intersections - actually doesn't work as well as it might seem on the surface.  I will try to post something about it later.)

Now we need to determine the radius of the circle.  If you're in a hurry, you just put your center at the point you calculated and indicate the radius to AutoCAD by clicking any one of the points you believe to be on the circle.  If you're in an even bigger hurry, you skip calculating the center and just eyeball the center based on your intersections as in the picture above.  On the other hand, a spreadsheet sure helps out a lot here if you don't know for sure which point to pull your radius out to (some might be more in error than others after all).

For this we need the equation of a circle.  Let's say our calculated center is given by (a, b).  Then the equation of the circle we are interested in is given by 

Your first guess might be to find the average distance between the center and each measured point and use that as the radius.  This actually turns out to be the case:

where the xi's and yi's are the points on the arc or circle and n is the number of points. See the development of the above formula here.

No comments: