非局部均值去噪（Non-Local Means，简称NLM）是一种用于图像处理的去噪算法，特别适用于去除图像中的高斯噪声。它的核心思想是考虑图像中的每个像素，并将其与图像中其他位置的相似区域进行比较。不同于传统的局部去噪方法，NLM算法利用了图像中更广泛区域的信息，从而更好地保持了图像的细节和结构。

NLM算法的基本步骤：

1. 选取搜索窗口：对于图像中的每个像素，选取一个围绕该像素的搜索窗口。
2. 计算相似性：对于搜索窗口内的每个像素，计算它与目标像素周围的小区域（局部邻域）的相似度。
3. 相似度加权：使用计算出的相似度作为权重，对搜索窗口内的像素进行加权平均，得到目标像素的去噪值。
4. 归一化处理：对所有的权重进行归一化处理，确保它们的总和为1。

NLM算法的优点：

• 保留细节：由于考虑了图像中更广泛的区域，NLM能够更好地保留图像细节和结构。
• 适用性广：适用于不同类型的噪声，特别是在处理高斯噪声方面效果显著。

缺点：

• 计算量大：NLM算法需要处理大量的像素和邻域比较，计算量较大，这在大尺寸图像或高分辨率图像处理时尤其明显。
• 参数选择：算法的效果很大程度上取决于参数的选择，如搜索窗口的大小、局部邻域的大小等。

对于一幅真实图像f(x)，假设图像收到的噪声干扰为加性高斯白噪声φ(x),最终结果为g(x)。

其中φ(x)为均值为零的高斯白噪声。

而加权平均是利用了图像中的自相似信息，根据像素之间的相似程度设置权值的大小。

非局部均值算法由此这两种算法结合而来。如何度量像素之间的相似性，或者构造权值函数，这是非局部均值算法的重点。

局部去噪虽然去噪效果好，但是在精细结构，细节信息和纹理方面明显不足。

非局部均值去噪则对于需要处理的目标像素点，用一个窗口区域遍历每个像素点，去计算窗口区域与目标像素点窗口区域的相似度，来确定遍历像素点的加权系数来估计目标像素点。

相似度采用高斯加权的欧氏距离来计算。即。

相似度大的像素点获得大的权值，相似度小的像素点获得小的权值。

#include "NLM.cuh"

__device__ float uchar3_dis(uchar3& a, uchar3& b) {

	return (a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y) + (a.z - b.z) * (a.z - b.z);

}


__global__ void nlmKernel(uchar3* input, uchar3* output, 
	int width, int height, int search_window_range,int patch_size, float h) {

	int x = blockIdx.x * blockDim.x + threadIdx.x;
	int y = blockIdx.y * blockDim.y + threadIdx.y;

	int half_pitch = patch_size / 2;
	int half_window = search_window_range / 2;
	
	if (x >= width || y >= height) { return; }
	float sumWeights = 0.0f;
	float3 result = make_float3(0.0f, 0.0f, 0.0f);

	//大循环当前像素点所对应的搜索窗口遍历
	for (int j = -half_window; j <= half_window; j++) {
		for (int i = -half_window; i <= half_window; i++) {
			float dis = 0.0f;
			float weight = 0.0f;

			//小循环
			for (int dy = -half_pitch; dy <= half_pitch; dy++) {
				for (int dx = -half_pitch; dx <= half_pitch; dx++) {
					int nx = x + i + dx;
					int ny = y + j + dy;
					int cx = x + dx;
					int cy = y + dy;
					if (nx >= 0 && nx < width && ny>=0 && ny < height) {
						if (cx >= 0 && cx < width && cy >= 0 && cy < height) {
							//计算欧氏距离
							dis += uchar3_dis(input[ny * width + nx], input[cy * width + cx]);
						}
					}
				}
			}
			//权重计算
			weight = exp(-dis / (h * h));
			//邻域窗口的中心像素点
			int nnx = x + i;
			int nny = y + j;
			if (nnx >= 0 && nnx < width && nny >= 0 && nny < height) {
				uchar3 nnpixel = input[nny * width + nnx];

				result.x += weight * nnpixel.x;
				result.y += weight * nnpixel.y;
				result.z += weight * nnpixel.z;
				sumWeights += weight;
			}
		}
	}
	if (sumWeights > 0) {
		result.x /= sumWeights;
		result.y /= sumWeights;
		result.z /= sumWeights;
	}

	output[y*width+x] = make_uchar3(static_cast<unsigned char>(result.x), static_cast<unsigned char>(result.y), static_cast<unsigned char>(result.z));
}


cv::Mat NLM_demo(cv::Mat& noise_img) {
	//一些参数设置   int search_window_range,int patch_size, float h
	//窗口尺寸
	int patch_size = 3;
	//搜索窗口大小
	int search_window_range = 7;
	//高斯指数函数衰减速度 h = c * sigma
	float c = 6.0;  //在0~10之间，可调
	float sigma = 25.0; //噪声方差的平方根，之前加噪声是25，这里设为5
	float h = c * sigma;

	//将图像上传GPU 
	cv::cuda::GpuMat d_img;
	d_img.upload(noise_img);

	cv::cuda::GpuMat d_out(d_img.size(),d_img.type());

	//核函数参数
	dim3 block(16, 16);
	dim3 grid((d_img.cols + block.x - 1) / block.x, (d_img.rows + block.y - 1) / block.y);

	//处理核函数
	nlmKernel << <grid, block >> > (d_img.ptr<uchar3>(), d_out.ptr<uchar3>(), d_img.cols, d_img.rows, search_window_range, patch_size, h);

	//处理后下载返回
	cv::Mat h_out;
	d_out.download(h_out);


	return h_out;
}

目前还存在的问题单通道图像的降噪可能效果好一些，彩色图像，可能因为我加噪声的方式的问题。

这里我是这么加的。

	Mat img = imread(file);
	//添加噪声
	Mat noise_img = Mat(img.size(), img.type());
	   //生成随机数
	random_device rd;  
	mt19937 gen(rd());
	normal_distribution<> d(0, 25); //均值为0，标准差为25，产生正态分布
	//#pragma omp parallel for
	for (int i = 0; i < img.rows; ++i) {
		for (int j = 0; j < img.cols; ++j) {
			for (int c = 0; c < img.channels(); ++c) {
				double noiseVal = d(gen);
				noise_img.at<Vec3b>(i, j)[c] = saturate_cast<uchar>(img.at<Vec3b>(i, j)[c] + noiseVal);
			}		
				//saturate_cast<uchar>确保像素值在合法范围内（0-255）
		}
	}

三个通道分别加高斯随机，也不知道这样加对不对。

修改过后的代码，虽然目前能跑，但还需要调节一些参数，并且仍未做优化，比如使用积分图来计算两个窗口之间的相似性，以及使用cuda编程中的各中内存加速。但基本上NLM算法的原理是实现了。知道算法原理，如何利用目前的硬件资源对算法加速还是一件很难的事，后期会继续更新。

这里主要使用了一个在核函数调用的设备函数，一个核函数，调用核函数的主机函数void，可以加上extern c关键字在.cpp文件中调用。