【老文新发】Otsu大津法详解及python实现

原文：A Threshold Selection Method from Gray-Level Histograms
A Fast Algorithm for Multilevel Thresholding

前言

大津法包含两个重要的概念：类间方差（between-class variance）和类内方差（within-class variance）

两者的详细关系推导可后文。

大津法又称为最大类间方差法是有原因的。因为这个算法的目的就是最大化类间方差，且这个最优阈值一定存在。

大津法作为阈值自动分割的经典算法，其思想很巧妙，值得学习。

大津法推导

如图所示，右边分割分明，其类间的差异大，区分明显，所以其类间方差更大。

大津法就是要实现这个过程。

我们先做图像的直方图，统计每个很小像素区间包含的像素个数。

即将图像像素值分为 [ 1 , 2 , 3 , …. , L] 个区间。用$n_i$ 表示各个水平像素值的像素个数，总像素个数与$n_i$关系为：
$$N=n_1+n_2+\dots +n_L$$
像元个数 $ni$ 比上总像元数 N 即可得到某个像素区间出现的频率，定义 $pi$
$$p_i=n_i/N_i,p_i\ge 0,\sum _{i=1}^L{p}_i=1$$
根据原文，是利用一个阈值 k 把图像分为两类$C0、C1。k\in [1,2,3,\dots .,L]$
我们分别求出这个阈值前后的局部频率之和，定义如下：
$$\begin{array}{c}{C}_0=[1,k]\\ C_1=[k+1,L]\\ w_0=\mathrm{Pr}\left(C_0\right)=\sum _{i=1}^k{p}_i=w(k)\\ w_1=\mathrm{Pr}\left(C_1\right)=\sum _{i=k+1}^L{p}_i=1-w(k)\end{array}$$
则灰度图像频率直方图的总的数学期望和 C0 、C1的数学期望如下：
$$\begin{array}{c}{u}_0=\sum _{i=1}^{k}i\ast \mathrm{Pr}\left(i\mid C_0\right)=\sum _{i=1}^{k}i\ast p_i/w_0=\frac{u(k)}{w(k)}\\ u_1=\sum _{i=k+1}^{L}i\ast \mathrm{Pr}\left(i\mid C_1\right)=\sum _{i=k+1}^{L}i\ast p_i/w_1=\frac{u_T-u(k)}{1-w(k)}\\ u_T=u(L)=\sum _{i=1}^{L}i\ast p_i\end{array}$$
上式各个变量之间的关系如下：
$$w_0{u}_0+w_1{u}_1=u_T{w}_0+w_1=1$$
期望有了，计算一下对应的方差：
$$\begin{array}{c}{\sigma }_0^2=\sum _{i=1}^k{\left(i-u_0\right)}^2\mathrm{Pr}\left(i\mid C_0\right)=\sum _{i=1}^k{\left(i-u_0\right)}^2{p}_i/w_0\\ {\sigma }_1^2=\sum _{i=k+1}^L{\left(i-u_1\right)}^2\mathrm{Pr}\left(i\mid C_1\right)=\sum _{i=k+1}^L{\left(i-u_1\right)}^2{p}_i/w_0\\ {\sigma }_T^2=\sum _{i=1}^L{\left(i-u_T\right)}^2{p}_i\end{array}$$
根据文献：Introduction to statistical pattern recognition，260-267。类内误差、类间误差、总误差有如下关系：
$$类内误差{\sigma }_w^2=w_0{\sigma }_0^2+w_1{\sigma }_1^2$$
$$类间误差{\sigma }_b^2=w_0{\left(u_0-u_T\right)}^2+w_1{\left(u_1-u_T\right)}^2$$

$$总误差{\sigma }_w^2+{\sigma }_b^2={\sigma }_T^2$$
注意：总误差是与阈值k无关的，但类间误差和类内误差是与阈值k相关的函数
$$\begin{array}{c}{\sigma }_b^2=w_0{\left(u_0-u_T\right)}^2+w_1{\left(u_1-u_T\right)}^2\\ =w_0{\left(u_0-\left(w_0{u}_0+w_1{u}_1\right)\right)}^2+w_1{\left(u_1-\left(w_0{u}_0+w_1{u}_1\right)\right)}^2\\ =w_0{w}_1{\left(u_1-u_0\right)}^2\end{array}$$
然后再分别把w0,u1,u0带入，可得到：
$${\sigma }_b^2=\frac{{\left[u_{T}w(k)-u(k)\right]}^2}{w(k)[1-w(k)]}$$

则求解最大类间误差为:
$${\sigma }_b^2\left(k^{\ast }\right)=\underset{1\le k<L}{\max }{\sigma }_b^2(k)$$

由上述 ${\sigma }_b^2$ 的分母可以发现， $w(k)$ 可以取到 1 也可以取到0，因此在边界上 ${\sigma }_b^2$ 可以无穷大，而在开基 $(0,1)$ 则类间方差有限，因此在定义域
$$S^{\ast }=k:w_0{w}_1=w(k)[1-w(k)]>0$$
因此必定存在一个阈值k使得两类类间方差最大。
以下是python代码实现：

def otsu(gray_img):
    n_count = gray_img.size

    gray_img_array = gray_img.flatten()
    index = np.flatnonzero(gray_img_array)
    gray_img_data = gray_img_array[index]
    
    threshold_t = 0
    max_g = 0
    
    t = np.linspace(start=-1, stop=1, num=256)
    # 遍历每一个灰度层
    for i in range(len(t)):
    	# 使用numpy直接对数组进行运算
        n0 = gray_img_data[np.where(gray_img_data < t[i])]
        n1 = gray_img_data[np.where(gray_img_data >= t[i])]
        w0 = len(n0) / n_count
        w1 = len(n1) / n_count
        u0 = np.mean(n0) if len(n0) > 0 else 0.
        u1 = np.mean(n1) if len(n0) > 0 else 0.
        
        g = w0 * w1 * (u0 - u1) ** 2
        if g > max_g:
            max_g = g
            threshold_t = t[i]
    print('类间方差最大阈值：', threshold_t)
    gray_img[gray_img < threshold_t] = 0
    gray_img[gray_img >= threshold_t] = 1
    return gray_img

这个在opencv中已经有实现，可以直接调用

import cv2
t, otsu = cv2.threshold(img, 0, 255, cv2>THRESH_BINARY + cv2.THRESH_OTSU)

多分类最大类间方差法

根据以上公式类推到多分类的最大类间方差法，假设有 $\mathrm{m}-1$ 个阈值 { t 1 , t 2 , … , t M − 1 } \{\mathrm{t} 1, \mathrm{t} 2, \ldots, \mathrm{tM}-1\} {t1,t2,…,tM−1} 将图像分为 $\mathrm{M}$ 类， $C_1$ ， $C_2\dots C_M$ 。 则存在一组阈值 { t 1 ∗ , t 2 ∗ , … , t M − 1 ∗ } \left\{t 1^*, \mathrm{t} 2^*, \ldots, \mathrm{tM}-1^*\right\} {t1∗,t2∗,…,tM−1∗} 使得
{ t 1 ∗ , t 2 ∗ , … , t M − 1 ∗ } = Arg ⁡ Max ⁡ { σ B 2 ( t 1 , t 2 , … , t M − 1 ) } , 1 ≤ t 1 < … < t M − 1 < L \begin{aligned} \left\{\mathrm{t}_1 *, \mathrm{t}_2 *, \ldots, \mathrm{t}_{\mathrm{M}-1} *\right\}= & \operatorname{Arg} \operatorname{Max}\left\{\sigma_{\mathrm{B}}{ }^2\left(\mathrm{t}_1, \mathrm{t}_2, \ldots, \mathrm{t}_{\mathrm{M}-1}\right)\right\}, \\ & 1 \leq \mathrm{t}_1<\ldots<\mathrm{t}_{\mathrm{M}-1}<\mathrm{L} \end{aligned} {t1​∗,t2​∗,…,tM−1​∗}=​ArgMax{σB​2(t1​,t2​,…,tM−1​)},1≤t1​<…<tM−1​<L​

成立
其中:
$$\begin{array}{rl}{\sigma }_{\mathrm{B}}{}^2& =\sum _{k=1}^{\mathrm{M}}{\omega }_{\mathrm{k}}{\left({\mu }_{\mathrm{k}}-{\mu }_{\mathrm{T}}\right)}^2\\ {\omega }_{\mathrm{k}}& =\sum _{\mathrm{i}\in \mathrm{Ck}}{\mathrm{p}}_{\mathrm{i}},\\ {\mu }_{\mathrm{k}}& =\sum _{\mathrm{i}\in \mathrm{Ck}}\mathrm{i}{\mathrm{p}}_{\mathrm{i}}/\omega (\mathrm{k}).\end{array}$$

因为
$$\begin{array}{c}\sum _{k=1}^{\mathrm{M}}{\omega }_{\mathrm{k}}=1\\ {\mu }_{\mathrm{T}}=\sum _{k=1}^{\mathrm{M}}{\omega }_{\mathrm{k}}{\mu }_{\mathrm{k}}.\end{array}$$
因此
$${\sigma }_{\mathrm{B}}{}^2\left({\mathrm{t}}_1,{\mathrm{t}}_2,\dots ,{\mathrm{t}}_{\mathrm{M}-1}\right)=\sum _{k=1}^{\mathrm{M}}{\omega }_{\mathrm{k}}{\mu }_{\mathrm{k}}^2-{\mu }_{\mathrm{T}}{}^2$$

$\mu \mathrm{T}$ 与间值无关，因此求上式的最大值可转为:
{ t 1 ∗ , t 2 ∗ , … , t M − 1 ∗ } = Arg ⁡ Max ⁡ { ( σ B ′ ) 2 { { t 1 , t 2 , … , t M − 1 } } 1 ≤ t 1 < … < t M − 1 < L ( σ B ) 2 = ∑ k = 1 M ω k μ k 2 \begin{array}{c} \left\{\mathrm{t}_{1}^{*}, \mathrm{t}_{2}^{*}, \ldots, \mathrm{t}_{\mathrm{M}-1}^{*}\right\}=\operatorname{Arg} \operatorname{Max}\left\{\left(\sigma_{\mathrm{B}}{ }^{\prime}\right)^{2}\left\{\left\{\mathrm{t}_{1}, \mathrm{t}_{2}, \ldots, \mathrm{t}_{\mathrm{M}-1}\right\}\right\}\right. \\ 1 \leq \mathrm{t}_{1}<\ldots<\mathrm{t}_{\mathrm{M}-1}<\mathrm{L} \\ \left(\sigma_{\mathrm{B}}\right)^{2}=\sum_{k=1}^{\mathrm{M}} \omega_{\mathrm{k}} \mu_{\mathrm{k}}{ }^{2} \end{array} {t1∗​,t2∗​,…,tM−1∗​}=ArgMax{(σB​′)2{{t1​,t2​,…,tM−1​}}1≤t1​<…<tM−1​<L(σB​)2=∑k=1M​ωk​μk​2​

类间方差、类内方差和总方差关系

https://blog.csdn.net/qq_42164483/article/details/119064535
https://blog.csdn.net/m0_38024332/article/details/104226806
以上这两篇博文讲的也很详细，内容有所参考，欢迎访问阅读。