圖像降噪算法——稀疏表達：K-SVD算法

• 圖像降噪算法——稀疏表達：K-SVD算法
• • 1. 基本原理
  • 2. python代碼
  • 3. 結論

圖像降噪算法——稀疏表達：K-SVD算法

為了完善下自己降噪算法知識版圖的完整性，我打算花一個周末的時間再了解下基于稀疏表達和低秩聚類這兩種原理實作的圖像降噪算法，因為學習的時間并不長，也沒有花太多時間去做實驗，是以對算法了解得可能比較膚淺，還願讀者見諒。

這裡我分享幾篇很優秀的部落格，我對稀疏表達概念的認識目前也就是來自于這幾篇部落格：

字典學習（Dictionary Learning, KSVD）詳解

從稀疏表示到K-SVD，再到圖像去噪

稀疏表達的意義在于？為什麼稀疏表達得到廣泛的應用？

1. 基本原理

這裡我們先來回答幾個問題：

（1）什麼是稀疏表達？

稀疏表達的定義是用較少的基本信号的線性組合來表達大部分或者全部的原始信号。在稀疏表達算法中有"字典"這樣一個概念，這個“字典”和我們現實生活中使用的字典作用是相同的，現實生活中，我們通過對詞條進行組合可以表達我們想要的資訊，同樣地，在稀疏表達算法中，我們對“原子”進行線性組合就可以還原我們想要的信号。稀疏表達在機器學習領域有着廣泛的應用，諸如壓縮感覺、特征提取等，圖像去噪隻是其應用的一個方面。

（2）稀疏表達為什麼可以進行圖像去噪？

噪聲圖像是原始圖像和噪聲合成的圖像，原始圖像認為是可稀疏的，即可以通過有限個"原子"進行表示，而噪聲是随機而不可稀疏的，是以通過提取圖像的稀疏成分，再利用稀疏成分來重構圖像，在這個過程中，噪聲被當做觀測圖像和重構圖像之間的殘差而被丢棄，進而起到降噪的作用。

在圖像去噪領域，盡管CNN的方法已經占據了半壁江山，但是CNN方法的一個弊端就是需要大量訓練資料，在難以采集大量訓練資料的場景中，例如某些醫學圖像，稀疏表達的方法仍然起到重要作用。

（3）稀疏表達和轉換域濾波方法的差別？

轉換域濾波指的是離散傅裡葉變換、小波變換這樣一些方法，這些方法獲得的是一系列正交基，而稀疏表達獲得的“字典”則是一些列非正交基。正交基往往隻能表示圖像的某一個特征而不能夠同時表示其他特征，是以正交基的稀疏性不及非正交基。

下面咱來開始讨論基于稀疏濾波的K-SVD算法，K-SVD算法中最核心的部分就是字典學習，字典學習的過程就是基于樣本集 Y = { y i } i = 1 M \mathbf{Y}=\{\mathbf{y}_i\}^M_{i=1} Y={yi​}i=1M​尋找一組“超完備”的基向量作為字典矩陣 D = { d i } i = 1 K \mathbf{D}=\{\mathbf{d}_i\}^{K}_{i=1} D={di​}i=1K​其中 M M M為樣本數， y i ∈ R N \mathbf{y}_{i} \in R^{N} yi​∈RN為單個樣本，是一個 N N N維的特征向量， M M M個單個樣本按列組成樣本集矩陣， d i ∈ R N \mathbf{d}_i \in R^{N} di​∈RN為基向量，次元也是 N N N維， K K K個基向量按列組成字典矩陣，樣本集中的任意樣本都可以根據字典求得起對應的稀疏表達形式。當 K > N K > N K>N時為超完備字典，即我們字典學習的通常情況。

字典學習的目标是學習一個字典矩陣 D \mathbf{D} D，使得 Y ≈ D ∗ X \mathbf{Y} \approx \mathbf{D} * \mathbf{X} Y≈D∗X同時要 X X X盡可能稀疏。這裡我們注意 D ∈ R N × K \mathbf{D} \in \mathbf{R}^{N \times K} D∈RN×K， Y ∈ R N × M \mathbf{Y} \in \mathbf{R}^{N \times M} Y∈RN×M，是以 X ∈ R K × M \mathbf{X} \in \mathbf{R}^{K \times M} X∈RK×M。如下圖所示

D \mathbf{D} D矩陣中不同顔色代表不同的基向量， X \mathbf{X} X矩陣中不同的顔色代表對應基向量的系數， Y \mathbf{Y} Y矩陣中某一樣本（灰色塊）就是由 D \mathbf{D} D矩陣的基向量通過 X \mathbf{X} X矩陣中黑框中的系數進行線性組合獲得的，由于是稀疏組合，是以中間有很多系數會是零，我們用白色塊表示。

這個問題用數學問題描述為 min ⁡ D , X ∥ Y − D X ∥ F 2 ,  s.t.  ∀ i , ∥ x i ∥ 0 ≤ T 0 \min _{\mathbf{D}, \mathbf{X}}\|\mathbf{Y}-\mathbf{D} \mathbf{X}\|_{F}^{2}, \quad \text { s.t. } \forall i,\left\|\mathbf{x}_{i}\right\|_{0} \leq T_{0} D,Xmin​∥Y−DX∥F2​, s.t. ∀i,∥xi​∥0​≤T0​或者 min ⁡ D , x ∑ i ∥ x i ∥ 0 ,  s.t.  min ⁡ D , x ∥ Y − D X ∥ F 2 ≤ ϵ \min _{\mathbf{D}, \mathbf{x}} \sum_{i}\left\|\mathbf{x}_{i}\right\|_{0}, \quad \text { s.t. } \min _{\mathbf{D}, \mathbf{x}}\|\mathbf{Y}-\mathbf{D} \mathbf{X}\|_{F}^{2} \leq \epsilon D,xmin​i∑​∥xi​∥0​, s.t. D,xmin​∥Y−DX∥F2​≤ϵ其中 x i \mathbf{x}_i xi​為稀疏矩陣 X \mathbf{X} X的行向量，代表字典矩陣的系數，和上面彩圖中是對應的。這裡值得注意的是 ∥ x i ∥ 0 \left\|\mathbf{x}_{i}\right\|_{0} ∥xi​∥0​是零階範數，它表示向量中不為0的數的個數。

上面兩個公式是帶有限制的優化問題，可以利用拉格朗日乘子法将其轉化為無限制優化問題 min ⁡ D , x ∥ Y − D X ∥ F 2 + λ ∥ x i ∥ 1 \min _{\mathbf{D}, \mathbf{x}}\|\mathbf{Y}-\mathbf{D} \mathbf{X}\|_{F}^{2}+\lambda\left\|\mathbf{x}_{i}\right\|_{1} D,xmin​∥Y−DX∥F2​+λ∥xi​∥1​其中主要是将 ∥ x i ∥ 0 \left\|\mathbf{x}_{i}\right\|_{0} ∥xi​∥0​用 ∥ x i ∥ 1 \left\|\mathbf{x}_{i}\right\|_{1} ∥xi​∥1​代替，這樣更加便于求解。這裡要優化 D \mathbf{D} D， X \mathbf{X} X兩個變量，而樣本集 Y \mathbf{Y} Y是已知的，我們做法是固定一個變量，優化另一個，交替進行。下面分開來說：

已知 D \mathbf{D} D，優化 X \mathbf{X} X：這個問題有許多經典算法可以進行求解，例如Lasso，OMP，我們之後再進行補充，K-SVD中比較有特點是下面這種情況。

已知 X \mathbf{X} X，優化 D \mathbf{D} D：我們需要逐列來更新字典，記 d k \mathbf{d}_k dk​為字典 D \mathbf{D} D的第 k k k列，其實就是某個基向量， x T k \mathbf{x}_{T}^{k} xTk​為稀疏矩陣 X \mathbf{X} X的第 k k k行，按照本文之前的定義，其實就是該基向量不同樣本對應的系數。那麼有：

∥ Y − D X ∥ F 2 = ∥ Y − ∑ j = 1 K d j x T j ∥ F 2 = ∥ ( Y − ∑ j ≠ k d j x T j ) − d k x T k ∥ F 2 = ∥ E k − d k x T k ∥ F 2 \begin{aligned}\|\mathbf{Y}-\mathbf{D} \mathbf{X}\|_{F}^{2} &=\left\|\mathbf{Y}-\sum_{j=1}^{K} \mathbf{d}_{j} \mathbf{x}_{T}^{j}\right\|_{F}^{2} \\ &=\left\|\left(\mathbf{Y}-\sum_{j \neq k} \mathbf{d}_{j} \mathbf{x}_{T}^{j}\right)-\mathbf{d}_{k} \mathbf{x}_{T}^{k}\right\|_{F}^{2} \\ &=\left\|\mathbf{E}_{k}-\mathbf{d}_{k} \mathbf{x}_{T}^{k}\right\|_{F}^{2} \end{aligned} ∥Y−DX∥F2​​=∥∥∥∥∥​Y−j=1∑K​dj​xTj​∥∥∥∥∥​F2​=∥∥∥∥∥∥​⎝⎛​Y−j​=k∑​dj​xTj​⎠⎞​−dk​xTk​∥∥∥∥∥∥​F2​=∥∥​Ek​−dk​xTk​∥∥​F2​​上式中定義殘差 E k = Y − ∑ j ≠ k d j x T j \mathbf{E}_{k}=\mathbf{Y}-\sum_{j \neq k} \mathbf{d}_{j} \mathbf{x}_{T}^{j} Ek​=Y−∑j​=k​dj​xTj​就是除了目前更新基向量外其它基向量組合與樣本集的誤差，我們的目的就是讓該殘差項最小，是以有 min ⁡ d k , x T k ∥ E k − d k x T k ∥ F 2 \min _{\mathbf{d}_{k}, \mathbf{x}_{T}^{k}}\left\|\mathbf{E}_{k}-\mathbf{d}_{k} \mathbf{x}_{T}^{k}\right\|_{F}^{2} dk​,xTk​min​∥∥​Ek​−dk​xTk​∥∥​F2​這個目标函數如何求呢？這裡采用的就是SVD方法，如下： E k ′ = U Σ V T \mathbf{E}_{k}^{\prime}=U \Sigma V^{T} Ek′​=UΣVT其中 Σ \Sigma Σ矩陣中的奇異值是從大到小排列的，我們要取最大奇異值對應的向量作為最優的 d k , x T k \mathbf{d}_{k}, \mathbf{x}_{T}^{k} dk​,xTk​，那麼就是取矩陣 U U U的第一個列向量 u 1 = U ( ⋅ , 1 ) \mathbf{u}_{1}=U(\cdot, 1) u1​=U(⋅,1)作為 d k \mathbf{d}_{k} dk​，取矩陣 V V V的第一個行向量與第一個奇異值的乘積作為 x T ′ k \mathbf{x}_{T}^{\prime k} xT′k​，将其更新到 x T k \mathbf{x}_{T}^{k} xTk​，這裡值得注意的是，這裡更新 x T k \mathbf{x}_{T}^{k} xTk​并不能保證其稀疏性，要保證稀疏性還是需要回到上一種情況使用OMP算法進行稀疏編碼

這裡為什麼可以通過SVD算法求得最優解呢？關于這一點之後再補充，那麼以上就完成了字典學習的過程，那麼來看一段代碼印證下KSVD的流程。

2. python代碼

這段代碼是我從部落格K-SVD字典學習及其實作（Python）扒過來的，該了一下opencv的接口，這裡要注意的是，如果直接按照我這個代碼把整張圖檔輸入的話其實降噪效果不理想的，為什麼呢？

#!/usr/bin/python
# -*- coding: UTF-8 -*-
import numpy as np
from sklearn import linear_model
import scipy.misc
from matplotlib import pyplot as plt
import cv2


class KSVD(object):
    def __init__(self, n_components, max_iter=30, tol=5000,
                 n_nonzero_coefs=None):
        """
        稀疏模型Y = DX，Y為樣本矩陣，使用KSVD動态更新字典矩陣D和稀疏矩陣X
        :param n_components: 字典所含原子個數（字典的列數）
        :param max_iter: 最大疊代次數
        :param tol: 稀疏表示結果的容差
        :param n_nonzero_coefs: 稀疏度
        """
        self.dictionary = None
        self.sparsecode = None
        self.max_iter = max_iter
        self.tol = tol
        self.n_components = n_components
        self.n_nonzero_coefs = n_nonzero_coefs

    def _initialize(self, y):
        """
        初始化字典矩陣
        """
        u, s, v = np.linalg.svd(y)
        self.dictionary = u[:, :self.n_components]
        print(self.dictionary.shape)

    def _update_dict(self, y, d, x):
        """
        使用KSVD更新字典的過程
        """
        for i in range(self.n_components):
            index = np.nonzero(x[i, :])[0]
            if len(index) == 0:
                continue

            d[:, i] = 0
            r = (y - np.dot(d, x))[:, index]
            u, s, v = np.linalg.svd(r, full_matrices=False)
            d[:, i] = u[:, 0].T
            x[i, index] = s[0] * v[0, :]
        return d, x

    def fit(self, y):
        """
        KSVD疊代過程
        """
        self._initialize(y)
        for i in range(self.max_iter):
            x = linear_model.orthogonal_mp(self.dictionary, y, n_nonzero_coefs=self.n_nonzero_coefs)
            e = np.linalg.norm(y - np.dot(self.dictionary, x))
            if e < self.tol:
                break
            self._update_dict(y, self.dictionary, x)

        self.sparsecode = linear_model.orthogonal_mp(self.dictionary, y, n_nonzero_coefs=self.n_nonzero_coefs)
        return self.dictionary, self.sparsecode


if __name__ == '__main__':
    im_ascent = cv2.imread("./input.png", 0).astype(np.float)
    print(im_ascent.shape)
    ksvd = KSVD(300)
    dictionary, sparsecode = ksvd.fit(im_ascent)

    output = dictionary.dot(sparsecode)
    output = np.clip(output, 0, 255)
    cv2.imwrite("./output.png", output.astype(np.uint8))
           

從代碼中可以看到，我們輸入的樣本集 Y \mathbf{Y} Y矩陣就是整張圖檔，這樣的話，單個樣本就是圖檔的每一列的資料，這樣合理嗎？當然不合理，我覺得這就是造成僅僅通過上面這段代碼降噪效果不理想的原因，合理的做法應該是從圖像中提取patch，将patch轉換成一列一列的樣本資料，組成樣本集矩陣，然後再利用K-SVD算法從樣本集矩陣學習字典和稀疏矩陣。感興趣的同學可以根據該思路嘗試改下代碼，或者參考代碼nel215/ksvd

3. 結論

1. KSVD隻是衆多基于稀疏表達降噪的方法中最經典的一種，類似的還有LSSC，NCSR等，KSVD算法是2006年提出的，是以效果有限，對這個方向感興趣的同學可以多關注後續的一些paper
2. 在Image Denoising Via Sparse and Redundant Representations Over Learned Dictionaries中展示了KSVD算法的效果，這篇論文中描述的KSVD算法就是基于patch實作的

   

   其中，Global Trained Dictionary指的是從一堆無噪聲的圖檔中抽取的patch學習的字典，Adaptive Dictionary指的是從原噪聲圖像patch中學習的字典，從結果看，Adaptive Dictionary是要優于Global Trained Dictionary的

這篇部落格沒有做太多實驗去驗證明驗效果，有問題歡迎交流

此外，這裡我寫一個各種算法的總結目錄圖像降噪算法——圖像降噪算法總結，對圖像降噪算法感興趣的同學歡迎參考