详解-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]两部分。然后，计算该最佳阈值来区分这两类像素，使得两类像素区分度最大。

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

特性：

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

如果图像有存在局部噪声，则会影响大津法的判断
当目标与背景的面积比例悬殊的时候，类间方差函数可能呈现双峰或者多峰，这个时候大津法的效果不好

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

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

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

代码实现：

基本所有的博客都是讲完原理，然后贴一整个代码，这样很多读者可能搞懂了原理，但是在阅读代码的时候突然又看不懂了，所以代码部分会按照上面的公式，推导编写，下面以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 256int 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;