C++数字图像处理（4）—均值滤波_c++分别采用3*3和5*5的窗口对噪声污染图像进行滤波

1、算法原理

在数字图像处理中，滤波是一个很重要的操作，许多算法其本质都是滤波操作。使用白话的形式对滤波定义：对于一个像素点，使用其领域像素（可以包含自身，也可不包含）的相关特性，计算出一个值，代替当前像素值。举个例子，3X3均值滤波，就是计算每个像素对应的3X3领域所有像素值的平均值，代替当前像素。滤波操作可分为线性滤波和非线性滤波两大种类，其中，常见的线性滤波操作有：均值滤波、高斯滤波、方差滤波（局部方差）等等，常见的非线性滤波操作有：中值滤波、最大值滤波、最小值滤波等等。

2、算法实现

（1）、边界处理

一个完善的滤波算法，边界处理是必要的。对于图像的边界像素，其领域是不存在的，比如坐标为（0,0）（y，x形式）的像素点，其3X3的领域坐标集合为

其中，所有坐标包含-1的像素点都是不存在的，这个时候滤波操作就需要对这些不存在的点进行处理。有几种边界处理方式：固定值填充法、镜像法、重叠法。所谓的固定值填充法即为：所有不存在的像素点的像素值认为是某一个固定的值（一般由算法调用者去设定，或者默认为0）；镜像法：以边界像素为对称轴，使用对称像素点的像素值代替；重叠法：使用边界像素值代替。示意图如下：

一般来说，如无特殊的应用场景，推荐大家使用边界重叠法进行边界处理。C++代码实现边界重叠法代码如下：

//函数名：makeRepeatBorder
//作用：边界填充（边界像素重叠法）
//参数：
//matInput：输入图像
//matOutput : 输出图像
//cnBorder : 边界尺寸
//返回值：无
//注：支持单通道8位灰度图像
void makeRepeatBorder(cv::Mat& matInput, cv::Mat& matOutput, const int& cnBorder)
{
	//构造输出图像
	matOutput = cv::Mat::zeros(matInput.rows + cnBorder * 2, matInput.cols + cnBorder * 2, matInput.type());
 
	unsigned char *pSrc = NULL, *pDst = NULL;
 
	//拷贝原始图像数据
	for (int m = 0; m < matInput.rows; m++)
	{
		//源地址
		pSrc = matInput.data + m * matInput.step[0];
		//目的地址
		pDst = matOutput.data + (m + cnBorder) * matOutput.step[0] + cnBorder;
		memcpy(pDst, pSrc, matInput.step[0]);
	}
	//边界处理
	for (int m = 0; m < cnBorder; m++)
	{
		//顶部
		pSrc = matOutput.data + cnBorder * matOutput.step[0];
		pDst = matOutput.data + m * matOutput.step[0];
		memcpy(pDst, pSrc, matOutput.step[0]);
		//底部
		pSrc = matOutput.data + (cnBorder + matInput.rows - 1) * matOutput.step[0];
		pDst = matOutput.data + (cnBorder + matInput.rows + m) * matOutput.step[0];
		memcpy(pDst, pSrc, matOutput.step[0]);
	}
 
	for (int m = 0; m < matOutput.rows; m++)
		for (int n = 0; n < cnBorder; n++)
		{
			//左
			matOutput.at<unsigned char>(m, n) = matOutput.at<unsigned char>(m, cnBorder);
			//右
			matOutput.at<unsigned char>(m, n + matInput.cols + cnBorder) = matOutput.at<unsigned char>(m, matOutput.cols - cnBorder - 1);
		}
}

效果如下：

其余两种方法，大家可以自行编码实现。

（2）、均值滤波

经过边界处理后，图像尺寸增大了。假如进行5*5的均值滤波，图像的尺寸增大了4，那么，在边界处理后的图像进行滤波操作的时候，滤波的起始、停止下标同样需要进行改变。一般我们滤波的窗口尺寸设为奇数，主要是考虑到窗口存在滤波中心（锚点），便于计算（个人理解，未经证实）。C++实现代码如下：

//函数名：meanFilter
//作用：均值滤波实现
//参数：
//matInput：输入图像
//matOutput : 输出图像
//cnWinSize : 滤波窗口尺寸
//返回值：无
//注：支持单通道8位灰度图像
void meanFilter(cv::Mat& matInput, cv::Mat& matOutput, const int& cnWinSize)
{
	//滤波窗口为奇数
	assert(cnWinSize % 2 == 1);
 
	//构造输出图像
	matOutput = cv::Mat::zeros(matInput.rows, matInput.cols, matInput.type());
 
	int nAnchor = cnWinSize / 2;
	//边界处理
	cv::Mat matMakeBorder;
	makeRepeatBorder(matInput, matMakeBorder, nAnchor);
 
	//使用行首地址存储法遍历滤波
	unsigned char ** pRow = new unsigned char*[matMakeBorder.rows];
	for (int r = 0; r < matMakeBorder.rows; r++)
		pRow[r] = matMakeBorder.data + r * matMakeBorder.step[0];
	
	//乘法代替除法
	float fNormal = 1.0f / (cnWinSize * cnWinSize);
 
	//滤波操作
	for (int r = 0; r < matMakeBorder.rows - cnWinSize + 1; r++)
	{
		unsigned char* pOut = matOutput.data + r * matOutput.step[0];
		for (int c = 0; c < matMakeBorder.cols - cnWinSize + 1; c++)
		{
			//求和
			float fSum = 0;
			for (int m = 0; m < cnWinSize; m++)
			for (int n = 0; n < cnWinSize; n++)
			{
				fSum += pRow[m + r][n + c];
			}
 
			//求均值输出
			pOut[c] = (unsigned char)(fSum * fNormal);
		}
	}
	delete[] pRow;
}

滤波结果下：

3、均值滤波的加速

经过上面的编码实现，我们可以粗略的计算下，该种方式实现均值滤波的算法复杂度。假设滤波图像尺寸为W*H，滤波窗口为M*N，不考虑边界处理的情况下，内存访问的次数为W*H*M*N。也就是说，算法的耗时与滤波窗口的尺寸相关。算法耗时情况如下表：

那么，是否有更加优雅的方式实现呢？我这里提供一种基于积分图加速办法（查看OpenCV源码可知，cpu处理下，还存在行列滤波加速、频域加速等的方法）：

（1）、积分图原理

积分图的首次出现在基于HAAR特征的人脸检测的框架中，原始论文为《Rapid object detection using a boosted cascade of simple features》。Viola 提出了一种利用积分图(integral image)快速计算 Haar 特征的方法, 这个方法使得图像的局部矩形求和运算的复杂度从 O(mn) 下降到了 O(4) 。图像是由一系列的离散像素点组成，因此图像的积分其实就是求和。 图像积分图中每个点的值是原图像中该点左上角的所有像素值之和。示意图如下：

 

在积分图像中，I(x,y)存储的是SAT(x，y)矩形区域对应原始图像的像素和值。那么，在给定积分图后，计算任意矩形区域(x,y,w,h)的像素的和值的公式为：

也就是说，得到积分图后，仅仅需要四个点就可以计算任意矩形区域的像素值之和。

（2）、积分图快速实现

积分图的C++实现有许多种方法，网络上也有很多博主写出了。我查看不少的博文，我向大家推荐这位博主的原理讲解https://blog.csdn.net/xiaowei_cqu/article/details/17928733。其实，大家去查看OpenCV的积分图实现源码就可以发现，积分图的尺寸是比原始图像的尺寸大1的。

这是为什么呢？某个点的积分图反应的是原图中此位置左上角所有像素之和，这里是的累加和是不包括这个点像素本身的，那么这样，对于原图的第一行和第一列的所有像素，其对应位置的积分图就应该是0，这样考虑到所有的像素，为了能容纳最后一列和最后一行的情况，最终的积分图就应该是 (W + 1) X (H + 1)大小。C++实现代码如下：

//函数名：integral
//作用：快速积分图计算实现
//参数：
//matInput：输入图像
//matOutput : 输出积分图（32位整型数据）
//返回值：无
//注：支持单通道8位灰度图像
void integral(cv::Mat& matInput, cv::Mat& matOutput)
{
	//构造积分图
	matOutput = cv::Mat::zeros(matInput.rows + 1, matInput.cols + 1, CV_32SC1);
	for (int r = 0; r < matInput.rows; r++)
	{
		int nLineACC = 0;
		for (int c = 0; c < matInput.cols; c++)
		{
			nLineACC += matInput.at<unsigned char>(r, c);
			matOutput.at<int>(r + 1, c + 1) = matOutput.at<int>(r, c + 1) + nLineACC;
		}
	}
}

（3）、基于积分图的快速均值滤波

//函数名：fastMeanFilter
//作用：基于积分图的快速均值滤波实现
//参数：
//matInput：输入图像
//matOutput : 输出图像
//cnWinSize : 滤波窗口尺寸
//返回值：无
//注：支持单通道8位灰度图像
void fastMeanFilter(cv::Mat& matInput, cv::Mat& matOutput, const int& cnWinSize)
{
	//滤波窗口为奇数
	assert(cnWinSize % 2 == 1);
 
	//构造输出图像
	matOutput = cv::Mat::zeros(matInput.rows, matInput.cols, matInput.type());
 
	int nAnchor = cnWinSize / 2;
	//边界处理
	cv::Mat matMakeBorder;
	makeRepeatBorder(matInput, matMakeBorder, nAnchor);
 
	//计算积分图
	cv::Mat matIntegral;
	integral(matMakeBorder, matIntegral);
 
	//使用行首地址存储法遍历滤波
	int ** pRow = new int*[matIntegral.rows];
	for (int r = 0; r < matIntegral.rows; r++)
		pRow[r] = (int*)(matIntegral.data + r * matIntegral.step[0]);
 
 
	float fNormal = 1.0f / (cnWinSize * cnWinSize);
	//滤波操作
	for (int r = 0; r < matMakeBorder.rows - cnWinSize + 1; r++)
	{
		unsigned char* pOut = matOutput.data + r * matOutput.step[0];
		for (int c = 0; c < matMakeBorder.cols - cnWinSize + 1; c++)
		{
			//求和
			float fSum = pRow[r][c] + pRow[r + cnWinSize][c + cnWinSize] - \
				pRow[r + cnWinSize][c] - pRow[r][c + cnWinSize];
			//求均值输出
			pOut[c] = (unsigned char)(fSum * fNormal);
		}
	}
 
	delete[] pRow;
}

经过积分图加速后，均值滤波算法耗时如下：

4、总结

（1）、均值滤波优势在于计算量小、有多种加速的方法，对于随机噪声、点噪声有一定的抑制作用；缺陷在于，算法没有保护图像细节，滤波的同时也把图像边缘细节进行了平滑。

（2）、使用积分图进行均值滤波加速，能够使算法耗时基本与滤波窗口尺寸无关。

（3）、博文中的算法绝对耗时的值是与计算机配置、编译器等相关的，但是，耗时的趋势基本是可信的。

技术交流合作QQ：3355138068