Separating Objects in an Image

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?

http://www.flickr.com/photos/60125138@N08/5483344701/

clear all;   
%Step 1 Image
I=imread('748C_A.jpg');
figure,imshow(I)

%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')

Answer 1

[Brett Shoelson]

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);
imshow(revim,[])
D = bwdist(revim);
D2 = imcomplement(D);
imshow(D2,[])
% Suppress shallow minima
D3 = imhmin(D2,1);
L = watershed(D3);
grainsBW(L == 0) = 0;
imshow(grainsBW)

Note that i just grabbed a quick screen capture o fthe image...hence the rgb2gray call. Also, your threshold might be different from mine.

Cheers,

Brett

[Alyssa]

Thank you, the code did work for separating the particles. I was wondering if there is a way to measure sphericity in Matlab with the SEM images that I have? I know there is other software out there with these capabilities but I do not have access to it.

[Sean de Wolski]

Alyssa,

You'll have to determine your measurement method for surface area and volume. For surface area do you want to use a marching squares/cubes type algorithm or just a pixel/voxel face count?

For the volume do you want just the count of pixel/voxels or do you want something that weights arrangements differently?

[Alyssa]

Which way of calculating surface area will be the most accurate but easier to code? I am not very proficient in Matlab.

[Brett Shoelson]

I bow to Sean de's obviously superior knowledge of sphericity calculations. But if you can assume that 2-dimensional circles are actually planar representations of spheres (i.e, that the objects are symmetric), then perhaps you can just calculate sphericity from eccentricity? The latter is easy to measure in 2-dimensions. (Try REGIONPROPS).

Cheers,

Brett

[Alyssa]

Brett,

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);
imshow(revim,[])
D = bwdist(revim);
D2 = imcomplement(D);
imshow(D2,[])
% Suppress shallow minima
D3 = imhmin(D2,1);
L = watershed(D3);
grainsBW(L == 0) = 0;
imshow(grainsBW)

   [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)...
/max([graindata(cnt).MajorAxisLength,graindata(cnt).MinorAxisLength]));
score2(cnt) = abs(1 - abs(pi*((mean([graindata(cnt).MajorAxisLength,...
graindata(cnt).MinorAxisLength]))/2)^2-graindata(cnt).Area)/graindata(cnt).Area);
score3(cnt) = abs(1 - abs(pi*(max([graindata(cnt).MajorAxisLength,...
graindata(cnt).MinorAxisLength]))-graindata(cnt).Perimeter)/graindata(cnt).Perimeter);
end
score = mean([score1;score2;score3]);

%This result can be displayed visually by placing the text over the image
imshow(grainsBW)
for cnt = 1:length(graindata)
    text(graindata(cnt).Centroid(1),graindata(cnt).Centroid(2),num2str(score(cnt)), 'color','red');
end

xlswrite('eccentricity_A.xls', score', 'Eccentricities', 'A1')

[Brett Shoelson]

Alyssa,

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)

text(pared_graindata(ii).Centroid(1),pared_graindata(ii).Centroid(2),...)

end

Make sense?

Brett

Answer 2

[Walter Roberson]

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.

Answer 3

[Sean de Wolski]

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!

Answer 4

[anto]

Use 'bwlabel' or 'bwlabeln'

[Image Analyst]

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.