Particle size distribution using image processing in MATLAB
Show older comments
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 Comments
Image Analyst
on 3 Mar 2013
Which color is the fat cells? Why are you doing marker-based watershed? Are you trying to follow Steve's demo: http://blogs.mathworks.com/steve/2006/06/02/cell-segmentation/
amrutha Priya
on 3 Mar 2013
Image Analyst
on 3 Mar 2013
I can't see the article because I don't have an account. Can you show your code so far? Are you able to adapt Steve's algorithm? It shouldn't be too hard.
amrutha Priya
on 3 Mar 2013
Edited: amrutha Priya
on 3 Mar 2013
amrutha Priya
on 3 Mar 2013
amrutha Priya
on 3 Mar 2013
Image Analyst
on 3 Mar 2013
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
on 4 Mar 2013
Accepted Answer
Image Analyst
on 4 Mar 2013
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 Comments
amrutha Priya
on 4 Mar 2013
Image Analyst
on 4 Mar 2013
You can use plot() if you want a line plot.
amrutha Priya
on 4 Mar 2013
Image Analyst
on 4 Mar 2013
If it's a particle size distribution, it's probably actually a log-normal distribution - you just don't have enough bins to see the shape accurately. Anything you measure with particle sizing is almost always log-normal, no matter what the measurement is. Here is the way your determine the parameters for it: http://en.wikipedia.org/wiki/Log-normal_distribution If you want MATLAB code, I'd suggest you start a new discussion and ask how to do a log normal fit, or do a search of Answers for log-normal.
Valeriy
on 13 Dec 2020
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
on 13 Dec 2020
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
on 2 May 2024
Yes, I answered you below in
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.
More Answers (3)
khaled soliman
on 10 Jul 2021
0 votes
Dear ImageAnalyst,
I want to determine the particular size of spray droplets which is attached
9 Comments
Image Analyst
on 11 Jul 2021
Sorry - the droplets are not individually resolvable.
Jordan Rayner
on 11 Nov 2021
Hi Image Analyst
Can the particle size, of the particles in focus, be determined?
Image Analyst
on 11 Nov 2021
@khaled soliman, this is not an answer to @amrutha Priya. Please start your own question. ALso search for "spray" since I've answered spray questions several times.
@Jordan Rayner you should try to get a better picture first.
Image Analyst
on 21 Mar 2024
It's a generic, general purpose demo of how to threshold an image to find blobs, and then measure things about the blobs, and extract certain blobs based on their areas or diameters.
Basically get the diameters or areas and take the histogram
props = regionprops(mask, 'Area', 'EquivDiameter');
histArea = histogram([props.Area])
histDiam = histogram([props.EquivDiameter])
Aditya Sai Deepak
on 30 Apr 2024
how can we find the size of the particles for few different images and make a histogram of all the images into one histogram graph stating like we found ---- this many particles around --- this size ??
Image Analyst
on 30 Apr 2024
@Aditya Sai Deepak, see the FAQ for code samples on how to process a sequence of files:
Inside the loop over files, accumulate your data for areas into one variable then pass it to histogram
% Set these to null before your loop begins though to initialize them!!!
allAreas = [];
allDiams = [];
for k = 1 : numberOfFiles
% Insert extra code here, to read in image, create mask, etc.
% Now measure particles.
props = regionprops(mask, 'Area', 'EquivDiameter');
% Display histograms for this one image.
histArea = histogram([props.Area])
histDiam = histogram([props.EquivDiameter])
% Accumulate
allAreas = [allAreas, [props.Area]];
allDiams = [allDiams, [props.EquivDiameter]];
end
Then after the loop display the histograms of the grand totals:
histogram(allAreas);
Aditya Sai Deepak
on 30 Apr 2024
suppose if i have a 500 nm images , how will I get the size of the particales in nm thorugh histogram data .? I am confused in converting this ,
The bin Values in histogram data are the actual sizes of the particle or should we use some formula to convert them into nm ??
Image Analyst
on 1 May 2024
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
on 2 May 2024
@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.
Maria Luisa Muñoz Leon
on 23 Apr 2023
0 votes
Hi
Can you give the name of the dataset please ?
2 Comments
Image Analyst
on 23 Apr 2023
What do you mean by dataset? The image that the original poster was using was attached to the original post. Just download it.
DGM
on 23 Apr 2023
Or you can look at the link that's next to the attached image in order to find the source.
Histological patterns in drug-induced liver disease
Witek Jachimczyk
37 minutes ago
Edited: Torsten
1 minute ago
This thread has a good accepted answer for plotting a histogram of droplet sizes from a binary image, but a few parts of the original question were not yet addressed. The script below runs on the liver image attached to the question and fills those gaps. It carries out the steps 3 and 4 that were asked about: after inverting the distance transform it sets the background to Inf, then extracts markers with imextendedmin (an h-minima transform built on imreconstruct) and uses them in a marker-controlled watershed. That separates touching fat droplets that bwlabel alone counts as one object (about 490 droplets instead of 403). It also discards droplets cut off by the image edge, since their true size is unknown and they would make the distribution come out too small. It fits the size distribution instead of only plotting it. The earlier suggestion of a log-normal distribution is implemented with fitdist, and because this image clearly contains two droplet populations (small droplets with a median of about 19 px and large ones with a median of about 39 px), the script also fits a two-component log-normal mixture and compares the two models with BIC. Finally, on the ~1.129 ratio @Valeriy reported: EquivDiameter is the diameter of a circle with the same area, so a square of side s gives 2s/sqrt(pi) = 1.128*s. The ratio comes from the square test particles, not from overlap. Set umPerPixel to your calibration to get sizes in micrometres. The threshold and hMin are the two parameters most worth adjusting for other images.
%% Fat droplet size distribution from a liver histology image
% Marker-controlled watershed separates touching droplets, then the
% equivalent diameters are fitted with a log-normal distribution.
%% Parameters
imageFile = 'temp.jpg'; % image attached to the original question
umPerPixel = 1; % spatial calibration (set to 1 to stay in pixels)
threshold = 205; % green-channel level separating droplets from tissue
minArea = 50; % smallest droplet to keep, in pixels
hMin = 2; % marker depth; raise if over-split, lower if under-split
nBins = 30;
%% 1. Read image and take the green channel (best white-vs-pink contrast)
rgb = imread(imageFile);
g = imgaussfilt(rgb(:,:,2), 1.5);
%% 2. Threshold and clean-up
% Droplets are near-saturated white; Otsu (~171) would include pale tissue
bw = g >= threshold;
bw = imopen(bw, strel('disk', 2));
bw = imfill(bw, 'holes');
bw = bwareaopen(bw, minArea);
%% 3. Inverted distance transform, background set to a high value
D = -bwdist(~bw);
D(~bw) = Inf;
%% 4. Markers by grayscale reconstruction (h-minima) and watershed
markers = imextendedmin(D, hMin);
D2 = imimposemin(D, markers);
L = watershed(D2);
L(~bw) = 0;
%% 5. Discard droplets cut off by the image edge (their size is unknown)
edgeLabels = unique([L(1,:), L(end,:), L(:,1)', L(:,end)']);
L(ismember(L, edgeLabels)) = 0;
%% 6. Measure each droplet
stats = regionprops(L, 'Area', 'EquivDiameter');
stats = stats([stats.Area] >= minArea);
diams = [stats.EquivDiameter]' * umPerPixel;
fprintf('Droplets measured: %d\n', numel(diams));
%% 7. Fit a log-normal distribution
pd = fitdist(diams, 'Lognormal');
fprintf('Log-normal fit: mu = %.3f, sigma = %.3f (median diameter = %.2f)\n', ...
pd.mu, pd.sigma, median(pd));
%% 8. Two-population fit (small + large droplets): log-normal mixture
gm = fitgmdist(log(diams), 2, 'Replicates', 10, 'Options', statset('MaxIter', 1000));
[~, order] = sort(gm.mu);
fprintf('Mixture: median %.2f, sigma %.3f, fraction %.2f\n', [exp(gm.mu(order))'; ...
sqrt(squeeze(gm.Sigma(order)))'; gm.ComponentProportion(order)]);
fprintf('BIC: single log-normal %.1f, mixture %.1f (lower is better)\n', ...
fitgmdist(log(diams), 1).BIC, gm.BIC);
%% 9. Display
if umPerPixel == 1, unitLabel = 'pixels'; else, unitLabel = '\mum'; end
figure('Name', 'Fat droplet sizes', 'Position', [100 100 1100 750]);
subplot(2,2,1); imshow(rgb); title('Image');
subplot(2,2,2); imshow(bw); title('Binary');
subplot(2,2,3); imshow(imoverlay(rgb, boundarymask(L), 'k')); title('Separated droplets');
subplot(2,2,4);
histogram(diams, nBins, 'Normalization', 'pdf'); hold on
x = linspace(0, max(diams), 300);
plot(x, pdf(pd, x), 'r-', 'LineWidth', 2);
x = x(2:end); % mixture pdf in diameter space: f(log d)/d
plot(x, pdf(gm, log(x'))' ./ x, 'k--', 'LineWidth', 2); hold off
xlabel(['Equivalent diameter (' unitLabel ')']); ylabel('Probability density');
legend('Measured', 'Log-normal fit', 'Two-population fit'); title('Size distribution');
Categories
Find more on Data Preprocessing in Help Center and File Exchange
Products
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
