Code covered by the BSD License

Highlights from fitcircle.m

4.66667
4.7 | 6 ratings Rate this file 33 Downloads (last 30 days) File Size: 39.3 KB File ID: #15060 Version: 1.0

fitcircle.m

Richard Brown (view profile)

20 May 2007 (Updated )

Fits circles to 2D data using nonlinear least squares to minimise geometric error

File Information
Description

Although a linear least squares fit of a circle to 2D data can be computed, this is not the solution which minimizes the distances from the points to the fitted circle (geometric error). The linear solution minimizes the algebraic error of a function something like
f(x) = ax'x + b'x + c = 0

Minising the geometric error is a nonlinear least squares problem. fitcircle allows you to compute either - it uses the algebraic fit as the initial guess for the geometric error minimization.

e.g.
x = randn(2, 10);
% Linear least squares fit
[z, r] = fitcircle(x, 'linear')

% True best fit (minimizing geometric error)
[z, r] = fitcircle(x)

This submission is based on the paper:
"Least-squares fitting of circles and ellipses", W. Gander, G. H. Golub, R. Strebel, BIT Numerical Mathematics, Springer 1994

A similar submission for ellipses should be forthcoming

MATLAB release MATLAB 7.1.0 (R14SP3)
Other requirements Should work for R14 onwards. Changes to make it work for earlier releases are probably fairly straightforward (e.g. replacing boolean false and true with 0, 1)
18 Jul 2016 Aaron Greenbaum

Aaron Greenbaum (view profile)

07 Sep 2012 Richard Brown

Richard Brown (view profile)

@Tolga, thanks -- and yes, I know. I deliberately wanted to avoid any toolbox dependencies.

I haven't done any comparisons of efficiency -- if you do you could let us know! But accuracy is simply dependent on the tolerance you specify.

Comment only
19 Aug 2012 Tolga Birdal

Tolga Birdal (view profile)

Richard,

I appreciate your effort and this great function. I would just want to remind that (and I am sure you already know) for the nonlinear part it is also possible to use:

options = optimset('Jacobian','on', 'Algorithm','levenberg-marquardt');
u = lsqnonlin(@sys,u,[],[],options);

Did you do any comparisons of accuracy? (Because Walter states that Gauss-Newton performed similarly, while being more efficient)

Cheers,

Comment only
21 May 2012 Richard Brown

Richard Brown (view profile)

My previous comment got lost -- the 2 norm residual is the sum of the squared perpendicular distance from each data point to the fitted circle

Comment only
21 May 2012 Richard Brown

Richard Brown (view profile)

Also, @Graeme, the accuracy of the fitted circle depends on the kinds of errors that are in your data. If you assume there is a true underlying set of parameters that you're trying to find, and that your data is normally distributed, then the accuracy will decrease with number of data points as something like s / sqrt(n), where s is the standard deviation of perpendicular distances to the fitted circle, and n the number of data points. The error will then probably follow some kind of t distribution.

Comment only
10 May 2012 Sharif

Sharif (view profile)

IGNORE last comment.
http://www.mathworks.com/matlabcentral/fileexchange/15125-fitellipse-m

Comment only
10 May 2012 Sharif

Sharif (view profile)

Hi Richard,

Thanks you for sharing your script.

i was wondering if it possible to modify your script to fit an ellipse. if so, do you have any suggestions?

thanks again,
sharif

Comment only
07 May 2012 Graeme Crouchley

Graeme Crouchley (view profile)

Hi all,

I was wondering what exactly the 2 norm residual is? Is it the sum of the residuals?

Also how would I use this value (or any other) to determine the accuracy of the fitted circle. i.e what accuracy has the radius of the circle been determined to?

Thanks,

Graeme

Comment only
04 Apr 2012 Daniel

Daniel (view profile)

thanks,

Comment only
01 Apr 2012 Richard Brown

Richard Brown (view profile)

@Daniel. Yes, it is very easy to do this, and you don't need my file (it has a simple solution).

Assume that you want to force the centre to be 0,0. (if something else, then you can just change coordinates to solve the problem, then change back afterward).

Let your points be in a 2xN array called x. Then the objective you're minimising is
f(r) = \sum_i (\norm{x_i} - r)^2
where x_i is the i-th column of x.

You can minimise using basic calculus
(1) f'(r) = -2\sum_i (\norm{x_i} - r) = -2\sum_i \norm{x_i} + 2Nr
(2) f''(r) = 2N

(1) has a unique minimiser at
r = \frac{\sum_i \norm{x_i}}{N}
and (2) confirms that the function is convex

i.e. in Matlab notation:
r = sum(hypot(x(1, :), x(2, :)) / N;

cheers,

Richard

Comment only
29 Mar 2012 Daniel

Daniel (view profile)

29 Mar 2012 Daniel

Daniel (view profile)

Great and simple. Any idea if it would be easy to force the center to zero,zero (or anything other coordinate)?

Comment only
19 Mar 2012 Akash Yalagach

Akash Yalagach (view profile)

18 Sep 2011 Richard Brown

Richard Brown (view profile)

Hi Brian

This code computes a least squares fit, which is not robust to large outliers - the results you got are not unexpected, nor are they wrong.

Because the objective being minimised is a function of the square of the distance, big outliers will skew the fit. If your data is known to have significant outliers, least squares is not a good choice of objective to minimise.

cheers,

Richard

Comment only
30 Aug 2011 Brian Wheeler

Brian Wheeler (view profile)

To test thus function I fed it data for a 360 point unit circle. As expected it returned a unit circle centered on zero.

Distorting one point of the circle to a radius of 6 gave expected results:
center [ 3.76743045176789e-17 0.0281775669475842]

Distorting the same point of the circle to a radius of 7 gave UNEXPECTED results:
center [1.25049755875451e-15 3.17478833187820]

Thoughts anyone?

Comment only
26 Apr 2010 SengChoy Liew

SengChoy Liew (view profile)

if i have a gray scale picture with circles of various diameter touching each other. How can i make use of this fitting to evaluate the distribution in histogram of circles in that pictures ?

I could not find any codes in image processing toolbox for it.

Anyone out there for help ?

Comment only
19 Oct 2009 Kristopher White

Kristopher White (view profile)

Great work.

23 Mar 2009 Levi Kilcher

Levi Kilcher (view profile)

This function is fantastic!

Great work, thanks.

23 May 2007 Ali Özgül

Documentation:
+Publis M-files application and theory (5)
+Graphical program interface (5)
+Workable program description (5)
+Works references (3)

Programming:
+M-files complexity (3)
+Mlint error message (5)
+Help problems (4)
+Filename of clashes (5)
+Clone the code line (5)
+Program's symbolic manuscript (3)
+Program samples and runing test (5)
+System requirements (5)
+Program running speed (4)

Files archive:
Seleceted file name (5)
Per files comment (5)
All files usefuly (4)

Thanks:
Very nice application and documentation.
Good works.