深度学习结合非局部均值滤波的图像去噪算法

其实这是半年之前完成的内容，一直懒着没有总结，今天看了看代码，发觉再不总结自己以后都看不懂了，故整理如下。

非局部均值是一种基于块匹配来确定滤波权值的。即先确定一个块的大小，例如7x7，然后在确定一个搜索区域，例如15x15，在15x15这个搜索区域中的每一个点，计算7x7的窗口与当前滤波点7x7窗口的相似性（使用绝对差和SAD，一般而言，窗口中各点的差值还需要乘以经高斯核生成的权重参数，离中心点越近，权重值越大一些），然后根据相似性值使用指数函数生成窗口中心点的权重参数，相似性越高，该中心点的权重越大，最后各中心点的加权平均就是最终滤波图像，能获得很好的视觉效果。

非局部均值的成功之处主要在于充分利用了块的相似性，而后续步骤由相似性计算对应权重值，按照经验使用指数函数，其参数h有着至关重要的作用，许多论文也是在h上面做改进。如果我们跳出加权平均和指数函数的思路，完全可以将含噪图像所有相邻点的像素值、相似性值、距离等做为输入送给深度学习网络，将原图像值作为输出进行训练啊，训练好的模型就可以直接用于滤波。

下面附一个简化版的python代码，经实测改进后的算法比原生的非局部均值滤波要好，里面的网络模型过于简单，想提升效果的自己修改调优吧。

注意使用的是python3环境

#coding:utf8
import cv2, datetime,sys,glob
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.cm as cm

from keras.models import Sequential, model_from_json
from keras.layers import Dense, Activation,Dropout,Flatten,Merge
from keras.callbacks import EarlyStopping
from keras.layers.convolutional import Convolution2D,Convolution3D

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 GetNlmData(I, templateWindowSize=4,  searchWindowSize=9):
    f = int(templateWindowSize / 2)
    t = int(searchWindowSize / 2)
    height, width = I.shape[:2]
    padLength = t + f
    I2 = np.pad(I, padLength, 'symmetric')
    I_ = I2[padLength - f:padLength + f + height, padLength - f:padLength + f + width]

    res = np.zeros((height, width, templateWindowSize+2, t+t+1, t+t+1))
    for i in range(-t, t + 1):
        for j in range(-t, t + 1):
            I2_ = I2[padLength + i - f:padLength + i + f + height, padLength + j - f:padLength + j + f + width]
            for kk in range(templateWindowSize):
                kernel = np.ones((2*kk+1, 2*kk+1))
                kernel = kernel/kernel.sum()
                res[:, :, kk, i+t, j+t] = cv2.filter2D((I2_-I_) ** 2, -1,  kernel)[f:f + height, f:f + width]
            res[:, :, -2, i+t, j+t] = I2_[f:f + height, f:f + width]-I
            res[:, :, -1, i+t, j+t] = np.exp(-np.sqrt(i**2+j**2))
    print(res.max(), res.min())
    return res

def zmTrain(trainX, trainY):
    model = Sequential()
    if 1:
        model.add(Dense(100, init='uniform', input_dim=trainX.shape[1]))
        model.add(Activation('relu'))
        model.add(Dense(50))
        model.add(Activation('relu'))
        model.add(Dense(1))
        model.compile(loss='mse', optimizer='adam', metrics=['accuracy'])
    else:
        with open('model.json', 'rb') as fd:
            model = model_from_json(fd.read())
            model.load_weights('weight.h5')
            model.compile(loss='msle', optimizer='adam', metrics=['accuracy'])

    early_stopping = EarlyStopping(monitor='val_loss', patience=5)
    hist =model.fit(trainX, trainY, batch_size=150, epochs=200, shuffle=True, verbose=2, validation_split=0.1
                    ,callbacks=[early_stopping])
    print(hist.history)

    res = model.predict(trainX)
    res = np.clip(np.round(res.ravel() * 255), 0, 255)
    print(psnr(res, trainY*255))
    return model
if __name__ == '__main__':
    sigma = 20.0
    if 1:                         #这部分代码用于训练模型
        trainX = None
        trainY = None

        for d in glob.glob('./img/_*'):
            I = cv2.imread(d,0)
            I1 = double2uint8(I + np.random.randn(*I.shape) *sigma)
            data = GetNlmData(I1.astype(np.double)/255)
            s = data.shape
            data.resize((np.prod(s[:2]), np.prod(s[2:])))
            
            if trainX is None:
                trainX = data
                trainY = ((I.astype(np.double)-I1)/255).ravel()
            else:
                trainX = np.concatenate((trainX, data), axis=0)
                trainY = np.concatenate((trainY,  ((I.astype(np.double)-I1)/255).ravel()), axis=0)
            
        
        model = zmTrain(trainX, trainY)
        with open('model.json', 'wb') as fd:
            #fd.write(model.to_json())
            fd.write(bytes(model.to_json(),'utf8'))
        model.save_weights('weight.h5')
    if 1:                       #滤波
        with open('model.json', 'rb') as fd:
            model = model_from_json(fd.read().decode())
            model.load_weights('weight.h5')
        I = cv2.imread('lena.jpg', 0)
        I1 = double2uint8(I + np.random.randn(*I.shape) * sigma)

        data= GetNlmData(I1.astype(np.double)/255)
        s = data.shape
        data.resize((np.prod(s[:2]), np.prod(s[2:])))
        res = model.predict(data)
        res.resize(I.shape)
        res = np.clip(np.round(res*255 +I1), 0, 255)
        print('nwNLM PSNR', psnr(res, I))
        res = res.astype(np.uint8)
        cv2.imwrite('cvOut.bmp', res)