大津阈值分割(OSTU)

文章首发于我的个人博客

大津法是一种灰度图像自适应阈值分割算法,是日本学者Ostu于1979年提出,又称类间方差阈值分割法。大津法根据图像的灰度分布将图像分为前景和背景两部分,前景是我们分割出来的部分。前景和背景的分割值就是我们要通过类间方差法求出的阈值。

一、算法原理

设图像灰度级为L,灰度级为i的像素点数为ni,那么直方图分布为:
p i = n i / N , ∑ i = 0 L − 1 p i = 1 p_i = n_i /N, \sum_{i=0}^{L-1}p_i=1 p i ​=n i ​/N ,i =0 ∑L −1 ​p i ​=1
按灰度级用阈值t划分为两类:C_0=(0,1,…,t)和C_1=(t+1,t+2…,L-1)。因此C_0类和C_1类的概率级均值分别由一下各式给出:
w 0 = P r ( C 0 ) = ∑ i = 0 k p i w_0=Pr(C_0)=\sum_{i=0}^{k}p_i w 0 ​=P r (C 0 ​)=i =0 ∑k ​p i ​

w 1 = P r ( C 1 ) = ∑ i = k + 1 L − 1 p i = 1 − w 0 w_1=Pr(C_1)=\sum_{i=k+1}^{L-1}p_i = 1-w_0 w 1 ​=P r (C 1 ​)=i =k +1 ∑L −1 ​p i ​=1 −w 0 ​

u 0 = ∑ i = 0 k i P r ( i ∣ C 0 ) = ∑ i = 0 k i p i / w 0 u_0=\sum_{i=0}^{k}iPr(i|C_0)=\sum_{i=0}^{k}ip_i/w_0 u 0 ​=i =0 ∑k ​i P r (i ∣C 0 ​)=i =0 ∑k ​i p i ​/w 0 ​

u 1 = ∑ i = k + 1 L − 1 i P r ( i ∣ C 1 ) = ∑ i = k + 1 L − 1 i p i / w 1 u_1=\sum_{i=k+1}^{L-1}iPr(i|C_1)=\sum_{i=k+1}^{L-1}ip_i/w_1 u 1 ​=i =k +1 ∑L −1 ​i P r (i ∣C 1 ​)=i =k +1 ∑L −1 ​i p i ​/w 1 ​

u T = ∑ i = 0 L − 1 i p i u_T=\sum_{i=0}^{L-1}ip_i u T ​=i =0 ∑L −1 ​i p i ​

可以看出,对任何t值,下式都能成立:
w o u 0 + w 1 u 1 = u T , w 0 + w 1 = 1 w_ou_0+w_1u_1=u_T,w_0+w_1=1 w o ​u 0 ​+w 1 ​u 1 ​=u T ​,w 0 ​+w 1 ​=1
C_0和C_1类的方差可由下式求得:
σ 0 2 = ∑ i = 0 t ( i − u 0 ) 2 p i / w 0 \sigma_0^2=\sum_{i=0}^{t}(i-u_0)^2p_i/w_0 σ0 2 ​=i =0 ∑t ​(i −u 0 ​)2 p i ​/w 0 ​

σ 1 2 = ∑ i = t + 1 L − 1 ( i − u 1 ) 2 p i / w 1 \sigma_1^2=\sum_{i=t+1}^{L-1}(i-u_1)^2p_i/w_1 σ1 2 ​=i =t +1 ∑L −1 ​(i −u 1 ​)2 p i ​/w 1 ​

由此便可定义,类内方差(within):
σ w 2 = w 0 σ 0 2 + w 1 σ 1 2 \sigma_w^2=w_0\sigma_0^2+w_1\sigma_1^2 σw 2 ​=w 0 ​σ0 2 ​+w 1 ​σ1 2 ​
类间方差(Between):
σ B 2 = w 0 ( u 0 − u T ) 2 + w 1 ( u 1 − u T ) 2 = w 0 w 1 ( u 1 − u 0 ) 2 \sigma_B^2=w_0(u_0-u_T)^2+w_1(u_1-u_T)^2=w_0w_1(u_1-u_0)^2 σB 2 ​=w 0 ​(u 0 ​−u T ​)2 +w 1 ​(u 1 ​−u T ​)2 =w 0 ​w 1 ​(u 1 ​−u 0 ​)2
总体方差:
σ T 2 = σ B 2 + σ w 2 \sigma_T^2=\sigma_B^2+\sigma_w^2 σT 2 ​=σB 2 ​+σw 2 ​
何为类间方差?我们将灰度级L按t非为两类,反应在直方图上就是两类灰度。我们要做的就是找到能使两边灰度差的最大的阈值,根据这个阈值对灰度图像进行分割。因此我们要做的就是找到能使类间方差最大的t值。

二、类间方差推导

σ B 2 = w 0 ( μ 0 − μ T ) 2 + w 1 ( μ 1 − μ T ) 2 = w 0 [ μ 0 − w 0 μ 0 − w 1 μ 1 ] 2 + w 1 [ μ 1 − w 0 μ 0 − w 1 μ 1 ] 2 = w 0 [ ( 1 − w 0 ) μ 0 − w 1 μ 1 ] 2 + w 1 [ ( 1 − w 1 ) μ 1 − w 0 μ 0 ] 2 = w 0 [ w 1 μ 0 − w 1 μ 1 ] 2 + w 1 [ w 0 μ 1 − w 0 μ 0 ] 2 = w 0 w 1 ( μ 0 − μ 1 ) 2 + w 0 w 1 ( μ 1 − μ 0 ) 2 = ( μ 1 − μ 0 ) 2 w 0 w 1 ( w 0 + w 1 ) = w 0 w 1 ( μ 1 − μ 0 ) 2 {\sigma_B}^2=w_0(\mu_0-\mu_T)^2+w_1(\mu_1-\mu_T)^2 \ = w_0[\mu_0-w_0\mu_0-w_1\mu_1]^2+w_1[\mu_1-w_0\mu_0-w_1\mu_1]^2 \ = w_0[(1-w_0)\mu_0-w_1\mu_1]^2+w_1[(1-w_1)\mu_1-w_0\mu_0]^2 \ =w_0[w_1\mu_0-w_1\mu_1]^2+w_1[w_0\mu_1-w_0\mu_0]^2 \ =w_0w_1(\mu_0-\mu_1)^2+w_0w_1(\mu_1-\mu_0)^2 \ =(\mu_1-\mu_0)^2w_0w_1(w_0+w_1)\ =w_0w_1(\mu_1-\mu_0)^2 σB ​2 =w 0 ​(μ0 ​−μT ​)2 +w 1 ​(μ1 ​−μT ​)2 =w 0 ​[μ0 ​−w 0 ​μ0 ​−w 1 ​μ1 ​]2 +w 1 ​[μ1 ​−w 0 ​μ0 ​−w 1 ​μ1 ​]2 =w 0 ​[(1 −w 0 ​)μ0 ​−w 1 ​μ1 ​]2 +w 1 ​[(1 −w 1 ​)μ1 ​−w 0 ​μ0 ​]2 =w 0 ​[w 1 ​μ0 ​−w 1 ​μ1 ​]2 +w 1 ​[w 0 ​μ1 ​−w 0 ​μ0 ​]2 =w 0 ​w 1 ​(μ0 ​−μ1 ​)2 +w 0 ​w 1 ​(μ1 ​−μ0 ​)2 =(μ1 ​−μ0 ​)2 w 0 ​w 1 ​(w 0 ​+w 1 ​)=w 0 ​w 1 ​(μ1 ​−μ0 ​)2

三、c++实现(核心代码)

    auto tempImage = originImage.clone();
    std::vector<int> countPixel(256, 0);

    std::vector<float> pixelProb;
    std::vector<int> pixel(256);
    for(int i = 0; i < 256; i++){
        pixel[i] = i;
    }

    for(int i = 0; i < tempImage.rows; i++){
        for(int j = 0; j < tempImage.cols; j++){
            countPixel[tempImage.at<uchar>(i, j)]++;
        }
    }

    float sum_pixel = tempImage.rows * tempImage.cols;
    for(auto iter = countPixel.begin(); iter != countPixel.end(); iter++){
        pixelProb.push_back(*iter / sum_pixel);
    }

    std::vector<float> class_within_variance_vec;

    std::vector<float> class_between_variance_vec;

    for(int k = 0; k < 256; k++){

        float class_0_prob = 0;

        float class_1_prob = 0;
        for(int i = 0; i  k; i++){
            class_0_prob += pixelProb[i];
        }
        class_1_prob = 1 - class_0_prob;

        float class_0_average = 0;

        float class_1_average = 0;
        for(int i = 0; i  k; i++){
            class_0_average += i * (pixelProb[i] / (class_0_prob+0.000001));
        }
        for(int i = k+1; i  255; i++){
            class_1_average += i * (pixelProb[i] / (class_1_prob+0.000001));
        }

        float class_0_variance = 0;

        float class_1_variance = 0;
        for(int i = 0; i  k; i++){
            class_0_variance += std::pow((i-class_0_average),2)*(pixelProb[i]/(class_0_prob+0.000001));
        }

        for(int i = k+1; i  255; i++){
            class_1_variance += std::pow((i-class_1_average),2)*(pixelProb[i]/(class_1_prob+0.000001));
        }

        float class_within_variance = class_0_prob*std::pow(class_0_variance,2) + class_1_prob*std::pow(class_1_variance,2);
        class_within_variance_vec.push_back(class_within_variance);

        float class_between_variance = class_0_prob * class_1_prob * std::pow((class_1_average-class_0_average),2);
        class_between_variance_vec.push_back(class_between_variance);
    }

    auto max_value_iter = std::max_element(class_between_variance_vec.begin(), class_between_variance_vec.end());

    int k_value = std::distance(class_between_variance_vec.begin(), max_value_iter);
    std::cout<<"k_value:"<< k_value<<"  max_var:"<<*max_value_iter<<std::endl;

    for(int i = 0; i < tempImage.rows; i++){
        for(int j = 0; j < tempImage.cols; j++){
            if(tempImage.at<uchar>(i, j)  k_value)
                tempImage.at<uchar>(i, j) = 0;
            else
                tempImage.at<uchar>(i, j) = 255;
        }
    }

四、效果对比

下图从左到右依次为,原灰度图、OpenCV提供OSTU算法分割结果、自己实现的OSTU分割效果。

大津阈值分割(OSTU)

下图为调用matlab函数实现的大津法阈值分割效果

大津阈值分割(OSTU)

必不可少的Lenna。

大津阈值分割(OSTU)

大津阈值分割(OSTU)

Original: https://blog.csdn.net/qq_42780025/article/details/113000940
Author: nachr
Title: 大津阈值分割(OSTU)

原创文章受到原创版权保护。转载请注明出处:https://www.johngo689.com/193724/

转载文章受原作者版权保护。转载请注明原作者出处!

(0)

大家都在看

免费咨询
免费咨询
扫码关注
扫码关注
联系站长

站长Johngo!

大数据和算法重度研究者!

持续产出大数据、算法、LeetCode干货,以及业界好资源!

2022012703491714

微信来撩,免费咨询:xiaozhu_tec

分享本页
返回顶部