# Content-Aware Compressive Sensing Recovery Using Laplacian Scale Mixture Priors and Side Information

Authors: Zhonghua Xie; Lihong Ma; Lingjun Liu
Journal: Mathematical Problems in Engineering
 (2018)

---

## Abstract

Nonlocal methods have shown great potential in many image restoration tasks including compressive sensing (CS) reconstruction through use of image self-similarity prior. However, they are still limited in recovering fine-scale details and sharp features, when rich repetitive patterns cannot be guaranteed; moreover the CS measurements are corrupted. In this paper, we propose a novel CS recovery algorithm that combines nonlocal sparsity with local and global prior, which soften and complement the self-similarity assumption for irregular structures. First, a Laplacian scale mixture (LSM) prior is utilized to model dependencies among similar patches. For achieving group sparsity, each singular value of similar packed patches is modeled as a Laplacian distribution with a variable scale parameter. Second, a global prior and a compensation-based sparsity prior of local patch are designed in order to maintain differences between packed patches. The former refers to a prediction which integrates the information at the independent processing stage and is used as side information, while the latter enforces a small (i.e., sparse) prediction error and is also modeled with the LSM model so as to obtain local sparsity. Afterward, we derive an efficient algorithm based on the expectation-maximization (EM) and approximate message passing (AMP) frame for the maximum a posteriori (MAP) estimation of the sparse coefficients. Numerical experiments show that the proposed method outperforms many CS recovery algorithms.


---

## Body

## 1. Introduction

Compressive sensing (CS) [1, 2] allows us to reconstruct high-dimensional data with only a small number of random samples or measurements, if the original signal can be sparsely represented by some given appropriate basis. Owing to the fact that image prior knowledge plays a critical role in the performance of compressive sensing reconstruction, much efforts have been made to develop an effective regularization term or signal model to reflect the image prior knowledge. Standard CS methods exploit the sparsity of signal in some domains, such as DCT [3], wavelets [4, 5], total variation (TV) [6, 7], and learned dictionary [8, 9]. Unfortunately, these methods are less appropriate for many imaging applications. The reason for this failure is that natural images do not have an exactly sparse representation in any above basis. These models favor piecewise constant image structures and hence tend to smooth much the image details.More recently, the concept of sparsity has evolved into various sophisticated forms, including group sparsity [10, 11], tree sparsity [12–14], and nonlocal sparsity [15–19], where higher-order dependency among sparse coefficients is exploited. Among them, nonlocal sparsity, which refers to the fact that a patch often has many nonlocal similar patches to it across the image, has been shown to be most beneficial to CS image recovery. In [15], a nonlocal total variation (NLTV) regularization model for CS image recovery is proposed by using the self-similarity property in gradient domain. In order to obtain an adaptive sparsity regularization term for CS image recovery process, a local piecewise autoregressive model is designed in [16]. In [17], similar patches are grouped to form a two-dimensional data matrix for characterizing the low-rank property, leading to a CS recovery method via nonlocal low-rank regularization (NLR-CS). In [18, 19], a probabilistic graphical model is established, which uses collaborative filtering [20] to promote sparsity of packed patches. Despite the steady progress in nonlocal methods, they still tend to smooth the detailed image textures, degrading the image visual quality, for the reason that the lack of self-repetitive structures and corruption for data is unavoidable.To deal with this issue, local and global priors are designed to soften and complement the nonlocal sparsity for irregular structures so as to preserve image details. More specifically, the nonlocal sparsity is only imposed on a set of patches with limited influence from neighboring pixels, while the global prior refers to a prediction used as a reference which integrates the outcomes at the independent processing stage and can maintain the entire consistency of image; moreover, a compensation-based constraint term of local patch is utilized to enforce a small (i.e., sparse) prediction error. Both local sparsity and nonlocal sparsity are represented by Laplacian scale mixture (LSM) [21, 22] models, which are adopted to force coefficients, that is, singular values of local patches and similar packed patches, to be sparse. Each coefficient is modeled as a Laplacian distribution with a variable scale parameter, resulting in weighted singular value minimization problems, where weights are adaptively assigned according to the signal-to-noise ratio. On the other hand, the reference image can be used as side information. Finally, we obtain a side information-aided LSM prior model for CS image reconstruction. To solve this model, the expectation-maximization (EM) [23] method is adopted, turning the CS recovery problem into a prior parameter estimation problem and a singular value minimization problem. In particular, owing to its promising performance and efficiency, we are motivated to apply the approximate message passing (AMP) algorithm [24, 25], which is an iterative algorithm that can be used in signal and image reconstruction by performing denoising at each iteration, to solve the latter. Experimental results on natural images show that our approach can achieve more accurate reconstruction than other competing approaches.


## 2. Background

### 2.1. CS Recovery Problem

The CS recovery problem aims to find the sparsest solutionx∈Cn from the underdetermined linear system y=Ax+w, where y∈Cm are the measurements, A∈Cm×n, m<n is the measurement matrix, and w denotes the additive noise. One can solve the following objective function: (1)x^=arg⁡minx⁡y-Ax222+λRx.The first term is the data fidelity term that represents the closeness of the solution to the measurements. The second term is a regularization term that represents a priori sparse information of the original signal. λ is a regularization parameter that balances the contribution of both terms. As mentioned in Introduction, CS recovery methods exploit the sparsity of signal in some domains, such as DCT [3], wavelets [4, 5], learned dictionary [8, 9], and total variation (TV) [6, 7], leading to various forms of R(x):DCT(x)p, Ψ(x)p, and Dic(x)p where p∈{0,1} and xTV, respectively.


### 2.2. Nonlocal Sparsity

The abundance of self-repeating patterns in natural image (as shown in Figure1) can be characterized by the nonlocal sparsity [16, 17]. As shown in Figure 2, for each local patch, we can find the first M most similar nonlocal patches to it. In practice, this can be done by Euclidean distance based block matching in a large enough local window. Let Rix (or xi) denote an exemplar patch located at the ith position. Patches that are similar to Rix including Rix itself are found to form a low-rank matrix Cix=[Ri1x,Ri2x,…,RiNx]: Cix∈RM×N. Here, we suppose that M≥N. An objective function that reflects the group sparsity of similar patches with a low-rank regularization term for CS recovery can be formulated as follows: (2)x^=arg⁡minx⁡y-Ax222+∑i=1GCix∗,where G is the total number of similar patch groups; Cix∗ is the nuclear norm of Cix, taking a sum value of its singular values, namely, Cix∗=∑jsi,j. si=[si,1,si,2,…,si,N] denotes the singular value vector of Cix, that is, Cix=UiΣiViT, and si is a vector that contains the diagonal elements of Σi.Figure 1
The nonlocal similar patches in natural images.Figure 2
Illustration for the low-rank matrices construction.


### 2.3. The Approximate Message Passing Algorithm

On the other hand, these minimizing problems (e.g., (1) and (2)) can be solved easily by iterative shrinkage/thresholding (IST) methods [26–28], alternating direction method of multipliers (ADMM) [29, 30], or Bregman iterative algorithms [31, 32]. The approximate message passing reconstruction algorithm defined by Donoho et al. [24] has recently become a popular algorithm for solving signal reconstruction problems in linear systems as defined in (1). It is based on the theory of belief propagation in graphical models. Unlike belief propagation that needs to calculate 2mn messages in each iteration, by employing quadratic approximation, the expression of each message can be simplified in the AMP algorithm, and the number of messages can be reduced to m+n [33]. The final alternating expressions to solve the objective function minx⁡y-Ax22/2+λΨ(x)1 are(3)xt+1=ηxt+A∗zt,(4)zt=y-Axt+zt-1η′xt-1+A∗zt-11m,where Ψ(·) denotes wavelet transforms; x(t) and z(t) are the estimates of x and the residual at iteration t. The iteration starts from x(0)=0 and z(0)=y. A∗ is the conjugate transpose of A. The functions η(·) and η′(·) are the wavelet threshold function and its first derivative, respectively. The last term in (4) is called the “Onsager reaction term” [24] in statistical physics. This Onsager reaction term helps improve the phase transition (trade-off between the measurement rate and signal sparsity) of the reconstruction process over existing IST algorithms [26–28]. We can summarize the AMP algorithm in three steps: a residue update step (i.e., (4)), a back-projection step to yield a noisy image q(t)=x(t)+A∗z(t), and a proximal denoising correction (i.e., (3)). These steps are identical to the ones of IST methods. As a result, the AMP algorithm can be viewed as a special IST method. One benefit of this interpretation is that if R(x) varies, only the proximal denoising operator needs to be altered, when the AMP algorithm is utilized to solve (1). Some AMP variants [13, 34, 35] have been proposed with various forms of R(x), such as total variation [34], a Cauchy prior in the wavelet domain [35], and tree sparsity [13]. This motivates us to find a more suitable prior for natural images to improve the AMP algorithm.


## 2.1. CS Recovery Problem

The CS recovery problem aims to find the sparsest solutionx∈Cn from the underdetermined linear system y=Ax+w, where y∈Cm are the measurements, A∈Cm×n, m<n is the measurement matrix, and w denotes the additive noise. One can solve the following objective function: (1)x^=arg⁡minx⁡y-Ax222+λRx.The first term is the data fidelity term that represents the closeness of the solution to the measurements. The second term is a regularization term that represents a priori sparse information of the original signal. λ is a regularization parameter that balances the contribution of both terms. As mentioned in Introduction, CS recovery methods exploit the sparsity of signal in some domains, such as DCT [3], wavelets [4, 5], learned dictionary [8, 9], and total variation (TV) [6, 7], leading to various forms of R(x):DCT(x)p, Ψ(x)p, and Dic(x)p where p∈{0,1} and xTV, respectively.


## 2.2. Nonlocal Sparsity

The abundance of self-repeating patterns in natural image (as shown in Figure1) can be characterized by the nonlocal sparsity [16, 17]. As shown in Figure 2, for each local patch, we can find the first M most similar nonlocal patches to it. In practice, this can be done by Euclidean distance based block matching in a large enough local window. Let Rix (or xi) denote an exemplar patch located at the ith position. Patches that are similar to Rix including Rix itself are found to form a low-rank matrix Cix=[Ri1x,Ri2x,…,RiNx]: Cix∈RM×N. Here, we suppose that M≥N. An objective function that reflects the group sparsity of similar patches with a low-rank regularization term for CS recovery can be formulated as follows: (2)x^=arg⁡minx⁡y-Ax222+∑i=1GCix∗,where G is the total number of similar patch groups; Cix∗ is the nuclear norm of Cix, taking a sum value of its singular values, namely, Cix∗=∑jsi,j. si=[si,1,si,2,…,si,N] denotes the singular value vector of Cix, that is, Cix=UiΣiViT, and si is a vector that contains the diagonal elements of Σi.Figure 1
The nonlocal similar patches in natural images.Figure 2
Illustration for the low-rank matrices construction.


## 2.3. The Approximate Message Passing Algorithm

On the other hand, these minimizing problems (e.g., (1) and (2)) can be solved easily by iterative shrinkage/thresholding (IST) methods [26–28], alternating direction method of multipliers (ADMM) [29, 30], or Bregman iterative algorithms [31, 32]. The approximate message passing reconstruction algorithm defined by Donoho et al. [24] has recently become a popular algorithm for solving signal reconstruction problems in linear systems as defined in (1). It is based on the theory of belief propagation in graphical models. Unlike belief propagation that needs to calculate 2mn messages in each iteration, by employing quadratic approximation, the expression of each message can be simplified in the AMP algorithm, and the number of messages can be reduced to m+n [33]. The final alternating expressions to solve the objective function minx⁡y-Ax22/2+λΨ(x)1 are(3)xt+1=ηxt+A∗zt,(4)zt=y-Axt+zt-1η′xt-1+A∗zt-11m,where Ψ(·) denotes wavelet transforms; x(t) and z(t) are the estimates of x and the residual at iteration t. The iteration starts from x(0)=0 and z(0)=y. A∗ is the conjugate transpose of A. The functions η(·) and η′(·) are the wavelet threshold function and its first derivative, respectively. The last term in (4) is called the “Onsager reaction term” [24] in statistical physics. This Onsager reaction term helps improve the phase transition (trade-off between the measurement rate and signal sparsity) of the reconstruction process over existing IST algorithms [26–28]. We can summarize the AMP algorithm in three steps: a residue update step (i.e., (4)), a back-projection step to yield a noisy image q(t)=x(t)+A∗z(t), and a proximal denoising correction (i.e., (3)). These steps are identical to the ones of IST methods. As a result, the AMP algorithm can be viewed as a special IST method. One benefit of this interpretation is that if R(x) varies, only the proximal denoising operator needs to be altered, when the AMP algorithm is utilized to solve (1). Some AMP variants [13, 34, 35] have been proposed with various forms of R(x), such as total variation [34], a Cauchy prior in the wavelet domain [35], and tree sparsity [13]. This motivates us to find a more suitable prior for natural images to improve the AMP algorithm.


## 3. Side Information-Aided LSM Prior Modeling

Image nonlocal self-similarity has been widely adopted in patch based CS image reconstruction methods. Despite the great success, most of the existing works exploit the nonlocal sparsity from the degraded image, which may cause the mismatching issue in the block matching stage. Moreover, they treat irregular and regular structures equally, resulting in an oversmoothed outcome. In fact, unlike large-scale edges, the fine-scale textures have much higher randomness in local structure and they are hard to characterize by using a local or nonlocal model. In this paper, we soften and complement the nonlocal sparsity for irregular structures by combining with local and global priors and propose a side information-aided LSM prior model for achieving detail-preserving and content-aware CS image reconstruction.


### 3.1. The Laplacian Scale Mixture Distribution

A random variableθi=κi-1ςi is a Laplacian scale mixture if ςi has a Laplacian distribution with scale 1; that is, p(ςi)=exp(-|ςi|)/2, and the multiplier variable κi is a positive random variable with probability p(κi) [21]. Supposing that ςi and κi are independent, conditioned on the parameter κi, the coefficient θi has a Laplacian distribution with inverse scale p(θi∣κi)=κiexp⁡(-κi|θi|)/2. The distribution over θi is therefore a continuous mixture of Laplacian distributions with different inverse scales: (5)pθi=∫0∞pθi∣κipκidκi=∫0∞κi2exp⁡-κiθipκidκi.The distribution in (5) is defined as a Laplacian scale mixture. Note that, for most choices of p(κi), we do not have an analytical expression for p(θi).


### 3.2. The Proposed Model

In this section, we formulate the side information-aided LSM prior model and apply the MAP theory to estimate the original signal. As discussed earlier, the global information of image refers to a prediction that is actually a reference image. We use this reference image denoted byu as side information which is supposed to be known in each iteration to assist the reconstruction. Using Bayesian formula, we might derive the following MAP estimation problem:(6)x^=arg⁡minx⁡px∣y,u=arg⁡minx⁡py,u∣xpxpy,u=arg⁡minx⁡py∣xpu∣xpxpy,u=arg⁡minx⁡py∣xpu∣xpx.The side information u and the measurements y are supposed to be independent. According to (6), we need to define three terms p(y∣x), p(u∣x), and p(x) for the MAP estimation.First, the additive noisew is assumed to be white Gaussian with variance σ2, that is, w~N(0,σ2). Thus, we have the following likelihood of CS:(7)py∣x=12πσ2m/2exp⁡-12σ2y-Ax22.Second, we characterize the nonlocal sparsity of image with the Laplacian scale mixture model. Recall that, in Section2, we define si=[si,1,si,2,…,si,N] as the singular value vector of the low-rank matrix Cix. For each coefficient si,j, we assign it a Laplacian distribution with a variable scale parameter,(8)psi,j=∫0∞psi,j∣γi,jpγi,jdγi,j=∫0∞γi,j2exp⁡-γi,jsi,jpγi,jdγi,j,with a Gamma distribution prior over the scale parameter, that is, γi,j~Gamma(a,b). Note that the mean of the Laplacian distribution is 0. Hence, the mean subtraction should be carried out as shown in Figure 2. The observation that the singular values can be modeled by a Laplacian distribution has been proposed and validated in [36]. Here, we extend this idea by viewing scale parameters as random variables for achieving a better spatial adaptation. Assume that the sparse coefficients are i.i.d., and then the LSM prior of the coefficients can be expressed as p(s)=∏i=1G∏j=1Np(si,j).Third, for irregular structures, the packed patches are less similar. We wish to preserve the relative difference among these patches. Supposing that the reference imageu is available, we wish to get a solution that has a small Euclidean distance fromu while being constrained by the self-similarity prior:(9)pu∣x=∏i=1Hpui∣xi=∏i=1Hexp⁡-12σu2xi-ui22,where ui and xi denote the patches located at the ith position and H is the total number of patches. In practice, the reference image u is produced by minimizing the data fidelity term with the gradient descent method [26]. Note that the resulting image u has fine details as well as a small amount of noise; thus the noise level might be a little higher than the estimated one. The image u which integrates the outcomes at the independent processing stage can represent the global information of image. This approach is similar to the one reported in [37], which is a denoising method. It reduces the local-global gap by encouraging the overlapping patches to reach an agreement before they merge their forces by the averaging.At last, when the packed patches have an obvious dissimilarity, a local constraint imposed on each local patch is added by enforcing a small (i.e., sparse) prediction errorei=xi-vi, where v is the noiseless version of the image u. To get v, the mean filtering is applied. The singular value vector of the prediction error of the ith patch is denoted by ri=[ri,1,ri,2,…,ri,K], which is imposed on a LSM distribution p(ri)=∏j=1Kp(ri,j), where(10)pri,j=∫0∞pri,j∣πi,jpπi,jdπi,j=∫0∞πi,j2exp⁡-πi,jri,jpπi,jdπi,j,with a Gamma distribution prior over the scale parameter, that is, πi,j~Gamma(c,d). K is the total number of singular values of the ith patch.The proposed hierarchical model can be summarized in Figure3. We will have an objective function that can be maximized with respect to x, if we observe the latent variable γ and π. The standard approach in machine learning when confronted with such a problem is the EM algorithm. Note that once the spare coefficients (i.e., s and r) are determined, the image x can be obtained by averaging all reconstructed patches.Figure 3
The proposed hierarchical model.


## 3.1. The Laplacian Scale Mixture Distribution

A random variableθi=κi-1ςi is a Laplacian scale mixture if ςi has a Laplacian distribution with scale 1; that is, p(ςi)=exp(-|ςi|)/2, and the multiplier variable κi is a positive random variable with probability p(κi) [21]. Supposing that ςi and κi are independent, conditioned on the parameter κi, the coefficient θi has a Laplacian distribution with inverse scale p(θi∣κi)=κiexp⁡(-κi|θi|)/2. The distribution over θi is therefore a continuous mixture of Laplacian distributions with different inverse scales: (5)pθi=∫0∞pθi∣κipκidκi=∫0∞κi2exp⁡-κiθipκidκi.The distribution in (5) is defined as a Laplacian scale mixture. Note that, for most choices of p(κi), we do not have an analytical expression for p(θi).


## 3.2. The Proposed Model

In this section, we formulate the side information-aided LSM prior model and apply the MAP theory to estimate the original signal. As discussed earlier, the global information of image refers to a prediction that is actually a reference image. We use this reference image denoted byu as side information which is supposed to be known in each iteration to assist the reconstruction. Using Bayesian formula, we might derive the following MAP estimation problem:(6)x^=arg⁡minx⁡px∣y,u=arg⁡minx⁡py,u∣xpxpy,u=arg⁡minx⁡py∣xpu∣xpxpy,u=arg⁡minx⁡py∣xpu∣xpx.The side information u and the measurements y are supposed to be independent. According to (6), we need to define three terms p(y∣x), p(u∣x), and p(x) for the MAP estimation.First, the additive noisew is assumed to be white Gaussian with variance σ2, that is, w~N(0,σ2). Thus, we have the following likelihood of CS:(7)py∣x=12πσ2m/2exp⁡-12σ2y-Ax22.Second, we characterize the nonlocal sparsity of image with the Laplacian scale mixture model. Recall that, in Section2, we define si=[si,1,si,2,…,si,N] as the singular value vector of the low-rank matrix Cix. For each coefficient si,j, we assign it a Laplacian distribution with a variable scale parameter,(8)psi,j=∫0∞psi,j∣γi,jpγi,jdγi,j=∫0∞γi,j2exp⁡-γi,jsi,jpγi,jdγi,j,with a Gamma distribution prior over the scale parameter, that is, γi,j~Gamma(a,b). Note that the mean of the Laplacian distribution is 0. Hence, the mean subtraction should be carried out as shown in Figure 2. The observation that the singular values can be modeled by a Laplacian distribution has been proposed and validated in [36]. Here, we extend this idea by viewing scale parameters as random variables for achieving a better spatial adaptation. Assume that the sparse coefficients are i.i.d., and then the LSM prior of the coefficients can be expressed as p(s)=∏i=1G∏j=1Np(si,j).Third, for irregular structures, the packed patches are less similar. We wish to preserve the relative difference among these patches. Supposing that the reference imageu is available, we wish to get a solution that has a small Euclidean distance fromu while being constrained by the self-similarity prior:(9)pu∣x=∏i=1Hpui∣xi=∏i=1Hexp⁡-12σu2xi-ui22,where ui and xi denote the patches located at the ith position and H is the total number of patches. In practice, the reference image u is produced by minimizing the data fidelity term with the gradient descent method [26]. Note that the resulting image u has fine details as well as a small amount of noise; thus the noise level might be a little higher than the estimated one. The image u which integrates the outcomes at the independent processing stage can represent the global information of image. This approach is similar to the one reported in [37], which is a denoising method. It reduces the local-global gap by encouraging the overlapping patches to reach an agreement before they merge their forces by the averaging.At last, when the packed patches have an obvious dissimilarity, a local constraint imposed on each local patch is added by enforcing a small (i.e., sparse) prediction errorei=xi-vi, where v is the noiseless version of the image u. To get v, the mean filtering is applied. The singular value vector of the prediction error of the ith patch is denoted by ri=[ri,1,ri,2,…,ri,K], which is imposed on a LSM distribution p(ri)=∏j=1Kp(ri,j), where(10)pri,j=∫0∞pri,j∣πi,jpπi,jdπi,j=∫0∞πi,j2exp⁡-πi,jri,jpπi,jdπi,j,with a Gamma distribution prior over the scale parameter, that is, πi,j~Gamma(c,d). K is the total number of singular values of the ith patch.The proposed hierarchical model can be summarized in Figure3. We will have an objective function that can be maximized with respect to x, if we observe the latent variable γ and π. The standard approach in machine learning when confronted with such a problem is the EM algorithm. Note that once the spare coefficients (i.e., s and r) are determined, the image x can be obtained by averaging all reconstructed patches.Figure 3
The proposed hierarchical model.


## 4. EM-AMP Algorithm for CS Recovery

In this section, we simultaneously learn the hidden parameters and do inference. To accomplish this task, we embed the AMP algorithm within an EM framework.


### 4.1. EM Learning of the Prior Parameters

We use Jensen’s inequality and obtain the following upper bound on the posterior likelihood, denoted byJ(Θ,x), where Θ={Q,D}, according to (6):(11)-log⁡px∣y,u≤-log⁡py∣x-log⁡pu∣x-∫γQγlog⁡ps,γQγdr-∫πDπlog⁡pr,πDπdπ,where(12)log⁡ps,γ=log⁡ps∣γ+log⁡pγ=∑i=1G∑j=1N-γi,jsi,j+log⁡γi,j2+log⁡pγ,log⁡pr,π=log⁡pr∣π+log⁡pπ=∑i=1H∑j=1K-πi,jri,j+log⁡πi,j2+log⁡pπ.Performing coordinate descent in the auxiliary function J(Θ,x) leads to the following updates that are usually called the E step and the M step:(13)Estep:Θt+1=arg⁡minΘ⁡JΘ,xt,(14)Mstep:xt+1=arg⁡minx⁡JΘt+1,x.The E step and the M step represent the prior learning and the signal reconstruction, respectively. First, we have equality in the Jensen inequality if Q(γ)=p(γ∣s) and D(π)=p(π∣r), which implies that the E step can be reduced to Q(t+1)=p(γ∣s(t)) and D(t+1)=p(π∣r(t)). Second, let τ and ε denote the expectation with respect to Q(γ) and D(π), respectively; that is, τi,j=∫γi,jγi,jQ(γi,j)dγi,j and εi,j=∫πi,jπi,jD(πi,j)dπi,j. The M Step (14) simplifies to(15)xt+1=arg⁡minx,s,r⁡y-Ax222σ2+∑i=1Hxi-ui222σu2+∑i=1G∑j=1Nτi,jt+1si,j+∑i=1H∑j=1Kεi,jt+1ri,j.Finally, we have the proposed objective function for CS recovery. These four terms are the data fidelity term, side information-aided term, the nonlocal regularization, and the local regularization, respectively.The Gamma distribution and Laplacian distribution are conjugate; that is, the posterior probability ofγi,j (or πi,j) given si,j (or ri,j) is also a Gamma distribution with parameters a+1 (or c+1) and b+|si,j| (or d+|ri,j|). Hence, the expectations are given by(16)τi,jt+1=∫γi,jγi,jpγi,j∣si,jtdγi,j=a+1b+si,jt,εi,jt+1=∫πi,jπi,jpπi,j∣ri,jtdπi,j=c+1d+ri,jt.


### 4.2. Image Recovery via the AMP Algorithm

Once the prior parameters are estimated, the composite regularization problem, that is, (15), can be solved with various CS reconstruction algorithms [29–32]. Owing to its promising performance, the AMP algorithm [24, 25] is employed. As discussed earlier, the proximal operator varies according to the regularization term R(x). Hence we have the following proximal operator l(q(t)) on the basis of (15):(17)xt+1=arg⁡minx,s,r⁡x-qt222σ2+∑i=1Hxi-ui222σu2+∑i=1G∑j=1Nτi,jt+1si,j+∑i=1H∑j=1Kεi,jt+1ri,j,where q(t)=x(t)+A∗z(t) denotes a noisy image yielded by back-projection and l(q(t)) is a denoising operation. To solve this composite regularization problem, that is, (17), we decompose it into two simpler regularization subproblems by using the composite splitting algorithm (CSA) [7], which includes three steps: (1) splitting variable x into two variables x1 and x2, (2) performing operator splitting to minimize the nonlocal regularization and the local regularization subproblems over x1 and x2, respectively, and (3) obtaining the solution x by linear combination of x1 and x2, that is, x=ω1x1+ω2x2. Subproblems can be written as(18)x1t+1=arg⁡minx,s⁡x-qt222σ2+∑i=1Hxi-ui222σu2+∑i=1G∑j=1Nτi,jt+1si,j=arg⁡minx,s⁡x-qt+σ2ut/σu2222σ2+∑i=1G∑j=1Nτi,jt+1si,j,(19)x2t+1=arg⁡minx,r⁡x-qt+σ2ut/σu2222σ2+∑i=1H∑j=1Kεi,jt+1ri,j.In order to solve (18), we transform the data from spatial domain to SVD domain for the unity of variable representation by (1) dividing the images into overlapped patches; (2) building the low-rank matrices Xi=Cix and Qi(t)=Ci(q(t)+σ2u(t)/σu2), and then we have x-(q(t)+σ2u(t)/σu2)22=∑i=1GXi-Qi(t)F2/α, where α is a scale factor for balancing the increasing energy caused by the repetitive computation because of the overlap division and ·F is the Frobenius norm [17]; (3) performing the singular value decomposition on Xi and Qi(t) to get the singular value vectors si=[si,1,si,2,…,si,N] and ρi=[ρi,1,ρi,2,…,ρi,N]. Based on the von Neumann trace inequality [38], we know that tr((Qi(t))TXi) achieves its upper bound ∑j=1Nρjsj. Then we have(20)Xi-QitF2=trQitTQit-2trQitTXi+trXiTXi≥∑j=1Nρj2-2∑j=1Nρjsj+∑j=1Nsj2.By incorporating (20) into (18), it yields(21)st+1=arg⁡mins⁡∑i=1G∑j=1Nρi,j2-2∑j=1Nρi,jsi,j+∑j=1Nsi,j22ασ2+∑j=1Nτi,jt+1si,j=arg⁡mins⁡∑i=1G∑j=1Nsi,j-ρi,j+ασ2τi,jt+12.Taking a derivative with respect to si,j, we have(22)si,jt+1=max⁡ρi,j-ασ2τi,jt+1,0,which is the global optimum of si,j. The noise variance σ2 is obtained by maximum likelihood estimation [39]; that is, σ2=z(t)22/m. Afterward, we obtain the matrix constructed by similar patches, that is, Xi(t+1)=Uidiag⁡(si(t+1))ViT, and then recover x1(t+1) by averaging all reconstructed patches.Similarly, we can derive the global optimum of another subproblem by incorporating (20) into (19):(23)ri,jt+1=max⁡ϕi,j-σ2πi,jt+1,0,where ϕi is the singular value vector of the ith patch. Then x2(t+1) can be recovered by aggregating all reconstructed patches. Because r is the singular value of the prediction error, the reconstructed image solved by this subproblem is the prediction error in fact. Thus, we add v to obtain the final image x2(t+1).To apply the AMP algorithm, computing the Onsager termz(t)l′(q(t))1/m, which involves computing the derivative l′(q(t)), is required. It is not easy to get l′(q(t)), since l(q(t)) do not have an explicit input-output relation. Thanks to the Monte Carlo (MC) method [40], we can simulate l′(q(t)) with random numbers. This method has been used to estimate the derivative of BM3D denoiser in [18]. More details can be found in [18].


### 4.3. Content-Aware Strategy

In order to preserve fine details, we distinguish irregular structures from regular structures through the similarity of packed patches. The local regularization and side information-aided term are only utilized for irregular structures. They are adopted to maintain differences among packed patches. We use the normalized mean squared error (NMSE) as the similarity measure:(24)E=∑jx~i-xj22x~i22,where x~i is an exemplar patch located at the ith position and xj denotes the jth patch in the similar patch group. The objective function (17) is then modified to(25)xt+1=arg⁡minx,s,r⁡x-qt222σ2+∑i=1HWi⊙xi-ui222σu2+∑i=1G∑j=1Nτi,jt+1si,j+∑i=1H∑j=1KW-i⊙εi,jt+1ri,j,where Wi and W-i denote the weights that correspond to the ith patch and ⊙ is the Hadamard product. Wi and W-i are determined by(26)Wi=0,ifE<TH11,otherwise,W-i=1,ifE>TH20,otherwise,where TH1 and TH2 (TH2>TH1) are two threshold parameters. First, if the patches are similar, that is, E<TH1, only the nonlocal regularization is utilized. Second, if the packed patches are not so similar, that is, TH1<E<TH2, the side information-aided term is added. At last, when the patches have an obvious dissimilarity, that is, E>TH2, the local regularization is superadded.In general, the proposed algorithm is summarized in Algorithm1, named as the AMP algorithm with side information-aided LSM priors (SI-LSM-AMP). We also propose an algorithm as a by-product, named as the AMP algorithm with nonlocal LSM prior (NL-LSM-AMP), which is a special form of SI-LSM-AMP by setting E<TH1 for all image patches. By comparing SI-LSM-AMP with NL-LSM-AMP in the experiments, one can validate the effectiveness of our content-aware strategy.Algorithm 1:CS image recovery via SI-LSM-AMP.
Input: y,A,T,a,b,c,d,ω1,ω2, TH1, TH2, x0=0, z0=y.
for  t=0  to  T-1  do
(a) Approximate the Onsager correction term via MC.
(b) Update the residual z(t)=y-Ax(t)+z(t-1)l′(q(t-1))1/m.
(c) Obtain the noisy image q(t)=x(t)+A∗z(t).
(d) Calculate the proximal operator l(q(t)) (i.e., solve (25))
for  i=1  to  G  do
(I) Construct the low-rank matrix Qi(t).
(II) Distinguish irregular structures from regular structures with the similarity measure E, and set
W
i and W-i via (26).
(III) Perform the SVD on Qi(t) to get the singular value vector [ρi,1,ρi,2,…,ρi,N].
(IV) Estimate the expectations of scale parameters τi,j(t+1) and εi,j(t+1) via (16), and the noise variance
σ
2
=
z
(
t
)
2
2
/
m.
(V) Compute the global optimums of coefficients si,j(t+1) and ri,j(t+1) via (22) and (23).
(VI) If i=G, recover the whole image x(t+1)=ω1x1(t+1)+ω2x2(t+1) by aggregating all recovered pixels.
end for
end for


## 4.1. EM Learning of the Prior Parameters

We use Jensen’s inequality and obtain the following upper bound on the posterior likelihood, denoted byJ(Θ,x), where Θ={Q,D}, according to (6):(11)-log⁡px∣y,u≤-log⁡py∣x-log⁡pu∣x-∫γQγlog⁡ps,γQγdr-∫πDπlog⁡pr,πDπdπ,where(12)log⁡ps,γ=log⁡ps∣γ+log⁡pγ=∑i=1G∑j=1N-γi,jsi,j+log⁡γi,j2+log⁡pγ,log⁡pr,π=log⁡pr∣π+log⁡pπ=∑i=1H∑j=1K-πi,jri,j+log⁡πi,j2+log⁡pπ.Performing coordinate descent in the auxiliary function J(Θ,x) leads to the following updates that are usually called the E step and the M step:(13)Estep:Θt+1=arg⁡minΘ⁡JΘ,xt,(14)Mstep:xt+1=arg⁡minx⁡JΘt+1,x.The E step and the M step represent the prior learning and the signal reconstruction, respectively. First, we have equality in the Jensen inequality if Q(γ)=p(γ∣s) and D(π)=p(π∣r), which implies that the E step can be reduced to Q(t+1)=p(γ∣s(t)) and D(t+1)=p(π∣r(t)). Second, let τ and ε denote the expectation with respect to Q(γ) and D(π), respectively; that is, τi,j=∫γi,jγi,jQ(γi,j)dγi,j and εi,j=∫πi,jπi,jD(πi,j)dπi,j. The M Step (14) simplifies to(15)xt+1=arg⁡minx,s,r⁡y-Ax222σ2+∑i=1Hxi-ui222σu2+∑i=1G∑j=1Nτi,jt+1si,j+∑i=1H∑j=1Kεi,jt+1ri,j.Finally, we have the proposed objective function for CS recovery. These four terms are the data fidelity term, side information-aided term, the nonlocal regularization, and the local regularization, respectively.The Gamma distribution and Laplacian distribution are conjugate; that is, the posterior probability ofγi,j (or πi,j) given si,j (or ri,j) is also a Gamma distribution with parameters a+1 (or c+1) and b+|si,j| (or d+|ri,j|). Hence, the expectations are given by(16)τi,jt+1=∫γi,jγi,jpγi,j∣si,jtdγi,j=a+1b+si,jt,εi,jt+1=∫πi,jπi,jpπi,j∣ri,jtdπi,j=c+1d+ri,jt.


## 4.2. Image Recovery via the AMP Algorithm

Once the prior parameters are estimated, the composite regularization problem, that is, (15), can be solved with various CS reconstruction algorithms [29–32]. Owing to its promising performance, the AMP algorithm [24, 25] is employed. As discussed earlier, the proximal operator varies according to the regularization term R(x). Hence we have the following proximal operator l(q(t)) on the basis of (15):(17)xt+1=arg⁡minx,s,r⁡x-qt222σ2+∑i=1Hxi-ui222σu2+∑i=1G∑j=1Nτi,jt+1si,j+∑i=1H∑j=1Kεi,jt+1ri,j,where q(t)=x(t)+A∗z(t) denotes a noisy image yielded by back-projection and l(q(t)) is a denoising operation. To solve this composite regularization problem, that is, (17), we decompose it into two simpler regularization subproblems by using the composite splitting algorithm (CSA) [7], which includes three steps: (1) splitting variable x into two variables x1 and x2, (2) performing operator splitting to minimize the nonlocal regularization and the local regularization subproblems over x1 and x2, respectively, and (3) obtaining the solution x by linear combination of x1 and x2, that is, x=ω1x1+ω2x2. Subproblems can be written as(18)x1t+1=arg⁡minx,s⁡x-qt222σ2+∑i=1Hxi-ui222σu2+∑i=1G∑j=1Nτi,jt+1si,j=arg⁡minx,s⁡x-qt+σ2ut/σu2222σ2+∑i=1G∑j=1Nτi,jt+1si,j,(19)x2t+1=arg⁡minx,r⁡x-qt+σ2ut/σu2222σ2+∑i=1H∑j=1Kεi,jt+1ri,j.In order to solve (18), we transform the data from spatial domain to SVD domain for the unity of variable representation by (1) dividing the images into overlapped patches; (2) building the low-rank matrices Xi=Cix and Qi(t)=Ci(q(t)+σ2u(t)/σu2), and then we have x-(q(t)+σ2u(t)/σu2)22=∑i=1GXi-Qi(t)F2/α, where α is a scale factor for balancing the increasing energy caused by the repetitive computation because of the overlap division and ·F is the Frobenius norm [17]; (3) performing the singular value decomposition on Xi and Qi(t) to get the singular value vectors si=[si,1,si,2,…,si,N] and ρi=[ρi,1,ρi,2,…,ρi,N]. Based on the von Neumann trace inequality [38], we know that tr((Qi(t))TXi) achieves its upper bound ∑j=1Nρjsj. Then we have(20)Xi-QitF2=trQitTQit-2trQitTXi+trXiTXi≥∑j=1Nρj2-2∑j=1Nρjsj+∑j=1Nsj2.By incorporating (20) into (18), it yields(21)st+1=arg⁡mins⁡∑i=1G∑j=1Nρi,j2-2∑j=1Nρi,jsi,j+∑j=1Nsi,j22ασ2+∑j=1Nτi,jt+1si,j=arg⁡mins⁡∑i=1G∑j=1Nsi,j-ρi,j+ασ2τi,jt+12.Taking a derivative with respect to si,j, we have(22)si,jt+1=max⁡ρi,j-ασ2τi,jt+1,0,which is the global optimum of si,j. The noise variance σ2 is obtained by maximum likelihood estimation [39]; that is, σ2=z(t)22/m. Afterward, we obtain the matrix constructed by similar patches, that is, Xi(t+1)=Uidiag⁡(si(t+1))ViT, and then recover x1(t+1) by averaging all reconstructed patches.Similarly, we can derive the global optimum of another subproblem by incorporating (20) into (19):(23)ri,jt+1=max⁡ϕi,j-σ2πi,jt+1,0,where ϕi is the singular value vector of the ith patch. Then x2(t+1) can be recovered by aggregating all reconstructed patches. Because r is the singular value of the prediction error, the reconstructed image solved by this subproblem is the prediction error in fact. Thus, we add v to obtain the final image x2(t+1).To apply the AMP algorithm, computing the Onsager termz(t)l′(q(t))1/m, which involves computing the derivative l′(q(t)), is required. It is not easy to get l′(q(t)), since l(q(t)) do not have an explicit input-output relation. Thanks to the Monte Carlo (MC) method [40], we can simulate l′(q(t)) with random numbers. This method has been used to estimate the derivative of BM3D denoiser in [18]. More details can be found in [18].


## 4.3. Content-Aware Strategy

In order to preserve fine details, we distinguish irregular structures from regular structures through the similarity of packed patches. The local regularization and side information-aided term are only utilized for irregular structures. They are adopted to maintain differences among packed patches. We use the normalized mean squared error (NMSE) as the similarity measure:(24)E=∑jx~i-xj22x~i22,where x~i is an exemplar patch located at the ith position and xj denotes the jth patch in the similar patch group. The objective function (17) is then modified to(25)xt+1=arg⁡minx,s,r⁡x-qt222σ2+∑i=1HWi⊙xi-ui222σu2+∑i=1G∑j=1Nτi,jt+1si,j+∑i=1H∑j=1KW-i⊙εi,jt+1ri,j,where Wi and W-i denote the weights that correspond to the ith patch and ⊙ is the Hadamard product. Wi and W-i are determined by(26)Wi=0,ifE<TH11,otherwise,W-i=1,ifE>TH20,otherwise,where TH1 and TH2 (TH2>TH1) are two threshold parameters. First, if the patches are similar, that is, E<TH1, only the nonlocal regularization is utilized. Second, if the packed patches are not so similar, that is, TH1<E<TH2, the side information-aided term is added. At last, when the patches have an obvious dissimilarity, that is, E>TH2, the local regularization is superadded.In general, the proposed algorithm is summarized in Algorithm1, named as the AMP algorithm with side information-aided LSM priors (SI-LSM-AMP). We also propose an algorithm as a by-product, named as the AMP algorithm with nonlocal LSM prior (NL-LSM-AMP), which is a special form of SI-LSM-AMP by setting E<TH1 for all image patches. By comparing SI-LSM-AMP with NL-LSM-AMP in the experiments, one can validate the effectiveness of our content-aware strategy.Algorithm 1:CS image recovery via SI-LSM-AMP.
Input: y,A,T,a,b,c,d,ω1,ω2, TH1, TH2, x0=0, z0=y.
for  t=0  to  T-1  do
(a) Approximate the Onsager correction term via MC.
(b) Update the residual z(t)=y-Ax(t)+z(t-1)l′(q(t-1))1/m.
(c) Obtain the noisy image q(t)=x(t)+A∗z(t).
(d) Calculate the proximal operator l(q(t)) (i.e., solve (25))
for  i=1  to  G  do
(I) Construct the low-rank matrix Qi(t).
(II) Distinguish irregular structures from regular structures with the similarity measure E, and set
W
i and W-i via (26).
(III) Perform the SVD on Qi(t) to get the singular value vector [ρi,1,ρi,2,…,ρi,N].
(IV) Estimate the expectations of scale parameters τi,j(t+1) and εi,j(t+1) via (16), and the noise variance
σ
2
=
z
(
t
)
2
2
/
m.
(V) Compute the global optimums of coefficients si,j(t+1) and ri,j(t+1) via (22) and (23).
(VI) If i=G, recover the whole image x(t+1)=ω1x1(t+1)+ω2x2(t+1) by aggregating all recovered pixels.
end for
end for


## 5. Experiments

In this section, we report the experimental results of the proposed CS recovery method SI-LSM-AMP.


### 5.1. Experiment Setup

We generate the CS measurements by randomly sampling the Fourier transform coefficients of test images; that is,A is partial Fourier transform with m rows and n columns. Thus, the sampling ratio is m/n. We follow the sampling strategy of previous works ([7, 12]), which randomly choose more Fourier coefficients from low frequency and less on high frequency and set the sampling ratio near to 0.2, as CS imaging is always interested in low sampling ratio cases. All measurements are mixed with Gaussian white noise with standard deviations 5 and 15, representing the environments with low noise and high noise, respectively. Peak Signal-to-Noise Ratio (PSNR) is used for quantitative evaluation.We conduct experiments on 22 natural images, as shown in Figure4. Besides the NL-LSM-AMP algorithm, we compare the SI-LSM-AMP algorithm with image reconstruction algorithms, including the original AMP method [24], two tree-based algorithms: WaTMRI [12] and Turbo-AMP [13], and two nonlocal sparsity algorithms: NLR-CS [17] and BM3D-AMP [19]. For fair comparisons, all codes are downloaded from the authors’ websites. All algorithms terminate after 50 iterations, except 20 iterations for Turbo-AMP and 260 iterations for NLR-CS. For the WaTMRI algorithm, in order to achieve the best result, the regularization parameters are tuned to 0.8 and 0.35. The original AMP uses Daubechies wavelet to decompose the images. Then the wavelet denoising is applied with a threshold 0.8σ. For the rest of algorithms, default settings in their codes are adopted. All experiments are on a desktop with 3.80 GHz AMD A10-5800K CPU. Matlab version is R2014a.Figure 4
The 22 test images.The main parameters of the proposed algorithms are set as follows: patch size is6×6; a total of 36 similar patches are selected for each exemplar patch. To reduce the computational complexity, we extract exemplar image patch in every 5 pixels along both horizontal and vertical directions. For better CS recovery performance, some parameters are tuned empirically, including (1)  σu2=3.3σ2, which implies that the proportion between the data fidelity term and the side information-aided term is 1 : 0.3; (2)  α=5.2 in (22); (3) the parameters of Gamma distribution a=c=0 and b=d=0.01 as suggested in [21]; (4) the combination parameters ω1=2/3 and ω2=1/3. Recall that the noise level should be a little higher (1.08 times in our experiments) than the estimated one because the reference image introduces a small amount of noise. The experimental results including average PSNR, visual quality, and runtime are present.


### 5.2. Average PSNR Evaluation

To reduce the randomness, we run each experiment 5 times for each image. The average PSNR results are shown in Figures5 and 6. Figure 5 shows the PSNR results in the environment with low noise, while Figure 6 shows the ones with high noise. From these figures, one can draw the following four conclusions. First, four nonlocal sparsity-based methods SI-LSM-AMP, NL-LSM-AMP, BM3D-AMP, and NLR-CS perform better than others, which implies that the assumption of the nonlocal sparsity structure is more appropriate for natural images when compared to the tree structure and the standard sparsity. Second, the highest PSNR results are achieved by the proposed algorithm SI-LSM-AMP. In fact, the PSNR gain of the SI-LSM-AMP algorithm over the next-best algorithm NL-LSM-AMP can be as much as 0.25 dB on average. This result validates the effectiveness of our content-aware strategy that treats irregular and regular structures differently. Third, through the use of the LSM prior model, the SI-LSM-AMP algorithm and the NL-LSM-AMP algorithm outperform the other nonlocal sparsity-based AMP method BM3D-AMP (in fact, by 0.88 dB and 0.63 dB on average). Therefore, we can conclude that the LSM prior model is more appropriate for representing the nonlocal sparsity of natural images. At last, the PSNR curves decline with the decrease in sampling ratio or with the increase in measurement noise. The SI-LSM-AMP algorithm always performs better with low sampling ratio or strong measurement noise. For better quality perception, experiments performed on Baboon and Boat images present the similar results in Figures 5(b) and 5(c) and Figures 6(b) and 6(c), respectively. In a word, these results validate the superiority of the proposed SI-LSM-AMP algorithm in objective quality and the effectiveness of the LSM prior model.Figure 5
Average PSNR at different sampling ratios with measurement noise with standard deviation 5. (a) Comparisons on 22 natural images; (b) comparisons on Baboon image; (c) comparisons on Boat image.
(a)
(b)
(c)Figure 6
Average PSNR at different sampling ratios with measurement noise with standard deviation 15. (a) Comparisons on 22 natural images; (b) comparisons on Baboon image; (c) comparisons on Boat image.
(a)
(b)
(c)


### 5.3. Visual Quality and Runtime Evaluation in Lower-Power Noise Environment

Figures7 and 8 show the visual comparisons of the reconstructed results on Baboon and Boat images with 20% sampling by different methods with measurement noise with standard deviation 5, while the corresponding iterative curves are given in Figures 9 and 10, respectively. From Figures 7 and 8, we can clearly see that four nonlocal sparsity-based methods are still better than others. Among them, the SI-LSM-AMP algorithm enjoys great advantages over the NL-LSM-AMP algorithm in producing clearer image, for example, on the area of hair (Figure 7) and beach (Figure 8). It can not only perfectly reconstruct large-scale sharp edges but also well recover small-scale fine structures. The images reconstructed by the NL-LSM-AMP algorithm and the BM3D-AMP algorithm are too smooth. These two methods all have a strong assumption of the nonlocal self-similarity structure. However, many images with irregular structures do not strictly follow this assumption. We could find that the SI-LSM-AMP algorithm has great superiority on the images with irregular structures, because we combine local and global sparsity, which soften and complement the nonlocal self-similarity structure assumption for irregular structures. We also compute the PSNR as well as the structural similarity index (SSIM), which better reflects the visual quality of the images. The SI-LSM-AMP algorithm achieves the highest quantitative results. The superiority of the proposed SI-LSM-AMP in visual quality could be demonstrated by these results.Figure 7
Visual comparisons on Baboon image at 20% sampling ratio with measurement noise with standard deviation 5. (a) The original image; (b) AMP (PSNR = 25.44 dB, SSIM = 0.6968); (c) Turbo-AMP (PSNR = 23.82 dB, SSIM = 0.6103); (d) WaTMRI (PSNR = 26.37 dB, SSIM = 0.7590); (e) NLR-CS (PSNR = 26.56 dB, SSIM = 0.7877); (f) BM3D-AMP (PSNR = 26.56 dB, SSIM = 0.7478); (g) NL-LSM-AMP (PSNR = 26.97 dB, SSIM = 0.7792); (h) SI-LSM-AMP (PSNR = 27.41 dB, SSIM = 0.8051).
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)Figure 8
Visual comparisons on Boat image at 20% sampling ratio with measurement noise with standard deviation 5. (a) The original image; (b) AMP (PSNR = 26.00 dB, SSIM = 0.7702); (c) Turbo-AMP (PSNR = 23.89 dB, SSIM = 0.6748); (d) WaTMRI (PSNR = 27.39 dB, SSIM = 0.8192); (e) NLR-CS (PSNR = 28.81 dB, SSIM = 0.8687); (f) BM3D-AMP (PSNR = 29.09 dB, SSIM = 0.8755); (g) NL-LSM-AMP (PSNR = 29.42 dB, SSIM = 0.8861); (h) SI-LSM-AMP (PSNR = 29.67 dB, SSIM = 0.8927).
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)Figure 9
Iterative curves on Baboon image at 20% sampling ratio with measurement noise with standard deviation 5. (a) Average PSNR to iterations; (b) average PSNR to CPU time.
(a)
(b)Figure 10
Iterative curves on Boat image at 20% sampling ratio with measurement noise with standard deviation 5. (a) Average PSNR to iterations; (b) average PSNR to CPU time.
(a)
(b)The CPU time and PSNR are traced in each iteration for each of the methods. Figures9 and 10 present the CPU time versus PSNR curves and iteration number versus PSNR curves, respectively. We can see that the SI-LSM-AMP algorithm achieves higher PSNR results after about 30 iterations in Figures 9(a) and 10(a) and after about 100 seconds in Figures 9(b) and 10(b). Note that, considering the fact that it is hard to distinguish irregular structures from regular structures at the beginning stage, we only use the nonlocal sparsity constraint to obtain an initial reconstructed image in the first 30 iterations of the SI-LSM-AMP algorithm. Therefore, the PSNR results of the SI-LSM-AMP algorithm in the first 30 iterations are identical to the ones of the NL-LSM-AMP algorithm. The SI-LSM-AMP algorithm is relatively slow. The main computational burdens are introduced by iteratively applying SVD on each of the patch groups. The computational time of our algorithm can be further reduced through the use of the parallel computing method or C language development. Besides, in the case of NLR-CS with more iterations, its results are every five iterations presented in Figures 9 and 10. These curves demonstrate that the SI-LSM-AMP algorithm can converge to a good reconstructed result in a reasonable amount of time.


### 5.4. Visual Quality and Runtime Evaluation in Higher-Power Noise Environment

In this subsection, we conduct similar experiments with noisy CS measurements subject to strong measurement noise to demonstrate the robustness of the proposed SI-LSM-AMP to noise. The standard derivation of additive Gaussian noise is 15. The subjective quality comparison results on Baboon image and Boat image are shown in Figures11 and 12, respectively. One can observe that (1) the quality of all reconstructed images degrades seriously as measurement noise increases; (2) compared to the reconstructed images in the presence of lower-power measurement noise, the ones in higher-power noise environment tend to be more noisy (e.g., AMP, Turbo-AMP, WaTMRI, and NLR-CS) or smooth (e.g., BM3D-AMP and NL-LSM-AMP); however, the SI-LSM-AMP algorithm can better remove the artifacts and preserve important image structures more effectively even when the measurement noise is high. The corresponding PSNR and SSIM results of these algorithms are also provided, from which we can see that the proposed algorithm achieves the highest quantitative results. This indicates that the SI-LSM-AMP algorithm is shown to be more robust to noise.Figure 11
Visual comparisons on Baboon image at 20% sampling ratio with measurement noise with standard deviation 15. (a) The original image; (b) AMP (PSNR = 24.16 dB, SSIM = 0.6047); (c) Turbo-AMP (PSNR = 23.27 dB, SSIM = 0.5405); (d) WaTMRI (PSNR = 23.99 dB, SSIM = 0.6342); (e) NLR-CS (PSNR = 25.33 dB, SSIM = 0.6972); (f) BM3D-AMP (PSNR = 25.32 dB, SSIM = 0.6755); (g) NL-LSM-AMP (PSNR = 25.55 dB, SSIM = 0.6962); (h) SI-LSM-AMP (PSNR = 25.82 dB, SSIM = 0.7159).
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)Figure 12
Visual comparisons on Boat image at 20% sampling ratio with measurement noise with standard deviation 15. (a) The original image; (b) AMP (PSNR = 24.34 dB, SSIM = 0.6725); (c) Turbo-AMP (PSNR = 23.24 dB, SSIM = 0.6235); (d) WaTMRI (PSNR = 24.25 dB, SSIM = 0.6793); (e) NLR-CS (PSNR = 26.00 dB, SSIM = 0.7398); (f) BM3D-AMP (PSNR = 26.17 dB, SSIM = 0.7682); (g) NL-LSM-AMP (PSNR = 26.58 dB, SSIM = 0.7895); (h) SI-LSM-AMP (PSNR = 26.84 dB, SSIM = 0.8021).
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)The corresponding CPU time versus PSNR curves and iteration number versus PSNR curves are given in Figures13 and 14. From these figures, we can see that all algorithms converge to reconstructed results more quickly in higher-power noise environment. Among them, the PSNR results of the proposed SI-LSM-AMP algorithm are higher than all other competing methods.Figure 13
Iterative curves on Baboon image at 20% sampling ratio with measurement noise with standard deviation 15. (a) Average PSNR to iterations; (b) average PSNR to CPU time.
(a)
(b)Figure 14
Iterative curves on Boat image at 20% sampling ratio with measurement noise with standard deviation 15. (a) Average PSNR to iterations; (b) average PSNR to CPU time.
(a)
(b)We could find that, compared to other nonlocal methods, the SI-LSM-AMP algorithm has great superiority in higher-power noise environment. These results are reasonable because the mismatching issue is worse when the measurement noise is higher, leading to the difficulty in finding similar patches. However, other nonlocal algorithms depend on self-similarity structure only, which makes them hard to perform well. In contrast, besides the self-similarity structure, the SI-LSM-AMP algorithm also uses the local and global sparse prior of natural images and thus has a soft assumption of the self-similarity structure. As a result, one advantage of the SI-LSM-AMP algorithm is that it is more appropriate for compressive imaging under severe environment, where compressive samples are subject to strong measurement noise.


## 5.1. Experiment Setup

We generate the CS measurements by randomly sampling the Fourier transform coefficients of test images; that is,A is partial Fourier transform with m rows and n columns. Thus, the sampling ratio is m/n. We follow the sampling strategy of previous works ([7, 12]), which randomly choose more Fourier coefficients from low frequency and less on high frequency and set the sampling ratio near to 0.2, as CS imaging is always interested in low sampling ratio cases. All measurements are mixed with Gaussian white noise with standard deviations 5 and 15, representing the environments with low noise and high noise, respectively. Peak Signal-to-Noise Ratio (PSNR) is used for quantitative evaluation.We conduct experiments on 22 natural images, as shown in Figure4. Besides the NL-LSM-AMP algorithm, we compare the SI-LSM-AMP algorithm with image reconstruction algorithms, including the original AMP method [24], two tree-based algorithms: WaTMRI [12] and Turbo-AMP [13], and two nonlocal sparsity algorithms: NLR-CS [17] and BM3D-AMP [19]. For fair comparisons, all codes are downloaded from the authors’ websites. All algorithms terminate after 50 iterations, except 20 iterations for Turbo-AMP and 260 iterations for NLR-CS. For the WaTMRI algorithm, in order to achieve the best result, the regularization parameters are tuned to 0.8 and 0.35. The original AMP uses Daubechies wavelet to decompose the images. Then the wavelet denoising is applied with a threshold 0.8σ. For the rest of algorithms, default settings in their codes are adopted. All experiments are on a desktop with 3.80 GHz AMD A10-5800K CPU. Matlab version is R2014a.Figure 4
The 22 test images.The main parameters of the proposed algorithms are set as follows: patch size is6×6; a total of 36 similar patches are selected for each exemplar patch. To reduce the computational complexity, we extract exemplar image patch in every 5 pixels along both horizontal and vertical directions. For better CS recovery performance, some parameters are tuned empirically, including (1)  σu2=3.3σ2, which implies that the proportion between the data fidelity term and the side information-aided term is 1 : 0.3; (2)  α=5.2 in (22); (3) the parameters of Gamma distribution a=c=0 and b=d=0.01 as suggested in [21]; (4) the combination parameters ω1=2/3 and ω2=1/3. Recall that the noise level should be a little higher (1.08 times in our experiments) than the estimated one because the reference image introduces a small amount of noise. The experimental results including average PSNR, visual quality, and runtime are present.


## 5.2. Average PSNR Evaluation

To reduce the randomness, we run each experiment 5 times for each image. The average PSNR results are shown in Figures5 and 6. Figure 5 shows the PSNR results in the environment with low noise, while Figure 6 shows the ones with high noise. From these figures, one can draw the following four conclusions. First, four nonlocal sparsity-based methods SI-LSM-AMP, NL-LSM-AMP, BM3D-AMP, and NLR-CS perform better than others, which implies that the assumption of the nonlocal sparsity structure is more appropriate for natural images when compared to the tree structure and the standard sparsity. Second, the highest PSNR results are achieved by the proposed algorithm SI-LSM-AMP. In fact, the PSNR gain of the SI-LSM-AMP algorithm over the next-best algorithm NL-LSM-AMP can be as much as 0.25 dB on average. This result validates the effectiveness of our content-aware strategy that treats irregular and regular structures differently. Third, through the use of the LSM prior model, the SI-LSM-AMP algorithm and the NL-LSM-AMP algorithm outperform the other nonlocal sparsity-based AMP method BM3D-AMP (in fact, by 0.88 dB and 0.63 dB on average). Therefore, we can conclude that the LSM prior model is more appropriate for representing the nonlocal sparsity of natural images. At last, the PSNR curves decline with the decrease in sampling ratio or with the increase in measurement noise. The SI-LSM-AMP algorithm always performs better with low sampling ratio or strong measurement noise. For better quality perception, experiments performed on Baboon and Boat images present the similar results in Figures 5(b) and 5(c) and Figures 6(b) and 6(c), respectively. In a word, these results validate the superiority of the proposed SI-LSM-AMP algorithm in objective quality and the effectiveness of the LSM prior model.Figure 5
Average PSNR at different sampling ratios with measurement noise with standard deviation 5. (a) Comparisons on 22 natural images; (b) comparisons on Baboon image; (c) comparisons on Boat image.
(a)
(b)
(c)Figure 6
Average PSNR at different sampling ratios with measurement noise with standard deviation 15. (a) Comparisons on 22 natural images; (b) comparisons on Baboon image; (c) comparisons on Boat image.
(a)
(b)
(c)


## 5.3. Visual Quality and Runtime Evaluation in Lower-Power Noise Environment

Figures7 and 8 show the visual comparisons of the reconstructed results on Baboon and Boat images with 20% sampling by different methods with measurement noise with standard deviation 5, while the corresponding iterative curves are given in Figures 9 and 10, respectively. From Figures 7 and 8, we can clearly see that four nonlocal sparsity-based methods are still better than others. Among them, the SI-LSM-AMP algorithm enjoys great advantages over the NL-LSM-AMP algorithm in producing clearer image, for example, on the area of hair (Figure 7) and beach (Figure 8). It can not only perfectly reconstruct large-scale sharp edges but also well recover small-scale fine structures. The images reconstructed by the NL-LSM-AMP algorithm and the BM3D-AMP algorithm are too smooth. These two methods all have a strong assumption of the nonlocal self-similarity structure. However, many images with irregular structures do not strictly follow this assumption. We could find that the SI-LSM-AMP algorithm has great superiority on the images with irregular structures, because we combine local and global sparsity, which soften and complement the nonlocal self-similarity structure assumption for irregular structures. We also compute the PSNR as well as the structural similarity index (SSIM), which better reflects the visual quality of the images. The SI-LSM-AMP algorithm achieves the highest quantitative results. The superiority of the proposed SI-LSM-AMP in visual quality could be demonstrated by these results.Figure 7
Visual comparisons on Baboon image at 20% sampling ratio with measurement noise with standard deviation 5. (a) The original image; (b) AMP (PSNR = 25.44 dB, SSIM = 0.6968); (c) Turbo-AMP (PSNR = 23.82 dB, SSIM = 0.6103); (d) WaTMRI (PSNR = 26.37 dB, SSIM = 0.7590); (e) NLR-CS (PSNR = 26.56 dB, SSIM = 0.7877); (f) BM3D-AMP (PSNR = 26.56 dB, SSIM = 0.7478); (g) NL-LSM-AMP (PSNR = 26.97 dB, SSIM = 0.7792); (h) SI-LSM-AMP (PSNR = 27.41 dB, SSIM = 0.8051).
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)Figure 8
Visual comparisons on Boat image at 20% sampling ratio with measurement noise with standard deviation 5. (a) The original image; (b) AMP (PSNR = 26.00 dB, SSIM = 0.7702); (c) Turbo-AMP (PSNR = 23.89 dB, SSIM = 0.6748); (d) WaTMRI (PSNR = 27.39 dB, SSIM = 0.8192); (e) NLR-CS (PSNR = 28.81 dB, SSIM = 0.8687); (f) BM3D-AMP (PSNR = 29.09 dB, SSIM = 0.8755); (g) NL-LSM-AMP (PSNR = 29.42 dB, SSIM = 0.8861); (h) SI-LSM-AMP (PSNR = 29.67 dB, SSIM = 0.8927).
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)Figure 9
Iterative curves on Baboon image at 20% sampling ratio with measurement noise with standard deviation 5. (a) Average PSNR to iterations; (b) average PSNR to CPU time.
(a)
(b)Figure 10
Iterative curves on Boat image at 20% sampling ratio with measurement noise with standard deviation 5. (a) Average PSNR to iterations; (b) average PSNR to CPU time.
(a)
(b)The CPU time and PSNR are traced in each iteration for each of the methods. Figures9 and 10 present the CPU time versus PSNR curves and iteration number versus PSNR curves, respectively. We can see that the SI-LSM-AMP algorithm achieves higher PSNR results after about 30 iterations in Figures 9(a) and 10(a) and after about 100 seconds in Figures 9(b) and 10(b). Note that, considering the fact that it is hard to distinguish irregular structures from regular structures at the beginning stage, we only use the nonlocal sparsity constraint to obtain an initial reconstructed image in the first 30 iterations of the SI-LSM-AMP algorithm. Therefore, the PSNR results of the SI-LSM-AMP algorithm in the first 30 iterations are identical to the ones of the NL-LSM-AMP algorithm. The SI-LSM-AMP algorithm is relatively slow. The main computational burdens are introduced by iteratively applying SVD on each of the patch groups. The computational time of our algorithm can be further reduced through the use of the parallel computing method or C language development. Besides, in the case of NLR-CS with more iterations, its results are every five iterations presented in Figures 9 and 10. These curves demonstrate that the SI-LSM-AMP algorithm can converge to a good reconstructed result in a reasonable amount of time.


## 5.4. Visual Quality and Runtime Evaluation in Higher-Power Noise Environment

In this subsection, we conduct similar experiments with noisy CS measurements subject to strong measurement noise to demonstrate the robustness of the proposed SI-LSM-AMP to noise. The standard derivation of additive Gaussian noise is 15. The subjective quality comparison results on Baboon image and Boat image are shown in Figures11 and 12, respectively. One can observe that (1) the quality of all reconstructed images degrades seriously as measurement noise increases; (2) compared to the reconstructed images in the presence of lower-power measurement noise, the ones in higher-power noise environment tend to be more noisy (e.g., AMP, Turbo-AMP, WaTMRI, and NLR-CS) or smooth (e.g., BM3D-AMP and NL-LSM-AMP); however, the SI-LSM-AMP algorithm can better remove the artifacts and preserve important image structures more effectively even when the measurement noise is high. The corresponding PSNR and SSIM results of these algorithms are also provided, from which we can see that the proposed algorithm achieves the highest quantitative results. This indicates that the SI-LSM-AMP algorithm is shown to be more robust to noise.Figure 11
Visual comparisons on Baboon image at 20% sampling ratio with measurement noise with standard deviation 15. (a) The original image; (b) AMP (PSNR = 24.16 dB, SSIM = 0.6047); (c) Turbo-AMP (PSNR = 23.27 dB, SSIM = 0.5405); (d) WaTMRI (PSNR = 23.99 dB, SSIM = 0.6342); (e) NLR-CS (PSNR = 25.33 dB, SSIM = 0.6972); (f) BM3D-AMP (PSNR = 25.32 dB, SSIM = 0.6755); (g) NL-LSM-AMP (PSNR = 25.55 dB, SSIM = 0.6962); (h) SI-LSM-AMP (PSNR = 25.82 dB, SSIM = 0.7159).
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)Figure 12
Visual comparisons on Boat image at 20% sampling ratio with measurement noise with standard deviation 15. (a) The original image; (b) AMP (PSNR = 24.34 dB, SSIM = 0.6725); (c) Turbo-AMP (PSNR = 23.24 dB, SSIM = 0.6235); (d) WaTMRI (PSNR = 24.25 dB, SSIM = 0.6793); (e) NLR-CS (PSNR = 26.00 dB, SSIM = 0.7398); (f) BM3D-AMP (PSNR = 26.17 dB, SSIM = 0.7682); (g) NL-LSM-AMP (PSNR = 26.58 dB, SSIM = 0.7895); (h) SI-LSM-AMP (PSNR = 26.84 dB, SSIM = 0.8021).
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)The corresponding CPU time versus PSNR curves and iteration number versus PSNR curves are given in Figures13 and 14. From these figures, we can see that all algorithms converge to reconstructed results more quickly in higher-power noise environment. Among them, the PSNR results of the proposed SI-LSM-AMP algorithm are higher than all other competing methods.Figure 13
Iterative curves on Baboon image at 20% sampling ratio with measurement noise with standard deviation 15. (a) Average PSNR to iterations; (b) average PSNR to CPU time.
(a)
(b)Figure 14
Iterative curves on Boat image at 20% sampling ratio with measurement noise with standard deviation 15. (a) Average PSNR to iterations; (b) average PSNR to CPU time.
(a)
(b)We could find that, compared to other nonlocal methods, the SI-LSM-AMP algorithm has great superiority in higher-power noise environment. These results are reasonable because the mismatching issue is worse when the measurement noise is higher, leading to the difficulty in finding similar patches. However, other nonlocal algorithms depend on self-similarity structure only, which makes them hard to perform well. In contrast, besides the self-similarity structure, the SI-LSM-AMP algorithm also uses the local and global sparse prior of natural images and thus has a soft assumption of the self-similarity structure. As a result, one advantage of the SI-LSM-AMP algorithm is that it is more appropriate for compressive imaging under severe environment, where compressive samples are subject to strong measurement noise.


## 6. Conclusion

We have proposed an effective structured AMP algorithm for CS image reconstruction. Our work has the following contributions. First, guided by structure sparsity theories, we introduce the Laplacian scale mixture distribution to model the nonlocal sparsity for higher-order sparse representation of natural images, and then it is used as a prior constraint for CS problem. By taking the scale parameter as a random variable, it makes the practical representation much more feasible for achieving a better spatial adaptation. Second, in order to maintain differences between packed patches for irregular structures, we combine local sparsity and global sparsity to soften and complement the nonlocal self-similarity structure assumption. It can substantially enhance the details of image. Afterward, an effective algorithm based on the EM and AMP frame is proposed in this paper to solve this model. The derived inference procedures are efficient to estimate both the sparse coefficients and the scale parameter. After learning the prior parameters, the AMP algorithm is utilized to solve the singular value minimization problem for achieving the accurate image reconstruction. Finally, our simulations and experiments on a variety of natural images demonstrate the superiority of the proposed algorithm to the original AMP algorithm and several tree-based and nonlocal sparsity-based algorithms or solvers in CS image recovery.



---

*Source: 7171352-2018-01-29.xml*