Copyright (c) 2016, Peter Bone
All rights reserved.
Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are
met:
* Redistributions of source code must retain the above copyright
notice, this list of conditions and the following disclaimer.
* Redistributions in binary form must reproduce the above copyright
notice, this list of conditions and the following disclaimer in
the documentation and/or other materials provided with the distribution
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
POSSIBILITY OF SUCH DAMAGE.
Rumen Georgiev (view profile)
jack (view profile)
hello Peter,
i wanted to ask whether your code works on n-dimensions or not?
If not then can you guide in achieving my purpose?
please advise.
Thanks in advance.
Juernjakob Dugge (view profile)
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.