File Exchange Pick of the Week

Fitting a Circle, easily. 20

Posted by Brett Shoelson,

You may have noticed some recent changes in the format of this blog. Here’s what to expect on a regular basis – two topics per week.

  • On Tuesdays Doug will provide MATLAB tutorials.
  • On Fridays guest bloggers Jiro, Brett and Bob will highlight File Exchange submissions.

A file need not be long to be useful. Brett's pick this week, Izhak Bucher's Circle Fit, is only 5 lines long, excluding comments. But I really like Izhak's entry, and have had many opportunities to use it in the several years since I downloaded it. Somehow, the requirement of fitting a circle to some points seems to occur with puzzling frequency in my work. Izhak's CIRCFIT makes it easy, and is quite robust in most cases.

So how do you use Izhak's function? Consider these noisy data created for illustrative purposes:

R = 10; x_c = 5; y_c = 8;
thetas = 0:pi/64:pi;
xs = x_c + R*cos(thetas);
ys = y_c + R*sin(thetas);

% Now add some random noise to make the problem a bit more challenging:
mult = 0.5;
xs = xs+mult*randn(size(xs));
ys = ys+mult*randn(size(ys));

Plot the points...

figure
plot(xs,ys,'b.')
axis equal

...and then calculate and display the best-fit circle

[xfit,yfit,Rfit] = circfit(xs,ys);
figure
plot(xs,ys,'b.')
hold on
rectangle('position',[xfit-Rfit,yfit-Rfit,Rfit*2,Rfit*2],...
    'curvature',[1,1],'linestyle','-','edgecolor','r');
title(sprintf('Best fit: R = %0.1f; Ctr = (%0.1f,%0.1f)',...
    Rfit,xfit,yfit));
plot(xfit,yfit,'g.')
xlim([xfit-Rfit-2,xfit+Rfit+2])
ylim([yfit-Rfit-2,yfit+Rfit+2])
axis equal

By the way, if it seems odd to you to use a RECTANGLE command to plot a circle, stay tuned for next week's Pick of the Week!

What other data-fitting problems do you commonly come up against, and how do you solve them?


Get the MATLAB code

Published with MATLAB® 7.6

20 CommentsOldest to Newest

Why not point to a modern version of the circle fitting algorithm #19094 in the FEX:

http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do?objectId=19094&objectType=file

Same algorithm as the code you mention, but with real help text, a documented algorithm, error checking, and data scaling for better accuracy.

Reply to ko: fitting an oval is fitting an ellipse. See #7012 in the File Exchange:

http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do?objectId=7012&objectType=FILE

Ko:
Assuming you were interested in fitting an ellipse (which is well defined in mathematical terms), try searching the File Exchange for “fit ellipse”; you’ll find lots of good code that will get you there:

http://www.mathworks.com/matlabcentral/fileexchange/loadFileList.do?objectType=search&criteria=fit%20ellipse&search_submit=fileexchange&query=fit%20ellipse

Duane:
If you search the File Exchange for “fit circle,” you get lots of options to chose from, including your own:

http://www.mathworks.com/matlabcentral/fileexchange/loadFileList.do?objectType=search&criteria=Fit%20circle

Izhak’s code may not be a paragon of beauty, but it’s well-vetted and robust, and it works. Plus, note from the date on Izhak’s submission that it’s been shared for years; I like the idea of recognizing the efforts of one whose code has been widely distributed, but who hasn’t been recognized in this forum yet.

This circle fit has been first published by P.Delogne and I.Kasa in the 1970s and is known as “Kasa method” in statistics. It works well when points cover a large part of the circle but is heavily biased when points are restricted to a small arc. Better fits were proposed by Pratt and Taubin. See my website for respective MATLAB codes.

What if you only have the points (x and y values) of an arc and then you have to find a best fit circle to ultimately determine the radius?

Bob, your circle fit solution looks very similar to Izhak’s. Am I missing something?

Nadine, your question is a bit vague; Izhak’s code (and several other submissions on the File Exchange) are all designed to do exactly what you’re after: fit a circle–including determining a radius–for points on an arc. If you are specifically looking for a solution optimized for the case when all points fall on or near a small arc of a larger circle, I’d suggest taking a look at Nikolai’s “Circle Fit (Pratt Method).”

Cheers,
Brett

Doh! My CIRCLE_FIT uses exactly the same solution as Izhak’s CIRCFIT. Had I looked under the hood at his implementation I would have seen that. Thanks to Brett for the (tactful) catch. For other (different) solutions be sure to read the review comments for Izhak’s submission. Also, Yuri wondered about applications for which arc/circle fitting is important. That submission I linked (6) is a soup-to-nuts example that not only provides such context, it also takes radial variation into account in the final result with measurement error.

question: Is there any way of calculating your circles diameter once it has been drawn? I assume regionprops cannot be used as this a plot graph and not a figure? Can our graph be converted to a figure maybe?

Thank you,
Lee

Bret,
My question was indeed vague, but Bob’s suggestion(fitcircle.m)solved my problem perfectly.
Thanks.

Lee, If I recall my geometry, one would calculate the diameter of the best-fit circle by doubling the best-fit radius. :) (Sorry…couldn’t resist.)

Do note that outputs from this (and, in general, other circle-fitting routines) include the radius of the circle, and the x-y- coordinates of its center.

Once you’ve plotted points in a MATLAB axes, you can save the axes in any of several different formats. Activate the figure, then look at the SAVE AS options from the FILE menu. (These formats include FIG, JPEG, BMP, EPS, PDF, ….)

If you are working from a PLOT in MATLAB, you can always recover the x- and y- data by examining the properties of the children of the axes. In the code above, for instance, you should find 3 children of the main axes. (a = get(gca,’children’). Examine them one by one (get(a(1)), get(a(2)), get(a(3))), and you’ll see that one gives the properties of the rectangle, one gives the properties of the line object that is the center marker, and one gives the properties of the line object (collection of points) that are the plotted x’s and y’s. Among those properties are XDATA and YDATA, from which you can recover the plotted values.

If your question is, if you only have a JPEG (for instance) of a circle, how do you calculate its diameter…that is indeed a different problem–and one for image processing approaches. First, you would have to determine a scale somehow. Then yes, you would likely want to use REGIONPROPS…perhaps after a morphological fill operation?

“If I recall my geometry, one would calculate the diameter of the best-fit circle by doubling the best-fit radius” – now thats just cheeky. funny though:)

Ah, if this outputs radius values then I shouldn’t have a problem. great!

But before i get that far, I’m attempting to determine diameter values from 2 overlapping circles. I am trying to use BWBOUNDARIES to extract arc values from the areas of the circles that are not overlapping. I then hope to use this circlefit.m to give a diameter estimate of each circle.

Can you help me with extracting the required values from bwboundaries to insert as my input to circlefit.m, if you get my drift? bwboundaries just outputs an array and not the raw data if im correct?

thank you:)

lee

BWBOUNDARIES operates on a binary image. If possible, I would avoid the route of converting your raw data to an image file if the only reason for doing so would be extracting and fitting the data for overlapping circles in preparation for re-plotting them.

I haven’t seen your data, but my initial thought is use some other (clustering?) approach on the raw data to associate points with one of the two circles, then fitting them independently.

If you do need to extract xs and ys from the output of BWBOUNDARIES (which returns a cell array of x-y- data), you can easily do so. For instance, the boundaries of the first object in the bw image IMG can be plotted with:

b = bwboundaries(img);
figure;
plot(b{1}(:,1),b{1}(:,2),’r.’);

Cheers,
Brett

CIRCLES
9 circles .
Each Circle lettered, A through I, has it’s own number value from 1 through 9.Sums of circles that overlap at that point.
for example… 5 is the sum of circles A and C. What are A and C’s Values?
How do you arrive at the values?

Quinno,
In a “Hough” world, one would have to parameterize the rectangle or polynomial and create the Hough implementation himself or herself. But the MATLAB function POLYFIT might be useful for you as well.

Hi
I have the same problem as nadine had,but the link which Bob was given is now corrupt & I can’t find a good answer.
would you help me please?!
tnx

Motahhare,
I don’t know why that link is no longer valid; perhaps the author pulled the file. However, if you search the File Exchange for “fit circle,” you should find several options that will do what you need.

@WDM:
I’m not sure I understand what your question is. In the example I used above, I only had several points on (or near) the arc of a circle. I used those points to calculate the best-fit circle through those arc points.
Brett

These postings are the author's and don't necessarily represent the opinions of MathWorks.