论文《气体泄漏被动式红外成像检测理论及方法研究》部分代码实现

论文名称：李家琨. 气体泄漏被动式红外成像检测理论及方法研究[D].

论文地址：http://cdmd.cnki.com.cn/Article/CDMD-10007-1015801355.htm

本文主要实现论文的4.2章节：基于各向异性扩散的气体泄漏红外图像增强算法

代码实现

废话不多说，直接上代码，基于opencv 4.3.0。

#include <opencv2/opencv.hpp>
#include <iostream>

using namespace std;
using namespace cv;

string	origin_name		= "原图像";
string	aniso_name		= "各向异性扩散图像";
string	equal_name		= "直方图均衡化后的图像";
string	integral_name		= "积分图像";
string	k_means_name		= "kmeans分割图像";
string	result_name		= "结果图像";
int		frame_num		= 0;	//输入帧的数量

//********各向异性扩散滤波参数********
Mat AnisoDiff(const Mat& src, int iter, int k, float lambda);
float   lambda			= 0.23;	//扩散速度控制，最大取0.25
int     kappa			= 10;	//热传导系数，取值范围待经验确认
int		iteration		= 1;	//迭代次数

//********帧间差分及积分参数********
Mat FrameDiff(const deque<Mat>& src, int difference_interval, int integral_interval);
Mat ThreeFrameDiff(const deque<Mat>& src, int difference_interval, int integral_interval);
int		difference_interval = 2;	//差分帧间间隔
int		integral_interval	= 5;	//积分帧间间隔
int		min_diff_threshold	= 60;

//********k-means图像分割参数********
Mat KMeans(const Mat& src, int k, TermCriteria criteria);
int		k				= 3;	//聚类的种类
TermCriteria criteria	= TermCriteria(TermCriteria::EPS | TermCriteria::COUNT, 10, 1);	//迭代终止条件

int main()
{
	//VideoCapture cap("C:/Users/admin/Desktop/论文/test.avi");
	VideoCapture cap(0);
	Mat frame, image;
	deque<Mat> image_deque;

	for (int i = 0; i < difference_interval + integral_interval; i++)
	{
		cap >> frame;
		frame_num++;
		if (frame.empty())
			break;
		imshow(origin_name, frame);

		frame.copyTo(image);
		cvtColor(image, image, COLOR_BGR2GRAY);
		//1.直方图均衡化
		equalizeHist(image, image);
		//2.各向异性扩散滤波
		image = AnisoDiff(image, iteration, kappa, lambda);
		image_deque.push_back(image);
	}
	
	for (;;)
	{
		Mat img_aniso = image_deque.at(difference_interval);
		imshow(aniso_name, img_aniso);

		//3.帧间差分和积分
		//Mat img_integral = FrameDiff(image_deque, difference_interval, integral_interval);
		Mat img_integral = ThreeFrameDiff(image_deque, difference_interval, integral_interval);
		threshold(img_integral, img_integral, min_diff_threshold, 255.0, cv::THRESH_BINARY);
		dilate(img_integral, img_integral, Mat());		//膨胀			
		erode(img_integral, img_integral, Mat());		//腐蚀
		imshow(integral_name, img_integral);
		createTrackbar("threshold", integral_name, &min_diff_threshold, 255, 0);	//创建过滤最小阈值的调节条

		4.k-means
		//Mat img_kmeans = KMeans(img_integral, k, criteria);
		//imshow(k_means_name, img_kmeans);
		
		//5.彩色渲染和融合
		Mat img_result = frame.clone();
		for(int i = 0; i < img_result.rows; i++)
			for (int j = 0; j < img_result.cols; j++)
			{
				if (img_integral.at<uchar>(i, j) != 0)
					img_result.at<Vec3b>(i, j) = Vec3b(64, 128, 255);
			}
		imshow(result_name, img_result);
		//6.读取下一帧，再次进行第1步、第2步
		cap >> frame;
		if (frame.empty())
			break;
		imshow(origin_name, frame);

		frame.copyTo(image);
		cvtColor(image, image, COLOR_BGR2GRAY);
		cout << "当前帧："<<frame_num++ << endl;
		image_deque.pop_front();
		equalizeHist(image, image);
		image = AnisoDiff(image, iteration, kappa, lambda);
		image_deque.push_back(image);

		char key = (char)waitKey(10);
		if (key == 27 || key == 'q' || key == 'Q') // 'ESC'
			break;
	}

	return 0;
}

//********各向异性扩散滤波（灰度图像版本）（已完成）********
//参数说明：src--原图像，iter--迭代次数，k--kappa（传导系数），lambda--扩散速度
Mat AnisoDiff(const Mat& src, int iter, int k, float lambda)
{
	float ei, si, wi, ni;	//东南西北4个梯度
	float ce, cs, cw, cn;	
	int rows = src.rows;
	int cols = src.cols;
	//Mat temp = src.clone();
	Mat temp = src;	
	Mat output = cv::Mat::zeros(rows, cols, CV_8UC1);

	if (temp.channels() != 1)
	{
		cvtColor(temp, temp, COLOR_BGR2GRAY);	//判断通道数，转换为灰度图像
		cout << "图像不是单通道灰度图，已转换为灰度图!" << endl;
	}
	//if(temp.type() != CV_8UC1)
	//	cvtColor(temp, temp, COLOR_BGR2GRAY);	//判断图像类型，转换为灰度图像

	for (int n = 0; n < iter; n++)
	{
		for(int i = 1; i < rows - 1; i++)
			for (int j = 1; j < cols - 1; j++)
			{
				uchar cur = temp.at<uchar>(i, j);
				ei = temp.at<uchar>(i - 1, j) - cur;
				si = temp.at<uchar>(i, j + 1) - cur;
				wi = temp.at<uchar>(i + 1, j) - cur;
				ni = temp.at<uchar>(i, j - 1) - cur;

				ce = exp(-ei * ei / (k * k));
				cs = exp(-si * si / (k * k));
				cw = exp(-wi * wi / (k * k));
				cn = exp(-ni * ni / (k * k));

				output.at<uchar>(i, j) = cur + lambda * (ce * ei + cs * si + cw * wi + cn * ni);
			}
		output.copyTo(temp);
	}
	return output;
}

//********帧间差分及积分（已完成）********
//参数说明：src--原图像，difference_interval--差分图像采样间隔，integral_interval--积分间隔
Mat FrameDiff(const deque<Mat>& src, int difference_interval, int integral_interval)
{
	int num_frame_stored = difference_interval + integral_interval;	//需要存储的帧的数量
	if (src.size() != num_frame_stored)
	{
		cout << "存储帧数量与需要帧数量不符，请检查程序代码！" << endl;
	}
	Mat difference_result = cv::Mat::zeros(src.at(0).rows, src.at(0).cols, src.at(0).type());	//差分结果
	Mat integral_result = cv::Mat::zeros(src.at(0).rows, src.at(0).cols, src.at(0).type());	//积分结果
	//Mat integral_result = cv::Mat::zeros(src.at(0).rows, src.at(0).cols, CV_8UC1);
	for (int r = 0; r < integral_interval; r++)
	{
		absdiff(src.at(r), src.at(r + difference_interval), difference_result);
		//去除微小变化（此步骤暂缺）
		integral_result = integral_result + difference_result;
	}
	
	return integral_result;
}

//********三帧差分及积分（已完成）********
//参数说明：src--原图像，difference_interval--差分图像采样间隔，integral_interval--积分间隔
Mat ThreeFrameDiff(const deque<Mat>& src, int difference_interval, int integral_interval)
{
	int num_frame_stored = difference_interval + integral_interval;	//需要存储的帧的数量
	if (src.size() != num_frame_stored)
	{
		cout << "存储帧数量与需要帧数量不符，请检查程序代码！" << endl;
	}
	Mat difference_result_1	= cv::Mat::zeros(src.at(0).rows, src.at(0).cols, src.at(0).type());	//差分结果1
	Mat difference_result_2 = cv::Mat::zeros(src.at(0).rows, src.at(0).cols, src.at(0).type());	//差分结果2
	Mat difference_result = cv::Mat::zeros(src.at(0).rows, src.at(0).cols, src.at(0).type());
	Mat integral_result		= cv::Mat::zeros(src.at(0).rows, src.at(0).cols, src.at(0).type());	//积分结果
	//Mat integral_result = cv::Mat::zeros(src.at(0).rows, src.at(0).cols, CV_8UC1);
	for (int r = difference_interval; r < integral_interval; r++)
	{
		int previous = r - difference_interval, future = r + difference_interval;
		absdiff(src.at(r), src.at(r - difference_interval), difference_result_1);
		absdiff(src.at(r), src.at(r + difference_interval), difference_result_2);
		bitwise_and(difference_result_1, difference_result_2, difference_result);
		integral_result = integral_result + difference_result;
	}

	return integral_result;
}

//********k-means聚类（已完成）********
//参数说明：src--原图像，k--聚类类别数量，criteria--迭代终止条件
Mat KMeans(const Mat& src, int k, TermCriteria criteria = TermCriteria(TermCriteria::EPS | TermCriteria::COUNT, 10, 1))
{

	Mat centers;	//centers：最终分割后每个聚类的中心位置
	Mat dst(src.size(), src.type());
	Mat src_array = src.clone();
	src_array.convertTo(src_array, CV_32FC1);
	src_array = src_array.reshape(0, src.rows * src.cols);	//转化为1维向量，通道数不变，reshape按行序列化（即按src的行取值，填入dst的行中）
	Mat_<int> labels(src_array.size(), CV_32SC1);	//labels：最终分类每个样本的标签
	
	去除大量0元素，提升速度（此步骤暂缺）
	//vector<float> src_vec;
	//for (int i = 0; i < src_array.rows; i++)
	//{
	//	if(src_array.at<float>(i, 0) != 0)
	//	{
	//		//cout << src_array.at<float>(i, 0) << endl;
	//		src_vec.push_back(src_array.at<float>(i, 0));
	//	}
	//}
	//src_array = Mat(src_vec, true);

	//kmeans
	double compactness = kmeans(src_array, k, labels, criteria, 1, cv::KMEANS_PP_CENTERS, centers);	
	//double compactness = kmeans(src_array, k, labels, criteria, 1, cv::KMEANS_RANDOM_CENTERS, centers);
	//将每一类的像素转化为聚类中心的像素值
	MatIterator_<uchar> itd = dst.begin<uchar>(), itd_end = dst.end<uchar>();	//这里还是人工设定了dst是灰度图，像素类型是uchar，不够“泛型”
	for (int i = 0; itd != itd_end; i++, itd++)
	{
		uchar color = centers.at<float>(labels.at<int>(i, 0), 0);
		(*itd) = saturate_cast<uchar>(color);
	}
	return dst;
}

解析

1. 各向异性扩散滤波

论文中提到的P-M方程实际上就是各向异性扩散滤波最基础的实现方式 ，原理代码实现参考了文章末尾的[1]。各向异性扩散可用于平滑图像，同时保留边缘信息。

2. 帧差法与积分

这部分非常好理解，差分就是当前帧与之前的某一帧做差，积分就是差分图像累计求和，论文中的公式为：
$$I_{out}^n=\sum _{r=0}^{q-1}(I^{n+r-p}-I^{n+r})$$
式中：$I^n$为当前帧，$p$ 表示差分间隔，$q$ 表示积分间隔。

根据公式，在程序实现中，在当前帧之前，需要存储的帧数量为 $q$，在当前帧之后（包含当前帧），需要存储的数量为 $q$，那么进行差分和积分操作，一共需要存储的帧数量为 $q+p$ ，之后更新过程就是删除存储帧的首元素，再在尾部插入新的帧，因为操作都在首尾，所以使用deque容器。

本文在实际应用中使用帧差法对背景变化非常敏感，因此选用了三帧差法，主要差别在：

	bitwise_and(difference_result_1, difference_result_2, difference_result);

两次差分结果按位进行“与”操作，即前后都产生变化的像素才会被计入最终的差分结果。

3. k-means聚类

k-means算法原理论文中有详细描述，网络上也有许多科普文章，本文不再赘述，此处说明一下opencv库自带的kmeans函数。

double cv::kmeans( 
	InputArray data, 
	int K, 
	InputOutputArray bestLabels,
	TermCriteria citeria, 
	int attempts, 
	int flags, 
	OutputArray centers
)

（1）data是需要聚类的样本集合，k-means按行组织样本，每一行是一个样本，列数只有1列。需要注意的是，kmeans函数只接受CV_32F的数据，即浮点数，灰度图类型为CV_8U，需要进行转换；

（2）K是聚类的类别数量；

（3）bestLabels是各个样本的类别标记，整型数字，同样只有1列。举例，样本被分为6类，则各个类别标记可表示为[0, 1, 2, 3, 4, 5]，如果data中第2行样本被分类至第[3]类，则bestLabels对应的在第二行存储[3]；

（4）citeria为迭代终止条件，主要包括迭代次数达到上限和误差满足精度要求两种，搜索TermCriteria类有详细的解释；

（5）attempts是进行k-means算法的次数，选取结果最佳的一次作为最终结果；

（6）flags是聚类中心初始化的方式，KMEANS_RANDOM_CENTERS 表示随机初始化聚类心KMEANS_PP_CENTERS 表示用kmeans++算法来初始化聚类中心，KMEANS_USE_INITIAL_LABELS 表示使用用户自定义值初始化聚类中心；

（7）centers用于存放聚类中心。

如何将图像像素转换为kmeans需要的数据？

灰度图像像素具有强度和位置两种属性，k-means用于图像分割时，将像素的强度作为样本值，利用reshape()函数将图像“碾平”为1维单列向量：

	src_array = src_array.reshape(0, src.rows * src.cols);

函数原型：

	C++: Mat Mat::reshape(int cn, int rows=0) const

（1）cn表示新Mat的通道数，如果为0则表示保持通道数不变；

（2）rows表示新Mat的行数，如果为0则表示保持行数不变。

reshape()函数可以看作将源图像的像素按行依次提取，再按行依次存储进reshape()设定的新图像中，前后像素总数相等。

kmeans完成之后？

获取聚类中心centers和类别标记bestLabels后，就可以根据bestLabels，将同一类的像素统一为相同的值了，这样就完成了图像分割的效果，相同值可以是聚类中心的值，也可以是自己设定的值。

问题与思考

1. 运行时间

整个程序的运行时间是无法达到实时检测的需求的，有大量时间花费于k-means步骤。k-means算法需要进行多次迭代，每次迭代都需要全部的像素与聚类中心计算欧式距离，而且opencv的kmeans()仅接受浮点数运算，计算量还是比较大的。

根据论文，k-means的目的在于去除微弱的背景信息，我自己的理解是：因为摄像头的拍摄质量，在视频中，除了漏气区域，背景也会有变化，因此在差分图像中也会存在背景变化的信息。进行3类聚类时，由于背景变化较为微弱，属于背景变化的非0像素值基本较小，会与0值像素分为同一类，统一同类像素值时这些微弱变化的像素可以被“抹去”，前景和背景就被分离开了。

但从实际应用上考虑，kmeans()的性价比确实有些低，本文尝试使用一种更简单的方法。获取积分图像后，进行二值化处理，并设定最小阈值，低于此阈值的像素值会被设置为0，本文同时设置了一个最低阈值的调节条：

	createTrackbar("threshold", integral_name, &min_diff_threshold, 255, 0)

通过滤除较小的像素值，也达成了一定的过滤微弱变化的背景的功能，速度可以达到实时了，但是不能实现论文中的两种气体浓度分割的效果了。

2. opencv函数Mat元素类型

刚学习opencv时，确实不太注意图像的类型，但其实还是需要注意的，比如代码中的直方图均衡函数equalizeHist()只接受CV_8UC1类型的图像，即8位1维图像，多维图像需要分通道处理；而kmeans()函数又只接受CV_32F类型的图像，两者结合使用时，如果不进行图像类型转换的话，可能就会出现如下的错误：

	Error: Assertion failed (data0.dims <= 2 && type == CV_32F && K > 0)

参考

[1]https://blog.csdn.net/bluecol/article/details/46690985

[2]https://blog.csdn.net/little_white__/article/details/88878567

[3]https://github.com/yuzsh/kmeans_segmentation