一、概要
本文利用图像灰度的方差以及灰度的平均值的物理意义，结合概率论，经过多次测试，构造出较为简单实用的空域局部降噪滤波器的函数解析式，最后在本文末给出代码演示。另外，本文构造的局部降噪滤波器适用范围是滤除图像的加性噪声干扰。
二、滤波器的设计
1、问题的背景
现有一幅原图像为F(x,y)，被加性噪声V(x,y)干扰后实际得到的图像为G(x,y),亦即G(x,y)=F(x,y)+V(x,y)。构造一个滤波器，使得图像G(x,y)通过该滤波器后得到的图像能够尽可能的接近F(x,y)。
2、均值滤波器的本质
设f(x,y)表示原始图像F(x,y)上某一点(x,y)的像素值，g(x,y)表示加噪图像G(x,y)上某一点(x,y)的像素值，v(x,y)表示加性噪声V(x,y)对图像上某一点(x,y)施加的干扰值。那么对于f(x,y)可以利用g(x,y)对其进行线性表示，f(x,y)的估计值如下所示：

图像灰度的均值反映的是一幅图像的亮暗程度，其方差反映的是图像的对比度，图像的细节。那么恢复一幅被加性噪声污染的图像可以基于以下考虑：

亦即恢复出来的图像从统计意义上来说必须是原始图像的无偏且一致的估计；从物理意义上来说，恢复出来的图像必须与原始图像的亮暗程度以及细节相吻合。事实上，我们无法准确的恢复图像F(x,y)中的每一个具体像素值f(x,y)，但是能够在统计意义上无限的接近原始图像，只要误差在可接受的范围内，我们就说这个图像是恢复成功的。
将g(x,y)对f(x,y)的线性估计关系式代入上述方程组，随后解得系数a，b的值如下所示：

最后得到的表达式如下所示：

从理论上来说，直接套用上式复原一幅被加性噪声干扰的图像，在图像的任意一个区域内都能在统计的意义上精确的复现原始图像的亮暗程度（灰度均值）以及图像细节（灰度方差）。但实际上无法准确的得知局部区域的噪声方差D[v(x,y)]及其均值E[v(x,y)]，因此考虑利用全局意义上的噪声方差以及噪声均值作为D[v(x,y)]的估计，亦即：

那么滤波器的表达式如下所示：

在代码实现中套用了上述公式对图像进行降噪，经过验证，发现降噪后的图像会出现错误的结果。经过分析得出造成该问题的原因是因为利用噪声在原始图像F(x,y)上的方差作为噪声在局部区域Sxy上方差的估计的这一近似方法并不准确，有时根号内的计算结果为负值（局部区域的噪声方差小于全局意义上的噪声方差）于是导致计算结果错误。另外，尽管也利用了噪声在原始图像F(x,y)上的均值作为噪声在局部区域Sxy上均值的估计，这样做虽然会带来一定的误差，但反映到图像中不外乎是图像的亮暗变化，而抑制噪声更关心的是噪声的方差对图像的影响，因为噪声的方差会对图像质量直接造成较大的影响（直接破坏图像的细节）
综上所述，需要对上式中有根号运算的系数做出修正。
3、修正滤波函数
将上述滤波函数改写为：

对上式求方差得

如上所示，现在问题转化为构造一个系数C，在抑制噪声方差D[v(x,y)]的同时，保留原始图像的方差D[f(x,y)]。即抑制噪声对该区域Sxy的影响的同时，保留该区域的细节信息。
为此，经过多次测试，外加理论推导，系数a的构造方法如下所示：

其中，参数A根据实际加性噪声干扰情况进行设定，一般来说是加性噪声越大，A值也应该越大。
综上，最终的局部降噪滤波器的函数解析式如下所示：

该滤波器函数在统计意义上保留了图像中某一区域的细节信息，并抑制了加性噪声，这在后面会说明其原理。
三、算法流程
算法流程如下所示：
1、利用一个尺寸为m·n的矩形窗口自左向右，自上而下的选取图像G(x,y)的局部区域Sxy
2、计算局部区域Sxy的统计特性，即平均灰度值E[g(x,y)]，方差D[g(x,y)]
3、代入如下公式，并计算局部区域Sxy中心点(x,y)的像素值g(x,y)对原始图像对应位置像素值f(x,y)的估计。

四、滤波效果演示及分析
以一幅泰勒斯威夫特的图像作为原始图像，输入图像如下所示：

随后加上均值为20，方差为400的高斯白噪声。加噪结果如下所示：

显然，加上高斯白噪声后的图像变得斑驳，这是噪声的方差对图像细节的干扰所致；同时图像的亮度变亮了，这是噪声的均值对图像影响的结果。
运用传统的低通滤波器对该图像进行滤波会导致图像细节成分的大量丢失，在这里以均值滤波器为例，其滤波后的结果如下所示：

显然均值滤波后的图像比较模糊，丢失了相当一部分的图像细节。
运用本文构造的局部降噪滤波器处理该图片得到的结果如下所示：

显然，该滤波器在对高斯白噪声进行抑制的同时保留了更多的图像细节，同时对于图像边缘有锐化效果。因此相比于均值滤波来说，滤波后的图像显得更为清晰。
下面分析本文构造的滤波器抑制噪声的同时保留图像细节的原因。滤波器解析式如下所示：

对上式求方差，并化简得：

并且有以下定义：

G的大小反映的是某一区域内的图像细节被增强的程度。显然G的取值与原始图像的方差D[f(x,y)]以及噪声方差D[v(x,y)]有关。D[f(x,y)]越大，降噪后的图像细节也就越明显。
1、当A满足下述关系时（图像的轮廓往往会满足这个条件），该滤波器能够抑制噪声并增强图像的细节。

2、当A满足下述关系时（图像平坦部分会满足这个条件），图像中的加性噪声被抑制。

这是因为图像某一区域的方差为0，在图像细节丢失的同时也消除了噪声方差对图像的影响，但仍保留了原始图像的亮暗信息，亦即：

综上。该滤波器在图像中细节较多的地方，进行降噪的同时能够增强图像细节；在细节较少的地方只做降噪处理，因此最后的视觉效果会比均值滤波等传统滤波器的滤波效果更好。
最后，经过多次验证，本文所给的滤波器在方差更大的加性噪声干扰下仍能保持较好的抑制噪声能力，并能突出图像细节。
四、代码演示

#include<iostream>
#include<opencv2/opencv.hpp>
#include<math.h>
using namespace cv;
using namespace std;
double generateGaussianNoise(double m, double sigma);
Mat addGaussianNoise(Mat &srcImag);
Mat Means(Mat src, int m, int n);//m是行，n是列
Mat Means_2(Mat src, int m, int n);//m是行，n是列
Mat Variance(Mat src, Mat input_1, Mat input_2, int m, int n);//m是行，n是列
int main()
{Mat srcImage = imread("T.jpg");Mat gray;cvtColor(srcImage, gray, CV_BGR2GRAY);imshow("gray", gray);Mat dstImage = addGaussianNoise(srcImage);imshow("添加高斯噪声后的图像", dstImage);Mat gray_src;cvtColor(dstImage, gray_src, CV_BGR2GRAY);imshow("灰度化图像", gray_src);Mat example;gray_src.copyTo(example);blur(example, example, Size(7, 7));imshow("均值滤波结果", example);//求局部区域的平均灰度值E(x)Mat means1 = Mat::zeros(gray_src.size(), CV_32FC1);means1 = Means(gray_src, 7, 7);//E[g(x,y)]//求图像局部区域的二阶原点矩E(x^2)Mat temp1 = Mat::zeros(gray_src.size(), CV_32FC1);temp1 = Means_2(gray_src, 7, 7);//E[g(x,y)^2]//求图像局部区域的方差Mat variance1 = Mat::zeros(gray_src.size(), CV_32FC1);variance1 = Variance(gray_src, temp1, means1, 7, 7);//D[g(x,y)]//滤波Mat rst = Mat::zeros(gray_src.size(), CV_32FC1);for (int row = 3; row < gray_src.rows - 3; row++){for (int col = 3; col < gray_src.cols - 3; col++){float a = variance1.at<float>(row, col) / 1500;float b = means1.at<float>(row, col) - a*means1.at<float>(row, col);rst.at<float>(row, col) = (a*gray_src.at<uchar>(row, col) + b)-20;//for (int i = -3; i <= 3; i++)//{//	for (int j = -3; j <= 3; j++)//	{//		rst.at<uchar>(row+i, col+j) = saturate_cast<uchar>(a*gray_src.at<uchar>(row+i, col+j) + b);//	}//}}}rst.convertTo(rst, gray_src.type());imshow("rst", rst);waitKey(0);system("pause");return 0;
}
Mat Means(Mat src, int m, int n)
{Mat means = Mat::zeros(src.size(), CV_32FC1);//求局部区域的平均灰度值E(x)for (int row = (m - 1) / 2; row < src.rows - (m - 1) / 2; row++){for (int col = (n - 1) / 2; col < src.cols - (n - 1) / 2; col++){float sum = 0;for (int i = -(m - 1) / 2; i <= (m - 1) / 2; i++){for (int j = -(n - 1) / 2; j <= (n - 1) / 2; j++){sum = sum + src.at<uchar>(row + i, col + j);}}means.at<float>(row, col) = sum / (m*n);}}return means;
}
Mat Means_2(Mat src, int m, int n)
{Mat temp = Mat::zeros(src.size(), CV_32FC1);for (int row = (m - 1) / 2; row < src.rows - (m - 1) / 2; row++){for (int col = (n - 1) / 2; col < src.cols - (n - 1) / 2; col++){float sum = 0;for (int i = -(m - 1) / 2; i <= (m - 1) / 2; i++){for (int j = -(n - 1) / 2; j <= (n - 1) / 2; j++){sum = sum + pow(src.at<uchar>(row + i, col + j), 2);}}temp.at<float>(row, col) = sum / (m*n);}}return temp;
}
Mat Variance(Mat src, Mat input_1, Mat input_2, int m, int n)
{Mat variance = Mat::zeros(src.size(), CV_32FC1);for (int row = (m - 1) / 2; row < src.rows - (m - 1) / 2; row++){for (int col = (n - 1) / 2; col < src.cols - (n - 1) / 2; col++){variance.at<float>(row, col) = input_1.at<float>(row, col) - pow(input_2.at<float>(row, col), 2);}}return variance;
}
//生成高斯噪声
double generateGaussianNoise(double mu, double sigma)
{//定义小值const double epsilon = numeric_limits<double>::min();static double z0, z1;static bool flag = false;flag = !flag;//flag为假构造高斯随机变量Xif (!flag)return z1 * sigma + mu;double u1, u2;//构造随机变量do{u1 = rand() * (1.0 / RAND_MAX);u2 = rand() * (1.0 / RAND_MAX);} while (u1 <= epsilon);//flag为真构造高斯随机变量z0 = sqrt(-2.0*log(u1))*cos(2 * CV_PI*u2);z1 = sqrt(-2.0*log(u1))*sin(2 * CV_PI*u2);return z0*sigma + mu;
}
//为图像添加高斯噪声
Mat addGaussianNoise(Mat &srcImag)
{Mat dstImage = srcImag.clone();int channels = dstImage.channels();int rowsNumber = dstImage.rows;int colsNumber = dstImage.cols*channels;//判断图像的连续性if (dstImage.isContinuous()){colsNumber *= rowsNumber;rowsNumber = 1;}for (int i = 0; i < rowsNumber; i++){for (int j = 0; j < colsNumber; j++){//添加高斯噪声int val = dstImage.ptr<uchar>(i)[j] +generateGaussianNoise(20, 20);if (val < 0)val = 0;if (val>255)val = 255;dstImage.ptr<uchar>(i)[j] = (uchar)val;}}return dstImage;
}