算法解剖系列-Otsu二值化原理及实现

基本原理

• Otsu是以统计决策为基础的；

• Otsu方法在类间方差最大的情况下获取最佳阈值，达到最佳分割效果；（即，在目标与背景之间差别最大时获取的阈值）

• Otsu方法是以图像的相对直方图（归一化直方图）为基础对一维阵列进行计算的；

推导大致过程：

对一幅大小为$M\times{N}$的数字图像：

令$L$表示灰度级数；

$n_i$表示灰度级为$i$的像素数，则图像中像素总数$MN = n_1+n_2+n_3+···+n_L；$

$p_i = n_i/MN,p_i$为相对直方图，即灰度级为$i$出现的概率。则$\sum_{i = 1}^L p_i = 1$；

设阈值$T(k) = k,1<k<L可以将图像分为C_1和C_2两类，灰度值范围分别为[1,k]和[k+1,L]$。

• 则被分到$C_1和C_2$的概率分别为：
  $P_1(k) = \sum_{i = 1}^k p_i \tag {1}$
  $P_2(k) = \sum_{i = k+1}^L p_i = 1-P_1 (k)\tag {2}$

• 则被分到$C_1和C_2$像素的灰度均值分别为：
  $m_1(k) =\sum_{i = 1}^k i P(i/C_1) = \sum_{i = 1}^k i P(C_1/i)P(i)/P(C_1)= \frac{1}{P_1(k)} \sum_{i = 1}^k ip_i \tag {3}$
  $m_2(k) =\sum_{i = k+1}^L i P(i/C_2)= \frac{1}{P_2(k)} \sum_{i = k+1}^L ip_i \tag {4}$

• 直至$k$级的累加均值：
  $m(k) =\sum_{i = 1}^k ip_i \tag{5}$

• 图像的灰度均值（即全局均值）：
  $m_G =\sum_{i = 1}^L ip_i \tag {6}$

由$（3）和（4）$得：
$P_1m_1+P_2m_2 = m_G\tag{7}$
类间方差：
$\sigma_B^2 = P_1(m_1-m_G)^2 +P_2(m_2-m_G)^2\tag{8}$
将$（7）中的m_G代入（8）中，以及由P_1+P_2 = 1$得：

$$\sigma_B^2 = P_1P_2(m_1-m_2)^2\tag {9}$$
最后由$(3)(5)得m_1 = \frac{1}{P_1}m，由（7）得m_2 = \frac{m_G-P_1m}{P_2},以及由P_1+P_2 = 1得$：
$$\sigma_B^2 = \frac{(m_GP_1-m)^2}{P_1(P_1-1)}$$
所以要获取类间方差，只需要计算全局均值$m_G、k级的累加均值m和累加概率P_1$。

设$\sigma_G^2 代表全局方差，\eta$表示可分性度量，表示可分割性：
$$\sigma_G^2 =\sum_{i = 1}^L (i-m_G)^2p_i$$

$$\eta = \frac{\sigma_B^2}{\sigma_G^2 }$$

最终的要计算的是：

$$\sigma_B^2 （k）= \frac{(m_GP_1(k)-m(k))^2}{P_1(k)(P_1(k)-1)}\tag {10}$$
$$\eta (k^*)= \frac{\sigma_B^2(k^*)}{\sigma_G^2 }\tag{11}$$
$k^*是最佳阈值，即\sigma_B^2 （k）为最大值的k值。$

算法基本步骤

---

1. 计算输入图像的相对直方图:根据$p_i = n_i/MN$；

2. 计算像素被分到$C_1 累计和P_1(k) ,根据上式(1)$；

3. 计算累计均值$m(k),根据上式(5)$；

4. 计算全局灰度均值$m_G，根据上式(6)$；

5. 计算类间方差$\sigma_B^2(k) ，根据上式(10)$ ；

6. 获取最大的$\sigma_B^2(k)中的k$值，即为最佳阈值；

7. 计算可分性度量$\eta^* ，根据上式(11)$。

---

算法实现

matlab 代码

close all;
clear all;
clc;

input = imread('R.png');%读图
input = rgb2gray(input);%灰度转换

L = 256;%给定灰度级
[ni, li] = imhist(input,L);  %ni-各灰度等级出现的次数；li-对应的各灰度等级
% figure,plot(xi,ni);%显示绝对直方图统计
% title('绝对直方图统计')

[M,N] = size(input);%获取图像大小
MN = M*N;%像素点总数

%%Step1 计算归一化直方图

pi = ni/MN;  %pi-统计各灰度级出现的概率
figure,plot(li,pi);%显示相对直方图统计
title('相对直方图统计')

%%Step2  计算像素被分到C1中的概率P1(k)

sum = 0;%用来存储各灰度级概率和
P1 = zeros(L,1);%用来存储累积概率和
for k = 1:L
    sum = sum +pi(k,1);
    P1(k,1) = sum;%累加概率
end

%%Step3  计算像素至K级的累积均值m(k)

sum1 = 0;%用来存储灰度均值
m = zeros(L,1);%用来存储累计均值
for k = 1:L
    sum1 = sum1+k*pi(k,1);
    m(k,1) = sum1;%累积均值
end

%%Step4  计算全局灰度均值mg

mg = sum1;

%%Step5 计算类间方差sigmaB2 
sigmaB2 = zeros(L,1);
for k = 1:L
    if(P1(k,1) == 0)
        sigmaB2(k,1) = 0;  %为了防止出现NaN
    else
        sigmaB2(k,1) = ((mg*P1(k,1)-m(k,1))^2)/(P1(k,1)*(1-P1(k,1)));
    end
end

%%Step6 得到最大类间方差以及阈值

[MsigmaB2,T] = max(sigmaB2);%获取最大类间方差MsigmaB2，以及所在位置（即阈值）
output = zeros(M,N);%定义二值化输出图像
for i = 1:M
    for j = 1:N
        if input(i,j)>T
            output(i,j) = 1;
        else
            output(i,j)=0;
        end
    end
end
figure,imshow(output);%显示结果

%%Step7 可分性度量eta

sigmaG2 = 0;%全局方差
for k = 1:L
    sigmaG2 = sigmaG2+(k-mg)^2*pi(k,1);
end
eta = MsigmaB2/sigmaG2;

效果图