Recognize centroids and radii of many circle

I have an image that contains many broken/slightly deformed circles because it is a cross section of an object. I need to be able to recognize the centroids and radii of them. I have attached the picture below. I have the image processing toolbox and have tried the method in the link below but it returns no radii or centroids. I have estimated the radius to be between 25 and 75. I also attached my code and picture I used.

<http://blogs.mathworks.com/pick/2008/05/23/detecting-circles-in-an-image/>

img = imread('img.png'); %Reads image from directory
red = img(:,:,2); %Use first (red) color plane
red = adapthisteq(red);  %Use adaptive histogram equalization
tic
[accum,circen,cirrad] = CircularHough_Grd(red,[15 75],20,13,1);
toc

if any(cirrad <= 0) %Eliminate zero radii
    inds = find(cirrad>0);
    cirrad = cirrad(inds);
    circen = circen(inds,:);
end
imshow(red);
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;

[Image Analyst]

I don't see the img.png that you attached. Are you sure you attached it? Did you use the image icon or the paper clip icon?

Answer 1

[Image Analyst]

I don't see your image. Can you use a method like in my Image Segmentation Tutorial : http://www.mathworks.com/matlabcentral/fileexchange/25157-image-segmentation-tutorial-blobsdemo?

[Michael]

Sorry, I posted. This does not go into too much detail on how (specific codes) to recognize the shapes and then find their centroids and radii.

[Image Analyst]

That image looks like it has bad chromatic aberration. Can you get a "blank shot" so that we can do some kind of background correction on it. Myabe you could try this: http://www.mathworks.com/matlabcentral/fileexchange/25619-image-segmentation-using-statistical-region-merging or this http://www.mathworks.com/products/demos/image/color_seg_k/ipexhistology.html

[Michael]

When you mean a "blank shot" is that an image of just the background without the circles in them? If I could get that shot how would you do the background correction? Using the first link I get these modified images which dont seem very useful

.

[Michael]

I was able to get these b&w images using the following code. Does this make it easier to recongize the centroids and radii?

A = imread('img.png');
I = rgb2gray(A);
imshow(I)
I2 = mat2gray(I);
subplot(2,2,1)
imshow(I)
subplot(2,2,2)
imshow(I2)

[Image Analyst]

Yes that looks sharper, I think. It looked like in your original image there was shear between the red, green, and blue color channels so that the centroids might not be in the same place for each color channel. If these centroids are okay, you can look at my BlobsDemo app which does that: http://www.mathworks.com/matlabcentral/answers/contributors/1343420-image-analyst/answers You might have to separate any touching blobs using something like watershed: http://blogs.mathworks.com/steve/2006/06/02/cell-segmentation/

[Michael]

I see how to find the centroids and radii from your BlobsDemo app but seperating the circles is still difficult. The watershed method you used doesnt seem to help with my image. I believe it is because there is no intense maximum of white within the circles. Is there another method of seperating them? I was able to trace around the B&W image (not that that helps because the circles are still connected).

[Image Analyst]

Attach the image as a file (use the paper clip icon). If I get time (I'm traveling in South America at the moment), I'll see what I can do. I think that it should work though if you just adjust some parameters.

[Michael]

I was just wonder, how does the program calculate the radius of each blob since they are not perfect circles?

[Image Analyst]

If a blob is not a perfect circle, you can get the "diameter" by asking regionprops() for "EquivDiameter" which is the diameter you'd get if you molded your blob into a perfect circle with the same number of pixels.

[Michael]

Ok thanks, I was just wondering. Have you been able to adjust any parameters to captures all the white blobs? From the last plot in the program it seems like all the blobs arent being recognized and having their radii and centroids calculated.

[Image Analyst]

No - I just haven't had time, and I'm about to leave town again in a few hours for a few days so I may not get to it soon, if ever. Sorry.

Answer 2

[Michael]

Here is the original image and my result images. Below is also the code I used to generate the centroids and radii. Let me know if there is anything you think I could tweak. Thanks for your help.

a =2; %rows
b =2; %columns
I = imread('test2.png');
I2 = rgb2gray(I);
originalImage = mat2gray(I2);
subplot(a, b, 1);
imshow(originalImage);
I_eq = adapthisteq(originalImage);
subplot(a, b, 2);
imshow(I_eq);
bw = im2bw(I_eq, graythresh(I_eq));
subplot(a, b, 3);
imshow(bw);
BW = imopen(bw, strel('square', 25));
subplot(a, b, 4);
imshow(BW);

%%From Blob
blobMeasurements = regionprops(BW, BW, 'all');   
numberOfBlobs = size(blobMeasurements, 1);

% bwboundaries() returns a cell array, where each cell contains the row/column coordinates for an object in the image.
% Plot the borders of all the coins on the original grayscale image using the coordinates returned by bwboundaries.
subplot(a, a, 4); imagesc(BW);
title('Outlines, from bwboundaries()'); axis square;
hold on;
boundaries = bwboundaries(BW);  
numberOfBoundaries = size(boundaries);
for k = 1 : numberOfBoundaries
  thisBoundary = boundaries{k};
  plot(thisBoundary(:,2), thisBoundary(:,1), 'g', 'LineWidth', 2);
end
hold off;

fontSize = 12;  % Used to control size of "blob number" labels put atop the image.
labelShiftX = -7;  % Used to align the labels in the centers of the coins.
blobECD = zeros(1, numberOfBlobs);
% Print header line in the command window.
fprintf(1,'Blob #      Mean Intensity  Area   Perimeter    Centroid       Diameter\n');
% Loop over all blobs printing their measurements to the command window.
for k = 1 : numberOfBlobs           % Loop through all blobs.
  % Find the mean of each blob.  (R2008a has a better way where you can pass the original image
  % directly into regionprops.  The way below works for all versions including earlier versions.)
    thisBlobsPixels = blobMeasurements(k).PixelIdxList;  % Get list of pixels in current blob.
    meanGL = mean(originalImage(thisBlobsPixels)); % Find mean intensity (in original image!)
  meanGL2008a = blobMeasurements(k).MeanIntensity; % Mean again, but only for version >= R2008a
  
  blobArea = blobMeasurements(k).Area;    % Get area.
  blobPerimeter = blobMeasurements(k).Perimeter;    % Get perimeter.
  blobCentroid = blobMeasurements(k).Centroid;    % Get centroid.
  blobECD(k) = sqrt(4 * blobArea / pi);          % Compute ECD - Equivalent Circular Diameter.
    fprintf(1,'#%2d %17.1f %11.1f %8.1f %8.1f %8.1f % 8.1f\n', k, meanGL, blobArea, blobPerimeter, blobCentroid, blobECD(k));
  % Put the "blob number" labels on the "boundaries" grayscale image.
  text(blobCentroid(1) + labelShiftX, blobCentroid(2), num2str(k), 'FontSize', fontSize, 'FontWeight', 'Bold');
end

MY RESULTS

Blob #      Mean Intensity  Area   Perimeter    Centroid       Diameter
# 1               0.6       169.0     48.0      7.0      7.0     14.7
# 2               0.7     19114.0    972.0     83.2    107.3    156.0
# 3               0.8    316607.0  13616.2    684.7    722.0    634.9
# 4               0.7      1380.0    158.4     12.2    504.2     41.9
# 5               0.7      7980.0    432.1     72.7    863.4    100.8
# 6               0.7      9832.0    506.0     55.7    980.7    111.9
# 7               0.7     37468.0   1681.8    232.2   1080.1    218.4
# 8               0.7      4647.0    266.4     50.4    250.2     76.9
# 9               0.7      4518.0    271.5     47.1    326.3     75.8
#10               0.7      9478.0    520.4     81.1    713.8    109.9
#11               0.7     10015.0    506.3     96.7    553.4    112.9
#12               0.7       937.0    145.1     81.3      8.3     34.5
#13               0.7      4563.0    264.1    128.6    333.5     76.2
#14               0.7     25228.0   1129.5    260.6    190.1    179.2
#15               0.7      2651.0    210.8    184.7     20.0     58.1
#16               0.7     49666.0   2108.4    321.3    882.9    251.5
#17               0.7     22384.0    964.1    416.7    384.9    168.8
#18               0.7      2518.0    222.3    347.7     16.6     56.6
#19               0.7     32146.0   1336.7    486.6    158.3    202.3
#20               0.7      4101.0    263.4    441.8     25.6     72.3
#21               0.7      5796.0    296.4    465.8    319.0     85.9
#22               0.8      3614.0    252.7    513.3   1123.7     67.8
#23               0.7      3795.0    248.9    552.6     30.3     69.5
#24               0.8     56637.0   2349.6    716.8    335.9    268.5
#25               0.8     15424.0    856.8    720.3   1105.5    140.1
#26               0.7      7792.0    428.0    656.4     37.9     99.6
#27               0.8      4549.0    276.0    704.6   1029.5     76.1
#28               0.7     16582.0    740.0    772.8    137.7    145.3
#29               0.8      4463.0    260.6    802.7    997.0     75.4
#30               0.7      2656.0    202.2    802.4     21.1     58.2
#31               0.7      4901.0    280.3    890.6     44.8     79.0
#32               0.8      5040.0    288.4    912.1   1107.8     80.1
#33               0.7     21588.0   1079.0   1011.0    312.4    165.8
#34               0.7      4266.0    255.1    925.5    198.8     73.7
#35               0.8      5091.0    281.1    947.0   1028.4     80.5
#36               0.7      4783.0    271.4    957.2    115.3     78.0
#37               0.7      7007.0    433.4    989.0     28.8     94.5
#38               0.8      5310.0    292.3   1003.7    407.1     82.2
#39               0.8      4529.0    265.4   1018.8   1114.8     75.9
#40               0.7      5055.0    277.7   1019.1    180.5     80.2
#41               0.8      4172.0    250.3   1028.7   1000.3     72.9
#42               0.7      5035.0    275.4   1071.3    523.9     80.1
#43               0.7      5025.0    276.3   1075.8     85.6     80.0
#44               0.8      3607.0    247.2   1118.3   1119.4     67.8
#45               0.7     11884.0    759.7   1169.4    484.4    123.0
#46               0.7     12108.0    683.1   1169.1    229.5    124.2
#47               0.7      4426.0    254.5   1159.5     65.8     75.1
#48               0.8      3580.0    238.0   1175.3    832.0     67.5
EDU>>