Compressive Imaging via Approximate Message Passing with Image Denoising

I Introduction

Compressed sensing (CS) has sparked a tremendous amount of research activity in recent years, because it performs signal acquisition and processing using far fewer samples than required by the Nyquist rate. Breakthroughs in CS have the potential to greatly reduce the sampling rates in numerous signal processing applications such as cameras , medical scanners, fast analog to digital converters , and high speed radar .

Compressed sensing has been used in compressive imaging, where the input signal is an image, and the goal is to acquire the image using as few measurements as possible. Acquiring images in a compressive manner requires less sampling time than conventional imaging technologies. Applications of compressive imaging appear in medical imaging , seismic imaging , and hyperspectral imaging .

I-B Related work

Many compressive imaging algorithms have been proposed in the literature. For example, Som and Schniter modeled the structure of the wavelet coefficients by a hidden Markov tree (HMT), and applied a turbo scheme that alternates between inference on the HMT structure with standard belief propagation and inference on the compressed sensing measurement structure with the generalized approximate message passing algorithm. He and Carin proposed a hierarchical Bayesian approach with Markov chain Monte Carlo (MCMC) for natural image reconstruction. Soni and Haupt exploited a hierarchical dictionary learning method and assumed that projecting images onto the learned dictionary will yield tree-sparsity, and therefore the nonzero supports of the dictionary can be identified and estimated accurately by setting an appropriate threshold.

However, existing compressive imaging algorithms may either not achieve good reconstruction quality or not be fast enough. Therefore, in this paper, we focus on a variation of a fast and effective algorithm called approximate message passing (AMP) to improve over the prior art. AMP is an iterative signal reconstruction algorithm that performs scalar denoising within each iteration, and proper selection of the denoising function used within AMP is needed to obtain better reconstruction quality. One challenge in applying image denoisers within AMP is that it may be hard to compute the so-called “Onsager reaction term” in the AMP iteration steps. The Onsager reaction term includes the derivative of the image denoising function, and thus if an image function does not have a convenient closed form, then the Onsager reaction term may be difficult to compute.

Dictionary learning is an effective technique that has attracted a great deal of attention in image denoising. Dictionary learning based methods generally achieve lower reconstruction error than wavelet-based methods. However, the learning procedure requires a large amount of training images, and may involve manual tuning. Owing to these limitations, our main focus in this paper is to integrate relatively simple and fast image denoisers into compressive imaging reconstruction algorithms.

Dabov et al. developed an image denoising strategy that employs collaborative filtering in a sparse 3-D transform domain, and they offered an efficient implementation that achieves favorable denoising quality. Other efficient denoising schemes include wavelet-based methods. A typical wavelet-based image denoiser proceeds as follows: (i) apply a wavelet transform to the image and obtain wavelet coefficients; (ii) denoise the wavelet coefficients; and (iii) apply an inverse wavelet transform to the denoised wavelet coefficients, yielding a denoised image. Two popular examples of denoisers that can be applied to the wavelet coefficients are hard thresholding and soft thresholding . Variations on the thresholding scheme can be found in ; other wavelet-based methods were proposed by Simoncelli and Adelson , Mıhçak et al. , and Moulin and Liu .

I-C Contributions

The objective of this paper is to develop compressive imaging algorithms that are fast and reconstruct well. We apply the “amplitude-scale-invariant Bayes estimator” (ABE) and the adaptive Wiener filter as image denoisers within AMP. Our numerical results with these denoisers are promising, showing that our AMP based algorithms run at least 3.5 times faster than the prior art algorithms. Moreover, with a proper choice of image denoiser within AMP, the proposed algorithm also outperforms the prior art algorithms in reconstruction quality. Although we only test for two wavelet-based image denoisers within AMP, we believe that other image denoisers would also work within the AMP framework, and these denoisers need not be wavelet-based.

The remainder of the paper is arranged as follows. We review AMP in Section II, and describe the image denoisers that we use within AMP in Section III. Numerical results are presented in Section IV, and the paper concludes with a discussion in Section V.

II Review of Approximate Message Passing

Before reviewing AMP, let us first model how measurements are obtained in a compressive imaging system.

Scalar channels: We now define scalar channels, and will describe in Section II-B how the matrix channel is converted to scalar channels in AMP. In scalar channels, the noisy observations obey

II-B Algorithmic framework

where $v^{t}_{i}\sim\mathcal{N}(0,\sigma_{t}^{2})$. The asymptotic performance of AMP can be characterized by a state evolution (SE) formalism:

where the random variables $W\sim\mathcal{N}(0,1)$ and $X\sim f_{X}$. Formal statements about SE appear in Bayati and Montanari . Note that SE (6) tracks the noise variance for AMP iterations, but the noise variance $\sigma_{t+1}^{2}$ need not necessarily be the smallest possible, unless in each iteration the denoiser $\eta_{t}(\cdot)$ achieves the minimum mean square error (MMSE). We note in passing that we have proposed a denoiser called “MixD” that achieves the MMSE of scalar channels (2), and thus applying MixD within AMP achieves the MMSE of matrix channels (1). On the other hand, it is unrealistic to expect existing image denoisers to achieve the MMSE, because the statistical distribution of natural images has yet to be determined. That said, running AMP with good image denoisers that achieve lower mean square error (MSE) may yield lower MSE in compressive imaging problems.

Finally, SE theoretically characterizes the noise variance $\sigma_{t}^{2}$ of the scalar channel at each iteration. However, the MSE performance of image denoisers cannot be characterized theoretically. Therefore, we must estimate the effective Gaussian noise level $\sigma^{2}_{t}$ empirically in each AMP iteration. The estimated noise variance $\widehat{\sigma}^{2}_{t}$ can be calculated as :

III Image Denoising within AMP

In this section, we describe how wavelet-based image denoisers are applied within AMP, and then outline two image denoisers that were proposed by Figueiredo and Nowak and Mıhçak et al. , respectively.

In image processing, one often computes the wavelet coefficients of images, applies some signal processing technique to the wavelet coefficients, and finally applies the inverse wavelet transform to the processed coefficients to obtain processed images. We now show how image denoising can be performed within AMP in the wavelet domain. Let us denote the wavelet transform by $\mathcal{W}$ and the inverse wavelet transform by $\mathcal{W}^{-1}$. By applying the wavelet transform to a vectorized image signal ${\bf x}$ (a 2-dimensional wavelet transform is used), we obtain the wavelet coefficient vector ${\bf\theta_{x}}=\mathcal{W}{\bf x}$. Conversely, ${\bf x}=\mathcal{W}^{-1}{\bf\theta_{x}}$. Therefore, the matrix channel (1) becomes ${\bf y=A}\mathcal{W}^{-1}{\bf\theta_{x}+z}$, where ${\bf A}\mathcal{W}^{-1}$ can be regarded as a new matrix in the matrix channel (1) and ${\bf\theta_{x}}$ as the corresponding input signal.

Let us express the AMP iterations (3, 4) for settings where the matrix is ${\bf A}\mathcal{W}^{-1}$,

Because the wavelet transform matrix is orthonormal, i.e., $\mathcal{W}\mathcal{W}^{T}={\bf I}=\mathcal{W}\mathcal{W}^{-1}$, it can be shown that $({\bf A}\mathcal{W}^{-1})^{T}=\mathcal{W}{\bf A}^{T}$. Therefore, the input of the denoiser $\eta_{t}(\cdot)$ (8) becomes

where ${\bf q}^{t}$ (5) is the noisy image at iteration $t$, and $\mathcal{W}{\bf q}^{t}$ is the wavelet transform applied to the noisy image.

With the above analysis of the modified AMP (8, 9), we formulate a compressive imaging procedure as follows. Let us denote the the wavelet transform of the scalar channel (5) by

where $\theta_{\bf q}^{t}=\mathcal{W}{\bf q}^{t}$, $\theta_{\bf x}=\mathcal{W}{\bf x}$, and $\theta_{\bf v}^{t}=\mathcal{W}{\bf v}^{t}$. First, ${\bf r}^{t}$ and ${\bf x}^{t}$ are initialized to all-zero vectors. Then, at iteration $t$ the algorithm proceeds as follows,

Calculate the residual term ${\bf r}^{t}$.

Calculate the noisy image ${\bf q}^{t}={\bf A}^{T}{\bf r}^{t}+{\bf x}^{t}$, and apply the wavelet transform $\mathcal{W}$ to the noisy image ${\bf q}^{t}$ to obtain wavelet coefficients $\theta_{\bf q}^{t}$, which are the inputs of the scalar denoiser $\eta_{t}(\cdot)$ in (8).

Apply the denoiser $\eta_{t}(\cdot)$ to the wavelet coefficients $\theta_{\bf q}^{t}$, and obtain denoised coefficients $\theta^{t+1}_{\bf x}$.

Apply the inverse wavelet transform $\mathcal{W}^{-1}$ to the coefficients $\theta_{\bf x}^{t+1}$ to obtain the estimated image ${\bf x}^{t+1}$, which is used to compute the residual term in the next iteration.

III-B Image denoisers

We choose to denoise the wavelet coefficients using scalar denoisers proposed by Figueiredo and Nowak and Mıhçak et al. , respectively, because these two denoisers are simple to implement while revealing promising numerical results (see Section IV). We call the algorithm where ABE is utilized within AMP “AMP-ABE,” and the algorithm where the adaptive Wiener filter is utilized “AMP-Wiener.” In both algorithms, the variance of the noise $\sigma_{t}^{2}$ in the noisy image is obtained using (7). Because we use an orthonormal wavelet transform, the noise variance in the wavelet domain is equal to that in the image domain. Although we only show how to employ two image denoisers within AMP, they serve as a proof of concept that other image denoisers could also be applied within AMP, possibly leading to further improvements in both image reconstruction quality and runtime.

Figueiredo and Nowak’s denoiser is an amplitude-scale-invariant Bayes estimator (ABE), and it is a scalar function. More specifically, for each noisy wavelet coefficient $\theta_{{\bf q},i}^{t}$ (11), the estimate of $\theta_{{\bf x},i}$ for the next iteration is

where $\sigma_{t}^{2}$ is the noise variance of the scalar channel (5) at the $t$-th AMP iteration, and $(\cdot)_{+}$ is a function such that $(u)_{+}=u$ if $u>0$ and $(u)_{+}=0$ if $u\leq 0$. Note that because the wavelet transform matrix $\mathcal{W}$ is orthonormal, the variance of the noise ${\bf v}^{t}$ (5) is equal to the variance of $\theta_{\bf v}^{t}$ (11).

The ABE function is continuous and differentiable except for two points ($\theta_{{\bf q},i}^{t}=\pm\sqrt{3}\sigma_{t}$), and we calculate the derivative of this denoising function numerically to obtain the Onsager reaction term in (4).

III-B2 Adaptive Wiener filter

Mıhçak et al. proposed a method to estimate the variances of the wavelet coefficients, and then apply the corresponding Wiener filter to each wavelet coefficient. The variance of the noisy wavelet coefficient $\theta_{{\bf q},i}^{t}$ is estimated from its neighboring coefficients. More specifically, a set of $3\times 3$ or $5\times 5$ neighboring coefficients $\mathcal{N}_{i}$ that is centered at $\theta_{{\bf q},i}^{t}$ is considered, and the variance of $\theta_{{\bf q},i}^{t}$ is estimated by averaging the sum of $(\theta_{{\bf q},k}^{t})^{2}$ where $k\in\mathcal{N}_{i}$. This method of averaging the neighboring coefficients can be regarded as first convolving a $3\times 3$ or $5\times 5$ mask of all $1$’s with the matrix of squared wavelet coefficients $\theta_{\bf q}^{t}$, and then dividing by the normalizing constant $9$ (for a $3\times 3$ mask) or $25$ (for a $5\times 5$ mask). Other masks can be applied to produce different and possibly better denoising results. For example, we have found that the mask

obtains lower MSE than other $5\times 5$ masks we have considered. Recall the scalar channel defined in (11) where the noise variance is ${\sigma}_{t}^{2}$; we estimate the variance of a noisy wavelet coefficient $\theta_{{\bf q},i}^{t}$ by $\widehat{\sigma}_{i}^{2}$, and the variance of the true wavelet coefficient $\theta_{{\bf x},i}^{t}$ by $\widehat{\sigma}_{i}^{2}-\sigma_{t}^{2}$.We use $\max\{\widehat{\sigma}_{i}^{2}-\sigma_{t}^{2},0\}$ to restrict the variance to be non-negative. Therefore, the scaling factor in the Wiener filter is given by $\frac{\widehat{\sigma}_{i}^{2}-\sigma_{t}^{2}}{(\widehat{\sigma}_{i}^{2}-\sigma_{t}^{2})+\sigma_{t}^{2}}$, and the adaptive Wiener filter being used as the denoising function can be expressed as follows,

Finally, the derivative of this denoising function with respect to $\theta_{{\bf q},i}^{t}$ is simply the scaling factor $\frac{\widehat{\sigma}_{i}^{2}-\sigma_{t}^{2}}{\widehat{\sigma}_{i}^{2}}$ of the Wiener filter, and so the Onsager reaction term in (4) can be obtained efficiently.

In standard AMP , the denoising function $\eta_{t}(\cdot)$ is separable, meaning that ${\theta}_{{\bf x},i}^{t+1}$ only depends on its corresponding noisy wavelet coefficient ${\theta}_{{\bf q},i}^{t}$. In the adaptive Wiener filter, however, the estimated variance $\widehat{\sigma}_{i}^{2}$ of each noisy wavelet coefficient depends on the neighboring coefficients of ${\theta}_{{\bf q},i}^{t}$, and so the denoising function in (13) implicitly depends on the neighboring coefficients of ${\theta}_{{\bf q},i}^{t}$. Therefore, the adaptive Wiener filter in (13) is not a strictly separable denoising function, and AMP-Wiener encounters convergence issues. Fortunately, a technique called “damping” solves for the convergence problem of AMP-Wiener. Specifically, damping is an extra step in the AMP iteration (3); instead of updating the value of ${\bf x}^{t+1}$ by the output of the denoiser $\eta_{t}({\bf A}^{T}{\bf r}^{t}+{\bf x}^{t})$, we assign a weighted sum of $\eta_{t}({\bf A}^{T}{\bf r}^{t}+{\bf x}^{t})$ and ${\bf x}^{t}$ to ${\bf x}^{t+1}$ as follows,

for some constant $0\leq\lambda<1$. It has been shown by Rangan et al. that sufficient damping ensures the convergence of AMP where the measurement matrix ${\bf A}$ is not i.i.d. Gaussian. However, we did indeed use i.i.d. Gaussian matrices in our numerical results in Section IV, and damping solved the convergence problem of AMP-Wiener, which suggests that damping may be an effective technique when various convergence issues arise in AMP based algorithms. We note in passing that other techniques such as SwAMP and ADMM-GAMP also solve for the convergence problem in AMP.

IV Numerical Results

Having described the AMP algorithm and two image denoisers , in this section we present the numerical results of applying these two denoisers within AMP.

We compare AMP-ABE and AMP-Wiener with three prior art compressive imaging algorithms, (i) Turbo-BG proposed by Som and Schniter ; (ii) Turbo-GM, also by Som and Schniter ; and (iii) a Markov chain Monte Carlo (MCMC) method by He and Carin . Both Turbo-BG and Turbo-GM are also message passing based algorithms. However, these two algorithms require more computation than AMP-ABE and AMP-Wiener, because they include two message passing procedures; the first procedure solves for dependencies between the wavelet coefficients and the second procedure is AMP. The performance metrics that we use to compare the algorithms are runtime and normalized MSE (NMSE), $\text{NMSE}({\bf x},\widehat{\bf x})=10\log_{10}(\|{\bf x-\widehat{x}}\|_{2}^{2}/\|{\bf x}\|_{2}^{2})$, where $\widehat{\bf x}$ is the estimate of the vectorized input image ${\bf x}$. In all simulations, we use the Haar wavelet transform .

Let us begin by contrasting the three prior art compressive imaging algorithms based on the numerical results provided in . Turbo-BG and Turbo-GM have similar runtimes; the NMSE of Turbo-GM is typically 0.5 dB better (lower) than the NMSE of Turbo-BG. At the same time, the NMSE of the MCMC algorithm is comparable to those of Turbo-BG and Turbo-GM, but MCMC is 30 times slower than the Turbo approaches of Som and Schniter . Other algorithms have also been considered for compressive imaging. For example, compressive sampling matching pursuit (CoSaMP) requires only half the runtime of Turbo-GM, but its NMSE is roughly 4 dB worse than that of Turbo-GM; and model based CS is twice slower than Turbo-GM and its NMSE is also roughly 4 dB worse. Therefore, we provide numerical results for Turbo-BG, Turbo-GM, MCMC, and our two proposed AMP based approaches.

Numerical setting: We downloaded 591 images from “pixel-wise labeled image database v2” at http://research. microsoft.com/en-us/projects/objectclassrecognition, and extracted image patches using the following two methods.

Method 1: A $192\times 192$ patch is extracted from the upper left corner of each image, and then the patch is resized to $128\times 128$; this image patch extraction method was used by Som and Schniter .

Method 2: A $192\times 192$ patch is extracted from the upper left corner of each image without resizing.

Result 1: Tables I and II show the NMSE and runtime averaged over the 591 image patches that are extracted by Methods 1 and 2, respectively. Runtime is measured in seconds on a Dell OPTIPLEX 9010 running an Intel(R) $\text{Core}^{\text{TM}}$ i7-860 with 16GB RAM, and the environment is Matlab R2013a. Figures 2 and 3 complement Tables I and II, respectively, by plotting the average NMSE over 591 images from iteration 1 to iteration 30.

It can be seen from Table I that the NMSE of AMP-Wiener is the best (lowest) among all the algorithms compared. At the same time, AMP-Wiener runs approximately 3.5 times faster than the Turbo approaches of Som and Schniter , and 120 times faster than MCMC . Although AMP-ABE does not outperform the competing algorithms in terms of NMSE, it runs as fast as AMP-Wiener.

Table I presents the runtimes of AMP-ABE and AMP-Wiener for image patches extracted by Method 1 with 30 iterations. However, we can see from Figure 2 that AMP-ABE and AMP-Wiener with fewer iterations already achieve NMSEs that are close to the NMSE shown in Table I. In Figure 2, the horizontal axis represents iteration numbers, and the vertical axis represents NMSE. It is shown in Figure 2 that the NMSE drops markedly from $-10$ dB to $-21$ dB for AMP-Wiener (solid line) and from $-5$ dB to $-19$ dB for AMP-ABE (dash-dot line), respectively. Note that the average NMSE is approximately $-21$ dB for AMP-Wiener and $-19$ dB for AMP-ABE around iteration 15. Therefore, we may halve the runtimes of AMP-ABE and AMP-Wiener (to approximately 1.7 seconds) by reducing the number of AMP iterations from 30 to 15.

The simulation for the larger image patches extracted by Method 2 is slow, and thus the results for Turbo-BG and MCMC have not been obtained for Table II. We believe that Turbo-BG is only slightly worse than Turbo-GM. At the same time, we did test for MCMC on several images, and found that the NMSEs obtained by MCMC were usually 0.5 dB higher than AMP-Wiener and the runtimes of MCMC usually exceeded 1,500 seconds. Similar to Figure 2, it can be seen from Figure 3 that the runtimes of our AMP based approaches could be further reduced by reducing the number of AMP iterations without much deterioration in estimation quality.

Result 2: As a specific example, Figure 4 illustrates one of the 591 image patches and the estimated patches using AMP-Wiener at iterations 1, 3, 7, 15, and 30. We also present the estimated patches using Turbo-GM and MCMC. It can be seen from Figure 4 that the estimated images using AMP-Wiener are gradually denoised as the number of iterations is increased, and the NMSE achieved by AMP-Wiener at iteration 15 already produces better reconstruction quality than Turbo-GM and MCMC.

IV-B Performance of scalar denoisers

Having seen that AMP-Wiener consistently outperforms AMP-ABE, let us now understand why AMP-Wiener achieves lower NMSE than AMP-ABE.

We test for ABE and the adaptive Wiener filter as scalar image denoisers in scalar channels as defined in (2). In this simulation, we use the 591 image patches extracted by Method 2, and add i.i.d. Gaussian noise $\mathcal{N}(0,\sigma^{2})$ to the image patches. The pixel values are normalized to be between and $1$, and we verified from the simulations for Table II that the estimated noise variances of the scalar channels in AMP iterations are typically between $1\times 10^{-4}$ and $1$. In Figure 5, the vertical axis represents NMSE, and the horizontal axis represents different noise variances varying from $1\times 10^{-4}$ to $1$ . It is shown in Figure 5 that the adaptive Wiener filter (solid line) consistently achieves lower NMSE than ABE (dash-dot line) for all noise variances, which suggests that AMP-Wiener outperforms AMP-ABE in every AMP iteration, and thus outperforms AMP-ABE when we stop iterating in iteration 30. Therefore, in order to achieve favorable reconstruction quality, it is important to select a good image denoiser within AMP. With this in mind, we include the NMSE of the image denoiser “block-matching and 3-D filtering” BM3D in Figure 5, and find that BM3D (dashed line) has lower NMSE than the adaptive Wiener filter, especially when the noise variance is large. Note that the NMSEs of ABE for difference noise variances are within 1 dB of the NMSEs of the adaptive Wiener filter, but this performance gap in scalar denoisers is amplified to more than 2 dB (refer to Table II) when applying the scalar denoisers within AMP. In other words, it is possible that applying BM3D within AMP could achieve better reconstruction quality than AMP-Wiener. However, one challenge of applying BM3D within AMP will be that it is not clear whether the Onsager reaction term in (4) can be computed in closed form or numerically, and thus an alternative way of approximating the Onsager reaction term may need to be developed. During the review process of our paper, Metzler et al. showed how to compute the Onsager correction term numerically, thus allowing to use different image denoisers within AMP.

IV-C Reconstruction quality versus measurement rate

Finally, we also evaluate the performance of each algorithm by plotting the NMSE (average NMSE over 591 images) versus the measurement rate $R=M/N$. The measurement matrix ${\bf A}$ is generated the same way as the numerical setting in Section IV-A.

Result: Figures 6 and 7 illustrate how the NMSEs achieved by AMP-Wiener and Turbo-GM vary when the measurement rate $R$ changes, where the horizontal axis represents the measurement rate $R=M/N$, and the vertical axis represents NMSE. Figures 6 shows the results for image patches extracted by Method 1, and the measurement rate $R$ varies from $0.1$ to $1$. Figure 7 shows the results for image patches extracted by Method 2. Because the simulation for $192\times 192$ image patches is relatively slow, we only show results for $R$ that varies from $0.1$ to $0.6$. It can be seen from Figures 6 and 7 that AMP-Wiener (solid line with pentagram markers) achieves lower NMSE than that of Turbo-GM (dash-dot line with asterisks) for all values of $R$.

V Discussion

In this paper, we proposed compressive imaging algorithms that apply image denoisers within AMP. Specifically, we used the “amplitude-scale-invariant Bayes estimator” (ABE) and an adaptive Wiener filter within AMP. Numerical results showed that AMP-Wiener achieves the lowest reconstruction error among all competing algorithms in all simulation settings, while AMP-ABE also offers competitive performance. Moreover, the runtimes of AMP-ABE and AMP-Wiener are significantly lower than those of MCMC and the Turbo approaches , and Figures 2 and 3 suggest that the runtimes of our AMP based algorithms could be reduced further if we accept a slight deterioration in NMSE.

Recall that the input of the denoising function $\eta_{t}$ in (3) is a noisy image with i.i.d. Gaussian noise, and so we believe that any image denoiser that deals with i.i.d. Gaussian noise can be applied within AMP. At the same time, in order to develop fast AMP based compressive imaging algorithms, the image denoisers that are applied within AMP should be fast. By comparing the denoising quality of ABE and the adaptive Wiener filter as image denoisers in scalar channels, we have seen that AMP with a better denoiser produces better reconstruction quality for compressive imaging problems. With this in mind, employing more advanced image denoisers within AMP may produce promising results for compressive imaging problems. The development of such compressive imaging algorithms is left for future work.

Acknowledgments

We thank Phil Schniter for providing information about the numerical settings evaluated in his work with Som ; Liyi Dai, Nikhil Krishnan, and Junan Zhu for useful discussions; and the reviewers for their careful evaluation of the manuscript.

References