灰度图像--图像分割阈值处理之OTSU阈值_ostu阈值分割提取清晰轮廓-CSDN博客

本文链接：https://blog.csdn.net/TonyShengTan/article/details/44132273

学习DIP第55天
转载请标明本文出处：***http://blog.csdn.net/tonyshengtan ***，出于尊重文章作者的劳动，转载请标明出处！文章代码已托管，欢迎共同开发：https://github.com/Tony-Tan/DIPpro
更多图像处理机器学习内容请访问最新网站www.tony4ai.com
#开篇废话
废话开始，今天介绍OTSU算法，本算法比前面给出的算法更能够给出数学上的最佳阈值，不需要任何输入附加参数、与同样不需要输入附加参数的迭代均值和均值阈值来比较，OTSU给出的阈值能使分类更加均匀。
阈值处理将灰度值分为两类，而对于分类问题，已有的一种最优闭合解–贝叶斯决策规则。
#贝叶斯决策规则
首先介绍下贝叶斯公式的形象化理解，考虑下图
这里写图片描述

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

$p(A)=\frac{8}{12}=\frac{2}{3} $..........(1) $p(B)=\frac{7}{12} $..........(2) $p(AB)=\frac{3}{12}=\frac{1}{4} $..........(3)

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

$p(B|A)=\frac{3}{8} $..........(4)

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

$p(A|B)=\frac{3}{7} $..........(5)

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

$p(AB)=p(A|B)*p(B)=p(B|A)*p(A) $..........(6)

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

$p(A)=\frac{2}{3} $

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

$p(B|A)=\frac{3}{8} $

拿出弹夹的概率：

$p(B)=\frac{7}{12} $

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

$p(A|B)=\frac{p(B|A)*p(A)}{p(B)} $

以上就是贝叶斯公式的一般形式，更复杂的形式会在后面的文章中详细介绍。（更复杂的形式是指盒子里有枪，子弹，弹夹，手榴弹。。。。。。）
#数学原理
OTSU算法可以基于直方图计算，考虑灰度级为{0，1，2…L-1}大小为$M \times N $的图像，设$ n_i $为灰度级为i的像素的总数量，那么:

$M \times N=\sum^{L-1}_{i=0}n_i $ $p(n_i)=\frac{n_i}{M \times N} $ $\sum^{L-1}_{i=0}p_i=1 $

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

$p_1(k)=\sum^{k}_{i=0}p_i $

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

$p_2(k)=\sum^{L-1}_{i=k+1}p_i $

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

$m_1(k)=\sum^{k}_{i=0}i*P(i|C_1)=\sum^{k}_{i=0}i*\frac{P(C_1|i)*P(i)}{P(C_1)} $

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

$m_1(k)=\frac{1}{P_1(k)} \sum^{k}_{i=0}i*p_i $

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

$m_2(k)=\frac{1}{P_2(k)} \sum^{L-1}_{k+1}i*p_i $

全图的均值

$m_G=\sum^{L-1}_{i=0}iP_i $

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

$P_1m_1+P_2m_2=m_G $ $P_1+P_2=1 $

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

$\eta=\frac{\delta_B^2}{\delta_G^2} $ $\delta_G^2=\sum^{L-1}_{i=0}(i-m_G)^2*p_i $

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

$\delta_B^2=P_1(m_1-m_G)^2+P_2(m_2-m_G)^2 $

公式还可以写成：

$\delta^2_B=P_1P_2(m_1-m_2)^2=\frac{P_1(m_1-m_G)^2}{1-P_1} $

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

$\delta^2_B(k^*)=max_{0\leq k \leq L-1}\delta^2_B(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);
}

#观察结果
原图：

![这里写图片描述](https://img-blog.csdn.net/20150308152743862) ![这里写图片描述](https://img-blog.csdn.net/20150308152907539)

加入1%的高斯噪声：

![这里写图片描述](https://img-blog.csdn.net/20150308152918272) ![这里写图片描述](https://img-blog.csdn.net/20150308152814765)

加入3%的高斯噪声：

![这里写图片描述](https://img-blog.csdn.net/20150308152827323) ![这里写图片描述](https://img-blog.csdn.net/20150308152836777)

加入5%的高斯噪声：

![这里写图片描述](https://img-blog.csdn.net/20150308153000860) ![这里写图片描述](https://img-blog.csdn.net/20150308153012264)

加入7%的高斯噪声：

![这里写图片描述](https://img-blog.csdn.net/20150308152915106) ![这里写图片描述](https://img-blog.csdn.net/20150308152932001)

加入9%的高斯噪声：

![这里写图片描述](https://img-blog.csdn.net/20150308153053245) ![这里写图片描述](https://img-blog.csdn.net/20150308153102277)