Otsu算法原理与python实现

[注：下面公式所涉及的是带权重的均值与方差，一开始我并不明白为什么要这样做，还去查了关于带权重与不带权重计算均值或方差的区别，后面发现，应该是因为该算法的计算是基于概率分布，而概率分布意味着这是一个可以根据概率，产生无数样本的population，而非sample，因此对概率分布求均值就相当于求期望，方差也同理，需要带上概率]

原理

一幅图片的所有像素有$L$个灰度级，$[1,2,\cdots ,L]$
处于灰度级$i$的像素数量表示为$n_i$，因此像素的总数为$N=n_1+n_2+\cdots +n_L$。后续为了简化讨论，将灰度直方图归一化，因而就能被视为一个概率分布：
$$p_i=\frac{n_i}{N},p_i\ge ,\sum _{i=1}^L{p}_i=1(1)$$
现在，将所有像素以灰度级$k$为阈值，划分为两类$C_0$和$C_1$，$C_0$代表处于灰度级$[1,\cdots ,k]$的像素，$C_1$代表处于灰度级$[k+1,\cdots ,L]$的像素。
某一类别出现的概率为：
$${\omega }_0=Pr(C_0)=\sum _{i=1}^k{p}_i=\omega (k)(2)$$
$${\omega }_1=Pr(C_1)=\sum _{i=k+1}^L{p}_i=1-\omega (k)(3)$$
某一类别的平均灰度级为（注意是条件概率）：
$${\mu }_0=\sum _{i=1}^{k}iPr(i|C_0)=\sum _{i=1}^k\frac{i{p}_i}{{\omega }_0}=\frac{\mu (k)}{\omega (k)}(4)$$
其中$$\mu (k)=\sum _{i=1}^{k}i{p}_i(5)$$
$\mu (k)$可以看成灰度直方图到k的一阶矩，同理，公式(2)中的$\omega (k)$可以看成灰度直方图到k的零阶矩。
$${\mu }_1=\sum _{i=k+1}^{L}iPr(i|C_1)=\sum _{i=k+1}^L\frac{i{p}_i}{{\omega }_1}=\frac{{\mu }_T-\mu (k)}{1-\omega (k)}(6)$$
其中$${\mu }_T=\mu (L)=\sum _{i=1}^{L}i{p}_i(7)$$
${\mu }_T$是原图像的整体平均灰度级。根据上述公式，对于任意的k，我们有：
$${\omega }_0{\mu }_0+{\omega }_1{\mu }_1={\mu }_T,{\omega }_0+{\omega }_1=1(8)$$
(ps:这可以理解为全概率公式，也是公式(4),(6)用条件概率的原因，比如考虑人患艾滋病的概率$P(S)$，现在你可以把人群分的细一些，比如（男性$M$，女性$W$），它们都有自己的发生概率。当考虑一个完全的分类结果时候，你就得到了全概率公式。在 这里你会得到：患病概率 = 男性患病概率 * 男性人口比例 + 女性患病概率 * 女性人口比例。用概率公式写出来就是$P(S)=P(S|M)P(M)+P(S|W)P(W)$)
某一类的灰度级方差为：
$${\sigma }_0^2=\sum _{i=1}^k(i-{\mu }_0{)}^{2}Pr(i|C_0)=\sum _{i=1}^k(i-{\mu }_0{)}^2\frac{p_i}{{\omega }_0}(9)$$
$${\sigma }_1^2=\sum _{i=k+1}^L(i-{\mu }_1{)}^{2}Pr(i|C_1)=\sum _{i=k+1}^L(i-{\mu }_1{)}^2\frac{p_i}{{\omega }_1}(10)$$
接下来，为了衡量阈值k选取的好坏，我们首先引进线性判别式分析（Linear Discriminant Analysis, LDA）里的思想，希望选定的阈值能使不同的类之间相差越大越好，相同的类则相差越小越好(即最大化$\lambda$)，Otsu论文中的公式如下(但是在LDA里，我只看到了$\lambda$，不知道$\kappa$和$\eta$是本文提出来的还是怎样)：
$$\lambda =\frac{{\sigma }_B^2}{{\sigma }_W^2},\kappa =\frac{{\sigma }_T^2}{{\sigma }_W^2},\eta =\frac{{\sigma }_B^2}{{\sigma }_T^2}(11)$$
其中${\sigma }_W^2$、${\sigma }_B^2$、${\sigma }_T^2$分别代表类内(within-class)方差，类间(between-class)方差与所有灰度级的总(total levels)方差。计算方法如下：
$${\sigma }_W^2={\omega }_0{\sigma }_0^2+{\omega }_1{\sigma }_1^2(12)$$
$${\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(13)$$
$${\sigma }_T^2=\sum _{i=1}^L(i-{\mu }_T{)}^2{p}_i(14)$$
又因为${\sigma }_W^2+{\sigma }_B^2={\sigma }_T^2$,因此$\kappa =\lambda +1$，$\eta =\frac{\lambda }{\lambda +1}$，都是关于$\lambda$的函数，而${\sigma }_W^2$涉及二阶统计量（方差），${\sigma }_B^2$只涉及一阶统计量，因此，我们只需要选取其中最简单的$\eta$求最大值即可。由于${\sigma }_T^2$与k的选取无关，因此最后的优化目标等同于最大化${\sigma }_B^2$,类间方差。

代码

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