【老文新发】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]$
我们分别求出这个阈值前后的局部频率之和，定义如下：
C 0 = [ 1 , k ] C 1 = [ k + 1 , L ] w 0 = Pr ⁡ ( C 0 ) = ∑ i = 1 k p i = w ( k ) w 1 = Pr ⁡ ( C 1 ) = ∑ i = k + 1 L p i = 1 − w ( k )

C0​=[1,k]C1​=[k+1,L]w0​=Pr(C0​)=i=1∑k​pi​=w(k)w1​=Pr(C1​)=i=k+1∑L​pi​=1−w(k)​
则灰度图像频率直方图的总的数学期望和 C0 、C1的数学期望如下：
u 0 = ∑ i = 1 k i ∗ Pr ⁡ ( i ∣ C 0 ) = ∑ i = 1 k i ∗ p i / w 0 = u ( k ) w ( k ) u 1 = ∑ i = k + 1 L i ∗ Pr ⁡ ( i ∣ C 1 ) = ∑ i = k + 1 L i ∗ p i / w 1 = u T − u ( k ) 1 − w ( k ) u T = u ( L ) = ∑ i = 1 L i ∗ p i

u0=∑i=1ki∗Pr⁡(i∣C0)=∑i=1ki∗pi/w0=u(k)w(k)u1=∑i=k+1Li∗Pr⁡(i∣C1)=∑i=k+1Li∗pi/w1=uT−u(k)1−w(k)uT=u(L)=∑i=1Li∗pi" role="presentation" style="position: relative;">u0=∑i=1ki∗Pr(i∣C0)=∑i=1ki∗pi/w0=u(k)w(k)u1=∑i=k+1Li∗Pr(i∣C1)=∑i=k+1Li∗pi/w1=uT−u(k)1−w(k)uT=u(L)=∑i=1Li∗pi$$\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}$$

u0​=i=1∑k​i∗Pr(i∣C0​)=i=1∑k​i∗pi​/w0​=w(k)u(k)​u1​=i=k+1∑L​i∗Pr(i∣C1​)=i=k+1∑L​i∗pi​/w1​=1−w(k)uT​−u(k)​uT​=u(L)=i=1∑L​i∗pi​​
上式各个变量之间的关系如下：
$$w_0{u}_0+w_1{u}_1=u_T{w}_0+w_1=1$$
期望有了，计算一下对应的方差：
σ 0 2 = ∑ i = 1 k ( i − u 0 ) 2 Pr ⁡ ( i ∣ C 0 ) = ∑ i = 1 k ( i − u 0 ) 2 p i / w 0 σ 1 2 = ∑ i = k + 1 L ( i − u 1 ) 2 Pr ⁡ ( i ∣ C 1 ) = ∑ i = k + 1 L ( i − u 1 ) 2 p i / w 0 σ T 2 = ∑ i = 1 L ( i − u T ) 2 p i

σ02=∑i=1k(i−u0)2Pr⁡(i∣C0)=∑i=1k(i−u0)2pi/w0σ12=∑i=k+1L(i−u1)2Pr⁡(i∣C1)=∑i=k+1L(i−u1)2pi/w0σT2=∑i=1L(i−uT)2pi" role="presentation" style="position: relative;">σ20=∑i=1k(i−u0)2Pr(i∣C0)=∑i=1k(i−u0)2pi/w0σ21=∑i=k+1L(i−u1)2Pr(i∣C1)=∑i=k+1L(i−u1)2pi/w0σ2T=∑i=1L(i−uT)2pi$$\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}$$

σ02​=i=1∑k​(i−u0​)2Pr(i∣C0​)=i=1∑k​(i−u0​)2pi​/w0​σ12​=i=k+1∑L​(i−u1​)2Pr(i∣C1​)=i=k+1∑L​(i−u1​)2pi​/w0​σT2​=i=1∑L​(i−uT​)2pi​​
根据文献：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相关的函数
σ b 2 = w 0 ( u 0 − u T ) 2 + w 1 ( u 1 − u T ) 2 = w 0 ( u 0 − ( w 0 u 0 + w 1 u 1 ) ) 2 + w 1 ( u 1 − ( w 0 u 0 + w 1 u 1 ) ) 2 = w 0 w 1 ( u 1 − u 0 ) 2

σb2​=w0​(u0​−uT​)2+w1​(u1​−uT​)2=w0​(u0​−(w0​u0​+w1​u1​))2+w1​(u1​−(w0​u0​+w1​u1​))2=w0​w1​(u1​−u0​)2​
然后再分别把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
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29

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

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

多分类最大类间方差法

根据以上公式类推到多分类的最大类间方差法，假设有 $\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

{t1​∗,t2​∗,…,tM−1​∗}=​ArgMax{σB​2(t1​,t2​,…,tM−1​)},1≤t1​<…<tM−1​<L​

成立
其中:
σ B 2 = ∑ k = 1 M ω k ( μ k − μ T ) 2 ω k = ∑ i ∈ C k p i , μ k = ∑ i ∈ C k i p i / ω ( k ) .

σB​2ωk​μk​​=k=1∑M​ωk​(μk​−μT​)2=i∈Ck∑​pi​,=i∈Ck∑​ipi​/ω(k).​

因为
∑ k = 1 M ω k = 1 μ T = ∑ k = 1 M ω k μ k .

k=1∑M​ωk​=1μT​=k=1∑M​ωk​μk​.​
因此
$${\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

{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
以上这两篇博文讲的也很详细，内容有所参考，欢迎访问阅读。