Saturday, April 5, 2014

Best Fit Circle: 3D considerations

Somebody recently searched for a 3D best fit circle and I realized there may be some special considerations if your circle is on a slope.  There would also be some tweaks needed if you wanted a best fit sphere instead.

Sphere

First, I should mention, it's possible that the searcher was really looking for a sphere because I have observed that many people don't remember specific terms and use an analogy from a 2D object that they do know the name of and prefix it with "3D" (there, all better now; Darren: *rolls eyes*). I'm not one to criticize analogies, though I sometimes wonder if people realize they are making one. (If you tell me you're making an analogy, I promise not to roll my eyes.  But tell me up front, so I can catch myself.)

Anyways, if it's a sphere you want a best fit of, the equations for the sum of the square of the error (SSE) and radius (r) become:
and

However, if you really want a circle on a sloped plane, I don't recommend using the spherical treatment as a "cheat".  I recommend the following method.

Circle on a Sloped Plane

Before you get carried away on this method, you need to decide if this is really what you need. If you want a shape which will make a circle in plan view (top view), then you need to project your points onto that plane (the xy-plane). This is the same as ignoring the elevation (z-axis). In practice, I'd guess this is really what most people would want. It's true that it will actually be an ellipse when viewed in the plane of the sloped surface; albeit, on only slightly sloping terrain, the difference between that ellipse and a circle probably does not affect the price of rice in China (as we say). If you think the plan view should show a circle, just lop off the z-coordinate and use a plane old 2D method—that is what you actually want.

But if the slope is considerable and you are sure that you want a best fit circle in the plane of that sloped (and basically planar) surface, there are some steps you can take to obtain transformed data to feed into a regular planar method.
  1. Find the best fit plane for your points.
  2. Determine a transform matrix (T) to rotate your plane into the xy-plane as well as its inverse (T–1).
    • You may use the normal vector of the best fit plane to determine this by using two rotation matrix transforms in composition.
    • You will also need to watch for translation. The plane, prior to rotation, should intersect the origin to work properly. It may be appropriate to use an augment (4x4 matrix) to incorporate translation into the transform as well—this way the inverse will also apply the translation as well, rather than having a separate step (only one point though, so no biggie, time-wise).
  3. Apply the transform T to your points to obtain points in the xy-plane.
    • The two rotation matrices can be premultiplied so that you only have to carry out the multiplication between a single transform matrix and the point data. (Recall that matrix multiplication is associative.)
  4. Ignore the z-coordinates of the transformed points. Set them to zero or produce (x,y) data by the trivial xy projection [(x,y,z) => (x,y)] only if this is expected by your 2D best fit circle function or method.
  5. Apply 2D methods to obtain a best fit circle on the xy-plane data.  See here for how to set this up in Excel.
  6. Apply the inverse transform (T–1) to your resulting best fit center to put it into the best fit plane.
Tip: Be watchful of whether your rotational transforms are assuming column vectors or row vectors. Many textbooks give column vector representations, but some programs may assume a row vector.  In the case of Maxima, it will assume a row vector if the vector is given first and a column vector if it is given last:

m1: matrix([c1,-s1,0],[s1,c1,0],[0,0,1]);
r: [a,b,c];
m1.r;            /* column vector representation */
r.transpose(m1); /* row vector representation */

(%o1) matrix([c1,-s1,0],[s1,c1,0],[0,0,1])
(%o2) [a,b,c]
(%o3) matrix([a*c1-b*s1],[a*s1+b*c1],[c])
(%o4) matrix([a*c1-b*s1,a*s1+b*c1,c])