# Statistical Prior Aided Separate Compressed Image Sensing for Green Internet of Multimedia Things

Authors: Shaohua Wu; Tiantian Zhang; Jian Jiao; Jingran Yang; Qinyu Zhang
Journal: Mobile Information Systems
 (2017)

---

## Abstract

In this paper, we aim to propose an image compression and reconstruction strategy under the compressed sensing (CS) framework to enable the green computation and communication for the Internet of Multimedia Things (IoMT). The core idea is to explore the statistics of image representations in the wavelet domain to aid the reconstruction method design. Specifically, the energy distribution of natural images in the wavelet domain is well characterized by an exponential decay model and then used in the two-step separate image reconstruction method, by which the row-wise (or column-wise) intermediates and column-wise (or row-wise) final results are reconstructed sequentially. Both the intermediates and the final results are constrained to conform with the statistical prior by using a weight matrix. Two recovery strategies with different levels of complexity, namely, the direct recovery with fixed weight matrix (DR-FM) and the iterative recovery with refined weight matrix (IR-RM), are designed to obtain different quality of recovery. Extensive simulations show that both DR-FM and IR-RM can achieve much better image reconstruction quality with much faster recovery speed than traditional methods.


---

## Body

## 1. Introduction

As a means of connecting the proliferating embedded devices to the Internet, the Internet of Things (IoT) has drawn extensive attention in recent years. IoT has the great potential to significantly influence our lives and the way we interact with devices including sensors, security/surveillance cameras, mobile phones, home automation/control devices, and vehicle state meters [1]. Different from the traditional wireless sensor networks (WSNs), IoT is envisioned to have a much larger deployment scale and cover a much broader range of services and applications. Therefore, the critical issue of efficient energy utilization in the design of wireless network transmission technologies, including the WSNs, LTE networks [2], cognitive radio networks [3], interference-alignment based networks [4], becomes much more severe for IoT, especially for IoT specialized for multimedia services and applications, which is referred to as “Internet of Multimedia Things” (IoMT) [5].The multimedia content, for example, image, audio, video, acquired from the physical environment possesses distinct characteristics as compared with the scalar data acquired by typical IoT devices, so the IoMT devices require much higher processing and memory resources. Moreover, the multimedia transmission is more bandwidth hungry as compared with the conventional scalar data traffic in IoT. Considering that the devices are usually low-computation capability devices and may run on battery with limited power, there is a great need for developing special multimedia data processing and transmission technologies that enable green computation and communications, that is, low-power and low-complexity signal processing, and energy efficient transmission. A core module of the multimedia data processing is compression, which aims to reduce the amount of data needed to be transmitted. Traditional compression techniques, such as the JPEG, JPEG 2000 algorithm for image compression, the MPEG, and H.26X algorithm for video compression, are generally too complex for low-cost IoMT devices. In this paper, we focus on the compressed sensing (CS) based image compression and reconstruction methods, aiming to design low-complexity and robust image transmission strategies for green IoMT applications.The CS technique was first proposed by Candes et al. [6, 7], in which the original data can be accurately reconstructed from only a portion of the sampled data, sampled at rates lower than the Nyquist rate. The key idea of CS is to make use of the sparse or compressible nature of the original data to reduce the size of the sampled/transmitted/stored data. Multimedia data, including images, generally holds this sparse or compressible nature. For example, images are approximately sparse (or compressible) in the wavelet domain. Since its invention, CS has sparked tremendous research interests from related areas such as medical imaging, wireless communication, and image/video compression. This technique lends itself naturally to the IoMT domain, and CS based image compression/reconstruction has become a very hot research topic. In [8], a novel CS based data aggregation mechanism tailored to the application of large-scale air-pollution monitoring with IoT devices was proposed. Their design exploits both the intra- and intersource correlations among air-pollution data using the framework of compressed sensing with side information. In [9], by applying the dictionary learning based denoising method within the approximate message passing (AMP) algorithm framework, a mechanism for efficient communication in IoMT was proposed.CS based image compression and reconstruction have already shown great potential for IoMP applications. However, there still remain two main problems. First, although the system complexity of “compression-reconstruction” has been greatly shifted from the encoder-end to the decoder-end, the encoding complexity is still too high to be handled by IoMT device when the size of the image is large. Second, although it is commonly assumed that the server node has abundant computation capability for the CS image reconstruction, there exist indeed situations where the decoder-end cannot afford to reconstruct the image, especially for real-time multimedia applications, for example, video surveillance using mobile handset to collect the IoMT packets. For the first problem, one solution is to divide the large image into small blocks and then compress each block sequentially. This method prominently reduces the computation complexity, but block processing destroys the boundary coherence between pixels, leading to obvious reconstruction quality degradation. To overcome this disadvantage, some new methods have been proposed, among which the separable image sensing encoder is the one that not only has affordable computation complexity, but also promises good reconstruction quality [10]. By separable sensing, the original 2D image is directly measured first in the row direction and then in the column direction, instead of being extended to a 1D signal. For reconstruction, one can transform the 2D image recovery problem to a 1D signal recovery by using of Kronecker product [11] or use the 2-dimensionality orthogonal matching pursuit (2D-OMP) algorithm proposed by [12]. However, both the Kronecker reconstruction method and the 2D-OMP algorithm are very computation-consuming; that is, the second problem we note above still exists with these methods.In our previous work, a separate-combine recovery has been proposed to yield a highly efficient recovery technique for CS image processing [13]. It can process large-sized images that can hardly be handled by existing methods such as the Kronecker reconstruction method or the 2D-OMP algorithm. The expense of the separate-combine recovery method is a slight degradation of the reconstruction quality. Moreover, its computation complexity is still high for practical real-time applications.In this paper, we aim to propose a reconstruction method that can further reduce the computation complexity and, meanwhile, increase the reconstruction quality. To achieve this goal, statistical information is considered in the recovery method design. By analyzing the image representations in the wavelet domain, the energy distribution of natural images is found to be well fitted by an exponential decay model, which is then used as statistical a priori information in the recovery method. Based on the separate compression image sensing, the recovery process is also composed of two steps: the row-wise (or column-wise) intermediates recovery and the column-wise (or row-wise) final results recovery. The aim of reconstruction in each step is to obtain a solution conform with the statistical prior by introducing a weight matrix. Two recovery strategies are, respectively, designed in this paper: the Direct Recovery strategy with Fixed weight Matrix (DR-FM) and the Iterative Recovery strategy with Refined weight Matrix (IR-RM). The DR-FM strategy is aimed to achieve extremely high computation efficiency by employing the same weight matrixes for the two recovery steps, while, for the IR-RM strategy, the weight matrix for the second step will be iteratively refined, aiming to achieve more accurate recovery results. Extensive simulations and experiments have been conducted, as results show that the proposed recovery method with DR-FM strategy has a much better computational efficiency than traditional methods, while the recovery quality is better; the recovery with IR-RM strategy can achieve the best quality of recovery at the expense of degradation in computational efficiency slightly, yet still faster than the traditional methods.The remainder of the paper is organized as follows. Section2 introduces the preliminaries and related work, including the CS basics and its applications to IoMT image transmission. Section 3 describes the separate image sensing and reconstruction methods. In Section 4, the proposed recovery method aided by statistical prior is described in detail. Section 5 presents the simulation results, followed by conclusions in the last section.


## 2. Preliminaries and Related Work

In this section, we summarize the preliminaries and related work of this paper, including the CS basics, the image transmission in IoMT, and prior aided CS image reconstruction.


### 2.1. The CS Basics

CS is applied for signal which is sparse or compressible. Under certain basesΨ=[ψ1,ψ2,…,ψn], a 1D signal x can be expressed as x=Ψs. We say x is c-sparse if it satisfies(1)x0=cor(2)s0=c.A signal is said to be nearly sparse or compressible if the largest c coefficients contribute to most of the signal energy, and the n-c others are small enough to be ignored. Signals in practical systems, including the images, are mostly compressible in a certain domain, so the CS theory can be used widely in natural signal processing.In CS framework, a spare or a nearly sparse signalx (size: n×1) is measured by a random matrix Φ (size: m×n, where m≪n) to generate a compressed result y (size: m×1):(3)y=Φx=ΦΨs=As,where A=ΦΨ.The reconstruction ofx corresponds to the problem of solving n unknown parameters from m linear equations, that is, a calculation for an underdetermined system of equations which contains innumerable solutions. Fortunately, s, the representation of the original signal under certain basis, is sparse enough, and the solution of the underdetermined system of equations can be achieved by solving the combinatorial optimization problem below:(4)s^=arg⁡min⁡s0,s.t.y=As.If the measurement matrixΦ is incoherent with the representation basis Ψ, especially when Φ is a random matrix and m=o(clog⁡(n/c)), A satisfies the restricted isometry property (RIP) with enormous probability [14]. The l0-norm optimization problem in (4) can be replaced by the l1-norm optimization problem which has the same unique solution; that is,(5)s^=arg⁡min⁡s1,s.t.y=As.The original signalx can then be obtained as(6)x^=Ψs.In conclusion, the crux of CS theory is to reconstruct the sparse representation of the original signal by making good use of the sparsity nature. There are mainly two kinds of reconstruction algorithms.


#### 2.1.1. Greedy Algorithms

They are a kind of recovery methods that can find out the locally optimal solution but not the globally optimal solution. Orthogonal matching pursuit (OMP) is a classical greedy method for CS recovery, and it is the base of many later proposed advanced methods, such as CoSaMP [15] and A∗OMP. The crux of greedy algorithms is to find out columns of the matrix A that make significant contributions to the compressed result iteratively.


#### 2.1.2. Convex Optimization Recovery

Different from the greedy algorithms, the convex optimization recovery methods aim to approach the globally optimal solution. The results by convex optimization recovery are usually more accurate, but the expense is more computation resources needed for the reconstruction. Convex optimization recovery can be used in situations where the recovery quality is preferred more than the recovery speed.


### 2.2. The Image Transmission in IoMT

The IoT will enable connections of a wide variety of things, ranging from small sensors to the cloud of servers (data storage and analysis). With the explosive growth of IoT, the device numbers will increase by several orders of magnitude, resulting in the great requirement for mechanisms that can reduce the computation, power, and communication loads for the end-devices. For the IoMT devices, this requirement becomes especially stringent due to large amount of multimedia data to be processed and transmitted. In addition, for some application scenarios such as security video surveillance, both the device-end and the server-end need to be capable of handling real-time processing.CS is one promising mechanism that can help to meet all the above requirements. First, CS is naturally a compression technique, so it can reduce the amount of data to be transferred, thus requiring less transmission power and bandwidth. Second, CS compression involves only linear operations, so the encoding complexity is extremely low compared with traditional nonlinear compression. For large image compression, by using methods such as separate sensing, the encoding complexity could further be reduced. Third, since the operation of CS compression is intrinsically global linear projection, every projected result has almost the same amount of information about the whole image. Therefore, CS compression is naturally robust to transmission loss. In other words, it can be regarded as a digital-fountain like erasure-correction codes. In all, the development of new and more efficient CS techniques for IoT is important, and our aim in this paper is to design reconstruction methods with not only low complexity but also high performance for real-time IoMT applications.


### 2.3. Prior Aided CS Image Reconstruction

In the basic CS framework, only the signal sparsity or compressibility is explored by the reconstruction. In fact, physical signals including IoMT images often show more features other than sparsity or compressibility, such as the structural and statistical features, which can be used as priors to further improve the reconstruction quality or efficiency. For example, the signal wavelet coefficients can be naturally organized into a tree structure, and those significant coefficients that contribute to most of the image energy often cluster along a few branches of this wavelet tree. The condensing sort and select algorithm (CSSA) [16] solving the optimal tree approximation is one of the good methods that can be employed by CS image reconstruction. Another example is the Block-Matching and 3D (BM3D) collaborative filtering model utilized to exploit the nonlocal self-similarity among pixels in global positioned patches [17].Besides structural prior exploration, statistical prior could also be used. In [18], recovery method adopting statistical information as prior knowledge was proposed to achieve a higher quality for the 1D sparse signals. Reference [19] proposed an innovative recovery method that achieves an estimation of the original signal by a series of iterations. In each iteration, a weight matrix which contains statistical a priori information of the original signal is used. However, the method involved numerous computation for an excess of iterations. A self-adapting statistical model of image expressed in wavelet domain is proposed by [20] to reduce noise in image storage. This method was not combined with compressed image sensing but still provides some inspiration for the next step research on CS. In [21], the Gaussian scale mixtures (GSM) model, which has been shown to effectively capture the local dependencies of intrascale wavelet coefficients [22], was incorporated into the proposed BLS-GSM recovery method using the Bayes least-square estimation. The BLS-GSM recovery has shown good performance, but the computation complexity is still too high to be affordable for IoMT devices. In summary, the development of low-complexity high-performance CS image reconstruction methods for IoMT image transmission systems remains an open problem requiring further investigation, which is just the aim of this paper.


## 2.1. The CS Basics

CS is applied for signal which is sparse or compressible. Under certain basesΨ=[ψ1,ψ2,…,ψn], a 1D signal x can be expressed as x=Ψs. We say x is c-sparse if it satisfies(1)x0=cor(2)s0=c.A signal is said to be nearly sparse or compressible if the largest c coefficients contribute to most of the signal energy, and the n-c others are small enough to be ignored. Signals in practical systems, including the images, are mostly compressible in a certain domain, so the CS theory can be used widely in natural signal processing.In CS framework, a spare or a nearly sparse signalx (size: n×1) is measured by a random matrix Φ (size: m×n, where m≪n) to generate a compressed result y (size: m×1):(3)y=Φx=ΦΨs=As,where A=ΦΨ.The reconstruction ofx corresponds to the problem of solving n unknown parameters from m linear equations, that is, a calculation for an underdetermined system of equations which contains innumerable solutions. Fortunately, s, the representation of the original signal under certain basis, is sparse enough, and the solution of the underdetermined system of equations can be achieved by solving the combinatorial optimization problem below:(4)s^=arg⁡min⁡s0,s.t.y=As.If the measurement matrixΦ is incoherent with the representation basis Ψ, especially when Φ is a random matrix and m=o(clog⁡(n/c)), A satisfies the restricted isometry property (RIP) with enormous probability [14]. The l0-norm optimization problem in (4) can be replaced by the l1-norm optimization problem which has the same unique solution; that is,(5)s^=arg⁡min⁡s1,s.t.y=As.The original signalx can then be obtained as(6)x^=Ψs.In conclusion, the crux of CS theory is to reconstruct the sparse representation of the original signal by making good use of the sparsity nature. There are mainly two kinds of reconstruction algorithms.


### 2.1.1. Greedy Algorithms

They are a kind of recovery methods that can find out the locally optimal solution but not the globally optimal solution. Orthogonal matching pursuit (OMP) is a classical greedy method for CS recovery, and it is the base of many later proposed advanced methods, such as CoSaMP [15] and A∗OMP. The crux of greedy algorithms is to find out columns of the matrix A that make significant contributions to the compressed result iteratively.


### 2.1.2. Convex Optimization Recovery

Different from the greedy algorithms, the convex optimization recovery methods aim to approach the globally optimal solution. The results by convex optimization recovery are usually more accurate, but the expense is more computation resources needed for the reconstruction. Convex optimization recovery can be used in situations where the recovery quality is preferred more than the recovery speed.


## 2.1.1. Greedy Algorithms

They are a kind of recovery methods that can find out the locally optimal solution but not the globally optimal solution. Orthogonal matching pursuit (OMP) is a classical greedy method for CS recovery, and it is the base of many later proposed advanced methods, such as CoSaMP [15] and A∗OMP. The crux of greedy algorithms is to find out columns of the matrix A that make significant contributions to the compressed result iteratively.


## 2.1.2. Convex Optimization Recovery

Different from the greedy algorithms, the convex optimization recovery methods aim to approach the globally optimal solution. The results by convex optimization recovery are usually more accurate, but the expense is more computation resources needed for the reconstruction. Convex optimization recovery can be used in situations where the recovery quality is preferred more than the recovery speed.


## 2.2. The Image Transmission in IoMT

The IoT will enable connections of a wide variety of things, ranging from small sensors to the cloud of servers (data storage and analysis). With the explosive growth of IoT, the device numbers will increase by several orders of magnitude, resulting in the great requirement for mechanisms that can reduce the computation, power, and communication loads for the end-devices. For the IoMT devices, this requirement becomes especially stringent due to large amount of multimedia data to be processed and transmitted. In addition, for some application scenarios such as security video surveillance, both the device-end and the server-end need to be capable of handling real-time processing.CS is one promising mechanism that can help to meet all the above requirements. First, CS is naturally a compression technique, so it can reduce the amount of data to be transferred, thus requiring less transmission power and bandwidth. Second, CS compression involves only linear operations, so the encoding complexity is extremely low compared with traditional nonlinear compression. For large image compression, by using methods such as separate sensing, the encoding complexity could further be reduced. Third, since the operation of CS compression is intrinsically global linear projection, every projected result has almost the same amount of information about the whole image. Therefore, CS compression is naturally robust to transmission loss. In other words, it can be regarded as a digital-fountain like erasure-correction codes. In all, the development of new and more efficient CS techniques for IoT is important, and our aim in this paper is to design reconstruction methods with not only low complexity but also high performance for real-time IoMT applications.


## 2.3. Prior Aided CS Image Reconstruction

In the basic CS framework, only the signal sparsity or compressibility is explored by the reconstruction. In fact, physical signals including IoMT images often show more features other than sparsity or compressibility, such as the structural and statistical features, which can be used as priors to further improve the reconstruction quality or efficiency. For example, the signal wavelet coefficients can be naturally organized into a tree structure, and those significant coefficients that contribute to most of the image energy often cluster along a few branches of this wavelet tree. The condensing sort and select algorithm (CSSA) [16] solving the optimal tree approximation is one of the good methods that can be employed by CS image reconstruction. Another example is the Block-Matching and 3D (BM3D) collaborative filtering model utilized to exploit the nonlocal self-similarity among pixels in global positioned patches [17].Besides structural prior exploration, statistical prior could also be used. In [18], recovery method adopting statistical information as prior knowledge was proposed to achieve a higher quality for the 1D sparse signals. Reference [19] proposed an innovative recovery method that achieves an estimation of the original signal by a series of iterations. In each iteration, a weight matrix which contains statistical a priori information of the original signal is used. However, the method involved numerous computation for an excess of iterations. A self-adapting statistical model of image expressed in wavelet domain is proposed by [20] to reduce noise in image storage. This method was not combined with compressed image sensing but still provides some inspiration for the next step research on CS. In [21], the Gaussian scale mixtures (GSM) model, which has been shown to effectively capture the local dependencies of intrascale wavelet coefficients [22], was incorporated into the proposed BLS-GSM recovery method using the Bayes least-square estimation. The BLS-GSM recovery has shown good performance, but the computation complexity is still too high to be affordable for IoMT devices. In summary, the development of low-complexity high-performance CS image reconstruction methods for IoMT image transmission systems remains an open problem requiring further investigation, which is just the aim of this paper.


## 3. Separate Image Sensing and the Reconstruction Methods

In this section, we will illustrate how the 2D images in IoMT applications can be compressed through the low-complexity separate sensing method and describe the existing corresponding reconstruction methods.


### 3.1. Separate Sensing for 2D Image Compression

To reduce the computation complexity of compressed image sensing, the 2D separate compression has been proved to be an effective method [10], especially for large-sized images. Let the original image signal be X, whose size is n×n, and let the random compression matrix be Φ, whose size is m×n (m<n). The original signal is measured column-wise by Φ and row-wise by ΦT, respectively, where ΦT is the transposition of Φ. So the compressed signal Y can be expressed as Y=ΦXΦT, as is shown in Figure 1. The size of compressed signal Y has been changed to m×m from n×n, after the compression.Figure 1
Illustration of the separate sensing method for image compression.Since the invention of separate compression, several corresponding recovery methods have been proposed in recent years, which will be described briefly in the next subsection.


### 3.2. Reconstruction Methods for Separate Sensing

#### 3.2.1. The Kronecker CS Method

Many existing methods of CS recovery are for 1D signals, not being able to be applied for 2D signals directly. A method of transforming the 2D image processing to 1D signal processing has been used for years by using the Kronecker product [11].The original image signalX at the receiver end usually has the equation form as follows:(7)Y=CSD,where S is the representation of X in the transform domain.LetY={y1,y2,…,yn}, X={x1,x2,…,xn}; then,(8)vec⁡Y=y1y2⋮yn,vec⁡X=x1x2⋮xn.Equation (7) can be transformed to(9)vec⁡Y=DT⊗Cvec⁡S,where DT⊗C denotes the Kronecker product of DT and C.The 2D signal recovery has been transformed to a 1D signal recovery by (9), so the method for 1D signals can be used on 2D signals directly. However, a main disadvantage with Kronecker CS is that the complexity has been increased sharply by the vector expansion operation. Therefore, the method involving Kronecker product cannot be used widely especially for real-time transmissions.


#### 3.2.2. The 2D-OMP Method

According to the separate compression, the signalsY=ΦXΦT obtained at the receiver can also be expressed as Y=ΦΨSΨTΦT=ASAT, where A=ΦΨ and S is the sparse representation of X.The crux of 2D-OMP [12] is to find a set of columns from A to represent Y as(10)Y=∑i=1n∑j=1nsi,jBi,j+R,where A=(a1,…,an), Bi,j=ai⊗ajT, and R denote the representation residual, which will be initialized as being equal to Y. By iterative procedures similar to the traditional 1D OMP, those columns of A that contribute significantly to R can be picked out sequentially, and the corresponding coefficients si,j will be estimated. The original signal X can finally be estimated after the sparse representation S of the original image is obtained. Compared with the Kronecker product based method, 2D-OMP is computationally more efficient. However, the computation complexity is still too high especially for large images.


#### 3.2.3. Separate-Combine Recovery

Separate recovery is a recovery method featured by its low computational complexity [13]. The compressed signal at the receiver end is(11)Y=ΦXΦT=ΦΨSΨTΦT,where X=ΨSΨT and S is the expression of X in wavelet domain. Let A=ΦΨ; we have(12)Y=ΦXΦT=ASAT.The recovery is composed of two steps: column-wise recovery and row-wise recovery. Let the intermediateS1=SAT; we have(13)Y=ΦXΦT=ASAT=AASTT=AS1.As is shown in Figure2, for the column-wise first recovery, the first step of the recovery is to obtain the intermediate solution S1 from the equation Y=AS1. Then, the second step of the recovery is to estimate the final solution S from the equation S1T=AST. The row-wise first recovery is just like the column-wise recovery. By combining the results of row-wise first recovery and column-wise first recovery together, a more accurate final result can be obtained.Figure 2
Illustration of separate recovery.Both of the two steps are to solve an underdetermined system of equations. The OMP algorithm can be chosen as the basic algorithm, by treatingS and S1 as sparse matrixes. The original signal X can be reconstituted on getting its sparse representation coefficients S.


## 3.1. Separate Sensing for 2D Image Compression

To reduce the computation complexity of compressed image sensing, the 2D separate compression has been proved to be an effective method [10], especially for large-sized images. Let the original image signal be X, whose size is n×n, and let the random compression matrix be Φ, whose size is m×n (m<n). The original signal is measured column-wise by Φ and row-wise by ΦT, respectively, where ΦT is the transposition of Φ. So the compressed signal Y can be expressed as Y=ΦXΦT, as is shown in Figure 1. The size of compressed signal Y has been changed to m×m from n×n, after the compression.Figure 1
Illustration of the separate sensing method for image compression.Since the invention of separate compression, several corresponding recovery methods have been proposed in recent years, which will be described briefly in the next subsection.


## 3.2. Reconstruction Methods for Separate Sensing

### 3.2.1. The Kronecker CS Method

Many existing methods of CS recovery are for 1D signals, not being able to be applied for 2D signals directly. A method of transforming the 2D image processing to 1D signal processing has been used for years by using the Kronecker product [11].The original image signalX at the receiver end usually has the equation form as follows:(7)Y=CSD,where S is the representation of X in the transform domain.LetY={y1,y2,…,yn}, X={x1,x2,…,xn}; then,(8)vec⁡Y=y1y2⋮yn,vec⁡X=x1x2⋮xn.Equation (7) can be transformed to(9)vec⁡Y=DT⊗Cvec⁡S,where DT⊗C denotes the Kronecker product of DT and C.The 2D signal recovery has been transformed to a 1D signal recovery by (9), so the method for 1D signals can be used on 2D signals directly. However, a main disadvantage with Kronecker CS is that the complexity has been increased sharply by the vector expansion operation. Therefore, the method involving Kronecker product cannot be used widely especially for real-time transmissions.


### 3.2.2. The 2D-OMP Method

According to the separate compression, the signalsY=ΦXΦT obtained at the receiver can also be expressed as Y=ΦΨSΨTΦT=ASAT, where A=ΦΨ and S is the sparse representation of X.The crux of 2D-OMP [12] is to find a set of columns from A to represent Y as(10)Y=∑i=1n∑j=1nsi,jBi,j+R,where A=(a1,…,an), Bi,j=ai⊗ajT, and R denote the representation residual, which will be initialized as being equal to Y. By iterative procedures similar to the traditional 1D OMP, those columns of A that contribute significantly to R can be picked out sequentially, and the corresponding coefficients si,j will be estimated. The original signal X can finally be estimated after the sparse representation S of the original image is obtained. Compared with the Kronecker product based method, 2D-OMP is computationally more efficient. However, the computation complexity is still too high especially for large images.


### 3.2.3. Separate-Combine Recovery

Separate recovery is a recovery method featured by its low computational complexity [13]. The compressed signal at the receiver end is(11)Y=ΦXΦT=ΦΨSΨTΦT,where X=ΨSΨT and S is the expression of X in wavelet domain. Let A=ΦΨ; we have(12)Y=ΦXΦT=ASAT.The recovery is composed of two steps: column-wise recovery and row-wise recovery. Let the intermediateS1=SAT; we have(13)Y=ΦXΦT=ASAT=AASTT=AS1.As is shown in Figure2, for the column-wise first recovery, the first step of the recovery is to obtain the intermediate solution S1 from the equation Y=AS1. Then, the second step of the recovery is to estimate the final solution S from the equation S1T=AST. The row-wise first recovery is just like the column-wise recovery. By combining the results of row-wise first recovery and column-wise first recovery together, a more accurate final result can be obtained.Figure 2
Illustration of separate recovery.Both of the two steps are to solve an underdetermined system of equations. The OMP algorithm can be chosen as the basic algorithm, by treatingS and S1 as sparse matrixes. The original signal X can be reconstituted on getting its sparse representation coefficients S.


## 3.2.1. The Kronecker CS Method

Many existing methods of CS recovery are for 1D signals, not being able to be applied for 2D signals directly. A method of transforming the 2D image processing to 1D signal processing has been used for years by using the Kronecker product [11].The original image signalX at the receiver end usually has the equation form as follows:(7)Y=CSD,where S is the representation of X in the transform domain.LetY={y1,y2,…,yn}, X={x1,x2,…,xn}; then,(8)vec⁡Y=y1y2⋮yn,vec⁡X=x1x2⋮xn.Equation (7) can be transformed to(9)vec⁡Y=DT⊗Cvec⁡S,where DT⊗C denotes the Kronecker product of DT and C.The 2D signal recovery has been transformed to a 1D signal recovery by (9), so the method for 1D signals can be used on 2D signals directly. However, a main disadvantage with Kronecker CS is that the complexity has been increased sharply by the vector expansion operation. Therefore, the method involving Kronecker product cannot be used widely especially for real-time transmissions.


## 3.2.2. The 2D-OMP Method

According to the separate compression, the signalsY=ΦXΦT obtained at the receiver can also be expressed as Y=ΦΨSΨTΦT=ASAT, where A=ΦΨ and S is the sparse representation of X.The crux of 2D-OMP [12] is to find a set of columns from A to represent Y as(10)Y=∑i=1n∑j=1nsi,jBi,j+R,where A=(a1,…,an), Bi,j=ai⊗ajT, and R denote the representation residual, which will be initialized as being equal to Y. By iterative procedures similar to the traditional 1D OMP, those columns of A that contribute significantly to R can be picked out sequentially, and the corresponding coefficients si,j will be estimated. The original signal X can finally be estimated after the sparse representation S of the original image is obtained. Compared with the Kronecker product based method, 2D-OMP is computationally more efficient. However, the computation complexity is still too high especially for large images.


## 3.2.3. Separate-Combine Recovery

Separate recovery is a recovery method featured by its low computational complexity [13]. The compressed signal at the receiver end is(11)Y=ΦXΦT=ΦΨSΨTΦT,where X=ΨSΨT and S is the expression of X in wavelet domain. Let A=ΦΨ; we have(12)Y=ΦXΦT=ASAT.The recovery is composed of two steps: column-wise recovery and row-wise recovery. Let the intermediateS1=SAT; we have(13)Y=ΦXΦT=ASAT=AASTT=AS1.As is shown in Figure2, for the column-wise first recovery, the first step of the recovery is to obtain the intermediate solution S1 from the equation Y=AS1. Then, the second step of the recovery is to estimate the final solution S from the equation S1T=AST. The row-wise first recovery is just like the column-wise recovery. By combining the results of row-wise first recovery and column-wise first recovery together, a more accurate final result can be obtained.Figure 2
Illustration of separate recovery.Both of the two steps are to solve an underdetermined system of equations. The OMP algorithm can be chosen as the basic algorithm, by treatingS and S1 as sparse matrixes. The original signal X can be reconstituted on getting its sparse representation coefficients S.


## 4. Statistical Prior Aided Low-Complexity Recovery

This paper is aiming to propose a low-complexity yet high-performance recovery method for compressed image sensing. For the encoder-end, the separate sensing method described in Section3 is adopted due to its low complexity. For the decoder-end, instead of using greedy iterations such as OMP algorithm in the two-stage separate reconstruction, a brand new algorithm is proposed to achieve performance improvement and complexity reduction. First, we investigate the statistical characteristics of natural images, establishing a statistical model. Then the statistical information is used to generate a weight matrix K to constrain the energy distribution of the estimates. The recovery strategy is named as Direct Recovery strategy with Fixed weight Matrix (DR-FM). Furthermore, an improved version which uses refined weigh matrix, that is, the Iterative Recovery strategy with Refined weight Matrix (IR-RM), is proposed to obtain a more accurate result.


### 4.1. Image Statistical Information

The image representation coefficient matrix in wavelet domain has significant characteristics, where low frequency components concentrate at top left corner and high frequency components concentrate at bottom right side. The matrixS is partitioned to blocks of different sizes, as Figure 3 shows. Sum the elements of S in each row; that is,(14)Ssum=∑js1,j∑js2,j⋮∑jsn,j;then ladder-like decaying distribution of Ssum can be observed. Take the image “Lena” for instance; as Figure 4 shows, the blue curve shows the average values of each row of the image in wavelet representation. The result shows that the average values conform to ladder-like decaying. The matrix Ssum has multiple energy levels, and the minimum length of one energy level is 8, then 8, 16, and 32, and so on. Table 1 shows the length of each energy level and the averaged absolute values of pixels contained in each energy level.Table 1
The length of each energy level and the average absolute values of pixels contained in each energy level.
Level
Length of level
The average absolute values of pixels
1
8
180.4892
2
8
42.2611
3
16
27.1487
4
32
17.0842
5
64
9.800
6
128
5.0789
7
256
2.5370
⋯
⋯
⋯Figure 3
Illustration of the wavelet coefficient blocks.Figure 4
The ladder-like decaying attributes of wavelet representations for images.As the energy level index goes up, the absolute values shows substantially decaying. And the exponential decaying can be used to fit it; as Figure5 shows, the result of fitting is(15)fx=aebx,where a=597.5 and b=-1.208. From experiments on wavelets representation of other images, similar decaying properties can be obtained. Therefore, the ladder-like decaying property of images in wavelet domain can be well characterized by the following model:(16)fn=aebn,where n is the energy level index of a certain point and a and b are two model parameters. From multiple investigations on different images, it is found that the values of parameter b mostly stay between −0.7 and −1.3. This decaying model can be used as the a priori statistical knowledge in the recovery strategy design.Figure 5
Curve fitting for the decaying of image wavelet representations.


### 4.2. Algorithm Design and the Analysis

To satisfy different application demands in practice, two recovery strategies are proposed in this paper, that is, the DR-FM strategy and the IR-RM strategy. In each strategy, the recovery process is composed of two steps, that is, the recovery for row-wise (or column-wise) intermediates and the recovery for column-wise (or row-wise) final results. Here, we take the column-wise first recovery followed by row-wise recovery as the example for description and analysis.


#### 4.2.1. The Direct Recovery Strategy with Fixed Weight Matrix (DR-FM)

Based on separate recovery, reconstructingS1 from Y1=A×S1 and reconstructing S2 from Y2=A×S2 are the two steps of the algorithm, as is shown in Figure 2. Each of the two matrix equations is an underdetermined equation system. For the statistical model established, the crux of each step is to identify the best match set based on the decaying law from candidate solutions of the underdetermined equation system. The two steps of the recovery can be expressed as follows:(17)Y=ASAT=AASTT=AS1Y1=AST,where S1=(AST)T and Y1=S1T.The two steps have a similar equation form. Take the first step for example, that is, reconstructS1 from Y=A×S1, and select a column of it, that is, y=A×z, where y (size: m×1) is a column vector in matrix Y and z (size: n×1) is a column vector in S1; let(18)z=z1z2⋮zn,A=a11…a1n⋮⋱⋮am1⋯amn.Based on the decaying model of the energy level, that is,f(n)=aebn, the weight matrix K can be constructed as(19)K=k10⋯00k2⋯0⋮⋮⋱⋮00⋯kn,where the first 8 values are in the first energy level, being calculated as k1 : 8=a×eb, and the 9 to 16 values are the second energy lever, being calculated as k9 : 16=a×e2×b. In general,(20)ki=a×exp⁡bifi=1,2,…,8a×exp⁡blog⁡i-1else,where ⌈∗⌉ denotes the round-up operation. The object of constraining the estimates to conform to the previous built decaying model can be easily formulated as the weighted least-square problem as follows:(21)z^=arg⁡min⁡z12k12+z22k22+⋯+zn2kn2.s.t.a11⋯a1n⋮⋱⋮am1⋯amn×z1⋮zn=yFor simplicity, we can normalize the elements of the weight matrixK by a×exp⁡(b); that is,(22)ki=1ifi=1,2,…,8exp⁡blog⁡i-2else.The result of (21) can be easily obtained as(23)z^=K2×AT×A×K2×AT-1×y.Since this method involves simply matrix multiplication, we can directly obtain the solution ofY=A×S1 accordingly as(24)S^1=K2×AT×A×K2×AT-1×Y.Similarly, this method can be applied to the second-step reconstruction fromY1=AST=AS2. The recovery result can be expressed as(25)S^2=K2×AT×A×K2×AT-1×Y1.Finally, the sparse representation of the image can be obtained asS^=S2T, and the original image can be reconstructed as(26)X^=ΨS^ΨT.The whole process of this method is based on matrix multiplication. Therefore, this method has extremely low complexity and can improve the operation speed effectively.


#### 4.2.2. The Iterative Recovery Strategy with Refined Weight Matrix (IR-RM)

By observing the differences between the recovery result ofS1 in the first step and the recovery result of S2 in the second-step, it can be seen that S1 fits the decaying model better than S2 since the latter is sparser. To get more accurate recovery results, the strategy with an extra iteration can be introduced during the second-step. Due to the extra iteration based on the first-step calculation, the result would be closer to the original signal.The energy distribution weight matrixK in the second-step recovery of DR-FM is derived by giving an original estimation based on the energy distribution of general signal. To obtain a result that is closer to the real signal based on the original estimation, the weight matrix K has to be refined. By comparison with the weight matrix K, it can be seen that the one-time estimation z^ is closer to the energy distribution of real signal, while there is still slight difference, as is shown in Figure 6.Figure 6
Comparison between the estimation results and the weight matrixK.
(a)
Comparison between results of the first iteration and real values
(b)
Comparison between results of the first iteration and the weight matrixKSimilar to DR-FM, the two steps of IR-RM can be expressed as(27)Y=ASAT=AASTT=AS1Y1=AST,where S1=(AST)T and Y1=S1T.The first step of IR-RM is the same as DR-FM, and the recovery can be completed in a linear step(28)S^1=K2×AT×A×K2×AT-1×Y.In the second-step recovery, letY1={y1(2),y2(2),…,yn(2)}, ST={z1(2),z2(2),…,zn(2)}. Then we have yj(2)=Azj(2), where j∈1,2,…,m.Let(29)zj2=z1j2z2j2⋮znj2,yj2=y1j2y2j2⋮ymj2,A=a11…a1n⋮⋱⋮am1⋯amn.The equation yj(2)=Azj(2) can be expressed as(30)y1j2y2j2⋮ymj2=a11…a1n⋮⋱⋮am1⋯amnz1j2z2j2⋮znj2.Similar to the DR-FM, the solution of the underdetermined equation system is(31)z^j2=K2×AT×A×K2×AT-1yj2.To obtain better recovery results, the weight matrix can be refined column by column, respectively, in the extra iteration of IR-RM as(32)K2,j=k1,j20⋯00k2,j2⋯0⋮⋮⋱⋮00⋯kn,j2.The elements on the main diagonal satisfyki(2)∝z^i. z^i is the element of the results from the method with one-time estimation z^; that is, z^=z^1,z^2,…,z^nT. Hence,(33)K2,j=ez^1j20⋯00ez^2j2⋯0⋮⋮⋱⋮00⋯ez^nj2.The parameter e is an arbitrary constant, which can be divided out:(34)K2,j=z^1j20⋯00z^2j2⋯0⋮⋮⋱⋮00⋯z^nj2.The reconstruction using the refined weight matrix can be formulated as(35)z^j3=arg⁡min⁡z1j32z1j22+z2j32z2j22+⋯+znj32znj22s.t.a11…a1n⋮⋱⋮am1⋯amn×z1j3z2j3⋮znj3=y1j2y2j2⋮ymj2,where z^j(3)=z^1j(3),z^2j(3),…,z^nj(3)T is the recovery result of the extra iteration of IR-RM.The solution of the above problem is(36)z^j3=K2,j2×AT×A×K2,j2×AT-1yj2.Since the results of the method with IR-RM derived from the second step are different from DR-FM’s, the solution with an extra iteration cannot be converted to the form of multiplication of matrixes and can only be solved vector by vector.


## 4.1. Image Statistical Information

The image representation coefficient matrix in wavelet domain has significant characteristics, where low frequency components concentrate at top left corner and high frequency components concentrate at bottom right side. The matrixS is partitioned to blocks of different sizes, as Figure 3 shows. Sum the elements of S in each row; that is,(14)Ssum=∑js1,j∑js2,j⋮∑jsn,j;then ladder-like decaying distribution of Ssum can be observed. Take the image “Lena” for instance; as Figure 4 shows, the blue curve shows the average values of each row of the image in wavelet representation. The result shows that the average values conform to ladder-like decaying. The matrix Ssum has multiple energy levels, and the minimum length of one energy level is 8, then 8, 16, and 32, and so on. Table 1 shows the length of each energy level and the averaged absolute values of pixels contained in each energy level.Table 1
The length of each energy level and the average absolute values of pixels contained in each energy level.
Level
Length of level
The average absolute values of pixels
1
8
180.4892
2
8
42.2611
3
16
27.1487
4
32
17.0842
5
64
9.800
6
128
5.0789
7
256
2.5370
⋯
⋯
⋯Figure 3
Illustration of the wavelet coefficient blocks.Figure 4
The ladder-like decaying attributes of wavelet representations for images.As the energy level index goes up, the absolute values shows substantially decaying. And the exponential decaying can be used to fit it; as Figure5 shows, the result of fitting is(15)fx=aebx,where a=597.5 and b=-1.208. From experiments on wavelets representation of other images, similar decaying properties can be obtained. Therefore, the ladder-like decaying property of images in wavelet domain can be well characterized by the following model:(16)fn=aebn,where n is the energy level index of a certain point and a and b are two model parameters. From multiple investigations on different images, it is found that the values of parameter b mostly stay between −0.7 and −1.3. This decaying model can be used as the a priori statistical knowledge in the recovery strategy design.Figure 5
Curve fitting for the decaying of image wavelet representations.


## 4.2. Algorithm Design and the Analysis

To satisfy different application demands in practice, two recovery strategies are proposed in this paper, that is, the DR-FM strategy and the IR-RM strategy. In each strategy, the recovery process is composed of two steps, that is, the recovery for row-wise (or column-wise) intermediates and the recovery for column-wise (or row-wise) final results. Here, we take the column-wise first recovery followed by row-wise recovery as the example for description and analysis.


### 4.2.1. The Direct Recovery Strategy with Fixed Weight Matrix (DR-FM)

Based on separate recovery, reconstructingS1 from Y1=A×S1 and reconstructing S2 from Y2=A×S2 are the two steps of the algorithm, as is shown in Figure 2. Each of the two matrix equations is an underdetermined equation system. For the statistical model established, the crux of each step is to identify the best match set based on the decaying law from candidate solutions of the underdetermined equation system. The two steps of the recovery can be expressed as follows:(17)Y=ASAT=AASTT=AS1Y1=AST,where S1=(AST)T and Y1=S1T.The two steps have a similar equation form. Take the first step for example, that is, reconstructS1 from Y=A×S1, and select a column of it, that is, y=A×z, where y (size: m×1) is a column vector in matrix Y and z (size: n×1) is a column vector in S1; let(18)z=z1z2⋮zn,A=a11…a1n⋮⋱⋮am1⋯amn.Based on the decaying model of the energy level, that is,f(n)=aebn, the weight matrix K can be constructed as(19)K=k10⋯00k2⋯0⋮⋮⋱⋮00⋯kn,where the first 8 values are in the first energy level, being calculated as k1 : 8=a×eb, and the 9 to 16 values are the second energy lever, being calculated as k9 : 16=a×e2×b. In general,(20)ki=a×exp⁡bifi=1,2,…,8a×exp⁡blog⁡i-1else,where ⌈∗⌉ denotes the round-up operation. The object of constraining the estimates to conform to the previous built decaying model can be easily formulated as the weighted least-square problem as follows:(21)z^=arg⁡min⁡z12k12+z22k22+⋯+zn2kn2.s.t.a11⋯a1n⋮⋱⋮am1⋯amn×z1⋮zn=yFor simplicity, we can normalize the elements of the weight matrixK by a×exp⁡(b); that is,(22)ki=1ifi=1,2,…,8exp⁡blog⁡i-2else.The result of (21) can be easily obtained as(23)z^=K2×AT×A×K2×AT-1×y.Since this method involves simply matrix multiplication, we can directly obtain the solution ofY=A×S1 accordingly as(24)S^1=K2×AT×A×K2×AT-1×Y.Similarly, this method can be applied to the second-step reconstruction fromY1=AST=AS2. The recovery result can be expressed as(25)S^2=K2×AT×A×K2×AT-1×Y1.Finally, the sparse representation of the image can be obtained asS^=S2T, and the original image can be reconstructed as(26)X^=ΨS^ΨT.The whole process of this method is based on matrix multiplication. Therefore, this method has extremely low complexity and can improve the operation speed effectively.


### 4.2.2. The Iterative Recovery Strategy with Refined Weight Matrix (IR-RM)

By observing the differences between the recovery result ofS1 in the first step and the recovery result of S2 in the second-step, it can be seen that S1 fits the decaying model better than S2 since the latter is sparser. To get more accurate recovery results, the strategy with an extra iteration can be introduced during the second-step. Due to the extra iteration based on the first-step calculation, the result would be closer to the original signal.The energy distribution weight matrixK in the second-step recovery of DR-FM is derived by giving an original estimation based on the energy distribution of general signal. To obtain a result that is closer to the real signal based on the original estimation, the weight matrix K has to be refined. By comparison with the weight matrix K, it can be seen that the one-time estimation z^ is closer to the energy distribution of real signal, while there is still slight difference, as is shown in Figure 6.Figure 6
Comparison between the estimation results and the weight matrixK.
(a)
Comparison between results of the first iteration and real values
(b)
Comparison between results of the first iteration and the weight matrixKSimilar to DR-FM, the two steps of IR-RM can be expressed as(27)Y=ASAT=AASTT=AS1Y1=AST,where S1=(AST)T and Y1=S1T.The first step of IR-RM is the same as DR-FM, and the recovery can be completed in a linear step(28)S^1=K2×AT×A×K2×AT-1×Y.In the second-step recovery, letY1={y1(2),y2(2),…,yn(2)}, ST={z1(2),z2(2),…,zn(2)}. Then we have yj(2)=Azj(2), where j∈1,2,…,m.Let(29)zj2=z1j2z2j2⋮znj2,yj2=y1j2y2j2⋮ymj2,A=a11…a1n⋮⋱⋮am1⋯amn.The equation yj(2)=Azj(2) can be expressed as(30)y1j2y2j2⋮ymj2=a11…a1n⋮⋱⋮am1⋯amnz1j2z2j2⋮znj2.Similar to the DR-FM, the solution of the underdetermined equation system is(31)z^j2=K2×AT×A×K2×AT-1yj2.To obtain better recovery results, the weight matrix can be refined column by column, respectively, in the extra iteration of IR-RM as(32)K2,j=k1,j20⋯00k2,j2⋯0⋮⋮⋱⋮00⋯kn,j2.The elements on the main diagonal satisfyki(2)∝z^i. z^i is the element of the results from the method with one-time estimation z^; that is, z^=z^1,z^2,…,z^nT. Hence,(33)K2,j=ez^1j20⋯00ez^2j2⋯0⋮⋮⋱⋮00⋯ez^nj2.The parameter e is an arbitrary constant, which can be divided out:(34)K2,j=z^1j20⋯00z^2j2⋯0⋮⋮⋱⋮00⋯z^nj2.The reconstruction using the refined weight matrix can be formulated as(35)z^j3=arg⁡min⁡z1j32z1j22+z2j32z2j22+⋯+znj32znj22s.t.a11…a1n⋮⋱⋮am1⋯amn×z1j3z2j3⋮znj3=y1j2y2j2⋮ymj2,where z^j(3)=z^1j(3),z^2j(3),…,z^nj(3)T is the recovery result of the extra iteration of IR-RM.The solution of the above problem is(36)z^j3=K2,j2×AT×A×K2,j2×AT-1yj2.Since the results of the method with IR-RM derived from the second step are different from DR-FM’s, the solution with an extra iteration cannot be converted to the form of multiplication of matrixes and can only be solved vector by vector.


## 4.2.1. The Direct Recovery Strategy with Fixed Weight Matrix (DR-FM)

Based on separate recovery, reconstructingS1 from Y1=A×S1 and reconstructing S2 from Y2=A×S2 are the two steps of the algorithm, as is shown in Figure 2. Each of the two matrix equations is an underdetermined equation system. For the statistical model established, the crux of each step is to identify the best match set based on the decaying law from candidate solutions of the underdetermined equation system. The two steps of the recovery can be expressed as follows:(17)Y=ASAT=AASTT=AS1Y1=AST,where S1=(AST)T and Y1=S1T.The two steps have a similar equation form. Take the first step for example, that is, reconstructS1 from Y=A×S1, and select a column of it, that is, y=A×z, where y (size: m×1) is a column vector in matrix Y and z (size: n×1) is a column vector in S1; let(18)z=z1z2⋮zn,A=a11…a1n⋮⋱⋮am1⋯amn.Based on the decaying model of the energy level, that is,f(n)=aebn, the weight matrix K can be constructed as(19)K=k10⋯00k2⋯0⋮⋮⋱⋮00⋯kn,where the first 8 values are in the first energy level, being calculated as k1 : 8=a×eb, and the 9 to 16 values are the second energy lever, being calculated as k9 : 16=a×e2×b. In general,(20)ki=a×exp⁡bifi=1,2,…,8a×exp⁡blog⁡i-1else,where ⌈∗⌉ denotes the round-up operation. The object of constraining the estimates to conform to the previous built decaying model can be easily formulated as the weighted least-square problem as follows:(21)z^=arg⁡min⁡z12k12+z22k22+⋯+zn2kn2.s.t.a11⋯a1n⋮⋱⋮am1⋯amn×z1⋮zn=yFor simplicity, we can normalize the elements of the weight matrixK by a×exp⁡(b); that is,(22)ki=1ifi=1,2,…,8exp⁡blog⁡i-2else.The result of (21) can be easily obtained as(23)z^=K2×AT×A×K2×AT-1×y.Since this method involves simply matrix multiplication, we can directly obtain the solution ofY=A×S1 accordingly as(24)S^1=K2×AT×A×K2×AT-1×Y.Similarly, this method can be applied to the second-step reconstruction fromY1=AST=AS2. The recovery result can be expressed as(25)S^2=K2×AT×A×K2×AT-1×Y1.Finally, the sparse representation of the image can be obtained asS^=S2T, and the original image can be reconstructed as(26)X^=ΨS^ΨT.The whole process of this method is based on matrix multiplication. Therefore, this method has extremely low complexity and can improve the operation speed effectively.


## 4.2.2. The Iterative Recovery Strategy with Refined Weight Matrix (IR-RM)

By observing the differences between the recovery result ofS1 in the first step and the recovery result of S2 in the second-step, it can be seen that S1 fits the decaying model better than S2 since the latter is sparser. To get more accurate recovery results, the strategy with an extra iteration can be introduced during the second-step. Due to the extra iteration based on the first-step calculation, the result would be closer to the original signal.The energy distribution weight matrixK in the second-step recovery of DR-FM is derived by giving an original estimation based on the energy distribution of general signal. To obtain a result that is closer to the real signal based on the original estimation, the weight matrix K has to be refined. By comparison with the weight matrix K, it can be seen that the one-time estimation z^ is closer to the energy distribution of real signal, while there is still slight difference, as is shown in Figure 6.Figure 6
Comparison between the estimation results and the weight matrixK.
(a)
Comparison between results of the first iteration and real values
(b)
Comparison between results of the first iteration and the weight matrixKSimilar to DR-FM, the two steps of IR-RM can be expressed as(27)Y=ASAT=AASTT=AS1Y1=AST,where S1=(AST)T and Y1=S1T.The first step of IR-RM is the same as DR-FM, and the recovery can be completed in a linear step(28)S^1=K2×AT×A×K2×AT-1×Y.In the second-step recovery, letY1={y1(2),y2(2),…,yn(2)}, ST={z1(2),z2(2),…,zn(2)}. Then we have yj(2)=Azj(2), where j∈1,2,…,m.Let(29)zj2=z1j2z2j2⋮znj2,yj2=y1j2y2j2⋮ymj2,A=a11…a1n⋮⋱⋮am1⋯amn.The equation yj(2)=Azj(2) can be expressed as(30)y1j2y2j2⋮ymj2=a11…a1n⋮⋱⋮am1⋯amnz1j2z2j2⋮znj2.Similar to the DR-FM, the solution of the underdetermined equation system is(31)z^j2=K2×AT×A×K2×AT-1yj2.To obtain better recovery results, the weight matrix can be refined column by column, respectively, in the extra iteration of IR-RM as(32)K2,j=k1,j20⋯00k2,j2⋯0⋮⋮⋱⋮00⋯kn,j2.The elements on the main diagonal satisfyki(2)∝z^i. z^i is the element of the results from the method with one-time estimation z^; that is, z^=z^1,z^2,…,z^nT. Hence,(33)K2,j=ez^1j20⋯00ez^2j2⋯0⋮⋮⋱⋮00⋯ez^nj2.The parameter e is an arbitrary constant, which can be divided out:(34)K2,j=z^1j20⋯00z^2j2⋯0⋮⋮⋱⋮00⋯z^nj2.The reconstruction using the refined weight matrix can be formulated as(35)z^j3=arg⁡min⁡z1j32z1j22+z2j32z2j22+⋯+znj32znj22s.t.a11…a1n⋮⋱⋮am1⋯amn×z1j3z2j3⋮znj3=y1j2y2j2⋮ymj2,where z^j(3)=z^1j(3),z^2j(3),…,z^nj(3)T is the recovery result of the extra iteration of IR-RM.The solution of the above problem is(36)z^j3=K2,j2×AT×A×K2,j2×AT-1yj2.Since the results of the method with IR-RM derived from the second step are different from DR-FM’s, the solution with an extra iteration cannot be converted to the form of multiplication of matrixes and can only be solved vector by vector.


## 5. Experimental Results for Multimedia Images

To demonstrate the efficacy of the proposed recovery strategy, simulation experiments have been performed with MATLAB in a 64-bit operating system with 8 GB memory and CPU clock speed of 3.30 GHz. Different sizes of Lena, Girl, and Man images are chosen in the experiments, respectively, and the results are evaluated by Peak Signal to Noise Ratio (PSNR) and reconstruction time.


### 5.1. Reconstruction Performance Evaluation

The original image and recovery images of Lena, Girl, and Man are shown in Figures7–9. The parameter values and recovery PSNR are shown in Table 2.Table 2
The parameters and recovery PSNR of Lena, Girl, and Man.
Image size
Size after compression
Parameterb
PSNR of DR-FM
PSNR of IR-RM
Lena
512
×
512
256
×
256
b
=
-
1.2
29.0912
29.4311
Girl
512
×
512
256
×
256
b
=
-
1.2
30.3810
30.5047
Man
512
×
512
256
×
256
b
=
-
1.2
23.8872
24.1672Figure 7
The simulation results of the image “Lena.”
(a)
The original image
(b)
Results by DR-FM
(c)
Results by IR-RMFigure 8
The simulation results of the image “Girl.”
(a)
The original image
(b)
Results by DR-FM
(c)
Results by IR-RMFigure 9
The simulation results of the image “Man.”
(a)
The original image
(b)
Results by DR-FM
(c)
Results by IR-RMBased on the results of the three images, DR-FM and IR-RM are both effective for image CS recovery. Compared with DR-FM, the IR-RM recovery quality has been slightly improved after the extra iteration. As is described in Section4.2, the quality improvement is due to the weight matrix refinement by the extra iteration, which aims to generate more accurate energy distribution constrains on the final results.


### 5.2. Comparison with Other Reconstruction Methods

Simulation on Lena images using methods in the literature, including the 2D-OMP [12] and the separate-combine recovery [13] on the same set of conditions, is also obtained and shown in Figures 10 and 11. The two figures show the comparison of reconstruction time and PSNR, respectively. The reconstruction time results are on logarithmic coordinate to show a clearer comparison. The image “Lena” is used for simulation, and the compression ratio is (m/n)2=(3/4)2=0.5625. Note that the results of using 2D-OMP to reconstruct images with dimensions above 512×512 are not obtained due to the extremely high complexity.Figure 10
Recovery time comparison between different methods.Figure 11
Recovery quality comparison between different methods.Compared with existing methods in the literature, that is, the separate-combine recovery and the 2D-OMP, the two proposed strategies with statistic information aids (DR-FM and IR-RM) are improved on both computational efficiency (less reconstruction time) and recovering accuracy (higher PSNR). The improvement of computational efficiency is mainly because the proposed strategies involve only linear operation of products between matrixes or vectors, while the existing methods are intrinsically greedy algorithms involving multiple nonlinear iterations. The improvement of reconstruction quality is due to the use of statistical information, in terms of the weigh matrix setting which aims to reflect the energy distribution of natural images in the wavelet domain. The algorithm with the highest accuracy is IR-RM, while the method with the lowest computing complexity is DR-FM. As is noted above, the performance difference between DR-FM and IR-RM is due to the weight matrix refinement by the extra iterations in the second step. In all, the two strategies proposed in this paper are more promising for IoMT applications, and they can be applied flexibly according to the specific requirements, namely, to achieve better recovering accuracy or to be more computation efficient.


### 5.3. Analysis of the Method Sensitiveness on Statistical Information

The method proposed in this paper relies on the statistics a lot. Although different images contain different statistical information, favorable result could be stably obtained as long as the fitting shows a fast decay. There is little influence on the final result, as the parameterb fluctuates, as shown in Figure 12. The simulated images are Lena, Girl, and Man. The compressing ratio of the simulation is (m/n)2=(1/2)2=0.25, and the size of the simulation images is 512×512.Figure 12
Recovery quality of Lena, Girl, and Man with different parameterb.Therefore, the proposed method is not sensitive to the parameterb of statistical information as long as the decaying is fast enough, which means that it can be used widely and stably.


## 5.1. Reconstruction Performance Evaluation

The original image and recovery images of Lena, Girl, and Man are shown in Figures7–9. The parameter values and recovery PSNR are shown in Table 2.Table 2
The parameters and recovery PSNR of Lena, Girl, and Man.
Image size
Size after compression
Parameterb
PSNR of DR-FM
PSNR of IR-RM
Lena
512
×
512
256
×
256
b
=
-
1.2
29.0912
29.4311
Girl
512
×
512
256
×
256
b
=
-
1.2
30.3810
30.5047
Man
512
×
512
256
×
256
b
=
-
1.2
23.8872
24.1672Figure 7
The simulation results of the image “Lena.”
(a)
The original image
(b)
Results by DR-FM
(c)
Results by IR-RMFigure 8
The simulation results of the image “Girl.”
(a)
The original image
(b)
Results by DR-FM
(c)
Results by IR-RMFigure 9
The simulation results of the image “Man.”
(a)
The original image
(b)
Results by DR-FM
(c)
Results by IR-RMBased on the results of the three images, DR-FM and IR-RM are both effective for image CS recovery. Compared with DR-FM, the IR-RM recovery quality has been slightly improved after the extra iteration. As is described in Section4.2, the quality improvement is due to the weight matrix refinement by the extra iteration, which aims to generate more accurate energy distribution constrains on the final results.


## 5.2. Comparison with Other Reconstruction Methods

Simulation on Lena images using methods in the literature, including the 2D-OMP [12] and the separate-combine recovery [13] on the same set of conditions, is also obtained and shown in Figures 10 and 11. The two figures show the comparison of reconstruction time and PSNR, respectively. The reconstruction time results are on logarithmic coordinate to show a clearer comparison. The image “Lena” is used for simulation, and the compression ratio is (m/n)2=(3/4)2=0.5625. Note that the results of using 2D-OMP to reconstruct images with dimensions above 512×512 are not obtained due to the extremely high complexity.Figure 10
Recovery time comparison between different methods.Figure 11
Recovery quality comparison between different methods.Compared with existing methods in the literature, that is, the separate-combine recovery and the 2D-OMP, the two proposed strategies with statistic information aids (DR-FM and IR-RM) are improved on both computational efficiency (less reconstruction time) and recovering accuracy (higher PSNR). The improvement of computational efficiency is mainly because the proposed strategies involve only linear operation of products between matrixes or vectors, while the existing methods are intrinsically greedy algorithms involving multiple nonlinear iterations. The improvement of reconstruction quality is due to the use of statistical information, in terms of the weigh matrix setting which aims to reflect the energy distribution of natural images in the wavelet domain. The algorithm with the highest accuracy is IR-RM, while the method with the lowest computing complexity is DR-FM. As is noted above, the performance difference between DR-FM and IR-RM is due to the weight matrix refinement by the extra iterations in the second step. In all, the two strategies proposed in this paper are more promising for IoMT applications, and they can be applied flexibly according to the specific requirements, namely, to achieve better recovering accuracy or to be more computation efficient.


## 5.3. Analysis of the Method Sensitiveness on Statistical Information

The method proposed in this paper relies on the statistics a lot. Although different images contain different statistical information, favorable result could be stably obtained as long as the fitting shows a fast decay. There is little influence on the final result, as the parameterb fluctuates, as shown in Figure 12. The simulated images are Lena, Girl, and Man. The compressing ratio of the simulation is (m/n)2=(1/2)2=0.25, and the size of the simulation images is 512×512.Figure 12
Recovery quality of Lena, Girl, and Man with different parameterb.Therefore, the proposed method is not sensitive to the parameterb of statistical information as long as the decaying is fast enough, which means that it can be used widely and stably.


## 6. Conclusions

To propose low-complexity CS image reconstruction methods for the green computation and communications of IoMT, this paper proposed a method with two reconstruction strategies based on the statistical prior information. Simulations show that it can reduce the complexity greatly by making good use of the prior information of image representations in the wavelet domain, even when the original image has been compressed heavily. The computation of the proposed DR-FM involves only linear matrix multiplications. The computation complexity of the proposed IR-RM is slightly heavier than DR-FM but still decrease a lot compared with traditional methods. The reconstruction quality of DR-FM is much better than traditional methods due to the use of statistical prior information, and the extra iteration of IR-RM further improves the reconstruction quality. The two strategies presented in this paper can be selected flexibly according to the practical requirements of different IoMT application situations, that is, higher reconstruction quality requirements or less reconstruction time requirements.



---

*Source: 2314062-2017-03-16.xml*