File Exchange Pick of the Week

Detecting Circles in an Image 52

Posted by Brett Shoelson,

A customer recently provided me with an image of cells that were roughly circular, but not very well defined, and often overlapping. He asked how we might use MATLAB and the Image Processing Toolbox to segment the cells in the presence of noise.

Many of us know the Hough transform functionality in the Image Processing Toolbox, and the ability of that function to detect lines in an image. With some modifications, the Hough transform can be used to find other shapes as well. In this case, I wanted to find circles; a quick search for "detect circles" on the MATLAB Central File Exchange revealed Tao Peng's implementation of the Circular Hough Transform. After an easy download, I was on my way to solving the problem. Tao's file is Brett's choice for this week's Pick of the Week.

(Thanks to Cem Girit at Biopticon Corp. for permission to use his image!)

Contents

Read, Display image

togglefig Original
% Note: togglefig is a utility function I wrote for managing figure
% windows, and is available
% <http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do?objectId=18220&objectType=file here>.
img = imread('NB1ln1.jpg');
imshow(img);

How would you segment this?

Segmentation is often a very challenging task, particularly with noisy images like this one. One might try simple threshold or edge detection routines, or more sophisticated watershed approaches, for instance. Unfortunately, none of these approaches was giving me the results I wanted for this particular image. I was at a bit of a loss--until I found CircularHough_Grd. Tao suggests some filtering operations in his very helpful comments. For this demonstration, I'm going to see how the algorithm can perform with just some very minor pre-processing.

Discarding color information

Tao's function works directly on grayscale images. Rather than converting the color image to grayscale with the Image Processing Toolbox's RGB2GRAY function, I elected to simply use the first (red) color plane, and to use adaptive histogram equalization:

togglefig('Red Plane, Adjusted')
red = img(:,:,2);
red = adapthisteq(red);
imshow(red)

Parameters for the segmentation

Before I segment, I needed to know how big the cells were in the image; CircularHough_Grd takes as an input a range of radii to search for. Using the IMDISTLINE function in the Image Processing Toolbox's IMTOOL, I estimated that the radii of interest range from about 5 to 25 pixels.

You can play around with the other input arguments to modify the function's sensitivity to less-than-perfect circles, for instance, or to change the way it deals with concentric circles. This makes the program pretty flexible!

Now for the segmentation...

tic;
[accum, circen, cirrad] = ...
    CircularHough_Grd(red, [5 25],...
    20, 13, 1);
toc
Elapsed time is 4.169272 seconds.

...and a bit of post-processing

Note to Tao: occasionally, your algorithm returns zero-radii "hits":

any(cirrad <= 0)
ans =
     1

This is easy to address (for instance, to keep the RECTANGLE command below from erroring), but might be an opportunity for enhancement.

if any(cirrad <= 0)
    inds = find(cirrad>0);
    cirrad = cirrad(inds);
    circen = circen(inds,:);
end

View the results

Now let's see how well the algorithm performed:

togglefig Results
imshow(img);
hold on;
plot(circen(:,1), circen(:,2), 'r+');
for ii = 1 : size(circen, 1)
    rectangle('Position',[circen(ii,1) - cirrad(ii), circen(ii,2) - cirrad(ii), 2*cirrad(ii), 2*cirrad(ii)],...
        'Curvature', [1,1], 'edgecolor', 'b', 'linewidth', 1.5);
end
hold off;

That's pretty remarkable, especially given the simplicity of my pre-processing. (Adaptive histogram equilization helped a lot; Tao's suggested filters improve the performance even more.)

Final comments

The time it takes to run this algorithm varies markedly, depending on the user settings. In this case, it took my computer approximately 4 seconds--but did a pretty amazing job of segmentating the image. Note how well it handled overlapping cells (circles), for instance:

togglefig Blowup
imshow(img)
xlims = [406 520];
ylims = [52 143];
set(gca,'xlim',xlims,'ylim',ylims)
inImageCircles = find(inpolygon(circen(:,1), circen(:,2), xlims([1 2 2 1]), ylims([1 1 2 2])));
for ii = 1 : numel(inImageCircles)
    rectangle('Position',...
        [circen(inImageCircles(ii),1) - cirrad(inImageCircles(ii)),...
        circen(inImageCircles(ii),2) - cirrad(inImageCircles(ii)),...
        2*cirrad(inImageCircles(ii)),...
        2*cirrad(inImageCircles(ii))],...
        'Curvature', [1,1], 'edgecolor', 'b', 'linewidth', 1.5);
end

If you do any image processing, you might recognize how Tao's function made a very challenging problem pretty manageable. If you're not an image processing expert, this might seem a bit arcane to you...but hopefully you'll find some value in it anyway. I'd love to hear your comments!


Get the MATLAB code

Published with MATLAB® 7.6

52 CommentsOldest to Newest

Hi Brett,

very nice post, with mentioning some useful tools like your togglefig, -
similar to a function I also wrote for myself :-) – adapthisteq and
imdistline.

However, in the end you say that the result is very good and overlapping
cells are handled well. However, in the image there are still at least
twenty or thirty cases where overlapping cells melt into a single detection.
Depending on the level of accuracy demanded by the application, a more
sophisticated processing might be needed.

Do you think a classification approach might do the trick? One could crop
the area inside and under the circles that you have plotted in blue and run
a classifier to decide if there is a single cell or two overlapping cells.
The objects are not that complicated, so a moderate training data base might
be sufficient. Which classification tools are provided by Matlab for such a
task?

Yours
Markus

Hi,

Doesn’t
red = img(:,:,2);
extract the green chanell of an RGB image (RGB=123)?

Looks like a usefull pick!

Best Regards,
Flemming

@Markus:
Regarding the specific cases of failure of the algorithm on this image–it is possible that there are different parameter levels that might give a better segmentation. Or, if you can show a different approach that gives better results, I would absolutely love to see it!

Regarding the classification of cells: right you are. In fact, the second half of this customer’s inquiry dealt with more sophisticated classification/clustering of detected circles. The MathWorks provides lots of tools for clustering, including some (kmeans, for instance) in the Statistics Toolbox, and the Neural Networks and Fuzzy Logic Toolboxes. My approach was lower-level for this problem; I simply looped over all of the circles, compiling information on the intensities and distributions of encircled pixels, and implemented some crude filters.

@Flemming:
Of course you’re right, Flemming. Consider this a typo! I went back between using red and green, and ended up with a mistake in my post. Thanks for point it out. For the record, given an RGB image: red = img(:,:,1), and green = img(:,:,2);

@Brett:
I am considering to make a small “home project” of this task, but I don’t know if I will find the time. However, I guess if I could come up with a better approach, it would be much more complicated than your solution in the post!

My first step would be to use a HSV color representation and extract the “pink” channel instead of the red channel.

By the way, are there more images like that available? I mean for training some sort of a classificator.

Markus

for serious work, use cellprofiler.org. it’s open source, written in matlab (and some java), and is capable of automatic classification of cells by types, or by selected examples (training). don’t bother writing your own classification software from scratch…

@User, Markus: I wasn’t familiar with cellprofiler. Thanks for sharing that link–nice to know it exists. (If you or others want to share their experiences with that tool, I’d be interested in hearing them.)

That said, I agree with Markus: while we don’t want to _have_ to recreate the wheel, sometimes doing so is fun. Markus, if you come up with a better (even if more complex) solution, please share it on the File Exchange!

Unfortunately, I can’t share this image; I received permission to use it in demos, but not to distribute it. You can likely find several useful images with a Google search. And the Image Processing Toolbox ships with a few that might be useful. (These come to mind: coins.png, blobs.png, circles.png, pears.png, pillsetc.png, snowflakes.png, and AT3_1M4_01.tif through AT3_1M4_10.tif.)

Hello Brett,

I have minimal experience with image analysis; however, I would like to know the reasoning behind your decision to use the red color plane and adaptive histogram equalization rather than the RGB2GRAY function. Thanks.

Michael

Hi Michael,

Whenever I start analyzing a color (rgb) image, the first thing I do is look at the color channels individually, alongside the RGB2GRAY output, to decide what information I want to work with. Often I find that RGB2GRAY gives the most visually appealing grayscale representation of a color image, but doesn’t necessarily provide the highest SNR for the subsequent analysis. In fact, I wrote (and recently posted to the File Exchange: http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do?objectId=19706&objectType=FILE) a utility function called ExploreRGB that facilitates side-by-side comparison of these different views of an RGB image. (My newer version, not yet posted, also includes HSV, L*a*b*, chrominance, etc.)

In the “Detecting Circles” case, I looked at all options and decided that I had the most information in the red or green planes, and proceeded from there.

Thanks for the question!
Brett

Hi Brett,

I am trying to run this program on simulated images generated in MATLAB but I don’t know what kind of filtering I am supposed to do on the image. I am drawing circles parametrically in one program and then I want to analyze the image that is generated to see if I can re-extract the locations of the circles.

What kind of filtering is necessary because I am getting some strange circle locations. For example, one time I ran the image through and there were hundreds of circles, all on the perimeter of the actual circle.

Thanks.
David

Hi David,
It’s hard to answer in the abstract…kind of like debugging your code without even seeing it. I do think that Tao’s algorithm has a few quirks to it, and that you might get better results just by changing the input parameter values. As far as filtering…Tao gives some good examples in which he discusses some filtering (a bit); I wonder if his suggested filter will help you get better results. Good luck!
Cheers,
Brett

That said, I agree with Markus: while we don’t want to _have_ to recreate the wheel, sometimes doing so is fun. Markus, if you come up with a better (even if more complex) solution, please share it on the File Exchange!

Hi Brett,
i was trying ur method of detecting circle but it seems like it nt working.. there hav been lots of varible undeclared n some functions cant be use like:
[accum, circen, cirrad] = …
CircularHough_Grd(red, [5 25],…
20, 13, 1);
do i need to declare or do something with the accum, circen, cirrad n the CircularHough_Grd??

*if it possible, can u send me the whole program code..
Thx in adv,
Andy

Hi Andy,
accum, circen, and cirrad are OUTPUTS of Tao’s (not my) program. You don’t need to declare them; function CircularHough_Grd does that for you. Tao’s program is “whole,” as far as I can tell. (And you download that from the File Exchange.)
Regards,
Brett

Hi Brett,i ran ur program..but im getting an error that
Function ‘-’ is not defined for values of class ‘single’..
pls do help

@Joel: Hi Joel, I was able to run Tao’s (not my) program. I’m afraid I can’t help you debug your usage of it. From the sound of the error message, I would suspect that something got corrupted in your copy of the mfile. You might try a fresh download.
Brett

@Joel – What version of MATLAB are you running? The MINUS function should have a single method, UNLESS you are running a really old copy of MATLAB. Non-double arithmetic was added to MATLAB in R14 (around 2004).

Brett,

I am using Tao’s code to detect circular shapes in Red blood cell images and its working quite well. I am attempting to add the diameter results to my figure, in simple text form.

I would like it to be interactive, so the user clicks on a circle and the desired diameter value appears. Is this doable?

Thanks in advance for your help
Lee

Hi Lee,
Of course! With MATLAB, all things are possible. :)
Here’s one approach using a nested function:

function ShowRectPos
axes(‘units’,'normalized’);
r(1) = rectangle(‘pos’,[0.2,0.3,0.2,0.2],…
‘facecolor’,'w’);
r(2) = rectangle(‘pos’,[0.5,0.6 0.3,0.3],…
‘facecolor’,'r’,'curvature’,[1 1]);

set(r,’buttondownfcn’,@showPos)
xlim([0 1]);ylim([0 1]);
axis equal

function showPos(obj, varargin)
delete(findobj(‘tag’, ‘tmptxt’));
tmp = get(obj,’pos’);
str = sprintf(‘Position = [%0.1f, %0.1f, %0.1f, %0.1f]‘,get(obj,’pos’));
title(str)
text(tmp(1),tmp(2),str,’tag’,'tmptxt’)
end

end

Cheers,
Brett

Thanks for the reply Brett.

Matlab really knows all doesn’t it?!

Quick question just, do I add this code in my command window, following my plotting code or do I add as an M File?

Thank you!
Lee

Lee, if you were to incorporate the original code I posted to show off Tao’s file into a function, and then add the SHOWPOS code block as a subfunction therein, you should be set. All you’d then have to do is set the BUTTONDOWNFCN property of any rectangle (circle) to call it. For instance, if R is a vector of handles to your rectangles:

set(r,’buttondownfcn’,@showPos)

Cheers,
Brett

Sorry Brett but im not sure what you mean. I have the following code to date, simply entered in my command window:

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Load image
rawimg = imread(‘clll40_005.tif’);
rawimg = rawimg(1:884,:);

%% Image binning
st_imgbin = immBinning( uint16(rawimg), [1,1] );
imgshrunk = uint8( st_imgbin.img / 4 );

%% Image smoothing using Gaussian filter
fltr4img = fspecial( ‘gaussian’, 5, 1.2 );
imgfltrd = imfilter( imgshrunk, fltr4img, ‘replicate’ );

%% Detection of circular features by Hough Transform
tic;
[accum, circen, cirrad] = …
CircularHough_Grd_SL( imgfltrd, [10 190], 6, 25, 0.6 );
toc;

%% Results display
figure(1); imagesc(accum); axis image;
title(‘Accumulation Array from Circular Hough Transform’);

figure(2); imagesc(imgfltrd); colormap(‘gray’); axis image;
hold on;
plot(circen(:,1), circen(:,2), ‘r+’);
for k = 1 : size(circen, 1),
DrawCircle(circen(k,1), circen(k,2), cirrad(k), 32, ‘b-’);
end
hold off;
title(['Raw Image with Circles Detected ', ...
'(center positions and radii marked)']);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Do I convert this into an m-file? and then add SHOWPROS within? Or do I add showpros as a separate m file?

When you say the SHOWPROS code block, is that everything you have added in the comment above?

And how do I go about knowing what to set r as? is r my required circle?

Overall, sorry im so confused. I’ve very little experience in calling functions to date so this seems to be going over my head so far.

Thanks you for any help you can give me.
Lee

Lee, Notice that in my original code, I didn’t use DrawCircle; I used the RECTANGLE command instead. I later modified my code to return handles to each rectangle so I could show you how to modify their BUTTONDOWNFCNs. If you want to use DRAWCIRCLE, you’ll have to modify it so that it returns a vector of handles to the plots it creates, then similary modify the buttondownfcns of those.
Regards,
Brett

Thank you for the help Brett. This is now working nicely for me. I have discovered a further issue however upon implementation of some alternative images I have been given. I have a collection of low res SEM images, similar to your example. Circle detection is not working however. I have altered Taos parameters significantly but to no avail. Any ideas?

See below
Original Image:
http://i44.tinypic.com/2liu3x5.jpg
Result: [Somewhat redundant:(]
http://i42.tinypic.com/2ik4i1c.jpg

Lee, I’m afraid I don’t have any better suggestions for you…except that maybe you can try posting your question to the newsgroup for some additional support. Good luck!
Brett

It seems that the function that Tao wrote has a bit of a minor problem, there is a function that is called filtfilt that is not defined. what would u recommend to overcome this problem?

Hi Suwei,
Thanks for pointing this out. Actually, FILTFILT is a function in the Signal Processing Toolbox that does zero-phase digital filtering. So Tao’s function should specify that as a pre-req. A workaround — besides buying the Signal Processing Toolbox ;) — would entail tweaking the code to avoid the dependency.

Cheers,
Brett

Hi Brett, I have tried Tao’s function, however I get around 25 bogus circles when I run with the default parameters (CircularHough_Grd(img, radrange, grdthres, fltr4LM_R, multirad, fltr4accum)) on the image link below.

Do you reckon it would work for this image: http://img15.imageshack.us/i/f142.png/ ?

If so, which input parameter shall I tweak?

Would be great if you could point me in the right direction.

Sincerely,
berkan

Hi Brett,

I’ve tested Tao’s program and yours too. They ran successfully. I’ve tried to applied your method above and Tao’s to detect sun features in an image that comes with text. Some of the suns are detected and some are not. Other features in the image were also detected such as the text.

I thought of removing the text so that it doesn’t interfere the program in detecting the circles in the image but not really sure how. Do you have any ideas how to start?

Many thanks!
Nickey

Sir,
I am trying to detect tires of a vehicle using Tao’s program. Before giving image to his program , I have done some preprocessing such as adaptive hist equalization, thresholding,dilation, erosion and edge detection on image. This was done in order to reduce false detection. But problem is that it is detecting circles in images at some places where no circle are there. Please suggest some method to solve the problem.

Nishtha,
I’m sorry I don’t have specific suggestions for how to improve your code. I’d play with the different input parameters to Tao’s code to see if you can improve your circle detection.
You might also consider posting a snippet of code to CSSM, along with a discussion of your results, and of your _desired_ results.
Brett

Hi Brett, you code looks kind of what I need for my “Coin Recognition” project!!!

Well I have an image which has a coin in it. This coin occupies no more than 20% (this is my region of interest, ROI) of the entire image; so what I need is:

• Firstly, the circle detection (I need you code here),
• Secondly, find out the radius of the circle (I have figured out how to do it already), and
• Thirdly, (x,y) coordinates detection (suggestions here are welcomed).

I am trying to make you code run but I get the following error message:

??? Undefined command/function ‘CircularHough_Grd’.

Error in ==> circdet at 14
[accum, circen, cirrad] = CircularHough_Grd(red, [5 25],20, 13, 1);

This error is driving me crazy; I don’t really know what to do next. I have installed all the toolboxes my Matlab (v 7.1.0.246) has. Please help me debugging it.

Thanks in advance for the help you can provide me with,
Luis

@Luis:
Hi Luis,
First, note that this is not my code…the circle detection code was Tao’s; I simply blogged about it herein.
Second, my guess–without seeing your image–is that the Image Processing Toolbox function REGIONPROPS will probably do just about everything you need for your project. Segment the objects in your image, then use REGIONPROPS to calculate areas, eccentricities, and centroids. Filter the image (read the doc for REGIONPROPS!) to find the object with the appropriate eccentricity and area. Then find the centroid of that object.
Cheers,
Brett

Thanks Brett, Im gonna do that during this week and if I have any question I will turn to you again.

Nice blog ;) 10/10

hello
i tested the code but it didnt work
some of the functions are missed like CircularHough_Grd
how can i have these functions?
thanks

Excellent post. Quick question: is there a particular reason why this page cannot be saved as a web archive? It’s the first Matlab blog post ever that gives me such an error. Thank you.

@arni:
Hi arni,
Can you clarify your question for me? What steps, exactly, are you trying that typically work but are failing here?
Thanks,
Brett

@luis:
CircularHough_grd is the Pick of the Week file, on which this blog post is focused. You can get it by following the links in the above verbiage.

I consistently get the following error:

??? Undefined function or method ‘filtfilt’ for input arguments of type ‘double’.

Error in ==> CircularHough_Grd at 606
SgnCv = filtfilt( fltr4SgnCv , [1] , SgnCv );
I run latest R2010b. Any ideas?

Thanks, Ev.

@Ev: filtfilt is a function from our Signal Processing Toolbox. Tao should have indicated that as a requirement–thanks for catching that.
Cheers,
Brett

i have downloaded a file to detect the circle and got the error.I have signal processing toolbox installed too but to no avail.Would you please help?

***************************************
Error in ==> detectCircle at 27
[accum, circen, cirrad] = CircularHough_Grd(red,[15 60]);
***************************************


img = imread('cell4.jpg');
imshow(img);

%% Discarding color information
% Tao's function works directly on grayscale images. Rather than converting
% the color image to grayscale with the Image Processing Toolbox's RGB2GRAY
% function, I elected to  simply use the first (red) color plane, and to
% use adaptive histogram equalization:
%togglefig('Red Plane, Adjusted')
red = img(:,:,2);
red = adapthisteq(red);
imtool(red)

%% Parameters for the segmentation
% Before I segment, I needed to know how big the cells were in the image;
% CircularHough_Grd takes as an input a range of radii to search for. Using
% the IMDISTLINE function in the Image Processing Toolbox's IMTOOL, I
% estimated that the radii of interest range from about 5 to 25 pixels. 
%
% You can play around with the other input arguments to modify the
% function's sensitivity to less-than-perfect circles, for instance, or to
% change the way it deals with concentric circles. This makes the program
% pretty flexible!

%% Now for the segmentation...
tic;
[accum, circen, cirrad] = CircularHough_Grd(red,[15 60]);
toc

%% ...and a bit of post-processing
% Note to Tao: occasionally, your algorithm returns zero-radii "hits":
any(cirrad <= 0)
%%
% This is easy to address (for instance, to keep the RECTANGLE command below
% from erroring), but might be an opportunity for enhancement.
if any(cirrad 0);
    cirrad = cirrad(inds);
    circen = circen(inds,:);
end

%% View the results
% Now let's see how well the algorithm performed: 
togglefig Results
imshow(img);
hold on;
plot(circen(:,1), circen(:,2), 'r+');
for ii = 1 : size(circen, 1)
    rectangle('Position',[circen(ii,1) - cirrad(ii), circen(ii,2) - cirrad(ii), 2*cirrad(ii), 2*cirrad(ii)],...
        'Curvature', [1,1], 'edgecolor', 'b', 'linewidth', 1.5);
end
hold off;
%%
% That's pretty remarkable, especially given the simplicity of my
% pre-processing. (Adaptive histogram equilization helped a lot; Tao's
% suggested filters improve the performance even more.)

%% Final comments
% The time it takes to run this algorithm varies markedly, depending on the
% user settings. In this case, it took my computer approximately 4
% seconds--but did a pretty amazing job of segmentating the image. Note how
% well it handled overlapping cells (circles), for instance:
togglefig Blowup
imshow(img)
xlims = [406 520];
ylims = [52 143];
set(gca,'xlim',xlims,'ylim',ylims)
inImageCircles = find(inpolygon(circen(:,1), circen(:,2), xlims([1 2 2 1]), ylims([1 1 2 2])));
for ii = 1 : numel(inImageCircles)
    rectangle('Position',...
        [circen(inImageCircles(ii),1) - cirrad(inImageCircles(ii)),...
        circen(inImageCircles(ii),2) - cirrad(inImageCircles(ii)),...
        2*cirrad(inImageCircles(ii)),...
        2*cirrad(inImageCircles(ii))],...
        'Curvature', [1,1], 'edgecolor', 'b', 'linewidth', 1.5);
end

Hi brett
I have a problem in opening mathwork site and can not pick togglefig function
pleas Email it to me
Thank you

Respected sir,
I am a M.Tech Student at Cochin University of Science & Technology,
Kerala. As a part of M.tech programme, i have to do a project. I select
the topic ” image Edge Detection Using Neural networks” . This project mainly used for underwater images.I wish to discuss the topic with
you through email.I kindly request you to help me doing this project.

Yours faithfully

NISHAD. V.H
M.TECH 3RD SEMESTER
DEPT. OF ELECTRONICS
CUSAT
COCHIN, KERALA
nishad.haneefa@cusat.ac.in
+91-9486107096

hai brett.. when i run that code,is show error..why

??? Undefined function or method ‘CircularHough_Grd’ for input arguments of type ‘uint8′.

Error in ==> code4 at 61
[accum, circen, cirrad] = CircularHough_Grd(red,[15 60]);

dear sir

i found the immBinning Function above .But i couldnt find the detail code of the above mentioned function. I am in need od the similar type of function to binn the Image Data.Would you please help me by providing the details of the above function

Sir, I am working on project in which i have to count the number of silkworm eggs which are near to circular in shape.Can u suggest me the diffrent ways to segment it and count the eggs.Here i cant upload that image.

hai..

What are the number 5,25,20,13,1 refer too. How to adjust the number when different picture been used.

CircularHough_Grd(red, [5 25], 20, 13, 1);

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