機器學習在影象處理中有非常多的應用，運用機器學習（包括現在非常流行的深度學習）技術，很多傳統的影象處理問題都會取得相當不錯的效果。今天我們就以機器學習中的字典學習（Dictionary Learning）為例，來展示其在影象去噪方面的應用。文中程式碼採用Python寫成，其中使用了Scikit-learn包中提供的API，讀者可以從【2】中獲得演示用的完整程式碼（Jupyter notebook）。

一、什麼是字典學習？

字典學習 (aka Sparse dictionary learning) is a branch of signal processing and machine learning. 特別地，我們也稱其為a representation learning method. 字典學習 aims at finding a sparse representation of the input data (also known as sparse coding or 字典) in the form of a linear combination of basic elements as well as those basic elements themselves。These elements are called atoms and they compose a dictionary。字典中，some training data admits a sparse representation。The sparser the representation, the better the dictionary。如下圖所示，現在我們有一組自然影象，我們希望學到一個字典，從而原圖中的每一個小塊都可以表示成字典中少數幾個atoms之線性組合的形式。

二、字典學習應用於影象去噪的原理

首先來看看影象去噪問題的基本模型，如下圖所示，對於一幅噪聲影象y，它應該等於原影象x 假設噪聲w，而在這個關係中，只有y是已知的，去噪的過程就是在此情況下推測x的過程。

在此基礎上，我們通常把去噪問題看成是一個帶約束條件的能量最小化問題，即對下面這個公式進行最下好從而求出未知項x。最小化公式中的第一項表示x和y要儘量接近，否則本來一幅貓噪聲的影象，降噪之後變成了狗，這樣的結果顯然不是我們所期望的。第二項則表示對x的一個約束條件，否則如果沒有這一項，那麼x就變成了噪聲影象y，那去噪也就失去了意義。這又變成了一個最大後驗估計(Maximum-a-Posteriori (MAP) Estimation) 問題，關於MAP的概念你可以參考【3】。

關於第二項到底應該符合什麼標準，不同學者都提出了各自的觀點。其中一派就認為，x應該滿足“稀疏”這個條件，因為這也是自然影象中所普遍存在的一個現實。

如果所x是一個m維的訊號，D=[d1,..., dn]是一組標準基向量（basis vectors），其大小為m×p，也就是我們所說的字典。D is "adapted" to y if it can represent it with a few basis vectors, that is, there exists a sparse vector a，這裡a是一個p維向量，such that y≈Da。這裡a 就是所謂的稀疏編碼。

為什稀疏對於去噪是好的？顯然，a dictionary can be good for representing a class of signals, but not for representing white Gaussian noise。直覺上，也即是說，用字典來近似表達y的時候，就會略掉噪聲。

於是乎，Elad and Aharon在2006年便提出了用字典學習來進行去噪的技術方案。首先，我們從y中提取所有的overlapping的矩形視窗（例如8×8），然後求解如下矩陣分解（matrix factorization）問題

基於這一基本原理，字典學習還可以應用於影象的inpaint（也就是影象修復），例如下面是修復前和修復後的影象對比

再來看看區域性效果，可見毫無違和感

三、Scikit-learn實現的基於字典學習的影象去噪例項

首先引入所需的各種包，並讀入一張Lena影象。

1. print(__doc__)

2. from PIL import Image

3. import matplotlib.pyplot as plt

4. import numpy as np

5. import scipy as sp

6. from sklearn.decomposition import MiniBatchDictionaryLearning

7. from sklearn.feature_extraction.image import extract_patches_2d

8. from sklearn.feature_extraction.image import reconstruct_from_patches_2d

9. %matplotlib inline

10. from keras.preprocessing.image import load_img

11. # load an image from file

12. image = load_img('lena_gray_256.tif')

13. from keras.preprocessing.image import img_to_array

14. # convert the image pixels to a numpy array

15. image = img_to_array(image)

16. image = image[:,:,0]

17. print("original shape", image.shape)

接下來將原本取值為0~255之間整數值的畫素灰度，轉換到0~1之間的浮點數。你也可以用imshow將影象顯示出來看看效果。

1. image = image.astype('float32')

2. image/=255

3. plt.imshow(image, cmap='gray')

然後是向其中加入高斯噪聲，我們把含噪聲的影象放在本文最後以便直觀地對比去噪效果。

1. noise = np.random.normal(loc=0, scale=0.05, size=image.shape)

2. x_test_noisy1 = image + noise

3. x_test_noisy1 = np.clip(x_test_noisy1, 0., 1.)

4. plt.imshow(x_test_noisy1, cmap='Greys_r')

接下來需要使用Scikit-learn中的Extract_pathches_2d來從影象中提取patches。通常來說，你可以從一組或者一個影象資料集中提取patches，並由此來訓練字典。你也可以從一張影象中提取patches，注意它們都是overlapping的。你還可以將二者結合起來。這裡我們的處理方式是從原圖（為加噪聲）中提取patches。你可能會疑惑，因為原圖在實際應用中是未知的。後面我們也會演示從噪聲影象中提取patches並去噪的效果。

1. # Extract all reference patches from the original image

2. print('Extracting reference patches...')

3. patch_size = (5, 5)

4. data = extract_patches_2d(image, patch_size)

5. print(data.shape)

這裡data的shape是一個三維向量[x, y, z]，x表示patches的數量，y和z是每個小patch的寬和高。但在進行運算時，我們需要把每個小矩形“拉直”，於是有

1. data = data.reshape(data.shape[0], -1)

2. print(data.shape)

3. data -= np.mean(data, axis=0)

4. data /= np.std(data, axis=0)

接下來使用MiniBatchDictionaryLearning函式和fit函式，就可以學到一個字典V。

1. # #############################################################################

2. # Learn the dictionary from reference patches

3. print('Learning the dictionary...')

4. dico = MiniBatchDictionaryLearning(n_components=144, alpha=1, n_iter=500)

5. V = dico.fit(data).components_

你也可以檢視一下字典V的形狀，或者把它顯示出來看看具體內容。這裡我們設計的字典中一共有144個atoms，所以分12排來列印，每排12個atoms。

1. print(V.shape)

2. plt.figure(figsize=(4.2, 4))

3. for i, comp in enumerate(V[:144]):

4. plt.subplot(12, 12, i + 1)

5. plt.imshow(comp.reshape(patch_size), cmap=plt.cm.gray_r,

6. interpolation='nearest')

7. plt.xticks(())

8. plt.yticks(())

9. plt.suptitle('Dictionary learned from patches\n', fontsize=16)

10. plt.subplots_adjust(0.08, 0.02, 0.92, 0.85, 0.08, 0.23)

字典如下圖所示。

接下來要做的是從噪聲影象中提取同樣大小的patches，然後從字典中找到一組最逼近這個小patch的一組原子的線性組合，並用這個組合來重構噪聲影象中的小patch。

1. # #############################################################################

2. # Extract noisy patches and reconstruct them using the dictionary

3. print('Extracting noisy patches... ')

4. data = extract_patches_2d(x_test_noisy1, patch_size)

5. data = data.reshape(data.shape[0], -1)

6. intercept = np.mean(data, axis=0)

7. data -= intercept

這裡所使用的演算法是OMP。函式recontruct_from_patches_2d用於重構影象。方法transform的作用是Encode the data as a sparse combination of the dictionary atoms

1. print('Orthogonal Matching Pursuit\n2 atoms' + '...')

2. reconstructions = x_test_noisy1.copy()

3. dico.set_params(transform_algorithm='omp', **{'transform_n_nonzero_coefs': 2})

4. code = dico.transform(data)

5. patches = np.dot(code, V)

6. patches += intercept

7. patches = patches.reshape(len(data), *patch_size)

8. reconstructions = reconstruct_from_patches_2d(patches, (256, 256))

9. plt.imshow(reconstructions, cmap='Greys_r')

去噪後的影象由文末對比結果中的左下角影象給出。

四、從噪聲影象中學習字典

現實中，噪聲影象的原影象是無法獲取的，否則就不需要去噪了。但我們可以直接從噪聲圖中學習字典。於是修改之前的程式碼，我們從直接從噪聲影象中獲取用於字典學習的patches。

1. # Extract all reference patches from noisy image

2. print('Extracting reference patches...')

3. patch_size = (5, 5)

4. data = extract_patches_2d(x_test_noisy1, patch_size)

5. print(data.shape)

同樣進行一些預處理和標準化之類的工作。

1. data = data.reshape(data.shape[0], -1)

2. print(data.shape)

3. data -= np.mean(data, axis=0)

4. data /= np.std(data, axis=0)

解下來進行字典學習，並將學到的字典顯示出來。

1. # #############################################################################

2. # Learn the dictionary from reference patches

3. print('Learning the dictionary...')

4. dico = MiniBatchDictionaryLearning(n_components=144, alpha=1, n_iter=500)

5. V = dico.fit(data).components_

6. print(V.shape)

7. plt.figure(figsize=(4.2, 4))

8. for i, comp in enumerate(V[:144]):

9. plt.subplot(12, 12, i + 1)

10. plt.imshow(comp.reshape(patch_size), cmap=plt.cm.gray_r,

11. interpolation='nearest')

12. plt.xticks(())

13. plt.yticks(())

14. plt.suptitle('Dictionary learned from patches\n', fontsize=16)

15. plt.subplots_adjust(0.08, 0.02, 0.92, 0.85, 0.08, 0.23)

下圖便是我們從噪聲影象中學習到的字典。

再次使用之前的程式碼。其中，set_params的作用是Set the parameters of this estimator。其中，transform_n_nonzero_coefs 是 Number of nonzero coefficients to target in each column of the solution.  也就是稀疏解的稀疏程度。預設情況下，它等於 0.1 * n_features（features的數量）。該引數 is only used by algorithm=’lars’ and algorithm=’omp’ and is overridden by alpha in the omp case。

1. # #############################################################################

2. # Extract noisy patches and reconstruct them using the dictionary

3. print('Extracting noisy patches... ')

4. data = extract_patches_2d(x_test_noisy1, patch_size)

5. data = data.reshape(data.shape[0], -1)

6. intercept = np.mean(data, axis=0)

7. data -= intercept

8. print('Orthogonal Matching Pursuit\n2 atoms' + '...')

9. reconstructions_frm_noise = x_test_noisy1.copy()

10. dico.set_params(transform_algorithm='omp', **{'transform_n_nonzero_coefs': 2})

11. code = dico.transform(data)

12. patches = np.dot(code, V)

13. patches += intercept

14. patches = patches.reshape(len(data), *patch_size)

15. reconstructions_frm_noise = reconstruct_from_patches_2d(patches, (256, 256))

16. plt.imshow(reconstructions_frm_noise, cmap='Greys_r')

重構的去噪效果如下圖中的右下角影象所示。

最後，順便補充一句，你可以使用下面的程式碼來儲存已經得到的影象結果。

1. imgs = (reconstructions_frm_noise * 255).astype(np.uint8)

2. Image.fromarray(imgs).save('lena_denoise_from_noise.png')

讀者還可以參考Scikit-learn官方文件中給出的例子【1】，以瞭解其他引數或演算法的使用。