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