Best Fit Circle: find the center using Excel

Finding the center of a best fit circle depends on minimizing the same function we were concerned with in finding the radius, except that we are going to view it as a function of the center instead of as a function of the radius:

where a and b are the x- and y-coordinates, respectively, of the center and r is given by

What we have here is a function of two variables. It looks like three variables until you realize that r is calculated in terms of the other two. So, we can do a three dimensional plot and see what the scoop is. I used Maxima to do this and obtained a very good view of the surface near the best fit center of the points I have been using in all of my investigations of this problem. Here is the 3D plot of SSE(a,b):

What we are most encouraged to see in this graph is that it looks very smooth and it looks like there is exactly one point that is the lowest point. This lowest point is where the SSE function is minimized and constitutes the best center of the circle. (It might be that there are a few local minima somewhat close together that we could see if we zoomed up really tight to the bottom and we are probably happy with any of these as the "answer". Welcome to numerics.)

These formulae can be used in Excel. Designate two cells for each of the values a and b. You don't know what these are, but start with some guesses for these. You will reference these guesses in your Excel formulae. Put your points in consecutive rows after the pattern (x, y, se, R) where se references the x and y for that line as well as the values for a, b, and r. R will only reference x, y for that line and a, b from above. r above is the average of all the R values in the rows (don't include the 1/n in the R)--you may want to create a cell to contain this average and reference it in your se columns. Use absolute references for a, b, and r (if you have a cell for it) so you can copy and paste the formula easily. Make a sum formula at the bottom of your SSE column and it represents your SSE function as above. You want to use the Excel solver now. The SSE cell is the cell you tell it to minimize and the a and b cells you designate as the cells to be modified. The solver will tweak with the a and b values in an attempt to make SSE as small as possible. (The instructions about r and R might seem circular until you actually implement them. Follow through to the end and you'll see it really isn't circular.) Don't try too hard to follow the instructions--try to do the likely intent (as always).

For a Maxima approach see here.

The Polygon Worksheet

I have put together a worksheet in Excel to demonstrate a technique for calculating the vertices of a regular polygon.  You can download the worksheet here.

Overview

It demonstrates an application of vector rotation (and, by the way, complex number multiplication makes a good way to remember how to do vector rotation).  The motivating principle is to demonstrate how you can reduce calls to trigonometric functions while calculating vertices of a regular polygon.  Note that the spreadsheet does not give such (time) savings because it uses both methods: (1) direct calls to trig functions for each vertex in columns B and C and (2) vector rotation in columns E and F.  Also, the error of the vector rotation results is given in columns H and I.  You will notice that the error becomes larger as the vertex number increases.  The error occurs in the vector rotation method, not the trig function method (which is our control, since we believe these results represent the best approximation we can achieve using the floating point precision available in Excel).

The only values you should need to adjust are the number of sides cell, the radius cell, and the initial angle cell.  Everything else follows from these values.  A few notes are in order:

1. The polygon that is determined, is inscribed inside of the circle with the given radius, and centered around the origin. 
2. All angles are specified and calculated in radian measure.
3. If fewer than 50 vertices are needed, the rows will repeat.
4. If you wanted the center somewhere other than the origin, just apply a translation to the points.  So, if you want center (a, b), then each point P = (x, y), becomes P' = (x + a, y + b).  If you've made it this far into this post, you are probably more than capable.

Why Try to Reduce Calls to Trig Functions?

As the user of calculation software and scientific calculators, you may be wondering why we would want to reduce our use of trig functions.  The main reason would be for time-critical applications.  This particular spreadsheet is obviously not one of those, but it demonstrates the accuracy of the method.  So, how do I know that trigonometric functions take a long time to perform.  First of all, let me debunk some myths:

1. Computers and calculators have tables of numbers in them that they use to determine values of "complicated" functions.  So wrong!  Computers and calculators compute/calculate these values!
2. The calculations that computers and calculators do are just based on tables and interpolation.  This is rarely true.  Sometimes a programmer will use a method like this when he knows that the precision needs of his application will be met this way and will be much faster.  But it is actually more programming effort and would only be done when you didn't need very accurate results but you needed them really, really fast.
3. Computers and calculators have hardware circuits that make the complex calculations almost as fast as the simple ones.  Not so.  Fact:  Multiplication is done with the fmult instruction.  Addition is done with the fadd instruction.  fmult takes way, way more computer cycles to compute than fadd.  fsincos takes way, way longer than fmult.  It is an unavoidable consequence of the nature of these calculations.

Recall that multiplication is dependant on addition.  Several additions actually.  So it should come as no surprise that by the time you put together a circuit that does the job that several additions and small multiplications can do, you end up with a circuit that takes longer to complete than one that just does addition.

The same is true of trig functions.  You know how you calculate them?  Here are the magic formulas:

The dots (...) indicate that it goes on and on and on and on and on, until the changes are small enough that you aren't changing the digits that mater to you. You may need more terms than are listed to get the accuracy you need. These formulas above are not written in the most computationally efficient manner, but you can rest assured that any circuitry or coding or combination thereof, which computes the values of the above functions will take longer than 4 multiplications and 2 addition/subtractions will (those are the operations involved in doing vector rotation).

Enter Vector Rotation

The key to how this method helps us, is that we can use a few trig function calls at the beginning and reuse the results several times over. We compute the direction vector for the central angle, θ, that each side covers. For a hexagon, this angle is 360°/6 = 60°, or π/3.

We need to determine the value of the central angle using direct calls to trig functions, but we will be able to reuse those values.  We also call trig functions to get our initial position.  From there we use our vector rotation method.  So, E14 and F14 use the values in E11 and F11 (renamed mx and my) and those of E13 and F13 to find their values:

E14 =E13*mx-F13*my

F14 =E13*my+F13*mx

Cells in rows below 14 do the same thing.  They still reference mx and my, which are the cosine and sine of our central angle, θ, respectively, and they reference the cells directly above them.  The concept for each row is:  rotate the point in the row above me around the origin by the central angle and tell me the x and y values for that rotated point.

Proof of the Formulas

Suppose we are given a point (x, y) and denote the distance from the origin as r and the angle it makes with the positive axis as i.  We wish to rotate (x, y) counter-clockwise about the origin by angle θ.  The diagram below illustrates:

Our new point (x', y') can be simplified:

                  

                      
                      
                 
                     
                     

If you look back at the formulas for E14 and F14, you may observe that mx and my correspond with sin θ and cos θ and E13 and F13 correspond with x and y.  That's all there is to it. 

Relationship with Complex Number Multiplication

If we represent the new point (x', y') as x' + y'j (where j is used to indicate the square root of -1), we can say:

This will produce an equivalent result to what we have above.