# Dictionary-Based Image Denoising by Fused-Lasso Atom Selection

Authors: Ao Li; Hayaru Shouno
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

We proposed an efficient image denoising scheme by fused lasso with dictionary learning. The scheme has two important contributions. The first one is that we learned the patch-based adaptive dictionary by principal component analysis (PCA) with clustering the image into many subsets, which can better preserve the local geometric structure. The second one is that we coded the patches in each subset by fused lasso with the clustering learned dictionary and proposed an iterative Split Bregman to solve it rapidly. We present the capabilities with several experiments. The results show that the proposed scheme is competitive to some excellent denoising algorithms.


---

## Body

## 1. Introduction

As an essential low-level image processing procedure, image denoising is studied extensively, which is also a classical inverse problem. The general observation with additive noise is modeled as:(1)y=x+n,
where y is the noisy observation and x and n present the original image and white Gaussian noise, respectively.With the degradation model in (1), many denoising algorithms are proposed during the past decades. From the earlier spatial and frequency domain filters [1, 2] to the lately developed wavelet and beyond wavelet based shrinkage methods [3–5]. Due to the the fact that the wavelet is not enough to present the local structure in images, an effective dictionary-learning algorithm called KSVD was proposed in [6] and achieved the good effect in image denoising task. In this method, a highly patch-based overcompleted dictionary is learned firstly, and then denosing is implemented by solving the sparse coding of each patch with assuming that the patches are sparsely representable under the learned dictionary. Foi et al. proposed the pointwise shape-adaptive DCT to be applied to the patch and its neighborhoods, which achieved the very sparse representation for noisy image and got the effective denoising results [7]. To lower the complexity, an orthogonal dictionary learning method is proposed in [8], which intends to train a global dictionary by collecting samples from degradation image randomly. Though it achieves good performance in image restoration, some progress can still be made in learning better dictionaries to present the patches accurately. To better present the image, a two-stage denoising method with PCA is proposed in [9], which trained the dictionary by the local neighboring patches.More recently, a set of approaches with nonlocal techniques (NL) is used for removing noise. The idea of NL can be tracked to [10], where the similar pixels are searched to weighted average the filtered pixel. But the weight in NL is determined only by the intensity distance between patches, so it cannot guarantee to search the patches with similar local geometric structure under the strong noise. Zhang et al. proposed a novel nonlocal mean method with application to MRI denoising [11]. To make full use of nonlocal similarity, Mairal et al. proposed the nonlocal sparse model for image restoration [12]. Also motivited by the NL, a collaborative denoising scheme called BM3D is proposed by Dabov et al. [13], where they took the patch matching to search the similar patches and grouped them into a 3D cube matrix. Then, the algorithm carried out the 3D sparse transform, such as the 3D Wavelet or 3D Curvelet, to the cube matrix and then removed the noise with Wiener filtering in the transform domain. Another effective method addressed the denoising problem under the kernel regression (KR) framework proposed by Takeda et al. [14]. Many spatial classical denoising algorithms, such as the bilaterial filtering [15] and nonlocal mean [16, 17], can be seen as the special cases in KR with different constraints.In this paper, we proposed a novel scheme for image denoising based on clustering dictionary learning. Firstly, we clustered the patches with similar geometric structure by taking a weight function as the feature. Secondly, we learned the patch-based dictionary by principle component analysis for each cluster. Lastly, we coded these patches by fused lasso and developed an iterative Split Bregman to solve it rapidly.The rest of the paper is outlined as follows. In Section2, we briefly review the kernel regression framework as we want to choose the weight function in KR as the feature for clustering. And then, we talk about how to learn a suitable dictionary to better describe the patches in each cluster in Section 3. An iterative Split Bregman to solve the fused lasso is proposed in Section 4, which is used to code the patches under the dictionary learned in Section 3 rapidly. Section 5 shows several experiments results compared to the current excellent algorithms and Section 6 concludes the paper with a summary.


## 2. Clustering with Weight Function in KR

Kernel regression is well studied in statistical signal processing. Recently, the KR is used to address many image restoration tasks, such as the denoising, interpolation, and deblurring [18]. The kernel regression expression can be written as follows:
(2)x^i=∑j=1NKh(pi-pj)yj∑j=1NKh(pi-pj),Kh(pi-pj)=1hK(pi-pjh),
where N is the number of pixels in neighborhood. xi presents the ith pixel in image, whose location is denoted by pi. K(·) is the local kernel to measure the distance between the center pixel and its neighborhoods, and h is the smoothing parameter to control the penalization strength.By (2), we can obtain that the key technique in KR was how to determine a effective form of kernel function, which was studied in literatures [19, 20]. Among these methods, the steering kernel is distinguished by producing the local regression weights, where local gradient was taken into consideration to analyze the similarity between the pixels in a neighborhood. The weights in steering kernel can be expressed as:
(3)ωij=det⁡(Cj)2πh2μj2exp⁡{(pi-pj)TCj(pi-pj)2h2μj2},
where ωij denotes the structure similarity between the ith center pixel with the jth pixel in neighborhood. Cj is the covariance matrix with the gradient of the jth pixel. Furthermore, the whole kernel consists of all weights in neighborhood.Introduced in [14], the ωij can present the underlying local structure of patch centered in ith pixel. In addition, Takeda et al. pointed out that the different locations patched with different intensities but similar underlying structure will still produce the similar kernel. Generally, the clustering is implemented by the Euclidean measurement of intensity, such as the denoising algorithm in NLM. But, distinguished from the regression algorithm, what we want here is to learn a dictionary to describe the patches with similar geometric structure. That is, we do not require them to have similar intensity simultaneously. So, we can take ωij as the feature to measure the structure similarity among the patches. The significant distinction between the general Euclidean measurement and KR weighted function is that the latter can obtain the patches with similar structure, but not the patches with similar intensity. To this end, we can take the weights ωi=(ωi1,ωi2,…,ωiN) formed from streeing kernel as the feature and it will show some advantage on learning the clustering dictionary to better describe the local structure in each cluster. Also, the l1 norm can be used to measure the distance between the features, which show the anisotropy property. Next, we will talk about how to obtain the weights of each patch.Conveniently, the matrixCj is divided into three components [14] and can be reformulated as follows:
(4)Cj=γjUθjΛjUθjT,
where Uθj=[cos⁡θjsinθj-sinθjcos⁡θj] and Λj=[σj00σj-1].By (4), what we need is to determine the two variables θj and σj. To this end, we calculated the local gradient matrix Gj of jth pixel as follows:
(5)Gj=[⋮⋮gh(pk)gv(pk)⋮⋮],pk∈Wj(k=1,2,…,M),
where gh(·) and gv(·) present the horizontal and vertical gradient operator, respectively; Wj is an analysis window around the jth pixel and M is the number of pixels in the window. With the singular value decomposition (SVD), we can obtain
(6)Gj=UjSjVjT.WithVj=(v1,v2)T and Sj=diag⁡{s1,s2}, we can calculate the parameters in (4) as:
(7)θj=arctan(v1v2),σj=s1+ηs2+η,γj=s1s2+λM,
where η is the regulation parameter, which is used to prevent the ratio from being degenerate. λ is used to keep γj from being zero.We summarized the calculation of weights of each patch in Algorithm1.Algorithm 1
Initialization: patchyi centered in ith pixel, N is the number of pixels in yi, η=1, λ=0.01, μi=1, K=1;
(1) Calculate the local gradient matrixGK of Kth pixel in yi;
(2) Implement the SVD toGK and computed the parameters by (7);
(3) Compute the covariance matrixCK by (4);
(4) Compute theωiK by (3);
(5) IfK<N, K=K+1; otherwise, output the feature vector ωi=(ωi1,ωi2,…,ωiN);Then, we can cluster the ordering overlapping patches into subsetΩk (k is the cluster indicator) by K-Means with l1 norm to measure the similarity between the samples and clustering centers.


## 3. Adaptive PCA Dictionary Learning

Once the clusters are formed, we can learn a dictionary with the principal component analysis suited to each cluster independently. To this end, according to the general dictionary learning algorithm, we need to solve the following minimization:(8)(D^k,A^k)=argminDk,Ak{∥Y-k-DkAk∥F2},
where Y-k={y-i,i∈Ωk} is the centered samples matrix which satisfied Y-k=Yk-y-(k). Yk is the samples matrix, whose columns are the patches in the kth cluster, and also we can present it as Yk={yi,i∈Ωk}. y-(k) is the mean vector of Yk and ∥·∥F denotes the Frobenius norm.As the minimization in (8), the numerical method of alternate minimization is used to estimate the two variables. That is to say, we should estimate one variable while the other is fixed. Separately, we rewrite (8) as follows:
(9)(D^k,A^k)=argminDk,Ak{∑i=1#Ωk∥y-i-Dkaki∥22},
where aki is the ith column of Ak. When the dictionary is fixed, the expression of aki is given by
(10)aki=(DkTDk)-1DkTy-i.Then, assuming theAk is fixed and taking formulation (10) into (9), we can obtain
(11)D^k=argminDk,Ak{∑i=1#Ωk∥y-i-Dk(DkTDk)-1DkTy-i∥22}.The patches in the same cluster have similar structure to each other, so we do not require the dictionary to be redundant enough. So, to simplify the problem (11), we add the constraint DkTDT=I to dictionary and formulation (11) can be changed to
(12)min⁡Dk{∑i=1#Ωk∥y-i-DkDkTy-i∥22},s.t.DkTDk=I.The minimization problem in (12) can be approximated by finding the first m(k) principle components of centered matrix Y-k with PCA framework. And m(k) satisfied
(13)m(k)=max⁡(m),∑j=m+1#Ωksj2≥rNσ2,
where N is the number of pixels in each patch, r is a constant, and σ is the noise standard deviation. sj is the singular value of Y-k and satisfied sj≥st≥0, if j<t. The reason why we take (13) to choose m(k) is that it can discard the principle components that present the variance arisen by noise.Our above learning method can train the dictionary with lower complexity. In additional, to make it more effective and compact, we also show a selection scenario in (13), which can tradeoff between the essential signal preservation and noise reduction.We summarized the dictionary learning with PCA in Algorithm2.Algorithm 2
Input: Centered samples matrixY-k (#Ωk presents the length of Y-k);
Output: DictionaryDk for the cluster Ωk;
(1) Compute the decomposition of1/#ΩkY-kY-kT : 1/#ΩkY-kY-kT=UkΛkUkT;
(2) Choose them(k) by (13) with Λk(sj=Λk(j,j));
(3) Output theDk=(u1k,u2k,…,um(k)k), ujk is thejth column of matrix Uk;


## 4. Patch-Based Coding by Fused Lasso

With the preparation work in Sections2 and 3, in this paragraph, we will study the sparse coding with fused lasso [21]. The reason why we take the fused lasso for sparse coding is that it not only constrains the sparsity of coefficients but also enforced the differences between the neighboring coefficients, which led to show some advantage in recovering the texture in image [22, 23]. So, we can recover the noisy image by minimizing the cost function of fused lasso, which is expressed as follows:
(14)f(β)=12∥y-Dβ∥22+λ∑i=1p|βi|+η∑i=2n|βi-βi-1|,
where y is the sample patch in image denoising task, D can be seen as the dictionary, β=(β1,β2,…,βn) is the sparse coding of y, and λ and η are regularization parameters.The reference [21] points out that solving fused lasso is computationally demanding. So we developed a novel scenario to solve the fused-lasso problem in (14) based on Split Bregman, which is originally proposed in [24] and used to solve the l1 minimization successfully [25, 26]. Conveniently, we rewrite (14) as follows:
(15)min⁡β{12∥y-Dβ∥22+λ∥β∥1+η∥Lβ∥1},
where L is a matrix, whose each row only contains two nonzero values −1 and 1 in the positions corresponding to the third term in (14). L is given by
(16)L=[-110⋯0-11⋯⋮⋮⋮⋯⋯⋯-11].To solve (15), we introduce the idea with two auxiliary variables based on Split Bregman, which was proposed in our previous work in [27], and change (15) as the following form:
(17)min⁡β{12∥y-Dβ∥22+λ∥u∥1+η∥v∥1}s.t.u=β,v=Lβ.
Note that the minimization with the augmentation Lagrangian function of (17) can be expressed as follows:
(18)min⁡β,u,v{12∥y-Dβ∥22+λ∥u∥1+η∥v∥1+〈p,β-u〉+〈q,Lβ-v〉+μ12∥β-u∥22+μ22∥Lβ-v∥2212}.
With the alternating direction method [28] (ADM), we can solve (18) with the following iterative algorithm:
(19)(βk+1,uk+1,vk+1)=argmin{12∥y-Dβ∥22+λ∥u∥1+η∥v∥1+〈pk,β-u〉+〈qk,Lβ-v〉+μ12∥β-u∥22+μ22∥Lβ-v∥2212∥y-Dβ∥22},(20)pk+1=pk+μ1(βk+1-uk+1),(21)qk+1=qk+μ2(Lβk+1-vk+1).Note that the initial idea to get the solution of (14) is presented in [21], which is called coordinate descent method (CDM). But, in practice, CDM is slow, complex, and difficult to implement to some specialized algorithms. In addition, it is not effective in some large-scale processing problems.Furthermore, by the iterative Split Bregman, we can divide (19) into the following three independent minimizations, whose convergence is guaranteed by [29]:
(22)βk+1=argminβ{12∥y-Dβ∥22+μ12∥β-uk+1μ1pk∥22+μ22∥Lβ-vk+1μ2qk∥22},(23)uk+1=argminu{μ12∥u-(βk+1+1μ1pk)∥22+λ∥u∥1},(24)vk+1=argminv{μ22∥v-(Lβk+1+1μ2qk)∥22+η∥v∥1}.As to the quadratic minimization in (22), we can obtain the solution by solving the following linear equation:
(25)(DTD+μ1I+μ2LTL)βk+1=DTy+μ1(uk-μ1-1pk)+μ2LT(vk-μ2-1qk),
where I is the identity matrix with the same size to LTL.While thel1 minimization in (23) and (24) can be solved by the shrinkage technique. So we can obtain the uk+1 and vk+1 by
(26)uk+1=Sμ1-1λ(βk+1+μ1-1pk),(27)vk+1=Sμ2-1η(Lβk+1+μ2-1qk),
where Sλ(x)=(tλ(x(1)),tλ(x(2)),…) and tλ(ω)=sgn⁡(ω)max⁡(0,|ω|-λ).We summarized the coding algorithm in Algorithm3.Algorithm 3
Initialization:β0=DTy, u0=v0=p0=q0=0, η=λ=0.2, μ1=μ2=0.1;
Fork=0 : MaxIternum
(1) Compute the βk+1 by  (25);
(2) Compute the uk+1 by  (26);
(3) Compute the vk+1 by  (27);
(4) Update the pk+1 and qk+1 by  (20) and  (21);
End
Output: the recovered patchx^=Dβ^.Now, with all the preparation works above, we can summarize the synthetic denoising algorithm in Algorithm4.Algorithm 4
Input: A set of overlapping patches from noisy imageY-={y-i}i=1P;
Output: The recovered imageX^;
(1) Fort=1:T do
(2)  Compute the weights vectorωi for y-i by Algorithm 1;
(3)  Cluster the patches intoΩk (k=1,2,…,K) by K-Means with
l1 norm measurement between ωi and the cluster center;
(4)  Fork=1:K do
(5)  Learn the dictionaryDk for cluster Ωk by Algorithm 2;
(6)  Recover all the patches inΩk by Algorithm 3;
(7)  End
(8)  Compute theX^ by averaging the recovered patches according to
their indexes. Then, updateY- by the recovered image X^.
(9) End


## 5. Experiment Results and Analysis

We conducted various experiments on image denoising to demonstrate the performance of our proposed algorithm. We degrade the images by adding artificial zero mean Gaussian noise with different standard deviations. The test images are shown in Figure1. The sizes of images in experiment are all 256 × 256, and the patch size is 8 × 8; that is, N is 64 in Algorithm 1. Empirically, we set the parameters r=1.0 and h=2.4 for all the experiment images in Algorithm 1. We fixed the window size in (5) to be 9 × 9, where we found the small size, such as 5 × 5, may not capture the local geometric structure of underlying image data sometimes. Surely, the window size can be tuned according to concrete experiment requirement. In particular, we extend the image boundary with the “symmetric” type according to the window size. The clustering number in K-Means is flexible, which is small in image with compact structure such as House in Figure 1, and large with complex structure such as Lena. According to the image underlying structure, in our experiments, the clustering number K is optional between 5 to 10. The maximum iteration T in Algorithm 4 is 10. In addition, the regularization parameters are set in the algorithms themselves.Test images.
(a)
Lena
(b)
House
(c)
Couple
(d)
ManWe compared our proposed algorithm to several current excellent denoising approaches, including the FGTV method in [30], denoising with the dictionary learned by KSVD in [6] (DKSVD), BM3D with Wavelet (BM3DW) in [13], the kernel regression method in [14], and the two-stage denoising method with PCA (TSPCA) in [9]. Due to the limited space, we only show the experiment results of Lena and Couple with noise standard deviation σ=20 in Figures 2 and 3, respectively. Furthermore, the PSNR results of all the recovered images are reported in Figure 4.Denoising results of Lena.
(a)
TV
(b)
BM3DW
(c)
KSVD
(d)
KR
(e)
TSPCA
(f)
ProposedDenoising results of Couple.
(a)
TV
(b)
BM3DW
(c)
KSVD
(d)
KR
(e)
TSPCA
(f)
ProposedPSNR results of test images.
(a)
PSNR results of Lena
(b)
PSNR results of House
(c)
PSNR results of Couple
(d)
PSNR results of ManFrom the denoising results, we can note that the FGTV method has the worst visualization among the compared methods. Because it only constrains the total variation but does not consider the local structure adaptively, so it lost many details in original image and also showed the disadvantage in PSNR results. KR algorithm generates many mottled artifacts in denoised image, but it indeed preserves some texture by capturing the local structure with kernel function, such as the curtain in Couple. But its results declined rapidly with the increasing noise standard deviation both in visual quality and in PSNR. TSPCA generated specific smoothness in the recovered image and was weak to present the texture region. In addition, the other factor resulting in bad performance in texture is that the PCA dictionary is learned with neighboring patches, which will show weaker results than the nonlocal scheme. We can see the PSNR of BM3DW, DKSVD, and the proposed method are approximate to each other in Figure4, although the BM3DW obtained the highest PSNR compared to other methods. But BM3DW and DKSVD have more distortions in part of texture regions than the proposed method. This is because in BM3DW, the wavelet is not a good representation for all types of images, such as the ones with some complex structures. Also, although the KSVD can obtain a dictionary learned by image itself to better present the structures in image, it produces a universal dictionary, which may be not effective to certain local structures. Compared to BM3DW and KSVD, the proposed method showed the advantage to present the underlying local geometric structures in image, such as the hairs in Lena and the face of woman in Couple.


## 6. Conclusion

In this paper, we propose a novel scheme for image denoising. To preserve the textures in image, we clustered the patches from the noisy image with the meaningful weights vector which can capture the underlying local structure. And then, we learned the dictionary to better present the patches for each cluster with PCA. Lastly, we coded the noisy patches with the learned dictionary by fused lasso and obtained recovery image. We compared our proposed scheme to some current excellent algorithms, and it can be seen that our method obtains the good performance both in visual quality and in PSNR among the compared methods. In addition, the dictionary learning and coding are performed independently in each cluster, so it can be easily parallelized by the processors with multiple cores when the image has been clustered. That means the proposed method can be used in the large-scale image denoising task and save more computational time effectively.



---

*Source: 368602-2014-08-28.xml*