How to calculate Spray Cone Angle

Dear everyone,

I have to calculate the spray cone angle based on the spray image that has been taken from my project experiment. I took reference from the Image Analyst code in separate section and the results is also attached. I changed the binaryImage value and also the thresholdValue to the value that representing the most similar actual spray boundary.

The question is, How can I know if the calculated value of spray angle is measured based on the spray origin? (marked by the red pen in attachment 1), and is there a way to fixed the spray origin into the specific position?

Because when I imported the image without spray, the code still calculated the cone angle is around 19 degree even though there were no spray in the images. (Please kindly refer to the attachment 2)

The image that will be processed is also attached.

Please kindly guide me to understand this matter.

Thank you again for your support.

Answer 1

[Image Analyst]

Try this:

% Demo to compute the cone angle of a liquid spray from a nozzle.
% Initialization steps.
clc;
close all;
workspace; % Make sure the workspace panel with all the variables is showing.
format longg;
format compact;
fontSize = 18;

%---------------------------------------------------------------------------
% Check that user has the Image Processing Toolbox installed.
hasIPT = license('test', 'image_toolbox');
if ~hasIPT
    % User does not have the toolbox installed.
    message = sprintf('Sorry, but you do not seem to have the Image Processing Toolbox.\nDo you want to try to continue anyway?');
    reply = questdlg(message, 'Toolbox missing', 'Yes', 'No', 'Yes');
    if strcmpi(reply, 'No')
        % User said No, so exit.
        return;
    end
end

%---------------------------------------------------------------------------
% Read in a standard MATLAB color demo image.
folder = pwd;
baseFileName = 'DE10-1200-500-30_00071.png';
% Get the full filename, with path prepended.
fullFileName = fullfile(folder, baseFileName);
if ~exist(fullFileName, 'file')
    % Didn't find it there.  Check the search path for it.
    fullFileName = baseFileName; % No path this time.
    if ~exist(fullFileName, 'file')
        % Still didn't find it.  Alert user.
        errorMessage = sprintf('Error: %s does not exist.', fullFileName);
        uiwait(warndlg(errorMessage));
        return;
    end
end
grayImage = imread(fullFileName);
% Get the dimensions of the image.  numberOfColorBands should be = 3.
[rows, columns, numberOfColorBands] = size(grayImage);
if numberOfColorBands > 1
    % If it's really color, then convert to gray scale by taking the green channel.
    grayImage = grayImage(:,:,2);
end

%---------------------------------------------------------------------------
% Display the original image.
hFig = figure;
subplot(2, 3, 1);
imshow(imadjust(grayImage)); % Use imadjust to increase contrast, but only for the displayed image.
impixelinfo;
axis on;
title('Original Image', 'FontSize', fontSize);
% Enlarge figure to full screen.
hFig.WindowState = 'maximized'
hFig.Name = 'Demo by Image Analyst'
hFig.NumberTitle = 'off'

%---------------------------------------------------------------------------
% Erase the circular part surrounding the image.  It's perfectly centered of known, fixed radius.
grayImage = CircleMaskImage(grayImage);

%---------------------------------------------------------------------------
% Let's compute and display the histogram.  Optional -- just for fun/info.
[pixelCount, grayLevels] = imhist(grayImage);
% Suppress bins at 0 and last bin because they're so tall we can't see the rest of it.
pixelCount(1) = 0;
pixelCount(end) = 0;
subplot(2, 3, 2);
bar(grayLevels, pixelCount);
grid on;
title('Histogram of original image', 'FontSize', fontSize);
xlim([0 grayLevels(end)]); % Scale x axis manually.

%---------------------------------------------------------------------------
% Threshold the image
lowThreshold = 0;
highThreshold = 8800;
% Interactively and visually set a threshold on a gray scale image.
% https://www.mathworks.com/matlabcentral/fileexchange/29372-thresholding-an-image?s_tid=srchtitle
% [lowThreshold, highThreshold, lastThresholdedBand] = threshold(lowThreshold, highThreshold, grayImage);
binaryImage = (grayImage >= lowThreshold) & (grayImage <= highThreshold);
% Display the image.
subplot(2, 3, 3);
imshow(binaryImage);
axis on;
impixelinfo;
title('Initial Binary Image', 'FontSize', fontSize);

%---------------------------------------------------------------------------
% Do an closing to smooth out the edges.
se = strel('disk', 3, 0);
binaryImage = imclose(binaryImage, se);

%---------------------------------------------------------------------------
% Fill any holes that might be present
binaryImage = imfill(binaryImage, 'holes');

%---------------------------------------------------------------------------
% Extract the blob with the largest area.
numberToExtract = 1;
binaryImage = bwareafilt(binaryImage, numberToExtract);
% Display the image.
subplot(2, 3, 4);
imshow(binaryImage);
axis on;
title('Biggest Blob, Final Mask', 'FontSize', fontSize);
drawnow;

%---------------------------------------------------------------------------
% Go down the image line by line finding the width
widths = nan(1, rows);
leftEdge = nan(1, rows);
rightEdge = nan(1, rows);
for row = 1 : rows
    thisRow = binaryImage(row, :);
    leftIndex = find(thisRow, 1, 'first');
    if ~isempty(leftIndex)
        % There's something there
        leftEdge(row) = leftIndex;
        rightEdge(row) = find(thisRow, 1, 'last');
        widths(row) = rightEdge(row) - leftEdge(row);
    end
end

%---------------------------------------------------------------------------
% Find the top row.  It's the first non-nan
topRow = find(~isnan(widths), 1, 'first')
% Find the widest row.
[maxWidth, bottomRow] = max(widths)
% Plot the widths
subplot(2, 3, 5);
plot(1:rows, widths, 'b-', 'LineWidth', 2);
grid on;
title('Widths as a function of line number', 'FontSize', fontSize);
axis on;

%---------------------------------------------------------------------------
% Now find min width at top and max width.
% This could be an algorithm all by itself, but since this is just a simple demo,
% I'm going to hard code in values that I can see from the plot of the widths.
% I'm going to fit between rows 90 and 240
x = topRow : bottomRow;
% Show top and bottom as red horizontal lines.
subplot(2, 3, 4);
yline(topRow, 'Color', 'r', 'LineWidth', 2)
yline(bottomRow, 'Color', 'r', 'LineWidth', 2)
subplot(2, 3, 5);
xline(topRow, 'Color', 'r', 'LineWidth', 2)
xline(bottomRow, 'Color', 'r', 'LineWidth', 2)

%---------------------------------------------------------------------------
% Get the left angle.
y = leftEdge(x);
leftCoefficients = polyfit(x, y, 1);
leftAngle = atand(leftCoefficients(1))
subplot(2, 3, 6);
plot(x, y, 'b-', 'LineWidth', 2);
hold on;
y = rightEdge(x);
rightCoefficients = polyfit(x, y, 1);
rightAngle = atand(rightCoefficients(1))
plot(x, y, 'r-', 'LineWidth', 2);
grid on;
coneAngle = abs(leftAngle) + abs(rightAngle);

%---------------------------------------------------------------------------
% Plot the fitted lines
yLeftFit = polyval(leftCoefficients, x);
plot(x, yLeftFit, 'b-', 'LineWidth', 2);
yRightFit = polyval(rightCoefficients, x);
plot(x, yRightFit, 'r-', 'LineWidth', 2);
legend('left', 'right', 'Location', 'east');
caption = sprintf('Cone Angle = %.1f degrees', coneAngle);
title(caption, 'FontSize', fontSize);
message = sprintf('Done with demo.\nThe cone angle is %.1f degrees.', coneAngle);
uiwait(helpdlg(message));

%==========================================================================================================
% Function to erase the image outside a circle.
function grayImage = CircleMaskImage(grayImage)
% Get the size of the image.
[imageSizeY, imageSizeX, numberOfColorChannels] = size(grayImage);
% Create a logical image of a circle with specified
% diameter, center, and image size.
% First create the image.
[columnsInImage, rowsInImage] = meshgrid(1:imageSizeX, 1:imageSizeY);
% Next create the circle in the image.
centerX = imageSizeX / 2;
centerY = imageSizeY / 2;
radius = 152;
% Display circle in the overlay above the current axes.
viscircles([centerX, centerY], radius, 'Color', 'r');
% Here is where we actually create the mask.
circlePixels = (rowsInImage - centerY).^2 ...
    + (columnsInImage - centerX).^2 <= radius.^2;
% circlePixels is a 2D "logical" array.
% Now, display it.
% image(circlePixels) ;
% colormap([0 0 0; 1 1 1]);
% title('Binary image of a circle');
% Set image outside the circle to some very bright value 
% so it won't be selected when we do a threshold for dark stuff.
grayImage(~circlePixels) = 65535;
end

[I Komang Gede Tryas Agameru Putra]

Dear Image Analyst,

Thank you for your answer and it works like charm!

Thank you again for providing this code.

Best Regards.