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
