Try this:
% Demo by Image Analyst
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 = 22;
markerSize = 40;
%--------------------------------------------------------------------------------------------------------
% READ IN IMAGE
folder = pwd;
baseFileName = 'Probe_1030_80ms_02.bmp';
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
rgbImage = imread(fullFileName);
% Get the dimensions of the image.
% numberOfColorChannels should be = 1 for a gray scale image, and 3 for an RGB color image.
[rows, columns, numberOfColorChannels] = size(rgbImage)
if numberOfColorChannels > 1
% It's not really gray scale like we expected - it's color.
% Extract the blue channel.
grayImage = rgbImage(:, :, 3);
else
grayImage = rgbImage;
end
%--------------------------------------------------------------------------------------------------------
% Display the image.
subplot(2, 3, 1);
imshow(grayImage);
impixelinfo;
axis('on', 'image');
title('Original Gray Scale Image', 'FontSize', fontSize, 'Interpreter', 'None');
hold on
% Maximize window.
g = gcf;
g.WindowState = 'maximized'
drawnow;
%--------------------------------------------------------------------------------------------------------
% Get a histogram
subplot(2, 3, 2);
imhist(grayImage);
grid on;
title('Histogram of Original Gray Scale Image', 'FontSize', fontSize, 'Interpreter', 'None');
%--------------------------------------------------------------------------------------------------------
% Get the vertical profile
% verticalProfile = mean(grayImage, 2);
% subplot(2, 3, 3);
% plot(1:rows, verticalProfile, 'b-');
% grid on;
% xlabel('Row', 'FontSize', fontSize)
% Ylabel('Intensity', 'FontSize', fontSize)
% Assume camera and substrate are fixed and everything below line 426 we can just throw away.
bottomRow = 426;
grayImage(bottomRow:end, :) = 255;
% Take the histogram
counts = histcounts(grayImage(1:bottomRow, :));
%--------------------------------------------------------------------------------------------------------
% Set thresholds for each level.
lowThreshold = 0;
% Get the triangle threshold for the upper threshold
highThreshold = triangle_threshold(counts, 'L', false)
% 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] = threshold(lowThreshold, highThreshold, grayImage)
% Draw threshold line over histogram.
xline(highThreshold, 'Color', 'r', 'LineWidth', 2);
%--------------------------------------------------------------------------------------------------------
% Binarize the image to get a mask.
mask = grayImage >= lowThreshold & grayImage <= highThreshold;
mask = imfill(mask, 'holes');
% Take largest blobs.
mask = bwareafilt(mask, 2);
% Display mask image.
subplot(2, 3, 3);
imshow(mask);
impixelinfo;
axis('on', 'image');
drawnow;
caption = sprintf('Mask with Threshold = %d gray levels', highThreshold);
title(caption, 'FontSize', fontSize, 'Interpreter', 'None');
% Take the blob closest to the center of the image.
props = regionprops(mask, 'Centroid', 'BoundingBox');
indexOfMin = 1; % Initialize.
if length(props) >= 2
% Extract centroids from structure to arrays.
xy = vertcat(props.Centroid);
x = xy(:, 1);
y = xy(:, 2);
% Find distances from each blob's centroid to the center of the image.
distances = sqrt((x - columns/2) .^ 2 + (y - rows/2) .^ 2);
[minDistance, indexOfMin] = min(distances);
end
% Get the bounding box of the droplet nearest the center.
labeledImage = bwlabel(mask);
mask = ismember(labeledImage, indexOfMin);
bb = props(indexOfMin).BoundingBox;
hold on;
rectangle('Position', bb, 'EdgeColor', 'r');
% Crop the image and display it.
croppedImage = imcrop(grayImage, bb);
mask = imcrop(mask, bb);
[rows, columns] = size(mask);
% Display image.
subplot(2, 3, 4:6);
imshow(croppedImage, []);
impixelinfo;
axis('on', 'image');
drawnow;
caption = sprintf('Mask with Threshold = %d gray levels', highThreshold);
title(caption, 'FontSize', fontSize, 'Interpreter', 'None');
%--------------------------------------------------------------------------------------------------------
% Scan to get lowest boundary
topRows = nan(1, columns);
x = nan(1, columns);
for col = 1 : columns
t = find(mask(:, col), 1, 'first');
if ~isempty(t)
topRows(col) = t;
x(col) = col;
end
end
% Get rid of nans, which were where there was no white pixels in a column.
x = x(~isnan(x));
topRows = topRows(~isnan(topRows));
% We need to fit the top boundary to a polynomial.
coefficients = polyfit(x, topRows, 2)
yFit = polyval(coefficients, x);
hold on;
plot(x, topRows, 'g-', 'LineWidth', 2);
plot(x, yFit, 'r-', 'LineWidth', 2);
% yFit = coefficients(1) * x ^ 2 + coefficients(2) * x + coefficients(3)
% The derivative is
% Slope = 2 * coefficients(1) * x + coefficients(2);
% Get the slope on the left end where x = 1.
leftSlope = 2 * coefficients(1) * x(1) + coefficients(2)
% Get the slope on the right end where x = the last value.
rightSlope = 2 * coefficients(1) * x(end) + coefficients(2)
% Get the average of the slopes
meanSlope = mean([abs(leftSlope), abs(rightSlope)])
% Get the average angle.
meanAngle = atand(meanSlope)
caption = sprintf('The mean slope is %f and the mean angle is %f degrees', meanSlope, meanAngle);
title(caption, 'FontSize', fontSize, 'Interpreter', 'None');
% Tell the user
message = sprintf('Done!\nThe mean slope is %f and the mean angle is %f degrees', meanSlope, meanAngle)
uiwait(msgbox(message))
