大津阈值法(OTSU)的应用

大津算法参见 点击打开链接

最大类间方差法是由日本学者大进展之于1979年提出的,是一种自适应的阈值确定的方法,又叫大津法,简称OTSU。它是按图像的灰度特性,将图像分成背景和目标两部分。背景和目标之间的类间方差越大,说明构成图像的两部分的差别越大,当部分目标错分为背景或部分背景错分为目标都会导致两部分差别变小。因此,使类间方差最大的分割意味着错分概率最小。对于图像I(x,y),前景(即目标)和背景的分割阈值记作T,属于前景的像素点数占整幅图像的比例记为\omega_0 \,,其平均灰度\mu_0;背景像素点数占整幅图像的比例为\omega_1,其平均灰度为\mu_1。图像的总平均灰度记为\mu,类间方差记为g。假设图像的背景较暗,并且图像的大小为M*N,图像中像素的灰度值小于阈值T的像素个数记作N_0,像素灰度大于阈值T的像素个数记作N_1,则有:

 \omega_0=\frac{N_0}{M*N}  (1)
  \omega_1=\frac{N_1}{M*N}             (2)
  N_0+N_1=M*N           (3)
  \omega_0+\omega_1=1           (4)
  \mu=\omega_0*\mu_0+\omega_1*\mu_1             (5)
  g=\omega_0(\mu_0-\mu)^2+\omega_1(\mu_1-\mu)^2       (6)

将式(5)代入式(6),得到等价公式:

  g=\omega_0\omega_1(\mu_0-\mu_1)^2              (7)

采用遍历的方法得到使类间方差最大的阈值T,即为所求。

在前面的博文中提到利用Sobel算子得到图像的梯度信息。那么本文将用大津算法将梯度幅值分割为背景和前景两部分,那么前景部分即为图像的边缘信息。

OpenCV中源码如下:

  1. getThreshVal_Otsu_8u( const Mat& _src )  
  2. {  
  3.     Size size = _src.size();  
  4.     if( _src.isContinuous() )  
  5.     {  
  6.         size.width *= size.height;  
  7.         size.height = 1;  
  8.     }  
  9.     const int N = 256;  
  10.     int i, j, h[N] = {0};  
  11.     for( i = 0; i < size.height; i++ )  
  12.     {  
  13.         const uchar* src = _src.data + _src.step*i;  
  14.         j = 0;  
  15.         #if CV_ENABLE_UNROLLED  
  16.         for( ; j <= size.width - 4; j += 4 )  
  17.         {  
  18.             int v0 = src[j], v1 = src[j+1];  
  19.             h[v0]++; h[v1]++;  
  20.             v0 = src[j+2]; v1 = src[j+3];  
  21.             h[v0]++; h[v1]++;  
  22.         }  
  23.         #endif  
  24.         for( ; j < size.width; j++ )  
  25.             h[src[j]]++;  
  26.     }  
  27.   
  28.     double mu = 0, scale = 1./(size.width*size.height);  
  29.     for( i = 0; i < N; i++ )  
  30.         mu += i*(double)h[i];  
  31.   
  32.     mu *= scale;  
  33.     double mu1 = 0, q1 = 0;  
  34.     double max_sigma = 0, max_val = 0;  
  35.   
  36.     for( i = 0; i < N; i++ )  
  37.     {  
  38.         double p_i, q2, mu2, sigma;  
  39.   
  40.         p_i = h[i]*scale;  
  41.         mu1 *= q1;  
  42.         q1 += p_i;  
  43.         q2 = 1. - q1;  
  44.   
  45.         if( std::min(q1,q2) < FLT_EPSILON || std::max(q1,q2) > 1. - FLT_EPSILON )  
  46.             continue;  
  47.   
  48.         mu1 = (mu1 + i*p_i)/q1;  
  49.         mu2 = (mu - q1*mu1)/q2;  
  50.         sigma = q1*q2*(mu1 - mu2)*(mu1 - mu2);  
  51.         if( sigma > max_sigma )  
  52.         {  
  53.             max_sigma = sigma;  
  54.             max_val = i;  
  55.         }  
  56.     }  
  57.   
  58.     return max_val;  
  59. }  
getThreshVal_Otsu_8u( const Mat& _src )
{
    Size size = _src.size();
    if( _src.isContinuous() )
    {
        size.width *= size.height;
        size.height = 1;
    }
    const int N = 256;
    int i, j, h[N] = {0};
    for( i = 0; i < size.height; i++ )
    {
        const uchar* src = _src.data + _src.step*i;
        j = 0;
        #if CV_ENABLE_UNROLLED
        for( ; j <= size.width - 4; j += 4 )
        {
            int v0 = src[j], v1 = src[j+1];
            h[v0]++; h[v1]++;
            v0 = src[j+2]; v1 = src[j+3];
            h[v0]++; h[v1]++;
        }
        #endif
        for( ; j < size.width; j++ )
            h[src[j]]++;
    }

    double mu = 0, scale = 1./(size.width*size.height);
    for( i = 0; i < N; i++ )
        mu += i*(double)h[i];

    mu *= scale;
    double mu1 = 0, q1 = 0;
    double max_sigma = 0, max_val = 0;

    for( i = 0; i < N; i++ )
    {
        double p_i, q2, mu2, sigma;

        p_i = h[i]*scale;
        mu1 *= q1;
        q1 += p_i;
        q2 = 1. - q1;

        if( std::min(q1,q2) < FLT_EPSILON || std::max(q1,q2) > 1. - FLT_EPSILON )
            continue;

        mu1 = (mu1 + i*p_i)/q1;
        mu2 = (mu - q1*mu1)/q2;
        sigma = q1*q2*(mu1 - mu2)*(mu1 - mu2);
        if( sigma > max_sigma )
        {
            max_sigma = sigma;
            max_val = i;
        }
    }

    return max_val;
}

参见 最大类间方差法(大津法,OTSU)中的C++代码,与OpenCV自带的大津阈值函数进行比较如下。

原图:


测试源码如下:

  1. #include “stdafx.h”  
  2. #include <opencv2\highgui\highgui.hpp>  
  3. #include <cv.h>  
  4.   
  5. using namespace cv;  
  6.   
  7. int _tmain(int argc, _TCHAR* argv[])  
  8. {  
  9.     Mat src = imread(”E:\\VC++Projects\\road.jpg”);  
  10.     Mat gray;  
  11.     cvtColor(src, gray, CV_BGR2GRAY);  
  12.   
  13.     // 求得x和y方向的一阶微分  
  14.     Mat sobelx;  
  15.     Mat sobely;  
  16.     Sobel(gray, sobelx, CV_32F, 1, 0, 3);  
  17.     Sobel(gray, sobely, CV_32F, 0, 1, 3);  
  18.   
  19.     // 求得梯度和方向  
  20.     Mat norm;  
  21.     Mat dir;  
  22.     cartToPolar(sobelx, sobely, norm, dir);  
  23.   
  24.     // 转换为8位单通道图像  
  25.     double normMax;  
  26.     minMaxLoc(norm, NULL, &normMax);  
  27.     Mat grad;  
  28.     norm.convertTo(grad, CV_8UC1, 255.0/normMax, 0);  
  29.     namedWindow(”grad”);  
  30.     imshow(”grad”,grad);  
  31.   
  32.     // OpenCV自带大津阈值函数处理  
  33.     Mat cvOstu;  
  34.     threshold(grad, cvOstu, 0, 255, CV_THRESH_BINARY|CV_THRESH_OTSU);  
  35.     namedWindow(”cvOstu”);  
  36.     imshow(”cvOstu”, cvOstu);  
  37.   
  38.     // 自己实现方法  
  39.     Mat myOstu;  
  40.     int width = grad.cols;  
  41.     int height = grad.rows;  
  42.     int pixelCount[256] = {0};  
  43.     float pixelPro[256] = {0};  
  44.     int i, j, pixelSum = width * height, Threshold = 0;  
  45.   
  46.     // 统计直方图  
  47.     for(i=0; i<height; i++)  
  48.     {  
  49.         uchar* g_ptr = grad.ptr<uchar>(i);  
  50.         for(j=0; j<width; j++)  
  51.         {  
  52.             pixelCount[ g_ptr[j] ]++;  
  53.         }  
  54.     }  
  55.   
  56.     // 每个灰度值的比例  
  57.     for(i=0; i<256; i++)  
  58.     {  
  59.         pixelPro[i] = (float)pixelCount[i] / pixelSum;  
  60.     }  
  61.   
  62.     // 计算最大类间方差  
  63.     float w0, w1, u0tmp, u1tmp, u0, u1, deltaTmp, deltaMax = 0;    
  64.     for(i = 0; i < 256; i++)    
  65.     {    
  66.         w0 = w1 = u0tmp = u1tmp = u0 = u1 = deltaTmp = 0;    
  67.         for(j = 0; j < 256; j++)    
  68.         {    
  69.             if(j <= i)       //背景部分    
  70.             {    
  71.                 w0 += pixelPro[j];    
  72.                 u0tmp += j * pixelPro[j];    
  73.             }    
  74.             else            //前景部分    
  75.             {    
  76.                 w1 += pixelPro[j];    
  77.                 u1tmp += j * pixelPro[j];    
  78.             }    
  79.         }    
  80.         u0 = u0tmp / w0;    
  81.         u1 = u1tmp / w1;    
  82.         deltaTmp = (float)(w0 *w1* pow((u0 - u1), 2)) ;    
  83.         if(deltaTmp > deltaMax)    
  84.         {    
  85.             deltaMax = deltaTmp;    
  86.             Threshold = i;    
  87.         }    
  88.     }  
  89.   
  90.     // 阈值处理  
  91.     threshold(grad, myOstu, Threshold, 255, CV_THRESH_BINARY);  
  92.     namedWindow(”myOstu”);  
  93.     imshow(”myOstu”, myOstu);  
  94.       
  95.     waitKey(0);  
  96.   
  97.     return 0;  
  98. }  
#include "stdafx.h"




#include <opencv2\highgui\highgui.hpp> #include <cv.h> using namespace cv; int _tmain(int argc, _TCHAR* argv[]) { Mat src = imread("E:\\VC++Projects\\road.jpg"); Mat gray; cvtColor(src, gray, CV_BGR2GRAY); // 求得x和y方向的一阶微分 Mat sobelx; Mat sobely; Sobel(gray, sobelx, CV_32F, 1, 0, 3); Sobel(gray, sobely, CV_32F, 0, 1, 3); // 求得梯度和方向 Mat norm; Mat dir; cartToPolar(sobelx, sobely, norm, dir); // 转换为8位单通道图像 double normMax; minMaxLoc(norm, NULL, &normMax); Mat grad; norm.convertTo(grad, CV_8UC1, 255.0/normMax, 0); namedWindow("grad"); imshow("grad",grad); // OpenCV自带大津阈值函数处理 Mat cvOstu; threshold(grad, cvOstu, 0, 255, CV_THRESH_BINARY|CV_THRESH_OTSU); namedWindow("cvOstu"); imshow("cvOstu", cvOstu); // 自己实现方法 Mat myOstu; int width = grad.cols; int height = grad.rows; int pixelCount[256] = {0}; float pixelPro[256] = {0}; int i, j, pixelSum = width * height, Threshold = 0; // 统计直方图 for(i=0; i<height; i++) { uchar* g_ptr = grad.ptr<uchar>(i); for(j=0; j<width; j++) { pixelCount[ g_ptr[j] ]++; } } // 每个灰度值的比例 for(i=0; i<256; i++) { pixelPro[i] = (float)pixelCount[i] / pixelSum; } // 计算最大类间方差 float w0, w1, u0tmp, u1tmp, u0, u1, deltaTmp, deltaMax = 0; for(i = 0; i < 256; i++) { w0 = w1 = u0tmp = u1tmp = u0 = u1 = deltaTmp = 0; for(j = 0; j < 256; j++) { if(j <= i) //背景部分 { w0 += pixelPro[j]; u0tmp += j * pixelPro[j]; } else //前景部分 { w1 += pixelPro[j]; u1tmp += j * pixelPro[j]; } } u0 = u0tmp / w0; u1 = u1tmp / w1; deltaTmp = (float)(w0 *w1* pow((u0 - u1), 2)) ; if(deltaTmp > deltaMax) { deltaMax = deltaTmp; Threshold = i; } } // 阈值处理 threshold(grad, myOstu, Threshold, 255, CV_THRESH_BINARY); namedWindow("myOstu"); imshow("myOstu", myOstu); waitKey(0); return 0; }
测试结果,梯度幅值图:


OpenCV自带函数处理结果:


参见网上源码的处理结果:


因此可见大津法可以很好的对车道线进行分割的作用,同时也验证了参见源码的可靠性。



  • 2
    点赞
  • 4
    收藏
    觉得还不错? 一键收藏
  • 0
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值