Particle size distribution using image processing in MATLAB

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.

1. Apply adaptive thresholding on the gray level image.
2. Find distance transform on the thresholded image.
3. Invert the distance transform and assign a high value to the background.
4. Extract marker using gray scale reconstruction and assign zero value to marker.
5. Apply watershed transformation on the marker embedded inverted distance image to obtain the segmented image.
6. 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?

[Image Analyst]

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]

the fat cells are the white ones. I'm doing marker based watershed in order to remove the small particles from the background and get the boundaries of fat cells. I need to get the radius of cells in order to plot their size against their number. I'm following the paper in this link which is similar to my work. http://www.sciencedirect.com.argo.library.okstate.edu/science/article/pii/S0040609000013961

[Image Analyst]

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]

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

[amrutha Priya]

I can email you the article, so that you can get a clear idea about this question.

[amrutha Priya]

1) http://i.imgur.com/zMzW09X.jpg?1 2) http://i.imgur.com/owefODD.jpg?1 3) http://i.imgur.com/7VirMDn.jpg?1

these are the outputs that i need.

[Image Analyst]

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]

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]

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);

[amrutha Priya]

which function can i use to get the plot in the form of a curve instead of bars. Like in the third image.

[Image Analyst]

You can use plot() if you want a line plot.

[amrutha Priya]

http://i.imgur.com/VT70iYb.jpg

I got the following plot after all the simulation. Now, I want to fit it into exponential curve. Can you tell me how to do that.

[Image Analyst]

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]

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]

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]

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]

Dear ImageAnalyst,

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

[Image Analyst]

Sorry - the droplets are not individually resolvable.

[Jordan Rayner]

Hi Image Analyst

Can the particle size, of the particles in focus, be determined?

[Image Analyst]

@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]

@Aditya Sai Deepak See my Image Segmentation Tutorial in my File Exchange:

Image Analyst's File Exchange

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]

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]

@Aditya Sai Deepak, see the FAQ for code samples on how to process a sequence of files:

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]

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]

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]

@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.

Answer 3

[Maria Luisa Muñoz Leon]

Hi

Can you give the name of the dataset please ?

[Image Analyst]

What do you mean by dataset? The image that the original poster was using was attached to the original post. Just download it.

[DGM]

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

https://jcp.bmj.com/content/62/6/481