阈值分割：最大类间方差法
一、简介
最大类间方差法，又称为大津阈值法，或OTSU算法。是由日本学者大津在
1979
1979
1
9
7
9
年提出的一种非参数的、无监督的自动选择阈值的图像分割方法。
二、算法描述
2.1 公式推导
对于给定的一幅具有
L
L
L
个灰度级（
[
0
,
1
,
2
,
⋯
,
L
−
1
]
[0,1,2,\cdots,L-1]
[
0
,
1
,
2
,
⋯
,
L
−
1
]
）的灰度图像，有以下描述：
-
使用
ni
n_{i}
n
i
表示处于灰度级
ii
i
的像素块的数目； -
使用
N=
n
0
+
n
1
+
⋯
+
n
L
−
1
N = n_0 + n_1 + \cdots + n_{L-1}
N
=
n
0
+
n
1
+
⋯
+
n
L
−
1
表示该图像所有像素块的数目总和。
对该图像的灰度级直方图进行
标准化处理
可以得到一个概率分布：
p
i
=
n
i
N
，
p
i
≥
0
,
∑
i
=
0
L
−
1
p
i
=
1
(1)
p_i = \frac{n_i}{N}，\quad p_i \ge 0,\quad \sum_{i=0}^{L-1}{p_i} =1 \tag{1}
p
i
=
N
n
i
，
p
i
≥
0
,
i
=
0
∑
L
−
1
p
i
=
1
(
1
)
以灰度级
k
k
k
为
阈值
，可将该图像中的像素块二分为两大类
C
0
C_{0}
C
0
和
C
1
C_{1}
C
1
（背景和前景，反之亦然），其中：
-
C0
C_0
C
0
表示灰度级在范围
[0
,
1
,
⋯
,
k
−
1
]
[0,1,\cdots,k-1]
[
0
,
1
,
⋯
,
k
−
1
]
的像素集合； -
C1
C_1
C
1
表示灰度级在范围
[k
,
k
+
1
,
⋯
,
L
−
1
]
[k,k+1,\cdots,L-1]
[
k
,
k
+
1
,
⋯
,
L
−
1
]
的像素集合。
记
w
0
w_0
w
0
与
w
1
w_1
w
1
分别为类
C
0
C_0
C
0
和类
C
1
C_1
C
1
发生的概率，则显然有：
∑
i
=
1
L
−
1
p
i
=
w
0
+
w
1
=
1
(2)
\sum_{i = 1}^{L-1}{p_i} = w_0 + w_1 = 1 \tag{2}
i
=
1
∑
L
−
1
p
i
=
w
0
+
w
1
=
1
(
2
)
其中，
w
0
w_0
w
0
与
w
1
w_1
w
1
可表示为
零阶矩
的形式：
w
0
=
P
r
(
C
0
)
=
∑
i
=
0
k
−
1
p
i
=
w
(
k
)
(3)
w_0 = Pr(C_0) = \sum_{i = 0}^{k-1}{p_{i}} = w(k) \tag{3}
w
0
=
P
r
(
C
0
)
=
i
=
0
∑
k
−
1
p
i
=
w
(
k
)
(
3
)
w
1
=
P
r
(
C
1
)
=
∑
i
=
k
L
−
1
p
i
=
1
−
w
(
k
)
(4)
w_1 = Pr(C_1) = \sum_{i = k}^{L-1}{p_{i}} = 1 – w(k) \tag{4}
w
1
=
P
r
(
C
1
)
=
i
=
k
∑
L
−
1
p
i
=
1
−
w
(
k
)
(
4
)
类
C
0
C_0
C
0
与 类
C
1
C_1
C
1
发生的
期望
（或均值）可表示为
一阶矩
的形式：
μ
0
=
∑
i
=
0
k
−
1
i
P
r
(
i
∣
C
0
)
=
∑
i
=
0
k
−
1
i
p
i
w
0
=
1
w
0
∑
i
=
0
k
−
1
i
p
i
=
μ
(
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
−
1
i
P
r
(
i
∣
C
0
)
=
i
=
0
∑
k
−
1
i
w
0
p
i
=
w
0
1
i
=
0
∑
k
−
1
i
p
i
=
w
(
k
)
μ
(
k
)
(
5
)
μ
1
=
∑
i
=
k
L
−
1
i
P
r
(
i
∣
C
1
)
=
∑
i
=
k
L
−
1
i
p
i
w
1
=
1
w
1
∑
i
=
k
L
−
1
i
p
i
=
μ
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
−
1
i
P
r
(
i
∣
C
1
)
=
i
=
k
∑
L
−
1
i
w
1
p
i
=
w
1
1
i
=
k
∑
L
−
1
i
p
i
=
1
−
w
(
k
)
μ
T
−
μ
(
k
)
(
5
)
其中：
μ
(
k
)
=
∑
i
=
0
k
−
1
i
p
i
,
μ
T
=
∑
i
=
0
L
−
1
i
p
i
(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
−
1
i
p
i
,
μ
T
=
i
=
0
∑
L
−
1
i
p
i
(
6
)
显然，有：
μ
0
w
0
+
μ
1
w
1
=
μ
T
(7)
\mu_0 w_0 + \mu_1 w_1 =\mu_{T} \tag{7}
μ
0
w
0
+
μ
1
w
1
=
μ
T
(
7
)
类
C
0
C_0
C
0
与 类
C
1
C_1
C
1
发生的
方差
可以表示为
二阶矩
的形式：
σ
0
2
=
∑
i
=
0
k
−
1
(
i
−
μ
0
)
2
P
r
(
i
∣
C
0
)
=
∑
i
=
0
k
−
1
(
i
−
μ
0
)
2
p
i
w
0
(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}
σ
0
2
=
i
=
0
∑
k
−
1
(
i
−
μ
0
)
2
P
r
(
i
∣
C
0
)
=
i
=
0
∑
k
−
1
(
i
−
μ
0
)
2
w
0
p
i
(
8
)
σ
1
2
=
∑
i
=
k
L
−
1
(
i
−
μ
1
)
2
P
r
(
i
∣
C
1
)
=
∑
i
=
k
L
−
1
(
i
−
μ
1
)
2
p
i
w
1
(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}
σ
1
2
=
i
=
k
∑
L
−
1
(
i
−
μ
1
)
2
P
r
(
i
∣
C
1
)
=
i
=
k
∑
L
−
1
(
i
−
μ
1
)
2
w
1
p
i
(
9
)
为了评估所选阈值
k
k
k
的
优良
（goodness），引入以下三种 判别标准度量：
λ
=
σ
B
2
σ
W
2
,
κ
=
σ
T
2
σ
W
2
,
η
=
σ
B
2
]
σ
T
2
(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}
λ
=
σ
W
2
σ
B
2
,
κ
=
σ
W
2
σ
T
2
,
η
=
σ
T
2
σ
B
2
]
(
1
0
)
其中：
（1）
类内方差
（within-class variance，简记为
σ
w
i
t
h
2
\sigma_{with}^{2}
σ
w
i
t
h
2
或
σ
W
2
\sigma_{W}^{2}
σ
W
2
）满足：
σ
W
2
=
w
0
σ
0
2
+
w
1
σ
1
2
(11)
\sigma_{W}^{2} = w_0 \sigma_{0}^{2} + w_1 \sigma_{1}^{2} \tag{11}
σ
W
2
=
w
0
σ
0
2
+
w
1
σ
1
2
(
1
1
)
（2）
类间方差
（between-class variance，简记为
σ
B
e
t
w
e
e
n
2
\sigma_{Between}^{2}
σ
B
e
t
w
e
e
n
2
或
σ
B
2
\sigma_{B}^{2}
σ
B
2
）满足：
σ
B
2
=
w
0
(
μ
0
−
μ
T
)
2
+
w
1
(
μ
1
−
μ
T
)
2
=
w
0
w
1
(
μ
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}
σ
B
2
=
w
0
(
μ
0
−
μ
T
)
2
+
w
1
(
μ
1
−
μ
T
)
2
=
w
0
w
1
(
μ
1
−
μ
0
)
2
(
1
2
)
（3）
全局方差
（total variance，简记为
σ
T
o
t
a
l
2
\sigma_{Total}^{2}
σ
T
o
t
a
l
2
或
σ
T
2
\sigma_{T}^{2}
σ
T
2
）满足：
σ
T
2
=
∑
i
=
0
L
−
1
(
i
−
μ
T
)
2
p
i
(13)
\sigma_{T}^{2} = \sum_{i = 0}^{L-1}{(i – \mu_{T})^{2}p_{i}} \tag{13}
σ
T
2
=
i
=
0
∑
L
−
1
(
i
−
μ
T
)
2
p
i
(
1
3
)
合适的阈值会将图像分割为两类。反过来就是说，能在灰度水平上实现最佳分离的阈值将是最合适的阈值。
因此，利用引入的判别标准度量，可以将问题转化为一个优化问题：寻找一个合适的阈值
k
k
k
最大化某一个判别标准度量函数（
λ
\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
λ
(
1
4
)
并且以下基本关系式始终成立：
σ
W
2
+
σ
B
2
=
σ
T
2
(15)
\sigma_{W}^{2} + \sigma_{B}^{2} = \sigma_{T}^{2} \tag{15}
σ
W
2
+
σ
B
2
=
σ
T
2
(
1
5
)
也就是说，当三个判别标准度量中的任意一个达到最大时，另外两个都会达到最大值。其中：
-
σW
2
\sigma_{W}^{2}
σ
W
2
和
σB
2
\sigma_{B}^{2}
σ
B
2
均是阈值
kk
k
的函数； -
σT
2
\sigma_{T}^{2}
σ
T
2
与阈值
kk
k
无关； -
σW
\sigma_{W}
σ
W
需要计算二阶矩，而
σB
2
\sigma_{B}^{2}
σ
B
2
仅仅需要计算一阶矩。
所以，
η
\eta
η
是关于阈值
k
k
k
的最简单的判别标准度量。
因此，使用
η
\eta
η
作为评估阈值
k
k
k
优良 的判别标准。
现在，求解最大化判别标准
η
\eta
η
时的最佳阈值
k
∗
k^{*}
k
∗
。 分析
η
=
σ
B
2
σ
T
2
\eta = \frac{\sigma_{B}^{2}}{\sigma_{T}^{2}}
η
=
σ
T
2
σ
B
2
可知，最大化
η
\eta
η
即最大化
σ
B
2
\sigma_{B}^{2}
σ
B
2
。即：
η
(
k
)
=
σ
B
2
(
k
)
σ
T
2
(16)
\eta(k) = \frac{\sigma_{B}^{2}(k)}{\sigma_{T}^{2}} \tag{16}
η
(
k
)
=
σ
T
2
σ
B
2
(
k
)
(
1
6
)
并且有：
σ
B
2
(
k
)
=
[
μ
T
w
(
k
)
−
μ
(
k
)
]
2
w
(
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}
σ
B
2
(
k
)
=
w
(
k
)
[
1
−
w
(
k
)
]
[
μ
T
w
(
k
)
−
μ
(
k
)
]
2
(
1
7
)
最优阈值
k
∗
k^{*}
k
∗
可以表示为：
k
∗
=
arg
max
1
≤
k
≤
L
−
1
σ
B
2
(
k
)
(18)
k^{*} = \arg \underset{1 \le k \le L-1}{\max}{\sigma_{B}^{2}{(k)}} \tag{18}
k
∗
=
ar
g
1
≤
k
≤
L
−
1
max
σ
B
2
(
k
)
(
1
8
)
其中，阈值
k
k
k
的搜索范围可以表示为：
S
∗
=
{
k
∣
w
0
w
1
=
w
(
k
)
[
1
−
w
(
k
)
]
>
0
,
o
r
0
<
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
∣
w
0
w
1
=
w
(
k
)
[
1
−
w
(
k
)
]
>
0
,
o
r
0
<
w
(
k
)
<
1
}
(
1
9
)
该范围称为灰度直方图的有效范围。
从公式（12）中
σ
B
2
=
w
0
(
μ
0
−
μ
T
)
2
+
w
1
(
μ
1
−
μ
T
)
2
=
w
0
w
1
(
μ
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}
σ
B
2
=
w
0
(
μ
0
−
μ
T
)
2
+
w
1
(
μ
1
−
μ
T
)
2
=
w
0
w
1
(
μ
1
−
μ
0
)
2
可以看出：
-
当选定的阈值
k∈
S
−
S
∗
=
{
k
∣
w
(
k
)
=
0
o
r
1
}
k \in S – S^{*} = \{k ~|~w(k) = 0~or~1\}
k
∈
S
−
S
∗
=
{
k
∣
w
(
k
)
=
0
o
r
1
}
时，判别标准度量
η\eta
η
取最小值
00
0
。 -
当选定的阈值
k∈
S
∗
k \in S^{*}
k
∈
S
∗
时，判别标准度量
η\eta
η
取一个正的且有界的值。
因此，显而易见，判别标准的最大值始终存在。
2.2 算法分析
对于选定的阈值
k
∗
k^{*}
k
∗
：
w
0
∗
=
P
r
(
C
0
∗
)
=
∑
i
=
0
k
∗
−
1
p
i
=
w
(
k
∗
)
(20)
w_{0}^{*} = Pr(C_{0}^{*}) = \sum_{i = 0}^{k^{*} – 1}{p_i} = w(k^{*}) \tag{20}
w
0
∗
=
P
r
(
C
0
∗
)
=
i
=
0
∑
k
∗
−
1
p
i
=
w
(
k
∗
)
(
2
0
)
w
1
∗
=
P
r
(
C
1
∗
)
=
∑
i
=
k
∗
L
−
1
p
i
=
1
−
w
(
k
∗
)
(21)
w_{1}^{*} = Pr(C_1^{*}) = \sum_{i = k^{*}}^{L-1}{p_{i}} = 1 – w(k^{*}) \tag{21}
w
1
∗
=
P
r
(
C
1
∗
)
=
i
=
k
∗
∑
L
−
1
p
i
=
1
−
w
(
k
∗
)
(
2
1
)
分别表示了灰度图像按照阈值
k
∗
k^{*}
k
∗
所划分的两类的发生概率。
类
C
0
∗
C_{0}^{*}
C
0
∗
与类
C
1
∗
C_{1}^{*}
C
1
∗
发生的 期望 分别为：
μ
0
∗
=
∑
i
=
0
k
∗
−
1
i
P
r
(
i
∣
C
0
∗
)
=
∑
i
=
0
k
∗
−
1
i
p
i
w
0
∗
=
μ
(
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
∗
−
1
i
P
r
(
i
∣
C
0
∗
)
=
i
=
0
∑
k
∗
−
1
i
w
0
∗
p
i
=
w
(
k
∗
)
μ
(
k
∗
)
(
2
2
)
μ
1
∗
=
∑
i
=
k
∗
L
−
1
i
P
r
(
i
∣
C
1
∗
)
=
∑
i
=
k
∗
L
−
1
i
p
i
w
1
∗
=
μ
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
−
1
i
P
r
(
i
∣
C
1
∗
)
=
i
=
k
∗
∑
L
−
1
i
w
1
∗
p
i
=
1
−
w
(
k
∗
)
μ
T
−
μ
(
k
∗
)
(
2
3
)
将判别标准
η
\eta
η
的最大值
η
(
k
∗
)
\eta(k^{*})
η
(
k
∗
)
简记为
η
∗
\eta^{*}
η
∗
，可以用作评估灰度图像中类的可分性的标准。这是一个重要的度量，它在灰度尺度的放射变化（也就是说，对于任意的位移和扩张）下是不变的。
2.3 算法扩展
实际上，利用判别准则，可以直接将OTSU算法推广至
多阈值
的情形。例如：
在三阈值的情形下，可以选择两个阈值
0
≤
k
1
k
2
≤
1
0 \le k_1 k_2 \le 1
0
≤
k
1
k
2
≤
1
将原始灰度图像分化为三类。此时标准度量
η
\eta
η
存在两个参数
k
1
k_1
k
1
与
k
2
k_2
k
2
，最佳阈值
k
1
∗
,
k
2
∗
k_1^{*},k_{2}^{*}
k
1
∗
,
k
2
∗
可通过最大化
η
\eta
η
：
(
k
1
∗
,
k
2
∗
)
=
arg
max
1
≤
k
1
≤
k
2
≤
L
−
1
σ
B
2
(
k
1
,
k
2
)
(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}
(
k
1
∗
,
k
2
∗
)
=
ar
g
1
≤
k
1
≤
k
2
≤
L
−
1
max
σ
B
2
(
k
1
,
k
2
)
(
2
4
)
进行求解。
参考文献
[1] Ostu N , Nobuyuki O , Otsu N . A thresholding selection method from gray level histogram. 1979.