阈值分割：最大类间方差法

一、简介

最大类间方差法，又称为大津阈值法，或OTSU算法。是由日本学者大津在197919791979年提出的一种非参数的、无监督的自动选择阈值的图像分割方法。

二、算法描述

2.1 公式推导

对于给定的一幅具有 LLL 个灰度级（[0,1,2,⋯,L−1][0,1,2,\cdots,L-1][0,1,2,⋯,L−1]）的灰度图像，有以下描述：

使用 nin_{i}ni 表示处于灰度级 iii 的像素块的数目；
使用 N=n0+n1+⋯+nL−1N = n_0 + n_1 + \cdots + n_{L-1}N=n0+n1+⋯+nL−1 表示该图像所有像素块的数目总和。

对该图像的灰度级直方图进行 标准化处理 可以得到一个概率分布：
pi=niN，pi≥0,∑i=0L−1pi=1(1)p_i = \frac{n_i}{N}，\quad p_i \ge 0,\quad \sum_{i=0}^{L-1}{p_i} =1 \tag{1} pi=Nni，pi≥0,i=0∑L−1pi=1(1)

以灰度级 kkk 为阈值，可将该图像中的像素块二分为两大类 C0C_{0}C0 和 C1C_{1}C1 （背景和前景，反之亦然），其中：

C0C_0C0 表示灰度级在范围 [0,1,⋯,k−1][0,1,\cdots,k-1][0,1,⋯,k−1] 的像素集合；
C1C_1C1 表示灰度级在范围 [k,k+1,⋯,L−1][k,k+1,\cdots,L-1][k,k+1,⋯,L−1] 的像素集合。

记 w0w_0w0 与 w1w_1w1 分别为类 C0C_0C0 和类 C1C_1C1 发生的概率，则显然有：
∑i=1L−1pi=w0+w1=1(2)\sum_{i = 1}^{L-1}{p_i} = w_0 + w_1 = 1 \tag{2} i=1∑L−1pi=w0+w1=1(2)

其中，w0w_0w0 与 w1w_1w1 可表示为 零阶矩 的形式：
w0=Pr(C0)=∑i=0k−1pi=w(k)(3)w_0 = Pr(C_0) = \sum_{i = 0}^{k-1}{p_{i}} = w(k) \tag{3} w0=Pr(C0)=i=0∑k−1pi=w(k)(3)

w1=Pr(C1)=∑i=kL−1pi=1−w(k)(4)w_1 = Pr(C_1) = \sum_{i = k}^{L-1}{p_{i}} = 1 - w(k) \tag{4} w1=Pr(C1)=i=k∑L−1pi=1−w(k)(4)

类 C0C_0C0 与类 C1C_1C1 发生的期望（或均值）可表示为 一阶矩 的形式：
μ0=∑i=0k−1iPr(i∣C0)=∑i=0k−1ipiw0=1w0∑i=0k−1ipi=μ(k)w(k)(5)\mu_0 = \sum_{i=0}^{k-1}{iPr(i|C_0)} = \sum_{i=0}^{k-1}{i \frac{p_i}{w_0}} = \frac{1}{w_0}\sum_{i = 0}^{k-1}{ip_i} = \frac{\mu(k)}{w(k)} \tag{5} μ0=i=0∑k−1iPr(i∣C0)=i=0∑k−1iw0pi=w01i=0∑k−1ipi=w(k)μ(k)(5)

μ1=∑i=kL−1iPr(i∣C1)=∑i=kL−1ipiw1=1w1∑i=kL−1ipi=μT−μ(k)1−w(k)(5)\mu_1 = \sum_{i=k}^{L-1}{iPr(i|C_1)} = \sum_{i=k}^{L-1}{i \frac{p_i}{w_1}} = \frac{1}{w_1}\sum_{i = k}^{L-1}{ip_i}=\frac{\mu_{T} - \mu(k)}{1 - w(k)} \tag{5} μ1=i=k∑L−1iPr(i∣C1)=i=k∑L−1iw1pi=w11i=k∑L−1ipi=1−w(k)μT−μ(k)(5)

其中：
μ(k)=∑i=0k−1ipi,μT=∑i=0L−1ipi(6)\mu(k) = \sum_{i = 0}^{k-1}{ip_i},\quad\mu_{T} = \sum_{i = 0}^{L-1}{ip_i} \tag{6} μ(k)=i=0∑k−1ipi,μT=i=0∑L−1ipi(6)

显然，有：
μ0w0+μ1w1=μT(7)\mu_0 w_0 + \mu_1 w_1 =\mu_{T} \tag{7} μ0w0+μ1w1=μT(7)

类 C0C_0C0 与类 C1C_1C1 发生的方差可以表示为 二阶矩 的形式：
σ02=∑i=0k−1(i−μ0)2Pr(i∣C0)=∑i=0k−1(i−μ0)2piw0(8)\sigma_{0}^{2} = \sum_{i = 0}^{k-1}{(i - \mu_{0})^{2}Pr(i|C_0)} = \sum_{i = 0}^{k-1}{(i - \mu_0)^2\frac{p_i}{w_0}} \tag{8} σ02=i=0∑k−1(i−μ0)2Pr(i∣C0)=i=0∑k−1(i−μ0)2w0pi(8)

σ12=∑i=kL−1(i−μ1)2Pr(i∣C1)=∑i=kL−1(i−μ1)2piw1(9)\sigma_1^{2} = \sum_{i = k}^{L-1}{(i - \mu_1)^{2}Pr(i|C_1)} = \sum_{i = k}^{L-1}{(i - \mu_1)^{2}\frac{p_i}{w_1}} \tag{9} σ12=i=k∑L−1(i−μ1)2Pr(i∣C1)=i=k∑L−1(i−μ1)2w1pi(9)

为了评估所选阈值 kkk 的优良（goodness），引入以下三种判别标准度量：
λ=σB2σW2,κ=σT2σW2,η=σB2]σT2(10)\lambda = \frac{\sigma_{B}^{2}}{\sigma_{W}^{2}},\quad \kappa = \frac{\sigma_{T}^{2}}{\sigma_{W}^{2}},\quad \eta = \frac{\sigma_{B}^{2]}}{\sigma_{T}^{2}} \tag{10} λ=σW2σB2,κ=σW2σT2,η=σT2σB2](10)

其中：
（1）类内方差（within-class variance，简记为 σwith2\sigma_{with}^{2}σwith2 或 σW2\sigma_{W}^{2}σW2）满足：
σW2=w0σ02+w1σ12(11)\sigma_{W}^{2} = w_0 \sigma_{0}^{2} + w_1 \sigma_{1}^{2} \tag{11} σW2=w0σ02+w1σ12(11)

（2）类间方差（between-class variance，简记为 σBetween2\sigma_{Between}^{2}σBetween2 或 σB2\sigma_{B}^{2}σB2）满足：
σB2=w0(μ0−μT)2+w1(μ1−μT)2=w0w1(μ1−μ0)2(12)\sigma_{B}^{2} = w_{0}(\mu_{0} - \mu_{T})^{2} + w_1(\mu_{1} - \mu_{T})^{2} = w_{0}w_{1}(\mu_{1} - \mu_{0})^{2} \tag{12} σB2=w0(μ0−μT)2+w1(μ1−μT)2=w0w1(μ1−μ0)2(12)

（3）全局方差（total variance，简记为 σTotal2\sigma_{Total}^{2}σTotal2 或 σT2\sigma_{T}^{2}σT2）满足：
σT2=∑i=0L−1(i−μT)2pi(13)\sigma_{T}^{2} = \sum_{i = 0}^{L-1}{(i - \mu_{T})^{2}p_{i}} \tag{13} σT2=i=0∑L−1(i−μT)2pi(13)

合适的阈值会将图像分割为两类。反过来就是说，能在灰度水平上实现最佳分离的阈值将是最合适的阈值。

因此，利用引入的判别标准度量，可以将问题转化为一个优化问题：寻找一个合适的阈值 kkk 最大化某一个判别标准度量函数（λ\lambdaλ，κ\kappaκ 以及 η\etaη 中的某一个）。

那么选择哪一个判别标准度量最为合适呢？

实际上，最大化判别标准 λ\lambdaλ，κ\kappaκ 以及 η\etaη 是相互等价的。
因为，以 λ\lambdaλ 为单位可以分别表示 κ\kappaκ 以及 η\etaη：
κ=λ+1,η=λλ+1(14)\kappa = \lambda + 1,\quad \eta = \frac{\lambda}{\lambda + 1} \tag{14} κ=λ+1,η=λ+1λ(14)

并且以下基本关系式始终成立：
σW2+σB2=σT2(15)\sigma_{W}^{2} + \sigma_{B}^{2} = \sigma_{T}^{2} \tag{15} σW2+σB2=σT2(15)

也就是说，当三个判别标准度量中的任意一个达到最大时，另外两个都会达到最大值。其中：

σW2\sigma_{W}^{2}σW2 和 σB2\sigma_{B}^{2}σB2 均是阈值 kkk 的函数；
σT2\sigma_{T}^{2}σT2 与阈值 kkk 无关；
σW\sigma_{W}σW 需要计算二阶矩，而 σB2\sigma_{B}^{2}σB2 仅仅需要计算一阶矩。

所以，η\etaη 是关于阈值 kkk 的最简单的判别标准度量。

因此，使用 η\etaη 作为评估阈值 kkk 优良的判别标准。

现在，求解最大化判别标准 η\etaη 时的最佳阈值 k∗k^{*}k∗。分析 η=σB2σT2\eta = \frac{\sigma_{B}^{2}}{\sigma_{T}^{2}}η=σT2σB2 可知，最大化 η\etaη 即最大化 σB2\sigma_{B}^{2}σB2。即：
η(k)=σB2(k)σT2(16)\eta(k) = \frac{\sigma_{B}^{2}(k)}{\sigma_{T}^{2}} \tag{16} η(k)=σT2σB2(k)(16)

并且有：
σB2(k)=[μTw(k)−μ(k)]2w(k)[1−w(k)](17)\sigma_{B}^{2}(k) = \frac{\left[ \mu_{T} w(k) - \mu(k)\right]^{2}}{w(k)\left[1 - w(k)\right]} \tag{17} σB2(k)=w(k)[1−w(k)][μTw(k)−μ(k)]2(17)

最优阈值 k∗k^{*}k∗ 可以表示为：
k∗=arg⁡max⁡1≤k≤L−1σB2(k)(18)k^{*} = \arg \underset{1 \le k \le L-1}{\max}{\sigma_{B}^{2}{(k)}} \tag{18} k∗=arg1≤k≤L−1maxσB2(k)(18)

其中，阈值 kkk 的搜索范围可以表示为：
S∗={k∣w0w1=w(k)[1−w(k)]>0,or0<w(k)<1}(19)S^{*} = \{k~|~ w_0 w_1 = w(k)\left[1 - w(k)\right] > 0, ~ o r~ 0<w(k)<1\} \tag{19} S∗={k ∣ w0w1=w(k)[1−w(k)]>0, or 0<w(k)<1}(19)
该范围称为灰度直方图的有效范围。

从公式（12）中
σB2=w0(μ0−μT)2+w1(μ1−μT)2=w0w1(μ1−μ0)2\sigma_{B}^{2} = w_{0}(\mu_{0} - \mu_{T})^{2} + w_1(\mu_{1} - \mu_{T})^{2} = w_{0}w_{1}(\mu_{1} - \mu_{0})^{2} σB2=w0(μ0−μT)2+w1(μ1−μT)2=w0w1(μ1−μ0)2
可以看出：

当选定的阈值 k∈S−S∗={k∣w(k)=0or1}k \in S - S^{*} = \{k ~|~w(k) = 0~or~1\}k∈S−S∗={k ∣ w(k)=0 or 1} 时，判别标准度量 η\etaη 取最小值 000。
当选定的阈值 k∈S∗k \in S^{*}k∈S∗ 时，判别标准度量 η\etaη 取一个正的且有界的值。

因此，显而易见，判别标准的最大值始终存在。

2.2 算法分析

对于选定的阈值 k∗k^{*}k∗：
w0∗=Pr(C0∗)=∑i=0k∗−1pi=w(k∗)(20)w_{0}^{*} = Pr(C_{0}^{*}) = \sum_{i = 0}^{k^{*} - 1}{p_i} = w(k^{*}) \tag{20} w0∗=Pr(C0∗)=i=0∑k∗−1pi=w(k∗)(20)

w1∗=Pr(C1∗)=∑i=k∗L−1pi=1−w(k∗)(21)w_{1}^{*} = Pr(C_1^{*}) = \sum_{i = k^{*}}^{L-1}{p_{i}} = 1 - w(k^{*}) \tag{21} w1∗=Pr(C1∗)=i=k∗∑L−1pi=1−w(k∗)(21)

分别表示了灰度图像按照阈值 k∗k^{*}k∗ 所划分的两类的发生概率。

类 C0∗C_{0}^{*}C0∗ 与类 C1∗C_{1}^{*}C1∗ 发生的期望分别为：
μ0∗=∑i=0k∗−1iPr(i∣C0∗)=∑i=0k∗−1ipiw0∗=μ(k∗)w(k∗)(22)\mu_{0}^{*} = \sum_{i = 0}^{k^{*}-1}{iPr(i|C_{0}^{*})} = \sum_{i = 0}^{k^{*} - 1}{i \frac{pi}{w_{0}^{*}}} = \frac{\mu(k^{*})}{w(k^{*})} \tag{22} μ0∗=i=0∑k∗−1iPr(i∣C0∗)=i=0∑k∗−1iw0∗pi=w(k∗)μ(k∗)(22)
μ1∗=∑i=k∗L−1iPr(i∣C1∗)=∑i=k∗L−1ipiw1∗=μT−μ(k∗)1−w(k∗)(23)\mu_{1}^{*} = \sum_{i = k^{*}}^{L-1}{iPr(i|C_{1}^{*})} = \sum_{i = k^{*}}^{L-1}{i\frac{p_i}{w_1^{*}}} = \frac{\mu_{T} - \mu(k^{*})}{1 - w(k^{*})} \tag{23} μ1∗=i=k∗∑L−1iPr(i∣C1∗)=i=k∗∑L−1iw1∗pi=1−w(k∗)μT−μ(k∗)(23)

将判别标准 η\etaη 的最大值 η(k∗)\eta(k^{*})η(k∗) 简记为 η∗\eta^{*}η∗，可以用作评估灰度图像中类的可分性的标准。这是一个重要的度量，它在灰度尺度的放射变化（也就是说，对于任意的位移和扩张）下是不变的。

2.3 算法扩展

实际上，利用判别准则，可以直接将OTSU算法推广至 多阈值 的情形。例如：
在三阈值的情形下，可以选择两个阈值 0≤k1k2≤10 \le k_1 k_2 \le 10≤k1k2≤1 将原始灰度图像分化为三类。此时标准度量 η\etaη 存在两个参数 k1k_1k1 与 k2k_2k2，最佳阈值 k1∗,k2∗k_1^{*},k_{2}^{*}k1∗,k2∗ 可通过最大化 η\etaη：
(k1∗,k2∗)=arg⁡max⁡1≤k1≤k2≤L−1σB2(k1,k2)(24)(k_{1}^{*},k_{2}^{*})=\arg \underset{1 \le k_{1} \le k_{2} \le L-1}{\max}{\sigma_{B}^{2}(k_{1},k_{2})} \tag{24} (k1∗,k2∗)=arg1≤k1≤k2≤L−1maxσB2(k1,k2)(24)

进行求解。

参考文献

[1] Ostu N , Nobuyuki O , Otsu N . A thresholding selection method from gray level histogram. 1979.