详解-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]两部分。然后，计算该最佳阈值来区分这两类像素，使得两类像素区分度最大。
大津法阈值采用最大类间方差的原理，适合于图像灰度分布整体呈现“双峰”的情况。大津法会自动找出一个阈值，使得分割后的两部分类间方差最大。
特性：
-
大津法对噪音十分敏感，在处理之前应对图片进行去噪处理。
如果图像有存在局部噪声，则会影响大津法的判断
-
当目标与背景的面积比例悬殊的时候，类间方差函数可能呈现双峰或者多峰，这个时候 大津法的效果不好
双峰图像，目标和背景面积差距不大，可以很好的判断
当图像中的目标与背景的面积相差很大时，灰度直方图没有明显的双峰，或者两个峰的大小相差很大，分割效果不佳
目标像素点远小于背景像素点，分割效果不佳
代码实现：
基本所有的博客都是讲完原理，然后贴一整个代码，这样很多读者可能搞懂了原理，但是在阅读代码的时候突然又看不懂了，所以代码部分会按照上面的公式，推导编写，下面以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语言代码实现的主要内容，如果未能解决你的问题，请参考以下文章