calcCircle

Find the circle that passes through 3 non-collinear points.

Description

Find the circle that passes through 3 non-collinear points. Uses a fast analytical method based on finding the intersection of the bisectors of 2 of the line segments. Avoids divide by zero errors in all cases. Faster than determinant methods.
Tested under MATLAB 7.5 but should work on previous versions

• MATLAB release: MATLAB 7.5 (R2007b)

[John D'Errico]

This submission and its cousin, points2circle (FEX 19082) were submitted together, in the spirit of friendly competition. Lets take an extra careful look at these two, in that same spirit. How do they compare?

The most obvious question is, does this function work? I'll create a set of three points on a known circle, here with center [1 2], and radius 1.3.

center = [1,2];
theta = [.1 .2 .98]; % radians
radius = 1.3;
circ = @(center,rad,t) center + rad*[cos(t),sin(t)];
p1 = circ(center,radius,theta(1));
p2 = circ(center,radius,theta(2));
p3 = circ(center,radius,theta(3));

[c,r] = calc_circle(p1,p2,p3)
c =
            1 2
r =
          1.3

So, yes, it does indeed seem to work properly, returning the correct center and radius for my test points.

How fast is this code? I've used Steve's code, timeit (find it on the fex!) for the comparison. It gives a good estimate for the time required.

timeit(@() points2circle(p1,p2,p3))
ans =
   0.00047171

timeit(@() calc_circle(p1,p2,p3))
ans =
   0.00017641

Note that both codes are fast, but that calc_circle is nearly 3 times as fast. There really is not a huge amount of work to do, so the time should not be a serious factor, UNLESS you have this function call deep inside nested loops and you need maximum speed. In that event, calc_circle (this code) will give you almost 3x the throughput.

Pure speed is not all that matters however. Do you want the right answer slowly, or the wrong one fast? While speed is important, we also need to look at other factors. Does this code detect collinear data?

p1 = [0 0]; p2 = [1 1]; p3 = [3 3];
[c,r] = calc_circle(p1,p2,p3)
c =
     0 0
r =
    -1

As I might hope, it finds the problem and returns a result that indicates a problem, although no error or warning message was thrown. The user must test the radius to know there was a problem. IMHO, points2circle dealt with this in a more friendly way.

Next, I looked at the interface. The alternative code, points2circle, has a more flexible interface, allowing the user to provide their points in either of two logical forms. However, I'll admit that the calc_circle interface is entirely adequate here.

calc_circle has no error checks on its arguments. This might partly explain its speed. There are no size checks, no checks on the argument types. I did find an H1 line here.

The help for calc_circle is poor, almost useless in fact. Lets try it here:

help calc_circle
  Fit a circle to a set of 3 points
  Returns radius of -1 if points are collinear
 
  Author: Peter Bone
  Date: 6th March 2008

That is all of it! What information does this provide about how to use calc_circle? Almost none! What are the arguments? What type should they be? What sizes? Are these three points that lie in the (x,y) plane? Can it be in higher dimensions? What output(s) would we expect? There is no example of use. I don't want to know the author's name and the date as much as how to use the code. That is what help is for, and this author should know better.

An interesting point to note is whether this code survives a common problem. What happens if your data is not composed of doubles?

p1 = rand(1,2);p2 = rand(1,2);p3 = rand(1,2);

This one works:
[c,r] = calc_circle(p1,p2,p3)
c =
      0.41636 0.48853
r =
      0.34388

Single data is fine:
[c,r] = calc_circle(single(p1),single(p2),single(p3))
c =
      0.41636 0.48853
r =
      0.34388

But no test is done for non-floating point data. calc_circle just returns its silent non-answer.

[c,r] = calc_circle(uint8(p1),uint8(p2),uint8(p3))
c =
     0 0
r =
    -1
 
An often overlooked issue is, can calc_circle survive complex data? It does so nicely. Well done here, since calc_circle uses norm to compute the radius at the end.

[c,r] = calc_circle(p1+i,p2+i,p3+i)
c =
      0.41636 + 1i 0.48853 + 1i
r =
      0.34388

Overall, this code gets a rating of 2 for me, although I did consider a 3 rating. The lack of any useful help, examples, and error checks pushed me down. Your job as a programmer does not end when you write the last line of code. (The author should know this fact.) I liked the better speed of calc_circle, as well as the use of norm for the radius.

With the proper cleanup, I'd rate this code a 5.