学习DIP第55天
转载请标明本文出处：http://blog.csdn.net/tonyshengtan ，出于尊重文章作者的劳动，转载请标明出处！文章代码已托管，欢迎共同开发：https://github.com/Tony-Tan/DIPpro

开篇废话

废话开始，今天介绍OTSU算法，本算法比前面给出的算法更能够给出数学上的最佳阈值，不需要任何输入附加参数、与同样不需要输入附加参数的迭代均值和均值阈值来比较，OTSU给出的阈值能使分类更加均匀。
阈值处理将灰度值分为两类，而对于分类问题，已有的一种最优闭合解–贝叶斯决策规则。

贝叶斯决策规则

首先介绍下贝叶斯公式的形象化理解，考虑下图

上面的12幅图中有手枪和弹夹，只有弹夹和手枪出现在同一个盒子的时候才有杀伤力，也就是你拿到一个盒子，你不知道里面是什么，有可能是枪，有可能是弹夹，有可能同时有枪和弹夹。下面来从概率学角度分析
设盒子里有枪为事件A，那么A出现的概率设为$p(A)$。
设盒子里有弹夹为事件B，那么B出现的概率设为$p(B)$。
那么同时出现事件A和事件B的概率为$p(AB)$
看图可以知道

考虑我们随机抽出一个盒子，先拿出一个东西，比如先拿出一把枪，那么也就是事件A发生了，那么我们继续从盒子里拿东西，有可能拿到弹夹，也有可能啥也没有，那么拿到弹夹的概率就如下：

同理，如果先拿出来的是个弹夹，那么接下来拿出枪的概率是：

结合(1)(2)(3)(4)(5)，可以得到：

假设下面情形：
已知拿出枪的概率是：

拿出枪以后拿出弹夹的概率

拿出弹夹的概率：

求拿出弹夹以后拿出枪的概率

以上就是贝叶斯公式的一般形式，更复杂的形式会在后面的文章中详细介绍。（更复杂的形式是指盒子里有枪，子弹，弹夹，手榴弹。。。。。。）

数学原理

OTSU算法可以基于直方图计算，考虑灰度级为{0，1，2……..L-1}大小为$M \times N$的图像，设$n_i$为灰度级为i的像素的总数量，那么:

假设阈值为k将直方图分成两部分。
部分1$(C_1)$的概率为：

部分2$(C_2)$的概率为：

部分1$(C_1)$的平均数：

$P(C_1|i)$的值为1，因为$i$是属于$C_1$的，所以发生$i$以后发生$C_1$的概率是100%，所以

部分2$(C_2)$的平均数：

全图的均值

上面的式子可以由下面验证：

下面就是关键部分了，如何评价一个阈值的好坏，提出一个阈值，将像素灰度分为两类，通过以下的公式来评价阈值质量：

$\delta_B^2$是类间方差，其定义为：

公式还可以写成：

于是最佳阈值$k^*$由下面得出：

通过上式可以通过迭代计算出最佳的k值。使用k作为阈值，对图像进行处理。

代码实现

/*
 *OTSU 算法
 *otsu 算法使用贝叶斯分类原理得到最好聚类
 *
 *
 */
//归一化直方图

void setHist2One(double *hist_d,double *dst_hist_d){
    double sum=0.0;
    for(int i=0;i<GRAY_LEVEL;i++)
        sum+=hist_d[i];
    if(sum!=0)
        for(int i=0;i<GRAY_LEVEL;i++)
            dst_hist_d[i]=hist_d[i]/sum;

}
//计算公式中最大的deta，并返回直方图灰度
double findMaxDeta(double *hist_d){
    double max_deta=-1.0;
    double max_deta_location=0.0;
    double m_g=0.0;

    for(int i=0;i<GRAY_LEVEL;i++)
        m_g+=i*hist_d[i];


    for(int i=0;i<GRAY_LEVEL;i++){
        double p1=0.0;
        double m1=0.0;
        double deta=0.0;
        for(int j=0;j<=i;j++){
            p1+=hist_d[j];
            m1+=j*hist_d[j];
        }
        deta=p1*(m1-m_g)*(m1-m_g)/(1-p1);
        if(deta>max_deta){
            max_deta_location=i;
            max_deta=deta;
        }
    }
    return max_deta_location;
}
void OTSUThreshold(double *src,double *dst,int width,int height,int type){
    int hist[GRAY_LEVEL];
    double hist_d[GRAY_LEVEL];
    setHistogram(src, hist, width, height);
    Hist_int2double(hist, hist_d);
    setHist2One(hist_d, hist_d);
    double threshold=findMaxDeta(hist_d);
    Threshold(src, dst, width, height, threshold, type);
}

观察结果

原图：

加入1%的高斯噪声：

加入3%的高斯噪声：

加入5%的高斯噪声：

加入7%的高斯噪声：

加入9%的高斯噪声：

加入11%的高斯噪声：

lena:

baboon:

总结

OTSU算法产生的阈值是数学角度上的最佳分类，数学基础的贝叶斯公式，但应用也有一定的局限性，比如，前面说过最多的，对全局阈值，目标与背景的大小关系，当目标和背景大小相差很多时，或者噪声很大的时候，对OTSU产生影响较大。
待续。。。