Clear Filters
Clear Filters

Separating Objects in an Image

15 views (last 30 days)
Alyssa on 27 Feb 2011
I am trying to perform image analysis on an image where the objects are touching. I am using the watershed function but I am still having issues with it not separating all the objects. Can someone please help me try to separate the objects? I have attached the code and image.
How do you measure the sphericity of each individual object in this image?
clear all;
%Step 1 Image
%Step 2: Use the Gradient Magnitude as the Segmentation Function
hy = fspecial('sobel');
hx = hy';
Iy = imfilter(double(I), hy, 'replicate');
Ix = imfilter(double(I), hx, 'replicate');
gradmag = sqrt(Ix.^2 + Iy.^2);
figure, imshow(gradmag,[]), title('Gradient magnitude (gradmag)')
L = watershed(gradmag);
Lrgb = label2rgb(L);
figure, imshow(Lrgb), title('Watershed transform of gradient magnitude (Lrgb)')
%Step 3: Mark the Foreground Objects
se = strel('disk', 20);
Io = imopen(I, se);
figure, imshow(Io), title('Opening (Io)')
Ie = imerode(I, se);
Iobr = imreconstruct(Ie, I);
figure, imshow(Iobr), title('Opening-by-reconstruction (Iobr)')
Ioc = imclose(Io, se);
figure, imshow(Ioc), title('Opening-closing (Ioc)')
Iobrd = imdilate(Iobr, se);
Iobrcbr = imreconstruct(imcomplement(Iobrd), imcomplement(Iobr));
Iobrcbr = imcomplement(Iobrcbr);
figure, imshow(Iobrcbr), title('Opening-closing by reconstruction (Iobrcbr)')
fgm = imregionalmax(Iobrcbr);
figure, imshow(fgm), title('Regional maxima of opening-closing by reconstruction (fgm)')
I2 = I;
I2(fgm) = 255;
figure, imshow(I2), title('Regional maxima superimposed on original image (I2)')
se2 = strel(ones(5,5));
fgm2 = imclose(fgm, se2);
fgm3 = imerode(fgm2, se2);
fgm4 = bwareaopen(fgm3, 20);
I3 = I;
I3(fgm4) = 255;
figure, imshow(I3)
title('Modified regional maxima superimposed on original image (fgm4)')
%Step 4: Compute Background Markers
bw = im2bw(Iobrcbr, graythresh(Iobrcbr));
figure, imshow(bw), title('Thresholded opening-closing by reconstruction (bw)')
D = bwdist(bw);
DL = watershed(D);
bgm = DL == 0;
figure, imshow(bgm), title('Watershed ridge lines (bgm)')
%Step 5: Compute the Watershed Transform of the Segmentation Function.
gradmag2 = imimposemin(gradmag, bgm | fgm4);
L = watershed(gradmag2);
%Step 6: Visualize the Result
I4 = I;
I4(imdilate(L == 0, ones(3, 3)) | bgm | fgm4) = 255;
figure, imshow(I4)
title('Markers and object boundaries superimposed on original image (I4)')
Lrgb = label2rgb(L, 'jet', 'w', 'shuffle');
figure, imshow(Lrgb)
title('Colored watershed label matrix (Lrgb)')
figure, imshow(I), hold on
himage = imshow(Lrgb);
set(himage, 'AlphaData', 0.3);
title('Lrgb superimposed transparently on original image')

Answers (4)

Brett Shoelson
Brett Shoelson on 1 Mar 2011
I agree with Walter about the depth ambiguity, but I think you still might be able to get reasonably good results with this image. Try something along these lines:
grains = imread('c:\MATLAB_Answers_stuff\grains.png');
grains = rgb2gray(grains);
grainsBW = im2bw(grains,126/255);
revim = imcomplement(grainsBW);
D = bwdist(revim);
D2 = imcomplement(D);
% Suppress shallow minima
D3 = imhmin(D2,1);
L = watershed(D3);
grainsBW(L == 0) = 0;
Note that i just grabbed a quick screen capture o fthe image...hence the rgb2gray call. Also, your threshold might be different from mine.
Alyssa on 11 Mar 2011
Edited: Walter Roberson on 19 Jul 2012
Thank you for the advice on the eccentricity, it seems to be working. Currently we are displaying the calculated eccentricities overlayed on top of the image. Is there a way to limit the range of eccentricities displayed? Below is our code:
clc; clear all;
grains = imread('M65_A.png');
grainsBW = im2bw(grains,105/255);
revim = imcomplement(grainsBW);
D = bwdist(revim);
D2 = imcomplement(D);
% Suppress shallow minima
D3 = imhmin(D2,1);
L = watershed(D3);
grainsBW(L == 0) = 0;
[label num] = bwlabel(grainsBW,8); % label all the connected objects
graindata = regionprops(label,'MinorAxisLength', 'MajorAxisLength', 'Area', 'Perimeter','centroid');
%We can evaluate the "roundness" property using three parameters
%a. Comparing the longest and shortest diameter for the object,
%b. Comparing the Area with the formula pi*r^2.
%c. Comparing the Perimeter with the formula 2*pi*r
%The "scores" are normalized so that the output value represent "best fit"
for cnt = 1:length(graindata)
score1(cnt) = abs(1-(graindata(cnt).MajorAxisLength-graindata(cnt).MinorAxisLength)...
score2(cnt) = abs(1 - abs(pi*((mean([graindata(cnt).MajorAxisLength,...
score3(cnt) = abs(1 - abs(pi*(max([graindata(cnt).MajorAxisLength,...
score = mean([score1;score2;score3]);
%This result can be displayed visually by placing the text over the image
for cnt = 1:length(graindata)
text(graindata(cnt).Centroid(1),graindata(cnt).Centroid(2),num2str(score(cnt)), 'color','red');
xlswrite('eccentricity_A.xls', score', 'Eccentricities', 'A1')
Brett Shoelson
Brett Shoelson on 11 Mar 2011
First, note that REGIONPROPS will return eccentricity directly; your method of calculating it from Major and Minor axes lengths is not exactly correct. (If you want to calculate it directly, do so by computing the ratio of the minor axis length to the major axis length; it will be scaled on (0,1]). Second, for performance efficiency, consider replacing your call to BWLABEL with one to BWCONNCOMP. Finally, to limit the display of eccentricities, just index to find which are the ones you're interested in. For instance:
eccentricityMask = eccentricity > 0.2 & eccentricity < 0.8;
eccentricitiesToShow = eccentricity(eccentricityMask);
Use those same indices to select/filter the corresponding Centroid pairs. Then you can loop over the good eccentricities:
for ii = 1:numel(eccentricitiesToShow)
Make sense?

Sign in to comment.

Walter Roberson
Walter Roberson on 28 Feb 2011
Edited: Walter Roberson on 15 Jul 2012
I had a look at the image, and it is obvious to me that you cannot separate the objects from that single image alone. You have many objects that are conglomerates of smaller objects, such as can happen when independent crystals grow together. Thus, when two potentially different crystals touch, it is necessary to determine whether they are independent or conglomerated. Due to the lack of depth information and due to some details being hidden, it is not possible to determine that in all cases; I find particular examples of ambiguity in several places near the right hand side of the image.

Sean de Wolski
Sean de Wolski on 4 Mar 2011
Voxel surface area/count are both very easy. Marching Cubes is significantly more difficult and, to my knowledge, there's nothing on the file exchange that does it. It's also, according to many journal articles, the best method. You could also look into using an isosurface and calculating the surface area of it; there is a function to do this on the FEX.
It all really depends on what you need though. If you just want "sphericity", I would give the voxel faces method a try. Once a method is developed you can figure out what it returns as a "sphericity" for various objects and correlate that to your data. I did exactly this for my research and found that for an ideal voxelized sphere there was a sphericity asymptote at about 0.66. Knowing that for a true sphere, sphericity is 1, I could base my measurement off of this scale.
Also, of note, since you only have a 2-dimensional image it appears (Or are you trying to expand this to 3d?) The Image Processing Toolbox has a few built in functions that can do this for you (bwarea, bwperim). Marching squares, is also easier to program.
Good Luck!

anto on 15 Jul 2012
Use 'bwlabel' or 'bwlabeln'
  1 Comment
Image Analyst
Image Analyst on 15 Jul 2012
Why do you think that will work? The vast majority of particles are touching and overlapping so there won't be anywhere near the number of labels as there are particles. You can't use bwlabel successfully on that image until all the particles have been separated, which is really what his question was about.

Sign in to comment.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!