Separating Objects in an Image

30 vues (au cours des 30 derniers jours)
Alyssa
Alyssa le 27 Fév 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
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')

Réponses (4)

Brett Shoelson
Brett Shoelson le 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);
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
  6 commentaires
Alyssa
Alyssa le 11 Mar 2011
Modifié(e) : Walter Roberson le 19 Juil 2012
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
Brett Shoelson le 11 Mar 2011
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

Connectez-vous pour commenter.


Walter Roberson
Walter Roberson le 28 Fév 2011
Modifié(e) : Walter Roberson le 15 Juil 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 le 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
anto le 15 Juil 2012
Use 'bwlabel' or 'bwlabeln'
  1 commentaire
Image Analyst
Image Analyst le 15 Juil 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.

Connectez-vous pour commenter.

Community Treasure Hunt

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

Start Hunting!

Translated by