Get notified about new tutorials
RECEIVE NEW TUTORIALS

First 15 Minutes Free

326 sessions given

since Jul 09, 2014

since Jul 09, 2014

Likelihood of Reply:
85%

Response Time:
within an hour

Ray Phan

Feb 02, 2015

<p>There are a few issues with your code, and I'll outline where it's wrong:</p>
<ol>
<li>I'm simply nitpicking, but you can count the total number of pixels (<code>N</code>) by using <a href="http://www.mathworks.com/help/matlab/ref/numel.html" rel="nofollow"><code>numel</code></a>... it's just cleaner :)</li>
<li>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. </li>
<li>You also need to change your <code>sigma</code> declaration so that there are 256 elements, not 255. Remember, there are 256 possible intensities in your image.</li>
<li>Within your <code>for</code> 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.</li>
<li>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.</li>
<li>Your probability of class mean definitions are slightly inaccurate. You need to divide by <code>qb</code> and <code>qo</code>, not <code>N</code>.</li>
<li>You are using <code>m</code> to calculate then mean when you should be storing it in <code>muT</code>.</li>
<li>Finally, when you find the maximum of the variance between the object and background, you need to <strong>subtract</strong> by 1, as this will provide an intensity between 0 and 255. </li>
</ol>
<p>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.</p>
<pre><code>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
</code></pre>
<hr>
<p>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 <code>for</code> loops. What I would personally do is vectorize it, but I'll leave that to you as a learning exercise :)</p>
<hr>
<h1>Edit</h1>
<p>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 :)</p>
<pre><code>%// 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
</code></pre>
<hr>
<h1>Edits by Divakar</h1>
<p>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 <code>s_b</code>:</p>
<pre><code>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);
</code></pre>
<p>This tip was originally posted on <a href="http://stackoverflow.com/questions/25829725/Otsu's%20Thresholding%20Implementation%20not%20working%20properly/25829961">Stack Overflow</a>.</p>