[注:下面公式所涉及的是带权重的均值与方差,一开始我并不明白为什么要这样做,还去查了关于带权重与不带权重计算均值或方差的区别,后面发现,应该是因为该算法的计算是基于概率分布,而概率分布意味着这是一个可以根据概率,产生无数样本的population,而非sample,因此对概率分布求均值就相当于求期望,方差也同理,需要带上概率]
原理
一幅图片的所有像素有
L
L
L个灰度级,
[
1
,
2
,
⋯
,
L
]
[1,2,\cdots,L]
[1,2,⋯,L]
处于灰度级
i
i
i的像素数量表示为
n
i
n_i
ni,因此像素的总数为
N
=
n
1
+
n
2
+
⋯
+
n
L
N=n_1+n_2+\cdots+n_L
N=n1+n2+⋯+nL。后续为了简化讨论,将灰度直方图归一化,因而就能被视为一个概率分布:
p
i
=
n
i
N
,
p
i
≥
,
∑
i
=
1
L
p
i
=
1
(
1
)
p_i=\frac{n_i}{N},\quad p_i\ge,\sum_{i=1}^{L}p_i=1 \qquad(1)
pi=Nni,pi≥,i=1∑Lpi=1(1)
现在,将所有像素以灰度级
k
k
k为阈值,划分为两类
C
0
C_0
C0和
C
1
C_1
C1,
C
0
C_0
C0代表处于灰度级
[
1
,
⋯
,
k
]
[1,\cdots,k]
[1,⋯,k]的像素,
C
1
C_1
C1代表处于灰度级
[
k
+
1
,
⋯
,
L
]
[k+1,\cdots,L]
[k+1,⋯,L]的像素。
某一类别出现的概率为:
ω
0
=
P
r
(
C
0
)
=
∑
i
=
1
k
p
i
=
ω
(
k
)
(
2
)
\omega_0=Pr(C_0)=\sum_{i=1}^{k}p_i=\omega(k) \qquad(2)
ω0=Pr(C0)=i=1∑kpi=ω(k)(2)
ω
1
=
P
r
(
C
1
)
=
∑
i
=
k
+
1
L
p
i
=
1
−
ω
(
k
)
(
3
)
\omega_1=Pr(C_1)=\sum_{i=k+1}^{L}p_i=1-\omega(k) \qquad(3)
ω1=Pr(C1)=i=k+1∑Lpi=1−ω(k)(3)
某一类别的平均灰度级为(注意是条件概率):
μ
0
=
∑
i
=
1
k
i
P
r
(
i
∣
C
0
)
=
∑
i
=
1
k
i
p
i
ω
0
=
μ
(
k
)
ω
(
k
)
(
4
)
\mu_0=\sum_{i=1}^{k}iPr(i|C_0)=\sum_{i=1}^{k}\frac{ip_i}{\omega_0}=\frac{\mu(k)}{\omega(k)} \qquad(4)
μ0=i=1∑kiPr(i∣C0)=i=1∑kω0ipi=ω(k)μ(k)(4)
其中
μ
(
k
)
=
∑
i
=
1
k
i
p
i
(
5
)
\mu(k)=\sum_{i=1}^{k}ip_i\qquad(5)
μ(k)=i=1∑kipi(5)
μ
(
k
)
\mu(k)
μ(k)可以看成灰度直方图到k的一阶矩,同理,公式(2)中的
ω
(
k
)
\omega(k)
ω(k)可以看成灰度直方图到k的零阶矩。
μ
1
=
∑
i
=
k
+
1
L
i
P
r
(
i
∣
C
1
)
=
∑
i
=
k
+
1
L
i
p
i
ω
1
=
μ
T
−
μ
(
k
)
1
−
ω
(
k
)
(
6
)
\mu_1=\sum_{i=k+1}^{L}iPr(i|C_1)=\sum_{i=k+1}^{L}\frac{ip_i}{\omega_1}=\frac{\mu_T-\mu(k)}{1-\omega(k)} \qquad(6)
μ1=i=k+1∑LiPr(i∣C1)=i=k+1∑Lω1ipi=1−ω(k)μT−μ(k)(6)
其中
μ
T
=
μ
(
L
)
=
∑
i
=
1
L
i
p
i
(
7
)
\mu_T=\mu(L)=\sum_{i=1}^{L}ip_i\qquad(7)
μT=μ(L)=i=1∑Lipi(7)
μ
T
\mu_T
μT是原图像的整体平均灰度级。根据上述公式,对于任意的k,我们有:
ω
0
μ
0
+
ω
1
μ
1
=
μ
T
,
ω
0
+
ω
1
=
1
(
8
)
\omega_0\mu_0+\omega_1\mu_1=\mu_T,\quad\omega_0+\omega_1=1\qquad(8)
ω0μ0+ω1μ1=μT,ω0+ω1=1(8)
(ps:这可以理解为全概率公式,也是公式(4),(6)用条件概率的原因,比如考虑人患艾滋病的概率
P
(
S
)
P(S)
P(S),现在你可以把人群分的细一些,比如(男性
M
M
M,女性
W
W
W),它们都有自己的发生概率。当考虑一个完全的分类结果时候,你就得到了全概率公式。在 这里你会得到:患病概率 = 男性患病概率 * 男性人口比例 + 女性患病概率 * 女性人口比例。用概率公式写出来就是
P
(
S
)
=
P
(
S
∣
M
)
P
(
M
)
+
P
(
S
∣
W
)
P
(
W
)
P(S)=P(S|M)P(M)+P(S|W)P(W)
P(S)=P(S∣M)P(M)+P(S∣W)P(W))
某一类的灰度级方差为:
σ
0
2
=
∑
i
=
1
k
(
i
−
μ
0
)
2
P
r
(
i
∣
C
0
)
=
∑
i
=
1
k
(
i
−
μ
0
)
2
p
i
ω
0
(
9
)
\sigma_0^{2}=\sum_{i=1}^{k}(i-\mu_0)^2Pr(i|C_0)=\sum_{i=1}^{k}(i-\mu_0)^2\frac{p_i}{\omega_0}\qquad(9)
σ02=i=1∑k(i−μ0)2Pr(i∣C0)=i=1∑k(i−μ0)2ω0pi(9)
σ
1
2
=
∑
i
=
k
+
1
L
(
i
−
μ
1
)
2
P
r
(
i
∣
C
1
)
=
∑
i
=
k
+
1
L
(
i
−
μ
1
)
2
p
i
ω
1
(
10
)
\sigma_1^{2}=\sum_{i=k+1}^{L}(i-\mu_1)^2Pr(i|C_1)=\sum_{i=k+1}^{L}(i-\mu_1)^2\frac{p_i}{\omega_1}\qquad(10)
σ12=i=k+1∑L(i−μ1)2Pr(i∣C1)=i=k+1∑L(i−μ1)2ω1pi(10)
接下来,为了衡量阈值k选取的好坏,我们首先引进线性判别式分析(Linear Discriminant Analysis, LDA)里的思想,希望选定的阈值能使不同的类之间相差越大越好,相同的类则相差越小越好(即最大化
λ
\lambda
λ),Otsu论文中的公式如下(但是在LDA里,我只看到了
λ
\lambda
λ,不知道
κ
\kappa
κ和
η
\eta
η是本文提出来的还是怎样):
λ
=
σ
B
2
σ
W
2
,
κ
=
σ
T
2
σ
W
2
,
η
=
σ
B
2
σ
T
2
(
11
)
\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}} \qquad(11)
λ=σW2σB2,κ=σW2σT2,η=σT2σB2(11)
其中
σ
W
2
\sigma_W^{2}
σW2、
σ
B
2
\sigma_B^{2}
σB2、
σ
T
2
\sigma_T^{2}
σT2分别代表类内(within-class)方差,类间(between-class)方差与所有灰度级的总(total levels)方差。计算方法如下:
σ
W
2
=
ω
0
σ
0
2
+
ω
1
σ
1
2
(
12
)
\sigma_W^{2}=\omega_0\sigma_0^{2}+\omega_1\sigma_1^{2}\qquad(12)
σW2=ω0σ02+ω1σ12(12)
σ
B
2
=
ω
0
(
μ
0
−
μ
T
)
2
+
ω
1
(
μ
1
−
μ
T
)
2
=
ω
0
ω
1
(
μ
1
−
μ
0
)
2
(
13
)
\sigma_B^{2}=\omega_0(\mu_0-\mu_T)^{2}+\omega_1(\mu_1-\mu_T)^{2}=\omega_0\omega_1(\mu_1-\mu_0)^2\qquad(13)
σB2=ω0(μ0−μT)2+ω1(μ1−μT)2=ω0ω1(μ1−μ0)2(13)
σ
T
2
=
∑
i
=
1
L
(
i
−
μ
T
)
2
p
i
(
14
)
\sigma_T^{2}=\sum_{i=1}^{L}(i-\mu_T)^2p_i\qquad(14)
σT2=i=1∑L(i−μT)2pi(14)
又因为
σ
W
2
+
σ
B
2
=
σ
T
2
\sigma_W^{2}+\sigma_B^{2}=\sigma_T^{2}
σW2+σB2=σT2,因此
κ
=
λ
+
1
\kappa=\lambda+1
κ=λ+1,
η
=
λ
λ
+
1
\eta=\frac{\lambda}{\lambda+1}
η=λ+1λ,都是关于
λ
\lambda
λ的函数,而
σ
W
2
\sigma_{W}^2
σW2涉及二阶统计量(方差),
σ
B
2
\sigma_{B}^2
σB2只涉及一阶统计量,因此,我们只需要选取其中最简单的
η
\eta
η求最大值即可。由于
σ
T
2
\sigma_T^2
σT2与k的选取无关,因此最后的优化目标等同于最大化
σ
B
2
\sigma_B^2
σB2,类间方差。
代码
def otsu(gray):
pixel_number = gray.shape[0] * gray.shape[1]
mean_weigth = 1.0/pixel_number
# 发现bins必须写到257,否则255这个值只能分到[254,255)区间
his, bins = np.histogram(gray, np.arange(0,257))
final_thresh = -1
final_value = -1
intensity_arr = np.arange(256)
for t in bins[1:-1]: # This goes from 1 to 254 uint8 range (Pretty sure wont be those values)
pcb = np.sum(his[:t])
pcf = np.sum(his[t:])
Wb = pcb * mean_weigth
Wf = pcf * mean_weigth
mub = np.sum(intensity_arr[:t]*his[:t]) / float(pcb)
muf = np.sum(intensity_arr[t:]*his[t:]) / float(pcf)
#print mub, muf
value = Wb * Wf * (mub - muf) ** 2
if value > final_value:
final_thresh = t
final_value = value
final_img = gray.copy()
print(final_thresh)
final_img[gray > final_thresh] = 255
final_img[gray < final_thresh] = 0
return final_img