详解-OTUS(大津法-最大类间方差)原理及C语言代码实现

Posted Z小旋

tags:

篇首语:本文由小常识网(cha138.com)小编为大家整理,主要介绍了详解-OTUS(大津法-最大类间方差)原理及C语言代码实现相关的知识,希望对你有一定的参考价值。

灰度图二值化:

我们在对灰度图像进行处理的时候,为了便于观察和分析,经常需要将图像中的目标主体和背景分割开来,变成二值化图像(只有黑和白,黑白图像)

加菲猫--灰度图——二值化

 

而我们知道灰度图像是有256个灰度级 255代表全白,0表示全黑,那么在进行二值化的时候,是设定一个阈值,根据灰度值大于或小于阈值进行黑白显示,我们假设背景用白色0表示,目标物体用黑色1表示,不同阈值的选取,对于图像二值化的效果影响是非常大的

上图可以看出,阈值的选取对于灰度图二值化有着至关重要的作用,那么怎么选取一个合适的阈值,使得我们的背景和前景能够很好的分开,并且最大程度的减少误判呢?

灰度直方图

这就需要我们对图像的灰度直方图进行分析,

灰度直方图就是图像的像素点在不同灰度级的分布情况

看下面的图片,横轴代表灰度级,纵轴代表像素点数量

从上图中我们可以看出该图片共有两个尖峰,第一个尖峰表示人的像素灰度值分布在0~ 20之间,第二个尖峰表示背景部分像素灰度值分布在灰度100~180之间 ,这个时候做二值化,我们可以在30-100之间设定一个阈值,大于阈值的背景像素全为0(白色),小于阈值的目标像素全为1(黑色) 这就是图像二值化

那么阈值取多少合适呢?观察直方图我们可以得到一个大致的区间,但是实际上的取值不能确定,所以我们需要一个算法来找到最优阈值才行。而大津法(OTUS),就给我们提供了怎么找到这个最佳阈值的方法。

图一为阈值120 图二为阈值75

我们可以看到,图一的特征值更多,图二只有人像部分,其余背景过滤较好。

类间方差

首先我们要知道什么是方差,根据百度百科:方差是一组数据离散程度的度量,方差越大,表示离散程度越大

请注意上面这句话, 既然我们可以方差来表示系统的离散程度,我们图像二值化取最佳阈值时,背景应该与前景差别最大,这时方差也应该是最大的,那这样,是不是可以用方差来表示当前图像的二值化质量?

所以,我们假设图像阈值为T

小于T的像素为目标,数量占总图像比例为 w0 ,平均灰度值为u0

大于T的像素为背景,数量占图像比例为 w1,平均灰度值为u1

那么 w0+w1=1 (公式1) 概率和为1

 

图像的平均灰度值u为:(公式2)

那么类内方差g的定义为:(公式3)

(注意现在还是类内方差,两个类别的类内方差按比例求和)

 

目标部分和平均灰度值的方差,加上背景部分和平均灰度值的方差,这个值最大的时候,表示当前阈值为最佳阈值。

那么我们把上面的公式3简化一下,可以得到类间方差的定义:(公式4)

我们现在需要做的事,就是修改阈值T的值,求方差g的最大值。

推导过程:

大津法(OTUS)

这个时候来看大津法的概念:

最大类间方差是由日本学者大津(Nobuyuki Otsu)于1979年提出,是一种确定图像二值化分割阈值的算法。算法假设图像像素能够根据全局阈值,被分成背景[background]和目标[objects]两部分。然后,计算该最佳阈值来区分这两类像素,使得两类像素区分度最大。

大津法阈值采用最大类间方差的原理,适合于图像灰度分布整体呈现“双峰”的情况。大津法会自动找出一个阈值,使得分割后的两部分类间方差最大。

特性:

  1. 大津法对噪音十分敏感,在处理之前应对图片进行去噪处理。

    如果图像有存在局部噪声,则会影响大津法的判断

  2. 当目标与背景的面积比例悬殊的时候,类间方差函数可能呈现双峰或者多峰,这个时候 大津法的效果不好

双峰图像,目标和背景面积差距不大,可以很好的判断

当图像中的目标与背景的面积相差很大时,灰度直方图没有明显的双峰,或者两个峰的大小相差很大,分割效果不佳

目标像素点远小于背景像素点,分割效果不佳

代码实现:

基本所有的博客都是讲完原理,然后贴一整个代码,这样很多读者可能搞懂了原理,但是在阅读代码的时候突然又看不懂了,所以代码部分会按照上面的公式,推导编写,下面以C语言为例:

首先建立一个函数,定义需要用到的参数:

//------------------摄像头参数--------------//
uint8 dis_image[60][80];

大津法二值化//

/*! 

 *  @brief      大津法二值化0.8ms程序
 *  @date:   2018-9  
 *  @since      v1.2
 *  *image :图像地址
 *  width:  图像宽
 *  height:图像高
 *  @author     Z小旋
    */

uint8 image_threshold ;  //图像阈值

uint8 otsuThreshold(uint8 *image, uint16 col, uint16 row)
{
    #define GrayScale 256//定义256个灰度级
    uint16 width = col;   //图像宽度
    uint16 height = row;  //图像长度
    int pixelCount[GrayScale];  //每个灰度值所占像素个数
    float pixelPro[GrayScale]; //每个灰度值所占总像素比例
    int i, j;
    int sumPixel = width * height;//总像素点
    uint8 threshold = 0; //最佳阈值
    uint8* data = image;  //指向像素数据的指针

这里需要注意,因为我们定义了两个局部变量数组,局部变量处于堆栈区,数值随机,所以需要对数组做一个简单的初始化

初始化的方法有三种:

1用for循环赋值

for (i = 0; i < GrayScale; i++){
​        pixelCount[i] = 0;
​        pixelPro[i] = 0;}

2使用memset

memset(pixelCount, 0, GrayScale); //使用memset方法 

memset(pixelPro, 0, GrayScale); //使用memset方法 

3声明时,使用 {0} 初始化

int pixelCount[GrayScale] = {0};//每个灰度值所占像素个数
float pixelPro[GrayScale] = {0};//每个灰度值所占总像素比例

上面三种初始化,实测for循环浪费的时间最多,{0} 与memset 耗时差不多。

 
 
 

然后我们需要统计每个灰度级在整个图像中的个数,就是统计灰度值为1的点有多少多少个,灰度值为2的点有多少多少个…

      //统计灰度级中每个像素在整幅图像中的个数  
        for (i = 0; i < height; i++)
        {
            for (j = 0; j < width; j++)
            {
                pixelCount[(int)data[i * width + j]]++;  //将像素值作为计数数组的下标
                 //   pixelCount[(int)image[i][j]]++;    若不用指针用这个
            }
        }

 
 
 

计算每个灰度级像素数占图像总像素的比例,以及图像的总平均灰度u


    //遍历灰度级[0,255]  
        float w0=0, w1=0, u0Sum=0, u1Sum=0, u0=0, u1=0, u=0, variance=0, maxVariance = 0;
            //目标像素占总图像比例w0
            //背景像素占总图像比例w1 
            //目标平均灰度值u0 
            //背景平均灰度值u1
            //总平均灰度值u
            //类间方差 variance
            //最大类间方差 maxVariance
    
        //计算每个像素在整幅图像中的比例  
        for (i = 0; i < GrayScale; i++)
        {
            pixelPro[i] = (float)pixelCount[i] / sumPixel;
             u += i * pixelPro[i];  //总平均灰度
        }
     

接下来,就是我们假设图像阈值T从1-255遍历,然后求每次的目标像素数量占总图像比例w0,平均灰度值u0,背景像素数量占图像比例 w1,平均灰度值u1

然后比较256次情况下

(公式3)

公式g的最大值


    //计算每个像素在整幅图像中的比例  
    for (i = 0; i < GrayScale; i++)
    {
        pixelPro[i] = (float)pixelCount[i] / sumPixel;
         u += i * pixelPro[i];  //总平均灰度
    }

    for (i = 0; i < GrayScale; i++)     // i作为阈值 阈值从1-255遍历 
    {
    
        for (j = 0; j < GrayScale; j++)    //求目标前景和背景
        {
            if (j <= i)   //前景部分    
            {
                w0 += pixelPro[j];   
                u0Sum += j * pixelPro[j];
            }
            else   //背景部分  
            {
                w1 += pixelPro[j];
                u1Sum += j * pixelPro[j];
            }
        }
        u0 = u0Sum / w0;
        u1 = u1Sum / w1;
    
        variance = w0 * pow((u0 - u), 2) + w1 * pow((u1 - u), 2);  //类间方差计算公式
    
        if (variance > maxVariance)   //判断是否为最大类间方差
        {
            maxVariance = variance;
            threshold = i;
        }
    }
    
    return threshold;
    
    }

完整代码如下:

//------------------摄像头参数--------------//
uint8 dis_image[60][80];

大津法二值化//
uint8 image_threshold = 46;  //图像阈值

   //0x4D;0x18-0x1A;
uint8 otsuThreshold(uint8 *image, uint16 col, uint16 row)
{
    #define GrayScale 256//定义256个灰度级
    uint16 width = col;   //图像宽度
    uint16 height = row;  //图像长度
    int pixelCount[GrayScale];  //每个灰度值所占像素个数
    float pixelPro[GrayScale]; //每个灰度值所占总像素比例
    int i, j;
    int sumPixel = width * height;//总像素点
    uint8 threshold = 0; //最佳阈值
    uint8* data = image;  //指向像素数据的指针

    for (i = 0; i < GrayScale; i++)
    {
        pixelCount[i] = 0;
        pixelPro[i] = 0;
    }

    //统计灰度级中每个像素在整幅图像中的个数  
    for (i = 0; i < height; i++)
    {
        for (j = 0; j < width; j++)
        {
            pixelCount[(int)data[i * width + j]]++;  //将像素值作为计数数组的下标
             //   pixelCount[(int)image[i][j]]++;    若不用指针用这个
        }
    }

 //遍历灰度级[0,255]  
    float w0=0, w1=0, u0Sum=0, u1Sum=0, u0=0, u1=0, u=0, variance=0, maxVariance = 0;
        //目标像素占总图像比例w0
        //背景像素占总图像比例w1 
        //目标平均灰度值u0 
        //背景平均灰度值u1
        //总平均灰度值u
        //类间方差 variance
        //最大类间方差 maxVariance

    //计算每个像素在整幅图像中的比例  
    for (i = 0; i < GrayScale; i++)
    {
        pixelPro[i] = (float)pixelCount[i] / sumPixel;
         u += i * pixelPro[i];  //总平均灰度
    }
 
   
    for (i = 0; i < GrayScale; i++)     // i作为阈值 阈值从1-255遍历 
    {

        for (j = 0; j < GrayScale; j++)    //求目标前景和背景
        {
            if (j <= i)   //前景部分    
            {
                w0 += pixelPro[j];   
                u0Sum += j * pixelPro[j];
            }
            else   //背景部分  
            {
                w1 += pixelPro[j];
                u1Sum += j * pixelPro[j];
            }
        }
        u0 = u0Sum / w0;
        u1 = u1Sum / w1;

        variance = w0 * pow((u0 - u), 2) + w1 * pow((u1 - u), 2);  //类间方差计算公式

        if (variance > maxVariance)   //判断是否为最大类间方差
        {
            maxVariance = variance;
            threshold = i;
        }
    }

    return threshold;
    
}



image_threshold = otsuThreshold(image[0],COL,ROW);  //大津法计算阈值   

这就是OTUS的代码实现部分,当然上面的代码使用到了两个for循环判断,代码运行时间较长,实际上我们可以根据

w0+w1=1 (公式1) 概率和为1


(公式4)

公式1 和公式4 对代码来进行简化,将两个for循环进行叠加,最终简化代码如下,在STM32F4系列下实测用时为0.8ms左右



//------------------摄像头参数--------------//

uint8 image_threshold = 46;  //图像阈值
uint8 dis_image[60][80];

大津法二值化//

/*! 
 *  @brief      大津法二值化0.8ms程序
 *  @date:   2018-10  
 *  @since      v1.2
 *  *image :图像地址
 *  width:  图像宽
 *  height:图像高
 *  @author     Z小旋
 */

  
uint8 otsuThreshold(uint8 *image, uint16 width, uint16 height)
{
    #define GrayScale 256
    int pixelCount[GrayScale] = {0};//每个灰度值所占像素个数
    float pixelPro[GrayScale] = {0};//每个灰度值所占总像素比例
    int i,j;   
    int Sumpix = width * height;   //总像素点
    uint8 threshold = 0;
    uint8* data = image;  //指向像素数据的指针


    //统计灰度级中每个像素在整幅图像中的个数  
    for (i = 0; i < height; i++)
    {
        for (j = 0; j < width; j++)
        {
            pixelCount[(int)data[i * width + j]]++;  //将像素值作为计数数组的下标
          //   pixelCount[(int)image[i][j]]++;    若不用指针用这个
        }
    }
    float u = 0;  
    for (i = 0; i < GrayScale; i++)
    {
        pixelPro[i] = (float)pixelCount[i] / Sumpix;   //计算每个像素在整幅图像中的比例  
        u += i * pixelPro[i];  //总平均灰度
    }

      
    float maxVariance=0.0;  //最大类间方差
    float w0 = 0, avgValue  = 0;  //w0 前景比例 ,avgValue 前景平均灰度
    for(int i = 0; i < 256; i++)     //每一次循环都是一次完整类间方差计算 (两个for叠加为1个)
    {  
        w0 += pixelPro[i];  //假设当前灰度i为阈值, 0~i 灰度像素所占整幅图像的比例即前景比例
        avgValue  += i * pixelPro[i];        
        
        float variance = pow((avgValue/w0 - u), 2) * w0 /(1 - w0);    //类间方差   
        if(variance > maxVariance) 
        {  
            maxVariance = variance;  
            threshold = i;  
        }  
    } 


    return threshold;
    
}

使用:

 threshold = otsuThreshold(image[0],COL,ROW);  //大津法计算阈值  
 threshold = threshold > 50?50:threshold;


以上是关于详解-OTUS(大津法-最大类间方差)原理及C语言代码实现的主要内容,如果未能解决你的问题,请参考以下文章

大津法---OTSU算法

大津法(最大类间阈值法)

网易笔试题——计算机视觉_深度学习方向

Matlab图像处理问题---大津法。

图像分割最大类间方差法(otsu)图像分割

大津阈值法(OTSU)功能实现