# find circle in noisy data

I have an image that looks like the one below:

I'm trying to find the radius (or diameter) of the circle. I have tried using circular Hough transform (via matlab's imfindcircles(bw,[rmin rmax],'ObjectPolarity','bright')) , and by fitting to a circle or an ellipse (home made function that works pretty well for less noisy data, see below).

I've also tried some image processing to get a clearer circle, for example, see below:

se = strel('disk', 2);
bw = imdilate(bw, se);
bw = bwareaopen(bw,100000);
bw =  edge(bw);


However, when I feed the processed image to either techniques (Hough and circle\ellipse fitting) neither of them manage to detect the circle in a decent manner.

Here's a code snippet of the circle finder I wrote (matlab) [row col]=find(bw); contour = bwtraceboundary(bw, row(1), col(1)], 'N', connectivity, num_points);

    x = contour(:,2);
y = contour(:,1);

% solve for parameters a, b, and c in the least-squares sense by
% using the backslash operator
abc = [x y ones(length(x),1)] \ -(x.^2+y.^2);
a = abc(1); b = abc(2); c = abc(3);

% calculate the location of the center and the radius
xc = -a/2;
yc = -b/2;


Alternative approaches will be appreciated...

-
Hough transform looks for a circle, not a filled disc. you would need to do edge detection first to convert the filled disc to an empty circle. what are the properties of your circles? is the size constant? can they be ellipses? can the dots be distributed differently? –  endolith Nov 6 '12 at 22:16
I tried (see edited example), it is either too noisy or, not circular enough? Additionally, the size is constant, and it may have minute ellipticity (though in reality it is a perfectly circular window) due to camera angle error. –  bla Nov 6 '12 at 23:13
if the size and shape is constant, you might try something like cross-correlation of a filled disc template with the original dot image –  endolith Nov 7 '12 at 4:11
Beside my answer, I think that you might be trying to do this in a too later stage of your image processing pipe. Could you tell us more about the problem, and show some prior steps? –  Andrey Nov 7 '12 at 22:23

Here is my solution, it is close to @Yoda's idea, but I changed some steps.

• Mark all pixels such that there are at least 6 pixels in their 7x7 neighborhood
• Remove all blobs, but the largest
• Fill holes
• Apply edge detection
• Find circle using Hough transform

Here is the relevant Matlab code. I am using Hough transform for circles .m file in my code.

function FindCircle()
close all;
im = im(:,:,2);

ims = conv2(double(im), ones(7,7),'same');
imbw = ims>6;
figure;imshow(imbw);title('All pixels that there are at least 6 white pixels in their hood');

props = regionprops(imbw,'Area','PixelIdxList','MajorAxisLength','MinorAxisLength');
[~,indexOfMax] = max([props.Area]);

largestBlobIndexes  = props(indexOfMax).PixelIdxList;
bw = false(size(im));
bw(largestBlobIndexes) = 1;
bw = imfill(bw,'holes');
figure;imshow(bw);title('Leaving only largest blob and filling holes');
figure;imshow(edge(bw));title('Edge detection');

[~,maxIndex] = max(h(:));
[i,j,k] = ind2sub(size(h), maxIndex);
center.x = j;
center.y = i;

title('Final solution (Shown on edge image)');

title('Final solution (Shown on initial image)');

end

-
What is the Hough transform doing here that makes it solve and find the blue circle? Is it projecting many circles of different radii at different positions on the image and finding the one that best fits? –  Mohammad Nov 8 '12 at 2:26
@Mohammad, it is the usual circle detector. It uses binning and voting. –  Andrey Nov 8 '12 at 8:18
You can also use Fast Radial Symmetry Transform (FRST) after the first step in this answer. –  Geniedesalpages Oct 21 '14 at 15:28

It's fairly straightforward to do it using image processing. The following is a proof of concept in Mathematica. You'll have to translate it to MATLAB.

• First, trim the axes and keep only the image part of it. I call this variable img.
• Binarize the image and dilate it, followed by a filling transform. I also remove stray small components that are not connected to the main blob. It should give you something like the following:

filled = Binarize@img ~Dilation~ 3 // FillingTransform // DeleteSmallComponents


• Next, find the centroid of this blob and the equivalent disk radius of the blob (openCV, MATLAB all have equivalent commands to do this)

{center, radius} = 1 /. ComponentMeasurements[filled, {"Centroid", "EquivalentDiskRadius"}]

• That's it! Now plot the original image and a circle with the above center and radius to see how it fits:

Show[img, Graphics[{Red, Circle[center, radius]}]]


-
Awesome answer! Can you please expand on the dilation & filling transform? –  Mohammad Nov 7 '12 at 3:13
@Mohammad Dilation is a basic operation and would be easily explained by the wiki article. Filling transform fills in "holes" or in other words, sets of pixels that are surrounded by pixels of higher value. See the "More information" section here –  Lorem Ipsum Nov 7 '12 at 3:20
Ah sorry, I mis-typed. I am somewhat familiar with dilation transform, I was actually wondering if you could expand on the 'filling transform'. What rule set it is using exactly? I cant seem to find information related to that. Perhaps it goes by another name? –  Mohammad Nov 7 '12 at 3:44
@yoda, thank you for the answer, but if you read the question you'd notice I've tried dilation, and fitting. The image produced before I detect edges is similar to yours. I get some fit, it is not accurate. The same is for your fit, you can see that the top part of the fitted circle is too big, presumably because you take into account the noisy point on the top part above the circle. I have also tried to fit an ellipse (as stated in the question), the problem is that the fit isn't good enough. I think the maybe a better way would be to use the better part of the circle (an arc) to do the fit. –  bla Nov 7 '12 at 5:11
@nate I don't understand what you mean by "top part of the fitted circle" and "better part of the circle". You can use different metrics... bounding box, major axis length, minor axis length, mean distance from centroid, median distance from centroid, etc. It all depends on what you want. –  Lorem Ipsum Nov 7 '12 at 6:33