Finding a circle from a grayscale image in MATLAB
Here is an alternative solution to imfindcircles. Basically a threshold image by extending it with a disk structuring element and then, after finding the edges, apply the Hough transform to define the circle using the circle_hough
file sharing algorithm available here .
Here is the code:
clear
clc
close all
A = imread('CircleIm.jpg');
%// Some pre-processing. Treshold image and dilate it.
B = im2bw(A,.85);
se = strel('disk',2);
C = imdilate(B,se);
D = bwareaopen(C,10000);
%// Here imfill is not necessary but you might find it useful in other situations.
E = imfill(D,'holes');
%// Detect edges
F = edge(E);
%// circle_hough from the File Exchange.
%// This code is based on Andrey answer here:
%https://dsp.stackexchange.com/questions/5930/find-circle-in-noisy-data.
%// Generate range of radii.
radii = 200:10:250;
h = circle_hough(F, radii,'same');
[~,maxIndex] = max(h(:));
[i,j,k] = ind2sub(size(h), maxIndex);
radius = radii(k);
center.x = j;
center.y = i;
%// Generate circle to overlay
N = 200;
theta=linspace(0,2*pi,N);
rho=ones(1,N)*radius;
%Cartesian coordinates
[X,Y] = pol2cart(theta,rho);
figure;
subplot(2,2,1)
imshow(B);
title('Thresholded image (B)','FontSize',16)
subplot(2,2,2)
imshow(E);
title('Filled image (E)','FontSize',16)
subplot(2,2,3)
imshow(F);hold on
plot(center.x-X,center.y-Y,'r-','linewidth',2);
title('Edge image + circle (F)','FontSize',16)
subplot(2,2,4)
imshow(A);hold on
plot(center.x-X,center.y-Y,'r-','linewidth',2);
title('Original image + circle (A)','FontSize',16)
Which gives the following:
You can play with the parameters passed to the threshold, or expand to see how this affects the result.
Hope it helps!
source to share
Here's another approach to solving this problem. It is not based on the Hough transform as imfindcircles and the previous answer.
Basically:
- Image segment (threshold estimated using the Otsu method)
- Find the connected components and leave only 2% of them, starting with the ones whose area is larger.
- Expand the image with a small circle (disk)
- Find the newly connected components and take the component with the largest area as a circle
HT is sometimes slow, depending on input size and resolution. It may be helpful to compare the execution times of both approaches (HT, non HT).
The proposed method can also detect objects of a different shape (non-circular).
function[circle] = ipl_find_circle(I)
% NOTE: I is assumed to be a grayscale image
% Step 1: Segment image using Otsu´s method
t = graythresh(I);
BW = im2bw(I, t);
% Step 2: Leave just "big" components on binary image
[L, num] = bwlabel(BW);
stats = regionprops(L, 'Area', 'PixelIdxList');
area_vector = [stats(:).Area];
area_vector = sort(area_vector);
threshold_pos = floor(num * 0.98);
threshold = area_vector(threshold_pos);
for i=1:num
if(stats(i).Area < threshold)
BW(stats(i).PixelIdxList) = false;
end
end
% Step 3: Dilate image with a circle of small radius
str = strel('disk', 5);
BW = imdilate(BW, str);
% Step 4: Take component with biggest area as the circle
L = bwlabel(BW);
stats = regionprops(L, 'Area', 'BoundingBox', 'Centroid', 'EquivDiameter');
area_vector = [stats(:).Area];
[max_value, max_idx] = max(area_vector);
soi = stats(max_idx);
% Set output variable
circle = imcrop(I, soi.BoundingBox);
% Display results
radius = soi.EquivDiameter/2;
N = 1000;
theta = linspace(0, 2*pi, N);
rho = ones(1, N) * radius;
[X,Y] = pol2cart(theta, rho);
X = soi.Centroid(1) - X;
Y = soi.Centroid(2) - Y;
figure;
subplot(1,2,1);
imshow(I);
hold on;
plot(X, Y, '-r', 'LineWidth', 2);
title('Original graycale image + circle', 'FontSize', 12)
subplot(1,2,2);
imshow(circle);
title('Circle region', 'FontSize', 12);
end
source to share