Multimed Tools Appl DOI 10.1007/s11042-015-2975-9

Reconstruction algorithm for block-based compressed sensing based on mixed variational inequality Kaixiong Su 1 & Jian Chen 1 Lichao Su 2

1

& Weixing Wang &

Received: 16 January 2015 / Revised: 13 August 2015 / Accepted: 29 September 2015 # Springer Science+Business Media New York 2015

Abstract Block compressed sensing based on mixed variational inequality (BCS-MVI) is proposed to improve the performance of current reconstruction algorithms for block-based compressed sensing. In the measurement phase, an image is sampled block by block. In the recovery period, BCS-MVI takes the sparse regularization of the natural image as prior knowledge and approaches the target function within the entire image through the modified augmented Lagrange method (ALM) and alternating direction method (ADM) of multipliers. Moreover, for the reconstruction problem including two regularization terms, an adaptive weight (−AW) strategy based on the gray entropy of the initialized image is studied. BCS-MVI achieves an average PSNR gain of 0.5–2.0 dB and an SSIM gain of 0.02–0.05 over previous block-based compressed sensing methods, and the reconstructing time only slightly fluctuates with the sampling rate. The algorithm is suitable for applications in multimedia data processing with fixed transmission delays. Keywords Mixed variational inequality . Image reconstruction . Block-based compressed sensing . Alternating direction method

1 Introduction Compressed sensing (CS) theory [4, 6] states that as long as a signal is sparse in a certain transform domain, a high-dimensional signal can be projected into a low-dimensional space through an observation matrix that is incoherent with a transform basis. The signal can be precisely reconstructed using a small number of observations. However, for large data (such as images and video), direct compressive sampling [6, 24] is very computationally expensive and extremely memory intensive. * Jian Chen [email protected] 1

College of Physics and Information Engineering, Fuzhou University, Fuzhou, China

2

School of Information Science and Engineering, Xiamen University, Xiamen, China

Multimed Tools Appl

To reduce the burden of sampling, LuGan proposed block compressed sensing (BCS) for a natural image [10] in 2007. In 2008, Stankovic et al. introduced orthogonal matching pursuit (OMP) [23] to reconstruct sparse blocks independently [21], and applied the BCS to video coding [22]. Although the use of BCS greatly improved the efficiency of the image and video sampling, it did not consider blocking artifacts. In 2009, Mun et al. combined the smooth projection algorithm and directional wavelet transform and suggested an algorithm called block compressed sensing and smoothed projected Landweber reconstruction (abbreviated as BCS-SPL) [19], which considered both the sparsity of transformation coefficients and smooth features of an image. The blocking artifacts were eliminated due to the frame-based recovering, but the recovery performance was still below the given requirements. In 2011, multi-scale BCS based on SPL (abbreviated as MS-BCS-SPL) [8] was proposed by Fowler et al. Although the recovery quality improved to a certain degree, the multi-scale sampling was rather complex. Later, the original BCS-SPL and MS-BCS-SPL were further extended to video coding [9], but their common and potential problem of original BCS-SPL and MS-BCS-SPL is that, they do not take image information by itself into account in the filtering process, which may make the image overly smooth and affect the reconstruction quality of texture characteristics. Besides, their reconstruction speed fluctuates with the threshold value and sampling rate. Generally speaking, the review of the block-based compressed sensing is summarized in the Table 1. To improve the recovery performance of current BCS for more types of nature images without increasing the sample complexity, a reconstruction algorithm based on the MVI unified framework [12, 15, 16] for BCS (abbreviated as BCS-MVI) as well as an adaptive adjustment strategy for the weights of the objective functions are proposed in this paper to rebuild natural images in high fidelity with low complexity within a reasonable amount of time.

2 Related works 2.1 Block-based sampling strategies Various strategies, such as adaptive compressive sampling, have been proposed to improve the performance of BCS. In 2013, Wang et al. proposed the block-based adaptive compressed sensing of an image using the texture information about the gray entropy [25]. In the same year, Zhang et al. assigned sampling rates based on each block’s texture complexity with respect to the gradient variance [29]. The reconstruction quality and visual effect of the above two adaptive BCS approaches were obvious improvements over earlier methods, but the sampling complexity was increased owing to the additional step of detecting texture information. In addition, a good measurement matrix, such as a measurement matrix based on a regression model, also helps improve the recovery quality. In 2015, Han et al. studied a novel generator of measurement matrixes for BCS, which was produced using the regression model between the coordinates of pixels and their gray levels [13]. The algorithm was speedy and accurate with good compliance, but the generator seemed more complex compared to that for conventional BCS. In this paper, conventional uniform sampling by block is applied to decrease the computational complexity. Moreover, the idea of adaptive sampling is applied to the recovery process to adjust the regularization weights.

MS-BCS-SPL [8, 9]

BCS-SPL [9, 19]

multi-scale sampling and frame-based recovering

compressive sampling and reconstructing data both by frame acquiring and recovering data both by block acquiring by block and recovering by frame

CS [6, 24]

BCS [21, 22]

Feature

Subject

Table 1 Review of the block-based compressed sensing

improving recovery quality to a certain degree

filtering without considering image itself; recovery with low quality and fluctuant speed filtering without considering image itself; sampling rather complex

image processing with unpredictable delays and regardless of complexity

image or video with smooth feature and processing with unpredictable delays

image or video in any size

image or video in small size

very computationally expensive and extremely memory intensive for large data

precisely reconstructing signal using a small number of observations greatly reducing the burden of sampling eliminating blocking artifacts blocking artifacts

Application

Disadvantage

Advantage

Multimed Tools Appl

Multimed Tools Appl

2.2 Reconstruction algorithms At present, the popular reconstruction algorithms for CS are optimization, thresholding, and greedy methods [7]. Greedy pursuits—such as OMP [23] and compressive sampling matching pursuit (CoSaMp) [20]—are easy to implement, but their reconstruction speeds depend on sparse degrees. The performance of thresholding-type algorithms—such as two-step iterative shrinkage/thresholding (TwIST) [5] and iterative hard thresholding (IHT) [3]—is not affected by the sparse degrees, but is greatly influenced by the threshold. Classical convex optimization algorithms—such as SOCP (second-order cone programming) [11] and L1-MAGIC—can achieve high reconstruction quality, but their complexity is also very high. A variety of reconstruction algorithms for convex optimization based on total variation (TV) minimization—such as NESTA (Nesterov’s algorithm) [2], TVAL3 (TV minimization by augmented Lagrangian and alternating direction algorithms) [17, 18], RecPF (reconstruction from partial Fourier data) [27], and IADPM (inexact alternating directions proximal method) [26]— have been developed in recent years and have greatly increased the speed of the convex optimization algorithm. For example, the NESTA algorithm could converge at a rate of O(1/k2), but the requirements for the measurement matrix were quite rigid. The TVAL3 method was applicable to numerous types of measurement matrices, and its recovery speed was faster than that of the conventional TValgorithm; however, this method considered only the sparsity of the image gradient. The RecPF algorithm had the advantages of low iteration cost and fast solution speed, but it was suitable only for partial Fourier matrixes; moreover, the weight of each regularization was fixed. The IADPM algorithm— because of its similar performance to TVAL3—was found to be suitable for signals influenced by a small amount of additive noise. The above optimization problems based on minimizing single or multiple objective functions could be incorporated into a unified framework of mixed variational inequality (MVI) [12, 15, 16]. He et al. proposed the MVI unified framework and proved that the worst convergence rate of the corresponding reconstruction algorithm was O(1/k), which provided a theoretical basis for the convergence rate of the above algorithms. In this paper, MVI with convergence is introduced to the recovery model of BCS for generality. Furthermore, because the above algorithms are applicable to objective functions with fixed weights, BCS-MVI with adaptive weights is also considered based on image characteristics to enable better adjustments to different images.

3 Preliminaries 3.1 Compressed sensing According to CS theory [4, 6, 7], a signal (u) that is sparse in the transform domain can be projected to the measured value f with a set of measurement vectors Φ that are incoherent to the transform basis (ψ): f ¼ Φu s:t: u ¼ ψθ

ð1Þ

where u is an N × 1 vector, f is an M × 1 measured value matrix, and Φ is an M × N measurement matrix (M < < N). Gaussian, Bernoulli, and scramble block Hadamard ensemble (SBHE) have been shown to be good choices for the measurement matrix Φ. θ is the sparse coefficient vector corresponding to the transform basis ψ.

Multimed Tools Appl

The problem of signal reconstruction is to solve u from formula (1). However, this is an illposed problem because M is far less than N. The problem can be solved by calculating the sparse coefficient vector from formula (2)—namely, solving the optimization problem with constraints under the lp-norm:







ð2Þ min

θ

lp s:t: f ¼ Φψθ θ

3.2 Block compressed sensing To alleviate the huge computation and memory burdens, some scholars have introduced the block-based sampling scheme into CS [8–10, 19, 21, 22]. In BCS, each frame is divided into NB non-overlapping blocks uj (sized B × B; subscript j denotes the block indicator) and acquired using a suitable MB × B2 measurement matrix ΦB; then, the corresponding fj is f j ¼ ΦB u j ; u j ∈RBB ; f j ∈RM B ; ΦB ∈RM B B

2

ð3Þ

It is straightforward to see that formula (3) applied block by block to an image is equivalent to taking measurements via formula (1) with a structure constraint that Φ be block diagonal: 2 3 ⋯ 0 ΦB 0 60 ΦB ⋯ 0 7 7 Φ¼6 ð4Þ 4⋮ ⋱ ⋮5 0 ⋯ 0 ΦB When the sparsity transform ψ is also a block-based operator, the image can be reconstructed by block at the decoding end. In general, the block-independent reconstruction will produce severe blocking artifacts, thus rebuilding by frame via formula (2)—such as BCS-TV—prior to by block. Because BCS-TV reconstruction requires heavy computations, Fowler et al. combined BCS with back projection and Wiener filtering and proposed an iterative thresholding algorithm: BCS-SPL [9, 19]. BCS-SPL algorithms based on the sparsity of discrete wavelet, dual-tree discrete wavelet and contourlet transforms—namely, BCS-SPL-DWT, BCS-SPL-DDWT and BCS-SPL-CT—are well known, efficient recovery algorithms for block-based CS.

3.3 MVI unified framework (1) The classical convex optimization problem The generic convex minimization problem with linear constraints is expressed by minfθðxÞjAx ¼ b; x∈X g

ð5Þ

, b∈R , and X⊂R are convex, and θ: R → R is a closed convex but not where A∈R necessarily smooth function. The objective function in (5) is generic. Its Lagrangian function is M×N

M

N

N

Lðx; λÞ ¼ θðxÞ−λT ðAx−bÞ where λ∈RM. Assuming Ω to be a closed convex subset, the standard form of MVI [12, 15, 16] for the convex optimization problem with linear constraints is    T   θðxÞ−θ x* þ ω−ω* F ω* ≥0

∀ω; ω*∈Ω

ð6Þ

Multimed Tools Appl

 T    x −A λ and Ω=X×Rm. The superscript B*^ represents the where ω ¼ , F ðωÞ ¼ λ Ax−b optimal solution. The unified approach based on MVI contains two steps: prediction and correction [12]. The iteration problem at the (k + 1)th step is           k k T k k T k ~ ≥ v−~v F ω Q vk −~v θðxÞ−θ ~x þ ω−~ ω ∀ω∈Ω and   ð7Þ k vkþ1 ¼ vk −C vk −~v where Q is a positive semidefinite matrix, ν is a partial parameter of ω, and C is a nonsingular matrix needed to satisfy the convergence conditions [12]. The superscript represents the number of iterations, and the B~^ symbol indicates the estimated value. (2) The structural convex optimization problem The unified framework applies to the objective function not only with a single variable in the classical optimization problem but also with multiple variables in the structural optimization problem. The structural optimization model with two separable variables is

 where x ¼

minfθ1 ðx1 Þ þ θ2 ðx2 ÞjA1 x1 þ A2 x2 ¼ b; x1 ∈X 1 ; x2 ∈X 2 g

ð8Þ

0 0 1    x1 x1 x 2 , θ(x)=θ1(x1)+θ2(x2), F ðωÞ ¼ @ , ω ¼ @ x2 A, v ¼ x2 λ λ

1 −A1 T λ T A −A2 λ A1 x1 þ A2 x2 −b

and Ω=X1 ×X2 ×Rm. The classical augmented Lagrangian and alternating direction algorithms (ALM-ADM) are an efficient scheme for solving this problem, and the (k + 1)th iterations are 8 2   β > > x1 kþ1 ¼ argmin θ1 ðx1 Þ−λT A1 x1 þ A2 xk2 −b þ A1 x1 þ A2 xk2 −b2 > > 2 x1 < 2   β ð9Þ x2 kþ1 ¼ argmin θ2 ðx2 Þ−λT A1 xkþ1 þ A2 x2 −b þ A1 xkþ1 þ A2 x2 −b2 > 1 1 > 2 > x 2 >   : λkþ1 ¼ λk −β A1 xkþ1 þ A2 xkþ1 1 2 −b According to references [12, 14], the above scheme can be represented as a prediction and correction scheme—a unified approach—and its worst convergence rate is O(1/k).

3.4 Efficient total variation ALM-ADM can be viewed as a special case of a unified approach, such as TVAL3, which is an efficient TV reconstruction algorithm. TVAL3 introduces a variable wi into the TV model and adds the constraint term wi =Diu to determine the discrete gradient of image u in the i-th pixel, where D is an operator used to calculate the gradient; therefore, the regularization model can be represented as min wi ;u

X

kwi kTV

s:t:

Φu ¼ f ; wi ¼ Di u

ð10Þ

Multimed Tools Appl

The corresponding augmented Lagrangian problem is  X β μ kwi kTV −υi T ðDi u−wi Þ þ i kDi u−wi k2 −λT ðΦu− f Þ þ kΦu−f k2 LA ðwi ; uÞ ¼ 2 2 i where υ and λ are Lagrange multipliers and β and μ are penalty parameters. ALM-ADM is applied to the recovery process. For the target function, the one-step steepest descent scheme and shrinkage-like formula are used to solve two sub-problems. The u sub-problem at the (k + 1)th iteration is ukþ1 ¼ uk −αk d k

ð11Þ

where d represents the gradient direction of the objective function and α represents the falling step. Barzilai and Borwein [1] suggested an aggressive manner to choose step length for the method, which is called the BB step or BB method. As can be seen, the BB step utilizes the previous two iterates and achieves the super-linear convergence [1]. However, the one-step steepest descent is not able to offer two iterates, Li C.B et al. [17, 18] derived the BB-like step, which leads to  k k−1 T  k k−1  u −u u −u ð12Þ αk ¼   ðuk −uk−1 ÞT d k −d k−1 To validate the BB-like step, a non-monotone line search algorithm (NLSA) is integrated to check the non-monotone Armijo condition [28]. When the non-monotone Armijo condition is not satisfied, backtracking is performed as αk ¼ ραk ; ρ∈ð0; 1Þ

ð13Þ

The iterative formula for the w sub-problem is   ! !   k  kþ1 k k   kþ1 υki  1 kþ1 kþ1 υi wi ¼ shrink Di u ; υi ; β i ¼ max Di u − k − k ; 0 sgn Di u − k  βi  βi βi Subsequently, the Lagrange multipliers and penalty parameters are updated as follows:   ¼ υki −β ki Di ukþ1 −wkþ1 υkþ1 i i   λkþ1 ¼ λk −μk Φukþ1 −f β kþ1 ¼ ηβ ki ; μkþ1 ¼ ημk ; η > 1 i In the following, the direct reconstruction via TVAL3 for the block-based measurement, which can achieve a similar recovery quality as BCS-TV with much lower complexity, is abbreviated as BCS-TVAL3.

4 Proposed method In the recovery process, the MVI framework is introduced into BCS for the structural optimization problem with multiple separable variables. As an improved algorithm for BCSTVAL3, BCS-MVI is used to illustrate the application of the unified approach. The basic

Multimed Tools Appl

process, reconstruction algorithm, and performance analysis for BCS-MVI are expanded upon below.

4.1 Basic flow (1) Scan by block First, an image I is divided into NB non-overlapping blocks Ij (sized m × n, subscript j denotes the block indicator): I j ¼ decompose ðI Þ

ð14Þ

Then, each m × n two-dimensional image block Ij is transformed into an mn column vector uj (j = 1, 2,..., NB) via an appropriate scanning method:   ð15Þ u j ¼ scan I j

(2) Measure by block Second, each column vector of an image block is measured via a normalized random matrix ΦB (sized MB × mn), and NB columns of observations fj (j = 1, 2,..., NB) are obtained. f j ¼ ΦB u j

ð16Þ

 All observations are combined into a measured value matrix F ¼ f 1 ; :::; f N B , which is equivalent to taking measurements with the matrix of the original image U ¼ ½u1 ; :::; uN B ; specifically, F ¼ ΦU

(3) Recovery by frame Third, the sparsity of the gradient and wavelet domain of the entire image is taken as prior knowledge, and the image matrix is reconstructed via the function below: ~ ¼ recoverð F; ΦÞ U

(4) Iscan and integrate Then, each column vector of the restored image block is inversely scanned and converted into the image block Ĩ j:

 ~I j ¼ iscan u~ j

Multimed Tools Appl

Finally, they are integrated into a complete image Ĩ :

 ~I ¼ compose ~I j

4.2 Reconstruction model The function recover(F, Φ) represents the reconstruction algorithm of BCS-MVI. The mixed regularization model is constructed with the constraint condition concerning reconstruction accuracy as well as the characteristics of the local smoothness in the spatial domain and the sparsity in the wavelet transform domain of the image.   min τ W kDU kTV þ τ Z ψT U 1 st: ΦU ¼ F U

where ‖DU‖ TV is the abbreviation of ∑ kDi ukTV , ‖ψ T U‖ 1 is the 1-norm of the i

wavelet coefficient of a signal, and τ W and τZ represent the weights for the TV and DWT terms, respectively. Φ includes scanning and measuring by block using formulae (14)–(16), which require fewer computations and reduced memory burden compared to direct measurement using conventional CS. The parameters W and Z are introduced to substitute the variables in the model such that W=DU and Z=ψTU. The problem is thus transformed into an optimization problem with constraints min τ W kW kTV þ τ Z kZ k1 st:

W ;Z;U

ΦU ¼ F; W ¼ DU; Z ¼ ψT U

ð17Þ

There are three variables in formula (17): U, W and Z. Because W and Z are mutually independent variables, they can be viewed as one variable, namely, x2. By defining x1 = U and x2 = (WT, ZT)T, the problem can be simplified into a structural optimization problem with two separable variables. Compared to the linear-constrained convex optimization problem with two separable variables in the unified framework (6), the objective function and constraint parameter are 2 3 2 3 2 3 F Φ 0 0 θ1 ðx1 Þ ¼ 0; θ2 ðx2 Þ ¼ τ W kW kTV þ τ Z kZ k1 ; A1 ¼ 4 D 5; A2 ¼ 4 −I 0 5; b ¼ 4 0 5 0 ψT 0 −I After the modification from the original optimization problem to the standard form of the MVI frame, the relevant parameters could satisfy the convergence conditions and the corresponding unified approach would have a worst convergence rate. That’s why the modified ALM-ADM algorithm below, as a special case of the unified approach, can recovery image with more stable time. 2 3 2 3 βU λU Assuming β ¼ 4 β W 5; λ ¼ 4 λW 5, the (k + 1)th iterations of the modified ALMβZ λZ ADM via formula (9) based on the unified approach are listed as follows:

Multimed Tools Appl

(1) The objective problem, including x1 and x2 sub-problems 8 2  T  T   β  β > > U kþ1 ¼ argmin − λkU ðΦU −F Þ þ U kΦU − F k22 − λkW DU−W k þ W DU−W k 2 > > > 2 2 U > > > 2  T  T  β  > > < − λkZ ψ U −Z k þ Z ψT U −Z k 2 2 2   β    T > > W kþ1 ¼ argmin τ W kW kTV − λkW DU kþ1 −W þ W DU kþ1 −W 2 > > 2 > W > > 2  β   T  > > > Z k ¼ argmin τ Z kZ k1 − λkZ ψT U kþ1 −Z þ Z ψT U kþ1 −Z 2 : 2 Z

(2) The dual problem, including 3 components  kþ1  8 kþ1 k < λU ¼ λU −βU ΦU − F  λkþ1 ¼ λkW −β W DU kþ1 −W kþ1   : Wkþ1 λZ ¼ λkZ −β Z ψT U kþ1 −Z kþ1

ð18Þ

ð19Þ

4.3 Reconstruction algorithm The reconstruction algorithm based on the MVI frame is explained below. Algorithm 1 takes a fixed weight (FW) for the objective function. Algorithm 1 Function BCS-MVI-FW Input: F, Φ; Output: U Initialize: U0 = ΦTF Iterations: While (stopping criteria unsatisfied) Do the (k+1)th iterations 1) the objective problem (1)the U sub-problem: ① compute BB-like step αk via formula (12) ② verify and update step via formula (13) when the non-monotone Armijo condition is not satisfied compute ③ compute Uk+1 via formula (11) (2) the W and Z sub-problems:   ! !  k λk   1  kþ1 λW − W kþ1 ¼ max DU kþ1 − W ; 0 sgn DU −   βkW  βkW βkW   ! !  k λ 1 λk   Z kþ1 ¼ max ψT U kþ1 − Zk − k ; 0 sgn ψT U kþ1 − Zk   βZ βZ βZ

(3) update parameters: Update the objective function and gradient value 2) the dual problem Update the Lagrange multipliers via formula (19) 3) update the penalty parameters   βkþ1 ¼ min ηβk ; β max End Do

Multimed Tools Appl

The solutions of sub-problems U and W and the update of the penalty parameters are similar 2 k 3 2 3 ðβ max ÞU βU to the TVAL3 algorithm [17, 18], where β k ¼ 4 β kW 5 has upper limit β max ¼ 4 ðβ max ÞW 5; ðβmax ÞZ β kZ that is,   βkþ1 ¼ min ηβ k ; β max k kþ1 ¼ ηβkU ; β kþ1 ¼ ηβ kZ ; β < βmax β kþ1 U W ¼ ηβ W ; β Z ¼ kþ1 kþ1 kþ1 β U ¼ ðβ max ÞU ; β W ¼ ðβ max ÞW ; β Z ¼ ðβ max ÞZ β ≥β max

Algorithm 1 adds a regularization of the sparsity in the wavelet domain to TVAL3. This algorithm can achieve good results for textured images. However, for images with smooth features, the reconstruction quality improves to a minimal extent, whereas the complexity significantly increases after wavelet regularization is included. As is well known, the firstorder gray entropy of an image reflects its texture information; we find that the entropy of the initialized image can maintain the textural properties similar to those in the original. Hence, we estimate the texture characteristics by the gray entropy of the middle data, although it is not accurate. Based on this idea, Algorithm 2 takes the adaptive weight (AW) for the target function according to the gray entropy from the initialized U0. Algorithm 2 Function BCS-MVI-AW Input: F, Φ; Output: U Initialize: 1) initialize U: U0 = ΦTF 2) calculate the gray entropy: entropy = entropy(U0 ) Iterations: While (stopping criteria not satisfied) Do if β k <βmaxthen assign the weights: [τW,τZ]= [1,0] Do the same iterations as those in Algorithm 1 else assign the weights: [τW,τZ]= thresh(entropy) Do the same iterations as those in Algorithm 1 end End Do

When βk is small, the iteration results are not exact. We can set [τW,τZ] as [1,0] to alleviate the calculational burden because sub-problem W is easier to resolve than is sub-problem Z. The accuracy requirement increases with increasing βk. The weights for the two regularization terms are set via a threshold comparison method depending on the gray entropy of the initialized image U 0 . For example, given the selectable weights for τ W (1 > τ1 > τ2 > τ3 > 0) and thresholds of gray entropy (th1 < th2 < th3), the weights should be 8 ½1; 0 entropy < ¼ th1 > > <  ½τ 1 ; 1−τ 1  th1 < entropy <¼ th2 τ W ; τ Z ¼ threshðentropyÞ ¼ ½τ ; 1−τ 2  th2 < entropy <¼ th3 > > : 2 ½τ 3 ; 1−τ 3  entropy > th3

Multimed Tools Appl

The smaller the entropy, the smoother the image. In this case, the weight of the TV, τW, should be larger, and that of DWT, τZ, should be smaller, and vice versa.

4.4 Performance analysis Taking BCS-TVAL3 and BCS-SPL as references, the performance of BCS-MVI is analyzed in this section. Because the above algorithms adopt uniform sampling by block, the sampling complexities are equal. Their main differences lie in the reconstruction phase. Because BCSTVAL3 contains only the TV regularization, the method may obtain good results for images with the characteristic of piecewise smoothness. Moreover, owing to the simple operation of the transform and inverse transform of the gradient, its processing speed is high, but the effect is not as good for textured images. BCS-SPL-DWT would process threshold shrinkage in the wavelet domain and Wiener filtering in the pixel domain at each iteration, which considers both the texture information and smooth features of an image. However, the method does not take image information by itself into account in the filtering process. Instead, using a 3 × 3 Wiener filtering template, it may make the image overly smooth and affect the reconstruction quality. Its reconstruction speed fluctuates with the threshold value and sampling rate. BCSMVI-FW adds wavelet regularization to BCS-TVAL3 and improves the reconstructing effect for textured images. This alternately minimizes the augmented Lagrangian function of the deviation, gradient and wavelet regularization at each iteration; therefore, BCS-MVI-FW would achieve the best quality, but its speed is very low because it considers additional factors and has a higher complexity. Because the gradient operation is very simple, the reconstruction speed of BCS-MVI-FW depends on that of DWT and IDWT. To achieve a balance between quality and speed, BCSMVI-AW allocates the weight to each regularization of the objective function according to the estimated gray entropy in the iteration of BCS-MVI-FW. For a smooth image, BCS-MVI-AW sets τW close to 1, and its reconstruction speed is similar to that of BCS-TVAL3 at the lowest complexity. BCS-TVAL3 is a special case of BCS-MVI-FW when the weights are τW = 1 and τZ = 0. For a textured image, τZ should be set between 0 and 1; therefore, the restoration speed is between that of BCS-TVAL3 and BCS-MVI-FW. Of course, BCS-SPL-DWT has been expanded to include the dual-tree discrete wavelet and contourlet transform (−DDWT and -CT) to provide better results. Those transformations can also be introduced to the regularization model of BCS-MVI to further improve recovery performance.

5 Experiments To validate the performance of BCS-MVI, we performed various simulations using Matlab and measured the peak signal-to-noise ratio (PSNR) and structural similarity (SSIM) to compare the reconstruction quality and reconstruction time (Time) to determine the computational complexity. The original image was first divided into 32 × 32 blocks (namely, m = n = 32); scanned by column; measured via a normalized random Gaussian matrix in the block structure, where the sampling number of each block is MB = m × n × sub-rate; and finally reconstructed via BCSMVI-FW and BCS-MVI-AW. The original materials were taken from four classical gray images: Lenna (with piecewise smooth and various features), Couple (with coarse texture features), Barbara (with exquisite texture features) and Boat (with cartoon texture features).

Multimed Tools Appl

Each image is 256 × 256 pixels. The experimental platform was a Lenovo notebook configured with an Intel Core i5-2520 M CPU at dominant frequency 2.50 GHz and 3.05 GB memory. To enable a convenient comparison, we also applied three efficient algorithms—BCS-SPLDDWT, BCS-SPL-CT and BCS-TVAL3—as references. The two weights τW and τZ for BCSMVI-FW are both fixed at 0.5, and the weights for BCS-MVI-AW are adaptively adjusted according to the strategy described in Section 4.3. Figures 1, 2, 3, 4, 5, 6 and 7 compare the objective qualities, including PSNR, SSIM and time, as well as the subjective quality of the reconstructed images. In addition, Table 2 shows the average objective qualities for various materials, where the mean and gain of PSNR, SSIM and time for the BCS-SPL and BCS-MVI algorithms are calculated. Figures 1, 2 and 3 compare the recovery quality and time costs for five reconstruction algorithms under different measurement rates (namely, sub-sampling rates, abbreviated as subrates). We can see that the PSNRs of BCS-TVAL3 and BCS-SPL (−DWT or -CT) have their own advantages for smooth and textured images. The SSIM for BCS-TVAL3 is superior to that of BCS-SPL, whereas BCS-MVI (−FW or -AW) achieves the optimal PSNR and SSIM among the test images. BCS-TVAL3 requires the shortest time—usually 3–4 s—whereas the reconstruction time for the BCS-SPL algorithms vary greatly with the sampling rate. BCS-

Fig. 1 Performance comparison in terms of PSNR

Multimed Tools Appl Table 2 Mean and gain of objective performance Test Images Lenna

Recovery BCS-SPL- BCSPerformance DDWT SPL-CT

BCSBCSAverage MVI-FW MVI-AW BCS-SPL

Average BCS-MVI

PSNR (dB)

33.4239

33.4261

SSIM

0.89249

Time (s) Couple

PSNR (dB) SSIM Time (s) SSIM Time (s)

0.8797

8.0816

8.2847 35.1287

0.90546 9.8316 36.611

33.4283 0.90543 5.2917 36.4705

0.894

0.9189

6.9306

9.0069

10.7222

5.1823

30.3973

30.3785

31.0501

30.9535

0.84326 7.7465

PSNR (dB)

32.2347

35.2155 0.89388

Barbara PSNR (dB)

Boat

32.961

29.2274

0.83947 8.1927 28.9323

SSIM

0.82052

0.81169

Time (s)

8.0035

7.0885

0.86131 9.5174 31.1268 0.86718 10.151

0.91874

0.86015 7.224 31.058

32.59785 0.886095 8.18315 35.1721 0.89394 7.96875 30.3879 0.841365 7.9696 29.07985

0.905445 7.56165

Average Gain 0.82825 0.01935 −0.6215

36.54075

1.36865

0.91882

0.02488

7.95225 31.0018 0.86073 8.3707 31.0924

−0.0165 0.6139 0.019365 0.4011 2.01255

0.86439

0.816105

0.865785

0.04968

6.026

7.546

8.0885

0.5425

MVI-FW reconstructs images in the longest time but slowly fluctuates, whereas the reconstruction time of BCS-MVI-AW is between that of BCS-TVAL3 and BCS-MVI-FW and provides certain advantages for real-time processing and reconstruction time estimation. Some

(a) Lenna

(b) Couple

(c) Barbara

(d) Boat

Fig. 2 Performance comparison in terms of SSIM

Multimed Tools Appl

Fig. 3 Performance comparison in terms of recovery time

may wonder why BCS-MVI-AW achieves a performance similar to that of BCS-MVI-FW in a shorter recovery time. The reasons are as follows: Most natural images are piecewise smooth and sparse in both TV and DWT domains; thus, the reconstruction quality via the two domains will be similar. However, DWT (via multiplication) is more complex than TV (via addition and subtraction); thus, it should require a longer recovery time when more weights are allocated to DWT. For smooth images whose gray entropy is smaller, BCS-MVI-AW allocates more weights to TV than to DWT. Therefore, the performance levels of BCS-MVI-FW and BCSMVI-AW are very similar, but BCS-MVI-AW exhibits the shorter recovery time. Because 50 % is the intermediate value of the sub-rates ranging from 10 to 90 %, we focus on the recovered image at this sub-rate to estimate the average subjective performance of the proposed algorithm. Figures 4, 5, 6 and 7 show the comparison between the original and reconstructed images for the above four test images at a sub-rate of 50 %.. The above four pictures show that the reconstruction effect of BCS-TVAL3 is good for images with piecewise smoothness and simple texture features, such as Lenna, Couple and Boat; however, the details are dull for complex texture images such as Barbara. BCS-SPL (−DDWT and -CT) can recover images with better layering and contrast, but it struggles with some noise particles. The restored images for BCS-MVI (−FW and -AW) are gentler and layered with less noise, and the definition and quality levels in the subjective sense are better than those of other algorithms.

Multimed Tools Appl

Fig. 4 Subjective effect comparison (for Lenna)

To quantitatively compare the average objective performance between the popular BCSSPL and the proposed BCS-MVI, Table 2 lists the average PSNR, SSIM and recovery time for

Fig. 5 Subjective effect comparison (for Couple)

Multimed Tools Appl

Fig. 6 Subjective effect comparison (for Barbara)

two types of BCS-SPL and BCS-MVI algorithms at sub-rates from 10 to 90 %. Because BCSTVAL3 is a special case of BCS-MVI-FW, we omit it for brevity..

Fig. 7 Subjective effect comparison (for Boat)

Multimed Tools Appl

In the above table, the highest PSNR and SSIM as well as the shortest recovery time are in boldface, highlighting the optimal performances among the BCS-SPL-DDWT, BCS-SPL-CT, BCS-MVI-FW and BCS-MVI-AW algorithms. Whereas BCS-MVI-FW almost attains the highest PSNR and SSIM, the reconstruction quality of BCS-MVI-AW is similar to that of BCS-MVI-FW while requiring the shortest recovery time and thus is more efficient than the other methods. The average performance of BCS-SPL is the mean performance of BCS-SPL-DDWT and BCS-SPLCT, the average performance is BCS-MVI is that of BCS-MVI-FW and BCS-MVI-AW, and the average gain is the difference between them. It can be seen that BCS-MVI achieves an average PSNR gain of 0.5–2.0 dB and SSIM gain of 0.02–0.05 over the BCS-SPL. Furthermore, noisy images such as Peppers and Goldhill corrupted by Gaussian noise are also tested at a low sub-rate (e.g., 30 %), where the mean of the noise is that of the image gray, and the noise level represents its variance (e.g., 0.05). Figures 8 and 9 show that the recovery performance of BCS-MVI-FW and BCS-MVI-AW are still superior to that of BCS-SPL-DDWT and BCSSPL-CT in noisy environments. As a special case of BCS-MVI-FW, BCS-TVAL3 also performs well and provides the advantage of high speed. BCS-MVI-AW appears to be most valuable because it achieves trade-offs between quality and speed. The flowing two figures show preliminary evidence of the algorithm’s performance for noisy images. More experiments will be performed and discussed in future studies due to space limitations.

6 Conclusions and future work BCS-MVI is proposed as a strategy to improve the performance of current reconstruction algorithms for BCS under the condition of low complexity. For encoding, an image is divided

Fig. 8 Performance comparison (for Peppers)

Multimed Tools Appl

Fig. 9 Performance comparison (for Goldhill)

into several non-overlapping blocks, which are scanned by column and then measured by block to obtain column vectors. For decoding, the column vectors corresponding to the measured values are integrated into a matrix, and the augmented Lagrangian functions of the deviation, gradient and wavelet regularization are alternately minimized, assuming prior knowledge of the sparsity in the gradient and wavelet domains of the image. The regularization can be assigned either a fixed or an adaptive weight. This new approach introduces an MVI frame into the reconstruction model of BCS as well as a reconstruction algorithm for BCS based on the unified approach to ensure convergence. Furthermore, this approach can estimate the image texture according to the gray entropy from the initialized data and adaptively choose the weights for different regularizations to improve the recovery efficiency. In contrast to other reconstruction algorithms for BCS-SPL, our simulations show that the BCS-MVI algorithm achieves an average PSNR gain of 0.5– 2.0 dB and SSIM gain of 0.02–0.05 over BCS-SPL, and the reconstruction time only slightly fluctuates with the sampling rate; this is especially favorable for multimedia data processing with a fixed transmission delay. The use of the integral gray entropy of the intermediate image to estimate the texture is suitable only for images with consistent textures. Thus, for images with large texture variations among blocks, this algorithm benefits from the local optimization of different texture layers. Therefore, future research will focus on a BCS-based reconstruction algorithm for images with uneven textures. In addition, the thresholds of the gray entropy are currently empirical values and are empirically determined, which should be further improved. Certainly, there are other strategies for weight decisions, such as the use of gradient entropy, gradient variances, and gray-level cooccurrence matrixes. Other methods should be compared in the future.

Multimed Tools Appl Acknowledgments The paper is supported by the National Natural Science Foundation of China (No. 61170147 and 61471124), and the Natural Science Foundation of Fujian Province (No. 2013 J01234, 2014 J01234 and 2015 J01251). The wonderful lectures and patient help of Prof. Peng Zheng in the College of Mathematics and Computer Science at Fuzhou University is greatly appreciated, as are the wonderful lectures and scholarly communications from Prof. He Bing-sheng and Mrs. Tao Ming at Nanjing University.

References 1. Barzilai J, Borwein JM (1988) Two-point step size gradient methods. IMA J Numer Anal 3:141–148. 2. Becker BJ, Candes E (2009) NESTA: a fast and accurate first-order method for sparse recovery. SIAM J Imag Sci 4(1):1–39. 3. Blumensath T, Davies ME (2009) Iterative hard thresholding for compressed sensing. Appl Comput Harmon Anal 27(3):265–274. 4. Candes EJ, Romberg J, Tao T (2006) Robust uncertainty principles: exact signal reconstruction from highly incomplete frequency information. IEEE Trans Inf Theory 52(2):489–509. 5. Dias JB, Figueiredo M (2007) A new TwIST: two-step iterative shrinkage/thresholding algorithm for image restoration. IEEE Trans Image Process 16(12):2992–3004. 6. Donoho DL (2006) Compressed sensing. IEEE Trans Inf Theory 52(4):1289–1306. 7. Foucart S, Rauhut H (2013) A mathematical instruction to compressive sensing. Springer New York Press, USA. 8. Fowler JE, Mun S, Tramel EW (2011) Multiscale block compressed sensing with smoothed projected landweber reconstruction. In Proc. European Signal Processing Conference (EUSIPCO), Barcelona, Spain, 2011: 564–568 9. Fowler JE, Mun S, Tramel EW (2012) Block-based compressed sensing of images and video. Found Trends Signal Process (4):297–416. 10. Gan L (2007) Block compressed sensing of natural images. In Proc. Digital Signal Processing (DSP), Cardiff, UK, 2007:403–406 11. Goldfarb D, Yin W (2005) Second-order cone programming methods for total variation based image restoration. SIAM J Sci Comput 27(2):622–645. 12. Gu GY, He BS, Yuan XM (2014) Customized proximal point algorithms for linearly constrained convex minimization and saddle-point problems: a unified approach. Comput Optim Appl 59:135–161. 13. Han H, Gan L, Liu SJ, Guo YY (2015) A novel measurement matrix based on regression model for block compressed sensing. J Math Imaging Vision 51:161–170. 14. He BS, Yuan XM (2012) On the O(1/t) convergence rate of the douglas-rachford alternating direction method. SIAM J Numer Anal 50:700–709. 15. He BS, Liao LZ, Wang X (2012a) Proximal-like contraction methods for monotone variational inequalities in a unified framework I: effective quadruplet and primary methods. Comput Optim Appl 51:649–679. 16. He BS, Liao LZ, Wang X (2012b) Proximal-like contraction methods for monotone variational inequalities in a unified framework II: General methods and numerical experiments. Comput Optim Appl 51:681–708. 17. Li CB (2009) An effcient algorithm for total variation regularization with applications to the single pixel camera and compressive sensing. M.A. Rice University, USA. 18. Li CB (2011) Compressive sensing for 3D data processing tasks: applications, models and algorithms. USA. Ph.D. Rice University. 19. Mun S, Fowler JE (2009) Block compressed sensing of images using directional transforms. In Proc. Int. Conf. on Image Processing(ICIP), Cairo, Egypt, 2009 20. Needell D, Tropp J (2008) CoSaMP: Iterative signal recovery from incomplete and inaccurate samples. Appl Comput Harmon Anal 26(3):301–321. 21. Stankovic V, Stankovic L, Cheng S (2008a) Compressive sampling of binary images. In Proc. Congress on Image and Signal Processing (CISP), Sanya, Hainan, China, 2008 22. Stankovic V, Stankovic L, Cheng S (2008b) Compressive video sampling. In Proc. European Signal Processing Conference (EUSIPCO), Switzerland, Lausanne, 2008 23. Tropp J, Gilbert A (2007) Signal recovery from random measurements via orthogonal matching pursuit. IEEE Trans Inf Theory 53(12):4655–4666. 24. Wakin MB, Laska JN, Duarte MF, Baron D, Sarvotham S, Takhar D, Kelly KF, Baraniuk RG. (2006) Compressive imaging for video representation and coding. In Proc. Picture Coding Symposium (PCS), Beijing, China, 2006

Multimed Tools Appl 25. Wang RF, Jiao LC, Liu F, Yang SY D2013] Block- based adaptive compressed sensing of image using texture information. Chin J Electron 41D8]:1506–1514Din Chinese]. 26. Xiao YH, Song HN (2012) An inexact alternation directions algorithms for constrained total variation regularized compressive sensing problems. J Math Imaging Vision 4:114–127. 27. Yang J, Zang Y, Yin W (2010) A fast alternating direction method for TVL1-L2 signal reconstruction from partial Fourier data. IEEE J Sel Top Signal Process 4(2):288–297. 28. Zhang H, Hager WW (2004) A nonmonotone line search technique and its application to unconstrained optimization. SIAM J Optim 14:1043–1056. 29. Zheng HB, Zhu XC (2013) Sampling adaptive block compressed sensing reconstruction algorithm for images based on edge detection. J China Univ Posts Telecommun 20(3):97–103.

Kaixiong Su was born in 1959 in Fuzhou, Fujian, China. Professor at the College of Physics and Information Engineering, Fuzhou University. Doctoral tutor. He received his master degree in 1988 from the University of Science and Technology of China. His Interesting covers wireless communications, video broadcasting and embedded systems.

Jian Chen was born in 1981 in Fuzhou, Fujian, China. Lecturer at the College of Physics and Information Engineering, Fuzhou University. She received her master degree in communication and information system in 2007 from Fuzhou University,China. She is currently pursuing Ph.D. degree in Fuzhou university and Interesting covers image processing, video coding and compressed sensing.

Multimed Tools Appl

Weixing Wang was born in 1959 in Shuangfeng, Hunan, China. Professor at the College of Physics and Information Engineering, Fuzhou University. Doctoral tutor. He received his Ph.D. degree in 1997 from Royal Institute of Technology of Switzerland. His Interesting covers image processing, pattern recognition and machine vision.

Lichao Su was born in 1989 in Fuzhou, Fujian, China. He received his Master degree of engineering in computer science and technology in 2014 from Fujian Normal University, China. He is currently pursuing Ph.D. degree in Xiamen University and his research efforts are mainly focused on data mining, multimedia forensics.