NLM非局部均值滤波

• 非局部均值滤波(Non-Local Means,NLM)是Buades等人于2005年在论文“A non-local algorithm for image denoising”中提出的对传统邻域滤波方法的一种改进滤波，考虑到了图像的自相似性质，它充分利用了图像中的冗余信息，在去噪的同时能够最大程度的保持图像的细节特征。
• 该算法需要计算图像中所有像素与当前像素之间的相似性，考虑到这个计算量与效率的问题，一般会设定两个固定大小的窗口，一个大的搜索窗口(D×D)和一个小的邻域窗口(d×d)，邻域窗口在搜索窗口中进行滑动，根据邻域间的相似性来确定对应中心像素对当前像素的影响度，也就是权值。
  基本算法流程
  2018092917323925.png
  非局部均值滤波的基本原理与均值滤波类似，都是要取平均值，但是非局部均值滤波在计算中加入了每一个点的权重值，所以能够保证在相邻且相差很大的点在方框中求平均值时相互之间的影响减小，也就对图像边缘细节部分保留很多，这样图像看起来会更清晰。

非局部均值滤波的算法详细步骤

1. 首先在一个点A周围取一个大的框（搜索框），设边长为s，A在方框的中心，然后再在方框中取小的方框，即相似框，设边长为d
2. 那么在A周围也有一个边长为d的方框，然后在大方框中找到所有边长为d的小方框的组合（就是一个小正方形在一个大正方形中到处移动，记录小正方形中心点的坐标就行了），设小方框的中心点为B，分别于A周围的相似框求减法，并且加入高斯核计算得到的加权值，这样可以计算出一个二维数组，里面存放着各个点的差值乘以权重后的值，加入高斯核主要是因为距离中心点距离不同对中心点的影响大小也不同，而且高斯核的权重和是1，所以就不用再归一化了。
3. 然后将这个二维数组求和并平均，得到的值就是这个相似框的中心点B对于A的权重值。计算出A周围所有点的权重值，其实这个时候这个值和权重是成反比的，以A本身为例（以A为中心点的相似框），计算出来A对于A的所谓权重值是零。然后根据计算出来的值用一个指数减函数就得到了成正比的权重关系，具体的函数见下面的代码，w=exp(-d/h)，就是这个，其中d就是计算出来的值啦，代入后w就是成正比的权重关系啦，h是一个滤波百分比值。可以先固定为一个常数。 而且这个计算出来w就是一个(0,1)的值哦，自动归一化啦
4. 然后就是根据得到的权重值以及各个点本身的灰度值计算出非局部均值滤波后A点的灰度值。
5. 以此类推，可以计算出图中所有点经过非局部均值滤波后的值
   优缺点

• 优点

1. 可以既去除噪声，又保留图像边缘细节
2. 当然去噪声指的一般是高斯白噪声，因为高斯白噪声的均值是0，所以求和取平均会比较有效果

• 缺点

1. 计算起来很慢
2. 如果图像像素点比较多，而且计算的时候取的框还比较大的话，那么计算一般几分钟是要的了。

# coding:utf8
import cv2
import numpy as np
def psnr(A, B):
    return 10 * np.log(255 * 255.0 / (((A.astype(np.float) - B) ** 2).mean())) / np.log(10)
def double2uint8(I, ratio=1.0):
    return np.clip(np.round(I * ratio), 0, 255).astype(np.uint8)
def make_kernel(f):
    kernel = np.zeros((2 * f + 1, 2 * f + 1))
    for d in range(1, f + 1):
        kernel[f - d:f + d + 1, f - d:f + d + 1] += (1.0 / ((2 * d + 1) ** 2))
    return kernel / kernel.sum()
def NLmeansfilter(I, h_=10, templateWindowSize=5, searchWindowSize=11):
    f = templateWindowSize // 2
    t = searchWindowSize // 2
    height, width = I.shape[:2]
    padLength = t + f
    I2 = np.pad(I, padLength, 'symmetric')
    kernel = make_kernel(f)
    h = (h_ ** 2)
    I_ = I2[padLength - f:padLength + f + height, padLength - f:padLength + f + width]
    average = np.zeros(I.shape)
    sweight = np.zeros(I.shape)
    wmax = np.zeros(I.shape)
    for i in range(-t, t + 1):
        for j in range(-t, t + 1):
            if i == 0 and j == 0:
                continue
            I2_ = I2[padLength + i - f:padLength + i + f + height, padLength + j - f:padLength + j + f + width]
            w = np.exp(-cv2.filter2D((I2_ - I_) ** 2, -1, kernel) / h)[f:f + height, f:f + width]
            sweight += w
            wmax = np.maximum(wmax, w)
            average += (w * I2_[f:f + height, f:f + width])
    return (average + wmax * I) / (sweight + wmax)
if __name__ == '__main__':
    I = cv2.imread('123.jpg', 0)
    sigma = 5.0
    I1 = double2uint8(I + np.random.randn(*I.shape) * sigma)#np.random.randn()函数所产生的随机样本基本上取值主要在-1.96~+1.96之间
    print ('噪声图像PSNR', psnr(I, I1))#I1表示噪声图像
    R1 = cv2.medianBlur(I1, 5)
    print ('中值滤波PSNR', psnr(I, R1))
    R2 = cv2.fastNlMeansDenoising(I1, None, sigma, 5, 11)
    print ('opencv的NLM算法', psnr(I, R2))
    R3 = double2uint8(NLmeansfilter(I1.astype(np.float), sigma, 5, 11))
    cv2.imshow('R3',R3)
    cv2.waitKey(0)