× {{alert.msg}} Never ask again
Get notified about new tutorials RECEIVE NEW TUTORIALS

Otsu's Thresholding Implementation not working properly

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 &gt;= 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>
comments powered by Disqus