非局部均值滤波（一）

算法原理简介
图像之中的像素点并不是孤立存在的，某一点的像素与别处的像素点存在某种关联，可以概括为灰度相关性和几何结构相似性。但是相似像素并不局限于某个局部区域，自然图像含有丰富的冗余信息，所以可以采用能够描述图像结构特征的图像块在整幅图像上寻求相似块。在去噪的同时可以保留能最大程度地保持图像地细节特征。其基本思想是：当前像素的估计值由图像中与它具有相似邻域结构的像素加权平均得到。
以下图为例，进行详细说明：

以上为例，假定原图像大小为7x7（图中的橙色部分），由于需要对边缘像素进行处理，所以需要将图像进行填充，这种填充采用的是边缘像素的镜像。像素x是需要处理的原图像中的第一个像素，这里把搜索区域设置为7x7，也就是说像素x的估计值由图中绿色部分的49个像素加权平均得到，而权重是通过计算邻域块之间的相似度决定（图中以x和y为中心的蓝色图像块，块大小设置为3x3），如果计算出的块与块相似度越高，那么该像素对所求像素的贡献越大。
参考代码：

#include <iostream>
#include <opencv2/opencv.hpp>
#include <cmath>
#include <cstdlib>
#include <ctime>
#include <algorithm>
#include <vector>
using namespace cv;
using namespace std;


int Rows;//图像的行数
int Cols;//图像的列数

//添加高斯噪声
void addNoise(Mat& src, Mat& dest, double sigma)
{
	Mat s;
	src.convertTo(s, CV_16S);
	Mat n(s.size(), CV_16S);
	//randn()将n填充为高斯分布均值为0，标准差为sigma的随机数
	randn(n, 0, sigma);
	Mat temp = s + n;//将原图与高斯噪声叠加在一起输出
	temp.convertTo(dest, CV_8U);
}


//非局部均值滤波
void NlMeans(Mat &src, Mat &dst, double h = 10,int templateWindowSize = 7, int searchWindowSize = 21)
{
	//增加边界宽度
	int r1 = searchWindowSize/ 2;
	int r2 = templateWindowSize/ 2;
	int r3 = -1 * r2;
	int r = r1 + r2;
	Mat srcCopy;
	copyMakeBorder(src, srcCopy, r, r, r, r, BORDER_DEFAULT);
	//使用一个矩阵用来存储图像
	int hh = r * 2 + Rows;
	double **matrix = new double *[hh];
	for (int k = 0; k < hh; k++)
	{
		matrix[k] = new double[hh];
	}
	for (int x = 0; x < hh; x++)
	{
		for (int y = 0; y < hh; y++)
		{
			matrix[x][y] = srcCopy.at<uchar>(x, y);
		}
	}
	//srcCopy.convertTo(srcCopy, CV_64F);
	/*cout << srcCopy.rows << "  " << srcCopy.cols << endl;
	imshow("image", srcCopy);
	waitKey(0);*/
	//非局部均值滤波
	int i, j,m,n,p,q; //用于循环时像素的索引
	double similarity; //图像块的相似度
	double w;  //参考块相当于当前块的权重
	double sum; //权重之和
	double s;   //权重乘以像素值之和
	double sub; //像素差值
	//vector<double>weight(searchWindowSize*searchWindowSize);
	//double **weight = new double *[searchWindowSize];//建立一个weight数组，用来存储权重
	//for (int k = 0; k < searchWindowSize; k++)
	//{
	//	weight[k] = new double[searchWindowSize];
	//}
	double f = templateWindowSize * templateWindowSize;
	double delta = 10; //高斯噪声的方差
	double delta2 = 2 * delta*delta;
	Mat dstImage = src.clone();
	//遍历图像的每一个像素点
	for (i = r; i < Rows+r; ++i)
	{
		for (j = r; j < Cols+r; ++j)
		{
			//遍历搜索域
			sum = 0;
			s = 0;
			for (m = i - r1; m <= i + r1; ++m)
			{
				for (n = j - r1; n <= j + r1; ++n)
				{
					//计算每一个参考块
					similarity = 0;
					for (p = r3; p <= r2; ++p)
					{
						for (q = r3; q <= r2; ++q)
						{
							//计算当前块与搜索域中参考块之间的相似度
							//similarity += pow((srcCopy.at<uchar>(m + p, n + q) - srcCopy.at<uchar>(i + p, j + q)),2)/f;
							sub=matrix[m + p][n + q] - matrix[i + p][j + q];
							similarity += sub*sub;
						}
					}
					//计算参考块相对于当前块的权重
					similarity=similarity/f;
					w = exp(max((sbimilarity -delta2), 0.0) / (h*h));
					w = 1 / w;
					//cout << w << endl;
					sum += w;
					//s += w * srcCopy.at<uchar>(m, n);
					s += w *matrix[m][n];
				}
			}
			dst.at<uchar>(i - r, j - r) = s / sum;
		}
	}
	//释放二维数组
	for (i = 0; i < hh; i++)
	{
		delete[]matrix[i];
	}
	delete matrix;
}

int main()
{
	Mat srcImg = imread("E:/vsfiles/lena1.jpg");
	Mat grayImg;
	cvtColor(srcImg, grayImg, CV_RGB2GRAY);
	resize(grayImg, grayImg, Size(256,256));
	Rows = grayImg.rows;//行数
	Cols = grayImg.cols;//列数
	Mat imgWithNoise;
	addNoise(grayImg,imgWithNoise,10);
	Mat dstImg = imgWithNoise.clone();
	NlMeans(imgWithNoise, dstImg);
	//fastNlMeansDenoising(imgWithNoise, dstImg,10);
	imshow("image1", grayImg);
	imshow("image2", imgWithNoise);
	imshow("image3", dstImg);
	waitKey(0);
	return 0;
}