Particle size distribution using image processing in MATLAB

Question

amrutha Priya 2013-3-3

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/65623-particle-size-distribution-using-image-processing-in-matlab

评论： Image Analyst 2024-5-2

采纳的回答： Image Analyst

temp.jpg

http://jcp.bmj.com/content/62/6/481/F6.large.jpg

I want to plot the number of fat cells in the liver versus their size. For which I wanted to use the following steps.

Apply adaptive thresholding on the gray level image.
Find distance transform on the thresholded image.
Invert the distance transform and assign a high value to the background.
Extract marker using gray scale reconstruction and assign zero value to marker.
Apply watershed transformation on the marker embedded inverted distance image to obtain the segmented image.
Use gray scale reconstruction to find radius of each droplet and subsequently calculate area of each labelled droplet.

In step 3 how do I assign a high value to the background and can you help me with step 4?

8 个评论
显示 6更早的评论隐藏 6更早的评论

amrutha Priya 2013-3-3

编辑：amrutha Priya 2013-3-3

在 MATLAB Online 中打开

if true
  clc
clear all
close all
rgb=imread('C:\Users\Desktop\picture1.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',3);
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)')
D = bwdist(~bw); 
figure, imshow(D,[]), title('Distance transform of ~bw')
background = imopen(D,strel('disk',15));
figure, surf(double(background(1:8:end,1:8:end))),zlim([0 255]);
set(gca,'ydir','reverse');
I2 = D- background;
figure, imshow(I2);
marker = imsubtract(I2,2);
recon = imreconstruct(marker, mask);
end

Image Analyst 2013-3-3

I don't have time to write or debug it for you. But I did download your image and looked at its color channels and noticed that you'll get a lot better contract just using the green channel than using rgb2gray() because the red and blue channels are practically worthless and you don't want them to ruin your gray scale image.

amrutha Priya 2013-3-4

I just want to know how I can get the plot from the binary image after all the thresholding is done.

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Image Analyst 2013-3-4

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/65623-particle-size-distribution-using-image-processing-in-matlab#answer_77235

在 MATLAB Online 中打开

After you process everything and get to a binary image where you have distinct, separate blobs, you call bwlabel, and call regionprops, like this (untested):

labeledImage = bwlabel(binaryImage)
measurements = regionprops(labeledImage, 'EquivDiameter');
allDiameters = [measurements.EquivDiameter];
numberOfBins = 50; % Or whatever you want.
[diamDistribution binDiameters] = hist(allDiameters, numberOfBins);
bar(binDiameters, diamDistribution, 'BarWidth', 1.0);

7 个评论
显示 5更早的评论隐藏 5更早的评论

Valeriy 2020-12-13

在 MATLAB Online 中打开

Dear ImageAnalyst,

Thanks a lot for your, as usual, clear, efficient, and concise code. I try to apply it for particles' size analysis. To verify I have created artificial test image Arr:

xMax = 2592;
yMax = 1944;
Arr = zeros(yMax,xMax,'uint8');     % test Array
R = [1; 3; 7; 10];                  % "radius" of test particles
nR = [200; 100; 50; 10];            % number of test particles
for r = 1:length(R)                 % radius
    xx = randi(xMax,nR(r));         
    xx = xx(1,:);                   
    yy = randi(yMax,nR(r));
    yy = yy(:,1);                   
    for i = 1:nR(r)
%         I have omitted analysis for violation of array borders...
        Arr(yy(i)-R(r):yy(i)+R(r),xx(i)-R(r):xx(i)+R(r))=1;
    end                             % of for i = 1:nR(r)
end                                 % of for r = 1:length(R)

When I applied your above code, I received practically exact number of test particles. But there are some discrepancies between initial and obtained particles' sizes. Due to methods of test of object designing, initial diameters should be something like [3; 7; 15; 21]

For numberOfBins = 128; binDiameters (diamDistribution) = 3.46(198); 3.94(1); 7.91(100); 17(50); 23.6(10)

For numberOfBins = 1024; binDiameters (diamDistribution) = 3.4(198); 3.91(1); 7.9(100); 16.9(50); 23.7(10);

It seems that second obtained value is obtained due to overlapping of some test particles.

For last case numberOfBins = 1024 ratio calculatedBinDiameters / initialDiameters = 1.129 ± 0.0016.

Perhaps such rather stable value is due to difference between square shape of test particles and circular ones, which was supposed for above calculations.

Image Analyst 2020-12-13

在 MATLAB Online 中打开

test5.m

Not quite sure what you're saying but yes, some particle may overlap when you lay them down, and depending on the width of the histogram bins, different ECDs may be in the same bin.

Here is the full code with your code incorporated plus some improvements.

clc;    % Clear the command window.
clear all;
close all;
workspace;  % Make sure the workspace panel is showing.
format short g;
format compact;
fontSize = 15;
fprintf('Beginning to run %s.m ...\n', mfilename);
% Create image of specified spot sizes at random locations.
xMax = 2592;
yMax = 1944;
Arr = zeros(yMax,xMax,'logical');     % test Array
R = [1; 3; 7; 10];                  % "radius" of test particles
nR = [200; 100; 50; 10];            % number of test particles
totalParticleCount = 0;
for r = 1:length(R)                 % radius
	xx = randi(xMax,nR(r));
	xx = xx(1,:);
	yy = randi(yMax,nR(r));
	yy = yy(:,1);
	for i = 1:nR(r)
		% Make sure x and y are within the borders and not 0.
		xx(i) = min([xMax-R(r), xx(i)]);
		xx(i) = max([R(r)+1, xx(i)]);
		yy(i) = min([yMax-R(r), yy(i)]);
		yy(i) = max([R(r)+1, yy(i)]);
		Arr(yy(i)-R(r):yy(i)+R(r),xx(i)-R(r):xx(i)+R(r)) = true;
		% Log that we've added a particle.  Particle may overlap existing particles though.
		totalParticleCount = totalParticleCount + 1;
	end                             % of for i = 1:nR(r)
end                                 % of for r = 1:length(R)
subplot(1, 2, 1);
imshow(Arr, []);
drawnow;
axis('on', 'image');
% Enlarge figure to full screen.
set(gcf, 'Units', 'Normalized','OuterPosition',[0 0 1 1]);
% Particle analysis
% Ask user for one integer number.
defaultValue = 50;
titleBar = 'Enter an integer value';
userPrompt = 'Enter the integer';
dialogBoxWidth = 100;
caUserInput = inputdlg(userPrompt, titleBar, [1, dialogBoxWidth], {num2str(defaultValue)});
if isempty(caUserInput),return,end % Bail out if they clicked Cancel.
% Round to nearest integer in case they entered a floating point number.
numberOfBins = round(str2double(cell2mat(caUserInput)));
% Check for a valid integer.
if isnan(numberOfBins)
	% They didn't enter a number.
	% They clicked Cancel, or entered a character, symbols, or something else not allowed.
	numberOfBins = defaultValue;
	message = sprintf('I said it had to be an integer.\nTry replacing the user.\nI will use %d and continue.', numberOfBins);
	uiwait(warndlg(message));
end
binaryImage = Arr;
labeledImage = bwlabel(binaryImage);
fprintf('Measuring particles.\n');
measurements = regionprops(labeledImage, 'EquivDiameter');
numParticles = length(measurements);
fprintf('Found %d particles from the %d that we laid down.\n', numParticles, totalParticleCount);
if numParticles == totalParticleCount
	caption = sprintf('%d particles from %d laid down randomly (none overlap)', numParticles, totalParticleCount);
else
	caption = sprintf('%d particles from %d laid down randomly (some overlap)', numParticles, totalParticleCount);
end
title(caption, 'FontSize', fontSize);
allDiameters = [measurements.EquivDiameter];
subplot(1, 2, 2);
histObject = histogram(allDiameters, numberOfBins)
diamDistribution = histObject.Values;
binWidth = histObject.BinEdges(2) - histObject.BinEdges(1);
% Get a number that is the center value of the bin instead of the left edge.
binDiameters = histObject.BinEdges(1:end-1) + binWidth/2;
bar(binDiameters, diamDistribution, 'BarWidth', 1.0);
% Put labels atop the bars
for k = 1 : length(diamDistribution)
	x = binDiameters(k);
	y = diamDistribution(k) + 2;
	if diamDistribution(k) > 0
		caption = sprintf('%d particles\nECD = %.2f', y, x);
		text(x, y, caption, 'FontWeight', 'bold', 'FontSize', 15, 'Color', 'r', 'HorizontalAlignment', 'center');
		% Report how many particles are in this bin
		fprintf('The bin at %6.2f has %5d particles in it.\n', x, diamDistribution(k));
	end
end
fprintf('Found %d particles from the %d that we laid down.\n', numParticles, totalParticleCount);
if numParticles ~= totalParticleCount
	fprintf('Meaning that some overlapped when we laid them down.\n');
end
title('Histogram of Equivalent Diameters in Pixels', 'FontSize', fontSize);
xlabel('Equivalent Diameter in Pixels', 'FontSize', fontSize);
ylabel('Count (# of particles)', 'FontSize', fontSize);
grid on;
fprintf('Done running %s.m ...\n', mfilename);

The bin at 3.42 has 197 particles in it.

The bin at 4.77 has 1 particles in it.

The bin at 7.91 has 99 particles in it.

The bin at 16.89 has 49 particles in it.

The bin at 23.63 has 8 particles in it.

The bin at 24.08 has 1 particles in it.

The bin at 25.43 has 1 particles in it.

Found 356 particles from the 360 that we laid down.

Meaning that some overlapped when we laid them down.

Image Analyst 2024-5-2

Yes, I answered you below in

https://www.mathworks.com/matlabcentral/answers/65623-particle-size-distribution-using-image-processing-in-matlab#comment_3149981

Maybe you forgot you posted there. Please start your own discussion thread rather than posting comments in several different places in an 11 year old thread.

请先登录，再进行评论。

Answer 2

khaled soliman 2021-7-10

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/65623-particle-size-distribution-using-image-processing-in-matlab#answer_743978

GrayImage13.jpg

Dear ImageAnalyst,

I want to determine the particular size of spray droplets which is attached

9 个评论
显示 7更早的评论隐藏 7更早的评论

Image Analyst 2024-5-1

在 MATLAB Online 中打开

spatial_calibration_demo.m

Everything is in pixels. You need to know how many nm per pixel and then multiply the diameter in pixels by the number of nm/pixel.

diameterInNm = diameterInPixels * nmPerPixel; % Pixels cancel out and you're left with units of nm.

See attached demo.

Image Analyst 2024-5-2

在 MATLAB Online 中打开

@Aditya Sai Deepak please start your own discussion thread so we don't keep bugging @amrutha Priya with emails of activity on this thread. Attach your code and a couple of images, zipped up.

There I can help you if changing this

baseFileName = TifFiles(k).name;

to this

baseFileName = theFiles(k).name;

does not work.

请先登录，再进行评论。