字典学习在图像和信号处理中是一种重要的算法，常常用于图像去噪、分类等，其中图像去噪可以认为是一种无监督学习技术。接下来简单介绍字典学习原理，并使用Python进行灰度图像去噪。
1 字典学习
灰度图像可以认为是二维信号，可以使用冗余字典和该字典下的稀疏编码来表示。
字典学习就是根据已知的数据找到合适的字典和其对应的稀疏编码，使误差尽可能的小。矩阵使用冗余字典和稀疏编码表示如图
1
所示，理想情况下希望重构的误差为
0
。字典学习方法为何可以用于图像去噪呢？因为针对带噪图像的噪声可以认为是稀疏的信号，原始图像不是稀疏的信号，字典学习和稀疏编码可以去除稀疏的信号，保留不带噪声的清晰图像，从而实现图像去噪的目的。
图1 矩阵使用冗余字典稀疏编码示意图 结合图1在应用于图像去噪时，和Y相同尺寸的带噪图像Y’，我们希望找到矩阵D和A，通过Y＝D＊A来重新得到不带噪声的原始图像Y。 2 字典学习图像去噪实战
在Python中可以使用sklearn库decomposition模块中的函数进行字典学习算法的应用，接下来加载所需要的函数，使用Python进行字典学习对图像去噪：
## 输出高清图像%config InlineBackend.figure_format = 'retina'%matplotlib inlineimport numpy as npimport matplotlib.pyplot as pltimport sklearn.decomposition as skd## 将2D图像转化为块的集合from sklearn.feature_extraction.image import extract_patches_2d## 从所有的块中重构图像from sklearn.feature_extraction.image import reconstruct_from_patches_2dfrom skimage.io import imread ## 从skimage库中引入读取图片的函数from skimage.color import rgb2gray from skimage.util import random_noisefrom skimage.metrics import peak_signal_noise_ratio
## 输出高清图像
%config InlineBackend.figure_format = 'retina'
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import sklearn.decomposition as skd
## 将2D图像转化为块的集合
from sklearn.feature_extraction.image import extract_patches_2d
## 从所有的块中重构图像
from sklearn.feature_extraction.image import reconstruct_from_patches_2d
from skimage.io import imread ## 从skimage库中引入读取图片的函数
from skimage.color import rgb2gray
from skimage.util import random_noise
from skimage.metrics import peak_signal_noise_ratio
上面加载的函数中extract_patches_2d函数是将2D的图像数据通过滑块采样的方式提取图像特征，对图像进行分块；reconstruct_from_patches_2d函数是对extract_patches_2d进行反向操作，得到原始图像大小的数据；rgb2gray函数是将RGB图像转化为灰度图像；random_noise函数是对图像数据添加噪声的函数。加载所需要的函数后，针对灰度图像数据具体的去噪过程，可以使用如下数据表所示：
灰度图像去噪步骤： |
(1)针对干净原始图像X添加零均值的高斯白噪声(或者其他类型的噪声)V，得到带噪图像Y。(2)针对带噪声图像使用大小为n*n(如8*8、10*10、12*12等)的图像块过采样处理，并且将每一个图像块转化为一个向量对待，得到我们的训练模型数据集M_1。(3)对数据集M1中的每个图像块进行中心化处理,即M_2=M_1-M_mean,M_mean为每个图像块的均值。(4)使用相应的字典学习算法，学习字典D和相应的稀疏编码A。(5)使用D*A重建数据集,并且将每个图像块重新加上它们的均值，即M = D*A+M_mean。(6)将去噪后的数据集M重新转化为去噪后的灰度图像。(7)使用数据可视化技术，评价图像去噪的效果。 |
下面结合以上的步骤使用Python进行字典学习，对图像去噪。
(1)读取图像然后转化为灰度图像，并且对图像添加高斯噪声：
im = imread("莱娜.tiff")orimga = rgb2gray(im)mode='gaussian' ## 图像添加噪声sigma = 30 ## 噪声的标准差sigma2 = sigma**2 / (255**2)noimga = random_noise(orimga,mode=mode, var=sigma2,seed=1243,clip=True)print("noimga",noimga.min(),"~",noimga.max())noimga 0.0 ~ 1.0使用rgb2gray()函数将RGB图像转化为灰度图像，使用random_noise()函数添加高斯噪声，mode参数指定噪声的形式，这里使用的是高斯噪声，var参数指定噪声的方差，clip=True表示将添加噪声后的数据转化到0～1之间。可以发现经过处理后的图像取值范围为[0,1]。
(2)将图像转化为样本图像小块，从原始图像上滑动采样10*10大小的图像块，可以得到253009个小图像块，将图像块展开，得到一个253009*100的矩阵，并对数据进行中心化处理。
## 将图像转化为图像块样本patch_size = (10,10) ## 块的大小## 使用带噪声图像训练字典noimgadata = extract_patches_2d(noimga,patch_size)noimgadata = noimgadata.reshape(noimgadata.shape[0],-1)noimgadata = noimgadata.Tinterp = noimgadata.mean(axis=0) ## 计算每个图像块的均值print("interp.shape:",interp.shape)noimgadata = noimgadata - interp ## 图像块减去相应的均值noimgadata = noimgadata.T ## 再次转置 print("训练模型数据 shape is :",noimgadata.shape)print('带噪声图像取值为:',noimgadata.min(),"~",noimgadata.max())print("noimgadata.shape:",noimgadata.shape)interp.shape: (253009,)训练模型数据 shape is : (253009, 100)带噪声图像取值为: -0.8323655493823148 ~ 0.8489631118015595noimgadata.shape: (253009, 100)使用extract_patches_2d()函数，可以将一个矩阵转化为按照指定图像块大小(10*10)的多维数组，然后使用noimgadata.mean( axis =0) 函数转化成训练数据大小为253009*100的矩阵。计算二维数组每个列的均值，noimgadata.T得到转置后的数据集。
(3)定义字典原子数量，并进行字典学习。
m = 16n_components = m**2 ## 字典原子数目n_iter=100 ## 迭代次数code,dictionary=skd.dict_learning_online(noimgadata,n_components=n_components, alpha=0.2,n_iter=n_iter,return_code=True, method='lars',n_jobs=6,random_state = 1234)print("dictionary shape = ",dictionary.shape)print("sparse code shape = ",code.shape)dictionary shape = (256, 100)sparse code shape = (253009, 256)
上面的程序使用
skd.dict_learning_online
进行在线字典学习，得到的字典包含
256(16*16)
个原子，即字典为
256*100
的矩阵，稀疏编码为
253009*256
的矩阵。在
skd.dict_learning_online()
函数中，
noimgadata
表示要进行字典学习的带噪声的数据集；
n_components=n_components
表示字典原子的数目，这里为
256
；
alpha=0.2
表示要进行字典学习时使用的惩罚参数；
n_iter=n_iter
指定了模型的迭代次数；
return_code=True
表示输出的结果中字典的稀疏编码。
(
4
)使用字典和稀疏编码重构去噪后的图像。
## 重建无噪声图像denoimgadata = np.dot(code,dictionary) denoimgadata = denoimgadata.Tdenoimgadata = denoimgadata + interpdenoimgadata = denoimgadata.Tpatch = denoimgadata.reshape(len(denoimgadata),*patch_size)denoimga = reconstruct_from_patches_2d(patch,orimga.shape)print("denoimga.shape",denoimga.shape)denoimga.shape (512, 512)代码中np.dot(code,dictionary)表示使用字典与稀疏编码相乘来重构去除噪声后的训练数据集，再加上均值，得到去噪后的数据，最后使用reconstruct_from_patches_2d()函数将数据转化为原始图像大小，得到512*512的图像，该函数是extract_patches_2d()函数的逆操作。带噪数据的去噪工作即完成。
(5)将原始图像、带噪图像、去噪后的图像一起可视化，分析可视化的结果，检查去噪的效果。
position = [200,200] ## 图片位置 size = 80 ## 图像大小## 提取三个图像块orimpatch = orimga[position[0]:position[0]+size,position[1]:position[1]+size]noimpatch = noimga[position[0]:position[0]+size,position[1]:position[1]+size]denoimpatch = denoimga[position[0]:position[0]+size,position[1]:position[1]+size]## 查看6幅图像plt.figure(figsize=(12,8))ax1 = plt.subplot(231) ## 原始图像plt.imshow(orimga,cmap=plt.cm.gray)## 绘制图像块的位置rect = plt.Rectangle(position,size,size,edgecolor="r",facecolor="none")ax1.add_patch(rect)plt.axis("off");plt.title("Original Image")ax2 = plt.subplot(232) ## 噪声图像plt.imshow(noimga,cmap=plt.cm.gray)rect = plt.Rectangle(position,size,size,edgecolor="r",facecolor="none")ax2.add_patch(rect)plt.axis("off");plt.title("Noisy Image $\sigma = 30$")ax3 = plt.subplot(233) ## 去噪图像图像plt.imshow(denoimga,cmap=plt.cm.gray)rect = plt.Rectangle(position,size,size,edgecolor="r",facecolor="none")ax3.add_patch(rect)plt.axis("off");plt.title("Denoisy Image")## 绘制三个图像块plt.subplot(234) ## 原始图像块plt.imshow(orimpatch,cmap=plt.cm.gray)plt.axis("off");plt.title("Original Image patch")plt.subplot(235) ## 噪声图像块plt.imshow(noimpatch,cmap=plt.cm.gray)plt.axis("off");plt.title("Noisy Image patch")plt.subplot(236) ## 去噪图像图像块plt.imshow(denoimpatch,cmap=plt.cm.gray)plt.axis("off");plt.title("Denoisy Image patch")plt.subplots_adjust(wspace = 0.08,hspace = 0.08)plt.show()得到的图像如图2所示，共绘制了6副子图，每列分别对应原始图像、带噪图像和去噪后图像。第二行图像是第一行图像的局部放大图。
图2字典学习去噪效果
从图2对比可以发现，局部放大图中，去噪后图像和原始图像更加接近。即去噪后的图像更接近于原始图像，该算法的去噪效果很直观，对带噪图像质量有很大的提升。下面计算去噪前后的PSNR取值的大小，从数值上度量去噪效果。
## 计算去噪前后的PSNRprint("加噪后的PSNR:",peak_signal_noise_ratio(orimga,noimga),"dB")print("去噪后的PSNR:",peak_signal_noise_ratio(orimga,denoimga),"dB")加噪后的PSNR: 18.734043388447393 dB去噪后的PSNR: 28.46829312043542 dB
3
可视化字典原子和稀疏编码
可以使用下面的程序可视化字典原子的情况，运行程序后可获得如图3所示的图像。
## 可视化字典plt.figure(figsize=(10,10))for ii in range(256): plt.subplot(16,16,ii+1) plt.imshow(dictionary[ii,:].reshape(10,10),cmap=plt.cm.gray) plt.axis("off")plt.subplots_adjust(wspace=0.05,hspace=0.05)plt.show()
图3 256个字典原子可视化
使用下面的程序可以可视化稀疏编码的情况，只可视化前4个样本的稀疏编码，运行程序后可获得如图4所示的图像。可见每个稀疏编码只有几个非0值。
x = range(code.shape[1])plt.figure(figsize=(10,6))for ii in range(4): plt.subplot(4,1,ii+1) plt.stem(x,code[ii,:],use_line_collection = True) plt.grid()plt.tight_layout()plt.show()
图4 稀疏编码可视化
