A Dual Sparse Decomposition Method for Image Denoising

TL;DR

This paper introduces a dual sparse decomposition technique for image denoising, especially effective under strong noise conditions, by decomposing dictionaries based on atom occurrence frequency, leading to superior denoising performance.

Contribution

It presents a novel dual sparse decomposition approach that enhances image denoising by utilizing a new criterion for sub-dictionary selection based on atom frequency.

Findings

Outperforms state-of-the-art denoising methods in PSNR and SSIM.

Improves subjective visual quality of denoised images.

Effective under strong noise conditions.

Abstract

This article addresses the image denoising problem in the situations of strong noise. We propose a dual sparse decomposition method. This method makes a sub-dictionary decomposition on the over-complete dictionary in the sparse decomposition. The sub-dictionary decomposition makes use of a novel criterion based on the occurrence frequency of atoms of the over-complete dictionary over the data set. The experimental results demonstrate that the dual-sparse-decomposition method surpasses state-of-art denoising performance in terms of both peak-signal-to-noise ratio and structural-similarity-index-metric, and also at subjective visual quality.

Tables

Table 1. TABLE I: Denoising Algorithm Based on Dual Sparse Decomposition

Input: Image data $𝐗={\left\{{𝐱}_m\right\}}_{m=1}^M$ (Equ.(1)

Grouping: For patch ${𝐱}_i$, form group ${𝐳}_i$ according Equs. (14) or (15)

Dual sparse decomposition: For each group ${𝐳}_i,\forall i$ do

- Sparse decomposition: ${𝐳}_i\approxeq \mathrm{𝐃𝐀}$ by solving Equ.(2)

- Subspace decomposition: $𝐃={𝐃}^{\left(𝐒\right)}+{𝐃}^{\left(𝐍\right)}$ by Equs.(5, (8)-(11)

- Linear reconstruction on ${𝐃}^{\left(𝐒\right)}$: ${\widehat{𝐒}}_i={𝐃}_P^{\left(𝐒\right)}.{𝐀}_P^{\left(𝐒\right)}$ by Equs. (12)-(13)

Aggregate: to form denoised image $weight-average\{{\widehat{𝐒}}_i,\forall i\}\Rightarrow \widehat{𝐒}$

Output: Denoised image $\widehat{𝐒}$.

Equations

$${\mathbf{X}}_{N\times M}=\{{\mathbf{x}}_m{\}}_{m=1}^M$$

$$\{\mathbf{D},{\mathbf{\alpha }}_m\}=\underset{\mathbf{D},{\mathbf{\alpha }}_m}{\mathrm{argmin}}\parallel {\mathbf{\alpha }}_m{\parallel }_0+\parallel \mathbf{D}{\mathbf{\alpha }}_m-{\mathbf{x}}_m{\parallel }_2^2\le \epsilon ,$$

$$1\le m\le M$$

$${\mathbf{X}}_{N\times M}\approxeq {\mathbf{D}}_{N\times K}{\mathbf{A}}_{K\times M}$$

$${\mathbf{A}}_{K\times M}=\left[{\mathbf{\alpha }}_1\cdots {\mathbf{\alpha }}_k\cdots {\mathbf{\alpha }}_K\right]$$

$$\mathbf{S}=\left[{\mathbf{d}}_1{\mathbf{d}}_2\cdots {\mathbf{d}}_k\cdots {\mathbf{d}}_K\right]$$

$$\mathbf{A}$$

$$=\left[\begin{array}{c}{\alpha }_1(1){\alpha }_2(1)\cdots {\alpha }_M(1)\\ {\alpha }_1(2){\alpha }_2(2)\cdots {\alpha }_M(2)\\ ⋮⋮\ddots ⋮\\ {\alpha }_1(K){\alpha }_2(K)\cdots {\alpha }_M(K)\end{array}\right]=\left[\begin{array}{c}{\mathbf{\beta }}_1\\ {\mathbf{\beta }}_2\\ ⋮\\ {\mathbf{\beta }}_K\end{array}\right]$$

$${\mathbf{\beta }}_k=[{\alpha }_1(k){\alpha }_2(k)\dots {\alpha }_M(k)]\in {\mathbb{R}}^{1\times M}$$

$${\mathbf{A}}_{K\times M}={\left[{\mathbf{\beta }}_1^T\cdots {\mathbf{\beta }}_k^T\cdots {\mathbf{\beta }}_K^T\right]}^T$$

$${\mathbf{S}}_{N\times M}=\left[{\mathbf{d}}_1\cdots {\mathbf{d}}_k\cdots {\mathbf{d}}_K\right].{\left[{\mathbf{\beta }}_1^T\cdots {\mathbf{\beta }}_k^T\cdots {\mathbf{\beta }}_K^T\right]}^T$$

$$f_k\triangleq Frequency({\mathbf{d}}_k|\mathbf{X})=\Vert {\mathbf{\beta }}_k{\Vert }_0$$

$$\tilde{\mathbf{A}}\triangleq$$

$$\text{s.t.}\parallel {\mathbf{\beta }}_1^{\prime }{\parallel }_0\ge \parallel {\mathbf{\beta }}_2^{\prime }{\parallel }_0\ge \cdots \ge \parallel {\mathbf{\beta }}_K^{\prime }{\parallel }_0$$

$$\mathbf{D}\text{xLongrightarrow}sort\tilde{\mathbf{D}}=[{\text{d}}_1^{\prime },\cdots ,{\text{d}}_k^{\prime },\cdots ,{\text{d}}_K^{\prime }]$$

$${\mathbf{S}}_{N\times M}$$

$$=\left[{\mathbf{d}}_1^{\prime },\cdots ,{\mathbf{d}}_k^{\prime },\cdots ,{\mathbf{d}}_K^{\prime }\right].{\left[{\mathbf{\beta }}_1^{{}^{\prime }T}\cdots {\mathbf{\beta }}_k^{{}^{\prime }T}\cdots {\mathbf{\beta }}_K^{{}^{\prime }T}\right]}^T$$

$${\mathbf{D}}_P^{(S)}=span\{{\mathbf{d}}_1^{\prime },{\mathbf{d}}_2^{\prime },\cdots ,{\mathbf{d}}_P^{\prime }\}$$

$${\mathbf{D}}_{K-P}^{(N)}=span\{{\mathbf{d}}_{P+1}^{\prime },{\mathbf{d}}_{P+2}^{\prime },\cdots ,{\mathbf{d}}_K^{\prime }\}$$

$$P=\underset{k}{arg\; max\; Hist}(\Vert {\mathbf{\beta }}_k^{\prime }{\Vert }_0)$$

$$\hat{\mathbf{S}}$$

$$=\left[{\mathbf{d}}_1^{\prime },\cdots ,{\mathbf{d}}_k^{\prime },\cdots ,{\mathbf{d}}_P^{\prime }\right].{\left[{\mathbf{\beta }}_1^{{}^{\prime }T}\cdots {\mathbf{\beta }}_k^{{}^{\prime }T}\cdots {\mathbf{\beta }}_P^{{}^{\prime }T}\right]}^T$$

$$\text{SNR}=20\cdot {\log }_{10}\left[\max \{\mathbf{S}(i,j)\}\right]-10\cdot {\log }_{10}\left[\text{MSE}\right]$$

$$\text{MSE}=\frac{1}{IJ}{\sum }_{i=0}^{I-1}{\sum }_{j=0}^{J-1}{\left[\mathbf{S}(i,j)-\hat{\mathbf{S}}(i,j)\right]}^2$$

$$\text{SSIM}(\mathbf{S},\hat{\mathbf{S}})=\frac{(2u_{\mathbf{S}}{u}_{\hat{\mathbf{S}}}+c_1)(2{\sigma }_{\mathbf{S}\hat{\mathbf{S}}}+c_2)}{(u_{\mathbf{S}}^2+u_{\hat{\mathbf{S}}}^2+c_1)({\sigma }_{\mathbf{S}}^2+{\sigma }_{\hat{\mathbf{S}}}^2+c_2)}$$

$$\gamma ({\mathbf{x}}_i,{\mathbf{x}}_j)=\Vert {\mathbf{x}}_i-{\mathbf{x}}_j{\Vert }_0^2$$

$$\gamma ({\mathbf{x}}_i,{\mathbf{x}}_j)=(2L-1)\sum _k\log \sqrt{\frac{{\mathbf{y}}_i(k)}{{\mathbf{y}}_j(k)}}+\sqrt{\frac{{\mathbf{y}}_j(k)}{{\mathbf{y}}_i(k)}}$$

Full text

I Introduction

Two main issues are involved in the denoising problem. One is the filtering technique by signal analysis to identify the information underlying the noisy data. The other is grouping technique by clustering technique to provide homogeneous signals for filtering.

Almost all filtering techniques assume that the involved signal should be homogeneous. Therefore, a grouping procedure is generally required before filtering. Many edge detection and image segmentation techniques [1] are used in image denoising. Recently, a nonlocal self-similarity method [2] provides a potential breakthrough for data grouping, which is adopted in this paper.

The filtering technique is developed in the past 50 years or so from many diverse points of view, statistical estimation method, such as Viener filter, adaptive filter, etc. [3]; transform-domain method, such as Principal Components Analysis [4], wavelet shrinkage [5], etc., and so on. The underlying assumption of these filtering methods is that information in the noisy data has a property of energy concentration in a small linear subspace of the overall space of possible data vectors, whereas additive noise is typically distributed through the larger space isotropically.

However, in many practical cases, some components with low energy might actually be important because they carry information relative to the signal details. On the contrary, when dealing with noise with non-Gaussian statistics, it may happen that some noise components may have higher energies. Consequently, a major difficulty of filtering is to separate the information details from noise. A way to deal with this problem is cooperative filtering technique, such as Turbo iterative filter [6].

In recent years, sparse coding has attracted significant interest in the field of signal denoising [7] upon an over-complete dictionary. A sparse representation is a signal decomposition on a very small set of components (called atoms) which are adapted to the observational data. The sparse-decomposition based denoising is much better at the trade-off between the preservation of details and the suppression of noise. However, the sparse decomposition is adapted to noisy data so that separating details from noise still is at issue.

In this paper, we propose a dual sparse decomposition method for filtering. The first decomposition is to make an over-complete dictionary to reject some noises which really distributed through the larger space isotropically but to preserve the information details as much as possible. The second decomposition is to identify principal atoms to form a sub-dictionary which preserve well the weak information details and simultaneously suppress strong noises.

This article is organized as follows: Section 2 analyzes some limitations of the classical sparse decomposition for denoising. Section 3 presents the principle of the proposed dual sparse decomposition. Section 4 shows some experimental results and comparisons with state-of-art image denoising methods. Finally, we draw the conclusion in Section 5.

II Sparse Decomposition for Denoising

We start with a brief description of the classical sparse decomposition and analyze their limitations for denoising.

The sparse decomposition of $M$ observations $\{\mathbf{x}_{m}\in\mathbb{R}^{N}\}_{m=1}^{M}$ based on a dictionary $\mathbf{D}=\{\mathbf{d}_{k}\}_{k=1}^{K}\in\mathbb{R}^{N\times K}$. When $K>N$, the dictionary is said over-complete. $\mathbf{d}_{k}\in\mathbb{R}^{N}$ is a basis vector, also called an atom of the dictionary. They are not necessarily independent. With observational data set:

[TABLE]

a dictionary and the coefficients can be the solution of the following equation [8]:

[TABLE]

where $\|\bullet\|_{2}$ denotes $\ell^{2}$-norm and $\|\bullet\|_{0}$ denotes $\ell^{0}$-norm. In equation (2), $\bm{\alpha}_{m}=\left[\alpha_{m}(1)\;\alpha_{m}(2)\;\dots\;\alpha_{m}({K})\right]^{T}\in\mathbb{R}^{K\times 1}$ is the sparse code of the observation $\mathbf{x}_{m}$. The allowed error tolerance $\varepsilon$ can be chosen according to the standard deviation of the noise. The sparse decomposition can be written in matrix form as:

[TABLE]

where the matrix $\mathbf{A}$ of size $K\times M$ is composed of $M$ sparse column vectors $\bm{\alpha}_{m}$:

[TABLE]

An estimate of the underlying signal $\mathbf{S}$ embedded in the observed data set $\mathbf{X}$ would be:

[TABLE]

The over-complete dictionary $\mathbf{D}$ in sparse decomposition can effectively capture the information patterns and reject white Gaussian noise patterns. However, we note that the learning algorithm for dictionary $\mathbf{D}$ by equation (2) would fall into a dilemma of preserving weak derails and suppressing noise. On one hand, in order to suppress noise, the allowed error tolerance $\varepsilon$ in equation (2) should be small enough. As a result, certain weak details would be lost. On the other hand, in order to capture weak details, $\varepsilon$ cannot be too small. Otherwise some atoms would be so noisy that degrade the denoising performance. Fig. 1 shows an example to show this situation. Taking a noisy image degraded by white noise with standard deviation $\sigma=35$ (Fig. 1a), we make two different dictionaries $\textbf{D}_{H}$ and $\textbf{D}_{L}$ (Fig. 1b) by solving equation (2) with $\varepsilon=40$ and $\varepsilon=35$ respectively. We got two different retrieved images $\textbf{S}_{H}$ and $\textbf{S}_{L}$ (Fig. 1b) respectively by equation (4). Intuitively, the noise is well suppressed in $\textbf{S}_{H}$ but some information details are lost. On the contrary, more details are reserved in $\textbf{S}_{L}$ but it is rather noisy.

Considering the above limitation of sparse-decomposition-based denoising, our idea of dual sparse decomposition is to make a two-stop sparse decomposition: The first step is to make an over-complete dictionary by learning from the observational data with a lower allowed error tolerance $\varepsilon$ according to equation (2). Thereby, the obtained dictionary $\textbf{D}_{L}$ can capture more information details although it contains some noisy atoms. The second step is to make a sub-dictionary decomposition on $\textbf{D}_{L}$ to reject some atoms too noisy.

III Sparse Subspace Decomposition

To get the sub-dictionary, we introduce a novel criterion to the sparse subspace decomposition of a learned dictionary and a corresponding index of significance of the atoms.

III-A Occurrence Frequency of Atom

Atoms $\{\mathbf{d}_{k}\}_{k=1}^{K}$ in the sparse decomposition are prototypes of signal segments. This property allows us to take the atoms as a signal pattern. Thereupon, some important features of the signal pattern could be considered as a criterion to identify significant atoms. We note a common knowledge about the regularity of signal: A signal pattern must occur in meaningful signals with higher frequency even with a lower energy, such as the geometrical regularity of image structures like edges and textures. On the contrary, a noise pattern would hardly be reproduced in observed data even with a higher energy. Therefore, we propose to take the frequency of atoms appeared in the data set as the criterion to identify principal atoms [9]. In fact, the frequency of atoms is a good description of the signal texture [10].

We intend to find out a measurement of the frequency of atom from the sparse codes. Coefficient matrix $\mathbf{A}$ in the sparse representation (equation (3)) is composed by $M$ sparse column vectors $\bm{\alpha}_{m}$. Let us consider the row vectors $\{\bm{\beta}_{k}\}_{k=1}^{K}$ of coefficient matrix $\mathbf{A}$ :

[TABLE]

where

[TABLE]

Note that the row vector $\bm{\beta}_{k}$ is not necessarily sparse.

Thus, the coefficient matrix $\mathbf{A}$ can be written by $K$ row vectors as:

[TABLE]

Then equation (4) can be expressed by reordered dictionary and its coefficient as:

[TABLE]

Denoting $\|\bm{\beta}_{k}\|_{0}$ the $\ell^{0}$ zero pseudo-norm of $\bm{\beta}_{k}$, we find that $\|\bm{\beta}_{k}\|_{0}$ is just the number of occurrences of atom $\mathbf{d}_{k}$ over the data set $\{\mathbf{x}_{m}\}_{m=1}^{M}$. We can define the frequency of the atom $\mathbf{d}_{k}$ as $f_{k}$:

[TABLE]

III-B Subspace decomposition on Over-complete Dictionary

Taking vectors $\{\bm{\beta}_{k}\}_{k=1}^{K}$ from equation (5), we calculate their $\ell^{0}$-norms $\{\|\bm{\beta}_{k}\|_{0}\}_{k=1}^{K}$ and rank them in descending order:

[TABLE]

Corresponding to the order of $\{\bm{\beta}^{\prime}_{k}\}_{k=1}^{K}$, the reordered dictionary is written as:

[TABLE]

Equation (6) becomes as:

[TABLE]

Then, the first $P$ atoms of $\tilde{\mathbf{D}}$ can span a principal subspace $\mathbf{D}_{P}^{(S)}$ and the remaining atoms span a noise subspace $\mathbf{D}_{K-P}^{(N)}$ as:

[TABLE]

In practical application, $P$ is the threshold of $f_{k}$ to separate the principal sub-dictionary from the noise sub-dictionary. We set the maximum point of the histogram of $\{\|\bm{\beta}_{k}\|_{0}\}_{k=1}^{K}$ to $P$ as:

[TABLE]

An estimate of the underlying signal $\hat{\mathbf{S}}$ embedded in the observed data set $\mathbf{X}$ can be obtained on the principal sub-dictionary $\mathbf{D}^{(\mathbf{S})}$ simply by linear combination:

[TABLE]

Note that $P\ll K$.

We show an example of the proposed dual sparse decomposition in Fig. 1(c). The learned over-complete dictionary $\mathbf{D}$ is decomposed into a principal sub-dictionary $\mathbf{D}^{(S)}$ and a noise sub-dictionary $\mathbf{D}^{(N)}$ under the atom’s frequency criterion. The retrieved image $\hat{\mathbf{S}}^{(S)}$ by the dual sparse decomposition method has a super performance at preserving fine details and at suppressing strong noise. We note that the residual image $\mathbf{X}^{(N)}$ on the noise sub-dictionary $\mathbf{D}^{(N)}$ contains some information but very noisy. This is because the atoms of the over-complete dictionary are not independent. The information in the residue image $\mathbf{X}^{(N)}$ is also in existence in $\hat{\mathbf{S}}^{(S)}$.

III-C Application to Filtering

A major difficulty of filtering is to suppress noise Gaussian or non-Gaussian and to preserve information details simultaneously. We use the peak signal-to-noise ratio (PSNR) to assess the noise removal performance:

[TABLE]

and the structural similarity index metric (SSIM) between denoised image $\hat{\mathbf{S}}$ and the pure one $\mathbf{S}$ to evaluate the preserving details performance:

[TABLE]

where $u_{x}$ is the average of $x$, $\sigma_{x}^{2}$ is the variance of $x$, $\sigma_{xy}$ is the covariance of $x$ and $y$, and $c_{1}$ and $c_{2}$ are small variables to stabilize the division with weak denominator.

From the example shown in Fig. 1, the retrieved image $\mathbf{S}_{H}$ actually by the K-SVD filter [11] with the classical sparse decomposition has a high performance with $\textit{PSNR}=34.25$ and $\textit{SSIM}=0.82$ but some information details are obviously lost. On the contrary, the retrieved image $\mathbf{S}_{L}$ is noisier with $\textit{PSNR}=29.62$ and $\textit{SSIM}=0.78$ but more information details are reserved. Making a dictionary decomposition on $\mathbf{D}_{L}$ noisier but with more details, the retrieved image $\hat{\mathbf{S}}^{(S)}$ based on the principal sub-dictionary $\mathbf{D}^{(S)}$ has a higher performance with $\textit{PSNR}=35.82$ and $\textit{SSIM}=0.86$.

Fig. 2 shows an image filtering result based on the proposed dual sparse decomposition and a comparison with K-SVD algorithm. From the results, the dual sparse decomposition method outperforms K-SVD method by about $1dB$ in PSNR and by about $1\%$ in SSIM. In terms of subjective visual quality, we can see that the corner of mouth and the nasolabial fold with weak intensities are much better recovered by the dual sparse decomposition method.

Fig. 3 shows the despeckling results of simulated one-look SAR scenario with a fragment of Barbara image. From the result by a probabilistic patch based (PPB) filter [12] which can cope with non-Gaussian noise, we can see that PPB can well remove speckle noise. However, it also removes low-intensity details. The dual sparse decomposition method shows advantages at preserving fine details and at suppressing strong noise.

IV Application to Denoising

In practical applications, our images are generally with spatial complicated scene. On the other hand, the used filtering techniques are generally suitable to homogeneous images. For image denoising based on the sparse decomposition, the hypotheses of signal sparsity and component reproducibility mean also the condition of homogeneity. In order to make the involved signal homogeneous, we select homogeneous pixels before filtering by a self-similarity measure [2] $\gamma$. In applications of image denoising, $\gamma$ can be specified as Euclidean distance between the reference patch $\mathbf{x}_{i}$ and a given patch $\mathbf{x}_{j}$ as:

[TABLE]

The smaller $\gamma$ is, the more similar between $\mathbf{x}_{i}$ and $\mathbf{x}_{j}$ is. This self-similarity matches well the property of highly repetitive structures of images.

In applications of image despeckling, $\gamma$ becomes the probabilistic patch-based similarity proposed by [12] as:

[TABLE]

where $\mathbf{y}_{i}=\exp(\mathbf{x}_{i})$ and $L$ the equivalent number of looks.

For a given reference patch $\mathbf{x}_{i}$, we make grouping stacks with its $\Gamma$-most similar patches to form a group of data $\mathbf{z}_{i}$. In our experiments, we take $\Gamma=90$. Then we apply a filtering algorithm to each of the data groups $\mathbf{z}_{i},\forall i$. Our denoising algorithm is presented in Table I:

To compare with the state-of-art denoising algorithm, we take the BM3D method [13], one of the best method nowadays for image denoising. In the BM3D method, a block-matching grouping is also used before filtering. In the experiments, the used dictionaries $\mathbf{D}$s are of size are of size $64\times 256$ ($K=256$ atoms), designed to handle image patches $\mathbf{x}_{m}$ of size $N=64=8\times 8$ pixels.

Fig. 4 shows the results of denoising an image with a strong additive zero-mean white Gaussian noise and their performances of the dual-sparse-decomposition method and the BM3D method. Fig. 5 shows the results of despeckling a simulated one-look SAR image with non-Gaussian noise and their performances of the dual-sparse-decomposition method and the SAR-BM3D method [14]. The experimental results demonstrate some advantage of the dual-sparse-decomposition method at preserving fine details and at suppressing speckle noise, also with a better subjective visual quality over the BM3D method.

V Conclusion

This work present a new signal analysis method by a proposed dual sparse decomposition, leading to state-of-the-art performance for image denoising. The proposed method introduces a sub-dictionary decomposition on an over-complete dictionary learned under a lower allowed-error-tolerance. The principal sub-dictionary is identified under a novel criterion based on the occurrence frequency of atoms. The experimental results have demonstrated that the proposed dual-sparse-decomposition-based denoising method has some advantages both at preserving information details and at suppressing strong noise, as well as provides retrieved image with better subjective visual quality.

It is perfectly possible to straightforward extension the proposed dual-sparse-decomposition to application of feature extraction, inverse problems, or machine learning.

Acknowledgment

This work was supported by the National Natural Science Foundation of China (Grant No. 60872131).

The idea of the dual Sparse decomposition arises through a lot of deep discussions with Professor Henri Ma$\hat{i}$tre at Telecom-ParisTech. The mathematical expressions in this paper are corrected by Professor Didier Le Ruyet at CNAM-Paris.