Plot the s_b_squared with the imahe histogram

Hello!

So I have an image,and I would like to implement image segmentation on it. To be more specific, I'm performing Otsu's thresholding in order to create a binary image. In the process, I figured a way to show where the Otsu's threshold is present in the histogram of the image and have done it with success. Yet when I wanted to mathematically prove,but always based on MATLAB implementation, that the point that I've found is the right one using the formula of s_b_squared, I got a very different implementation of the threshold point.

Below I'll leave to you, both the code that I've written so far, the image that I want to process and the produced outcomes:

clear 
clc

%Ερώτημα Α)

%Αποτύπωση των εικόνων και των παραγόμενων, από την κατωφλιοποίηση εικόνων

figure(1) 

img_1=imread('BP1.jpg');
img_2=imread('BP2.jpg');

level1=graythresh(img_1);
level2=graythresh(img_2);

BW1=imbinarize(img_1,level1);
BW2=imbinarize(img_2,level2);

subplot(1,2,1)
imshowpair(img_1,BW1,'montage')
title("Αποτέλεσμα της κατωφλιοποίησης για την εικόνα 1")
subplot(1,2,2)
imshowpair(img_2,BW2,'montage');
title("Αποτέλεσμα της κατωφλιοποίησης για την εικόνα 2")

figure(2) %Ιστόγραμμα με υπόδειξη του σημείου κατωφλίου

subplot(1,2,1)


imhist(img_1)

hold on 
plot(256*[level1 level1],ylim,'r','LineWidth',2);
legend({'Ιστόγραμμα','Κατώγλι κατά Otsu'},'Location','northwest')
title("Εμφάνιση του σημείου κατωφλιοποίησης,κατά Otsu,στην εικόνα 1")
hold off

subplot(1,2,2)


imhist(img_2)

hold on 

plot(256*[level2 level2],ylim,'r','LineWidth',2); %Κανονικοποίηση του level1 ως προς το ιστόγραμμα
legend({'Ιστόγραμμα','Κατώγλι κατά Otsu'},'Location','northwest')
title("Εμφάνιση του σημείου κατωφλιοποίησης κατά,Otsu,στην εικόνα 2")

hold off

figure(3)

counts=imhist(img_1);       %Ανάθεση του ιστογράμματος σε μια μεταβλητή
L=length(counts);           %Εναπόθεση της έκτασης της μεταβλητής counts στην μεταβλητή L
p=counts/sum(counts);       %Υπολογισμός της πιθανότητας 
Propability_1=cumsum(p);    %Υπολογισμός της πιθανότητας της κλάσης 1 

%Για τον υπολογισμό της μέσης έντασης μέχρι το επίπεδο του κατωφλίου
%(average intensity) 

m_k=cumsum(p.*(1:L));

%Για τον υπολογισμό της μέσης έντασης ολόκληρης της εικόνας (global mean)

m_g=m_k(end);

sigma_b_squared=(((m_g * Propability_1) - m_k).^2) ./ (Propability_1 .* (1-Propability_1));

yyaxis left
plot(counts);
ylabel('Ιστόγραμμα')

yyaxis right
plot(sigma_b_squared);
ylabel('\sigma_B^2');

xlim([1 256])

[~,k] = max(sigma_b_squared);

hold on
plot([k,k],'LineWidth',5)
hold off

Alsoo this is a page that I've found to help me implement the final part:https://blogs.mathworks.com/steve/2016/06/14/image-binarization-otsus-method/

Can you please help me?

Thank you in advance

P.S. There's another image that I want to process, so don't mind the second one

Answer 1

[Sudarsanan A K]

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:

https://mathworks.com/help/images/ref/otsuthresh.html

https://mathworks.com/help/images/ref/graythresh.html

I hope this helps!