算法001:大津法OTSU学习记录

OTSU算法

一、大津法主要的工作是什么？

大多数时候，我们需要获取到一幅图像中的特定目标。如果可以根据像素值将图像进行合理的分割，例如全局阈值分割那样，找到一个阈值TTT,大于阈值TTT的赋予一个像素值,小于阈值TTT的赋予一个像素值,那么就可以很容易的将目标与背景分离开来,如下公式所述:

g(x,y)={1,f(x,y)>T0,f(x,y)≤Tg(x,y) = \begin{cases} 1, f(x,y) > T\\ 0, f(x,y) \leq T\\ \end{cases} g(x,y)={1,f(x,y)>T0,f(x,y)≤T

如上所示:当TTT是一个适用于整幅图像的常数时, 该处理方式为全局阈值分割

在此基础上,出现了用OTSU方法实现的最佳全局阈值处理,通过该算法,可以找到最佳的全局阈值
总结起来: OTSU大津法是一个可以找到一幅图像的最佳全局阈值TTT的算法

二、使用算法的时候我们想要达到什么目的？如何才能说明求得阈值就是最佳阈值呢？

我们知道，经过阈值处理之后的图像其实是经过分类的，以某个阈值作为分类的**“分界点”**，理论上来讲，如果这个阈值足够好，那么像素被分成两组或者多组的过程中所引入的平均误差最小；从另一个角度来看，这两个组的平均灰度应该相差越大越好。
OTSU算法引入了类间方差的概念，然后当类间方差取到最大的时候，两个组的平均灰度相差最大，此时的阈值就是最佳阈值
总结起来：
目的就是可以是分成的两个类之间的差距尽可能的大，然后这个差距具体表现在平均灰度

最佳阈值就是可以让更好的让我们达成上述目的的阈值

三、推导过程

1.假设一幅大小为M×NM\times NM×N像素的图像,共有LLL个不同的灰度级,因此其灰度区间可以表示为[0,⋯,L−1][0,\cdots,L-1][0,⋯,L−1]。nin_ini表示灰度级为iii的像素的个数,那么就有:
MN=n0+n1+n2+⋯+nL−1MN = n_0+n_1+n_2+\cdots+n_{L-1} MN=n0+n1+n2+⋯+nL−1

像素为iii出现的概率约等于频率pi=niMNp_i =\frac {n_i}{MN}pi=MNni,于是可得:

∑i=0L−1pi=1,pi≥0\sum_{i=0}^{L-1} {p_i} = 1, p_i\geq0 i=0∑L−1pi=1,pi≥0
2. 如果先假定一个阈值T(k)=k,0<k<L−1T(k) = k,0<k<L-1T(k)=k,0<k<L−1,他会将图像阈值分割为C1C_1C1和C2C_2C2两类,其中C1C_1C1中的像素位于区间[0,k][0,k][0,k],C2C_2C2中的像素位于[k+1,L−1][k+1,L-1][k+1,L−1].

那么像素被分到类C1C_1C1中的概率为:

P1(k)=p0+p1+p2+⋯+pk=∑i=0kpiP_1(k) = p_0 + p_1 + p_2+\cdots +p_k =\sum_{i = 0}^{k} {p_i} P1(k)=p0+p1+p2+⋯+pk=i=0∑kpi

同理,像素被分到类C2C_2C2中的概率为:

P2(k)=pk+1+pk+2+pk+3+⋯+pL−1=∑i=k+1L−1piP_2(k) = p_{k+1} + p_{k+2} + p_{k+3} +\cdots +p_{L-1} =\sum_{i = k+1}^{L-1} {p_i} P2(k)=pk+1+pk+2+pk+3+⋯+pL−1=i=k+1∑L−1pi

很显然:

P1(k)+P2(k)=p0+p1+p2+⋯+pk+pk+1+pk+2+pk+3+⋯+pL−1=∑i=0kpi+∑i=k+1L−1pi=∑i=0L−1pi=1P_1(k)+P_2(k) =p_0 + p_1 + p_2+\cdots +p_k+p_{k+1} + p_{k+2} + p_{k+3} +\cdots +p_{L-1} =\sum_{i = 0}^{k} {p_i}+\sum_{i = k+1}^{L-1} {p_i}=\sum_{i= 0}^{L-1}p_i = 1 P1(k)+P2(k)=p0+p1+p2+⋯+pk+pk+1+pk+2+pk+3+⋯+pL−1=i=0∑kpi+i=k+1∑L−1pi=i=0∑L−1pi=1

简言之:

P1(k)+P2(k)=1P_1(k) + P_2(k) = 1 P1(k)+P2(k)=1

3. 在继续下面的推导之前,先引入下面的公式和概念（下述的是定义，当然也比较好理解）

某幅图像的平均灰度mmm为
m=∑i=0L−1riprim = \sum_{i = 0}^{L-1}r_i p_{r_i} m=i=0∑L−1ripri
如果这里的知识遗忘的话,可以回忆一下考研概论中的期望简述如下:
E(x)={∑xp(x),离散情况∫xf(x)dx,连续情况E(x) = \begin{cases} \sum_{}^{} {xp(x)},\,\,\,\,\,\,\,\,\,\,\text{离散情况}\\ \int{xf(x)}dx,\,\,\,\,\,\,\text{连续情况}\\ \end{cases} E(x)={∑xp(x),离散情况∫xf(x)dx,连续情况
某幅图像的灰度方差σ2\sigma^2σ2(二阶矩),也可以记为μ2(r)\mu_2(r)μ2(r)
μ2(r)=∑i=0L−1(ri−m)2pi\mu_2(r) = \sum_{i= 0}^{L-1}{(r_i-m)^2p_i} μ2(r)=i=0∑L−1(ri−m)2pi
其中rrr是像素值。

4.分到类C1C_1C1中像素的平均灰度值为:
m1(k)=∑i=0kip(i∣C1)m_1(k) = \sum_{i = 0}^{k}{ip(i|C_1)} m1(k)=i=0∑kip(i∣C1)
由贝叶斯公式:
P(A∣B)=P(B∣A)P(A)P(B)P(A|B) = \frac {P(B|A)P(A)}{P(B)} P(A∣B)=P(B)P(B∣A)P(A)

得到:
m1(k)=∑i=0kiP(i∣C1)=∑i=0kiP(C1∣i)P(i)P(C1)m_1(k) = \sum_{i = 0}^{k}{iP(i|C_1)} =\sum_{i=0}^{k}{i\frac {P(C_1|i)P(i)}{P(C_1)}} m1(k)=i=0∑kiP(i∣C1)=i=0∑kiP(C1)P(C1∣i)P(i)
又因为P(C1∣i)=1P(C_1|i) = 1P(C1∣i)=1(因为若像素其值为iii,则必定属于C1C_1C1),且P(C1)P(C_1)P(C1)不就是上文所述的P1(k)P_1(k)P1(k)嘛,因此可得下式:
m1(k)=∑i=0kiP(i)P(C1)=1P1(k)∑i=0kiPim_1(k) = \sum_{i = 0}^{k}{i\frac{P(i)}{P(C_1)}} = \frac {1}{P_1(k)}\sum_{i=0}^{k}{iP_i} m1(k)=i=0∑kiP(C1)P(i)=P1(k)1i=0∑kiPi
同理可得分到类C2C_2C2中像素的平均灰度值为:
m2(k)=1P2(k)∑i=k+1L−1iPim_2(k) = \frac{1}{P_2(k)}\sum_{i = k+1}^{L-1}{iP_i} m2(k)=P2(k)1i=k+1∑L−1iPi

5.对于kkk个灰度等级进行累加求均值时,表示为m(k)m(k)m(k):
m(k)=∑i=0kipim(k) = \sum_{i=0}^{k}{ip_i} m(k)=i=0∑kipi
于是—>全局的(也就是整个图像的)灰度均值mGm_GmG为:
mG＝∑i=0L−1ipim_G　＝\sum_{i = 0}^{L-1}{ip_i} mG　＝i=0∑L−1ipi
为什么呢?因为整幅图像的灰度级是[0,L−1][0,L-1][0,L−1]

6.经过上述的公式我们可以得到一下结论:(这里暂时忽略kkk)
mG=P1m1+P2m2m_G = P_1m_1 + P_2m_2 mG=P1m1+P2m2
mG−m=∑i=k+1L−1ipi=m2P2m_G-m = \sum_{i = k+1}^{L-1} {ip_i} = m_2P_2\\ mG−m=i=k+1∑L−1ipi=m2P2
P1+P2=1P_1 + P_2 = 1\\ P1+P2=1

7.我们知道OTSU算法又叫做最大类间方差法,所以自然少不了所谓的类间方差,其定义为:
σB2=P1(m1−mG)2+P2(m2−mG)2\sigma^2_B = P_1(m_1-m_G)^2+P_2(m_2-m_G)^2 σB2=P1(m1−mG)2+P2(m2−mG)2
上述的就是类间方差的定义,P1P_1P1、P2P_2P2可以理解为"权重",mGm_GmG是全局均值

当然还有一个重要的概念:全局方差σG2\sigma^2_GσG2(回顾前文中的灰度方差):

σG2=∑i=0L−1(i−mG)2pi\sigma^2_G = \sum_{i = 0}^{L-1}{(i-m_G)^2p_i} σG2=i=0∑L−1(i−mG)2pi

现在还是回到类间方差σB2\sigma^2_BσB2,经过推导我们发现σB2\sigma^2_BσB2还等价于下面两种形式:

σB2=P1P2(m1−m2)2=(mGP1−m)2P1(1−P1)\sigma^2_B = P_1P_2(m_1-m_2)^2=\frac{(m_GP_1-m)^2}{P_1(1-P_1)} σB2=P1P2(m1−m2)2=P1(1−P1)(mGP1−m)2

具体的推导过程见下图:

8.根据第7步的公式我们可以看出两个均值m1m_1m1和m2m_2m2彼此相差越大,σB2\sigma^2_BσB2就越大,还记得最开始的时候我们说过目的:就是可以是分成的两个类之间的差距尽可能的大，然后这个差距具体表现在平均灰度,在这里表现为均值m1m_1m1和m2m_2m2彼此之差。

也就是说，如果我们的阈值kkk选的最优，那么m1m_1m1和m2m_2m2的差距应该最大。另外因为不管阈值怎么改变全局的平均灰度σG2\sigma^2_GσG2都是不受影响的,也即是个常数。

此时再引入kkk,用类间方差σB2\sigma^2_BσB2和全局方差σG2\sigma^2_GσG2的商η\etaη来表征:
η(k)=σB2(k)σG2(k)\eta(k) = \frac{\sigma^2_B(k)}{\sigma^2_G(k)} η(k)=σG2(k)σB2(k)

σB2(k)=[mG(k)P1(k)−m(k)]2P1(k)[1−P1(k)]\sigma^2_B(k) = \frac{[m_G(k)P_1(k)-m(k)]^2}{P_1(k)[1-P_1(k)]} σB2(k)=P1(k)[1−P1(k)][mG(k)P1(k)−m(k)]2

如果k∗k^*k∗是最佳阈值,那么σB2(k∗)\sigma^2_B(k^*)σB2(k∗)就是最大化的σB2(k)\sigma^2_B(k)σB2(k).

9.至此,只要我们找出可以使得σB2(k)\sigma^2_B(k)σB2(k)最大化的k∗k^*k∗就算是找到了最优阈值

一般的寻找方式是:取kkk所有可能取到的整数值-----位于[0,L−1][0,L-1][0,L−1]中,挨个选取,求出使得σB2(k)\sigma^2_B(k)σB2(k)最大的kkk值即可

如果存在多个kkk值使得σB2(k)\sigma^2_B(k)σB2(k)最大,那么求kkk的均值即可.

之后将kkk作为全局阈值,进行图像分割:
g(x,y)={1,f(x,y)>k0,f(x,y)≤kg(x,y) = \begin{cases} 1, f(x,y) > k\\ 0, f(x,y) \leq k\\ \end{cases} g(x,y)={1,f(x,y)>k0,f(x,y)≤k