Otsu's Thresholding Implementation not working properly

Question

I have implemented otsu's threshold implementation which segments an image into foreground and background image.The output of my implementation seems off from the desired one. Any thoughts on it? Thanks in advance! Would really appreciate it someone can tell me how do i resolve the issue.

My output: For image 1-

For image 2-

My Code:

im1=imread('D:\root-image.pgm');
% im1=rgb2gray(im1);
[n,m]=size(im1);
 hst=imhist(im1);
mu=zeros(255,1);
N=0;
for i=1:255
    N=N+hst(i);
end
% The total mean level of the original image
for i=1:255
    mu(i)=mu(i)+((i.*hst(i))./N);
end


for T=1:254
    qb=0;
    muT=0;
    qo=0;
    for i=1:T
        qb=qb+(hst(i)./N); % probability of class occurence (background)
        m=m+((i.*hst(i))./N);% probability of class mean (background)
    end
    for i=T+1:255
        qo=qo+(hst(i)./N);% probability of class occurence (object)
    end
    sigma(T)=((mu(T)-(qb*muT))^2)/(qb*qo);
end
[Y,T] = max(sigma);

[n,m]=size(im1);
for i=1:n
    for j=1:m
        if im1(i,j)>T
            im(i,j)=1;
        else
            im(i,j)=0;
        end
    end
end
figure(1);
subplot(1,2,1);
imshow(im1);
% subplot(1,3,2);
% imhist(im1);
subplot(1,2,2);
imshow(im);

Answer 1

There are a few issues with your code, and I'll outline where it's wrong:

1. I'm simply nitpicking, but you can count the total number of pixels (N) by using numel... it's just cleaner :)
2. In your original code, you are checking for the right threshold between 1 and 254. You really should be checking from 0 to 255, as there are 256 possible intensities in your image.
3. You also need to change your sigma declaration so that there are 256 elements, not 255. Remember, there are 256 possible intensities in your image.
4. Within your for loop for checking each intensity, when you are calculating the probability of class occurrences, you need to check for intensity 0 as well. Because of the fact that MATLAB starts indexing arrays at 1, you'll need to offset your access index so that you're starting at 1.
5. Your definition of the variance between the object and background is slightly off. You need to also calculate the probability of the class mean for the object as well. You can check the code for more details.
6. Your probability of class mean definitions are slightly inaccurate. You need to divide by qb and qo, not N.
7. You are using m to calculate then mean when you should be storing it in muT.
8. Finally, when you find the maximum of the variance between the object and background, you need to subtract by 1, as this will provide an intensity between 0 and 255.

As such, this is what your code looks like. Note that I have omitted the code that thresholds your image. I'm only providing the code that calculates the threshold for your image.

hst=imhist(im1);
sigma = zeros(256,1); %// Change
N = numel(im1); %// Change

for T=0:255 %// Change
    qb=0;
    muT=0;
    qo=0;
    muQ=0; %// Change
    for i=0:T %// Change
        qb=qb+(hst(i+1)./N); % probability of class occurence (background)
    end    
    for i=T+1:255 
        qo=qo+(hst(i+1)./N);% probability of class occurence (object)
    end
    for i=0:T%// Change        
        muT=muT+((i.*hst(i+1))./qb);% probability of class mean (background)    
    end
    for i=T+1:255 %// Change
        muQ=muQ+((i.*hst(i+1))./qo);% probability of class mean (object)        
    end   
    sigma(T+1) = qb*qo*((muT-muQ)^2); %// Change
end
[Y,T] = max(sigma);
T = T-1; %// Change - For 0 to 255

---

This code should now work. I ran this code with my own implementation of Otsu, and I get the same calculated threshold. To be honest, I find this code to be rather inefficient due to the many for loops. What I would personally do is vectorize it, but I'll leave that to you as a learning exercise :)

---

Edit

OK, I'll give in. Here's some code that I wrote for Otsu that is more vectorized. This essentially performs what you are doing above, but in a more vectorized manner. You're more than welcome to use it for your own purposes, but do cite me if you intend to use it of course :)

%// Function that performs Otsu's adaptive bimodal thresholding
%// Written by Raymond Phan - Version 1.0
%// Input - im - Grayscale image
%// Output - out - Thresholded image via Otsu

function [out] = otsu(im)

%// Compute histogram
J = imhist(im);

%// Total number of intensities
L = length(J);

%// Some pre-processing stuff
%// Determine total number of pixels
num_pixels = sum(J(:));

%// Determine PDF
pdf = J / num_pixels;

%// Storing between-class variances for each intensity
s_b = zeros(1,L);

for idx = 1 : L
    %// Calculate w_0
    w_0 = sum(pdf(1:idx));

    %// Calculate w_1
    w_1 = 1 - w_0;

    %// Calculate u_0
    u_0 = sum((0:idx-1)'.*pdf(1:idx)) / w_0;

    %// Calculate u_1
    u_1 = sum((idx:L-1)'.*pdf(idx+1:L)) / w_1;

    % // Calculate \sigma_b^{2}
    s_b(idx) = w_0*w_1*((u_1 - u_0)^2);    
end

%// Find intensity that provided the biggest variance
[max_s_b T] = max(s_b);
T = T - 1; %// Must subtract by 1, since the index starts from 1

%// Now create an output image that thresholds the input image
out = im >= T;

end

---

Edits by Divakar

Divakar (thanks!) has created vectorized code to replace the loop portion of the above function code and this essentially gets rid of the pre-allocation for s_b:

w_0 = cumsum(pdf);
w_1 = 1 - w_0;
u_0 = cumsum((0:L-1)'.*pdf)./w_0;
u_1 = flipud([0 ; cumsum((L-1:-1:1)'.*pdf((L:-1:2)))])./w_1;
s_b = w_0.*w_1.*((u_1 - u_0).^2);