阈值分割:最大类间方差法
一、简介
最大类间方差法,又称为大津阈值法,或OTSU算法。是由日本学者大津在 1979 1979 1979年提出的一种非参数的、无监督的自动选择阈值的图像分割方法。
二、算法描述
2.1 公式推导
对于给定的一幅具有 L L L 个灰度级( [ 0 , 1 , 2 , ⋯ , L − 1 ] [0,1,2,\cdots,L-1] [0,1,2,⋯,L−1])的灰度图像,有以下描述:
- 使用 n i n_{i} ni 表示处于灰度级 i i i 的像素块的数目;
- 使用 N = n 0 + n 1 + ⋯ + n L − 1 N = n_0 + n_1 + \cdots + n_{L-1} N=n0+n1+⋯+nL−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}
pi=Nni,pi≥0,i=0∑L−1pi=1(1)
以灰度级 k k k 为 阈值,可将该图像中的像素块二分为两大类 C 0 C_{0} C0 和 C 1 C_{1} C1 (背景和前景,反之亦然),其中:
- C 0 C_0 C0 表示灰度级在范围 [ 0 , 1 , ⋯ , k − 1 ] [0,1,\cdots,k-1] [0,1,⋯,k−1] 的像素集合;
- C 1 C_1 C1 表示灰度级在范围 [ k , k + 1 , ⋯ , L − 1 ] [k,k+1,\cdots,L-1] [k,k+1,⋯,L−1] 的像素集合。
记
w
0
w_0
w0 与
w
1
w_1
w1 分别为类
C
0
C_0
C0 和类
C
1
C_1
C1 发生的概率,则显然有:
∑
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−1pi=w0+w1=1(2)
其中,
w
0
w_0
w0 与
w
1
w_1
w1 可表示为 零阶矩 的形式:
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}
w0=Pr(C0)=i=0∑k−1pi=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} w1=Pr(C1)=i=k∑L−1pi=1−w(k)(4)
类
C
0
C_0
C0 与 类
C
1
C_1
C1 发生的 期望 (或均值)可表示为 一阶矩 的形式:
μ
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−1iPr(i∣C0)=i=0∑k−1iw0pi=w01i=0∑k−1ipi=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−1iPr(i∣C1)=i=k∑L−1iw1pi=w11i=k∑L−1ipi=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−1ipi,μT=i=0∑L−1ipi(6)
显然,有:
μ
0
w
0
+
μ
1
w
1
=
μ
T
(7)
\mu_0 w_0 + \mu_1 w_1 =\mu_{T} \tag{7}
μ0w0+μ1w1=μT(7)
类
C
0
C_0
C0 与 类
C
1
C_1
C1 发生的 方差 可以表示为 二阶矩 的形式:
σ
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}
σ02=i=0∑k−1(i−μ0)2Pr(i∣C0)=i=0∑k−1(i−μ0)2w0pi(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} σ12=i=k∑L−1(i−μ1)2Pr(i∣C1)=i=k∑L−1(i−μ1)2w1pi(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}
λ=σW2σB2,κ=σW2σT2,η=σT2σB2](10)
其中:
(1)类内方差(within-class variance,简记为
σ
w
i
t
h
2
\sigma_{with}^{2}
σwith2 或
σ
W
2
\sigma_{W}^{2}
σW2)满足:
σ
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}
σW2=w0σ02+w1σ12(11)
(2)类间方差(between-class variance,简记为
σ
B
e
t
w
e
e
n
2
\sigma_{Between}^{2}
σBetween2 或
σ
B
2
\sigma_{B}^{2}
σB2)满足:
σ
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}
σB2=w0(μ0−μT)2+w1(μ1−μT)2=w0w1(μ1−μ0)2(12)
(3)全局方差(total variance,简记为
σ
T
o
t
a
l
2
\sigma_{Total}^{2}
σTotal2 或
σ
T
2
\sigma_{T}^{2}
σT2)满足:
σ
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}
σT2=i=0∑L−1(i−μT)2pi(13)
合适的阈值会将图像分割为两类。反过来就是说,能在灰度水平上实现最佳分离的阈值将是最合适的阈值。
因此,利用引入的判别标准度量,可以将问题转化为一个优化问题:寻找一个合适的阈值 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λ(14)
并且以下基本关系式始终成立:
σ
W
2
+
σ
B
2
=
σ
T
2
(15)
\sigma_{W}^{2} + \sigma_{B}^{2} = \sigma_{T}^{2} \tag{15}
σW2+σB2=σT2(15)
也就是说,当三个判别标准度量中的任意一个达到最大时,另外两个都会达到最大值。其中:
- σ W 2 \sigma_{W}^{2} σW2 和 σ B 2 \sigma_{B}^{2} σB2 均是阈值 k k k 的函数;
- σ T 2 \sigma_{T}^{2} σT2 与阈值 k k k 无关;
- σ W \sigma_{W} σW 需要计算二阶矩,而 σ B 2 \sigma_{B}^{2} σB2 仅仅需要计算一阶矩。
所以, η \eta η 是关于阈值 k k k 的最简单的判别标准度量。
因此,使用 η \eta η 作为评估阈值 k k k 优良 的判别标准。
现在,求解最大化判别标准
η
\eta
η 时的最佳阈值
k
∗
k^{*}
k∗。 分析
η
=
σ
B
2
σ
T
2
\eta = \frac{\sigma_{B}^{2}}{\sigma_{T}^{2}}
η=σT2σB2 可知,最大化
η
\eta
η 即最大化
σ
B
2
\sigma_{B}^{2}
σB2。即:
η
(
k
)
=
σ
B
2
(
k
)
σ
T
2
(16)
\eta(k) = \frac{\sigma_{B}^{2}(k)}{\sigma_{T}^{2}} \tag{16}
η(k)=σT2σB2(k)(16)
并且有:
σ
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}
σB2(k)=w(k)[1−w(k)][μTw(k)−μ(k)]2(17)
最优阈值
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∗=arg1≤k≤L−1maxσB2(k)(18)
其中,阈值
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 ∣ w0w1=w(k)[1−w(k)]>0, or 0<w(k)<1}(19)
该范围称为灰度直方图的有效范围。
从公式(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}
σB2=w0(μ0−μT)2+w1(μ1−μT)2=w0w1(μ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 or 1} 时,判别标准度量 η \eta η 取最小值 0 0 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}
w0∗=Pr(C0∗)=i=0∑k∗−1pi=w(k∗)(20)
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} w1∗=Pr(C1∗)=i=k∗∑L−1pi=1−w(k∗)(21)
分别表示了灰度图像按照阈值 k ∗ k^{*} k∗ 所划分的两类的发生概率。
类
C
0
∗
C_{0}^{*}
C0∗ 与类
C
1
∗
C_{1}^{*}
C1∗ 发生的 期望 分别为:
μ
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∗−1iPr(i∣C0∗)=i=0∑k∗−1iw0∗pi=w(k∗)μ(k∗)(22)
μ
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−1iPr(i∣C1∗)=i=k∗∑L−1iw1∗pi=1−w(k∗)μT−μ(k∗)(23)
将判别标准 η \eta η 的最大值 η ( k ∗ ) \eta(k^{*}) η(k∗) 简记为 η ∗ \eta^{*} η∗,可以用作评估灰度图像中类的可分性的标准。这是一个重要的度量,它在灰度尺度的放射变化(也就是说,对于任意的位移和扩张)下是不变的。
2.3 算法扩展
实际上,利用判别准则,可以直接将OTSU算法推广至 多阈值 的情形。例如:
在三阈值的情形下,可以选择两个阈值
0
≤
k
1
k
2
≤
1
0 \le k_1 k_2 \le 1
0≤k1k2≤1 将原始灰度图像分化为三类。此时标准度量
η
\eta
η 存在两个参数
k
1
k_1
k1 与
k
2
k_2
k2,最佳阈值
k
1
∗
,
k
2
∗
k_1^{*},k_{2}^{*}
k1∗,k2∗ 可通过最大化
η
\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}
(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.