NLM去噪算法

一、算法原理

非局部均值滤波(Non-Local Means,NLM)是Buades等人于2005年在论文“A non-local algorithm for image denoising”中提出的对传统邻域滤波方法的一种改进滤波，考虑到了图像的自相似性质，它充分利用了图像中的冗余信息，在去噪的同时能够最大程度的保持图像的细节特征。【论文及源码地址】

该算法需要计算图像中所有像素与当前像素之间的相似性，考虑到这个计算量与效率的问题，一般会设定两个固定大小的窗口，一个大的搜索窗口(D×D)和一个小的邻域窗口(d×d)，邻域窗口在搜索窗口中进行滑动，根据邻域间的相似性来确定对应中心像素对当前像素的影响度，也就是权值。

下图是NLM算法执行过程，大窗口是以目标像素为中心的搜索窗口，两个灰色小窗口分别是以x，y为中心的邻域窗口。其中以y为中心的邻域窗口在搜索窗口中滑动，通过计算两个邻域窗口间的相似程度为y赋以权值w(x,y) 。


基本原理

1. 该算法需要遍历整个原图像；首先取出一个原图像pixel，以该pixel坐标为中心，圈出一大、一小两个矩形。大的矩形表示纹理替换搜索区域R，小的矩形表示待处理pixel的纹理区域L。
2. 将R矩形区域，分成若干个和R矩形区域一样的大小的矩形L1。计算出每块L1和L之间的权重W。
3. 将L的像素值都设置为0，然后根据L与每块L1的权重W，叠加当前L的像素值为：w * L1。
4. 将L和每块L1之间的权重W，同样累加起来到 $W_{all}$。
5. 所有L1遍历完了之后，$L=\frac{L}{W_{all}}$。就得到了经过去噪处理之后的L区域像素。

二、代码实现

import numpy as np
import cv2
import time


def nlm(array,N,K,sigma):
    height = array.shape[0]
    width  = array.shape[1]

    pad_len = N+K
    arraypad = np.pad(array,pad_len,'constant',constant_values=0)
    yy = np.zeros(array.shape)
    B = np.zeros((height,width))
    for ny in range(-N,N+1):
        for nx in range(-N,N+1):
            ssd = np.zeros((height,width))
            for ky in range(-K,K+1):
                for kx in range(-K, K + 1):
                    ssd += np.square(arraypad[pad_len+ny+ky:height+pad_len+ny+ky,pad_len+nx+kx:width+pad_len+nx+kx] - arraypad[pad_len+ky:height+pad_len+ky,pad_len+kx:width+pad_len+kx])
            ex = np.exp(-ssd/(2*sigma*sigma))
            B +=ex
            yy += ex*arraypad[pad_len+ny:height+pad_len+ny,pad_len+nx:width+pad_len+nx]
    return yy/B

if __name__ == '__main__':
    t1 = time.time()
    noise_img = cv2.imread("./3.png",0)
    noise_img = noise_img/255.0
    img_nlm = nlm(noise_img,10,4,0.6)
    t2 = time.time()
    print("Cost time :%f s"%(t2-t1))
    cv2.imshow("noise_img", noise_img)
    cv2.imshow("nlm",img_nlm)
    cv2.waitKey(0)
    cv2.destroyAllWindows()

三、加速代码

运行速度提高了约5倍

#%%cython --a --compile-args=/openmp
import numpy as np
cimport numpy as np
cimport cython
from libc.math cimport exp
from cython.parallel import parallel, prange



@cython.boundscheck(False)
@cython.wraparound(False)
def nlm(np.ndarray[np.float64_t, ndim=2] array,int N,int K,float sigma):
    cdef:
        Py_ssize_t height = array.shape[0]
        Py_ssize_t width  = array.shape[1]

        int pad_len = N+K
        int pad_H = pad_len+height
        int pad_W = pad_len+width

        double value = 2*sigma*sigma
        Py_ssize_t ny,nx,ky,kx
        Py_ssize_t h=0, w=0, h1=0, w1=0, h2=0, w2=0, h3=0, w3=0

    arraypad = np.zeros((pad_H+pad_len,pad_W+pad_len),dtype=np.float64)
    arraypad[pad_len:pad_H,pad_len:pad_W] = array
    Y_array = np.zeros((height,width),dtype=np.float64)
    X_array = np.zeros((height,width),dtype=np.float64)

    ex = np.zeros((height,width),dtype=np.float64)
    yy = np.zeros((height,width),dtype=np.float64)
    B = np.zeros((height,width),dtype=np.float64)
    ssd = np.zeros((height,width),dtype=np.float64)
    ppd = np.zeros((height,width),dtype=np.float64)

    cdef:
        double [:,::1] pad_view = arraypad
        double [:,::1] Y_view = Y_array
        double [:,::1] X_view = X_array
        double [:,::1] yy_view = yy
        double [:,::1] ex_view = ex
        double [:,::1] B_view = B
        double [:,::1] ppd_view = ppd
        double [:,::1] ssd_view = ssd
        
    with nogil,parallel():
        for ny in prange(-N,N+1):
            for nx in range(-N,N+1):
                ssd_view[:,::1] = 0
                for ky in range(-K,K+1):
                    for kx in range(-K, K + 1):
                        Y_view[:,::1] = pad_view[pad_len+ky+ny:pad_H+ky+ny,pad_len+kx+nx:pad_W+kx+nx]
                        X_view[:,::1] = pad_view[pad_len+ky:   pad_H+ky,   pad_len+kx:   pad_W+kx]
                        for h in range(height):
                            for w in range(width):
                                ssd_view[h,w] +=(Y_view[h,w] - X_view[h,w])*(Y_view[h,w] - X_view[h,w])
                for h1 in range(height):
                    for w1 in range(width):
                        ex_view[h1,w1] = exp(-(ssd_view[h1,w1]/value))
                for h2 in range(height):
                    for w2 in range(width):
                        B_view[h2,w2] += ex_view[h2,w2]
                ppd_view[:,::1] = pad_view[pad_len+ny:pad_H+ny,pad_len+nx:pad_W+nx]
                for h3 in range(height):
                    for w3 in range(width):
                        yy_view[h3,w3] += ex_view[h3,w3] * ppd_view[h3,w3]
    return yy/B

test.py

import cv2
import time
from nlm import nlm




if __name__ == '__main__':
    t1 = time.time()
    noise_img = cv2.imread("./3.png",0)
    noise_img = noise_img/255.0
    img_nlm = nlm(noise_img,10,4,0.6)
    t2 = time.time()
    print("Cost time :%f s"%(t2-t1))
    cv2.imshow("noise_img",noise_img)
    cv2.imshow("img_nlm", img_nlm)
    cv2.waitKey(0)
    cv2.destroyAllWindows()

测试效果图：我查看了一些运行数据，exp与np.exp在浮点数运算效果有些微差别导致图片效果不佳。待改进。

【参考】
https://blog.csdn.net/weixin_30570101/article/details/97913312
https://max.book118.com/html/2016/1206/68911778.shtm