Otsu 的阈值实现无法正常工作

Question

我已经实现了 otsu 的阈值实现，它将图像分割成前景和背景图像。我的实现的输出似乎与期望的不符。有什么想法吗？提前致谢！真的很感激有人能告诉我如何解决这个问题。

我的输出： 对于图像 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);

Answer 1

您的代码存在一些问题，我将概述其中的错误：

1. 我只是在吹毛求疵，但您可以使用 numel 来计算像素总数 (N)...这样更简洁:)
2. 在您的原始代码中，您正在检查 1 到 254 之间的正确阈值。您确实应该检查 0 到 255，因为图像中有 256 种可能的强度。
3. 您还需要更改您的 sigma 声明，以便有 256 个元素，而不是 255 个。请记住，您的图像中有 256 种可能的强度。
4. 在检查每个强度的 for 循环中，当您计算类出现的概率时，您还需要检查强度 0。由于 MATLAB 从 1 开始索引数组，因此您需要偏移访问索引以便从 1 开始。
5. 您对对象和背景之间差异的定义略有偏差。您还需要计算对象的类均值的概率。您可以查看代码以获取更多详细信息。
6. 您的类均值定义概率略有不准确。您需要除以qb 和qo，而不是N。
7. 您正在使用m 计算然后意味着您应该将其存储在muT 中。
8. 最后，当您找到对象和背景之间的最大方差时，您需要减去 1，因为这将提供介于 0 和 255 之间的强度。

因此，这就是您的代码的样子。请注意，我省略了阈值图像的代码。我只提供计算图像阈值的代码。

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 循环，我发现这段代码效率很低。我个人会做的是矢量化它，但我会把它留给你作为学习练习:)

---

编辑

好的，我会放弃的。这是我为 Otsu 编写的一些代码，它更加矢量化。这基本上执行了您在上面所做的事情，但是以更加矢量化的方式。非常欢迎您将它用于您自己的目的，但如果您当然打算使用它，请引用我：)

%// 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

---

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