数字图像处理之高斯滤波
在前面的文章中,我们讲过了均值滤波的原理与实现,讲高斯滤波之前,我们先回顾一下均值滤波,其核心思路是取每一个像素点邻域的矩形窗口,计算矩形窗口内所有像素点的像素平均值,作为该点滤波之后的像素值。
比如对于3*3窗口,如上图所示,点P(x, y)滤波之后的像素值为:
更广泛的,对于(2n+1)*(2n+1)窗口,点P(x, y)滤波之后的像素值为:
我们可以把上式看作是矩形窗口区域内所有像素点的像素值加权和,每个点的权重为W(i, j)=1/(2n+1)2,所以把上式进一步普遍化,则有:
下面我们开始讲高斯滤波,高斯滤波与均值滤波类似,都是矩形窗口内所有像素点的像素值加权和,只不过其权重与均值滤波不一样。高斯滤波的权重服从二维正态分布,越靠近窗口中心点(也即当前滤波点),权重越大。对于(2n+1)*(2n+1)窗口,其权重计算公式如下,其中σ为标准差,σ越大,权重分布越均匀,滤波效果越好,同时图像越模糊,反之,σ越小,权重分布越偏向于窗口中心点,滤波效果越差,同时图像越能保留其原有清晰度。
实际上,为了使加权和之后的值不超过像素值本来的最大范围,需要对上式做一下归一化,使0<w(i,j)<1。如下式:
举个例子,画出当σ=2.0时15*15窗口的权重分布三维图像,如下图,可以看到越接近窗口的中心点(7, 7),权重值越大:
好了,讲完原理,接着我们基于C++与Opencv来实现高斯滤波,最原始的实现思路是:1. 计算滤波权重;2. 像实现均值滤波一样,对图像进行边缘扩充,使得图像原有的所有像素点都可以取到邻域(2n+1)*(2n+1)矩形窗口;3. 行循环与列循环,遍历原有的每一个像素点,计算其矩形窗口内像素点的加权和作为滤波值。
计算权重的代码如下:
#define Gaussian_Size 17 //高斯滤波的窗口尺寸
#define Gaussian_Size_2 (Gaussian_Size>>1) //窗口尺寸的一半
float Gaussian_Ker[Gaussian_Size][Gaussian_Size]; //保存权重的全局数组
void getGaussianArray(float sigma) //计算高斯滤波的卷积核
{
float sum = 0.0f;
const float sigma_2 = sigma*sigma; //计算σ*σ
const float a = 1.0 / (2 * 3.14159*sigma_2); //计算1/(2*π*σ*σ)
for (int i = 0; i < Gaussian_Size; i++)
{
float dy = i - Gaussian_Size_2; //计算i-n
for (int j = 0; j < Gaussian_Size; j++)
{
float dx = j - Gaussian_Size_2; //计算j-n
//计算权重
Gaussian_Ker[i][j] = a*exp(-((dx*dx + dy*dy) / (2.0*sigma_2)));
//计算所有权重的和,用于归一化
sum += Gaussian_Ker[i][j];
}
}
sum = 1.0 / sum;
for (int i = 0; i < Gaussian_Size; i++)
{
for (int j = 0; j < Gaussian_Size; j++)
{
Gaussian_Ker[i][j] *= sum; //对权重进行归一化
}
}
}
那么高斯滤波的代码已经呼之欲出了:
void gaussian_fiter(Mat src, Mat &dst)
{
Mat src_board;
//扩充边缘
copyMakeBorder(src, src_board, Gaussian_Size_2, Gaussian_Size_2, Gaussian_Size_2, Gaussian_Size_2, BORDER_REFLECT); //扩充边缘
dst = Mat::zeros(src.size(), CV_8UC1);
for (int i = Gaussian_Size_2; i < src_board.rows - Gaussian_Size_2; i++) //行循环
{
uchar *dst_p = dst.ptr<uchar>(i - Gaussian_Size_2);
for (int j = Gaussian_Size_2; j < src_board.cols - Gaussian_Size_2; j++) //列循环
{
float sum = 0.0;
for (int y = 0; y < Gaussian_Size; y++)
{
uchar *p = src_board.ptr<uchar>(i - Gaussian_Size_2 + y);
for (int x = 0; x < Gaussian_Size; x++)
{
sum += p[j- Gaussian_Size_2 +x] * Gaussian_Ker[y][x]; //计算窗口内像素点的加权和
}
}
dst_p[j - Gaussian_Size_2] = (uchar)sum;
}
}
}
使用不同大小的窗口与σ值,运行上述代码,分别对Lena图像进行高斯滤波,得到结果如下。由以下结果可知:σ值越大,滤波效果越好,图像越模糊;窗口取得越大,滤波效果越好,同样图像也会变得越模糊。
原图
3*3窗口,σ=0.5
3*3窗口,σ=1.5
3*3窗口,σ=3.5
5*5窗口,σ=0.5
5*5窗口,σ=1.5
5*5窗口,σ=3.5
17*17窗口,σ=0.5
17*17窗口,σ=1.5
17*17窗口,σ=3.5
上述实现方法,需要执行4层循环,外面两层的循环次数由图像的行与列决定,内两层则由窗口的行与列决定。特别的,当窗口尺寸增大时,计算耗时将大幅度增加,比如对于以上496*472的Lena图像,取3*3窗口,耗时约1.3 ms,取5*5窗口时,耗时约3.18 ms,取17*17窗口时,耗时约226 ms。所以当为了提高滤波效果增大窗口尺寸时,耗时也大幅度增加,严重影响计算的实时性。根据高斯滤波权重公式的可分离性(如下式),我们可以把加权和的计算过程分成两步:1. 先计算行向的加权和;2. 然后再对行向的加权和计算列向加权和。这样一来,相当于把二维的加权和转换为一维加权和,减小了计算的复杂度。
根据上述分离原理,对于(2n+1)*(2n+1)的窗口,图像中任一点P(x, y)的高斯滤波值可按下式计算,其中Conx(x,y)为行方向的一维加权和,Cony(x,y)为列方向的一维加权和,也即高斯滤波的结果。
上式中W(i)为一维权重,其计算与二维权重类似,i取值0~2n:
同样,讲完原理,上代码。一维权重的计算代码如下:
#define Gaussian_Size 17 //保存高斯滤波窗口尺寸的全局变量
#define Gaussian_Size_2 (Gaussian_Size>>1)
float Gaussian_Ker_XY[Gaussian_Size];
void getGaussianArray_XY(float sigma) //计算高斯滤波的卷积核
{
float sum = 0.0f;
const float sigma_2 = sigma*sigma;
const float a = 1.0 / (2 * 3.14159*sigma_2);
for (int i = 0; i < Gaussian_Size; i++)
{
float dx = i - Gaussian_Size_2;
Gaussian_Ker_XY[i] = a*exp(-dx*dx / (2.0*sigma_2));
sum += Gaussian_Ker_XY[i];
}
sum = 1.0 / sum;
for (int i = 0; i < Gaussian_Size; i++)
{
Gaussian_Ker_XY[i] *= sum;
}
}
分离的高斯滤波实现代码:
void gaussian_fiter_XY(Mat src, Mat &dst)
{
Mat src_board;
//边缘扩充
copyMakeBorder(src, src_board, Gaussian_Size_2, Gaussian_Size_2, Gaussian_Size_2, Gaussian_Size_2, BORDER_REFLECT); //扩充边缘
Mat dstx = Mat::zeros(src_board.size(), CV_8UC1);
//x方向一维卷积
for (int i = 0; i < src_board.rows; i++)
{
for (int j = Gaussian_Size_2; j < src_board.cols - Gaussian_Size_2; j++)
{
float sum = 0.0;
for (int x = 0; x < Gaussian_Size; x++)
{
sum += src_board.ptr<uchar>(i)[j - Gaussian_Size_2 + x] * Gaussian_Ker_XY[x];
}
dstx.ptr<uchar>(i)[j] = (uchar)sum;
}
}
dst = Mat::zeros(src.size(), CV_8UC1);
//y方向一维卷积
for (int i = Gaussian_Size_2; i < src_board.rows - Gaussian_Size_2; i++)
{
for (int j = Gaussian_Size_2; j < src_board.cols - Gaussian_Size_2; j++)
{
float sum = 0.0;
for (int y = 0; y < Gaussian_Size; y++)
{
sum += dstx.ptr<uchar>(i - Gaussian_Size_2 + y)[j] * Gaussian_Ker_XY[y];
}
dst.ptr<uchar>(i - Gaussian_Size_2)[j - Gaussian_Size_2] = (uchar)sum;
}
}
}
运行上述代码,分别使用3*3、5*5、17*17大小的窗口对496*472的Lena图像进行高斯滤波,耗时分别为1.33 ms、1.91 ms、5.55 ms左右,相比原来的1.3 ms、3.18 ms、226 ms,可以看到,窗口尺寸越大,加速效果越明显。把二维加权和运算分离成一维加权和,本身可以减小计算复杂度,减少耗时,在此基础上,还便于作进一步的优化,比如SSE指令优化、CUDA优化,在下一篇文章中,我们将详细讲解高斯滤波的SSE指令优化与CUDA优化,敬请期待哟~