Otsu的阈值实施不能正常工作 [英] Otsu's Thresholding Implementation not working properly

问题描述

我已经实现了otsu的阈值实现，将图像分割为前景和背景图像。我的实现的输出似乎偏离了所需的。有什么想法吗？提前致谢！非常感谢有人能告诉我如何解决这个问题。

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.

我的输出：
对于图片1 -

My output: For image 1-

>

对于图片2 -

我的代码：

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);

推荐答案

您的代码有几个问题，

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

1. 我只是简单地选择，但你可以计算总像素数（ N ），请使用 numel ...它只是更干净：）


2. 在原始代码中，您正在检查1到254之间的正确阈值。从0到255检查，因为您的图像中有256种可能的亮度。


3. 您还需要更改 sigma 声明，以便有256个元素，而不是255.记住，256个可能的强度您的图片。


4. 在您的中为循环检查每个强度，当您计算类出现的概率时，您需要检查对于强度0。因为MATLAB在1开始索引数组，所以您需要偏移访问索引，以便从1开始。


5. 您定义的对象和背景稍微偏离。您还需要计算对象的类平均值的概率。您可以检查代码以了解更多详情。


6. 您的类均值定义的概率略有不准确。您需要除以 qb 和 qo ，而不是 N 。


7. 您现在使用 m 来计算当您应该将其存储在 muT 。


8. 最后，当您找到对象和背景之间的最大差异时，您需要将减1 强度介于0到255之间。

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

---

代码现在工作。我运行这个代码与我自己的实现Otsu，我得到相同的计算阈值。说实话，我发现这个代码是非常低效的，因为许多 for 循环。我个人会做的是将它矢量化，但我会留给你作为一个学习练习：）

---

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 :)

好的，我会赠送。这里是一些我为Otsu写的更加矢量化的代码。这基本上执行你在上面做的，但是以更矢量化的方式。你非常欢迎使用它为您自己的目的，但如果你打算使用它，当然：）

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

---

by Divakar

Divakar（感谢！）创建了用于替换上述函数代码的循环部分的矢量化代码，并且这基本上摆脱了 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); 

这篇关于Otsu的阈值实施不能正常工作的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！