【OpenCV学习笔记】4.阈值分割

阈值分割是一种基于区域的、简单的通过灰度信息提取形状的技术，实现简单、计算量小、性能稳定，因此应用广泛。阈值分割后的输出图像只有两种灰度值：255和0，因此阈值分割处理又被称为图像二值化。

1.方法概述

全局阈值分割

全局阈值分割是将灰度值大于thresh（阈值）的像素设为白色，小于或者等于thresh的像素设为黑色；或者反过来。

假设输入图像为 I I I，高为 H H H，宽为 W W W， I ( r , c ) I(r,c) I(r,c)代表 I I I的第 r r r行第 c c c列的灰度值， 0 ≤ r < H , 0 ≤ c < W , 0 \leq r<H,0\leq c<W, 0≤r<H,0≤c<W,全局阈值处理后的输出图像为 O O O，则：

O ( r , c ) = { 255 , I ( r , c ) > t h r e s h 0 , I ( r , c ) ≤ t h r e s h O(r,c)= \begin{cases}255,&I(r,c)>thresh\\ 0,&I(r,c)\leq thresh \end{cases} O(r,c)={255,0,I(r,c)>threshI(r,c)≤thresh

或者

O ( r , c ) = { 0 , I ( r , c ) > t h r e s h 255 , I ( r , c ) ≤ t h r e s h O(r,c)= \begin{cases}0,&I(r,c)>thresh\\ 255,&I(r,c)\leq thresh \end{cases} O(r,c)={0,255,I(r,c)>threshI(r,c)≤thresh

threshold(InputArray src, OutputArray dst, double thresh, double maxval, int type)
//src-输入矩阵，数据类型为CV_8U或者CV_32F
//dst-输出矩阵
//thresh-阈值
//maxval-图像二值化时，一般为255
//type-类型

当type=THRESH_BINARY,

d s t ( r , c ) = { m a x v a l , s r c ( r , c ) > t h r e s h 0 , s r c ( r , c ) ≤ t h r e s h dst(r,c)= \begin{cases}maxval,&src(r,c)>thresh\\ 0,&src(r,c)\leq thresh \end{cases} dst(r,c)={maxval,0,src(r,c)>threshsrc(r,c)≤thresh

当type=THRESH_BINARY_INV,

d s t ( r , c ) = { 0 , s r c ( r , c ) > t h r e s h m a x v a l , s r c ( r , c ) ≤ t h r e s h dst(r,c)= \begin{cases}0,&src(r,c)>thresh\\ maxval,&src(r,c)\leq thresh \end{cases} dst(r,c)={0,maxval,src(r,c)>threshsrc(r,c)≤thresh

当type为THRESH_OTSU或THRESH_TRIANGLE时，src只支持uchar类型，此时thresh作为输出参数，即通过Otsu和TRIANGLE算法自动计算。

//examples：
threshold(src, dst, 150, 255, THRESH_BINARY);
triThe = threshold(src, dst_tri, 0, 255, THRESH_TRIANGLE+THRESH_BINARY);

手动设置阈值效果不是很好，直方图技术法、Otsu算法、熵算法是比较流行的自动选取全局阈值的算法。

局部阈值分割

在有些情况下，如光照不均，全局阈值分割效果不理想，使用局部阈值（自适应阈值）可以产生较好的效果。

O ( r , c ) = { 255 , I ( r , c ) > t h r e s h ( r , c ) 0 , I ( r , c ) ≤ t h r e s h ( r , c ) O(r,c)= \begin{cases}255,&I(r,c)>thresh(r,c)\\ 0,&I(r,c)\leq thresh(r,c) \end{cases} O(r,c)={255,0,I(r,c)>thresh(r,c)I(r,c)≤thresh(r,c)

或者

O ( r , c ) = { 0 , I ( r , c ) > t h r e s h ( r , c ) 255 , I ( r , c ) ≤ t h r e s h ( r , c ) O(r,c)= \begin{cases}0,&I(r,c)>thresh(r,c)\\ 255,&I(r,c)\leq thresh(r,c) \end{cases} O(r,c)={0,255,I(r,c)>thresh(r,c)I(r,c)≤thresh(r,c)

局部阈值分割的核心也是计算阈值矩阵，常用的方法是自适应阈值算法（移动平均值算法）。

2.直方图技术法

一幅含有一个与背景程明显对比的图像具有包含双峰的直方图，两个峰值对应于物体内部和外部较多数目的点，两个峰值之间的波谷对应于物体边缘附近相对较少数目的点。直方图技术法就是取两个峰值之间的波谷位置的灰度值为阈值。下面是一种自动选取波峰和波谷的方法。

【OpenCV学习笔记】4.阈值分割

假设输入图像为 I I I，高为 H H H，宽为 W W W， h i s t o g r a m I histogram_I histogramI代表其对应的灰度直方图， h i s t o g r a m I ( k ) histogram_I(k) histogramI(k)代表灰度值为 k k k的像素点个数，其中 0 ≤ k ≤ 255 0\leq k\leq 255 0≤k≤255。

第一步：找到灰度值第一个峰值，即灰度直方图最大值，并找到其对应灰度值，用 f i r s t P e a k firstPeak firstPeak表示。

第二步：找到灰度值第二个峰值，并找到其对应灰度值。第二个峰值不一定是直方图第二大值，因为它可能出现在第一个峰值附近。可通过公式计算：

s e c o n d P e a k = a r g k m a x { ( k − f i r s t P e a k ) 2 ∗ h i s t o g r a m I ( k ) } secondPeak=arg_kmax\{(k-firstPeak)^2*histogram_I(k)\} secondPeak=argkmax{(k−firstPeak)2∗histogramI(k)}

或 s e c o n d P e a k = a r g k m a x { ∣ k − f i r s t P e a k ∣ ∗ h i s t o g r a m I ( k ) } secondPeak=arg_kmax\{|k-firstPeak|*histogram_I(k)\} secondPeak=argkmax{∣k−firstPeak∣∗histogramI(k)}

第三步：找到这两个峰值之间的波谷，若出现两个或多个波谷，取左侧波谷，其对应灰度值即为阈值。

int threshTwoPeaks(const Mat & image, Mat & thresh_out)
{
	//计算灰度直方图
	Mat histogram = calcGrayHist(image);//函数调用
	//找到灰度直方图最大峰值对应的灰度值
	Point firstPeakLoc;
	minMaxLoc(histogram, NULL, NULL, NULL, &firstPeakLoc);
	int firstPeak = firstPeakLoc.x;
	//寻找灰度直方图第二个峰值对应的灰度值
	Mat measureDists = Mat::zeros(Size(256, 1), CV_32FC1);
	for (int k = 0; k < 256; k++)
	{
		int hist_k = histogram.at<int>(0, k);
		measureDists.at<float>(0, k) = pow(float(k - firstPeak), 2)*hist_k;
	}
	Point secondPeakLoc;
	minMaxLoc(measureDists, NULL, NULL, NULL, &secondPeakLoc);
	int secondPeak = secondPeakLoc.x;
	//找到两个峰值之间最小值对应的灰度值，作为阈值
	Point threshLoc;
	int thresh = 0;
	if (firstPeak < secondPeak)//第一个峰值在第二个峰值左侧
	{
		minMaxLoc(histogram.colRange(firstPeak, secondPeak), NULL, NULL, &threshLoc);
		thresh = firstPeak + threshLoc.x + 1;
	}
	else//第一个峰值在第二个峰值右侧
	{
		minMaxLoc(histogram.colRange(secondPeak, firstPeak), NULL, NULL, &threshLoc);
		thresh = secondPeak + threshLoc.x + 1;
	}
	//阈值分割
	threshold(image, thresh_out, thresh, 255, THRESH_BINARY);
	return thresh;
}

Mat calcGrayHist(const Mat & image)
{
	Mat histogram = Mat::zeros(Size(256, 1), CV_32SC1);//存储256个灰度级的像素个数
	int rows = image.rows;
	int cols = image.cols;//图像的宽和高
	//计算灰度级个数
	for (int r = 0; r < rows; r++)
	{
		for (int c = 0; c < cols; c++)
		{
			int index = int(image.at<uchar>(r, c));
			histogram.at<int>(0, index) += 1;
		}
	}
	return histogram;
}

3.熵算法

信息熵的概念来源于信息论，假设信源符号 u u u有 N N N种取值，记为

u 1 , u 2 , … … , u N u_1,u_2,……,u_N u1,u2,……,uN

且每一种信源符号出现的概率，记为

p 1 , p 2 , … … , p N p_1,p_2,……,p_N p1,p2,……,pN

那么该信源符号的信息熵，记为

e n t r o p y ( u ) = − ∑ i = 1 N p i l o g p i entropy(u)=-\sum_{i=1}^{N}p_ilogp_i entropy(u)=−∑i=1Npilogpi
图像也可看做一种信源，假设输入图像为 I I I， n o r m H i s t I normHist_I normHistI代表归一化的灰度图像直方图，那么对于8位图可以看作由256个灰度符号，且每个符号出现的概率为 n o r m H i s t I ( k ) normHist_I(k) normHistI(k)组成的信源，其中 0 ≤ k ≤ 255 0\leq k\leq 255 0≤k≤255。
利用熵计算阈值步骤如下：

第一步：计算 I I I的累加概率直方图，又称零阶累积矩，记为

c u m u H i s t ( k ) = ∑ i = 0 k n o r m H i s t I ( i ) , k ∈ [ 0 , 255 ] cumuHist(k)=\sum_{i=0}^{k}normHist_I(i),k\in [0,255] cumuHist(k)=∑i=0knormHistI(i),k∈[0,255]

第二步：计算各个灰度级的熵，记为

e n t r o p y ( u ) = − ∑ k = 0 t n o r m H i s t I ( k ) l o g ( n o r m H i s t I ( k ) ) , t ∈ [ 0 , 255 ] entropy(u)=-\sum_{k=0}^{t}normHist_I(k)log(normHist_I(k)),t\in [0,255] entropy(u)=−∑k=0tnormHistI(k)log(normHistI(k)),t∈[0,255]

第三步：计算使 f ( t ) = f 1 ( t ) + f 2 ( t ) f(t)=f_1(t)+f_2(t) f(t)=f1(t)+f2(t)最大化的 t t t值，该值即为得到的阈值，即 t h r e s h = a r g t m a x ( f ( t ) ) thresh=arg_tmax(f(t)) thresh=argtmax(f(t))，其中

f 1 ( t ) = e n t r o p y ( t ) e n t r o p y ( 255 ) l o g ( c u m u H i s t ( t ) ) l o g ( m a x { c u m u H i s t ( 0 ) , c u m u H i s t ( 1 ) , … … , c u m u H i s t ( t ) } ) f_1(t)=\frac{entropy(t)}{entropy(255)}\frac{log(cumuHist(t))}{log(max\{cumuHist(0),cumuHist(1),……,cumuHist(t)\})} f1(t)=entropy(255)entropy(t)log(max{cumuHist(0),cumuHist(1),……,cumuHist(t)})log(cumuHist(t))

f 2 ( t ) = ( 1 − e n t r o p y ( t ) e n t r o p y ( 255 ) ) l o g ( 1 − c u m u H i s t ( t ) ) l o g ( m a x { c u m u H i s t ( t + 1 ) , … … , c u m u H i s t ( 255 ) } ) f_2(t)=(1-\frac{entropy(t)}{entropy(255)})\frac{log(1-cumuHist(t))}{log(max\{cumuHist(t+1),……,cumuHist(255)\})} f2(t)=(1−entropy(255)entropy(t))log(max{cumuHist(t+1),……,cumuHist(255)})log(1−cumuHist(t))

4.Otsu阈值处理

Otsu是一种最大方差法，该算法是在判别分析最小二乘法原理的基础上推导的。
假设输入图像为 I I I，高为 H H H，宽为 W W W， h i s t o g r a m I histogram_I histogramI代表其对应的灰度直方图， h i s t o g r a m I ( k ) histogram_I(k) histogramI(k)代表灰度值为 k k k的像素点个数，其中 0 ≤ k ≤ 255 0\leq k\leq 255 0≤k≤255。算法详细步骤如下：

第一步：计算灰度直方图的零阶累积矩

z e r o C u m u M o m e n t ( k ) = ∑ i = 0 k h i s t o g r a m I ( i ) , k ∈ [ 0 , 255 ] zeroCumuMoment(k)=\sum_{i=0}^{k}histogram_I(i),k\in [0,255] zeroCumuMoment(k)=∑i=0khistogramI(i),k∈[0,255]

第二步：计算灰度直方图的一阶累积矩

o n e C u m u M o m e n t ( k ) = ∑ i = 0 k ( i ∗ h i s t o g r a m I ( i ) ) , k ∈ [ 0 , 255 ] oneCumuMoment(k)=\sum_{i=0}^{k}(i*histogram_I(i)),k\in [0,255] oneCumuMoment(k)=∑i=0k(i∗histogramI(i)),k∈[0,255]

第三步：计算图像 I I I总体的灰度平均值 m e a n mean mean，其实就是 k = 255 k=255 k=255时的一阶累积矩，即

m e a n = o n e C u m u M o m e n t ( 255 ) mean=oneCumuMoment(255) mean=oneCumuMoment(255)

第四步：计算每一个灰度级作为阈值时，前景区域的平均灰度、背景区域的平均灰度与整幅图像的平均灰度的方差。对方差的衡量采用以下度量：

σ 2 ( k ) = ( m e a n ∗ z e r o C u m u M o m e n t ( k ) − o n e C u m u M o m e n t ( k ) ) 2 z e r o C u m u M o m e n t ( k ) ∗ ( 1 − z e r o C u m u M o m e n t ( k ) ) , k ∈ [ 0 , 255 ] \sigma^2(k)=\frac{(mean*zeroCumuMoment(k)-oneCumuMoment(k))^2}{zeroCumuMoment(k)*(1-zeroCumuMoment(k))},k\in [0,255] σ2(k)=zeroCumuMoment(k)∗(1−zeroCumuMoment(k))(mean∗zeroCumuMoment(k)−oneCumuMoment(k))2,k∈[0,255]

第五步：找到上述最大的 σ 2 ( k ) \sigma^2(k) σ2(k)，对应的 k k k即为Otsu自动选取的阈值，即

t h r e s h = a r g k ∈ [ 0 , 255 ] m a x ( σ 2 ( k ) ) thresh=arg_{k\in [0,255]}max(\sigma^2(k)) thresh=argk∈[0,255]max(σ2(k))

int otsu(const Mat & image, Mat & OtsuThreshImage)
{
	//计算灰度直方图
	Mat histogram = calcGrayHist(image);//函数调用
	//归一化灰度直方图
	Mat normHist;
	histogram.convertTo(normHist, CV_32FC1, 1.0 / (image.rows*image.cols), 0.0);
	//计算零阶累积矩和一阶累积矩
	Mat zeroCumuMoment = Mat::zeros(Size(256, 1), CV_32FC1);
	Mat oneCumuMoment = Mat::zeros(Size(256, 1), CV_32FC1);
	for (int i = 0; i < 256; i++)
	{
		if (i == 0)
		{
			zeroCumuMoment.at<float>(0, i) = normHist.at<float>(0, i);
			oneCumuMoment.at<float>(0, i) = i * normHist.at<float>(0, i);
		}
		else
		{
			zeroCumuMoment.at<float>(0, i) = zeroCumuMoment.at<float>(0, i - 1) + normHist.at<float>(0, i);
			oneCumuMoment.at<float>(0, i) = oneCumuMoment.at<float>(0, i - 1) + i * normHist.at<float>(0, i);
		}
	}
	//计算类方间差
	Mat variance = Mat::zeros(Size(256, 1), CV_32FC1);
	//总平均值
	float mean = oneCumuMoment.at<float>(0, 255);
	for (int i = 0; i < 255; i++)
	{
		if (zeroCumuMoment.at<float>(0, i) == 0 || zeroCumuMoment.at<float>(0, i) == 1)
			variance.at<float>(0, i) = 0;
		else
		{
			float cofficient = zeroCumuMoment.at<float>(0, i)*(1.0 - zeroCumuMoment.at<float>(0, i));
			variance.at<float>(0, i) = pow(mean*zeroCumuMoment.at<float>(0, i) - oneCumuMoment.at<float>(0, i), 2.0) / cofficient;
		}
	}
	//找到阈值
	Point maxLoc;
	minMaxLoc(variance, NULL, NULL, NULL, &maxLoc);
	int otsuThresh = maxLoc.x;
	//阈值处理
	threshold(image, OtsuThreshImage, otsuThresh, 255, THRESH_BINARY);
	return otsuThresh;
}

Mat calcGrayHist(const Mat & image)
{
	Mat histogram = Mat::zeros(Size(256, 1), CV_32SC1);//存储256个灰度级的像素个数
	int rows = image.rows;
	int cols = image.cols;//图像的宽和高
	//计算灰度级个数
	for (int r = 0; r < rows; r++)
	{
		for (int c = 0; c < cols; c++)
		{
			int index = int(image.at<uchar>(r, c));
			histogram.at<int>(0, index) += 1;
		}
	}
	return histogram;
}

5.自适应阈值

在图像进行平滑处理时，均值平滑、高斯平滑、中值平滑用不同规则计算出以当前像素为中心的邻域内灰度的平均值，所以可以使用平滑处理后的输出结果作为每个像素设置阈值的参考值。
在自适应阈值处理中，平滑算子的尺寸决定了分割出来的物体的尺寸，如果滤波器尺寸太小，那么估计出的局部阈值将不理想。一般平滑算子的宽度必须大于被识别物体的宽度，平滑算子尺寸越大，平滑后的结果越能更好作为每个像素阈值的参考，当然不能无限大。
假设输入图像为 I I I，高为 H H H，宽为 W W W，平滑算子尺寸记为 H × W H\times W H×W，其中 W W W和 H H H均为奇数，自适应阈值分割步骤如下：

第一步：对图像进行平滑处理，平滑结果记为 f s m o o t h ( I ) f_{smooth}(I) fsmooth(I)，可以是均值平滑、高斯平滑、中值平滑。

第二步：自适应阈值矩阵 T h r e s h = ( 1 − r a t i o ) ∗ f s m o o t h ( I ) Thresh=(1-ratio)*f_{smooth}(I) Thresh=(1−ratio)∗fsmooth(I),一般 r a t i o = 0.15 ratio=0.15 ratio=0.15。

第三步：利用局部阈值分割规则进行阈值分割。

O ( r , c ) = { 255 , I ( r , c ) > t h r e s h ( r , c ) 0 , I ( r , c ) ≤ t h r e s h ( r , c ) O(r,c)= \begin{cases}255,&I(r,c)>thresh(r,c)\\ 0,&I(r,c)\leq thresh(r,c) \end{cases} O(r,c)={255,0,I(r,c)>thresh(r,c)I(r,c)≤thresh(r,c)

或者

O ( r , c ) = { 0 , I ( r , c ) > t h r e s h ( r , c ) 255 , I ( r , c ) ≤ t h r e s h ( r , c ) O(r,c)= \begin{cases}0,&I(r,c)>thresh(r,c)\\ 255,&I(r,c)\leq thresh(r,c) \end{cases} O(r,c)={0,255,I(r,c)>thresh(r,c)I(r,c)≤thresh(r,c)

void adaptiveThreshold(InputArray src, OutputArray dst, double maxValue,
         int adaptiveMethod, int thresholdType, int blockSize, double C)
//src-输入矩阵，数据类型为CV_8U
//dst-输出矩阵
//maxValue-一般为255
//adaptiveMethod-ADAPTIVE_THRESH_MEAN_C采用均值平滑，
//               ADAPTIVE_THRESH_GAUSSIAN_C采用高斯平滑
//thresholdType-THRESH_BINARY,THRESH_BINARY_INV
//blockSize-平滑算子尺寸，奇数
//C-比例系数

【OpenCV学习笔记】4.阈值分割

1.方法概述

2.直方图技术法

3.熵算法

4.Otsu阈值处理

5.自适应阈值

继续阅读

C语言第四章自述2第四章选择结构程序设计

面试题:vector和map的区别，异同。空间分布，100万数据存哪个比较合适。一、迭代器区别二、vector三、Map、Set四、vector_map 为什么比map效率高五、如何选择六、容器选择原则七、效率对比

C++ 多线程用条件变量确定线程的执行顺序而不是使用 sleep(1)

POJ 1284 Primitive Roots (欧拉函数&原根定理)

CQ V1.0分词bates(基于双数组tire树)—应该是目前最快的中文分词算法

成员函数初始化列表

2021-08-13c++——类之操作符重载

swmm与lisflood-fp源码如何一起编译 CMake命令

Windows下VS开发环境环境安装工程项目设置关于Debug和Release的提示

一文看懂字符串的加减乘除

C++ 第十五周报告1--《冒泡法排序》

C++实现简单顺序表

C经典书籍笔记——C陷阱与缺陷②(语法陷阱之优先级)一、错误案列二、优先级规律

线性表之顺序表的实现

C++判断素数、求最大公约数代码判断一个数是否为素数求两个数的最大公约数

SequoiaDB巨杉数据库C++驱动概述