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。也就是说，算法的耗时与滤波窗口的尺寸相关。算法耗时情况如下表：
| 滤波尺寸 | 3 | 5 | 7 | 9 | 11 |
| 耗时(ms) | 2.431 | 5.296 | 9.949 | 17.06 | 27.37 |
那么，是否有更加优雅的方式实现呢？我这里提供一种基于积分图加速办法（查看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;
}经过积分图加速后，均值滤波算法耗时如下：
| 滤波尺寸 | 3 | 5 | 7 | 9 | 11 |
| 耗时（ms） | 2.246 | 2.456 | 2.658 | 2.954 | 3.155 |
4、总结
（1）、均值滤波优势在于计算量小、有多种加速的方法，对于随机噪声、点噪声有一定的抑制作用；缺陷在于，算法没有保护图像细节，滤波的同时也把图像边缘细节进行了平滑。
（2）、使用积分图进行均值滤波加速，能够使算法耗时基本与滤波窗口尺寸无关。
（3）、博文中的算法绝对耗时的值是与计算机配置、编译器等相关的，但是，耗时的趋势基本是可信的。
技术交流合作QQ：3355138068
4万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
