其实这是半年之前完成的内容，一直懒着没有总结，今天看了看代码，发觉再不总结自己以后都看不懂了，故整理如下。
非局部均值是一种基于块匹配来确定滤波权值的。即先确定一个块的大小，例如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)
