论文名称：A non-local algorithm for image denoising

论文下载：https://www.researchgate.net/profile/Bartomeu-Coll/publication/4156453_A_non-local_algorithm_for_image_denoising/links/0f317534c2cac194e4000000/A-non-local-algorithm-for-image-denoising.pdfhttps://www.researchgate.net/profile/Bartomeu-Coll/publication/4156453_A_non-local_algorithm_for_image_denoising/links/0f317534c2cac194e4000000/A-non-local-algorithm-for-image-denoising.pdf

        05年的论文，比较老的算法，但也比较有代表性。 论文比较好理解，就是认为图像内部有很多相似的块，只要进行块间相似性匹配，如果差异比较小，权重就比较大，差异比较大，权重就比较小，根据相似性进行加权平均，得到滤波后的结果。算法对均值为0的高斯噪声效果较好。由于好理解，就直接贴论文核心部分。

        用欧式距离进行块间匹配，差异越小，权重越大，由于图像是有噪声的，所以即使没有噪声的两个块是完全一样的，由于噪声的存在，导致计算出的欧式距离不为0，其期望大概为2*σ2。所以最后在算权重时，需要多减去2*σ2。

参该网站上提供的代码，用Python重写了，运行速度很慢，用分辨率低的图做测试比较好。IPOL Journal · Non-Local Means Denoisinghttps://www.ipol.im/pub/art/2011/bcm_nlm/ 

import cv2
import os
import numpy as np


def AddGaussNoise(img, sigma, mean=0):
    # 大概率abs(noise) < 3 * sigma
    noise = np.random.normal(mean, sigma, img.shape)
    img = img.astype(np.float)
    img = img + noise
    img = np.clip(img, 0, 255)
    img = img.astype(np.uint8)
    return img

def AddGaussNoiseGray(img, sigma, mean=0):
    lab = cv2.cvtColor(img, cv2.COLOR_BGR2Lab)
    noise = np.random.normal(mean, sigma, lab[:, :, 0].shape)
    lab = lab.astype(np.float)
    lab[:, :, 0] = lab[:, :, 0] + noise
    lab[:, :, 0] = np.clip(lab[:, :, 0], 0, 255)
    lab = lab.astype(np.uint8)
    img = cv2.cvtColor(lab, cv2.COLOR_Lab2BGR)
    return img

def GetMeandiff(patch1, patch2):
    # patch1 = patch1.astype(float)
    # patch2 = patch2.astype(float)
    diff = patch1 - patch2
    diff = diff.flatten()
    diff = diff * diff
    diffmean = diff.mean()
    return diffmean

def CalculateWeightLut(sigma, h):
    weightLut = np.zeros((256 * 256), np.float)
    sigma2 = sigma * sigma
    h2 = h * h
    for i in range(256 * 256):
        tmp = -max(i - 2 * sigma2, 0.0) / h2
        weightLut[i] = np.exp(tmp)
        if weightLut[i] < 0.0001:
            break
    return weightLut

def NonLocalMeansColor(image, sigma, h, templateWindowSize, searchWindow):
    height, width = image.shape[0], image.shape[1]
    patchWin = templateWindowSize // 2
    searchWin = searchWindow // 2

    # Padding the image
    padLength = patchWin + searchWin
    img = cv2.copyMakeBorder(image, padLength, padLength, padLength, padLength, cv2.BORDER_CONSTANT, value=255)

    img = img.astype(np.float)
    tmpSum = np.zeros((height + 2 * padLength, width + 2 * padLength, 3), np.float)
    count = np.zeros((height + 2 * padLength, width + 2 * padLength, 3), np.int)

    weightLut = CalculateWeightLut(sigma, h)

    for j in range(height):
        for i in range(width):
            padj = j + padLength
            padi = i + padLength
            centerPatch = img[padj - patchWin: padj + patchWin + 1, padi - patchWin: padi + patchWin + 1, :]
            # print(centerPatch.shape)
            sumWeight = 0
            templatePixel = np.zeros((templateWindowSize, templateWindowSize, 3), np.float)
            for r in range(padj - searchWin, padj + searchWin):
                for c in range(padi - searchWin, padi + searchWin):
                    otherPatch = img[r - patchWin: r + patchWin + 1, c - patchWin: c + patchWin + 1, :]
                    diff = GetMeandiff(centerPatch, otherPatch)
                    diff = (int)(diff)
                    curWeight = weightLut[diff]
                    sumWeight += curWeight
                    templatePixel += otherPatch * curWeight

            if sumWeight > 0.0001:
                tmpSum[padj - patchWin: padj + patchWin + 1, padi - patchWin: padi + patchWin + 1, :] += templatePixel / sumWeight
                count[padj - patchWin: padj + patchWin + 1, padi - patchWin: padi + patchWin + 1] += 1

    outImg = tmpSum[padLength:padLength + height, padLength:padLength + width, :]
    outCnt = count[padLength:padLength + height, padLength:padLength + width, :]

    mask = outCnt > 0
    outImg = outImg * mask + image * (1 - mask)

    outCnt = outCnt * mask + 1 - mask

    outImg = outImg / outCnt
    outImg = np.clip(outImg, 0, 255)
    outImg = outImg.astype(np.uint8)

    return outImg

if __name__ == '__main__':
    img = cv2.imread('test3.jpg', 1)

    # img = cv2.resize(img, (600, 750//2), interpolation=cv2.INTER_AREA)
    print(img.shape)
    noiseImg = AddGaussNoise(img, 20, 0)

    denoise = NonLocalMeansColor(noiseImg, 20, 8, 5, 15)

    cv2.imwrite('test_gauss_noise_color.jpg', noiseImg)
    cv2.imwrite('test_gauss_denoise_color.jpg', denoise)

    denoisecv = cv2.fastNlMeansDenoisingColored(noiseImg, None, 10, 10, 5, 15)
    cv2.imwrite('test_gauss_noise_colorcv.jpg', denoisecv)

代码里写的只能处理三通道的彩色图，处理灰度图的话，需要仿造写个函数。运行效果如下，sigma为20的高斯噪声。

噪声图像，sigma=20

去噪效果

fastNlMeansDenoisingColored运行的结果

可以看到，对高斯噪声，已知sigma的情况下，去噪效果还是比较好的，OpenCV自带的部分区域噪声没有去干净，可能和参数有关，他们用的是快速算法，参数也不一样，下次再研究一下。