Hi Nick,
The underlying algorithm implemented in the function "graythresh()" is indeed the Otsu thresholding, where it calculates the optimal threshold value to separate an image into foreground and background by maximizing the inter-class variance ("sigma_b_squared"). You can verify the implementation of the algorithm by opening the function definition of "otsuthresh()" (which is used in "graythresh()") with the following command in the command window:
open otsuthresh
I have made some corrections to your code in the calculation of threshold ("level1") using the inter-class variance as follows:
% Load and display the image
img_1 = imread('cameraman.jpg');
% Convert the image to grayscale if necessary
if size(img_1, 3) > 1
img_1 = rgb2gray(img_1);
end
figure('Name', 'Otsu Thresholding', 'Position', [100, 100, 1200, 400])
% Original Image
subplot(1, 3, 1)
imshow(img_1)
title('Original Image')
% Compute histogram
counts = imhist(img_1, 256);
num_bins = numel(counts);
% Variables names are chosen according to your code snippet
p = counts / sum(counts);
Propability_1 = cumsum(p);
m_k = cumsum(p .* (1:num_bins)');
m_g = m_k(end);
sigma_b_squared = (m_g * Propability_1 - m_k).^2 ./ (Propability_1 .* (1 - Propability_1));
% Find the location of the maximum value of sigma_b_squared.
% The maximum may extend over several bins, so average together the
% locations. If maxval is NaN, meaning that sigma_b_squared is all NaN,
% then return 0.
maxval = max(sigma_b_squared);
isfinite_maxval = isfinite(maxval);
if isfinite_maxval
idx = mean(find(sigma_b_squared == maxval));
% Normalize the threshold to the range [0, 1].
level1 = (idx - 1) / (num_bins - 1);
else
level1 = 0.0;
end
% Plot histogram with threshold indication
subplot(1, 3, 2)
bar(0:255, counts)
hold on
plot(256*[level1 level1], ylim, 'r', 'LineWidth', 2)
hold off
title('Histogram with Threshold')
xlabel('Intensity')
ylabel('Frequency')
axis tight % Adjust the y-axis limits
% Binarize the image using the optimal threshold
BW1 = imbinarize(img_1, level1);
% Display the binary image
subplot(1, 3, 3)
imshow(BW1)
title('Binary Image after Otsu''s Thresholding')
% Adjust the spacing between subplots
sgtitle('Otsu Thresholding')
Feel free to use your own images as inputs to the algorithm.
To know more about the Otsu thresholding implementation in MATLAB, I suggest you refer the MathWorks documentation pages:
I hope this helps!
