Problem Corner: Circle Fitting Quest
Professional Surveyor Magazine - July 2011
It’s been about ten years now since the Problem Corner started; the first two problems were published in the double July/August 2001 issue. Back then you had to wait for the next issue of
the magazine for the solutions, so your wait was two months. This article celebrates the anniversary, and it includes the solution with the problem.
While searching for another unique problem for the Problem Corner, I found that I had never submitted a leastsquares-adjusted circle-fitting problem. It’s certainly not something you see or
do everyday, but I think it would be good to know how to do it, especially without a lot of complex mathematics.
I didn’t want to be confined by holding to the tangents of the circle. I didn’t want to be forced to use some integer foot number or whole degree central angle. I wanted it to be practical enough for the average surveyor to use if and when the need arose.
I was sure I could just look in my survey texts and find a good example with an easy explanation. Well, to my surprise, there are no problems involving circle fitting in modern surveying texts that do not involve iterative solutions or matrices. Every good surveyor should know matrix algebra, but I’m sure most don’t, so I elected to leave it out. There are examples of fitting a circle to three points, but I consider that trivial because any good surveyor can fit a circle to three points by intersecting the perpendicular bisectors of the two chords. There are examples of fitting a circle to three points by solving simultaneous equations of a circle using
Even the gold standard, Adjustment Computations, Spatial Data Analysis
by Ghilani and Wolf, didn’t have anything other than fitting a parabola to a few points. The drawback of that problem is that one of the variables is held fixed. I wanted a problem that let both variables have some inherent error, just like in the surveying world.
So, like anyone, I went to the internet and typed “circle fitting” in a search engine. After four hours or so I came across a November 18, 2008 23-page paper by N. Chernov and C. Lesort of the University of Alabama at Birmingham Department of Mathematics titled “Least Squares Fitting of Circles.
” I called John Nolton of Tombstone, Arizona, who was helping me research this problem in his vast library of survey literature, to let him know of my latest find. He called back the next day and asked if I had checked all
the references given in that paper, especially the authors’ website (www.math.uab.edu/cl/cl1
). He said he was going to copy and send to me one of the papers published there. Well, two days later I got a 254-page(!) paper by Nikolai Chernov titled “Circular and Linear Regression: Fitting circles and lines by least squares
This paper explains more than you will ever want to know about circle fitting. It does, however, have just what I wanted: a straightforward solution to fitting a circle to points and all the pros and cons of doing so. I already had the problem drawn up by this time (Figure 1
) and even sent it out to a Problem Corner fan and some friends for constructive criticism.
The Kasa Method
The solution I choose to highlight here, the Kasa method, definitely has some weak points. (I will explain them later.) It does have a closed, non-iterative solution that can be solved by anyone with an algebra background. It is fast and easy. To be absolutely correct, the sum of the residuals should be calculated as
One of the weaknesses of the Kasa method is that the original points must be “close” to the circle to start. Because the seven points in Figure 1 are on opposite sides of the street, they must be transferred to the centerline. To move them radially towards the centerline, you must estimate a radius point. Didn’t I just say any good surveyor could do that?
Using points a, b, and d, calculate the midpoints of each chord by averaging the coordinates of each end: the midpoint of chord a-b is North 887.38, East 798.39 and the midpoint of chord b-d is North 1138.55, East 905.395. The bearing of chord a-b is N 11°04’00” E and the bearing of chord b-d is N 30°16’04” E. Intersecting the perpendiculars to the midpoints gives a radius point of N 729.29, E 1606.66 with a radius of 830.12 through these three points.
If you would rather use the equation of a circle (remembering that Northing=y and Easting=x), x² + y² + 2dx + 2ey + f = 0 for a circle centered at (-d,-e) with radius equal to
Subtracting 700 from all northings and eastings and using points a, b and d again:
78.44² + 85.38² + (2)(78.44)d + (2)(85.38)e + f = 0
118.34² + 289.38² + (2)(118.34)d + (2)(289.38)e + f = 0
292.45² + 587.72² + (2)(292.45)d + (2)(587.72)e + f = 0
from which d = –906.67, e = –29.29 and f = 133,796.995 and the radius of 830.12 is centered at North 729.29, East 1606.67.
With the approximate radius just calculated, move the points 33.000 feet towards the centerline. These are the values we will use in the Kasa method.
a to a’: S 86°07’32” E 33.000 = N 783.150, E 811.365
b to b’: S 71°44’27” E 33.000 = N 979.040, E 849.678
c to c’: S 54°38’20” E 33.000 = N 1190.572, E 956.642
d to d’: S 47°43’24” E = N 1265.521, E 1016.867
e to e’: N 57°39’45” W 33.000 = N 1155.652, E 933.198
f to f’: N 68°11’51” W 33.000 = N 1025.346, E866.560
g to g’: N 80°26’39” W 33.000 = N 861.608, E 820.658
Be careful to note the equations are set up for (x, y) notation.
You can grind out the multiplications, additions, and divisions on your basic calculator or you can enter them into a spreadsheet and let it do all the work. I subtracted 700 from all the northings and 800 from all the eastings. Figure 2
is my printout. The top line explains what is in the columns. The bottom line is the coefficients for the linear equations:
\B(13,823.396) + C(42,854.233) + D(93.567) = –23,296,374.337
B(42,854.233) + C(140,691.149) + D(337.270) = –71,916,197.444
B(93.567) + C(337.270) + D = –154,514.545
When you get all of the equations set up, you can solve them by yourself, but why do that when the internet is available to do it for you? Typing “3 x 3 equation solver” into the search boxwill give you many choices. I used the one from Penn State Erie just because it was simple and came up first.
The final values of B = –1,613.582, C = –58.449 and D = 16,176.523 can be converted by
to a = 806.791, b = 29.224 and R = 797.238.
Then adding back the 700 to the northing and the 800 to the easting gives a radius point of N 729.224, E 1,606.791.
The final residuals, using R and the radius point fixed:
a: 830.252 – 33 – 797.238 = +0.014
b: 830.263 – 33 – 797.238 = +0.025
c: 830.205 – 33 – 797.238 = –0.033
d: 830.261 – 33 – 797.238 = +0.023
e: 764.225 + 33 – 797.238 = –0.013
f: 764.264 + 33 – 797.238 = +0.026
g: 764.202 + 33 – 797.238 = –0.036
One other weak link of this method is that small portions of arcs, typically less than 30°, will yield ridiculously wild radii. Just try entering the raw coordinate data into the spreadsheet without reducing to centerline and you will see what I mean: a short radius circle nestled amongst all the points. We are tempting fate with 40° of arc.
Note also that the sums of squares and cross multiples of the coordinates can give you extremely large numbers. These “moments” as they are called can get unwieldy. That was the reason I subtracted 700 from the northings and 800 from the eastings.
» Back to our July 2011 Issue