Sage Journals: Discover world-class research

Abstract

As the most critical part of compressive sensing theory, reconstruction algorithm has an impact on the quality and speed of image reconstruction. After studying some existing convex optimization algorithms and greedy algorithms, we find that convex optimization algorithms should possess higher complexity to achieve higher reconstruction quality. Also, fixed atomic numbers used in most greedy algorithms increase the complexity of reconstruction. In this context, we propose a novel algorithm, called variable atomic number matching pursuit, which can improve the accuracy and speed of reconstruction. Simulation results show that variable atomic number matching pursuit is a fast and stable reconstruction algorithm and better than the other reconstruction algorithms under the same conditions.

Keywords

Image reconstruction compressive sensing convex optimization algorithm greedy algorithm variable atomic number matching pursuit

Introduction

Compressive sensing (CS)¹ is a new sampling framework for acquisition and processing of signals. It has been widely used in the field of image processing, such as medical imaging, remote sensing, hyperspectral image processing, and image compression. Different from the traditional image processing methods, compressive sensing samples use information directly, which can save a lot of calculation, transmission, and storage resources. The theory also points out that any sparse signal (or compressible signal) can be exactly recovered from only a small amount of measured data.²

The successful implementation of compressive sensing requires three important steps. First, make signal sparse; then, measure sparse signal; and finally, recover original signal accurately from the observation vector. One of the most important steps is to select an appropriate reconstruction algorithm to accurately recover the original signal with high quality from measurement vector. Current reconstruction algorithms mainly include two categories.³ One is convex optimization algorithm, such as basis pursuit (BP),⁴ gradient projection (GP),⁵ iterative hard threshold (IHT), and iterative reweighted least square (IRLS). IRLS algorithm is to solve certain optimization problems by an iterative method, each step of which involves solving a weighted least squares problem. IHT algorithm solves the l₀ norm problem of NP-hard and considers the number of non-zero components of sparse signal directly to find the K support very approximating the sparse signal. The other is to use the greedy algorithm to approximate the minimum l₁ norm. Greedy algorithm include matching pursuit (MP),⁶ orthogonal matching pursuit (OMP), compression sampling matching pursuit (CoSaMP), and subspace tracking (SP).⁷ According to the ideas of MP, original signal is approximated by local optimal solution step by step.

The summary and analysis of the existing reconstruction algorithms show that the convex optimization algorithms in image reconstruction require a higher complexity in order to achieve high-quality reconstruction, the greedy algorithm select fixed number of atoms in each iteration, and reconstruction accuracy is not high. In order to solve these shortcomings, we propose a new algorithm, i.e. variable atomic number matching pursuit (VANMP). The algorithm is mainly based on greedy algorithm. Compared with the existing greedy algorithms, the number of atoms in each iteration of our new algorithm is variable, which greatly increases the probability of selecting the optimal atom. And the experiment results conclude that our new algorithm is much better than the other algorithms.

The basic theory of compressive sensing reconstruction

Compressive sensing reconstruction is to reconstruct the high dimensional original signal x with the length of $N (M << N)$ from measurement vector y of length M. Vector y is obtained by sampling. The essence of this problem is to solve the minimum l₀ norm of a non-convex optimization⁸

min x ‖ x ‖ 0 s . t . y = Φ x

(1)

where

Φ

is the measurement matrix and

s . t .

denotes ‘subject to’. But miserably solving equation (1) is an NP-hard problem that requires an exhaustive enumeration of all

C_{N}^{K}

possible combinations for the locations of the non-zero entries in x. Under normal circumstances, it can be transformed into the l₁ norm

min x ‖ x ‖ 1 s . t . y = Φ x

(2)

By this way, the non-convex optimization problem is transformed into a convex optimization problem. Candès⁹ also proved that when the random measurement matrix satisfies the restricted isometry property (RIP), equation (2) is equivalent to equation (1).

According to the compressive sensing mathematical model $θ = Ψ T x$ ( $Ψ$ is sparse matrix), we can conclude that reconstruction problem is the process to find the estimation of sparse vector $θ$

\overset{\land}{θ} = \arg min θ ‖ θ ‖ 1, s . t . y = A θ

(3)

where sensing matrix

A = Φ Ψ

A \in R M \times N

. Then we can get an estimation of x:

\overset{\land}{x} = Ψ \overset{\land}{θ}

according to

x = Ψ θ

. Algorithms for solving sparse vector

\overset{\land}{θ}

mainly consist of two major categories: convex optimization algorithm and greedy algorithm.

Overview of compressive sensing reconstruction algorithm

Iteratively reweighted least squares minimization (IRLS)

The specific algorithm steps are as follows:

Input: Sensing matrix $A = Φ Ψ$ , measurement vector y, signal sparsity K.

Output: The estimation of signal sparsity representation coefficient $\overset{\land}{θ}$ .

Step 1: Initialization: let $θ 0 = 0$ ( $θ$ is the least square solution of $y = A θ$ ), $ɛ = 1$ , iterations $t = 1$ ;

Step 2: The least square solution of each iteration can be obtained by Lagrange multiplication: $θ t = Q t A T (AQ t A T) - 1 y$ . In which, $Q t$ is the diagonal matrix with components $1 / W t$ , weight $W t = ((θ t - 1) 2 + ɛ) p / 2 - 1$ , then $Q t = diag (1 . / W t, 0)$ ;

Step 3: If $‖ θ t - θ t ‖ 2 ≺ (ɛ / 100) 1 / 2$ , then each iteration of ɛ down to the original 1/10, i.e. $ɛ = ɛ / 10$ ;

Step 4: If $ɛ ≻ 10 - 9$ and iterations $t < K$ , let $t = t + 1$ , then go to step 2, otherwise halt iteration and go to next step;

Step 5: The result $\overset{\land}{θ}$ of reconstruction is the result $θ t$ of the last iteration.

Note: p denotes l_p norm in step 2. Yin et al.⁵ proved that a suitable IRLS¹⁰ algorithm is convergent for $1 < p < 3$ . The first section of this paper also points out that the problem is translated into l₁ norm when $p = 1$ . After getting $\overset{\land}{θ}$ , the reconstructed signal $\overset{\land}{x}$ can be obtained by $\overset{\land}{x} = Ψ \overset{\land}{θ}$ .

IHT

The specific algorithm steps are as follows:

Input: Sensing matrix A = Φ Ψ, measurement vector y, signal sparsity K.

Output: Estimation of signal sparsity representation coefficient $\overset{\land}{θ}$ .

Step 1: Initialization: make $θ 0 = 0$ , residual $r 0 = y$ , iterations $t = 1$ ;

Step 2: Compute $a t = θ t - 1 + A T (y - A θ t - 1)$ ;

Step 3: Select the largest K components: $θ t = H K (a t)$ ;

Step 4: If the iteration condition is satisfied, then halt iteration, otherwise let $t = t + 1$ and go to step 2.

Note: $H K (a)$ is a nonlinear hard threshold operator, which sets all components to zero except the first K components of the largest amplitude in vector a.¹¹ After getting $\overset{\land}{θ}$ , reconstructed signal $\overset{\land}{x}$ can be obtained by $\overset{\land}{x} = Ψ \overset{\land}{θ}$ .

OMP

The specific algorithm steps are as follows¹²:

Input: Sensing matrix $A = Φ Ψ$ , measurement vector y, signal sparsity K;

Output: The estimation of signal sparsity representation coefficient $\overset{\land}{θ}$ , residual $r k = y - A k \overset{\land}{θ}$ .

Step 1: Initialization: residual $r 0 = y$ , index set $Λ 0 = Ø$ , sensing matrix $A 0 = Ø$ , iterations $t = 1$ ;

Step 2: Find index $λ t = \arg max j = 1, 2, \dots, N | á r t - 1, a j ñ |$ , $λ t$ is the corresponding position of inner product’s maximum value of column a_j ( $j = 1, 2, \dots, N$ ) of sensing matrix $A$ and residual r, then obtain $a λ t$ by $λ t$ ;

Step 3: Update index set $Λ t = Λ t - 1 \cup {λ t}$ , record the reconstructed atom set found in sensing matrix: $A t = A t - 1 \cup {a λ t}$ ;

Step 4: Compute the least square solution of $y = A t θ t$ : $\overset{\land}{θ} t = \arg min θ t ‖ y - A t θ t ‖ = (A_{t}^{T} A t) - 1 A_{t}^{T} y$ ;

Step 5: Update residual $r t = y - A t \overset{\land}{θ} t$ ;

Step 6: Let $t = t + 1$ , if $t \leq K$ , then go to step 2, otherwise halt iteration and go to next step;

Step 7: Reconstruction of the resulting $\overset{\land}{θ}$ in $Λ t$ has non-zero terms, its value is the result of last iteration.

Note: After getting $\overset{\land}{θ}$ , reconstructed signal $\overset{\land}{x}$ can be obtained by $\overset{\land}{x} = Ψ \overset{\land}{θ}$ .

CoSaMP

The specific algorithm steps are as follows:

Input: Sensing matrix $A = Φ Ψ$ , measurement vector y, signal sparsity K.

Output: The estimation of signal sparsity representation coefficient $\overset{\land}{θ}$ , residual $r k = y - A k \overset{\land}{θ}$ .

Step 1: Initialization: residual $r 0 = y$ , index set $Λ 0 = Ø$ , sensing matrix $A 0 = Ø$ , iterations $t = 1$ ;

Step 2: Compute $u = {u j | u j = | á r t - 1, a j ñ |}$ , ( $j = 1, 2, \dots, N$ ). Select $2 K$ maximum value in u, record the column numbers of $A$ corresponding to these values into set S_k;

Step 3: Update index set $Λ t = Λ t - 1 \cup {S k}$ , record the reconstructed atom set found in sensing matrix: $A t = A t - 1 \cup {a λ t}$ , (for all $a λ t \in S k$ );

Step 4: Compute the least square solution of $y = A t θ t$ : $\overset{\land}{θ} t = \arg min θ t ‖ y - A t θ t ‖ = (A_{t}^{T} A t) - 1 A_{t}^{T} y$ ;

Step 5: Select maximum K items of absolute value of $\overset{\land}{θ} t$ and record them as $\overset{\land}{θ} t k$ , which is corresponding to $A t k$ in $A t$ , record the index set as $Λ t k$ , then update the index set $Λ t = Λ t k$ ;

Step 6: Update residual $r t = y - A t k \overset{\land}{θ} t k$ ;

Step 7: Let $t = t + 1$ , if $t \leq K$ , then go to step 2, if $t > K$ or $r = 0$ , then halt iteration and go to next step;

Step 8: Reconstruction of the resulting $\overset{\land}{θ}$ in $Λ t$ has non-zero terms, its value is the result of last iteration.

Note: a_j( $j = 1, 2, \dots, N$ ) is the column vector of sensing matrix $A$ .¹³ After getting $\overset{\land}{θ}$ , reconstructed signal $\overset{\land}{x}$ can be obtained by $\overset{\land}{x} = Ψ \overset{\land}{θ}$ .

VANMP

Principle of VANMP

Each iteration of the OMP algorithm selects the corresponding position of maximum absolute value of the inner product of residual r and the column a_j( $j = 1, 2, \dots, N$ ) of sensing matrix $A$ fixedly, and stores them in the candidate set. So the algorithm chooses an atom which can best match the structure of the target signal in each iteration. The residual of OMP algorithm is always orthogonal to the selected atoms, which means that an atom cannot be selected two times, the result will be converged in a few steps. Thus, OMP has a quick computing speed. However, OMP selects only one optimal atom in each iteration, the degree of reconstruction accuracy is not very high. Compared with OMP, CoSaMP algorithm selects the $2 K$ corresponding position of maximum value and stores them into candidate set, then selects the K components of the absolute value as the optimal approximation of target signal. In this context, the reconstruction accuracy of CoSaMP algorithm is higher than that of OMP algorithm. But each iteration of the algorithm needs to select $2 K$ atoms; the reconstruction speed is much slower than that of OMP.¹⁴

Through the analysis and comparison of OMP and CoSaMP algorithms, we find that these algorithms have the disadvantages of selecting atoms too much or too less, and the number of selected atoms is fixed in each iteration. We assume that there exists such an algorithm, where each iteration selects the number of atoms which is variable, and it is better to select from an atom to multiple atoms in each iteration. By this way, we can overcome the shortcomings of OMP and CoSaMP. Based on this assumption, we propose a new reconstruction algorithm named VANMP. At the time of the first iteration, VANMP computes $u = {u j | u j = | á r t - 1, a j ñ |}$ , ( $j = 1, 2, \dots, N$ ), then selects the maximum value of u, and stores the column index number of $A$ corresponding to the value in the candidate set. In the next iteration, VANMP selects the number of maximum value in u is increased by 1, and the number of the maximum value in each iteration is not more than the sparsity K. Atomic number selected by VANMP algorithm in each iteration is variable, so it overcomes the problem of too many atoms or too few atoms resulted from selecting atomic number fixedly, so as to ensure the reconstruction speed and accuracy.

Regularization of candidate set

In this paper, we design a VANMP algorithm to select atoms by the method of variable atomic number and store them in candidate set S_k, then regularization S_k for secondary screening, once again raise the support set of reliability and signal reconstruction accuracy and speed. The purpose of regularization is to choose a subset with maximum energy as the new candidate set $S k 0$ , whose column vector and residual’s inner product absolute maximum value cannot be lesser than minimum value more than twice.

The main steps of regularization are as follows¹⁵:

Input: Candidate set S_k, sparsity K.

Output: The candidate set $S k 0$ after regularization.

Step 1: Sort the elements of candidate set S_k from large to small, and find the number of non-zero elements M in the set;

Step 2: Compare sparsity K and number of non-zero elements M, let the number of elements N is equal to that of the small one.

Step 3: Start from the elements $u (i)$ ( $i = 1 : N$ ), find all the elements $u (j)$ that satisfy the condition $| u (i) | \leq | u (j) |$ ( $j = i + 1 : N$ ), together with the elements $u (i)$ to form a subset S_i.

Step 4: Compute energy $‖ u ‖ 2 2$ ( $u \in S i$ ) of S_i, and find the subset $S k 0$ with maximum energy from S_i.

Concrete process of VANMP algorithm

The specific algorithm steps are as follows:

Input: Sensing matrix $A = Φ Ψ$ , measurement vector y, signal sparsity K.

Output: Estimation of signal sparsity representation coefficient $\overset{\land}{θ}$ , residual $r k = y - A k \overset{\land}{θ}$ .

Step 1: Initialization: residual $r 0 = y$ , index set $Λ 0 = Ø$ , sensing matrix $A 0 = Ø$ , atomic number $I = 1$ , iterations $t = 1$ ;

Step 2: Compute $u = {u j | u j = | á r t - 1, a j ñ |}$ , ( $j = 1, 2, \dots, N$ ). Select I maximum value in u, record the column numbers of $A$ corresponding to these values into set S_k, and then get $S k 0$ by regularizing S_k;

Step 3: Update index set $Λ t = Λ t - 1 \cup {S k 0}$ , record the reconstructed atom set found in sensing matrix $A t = A t - 1 \cup {a j}$ , (for all $a j \in S k 0$ );

Step 4: Compute the least square solution of $y = A t θ t$ : $\overset{\land}{θ} t = \arg min θ t ‖ y - A t θ t ‖ = (A_{t}^{T} A t) - 1 A_{t}^{T} y$ ;

Step 5: Select maximum I items of absolute value of $\overset{\land}{θ} t$ and record them as $\overset{\land}{θ} t I$ , which is corresponding to $A t I$ in $A t$ , then record the index set $Λ t I$ ;

Step 6: Update residual $r t = y - A t I \overset{\land}{θ} t I$ , if $r t = 0$ , then halt iteration and go to step 8;

Step 7: If $‖ r t ‖ 2 \geq ‖ r t - 1 ‖ 2$ , then update atomic number $I = I + 1$ and iterations $t = t + 1$ , go to step 2. If $‖ r t ‖ 2 < ‖ r t - 1 ‖ 2$ or $I > K$ , then halt iteration and go to next step;

Step 8: Reconstruction of the resulting $\overset{\land}{θ}$ in $Λ t I$ has non-zero terms, its value is the result of last iteration.

Note: After getting $\overset{\land}{θ}$ , reconstructed signal $\overset{\land}{x}$ can be obtained by $\overset{\land}{x} = Ψ \overset{\land}{θ}$ .

VANMP vs. existing algorithms

The most remarkable feature of VANMP is that it selects atomic number in each iteration is variable, which is distinguished from other iterative greedy algorithm. It can also be said that VANMP provides a generalized framework for the OMP and CoSaMP. Note that when the number of atoms selected by VANMP in each iteration is fixed to 1, VANMP becomes OMP. Similarly, when the number of atoms is fixed to $2 K$ , VANMP becomes CoSaMP. In this case, VANMP is always more accurate than OMP, and more faster than CoSaMP.

IRLS is considered to be a convex optimization algorithm with high reconstruction accuracy than most of greedy algorithms.¹⁶ Its complexity is $O (pN 3)$ (p denotes l_p norm, N denotes signal length). Our improved greedy algorithm retains the orthogonal atomic selection criteria of OMP algorithm, which guarantees the optimality of iterative structure. At the same time, we introduce regularization principle of ROMP into VANMP algorithm to carry out two screening of support set, which can realize fast and efficient selection of atoms. The complexity of VANMP is $O (K \log 2 N)$ (K denotes sparsity of signal), which is much lower than that of IRLS algorithm.

Simulation results

The simulation experiment of this paper is done in the environment of Windows 7 operating system, CPU is E4500 2.2 GHz Intel, memory is 2 GB, and software is R2012a matlab.

We know that most of the natural images themselves are non-sparse, but by transformation, they are sparse or compressible in a certain sparse matrix. Therefore, the selection of sparse matrix and measurement matrix has a great influence on the quality of the whole compressive sensing image reconstruction. In order to reflect the reconstruction effect of algorithms in the same conditions, we use $256 \times 256$ Lena gray image as test image, wavelet transform matrix as sparse matrix, and Gauss random matrix as measurement matrix. In order to test the reconstruction effect of the above five algorithms, we have run 100 times for each algorithm in different sampling rates of Lena image, and compared their peak signal-to-noise ratio (PSNR), reconstruction time, and reconstruction error. The results are shown in Figures 1 to 3.

Figure 1.

PSNR at different sampling rates.

From Figures 1 and 2, we can find that the reconstruction time of IHT algorithm is the shortest but the performance is not ideal, the reconstruction performance of IRLS algorithm is the best but it takes too long time. PSNR obtained by CoSaMP algorithm is slightly higher than the OMP algorithm, but with the increase of sampling rate, the reconstruction time of CoSaMP is growing faster. Time consumed by VANMP algorithm is almost the same as OMP algorithm, but the reconstruction performance is better than that of the OMP algorithm. And with the increase of the sampling rate, VANMP’s reconstruction time is almost not increased; with the decrease of the sampling rate, the decrease of PSNR is smaller than that of other algorithms. VANMP algorithm to reconstruct the image with PSNR is slightly lesser than IRLS, but time is far lesser than the IRLS algorithm. In addition, the reconstruction effect of VANMP algorithm without regularization is also shown in Figure 1. From Figure 3, we can find that the reconstruction error of VANMP algorithm is the smallest of the five algorithms. In order to measure the dispersion degree of PSNR of each algorithm, we plot the standard deviation of them in Figure 4. Through the comparison of these five algorithms, we can conclude that VANMP algorithm is a better compressive sensing reconstruction algorithm with high quality and fast speed.

Figure 2.

Time at different sampling rates.

Figure 3.

Reconstruction error at different sampling rates.

Figure 4.

Standard deviation of PSNR at different sampling rates.

When sampling rate is 0.74, the reconstructed images of each reconstruction algorithm are shown in Figure 5, and the reconstruction parameters are shown in Table 1.

Table 1.

The reconstruction parameters of different algorithms when sampling rate is 0.74.

Algorithms	IRLS	IHT	OMP	CoSaMP	VANMP
Time(s)	381.0625	2	15.0156	93.2031	15.8281
PSNR(dB)	38.6590	23.0817	34.6444	35.5875	38.6590
Error	5.9206e-04	0.2007	7.1901e-04	6.7015e-04	6.0045e-04

Figure 5.

The reconstructed image of different algorithms when sampling rate is 0.74.

As can be seen from Table 1, the construction time of VANMP algorithm is very close to that of OMP in the sampling rate of 0.74, both spend much less time. And its PSNR is equivalent to IRLS, significantly higher than the other algorithms, although their reconstruction effects are very good at this time. As can be seen from Figure 5, the details of the image reconstructed by VANMP algorithm are also relatively good. Through the comparison of two convex optimization algorithms and the three greedy algorithms, we can find that VANMP is a kind of algorithm with good robustness, low complexity, and strong practicability.

Conclusions

Based on the in-depth study of the compressive sensing theory and some classical reconstruction algorithms, we propose a novel matching pursuit algorithm with variable number of atoms; for the greedy algorithms, a fixed number of atoms at each iteration should be selected which leads to the shortcomings of choosing atoms multiple or less. And the simulation results show that our new algorithm is superior to other reconstruction algorithms in image reconstruction. Since we use simple sparse matrix and measurement matrix, leading to the need for a higher sampling rate to enhanced reconstruction precision, our focus will turn to search for a better sparse matrix and measurement matrix, and better to use the compressive sensing in practice.

Footnotes

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work is supported by Doctoral Fund of Henan Polytechnic University (B2012-104).

References

Richard

. Compressive sensing. IEEE Signal Process Mag 2007; 24: 118–121.

Candè

Michael

. An introduction to compressive sampling. IEEE Signal Process Mag 2008; 25: 21–30.

Rivenson

Adrian

Bahram

. Overview of compressive sensing techniques applied in holography [Invited]. Appl Opt 2013; 52: A423–A432.

Jenila

Sivanantha Raja

. Reduction of data acquisition time in Raman spectroscopy imaging using structure based compressive sampling algorithm. Opt Quant Electron 2015; 47: 3855–3862.

Yin

X-Y

Xiao-Jing

Zhao-Yun

. Pre-stack basis pursuit seismic inversion for brittleness of shale. Petrol Sci 2015; 12: 618–627.

Miandji

Joel

Jonas

. Compressive image reconstruction in reduced union of subspaces. Comput Graph Forum 2015; 34: 33–44.

François-Xavier D and Sandrine A. Generalized subspace pursuit and an application to sparse Poisson denoising. In: IEEE international conference on image processing (ICIP), Paris, France, 2014.

Zhang

. Image compressive sensing recovery using adaptively learned sparsifying basis via L0 minimization. Signal Process 2014; 103: 114–126.

Candès

. The restricted isometry property and its implications for compressed sensing. Comptes Rendus Math 2008; 346: 589–592.

10.

Daubechies

. Iteratively reweighted least squares minimization for sparse recovery. Commun Pure Appl Math 2010; 63: 1–38.

11.

Karahanoglu

Hakan

. Compressed sensing signal recovery via forward-backward pursuit. Digital Signal Process 2013; 23: 1539–1548.

12.

Chen

Defang

Jiashu

. Iterative gradient projection algorithm for two-dimensional compressive sensing sparse image reconstruction. Signal Process 2014; 104: 15–26.

13.

Davenport

Deanna

Michael

. Signal space CoSaMP for sparse recovery with redundant dictionaries. IEEE Transac Inform Theor 2013; 59: 6820–6829.

14.

Blumensath

Richard

. Non-convexly constrained image reconstruction from nonlinear tomographic X-ray measurements. Philos Transac R Soc Lond A 2015; 373: 20140393.

15.

Sajjad

Irfan

Sung

. Sparse coded image super-resolution using K-SVD trained dictionary based on regularized orthogonal matching pursuit. Bio Med Mater Eng 2015; 26: 1399–1407.

16.

Mousavi

Vishal

Trac

. Iterative convex refinement for sparse recovery. IEEE Signal Process Lett 2015; 22: 1903–1907.

Image reconstruction algorithm based on variable atomic number matching pursuit

Abstract

Keywords

Introduction

The basic theory of compressive sensing reconstruction

Overview of compressive sensing reconstruction algorithm

Iteratively reweighted least squares minimization (IRLS)

IHT

OMP

CoSaMP

VANMP

Principle of VANMP

Regularization of candidate set

Concrete process of VANMP algorithm

VANMP vs. existing algorithms

Simulation results

Conclusions

Footnotes

Declaration of conflicting interests

Funding

References