Particle size calculation from given image

I want to calculate sizes of each particle from the input image using image processing technique.

Answer 1

[Thiago Henrique Gomes Lobato]

The size for each particle is a task which is extremely hard to automatize, you will always get some error due to different light conditions, particle distribution, overlapping particles etc. A problem that is way more tractable and realistic is to calculate the rock size distribution, which will give you how the sizes are statistically distributed in the image. I took some inspirations from Image Analyst tutorial and made a code that could be used as a start for you to fine tune it. In this way you could at least get a estimation of particle size distribution:

I = imread('image.jpeg');
Ibinary = rgb2gray(I) > 150; % Try different values for your problem
Ibinary = imfill(Ibinary, 'holes');


C = 0.008; % Conversion factor m/pixel. This was a guess, you should be able to find the real one

[labeledImage, numberOfBlobs] = bwlabel(Ibinary);
blobMeasurements = regionprops(labeledImage, 'Centroid','EquivDiameter'); % I believe equivalent diameter is a good way to measure size,
                                                                          % you can also check different metrics
% Get equivalent diameter
EquivDiameter = [blobMeasurements.EquivDiameter]; 
ValidDia = find(EquivDiameter>10); % Some empiric threshold
blobMeasurements = blobMeasurements(ValidDia);
EquivDiameter         = EquivDiameter(ValidDia);
% Probability estimation of size
figure,histogram(EquivDiameter*C,'Normalization','probability'),xlabel('Size[m - uncalibrated]'),ylabel('Probability')

% Display some particles
figure,imshow(Ibinary); hold on
% Display areas on image
for idx = 1:2: length(ValidDia)           % Loop through all blobs.
  Centroid = [blobMeasurements(idx).Centroid(1), blobMeasurements(idx).Centroid(2)];
  DiamSize = ['Diam. = ', num2str((EquivDiameter(idx)*C))];
  text(Centroid(1), Centroid(2), DiamSize, 'Color', 'c');
  viscircles(Centroid,EquivDiameter(idx)/2);
end

[Vikrant Pratap]

Thank you for your answer Thiago.

Particle size distribution will also serve my purpose.

I have used the following code. The blobs generated are incorrect and does not cover all the particles.

Please suggest the changes in the code.

Also can u suggest a way to measure the longest side of particle instead of diameter.

clear all
close all
rgb=imread('image.jpg');
imshow(rgb)
I = rgb2gray(rgb);
imshow(I)
I=adapthisteq(I);
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)')
se = strel('disk',2);
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)')
bw = im2bw(Iobrcbr, graythresh(Iobrcbr));
figure, imshow(bw), title('Thresholded opening-closing by reconstruction (bw)')

[labeledImage, numberOfBlobs] = bwlabel(bw);
blobMeasurements = regionprops(labeledImage, 'Centroid','EquivDiameter');

EquivDiameter = [blobMeasurements.EquivDiameter]; 
ValidDia = find(EquivDiameter>10); % Some empiric threshold
blobMeasurements = blobMeasurements(ValidDia);
EquivDiameter         = EquivDiameter(ValidDia);
% Probability estimation of size
figure,histogram(EquivDiameter,'Normalization','probability'),xlabel('Size[m - uncalibrated]'),ylabel('Probability')
% Display some particles
figure,imshow(bw); hold on
% Display areas on image
for idx = 1:2: length(ValidDia)           % Loop through all blobs.
  Centroid = [blobMeasurements(idx).Centroid(1), blobMeasurements(idx).Centroid(2)];
  DiamSize = ['Diam. = ', num2str((EquivDiameter(idx)))];
  text(Centroid(1), Centroid(2), DiamSize, 'Color', 'c');
  viscircles(Centroid,EquivDiameter(idx)/2);
end

[Thiago Henrique Gomes Lobato]

The longest size of the particle can be achieve with the MajorAxisLength property

regionprops(labeledImage, 'Centroid','EquivDiameter','MajorAxisLength');

For getting the distribution one nice thing you can consider is that you don't need to get all particles right, but rather only a subset of them. As I said in the initial answer it will be extremely difficult to find all and I don't believe this shouldn't be your goal.

In order to improve the results you can consider that the main problem of the rocks is that they overlap too much. So threshold methods like Otsu (graythresh(Iobrcbr)) will mostly likely not work, neither trying to close too much of the holes. I believe the best way to get this, in a binary segmentation paradigm, is to use a very high binary treshold or high amount of erosion in order to separate the intersection as much as possible and them rely in the subset assumptions, so even though you may remove some rocks the result will still remain reasonable. Also note the idx = 1:2: length(ValidDia). To plot all blobs remove this 2, I initially let it there so the image wouldn't be too much filled.

[Vikrant Pratap]

Please help with plotting MajorAxisLength on the binary image.

Answer 2

[Image Analyst]

Since your particles can't be segmented by itnensity, I think you need to look at granulometry. Wikipedia Granulometry