Image Processing for Grain Size Analysis

22 views (last 30 days)
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)?

Answers (2)

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.
  1 Comment
Christopher Atkinson
Christopher Atkinson on 23 Sep 2022
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.

Sign in to comment.


Image Analyst
Image Analyst on 23 Sep 2022
Edited: Image Analyst on 23 Sep 2022
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');
  5 Comments
Christopher Atkinson
Christopher Atkinson on 30 Oct 2022
Edited: Christopher Atkinson on 30 Oct 2022
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));

Sign in to comment.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!