Image Processing for Grain Size Analysis

I'm interesting in calculating wildfire ashes like this .jpeg example below,

The idea is the use a coin as a measurement reference in order to determine the size of the ash. I'm very new to image processing in general as I mainly just know how to read in the image. Is there a way for me to calculate the size of the ash or even code it to where a user can point and click to measure the intermidiate axix of each ash (or something to measure the ash in general)?

Answer 1

[Constantino Carlos Reyes-Aldasoro]

You need to define your task first. What is it that you want to "calculate"? I am going to make an educated guess that it is the size of the ashes, and another is that the ashes are the bright regions. If these two assumptions are correct, then you would start by segmenting these, probably a simple thresholding algorithm based on intensity would be enough to start. Then you would segment the coin (this time based on the known colour) and measure the size to callibrate the size of the ashes. None of these steps are too difficult to implement in Matlab, but what you need to do is grab a book so that you are able to understand the image processing steps as implemented in Matlab. There are lots of books that will help, one of the classic ones is by Gonzalez and Woods and one of the editions has Matlab, and there are others.

[Christopher Atkinson]

Thank you for the response and suggestions! The size of the ashes is definitely one. Finding the bright regions (as ashes) could be useful too since I'll have different pictures with different backgrounds. I need to look into more image processing methods via Matlab since my coding skills are definitely beginner level. Then again, the closest example code I've seen is grain size analysis (although that is slightly different). Just need to find the right book/information that talks about its implementation in Matlab.

Answer 2

[Image Analyst]

See my Image Processing Tutorial in my File Exchange:

Image Analyst's File Exchange

Is the ash black or white or both? Can you put the ash on a different background, one smoother and more uniform and a different color?

You can find the cent coin with the Color Thresholder on the Apps tab of the tool ribbon. Then use the created mask and regionprops to find the diameter. Then compute your spatial calibration factor.

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

[Christopher Atkinson]

Ash is mostly black but could be both. I am unable to change the background. That is one problem I may face as all pieces won't be smooth and uniform.

I will try to use the Color Thresholder this weekend.

[Image Analyst]

If ash is black on a black background it will be very difficult to find by thresholding. In that case you might have to do some spatial processing like edge detection to find the black particles, or stdfilt to find highly textured or changing particles. Or else use deep learning to try to find them.

Try this:

clc;    % Clear the command window.
close all;  % Close all figures (except those of imtool.)
workspace;  % Make sure the workspace panel is showing.
format long g;
format compact;
fontSize = 16;
%===============================================================================
% Get the name of the image the user wants to use.
baseFileName = 'ash.jpeg';
folder = pwd;
fullFileName = fullfile(folder, baseFileName);
% Check if file exists.
if ~exist(fullFileName, 'file')
    % 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 demo image.
rgbImage = imread(fullFileName);
% Get the dimensions of the image.
[rows, columns, numberOfColorChannels] = size(rgbImage)
% Crop off screenshot stuff
% rgbImage = rgbImage(132:938, 352:1566, :);
% Display image.
subplot(2, 1, 1);
imshow(rgbImage, []);
impixelinfo;
axis on;
caption = sprintf('Original Color Image\n%s', baseFileName);
title(caption, 'FontSize', fontSize, 'Interpreter', 'None');
hp = impixelinfo(); % Set up status line to see values when you mouse over the image.
% Set up figure properties:
% Enlarge figure to full screen.
set(gcf, 'Units', 'Normalized', 'OuterPosition', [0 0.05 1 0.95]);
% Get rid of tool bar and pulldown menus that are along top of figure.
% set(gcf, 'Toolbar', 'none', 'Menu', 'none');
% Give a name to the title bar.
set(gcf, 'Name', 'Demo by ImageAnalyst', 'NumberTitle', 'Off')
drawnow;
%=======================================================================================
% Threshold the image to get the disk.
[mask,maskedRGBImage] = createMask(rgbImage);
% Take only the largest blob.
mask = bwareafilt(mask, 1);
% Fill holes
mask = imfill(mask, 'holes');
% Find out the size of all the blobs so we know what size to use when filtering with bwareaopen.
props = regionprops(mask, 'Area', 'EquivDiameter', 'Centroid');
allAreas = sort([props.Area]);
% Display the mask image.
subplot(2, 1, 2);
imshow(mask, []);
impixelinfo;
axis('on', 'image');
drawnow;
pennyRadiusInPixels = props.EquivDiameter / 2;
pennyRadiusInMm = 19.05 / 2;
% Overlay the circle over the image.
viscircles([props.Centroid(1), props.Centroid(2)], pennyRadiusInPixels, 'Color', 'b')
% Compute the spatial calibration factor.
mmPerPixel = pennyRadiusInMm / pennyRadiusInPixels
subplot(2, 1, 2);
caption = sprintf('The spatial calibration is %f mm per pixel', mmPerPixel)
title(caption, 'FontSize', fontSize, 'Interpreter', 'None');
% Plot the borders of all the blobs in the overlay above the original grayscale image 
% using the coordinates returned by bwboundaries().
% bwboundaries() returns a cell array, where each cell contains the row/column coordinates for an object in the image.
% Here is where we actually get the boundaries for each blob.
boundaries = bwboundaries(mask);
% boundaries is a cell array - one cell for each blob.
% In each cell is an N-by-2 list of coordinates in a (row, column) format.  Note: NOT (x,y).
% Column 1 is rows, or y.    Column 2 is columns, or x.
numberOfBoundaries = size(boundaries, 1); % Count the boundaries so we can use it in our for loop
% Here is where we actually plot the boundaries of each blob in the overlay.
subplot(2, 1, 1);
hold on; % Don't let boundaries blow away the displayed image.
for k = 1 : numberOfBoundaries
    thisBoundary = boundaries{k}; % Get boundary for this specific blob.
    x = thisBoundary(:,2); % Column 2 is the columns, which is x.
    y = thisBoundary(:,1); % Column 1 is the rows, which is y.
    plot(x, y, 'r-', 'LineWidth', 2); % Plot boundary in red.
end
hold off;
caption = sprintf('Original image with %d Outlines, from bwboundaries()', numberOfBoundaries); 
title(caption, 'FontSize', fontSize); 
axis('on', 'image'); % Make sure image is not artificially stretched because of screen's aspect ratio.
message = sprintf('Done!\nThe spatial calibration is %f mm per pixel', mmPerPixel)
uiwait(helpdlg(message));
%=========================================================================================================
function [BW,maskedRGBImage] = createMask(RGB)
%createMask  Threshold RGB image using auto-generated code from colorThresholder app.
%  [BW,MASKEDRGBIMAGE] = createMask(RGB) thresholds image RGB using
%  auto-generated code from the colorThresholder app. The colorspace and
%  range for each channel of the colorspace were set within the app. The
%  segmentation mask is returned in BW, and a composite of the mask and
%  original RGB images is returned in maskedRGBImage.
% Auto-generated by colorThresholder app on 23-Sep-2022
%------------------------------------------------------
% Convert RGB image to chosen color space
I = rgb2hsv(RGB);
% Define thresholds for channel 1 based on histogram settings
channel1Min = 0.000;
channel1Max = 1.000;
% Define thresholds for channel 2 based on histogram settings
channel2Min = 0.000;
channel2Max = 1.000;
% Define thresholds for channel 3 based on histogram settings
channel3Min = 0.790;
channel3Max = 1.000;
% Create mask based on chosen histogram thresholds
sliderBW = (I(:,:,1) >= channel1Min ) & (I(:,:,1) <= channel1Max) & ...
    (I(:,:,2) >= channel2Min ) & (I(:,:,2) <= channel2Max) & ...
    (I(:,:,3) >= channel3Min ) & (I(:,:,3) <= channel3Max);
BW = sliderBW;
% Initialize output masked image based on input image.
maskedRGBImage = RGB;
% Set background pixels where BW is false to zero.
maskedRGBImage(repmat(~BW,[1 1 3])) = 0;
end

[Christopher Atkinson]

Hi! Sorry for the delayed response and thank you for the test code. I was trying a few things and had no luck. This works pretty well for the coin. So to detect the ash in this image, I think you suggested edge detection or stdfilt since I can't change the background or consider the ash to be smooth etc. right?

[Image Analyst]

Right. You can try adapthisteq or stdfilt and then do a global threshold with imbinarize

[Christopher Atkinson]

I think I got it. But I ran into another issue. So the size of the coin isn't always the smallest or largest part of the image and I'm having trouble specifying this in my code. I know there is a way where I can use the function "ginput" so I can have the user iput/choose the reference point (the coin). Assuming this has to do with my props function (after line 85). This is important to me cause I want to also keep the measurement value (as you'll see in lines 102 and 103) so I can use that to measure other pieces of the image.

Will post the code below and appreciate your help as always. Also the image I'm using is below since I felt this was an easier one to use.

close all;
clear all;

%% Read in Image
RGB = imread('20210717_Tamarack_1524.jpg');

Error using imread>get_full_filename
File "20210717_Tamarack_1524.jpg" does not exist.

Error in imread (line 372)
fullname = get_full_filename(filename);

% Display Image
figure(1);
clf
imshow(RGB)

%% Convert Truecolor (RGB) Image into Grayscale Image
I = im2gray(RGB);
% Use Median Filter
I = medfilt2(I,[10 10]);
figure(2);
clf
imshow(I)

%% Calculate the Gradient Image and Apply a Threshold
% Use edge and Sobel operator to calculate the threshold value. Then  Tune the threshold value and use edge again to obtain a binary mask that contains the segmented cell
[~,threshold] = edge(I,'sobel');
fudgeFactor = 1;
BWs = edge(I,'sobel',threshold * fudgeFactor);
% Display resulting binary gradient mask
figure(3);
clf
imshow(BWs)

%% Dilate the Image.
% Create two perpendicular linear structuring elements by using strel function.
se90 = strel('line',3,90);
se0 = strel('line',3,0);
% Dilate the binary gradient mask using the vertical structuring element followed by the horizontal structuring element. The imdilate function dilates the image.
BWsdil = imdilate(BWs,ones(15,15));
imshow(BWsdil)
% Fill Interior Gaps
% Fill remaining holes using the imfill function
BWdfill = imfill(BWsdil,'holes');
figure(4);
clf
imshow(BWdfill)

%% Remove Connected Objects
% Remove any objects that are connected to the border of the image using
% imclearborder function. We use 4 as connectivity to remove 2D diagonal
% connections
BWnobord = imclearborder(BWdfill,4);
figure(5);
clf
imshow(BWnobord)

%% Smooth the Object
% We create the diamond structuring element using the strel function in
% order to make the object look natual/smooth
seD = strel('diamond',1);
BWfinal = imerode(BWnobord,seD);
BWfinal = imerode(BWfinal,seD);
figure(6);
clf
imshow(BWfinal)

%% Visualize the Segmentation
% Labelloverlay function allows us to display the mask over the original
% image
figure(7);
clf
imshow(labeloverlay(I,BWfinal))
%% bwareafilt function allows us to filter image, retaining (10) or n objects with
% the largest area
bwsizefilt = bwareafilt(BWfinal, 5);
% bwlabel returns the label matrix L that contains labels for the
% objects found in bwsizefilt
L = bwlabel(bwsizefilt);
figure(8);
clf
% Retains plot
hold on
% Creates pseudocolor plot using the values in matrix L
pcolor(L);
set(gca,'ydir','reverse')
% Sets color shading properties 
shading flat
% Measure properties of image regions using regionprops
% Calculate centroids and area for connected components in the image.
props = regionprops(bwsizefilt, 'Centroid', 'Area', 'EquivDiameter');
for i = 1:length(props)
    plot(props(i).Centroid(1),props(i).Centroid(2), '*k')
end

%% Find out the size of all the blobs so we know what size to use when filtering with bwareaopen.
% props = regionprops(mask, 'Area', 'EquivDiameter', 'Centroid');
allAreas = sort([props.Area]);
% Display the mask image.
% subplot(2, 1, 2);

% Coin radius, can be changed based on reference

% uiwait(msgbox('Choose Reference Point'));

% A = ginput(1);

pennyRadiusInPixels = props(3).EquivDiameter / 2;
pennyRadiusInMm = 19.05 / 2;

% Identify Points
%[x,y] = ginput(1)
%centers = [x,y];
%viscircles(centers,pennyRadiusInMm, 'Color', 'b')

% Overlay the circle over the image.

viscircles([props(3).Centroid(1), props(3).Centroid(2)], pennyRadiusInPixels, 'Color', 'b')
return
% Compute the spatial calibration factor.
mmPerPixel = pennyRadiusInMm / pennyRadiusInPixels
subplot(2, 1, 2);
caption = sprintf('The spatial calibration is %f mm per pixel', mmPerPixel)
title(caption, 'FontSize', fontSize, 'Interpreter', 'None');
% Plot the borders of all the blobs in the overlay above the original grayscale image 
% using the coordinates returned by bwboundaries().
% bwboundaries() returns a cell array, where each cell contains the row/column coordinates for an object in the image.
% Here is where we actually get the boundaries for each blob.
boundaries = bwboundaries(mask);
% boundaries is a cell array - one cell for each blob.
% In each cell is an N-by-2 list of coordinates in a (row, column) format.  Note: NOT (x,y).
% Column 1 is rows, or y.    Column 2 is columns, or x.
numberOfBoundaries = size(boundaries, 1); % Count the boundaries so we can use it in our for loop
% Here is where we actually plot the boundaries of each blob in the overlay.
subplot(2, 1, 1);
hold on; % Don't let boundaries blow away the displayed image.
for k = 1 : numberOfBoundaries
    thisBoundary = boundaries{k}; % Get boundary for this specific blob.
    x = thisBoundary(:,2); % Column 2 is the columns, which is x.
    y = thisBoundary(:,1); % Column 1 is the rows, which is y.
    plot(x, y, 'r-', 'LineWidth', 2); % Plot boundary in red.
end
hold off;
caption = sprintf('Original image with %d Outlines, from bwboundaries()', numberOfBoundaries); 
title(caption, 'FontSize', fontSize); 
axis('on', 'image'); % Make sure image is not artificially stretched because of screen's aspect ratio.
message = sprintf('Done!\nThe spatial calibration is %f mm per pixel', mmPerPixel)
uiwait(helpdlg(message));