Estimation of particle size in the image

Hi,

I am trying to estimate size of the particle in the attached image. As this image has both the particles and other shapes, hence, it is bit difficult to estimate the size of the particle alone. Could someone help with this?

Answer 1

[Image Analyst]

Yes, it's pretty easy. Would be better if you could also post an image with just the background in there. Can we assume that the stream is always in a v ertical line along the same column numbers in the image? Might not need the background image if we can just mask out a thin strip in the middle.

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.

If you still need help, post your original image -- not the one with a huge white surrounding frame on it.

[Turbulence Analysis]

Hi Image Analyst,

Thanks very much for the suggestion. I indeed tried your Image segmentation tutorial, it's very helpful. Infact I had used it in the past for other applications.

In the current application, the issue I am facing is background correction, for e.g. In the image there are some inherent particles present which is actually originaed from the dust in the lens. This shouldn't be accounted while calculating the size.

I am attaching you here the actual and background images. Moreover, the stream is relatively verticaly, but not always in the same line (for e.g. please refer Image 2).

[Turbulence Analysis]

Sorry forgot to mention, another major issue is presence of thin strip in the middle.

Answer 2

[Image Analyst]

Try this:

% Demo by Image Analyst
% Initialization steps:
clc;    % Clear the command window.
close all;  % Close all figures (except those of imtool.)
clear;  % Erase all existing variables. Or clearvars if you want.
workspace;  % Make sure the workspace panel is showing.
format long g;
format compact;
fontSize = 16;

%--------------------------------------------------------------------------------------------------------
%    READ IN BACKGROUND IMAGE
folder = pwd;
baseFileName = "BACKGROUND.bmp";
fullFileName = fullfile(folder, baseFileName);
% Check if file exists.
if ~isfile(fullFileName)
    % The file doesn't exist -- didn't find it there in that folder.
    % Check the entire search path (other folders) for the file by stripping off the folder.
    fullFileNameOnSearchPath = baseFileName; % No path this time.
    if ~exist(fullFileNameOnSearchPath, 'file')
        % Still didn't find it.  Alert user.
        errorMessage = sprintf('Error: %s does not exist in the search path folders.', fullFileName);
        uiwait(warndlg(errorMessage));
        return;
    end
end
% Read in image file.
backgroundImage = imread(fullFileName);
% Get size
[rows, columns, numberOfColorChannels] = size(backgroundImage)
if numberOfColorChannels == 3
    % Convert to gray scale if needed.
    backgroundImage = rgb2gray(backgroundImage);
end
% Display the image.
subplot(2, 3, 1);
imshow(backgroundImage);
axis('on', 'image');
impixelinfo;
title('Background Image', 'FontSize', fontSize, 'Interpreter', 'None');
% Maximize window.
g = gcf;
g.WindowState = 'maximized';
g.Name = 'Demo by Image Analyst';
g.NumberTitle = 'off';
drawnow;

%--------------------------------------------------------------------------------------------------------
%    READ IN TEST IMAGE
folder = pwd;
baseFileName = "image 1.bmp";
fullFileName = fullfile(folder, baseFileName);
% Check if file exists.
if ~isfile(fullFileName)
    % The file doesn't exist -- didn't find it there in that folder.
    % Check the entire search path (other folders) for the file by stripping off the folder.
    fullFileNameOnSearchPath = baseFileName; % No path this time.
    if ~exist(fullFileNameOnSearchPath, 'file')
        % Still didn't find it.  Alert user.
        errorMessage = sprintf('Error: %s does not exist in the search path folders.', fullFileName);
        uiwait(warndlg(errorMessage));
        return;
    end
end
% Read in image file.
grayImage = imread(fullFileName);
% Get size
[rows, columns, numberOfColorChannels] = size(grayImage)
if numberOfColorChannels == 3
    % Convert to gray scale if needed.
    grayImage = rgb2gray(grayImage);
end
% Display the image.
subplot(2, 3, 2);
imshow(grayImage);
axis('on', 'image');
impixelinfo;
title('Original Image', 'FontSize', fontSize, 'Interpreter', 'None');

%--------------------------------------------------------------------------------------------------------
% DIVIDE THE TEST IMAGE BY THE BACKGROUND IMAGE.
backgroundCorrectedImage = double(grayImage) ./ double(backgroundImage);
% Crop out middle part
backgroundCorrectedImage = backgroundCorrectedImage(:, 400:600);
subplot(2, 3, 3);
imshow(backgroundCorrectedImage, []);
axis('on', 'image');
impixelinfo;
title('Background corrected Image', 'FontSize', fontSize, 'Interpreter', 'None');
subplot(2, 3, 4);
% Show histogram.
histogram(backgroundCorrectedImage)
grid on;
title('Histogram of Background Image', 'FontSize', fontSize, 'Interpreter', 'None');

%--------------------------------------------------------------------------------------------------------
% CREATE MASK
% It's where the image is not near 1.
tolerance = 0.4; % Adjust as needed.
mask = abs(backgroundCorrectedImage - 1) > tolerance;
subplot(2, 3, 5);
imshow(mask);
axis('on', 'image');
impixelinfo;
title('Initial Mask Image', 'FontSize', fontSize, 'Interpreter', 'None');

%--------------------------------------------------------------------------------------------------------
% Fill any holes.
mask = imfill(mask, 'holes');
% Delete top stream and any blobs touching the edge of the image (because they are partial droplets not full droplets).
mask = imclearborder(mask);
% Label each blob with 8-connectivity, so we can make measurements of it
[labeledImage, numberOfBlobs] = bwlabel(mask, 8);
% Apply a variety of pseudo-colors to the regions.
coloredLabelsImage = label2rgb (labeledImage, 'hsv', 'k', 'shuffle');
% Display the pseudo-colored mask image.
subplot(2, 3, 6);
imshow(coloredLabelsImage);
axis('on', 'image');
impixelinfo;

%--------------------------------------------------------------------------------------------------------
% GET AREAS AND DIAMETERS
blobMeasurements = regionprops(mask, 'Area', 'EquivDiameter');
numberOfDroplets = size(blobMeasurements, 1)
allAreas = [blobMeasurements.Area]
allDiameters = [blobMeasurements.EquivDiameter]
caption = sprintf('Final Mask Image with %d Droplets', numberOfDroplets);
title(caption, 'FontSize', fontSize, 'Interpreter', 'None');

numberOfDroplets =

5

allAreas =

586 561 447 569 55

allDiameters =

Columns 1 through 4

27.3151674571987 26.7261554398762 23.8566149421207 26.9160417029382

Column 5

8.368283871884

Answers are in pixels. Multiply by microns per pixel to get your answers in microns or square microns.

[Turbulence Analysis]

Thank you very much! Is that cropping in x axis nescessary in the 'Background corrected Image'?

Can we get the position (x, y) and eccentricity of the particles?

[Turbulence Analysis]

Hello Again.

Instead of

blobMeasurements = regionprops(mask, 'Area', 'EquivDiameter');

I tried 'all'

blobMeasurements = regionprops(mask, 'all');

with this function I am getting both 'centroid', 'Eccentricity'.

Thanks!

[Image Analyst]

No you don't need to crop though I just did that to make sure I didn't have any noise blobs included. You could just blacken the mask in those regions, if necessary, and keep the image the same size.

mask(:, 1:400) = false;
mask(:, 600:end) = false;

[Turbulence Analysis]

Thanks!

Just a small quick follow up question,

Is it possible to convert the image in pixels to real coridnates (i.e. mm) before the processing?. For the size I can do this at the end by using pixel scale factor, however for centroids, the conversion to mm is required to represent its position on negative and positive side of the axis.

The below code working to displaying the image in real coridinates.

SizeX = 1024;  
SizeY = 1024;
Ymm = zeros(SizeX,SizeY);
Xmm = zeros(SizeX,SizeY);
MF = (10/178); % Pixel Resolution (scale factor)
Xref = 500; 
Yref = 26; 
count_pixel = 0;
for i = 1:SizeX
    for j = 1:SizeY
        count_pixel = count_pixel + 1;
        Xmm(i,j) = MF*(i-Xref);  
        Ymm(i,j) = MF*(j-Yref); 
        
    end
end

I = imread('Image 1.bmp');
    Y = Ymm (1,:)*(1);
    X = Xmm(:,1);
I=double(I);
imagesc(X,Y,I);

[Turbulence Analysis]

Sorry, I find a another way to sort this out.

Please Never mind!

Thanks again!

[Turbulence Analysis]

Sorry for the inconvienence. Still I have issue in representing the centroid values correctly.

For example, in the below shown plot I have plotted X cordinates of the centroid, here 0 mm lies at 350 pixels. This needs to be presented in mm cordinates. The X cordinates of the centroid is attached here

[Image Analyst]

Your image is 1024 pixels wide, I guess. What is the field of view in mm?

What is the centroid? Each image will have several droplets and thus several centroids. Which one(s) are you plotting?

[Turbulence Analysis]

The scale factor is 17.5 pixels/mm. The above plot is the centroid (only x cordinate) of all the droplets detected in the all the images (in this case 5000 images). The idea is is to illustrate the x location of all the detected droplets

Maybe you can try this in the single image i.e. Image 1 shared yesterday...

[Image Analyst]

So you just do this (for one image):

blobMeasurements = regionprops(mask, 'Area', 'EquivDiameter', 'Centroid', 'Eccentricity');
xy = vertcat(blobMeasurements.Centroid)
xInPixels = xy(:, 1);
xInMm = xInPixels / 17.5;
figure;
subplot(1, 2, 1)
histogram(xInMm);
grid on;
title('Histogram of x values');
subplot(1, 2, 2);
y = 1 : numel(xInMm);
plot(x, y, 'b.', 'MarkerSize', 20)
grid on;
xlabel('x location in mm')
ylabel('Droplet number')

for multiple images, accumulate all the x of each individual image into one vector and show that.

[Turbulence Analysis]

Hi,

I did this already. However, as shown in the image here (2.bmp) particles are sometime in the negative side of x axis.

The histogram obtained with simple scaling shows all the particles in the positive x values as it just divides the pixel coridinates by psitive value of 17.5. Is there any other way?

[Image Analyst]

Attach your m-file. It looks like you maybe used 'XData' and 'YData' options in imshow(). The x values should be the distance from the left side of the image unless you somehow introduced a shift in the axis values.

What do you want? Do you want both negative and positive coordinates, or do you just want positive coordinates? Or do you want displacement from the centerline of the stream?

[Turbulence Analysis]

The below is the code I used to display image in real cordinates in mm.

SizeX = 1024;  
SizeY = 1024;
Ymm = zeros(SizeX,SizeY);
Xmm = zeros(SizeX,SizeY);
MF = (10/178); % Pixel Resolution (scale factor)
Xref = 500; 
Yref = 26; 
count_pixel = 0;
for i = 1:SizeX
    for j = 1:SizeY
        count_pixel = count_pixel + 1;
        Xmm(i,j) = MF*(i-Xref);  
        Ymm(i,j) = MF*(j-Yref); 
        
    end
end
    grayImage = imread('1.bmp');
    backgroundImage=imread('background.bmp');
    backgroundCorrectedImage = double(grayImage) ./ double(backgroundImage);
    backgroundCorrectedImage= imcrop(backgroundCorrectedImage,[150 1 700 1024]);
    Y = Ymm (1,:)*(1);
    X = Xmm(:,1);
    I=double(backgroundCorrectedImage);
II = imagesc(X,Y,I);
axis on
colormap(gray)

Yes, I want both positive and negative cordinates of 'x' values of the centroid function calculated from the below script. As we are implementing the below script on the images in the pixel cordinates, so it is difficult to obtain x value of centroid in negative and positive cordinates.

Finally, I just want to do scatter plot of only x values of the centroid.

blobMeasurements = regionprops(mask, 'Area', 'EquivDiameter', 'Centroid', 'Eccentricity');

[Image Analyst]

You forgot to attach the m-file like I asked. You attached just a snippet from it. Attach the whole thing.

Since different images have different droplets all over the place, where is your global x=0 location? Or does it vary image to image, like each image's origin depends on where the central stream is?

[Turbulence Analysis]

Hi,

Could you please help with this...

[Image Analyst]

Did you try subtracting 350?

xy = vertcat(blobMeasurements.Centroid)
xInPixels = xy(:, 1) - 350;
xInMm = xInPixels / 17.5;

[Turbulence Analysis]

Yes I tried that too, unofrtunatelt it doesn't help.

[Image Analyst]

Did you attach the correct image? I get

Error using imread>get_full_filename (line 579)

Unable to find file "0_40ml_1_C001H001S0001000001.bmp".

Error in imread (line 372)

fullname = get_full_filename(filename);

Error in Particle_size (line 27)

grayImage = imread(fullFileName);

[Turbulence Analysis]

For me it works well. Please try with the attached code and images

[Turbulence Analysis]

Hi,

You previous suggestion of substracting the values by - 350 indeed working. Prevously I did some mistake with the scale factor. Thanks for all the suggestions and sorry for the inconvinience.

[Turbulence Analysis]

Hi, Just a quick followup, as shown in the below attached image, Is it possible to find the both the X and Y values of the dicontinuity point. Right now the code provides only Y cordinate values (i.e. height).

Thanks!!

[Image Analyst]

Assuming you have the bottom y (row), you can get the x value(s) like this

lastRow = mask(row, :);
leftMostX = find(lastRow, 1, 'first')
rightMostX = find(lastRow, 1, 'last')
allX = find(lastRow)
meanX = mean(allX)

[Turbulence Analysis]

Hi,

Sorry, I don't have y (row) in the code.

[Image Analyst]

Then how did you draw the bounding box on the upper part of the stream? If you drew that red box with something like rectangle then you must know the row of the bottom y value.

[Turbulence Analysis]

with command props = regionprops(binaryImage , 'BoundingBox'); I just got the below shown values, in which column 4 representing the y bottom value and the value in the column 1 may be representing the x position, but I looked on the attached image, unfortunately although Y position is correct but not the x value.

[Turbulence Analysis]

Finally, I found it. The X values are in props.Extrema;

The below command would do

props = regionprops(binaryImage , 'all');
B1 = props.Extrema;
X = mean (B1(:,1);