NLM 图像降噪算法以及 Python 实现

Edvard 的视界深度之旅

已于 2023-05-31 09:14:23 修改

阅读量1.9k

点赞数 1

分类专栏：图形图像算法 # 图像降噪文章标签： python 论文阅读算法计算机视觉图像处理

于 2023-05-04 23:21:52 首次发布

本文链接：https://blog.csdn.net/edvard_hua/article/details/130497283

版权

图形图像算法同时被 2 个专栏收录

9 篇文章 3 订阅

订阅专栏

图像降噪

4 篇文章 0 订阅

订阅专栏

Non-Local Means Denoising, IPOL 2011

论文链接：https://www.ipol.im/pub/art/2011/bcm_nlm/article.pdf

请添加图片描述

是啥

基于相似性原理的图像去噪算法，与传统的局部均值滤波（平均，中值）不同，NLM 可以处理复杂的噪声模型，更好的保留图像细节

总得来说，NLM 算法的优缺点如下

优点：

去除噪音的同时，能够有效的保留图像的细节和纹理，降噪效果较好
不需要先验知识，可以处理多种不同类型的图像噪声
原理简单，易于实现

缺点：

复杂度高，运算速度较慢
对噪声强度和分布敏感，需要对参数进行合理设置才有好效果
处理低对比度和均匀纹理时，效果不佳

怎么做

先来看看整个算法流程

算法流程

在这里插入图片描述

**参考窗口（Reference Patch）**为像素块的大小，为奇数，如 5，9

**搜索窗口（Search Window）**为邻域窗口移动的范围，为奇数，如 11，13

算法是逐像素点处理的，结合上面的流程图，可以分为以下几个步骤

生成窗口：以像素点为中心生成参考窗口大小的像素块，并以滑动窗口的方式在搜索窗口内生成候选像素块
相似度匹配：计算参考窗口与候选窗口之间的相似度，获得 Top-N 个相似窗口
加权平均：将相似像素块进行加权平均，再使用高斯核获得最终降噪后的像素点的新值
重构：新值替代第一步原图中心像素点，对图像所有像素都重复上面三步，得到降噪后的图片

NLM 采用滑动窗口来寻找相似块，下图以 y 为中心的像素块，它会从左向右，从上往下寻找
在这里插入图片描述

局部和非局部

局部滤波（Local Filter）是使用当前像素的临近像素值来对图像进行降噪

非局部滤波（Non-local Filter）在这里的意思是寻找了相似块来进行降噪

相似块如何降噪

如下图所示，图像内部是存在相似的图像块，如右上角三个蓝色的矩形框所示，这三个矩形框的像素值我们可以用 P0，P1，P2 表示，那么

P0 = v0 + n0

P1 = v1 + n1

P2 = v2 + n2

其中 v 表示像素值，而 n 表示噪音且服从独立高斯分布，我们找到这些相似图像块并且计算像素均值，那么随机噪声就会被削弱，而本身的像素值 v 会被保留下来，从而达到去除噪声的同时保留图像细节
请添加图片描述

非局部均值滤波（Non-local Means Filter）

欧几里得距离，用来度量相似的块

**相似块权重，**卷积后接幂函数获得

可用以下两个公式来描述非局部均值滤波的计算

在这里插入图片描述

u 表示降噪后的图像，v 表示带噪音的图像，而 w(x, xi) 则表示权重值

在计算权重的公式中，B 表示两个像素块的欧几里得距离，Z 表示归一化系数，为所有权重之和

h 越小，权重衰减速度越快，更倾向于叠加那些误差小的块，肉眼来看就是降噪强度低，细节保留较多

h 越大，权重衰减速度约慢，近似于等权叠加，肉眼来看就是降噪强度高，细节保留较少

获得相似块后，不是直接加权叠加获得新块，而是只有每个块的中心像素点才参与计算，因为就算是相似快，块内还是会有差别，叠加会让图像变得混乱

Python 代码实现

详细的可以看下面代码的实现

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]
            # TODO: filter2D 的作用
            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])
    # 原图像需要乘于最大权重参与计算
    average += (wmax*I)
    # sweight 为 weight 之和，用于计算 weight 的归一化
    sweight += wmax
    return average / sweight

if __name__ == '__main__':
    I = cv2.imread('lena.jpg', 0)
    sigma = 20.0
    I1 = double2uint8(I + np.random.randn(*I.shape) *sigma)
    print('噪声图像PSNR',psnr(I, 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))
    print('NLM PSNR',psnr(I, R3))

其中，关键的代码是下面这行

w = np.exp( # 使用 exp 来约束 w 取值范围为 0-1 之间
    # (I2_ - I_) ** 2 为计算欧几里得距离
    -cv2.filter2D((I2_ - I_) ** 2, -1, kernel) / h
)[f:f + height, f:f + width]