1. 非局部均值濾波

非局部均值濾波（Non-Local Means，NL-Means）是一種 非線性 的影像去噪演算法。它基於影像中的像素具有相似結構這一假設，利用影像的全域資訊來對影像進行去噪。

1.1 全域演算法 VS 局部演算法

非局部均值濾波在計算每個像素點的估計值時，會考慮影像中所有與該像素點具有 相似 鄰域結構的像素點。因此，非局部均值濾波是一種 全域演算法 。

那麽相對於全域演算法的 局部演算法 是什麽呢？ 局部演算法 是指僅利用影像局部資訊進行處理的演算法。例如其鄰域視窗內的資訊，來計算該像素點的估計值。常用的局部演算法包括：

局部演算法的優點是計算量小，速度快。但其缺點是去噪效果有限，容易造成影像細節遺失。

1.2 均值濾波和非局部均值濾波的區別

均值濾波：對於影像中的每個像素點，其濾波值是其鄰域視窗內所有像素點的平均值。

非局部均值濾波：該演算法需要計算影像中 所有像素 與當前像素之間的相似性。首先需要定義一個 大的搜尋視窗 和 一個小的鄰域視窗 。搜尋視窗用於搜尋與當前像素點具有相似鄰域結構的像素點，鄰域視窗用於計算像素點之間的相似度。鄰域視窗在搜尋視窗中滑動，對於搜尋視窗內的每個像素點，計算其與當前像素點的鄰域視窗的相似度，並將其作為 權重 。相似度越大則權重越大。

非局部均值濾波的基本原理與均值濾波類似都要取平均值，但是非局部均值濾波在計算中加入了每一個點的權重。

非局部均值濾波是一種比均值濾波更先進的影像去噪方法，但其計算量也更大。

1.3 非局部均值濾波的原理

非局部均值濾波的公式如下：

其中，w(x,y) 是一個權重，表示在原始影像 v 中，像素 x 和像素 y 的相似度。是像素 x 的搜尋視窗。是濾波後的影像。

常用的 相似度度量 方法包括：歐式距離、高斯相似度、L1 範數、L2 範數、MSE 等等。

衡量兩個影像塊的相似度最常用的方法是計算他們之間的歐氏距離：

其中：

非局部均值濾波的復雜度跟 影像的大小、搜尋視窗的大小、相似度計算方法、權重計算方法 密切相關。

2. 非局部均值濾波的實作

下面的例子，是在影像中添加斑點雜訊，然後用非局部均值濾波消除雜訊。

#include<opencv2/opencv.hpp>
#include<opencv2/core.hpp>
#include<opencv2/highgui.hpp>
#include<random>
usingnamespacestd;
usingnamespace cv;
voidaddSpeckleNoise(Mat& image, double scale, Mat &dst){
dst = image.clone();
RNG rng;
dst.forEach<Pixel>([&](Pixel &p, constint * position) -> void {
int row = position[0];
int col = position[1];
double random_value = rng.uniform(0.0, 1.0);
double noise_intensity = random_value * scale;
dst.at<Vec3b>(row, col) = dst.at<Vec3b>(row, col) + Vec3b(noise_intensity * 255, noise_intensity * 255, noise_intensity * 255);
});
}
//NL-means 演算法的實作
voidnonlocalMeansFilter(Mat& src, Mat& dst, int templeteWindowSize, int searchWindowSize, double h, double sigma = 0.0){
//鄰域的大小不能超過搜尋視窗的大小
if (templeteWindowSize > searchWindowSize){
cout << "searchWindowSize should be larger than templeteWindowSize" << endl;
return;
}
if (dst.empty())
dst = Mat::zeros(src.size(), src.type());
constint tr = templeteWindowSize >> 1;//tr為鄰域的中心位置
constint sr = searchWindowSize >> 1; //sr為搜尋域的中心位置
constint bb = sr + tr;//需增加的邊界寬度
constint D = searchWindowSize*searchWindowSize;//搜尋域中的元素個數
constint H = D / 2 + 1;//搜尋域中的中心點位置
constdouble div = 1.0 / (double)D;//均勻分布時，搜尋域中的每個點的權重大小
constint tD = templeteWindowSize*templeteWindowSize;//鄰域中的元素個數
constdouble tdiv = 1.0 / (double)(tD);//均勻分布時，搜尋域中的每個點的權重大小
//擴充邊界
Mat boardSrc;
copyMakeBorder(src, boardSrc, bb, bb, bb, bb, cv::BORDER_DEFAULT);
//weight computation;
vector<double> weight(256 * 256 * src.channels());
double* w = &weight[0];
constdouble gauss_sd = (sigma == 0.0) ? h : sigma;//高斯標準差
double gauss_color_coeff = -(1.0 / (double)(src.channels())) * (1.0 / (h*h));//高斯顏色系數
int emax=0;
//w[i]保存變異數，即鄰域平均歐氏距離對應的高斯加權權重，供後面計算出歐式距離後呼叫
for (int i = 0; i < weight.size(); i++){
double v = std::exp(max(i - 2.0*gauss_sd*gauss_sd, 0.0)*gauss_color_coeff);
w[i] = v;
if (v<0.001){
emax = i;
break;
}
}
for (int i = emax; i < weight.size(); i++)
w[i] = 0.0;
int height = src.rows;
int width = src.cols;
if (src.channels() == 3){
constint cstep = (int)boardSrc.step - templeteWindowSize * 3;
constint csstep = (int)boardSrc.step - searchWindowSize * 3;
#pragma omp parallel for
for (int j = 0; j<height; j++){
uchar* d = dst.ptr(j);
int* ww = newint[D];//D 為搜尋域中的元質數量，ww用於記錄搜尋域每個點的鄰域變異數
double* nw = newdouble[D];//根據變異數大小高斯加權歸一化後的權重
for (int i = 0; i<width; i++){
double tweight = 0.0;
//search loop
uchar* tprt = boardSrc.data + boardSrc.step * (sr + j) + 3 * (sr + i);
uchar* sptr2 = boardSrc.data + boardSrc.step * j + 3 * i;
for (int l = searchWindowSize, count = D - 1; l--;){
uchar* sptr = sptr2 + boardSrc.step * (l);
for (int k = searchWindowSize; k--;){
//templete loop
int e = 0;
uchar* t = tprt;
uchar* s = sptr + 3 * k;
for (int n = templeteWindowSize; n--;){
for (int m = templeteWindowSize; m--;){
// computing color L2 norm
e += (s[0] - t[0])*(s[0] - t[0]) + (s[1] - t[1])*(s[1] - t[1]) + (s[2] - t[2])*(s[2] - t[2]);//L2 norm
s += 3;
t += 3;
}
t += cstep;
s += cstep;
}
constint ediv = e*tdiv;
ww[count--] = ediv;
//get weighted Euclidean distance
tweight += w[ediv];
}
}
//weight normalization
if (tweight == 0.0){
for (int z = 0; z<D; z++) nw[z] = 0;
nw[H] = 1;
}else{
double itweight = 1.0 / (double)tweight;
for (int z = 0; z<D; z++) nw[z] = w[ww[z]] * itweight;
}
double r = 0.0, g = 0.0, b = 0.0;
uchar* s = boardSrc.ptr(j + tr); s += 3 * (tr + i);
for (int l = searchWindowSize, count = 0; l--;){
for (int k = searchWindowSize; k--;)
{
r += s[0] * nw[count];
g += s[1] * nw[count];
b += s[2] * nw[count++];
s += 3;
}
s += csstep;
}
d[0] = saturate_cast<uchar>(r);
d[1] = saturate_cast<uchar>(g);
d[2] = saturate_cast<uchar>(b);
d += 3;
}//i
delete[] ww;
delete[] nw;
}//j
} elseif (src.channels() == 1){
constint cstep = (int)boardSrc.step - templeteWindowSize;//在鄰域比較時，從鄰域的上一行末尾跳至下一行開頭
constint csstep = (int)boardSrc.step - searchWindowSize;//搜尋域迴圈中，從搜尋域的上一行末尾跳至下一行開頭
#pragma omp parallel for
for (int j = 0; j<height; j++){
uchar* d = dst.ptr(j);
int* ww = newint[D];
double* nw = newdouble[D];
for (int i = 0; i<width; i++){
double tweight = 0.0;
uchar* tprt = boardSrc.data + boardSrc.step * (sr + j) + (sr + i);
uchar* sptr2 = boardSrc.data + boardSrc.step * j + i;
for (int l = searchWindowSize, count = D - 1; l--;){
uchar* sptr = sptr2 + boardSrc.step * (l);
for (int k = searchWindowSize; k--;){
int e = 0;
uchar* t = tprt;
uchar* s = sptr + k;
for (int n = templeteWindowSize; n--;){
for (int m = templeteWindowSize; m--;){
e += (*s - *t)*(*s - *t);
s++;
t++;
}
t += cstep;
s += cstep;
}
constint ediv = e*tdiv;
ww[count--] = ediv;
tweight += w[ediv];
}
}
if (tweight == 0.0){
for (int z = 0; z<D; z++) nw[z] = 0;
nw[H] = 1;
}else{
double itweight = 1.0 / (double)tweight;
for (int z = 0; z<D; z++) nw[z] = w[ww[z]] * itweight;
}
double v = 0.0;
uchar* s = boardSrc.ptr(j + tr); s += (tr + i);
for (int l = searchWindowSize, count = 0; l--;){
for (int k = searchWindowSize; k--;){
v += *(s++)*nw[count++];
}
s += csstep;
}
*(d++) = saturate_cast<uchar>(v);
}//i
delete[] ww;
delete[] nw;
}//j
}
}
intmain(){
Mat src = imread(".../girl.jpg");
imshow("src", src);
Mat result;
Mat dst4;
addSpeckleNoise(src,0.5,dst4);
imshow("addSpeckleNoise", dst4);
nonlocalMeansFilter(dst4, result, 3, 15, 40,40);
imshow("removeSpeckleNoise", result);
waitKey(0);
return0;
}





















OpenCV 提供了非局部均值濾波演算法，並對其進行了加速。

intmain(){
Mat src = imread(.../girl.jpg");
imshow("
src", src);
Mat result;
Mat dst4;
addSpeckleNoise(src,0.5,dst4);
imshow("

addSpeckleNoise", dst4);
fastNlMeansDenoisingColored(dst4,result,40,40);
imshow("
removeSpeckleNoise2", result);
waitKey(0);
return 0;
}

3. 總結

非局部均值濾波能夠有效地去除影像中的各種雜訊，包括高斯雜訊、椒鹽雜訊、紋理雜訊等。非局部均值濾波能夠較好地保留影像的細節，對雜訊的型別和分布不敏感，具有較強的魯棒性。

當然，非局部均值濾波的缺點也很明顯：計算量大，容易受到影像邊緣的影響等等。

非局部均值濾波的計算量大、速度慢是可以透過減少搜尋視窗大小、采用快速相似度計算方法、利用影像的稀疏性、並列化計算、利用硬體加速等方法來加速。

推薦閱讀