Image Reconstruction using Matched Wavelet Estimated from Data Sensed Compressively using Partial Canonical Identity Matrix

TL;DR

This paper introduces a joint framework for estimating matched wavelets directly from compressively sensed images to improve reconstruction quality and speed, utilizing a novel multi-level wavelet decomposition strategy.

Contribution

It presents a method to design matched wavelets from compressively sensed data, reducing sensing time and enhancing reconstruction performance over existing techniques.

Findings

Faster image reconstruction compared to standard wavelet methods.

Improved reconstruction quality with the proposed wavelet and decomposition strategy.

Reduced sensing and processing time without performance loss.

Abstract

This paper proposes a joint framework wherein lifting-based, separable, image-matched wavelets are estimated from compressively sensed (CS) images and used for the reconstruction of the same. Matched wavelet can be easily designed if full image is available. Also matched wavelet may provide better reconstruction results in CS application compared to standard wavelet sparsifying basis. Since in CS application, we have compressively sensed image instead of full image, existing methods of designing matched wavelet cannot be used. Thus, we propose a joint framework that estimates matched wavelet from the compressively sensed images and also reconstructs full images. This paper has three significant contributions. First, lifting-based, image-matched separable wavelet is designed from compressively sensed images and is also used to reconstruct the same. Second, a simple sensing matrix is employed to sample data at sub-Nyquist rate such that sensing and reconstruction time is reduced considerably without any noticeable degradation in the reconstruction performance. Third, a new multi-level L-Pyramid wavelet decomposition strategy is provided for separable wavelet implementation on images that leads to improved reconstruction performance. Compared to CS-based reconstruction using standard wavelets with Gaussian sensing matrix and with existing wavelet decomposition strategy, the proposed methodology provides faster and better image reconstruction in compressive sensing application.

Tables

Table 1. TABLE I : Coefficients of analysis filters designed with image ‘Lena’

Image matched 5/3 filters
	Predict Stage filter $T\left(z\right)$	Update Stage filter $S\left(z\right)$	Filter coefficients
Column Matched	[0.5028 0.4941]	[0.2858 0.2790]	$h_0\left[n\right]$=[-0.1412 0.2858 0.7185 0.279 -0.1403]
			$h_1\left[n\right]$=[-0.4941 1.0000 -0.5028]
Row Matched	[0.4959 0.5084]	[0.2775 0.2871]	$h_0\left[n\right]$=[-0.1411 0.2775 0.7164 0.2871 -0.1424]
			$h_1\left[n\right]$=[-0.5084 1.0000 -0.4959 ]

Table 2. TABLE II : Reconstruction accuracy on CS-based image reconstruction with standard wavelets and with image-matched wavelets designed from compressively sensed data with 5/3 length filters. PCI sensing matrix and the proposed wavelet decomposition strategy have been used to generate these results.

Sampl- ing ratio	Wavelet used	Beads	Boat	Lena	Ballo- ons	Bar- bara	Pepp- ers	Mand- ril	House	Build- ing	Camer- aman	Average PSNR
90	db2	36.11	38.42	40.20	49.72	35.79	38.30	30.41	52.52	34.90	43.55	39.99
	db4	37.23	39.47	41.47	51.00	37.77	38.61	30.71	55.72	35.03	45.53	41.25
	bior 5/3	35.87	37.77	39.64	49.38	35.05	37.68	29.67	52.66	34.42	43.16	39.53
	Matched	41.02	41.07	43.35	52.41	36.40	38.02	31.63	60.52	37.30	50.54	43.23
80	db2	32.35	34.55	36.47	46.30	31.82	35.02	27.09	47.55	31.43	38.96	36.15
	db4	33.51	35.54	37.69	47.19	33.52	35.27	27.35	50.77	31.54	40.91	37.33
	bior 5/3	32.10	33.93	35.84	45.67	31.12	34.40	26.35	48.02	30.96	38.67	35.71
	Matched	37.16	37.47	39.67	48.96	32.43	35.05	28.29	55.84	33.76	46.09	39.47
70	db2	29.87	31.98	33.97	43.46	29.20	32.90	25.04	44.10	29.30	35.74	33.56
	db4	30.93	32.82	35.10	44.50	30.67	33.14	25.27	46.94	29.32	37.51	34.62
	bior 5/3	29.51	31.31	33.33	43.12	28.51	32.30	24.28	44.43	28.80	35.22	33.08
	Matched	34.57	35.08	37.20	46.52	29.71	33.27	26.23	52.37	31.51	42.53	36.90
60	db2	27.83	29.96	31.94	41.19	27.16	31.17	23.50	40.98	27.53	33.22	31.45
	db4	28.65	30.65	32.98	42.10	28.36	31.41	23.68	43.26	27.48	34.55	32.31
	bior 5/3	27.36	29.23	31.26	40.93	26.49	30.61	22.77	41.11	27.03	32.39	30.92
	Matched	32.46	33.22	35.29	44.63	27.55	32.00	24.68	49.08	29.77	39.73	34.84
50	db2	25.96	27.98	29.97	39.12	25.48	29.70	22.17	37.90	25.90	30.67	29.48
	db4	26.70	28.60	30.85	39.75	26.29	29.91	22.32	39.64	25.81	31.74	30.16
	bior 5/3	25.52	27.41	29.29	38.90	24.83	29.22	21.51	37.94	25.42	29.76	28.98
	Matched	30.59	31.52	33.47	42.73	25.86	30.94	23.38	45.73	28.34	37.00	32.96
40	db2	24.14	26.07	28.00	36.48	23.87	28.03	20.95	34.71	24.15	28.14	27.45
	db4	24.75	26.51	28.75	36.64	24.45	28.25	21.04	35.57	24.04	28.98	27.90
	bior 5/3	23.75	25.57	27.50	36.57	23.36	27.71	20.37	34.72	23.77	27.17	27.05
	Matched	28.78	29.82	31.81	40.34	24.43	29.93	22.25	42.33	27.00	34.36	31.11
30	db2	21.93	24.04	25.88	33.33	22.18	26.25	19.69	30.69	22.24	25.39	25.16
	db4	22.51	24.34	26.41	32.68	22.53	26.39	19.69	30.96	22.06	25.99	25.36
	bior 5/3	21.76	23.66	25.53	33.66	21.80	25.97	19.23	30.72	21.90	24.71	24.89
	Matched	27.05	28.22	30.06	37.06	23.21	28.93	21.18	38.50	25.64	31.75	29.16
20	db2	18.84	21.07	22.85	28.74	19.66	23.00	17.81	25.36	19.56	21.84	21.87
	db4	19.59	21.42	23.27	28.01	20.03	23.60	17.89	25.83	19.61	22.32	22.16
	bior 5/3	18.82	20.81	22.65	29.01	19.43	23.03	17.51	25.37	19.27	21.46	21.74
	Matched	24.98	26.32	28.25	33.05	22.18	27.70	20.17	34.41	24.03	28.80	26.99
10	db2	10.75	9.79	10.03	13.47	9.74	8.20	8.95	9.60	10.24	10.35	10.11
	db4	11.59	10.61	11.06	14.53	10.55	9.28	9.59	10.73	10.87	11.28	11.01
	bior 5/3	11.01	10.19	10.46	13.77	10.15	8.62	9.35	10.08	10.67	10.61	10.49
	Matched	20.00	22.29	23.74	30.95	20.08	23.04	18.08	27.69	20.26	23.03	22.92

Equations

$$\mathbf{x}=\mathbf{\Psi s}$$

$$s_j=C{j}^{-q},j=1,2,...,N,$$

$$y[i]=\sum _{j=0}^{N-1}{\varphi }_{i,j}x[j],i=0,1,...,M-1,\text{vspace}-0.5em$$

$${\mathbf{y}}_{M\times 1}=$$

$$=$$

$$M\ge CK\log (N/K)$$

$$\tilde{\mathbf{s}}=\underset{\mathbf{s}}{\min }||\mathbf{s}|{|}_0\text{ subject to: }\mathbf{y}=\mathbf{As}.$$

$$\tilde{\mathbf{s}}=\underset{\mathbf{s}}{\min }||\mathbf{s}|{|}_1\text{ subject to: }\mathbf{y}=\mathbf{As},$$

$$\mathbf{y}$$

$$=\mathbf{\Phi x}$$

$$H_1^{new}(z)=H_1(z)-H_0(z)T(z^2),$$

$$F_0^{new}(z)=F_0(z)+F_1(z)T(z^2).$$

$$H_0^{new}(z)=H_0(z)+H_1(z)S(z^2),$$

$$F_1^{new}(z)=F_1(z)-F_0(z)S(z^2).$$

$$\mathbf{\Phi}^{p}_{i,j}=\left\{\begin{array}[]{ll}1\quad\mbox{ if }i\in\{1,2,...,M\}\mbox{ and }j\in\Omega\\ 0\quad\mbox{otherwise}\end{array}\right.$$

$$\text{PSNR (in dB)}=10{\log }_{10}\left(\frac{\sum _{i=0}^{m-1}\sum _{j=0}^{n-1}|\text{I}(i,j)-\hat{\text{I}}(i,j)|}{\sum _{i=0}^{m-1}\sum _{j=0}^{n-1}|\text{I}(i,j){|}^2}\right),$$

$$\tilde{\mathbf{s}}=\underset{\mathbf{s}}{\min }||\mathbf{s}|{|}_1\text{ subject to: }\mathbf{y}=\mathbf{\Phi \Psi }\mathbf{s}.$$

$${\tilde{d}}_{-1}[n]$$

$$={\tilde{x}}_o[n]-p[n]$$

$$=\tilde{x}[2n+1]-p[n]$$

$$p[n]=\sum _{k=0}^{L_t-1}t[k]{\tilde{x}}_e[n-k],$$

$$T(z)=z^{-(\frac{L_t}{2}-1)}\sum _{i=0}^{L_t-1}t[i]z^i,$$

$$p[n]=\sum _{k=0}^{L-1}t[k]\tilde{x}[2n-L_t+2+2k].$$

$$p[n]=$$

$$+t[\frac{L_t}{2}]\tilde{x}[2n+2]+...+t[L_t-1]\tilde{x}[2n+L_t].$$

$${\tilde{\mathbf{d}}}_{-1}={\mathbf{d}}_{-1}+{\mathbf{\eta }}_1,$$

$${\hat{\mathbf{d}}}_{-1}={\mathbf{d}}_{-1}+{\mathbf{\eta }}_2,$$

$${\hat{\mathbf{d}}}_{-1}={\tilde{\mathbf{d}}}_{-1}+\mathbf{\eta },$$

$${\hat{\mathbf{d}}}_{-1}=\mathbf{At}+\mathbf{\eta },$$

$$a_{-1}[n]=x_e[n]+{\tilde{d}}_{-1}[n]\ast s[n],$$

$$S(z)=z^{(\frac{L_s}{2}-1)}\sum _{i=0}^{L_s-1}s[i]z^{-i},$$

$${\tilde{d}}_{-1}[n]=-t[0]\tilde{x}[2n]+\tilde{x}[2n+1]-t[1]\tilde{x}[2n+2]$$

Full text

I Introduction

Classical signal acquisition technique in signal processing branch involves sensing the full signal at or above the Nyquist rate. In general, this signal is transformed to a domain where it is compressible. Only some of the largest coefficients of transformed signal having sufficient amount of energy are stored and transmitted to the receiver end along with the position information of the transmitted coefficients. The transmitted signal is decoded at the receiver end to recover the original signal. Thus, this process involves sensing the full signal, although most of the samples in the transformed domain are to be discarded.

In [1, 2], researchers proposed compressed sensing (CS) method that combines sensing and compression in one stage, where instead of sampling the signal sample-wise above Nyquist rate, it’s projections are captured via the measurement basis. These samples are very few as compared to those sampled at Nyquist rate. If the signal is sparse in some transform domain and also if the measurement basis are incoherent with the sparsifying basis, then the whole signal can be reconstructed with a good probability from a very few projections of the signal.

In this context, wavelets are extensively used as sparsifying basis in compressed sensing problem [2]. One of the advantage with wavelets is that there is no unique basis unlike Fourier transform. One may choose the set of basis depending upon the type of application. Since the wavelet basis are not unique, it is better to design wavelet that are matched to a given signal in a particular application. The designed wavelet basis are called signal-matched wavelet [3, 4, 5].

Motivated with the above discussion, this paper proposes to design matched wavelets for compressive sensing application. Unlike previous works on matched wavelets [3, 4, 5, 6, 7, 8, 9] where wavelets are designed from fully sampled signal, this paper proposes a novel method of designing signal matched wavelets from compressively sensed images at the receiver. The proposed method employs lifting based framework, also called second generation wavelet framework, to design image-matched wavelet [10].

Many researchers have designed wavelets using lifting [11, 8, 12, 13, 14, 15, 16][5, 17] that requires design of predict and update stage filters. For images-matched wavelet design, these filters are required to be designed from the given image itself. Although it is easy to design predict stage filters, design of update stage filters offer a real challenge. One of the criteria used in the literature to design the update stage filter is to minimize the reconstruction error[11]. In [8], update step is chosen based on the local gradient of the signal. In [12], authors use non-separable lifting based approach of wavelet implementation on images with regularity conditions imposed. In [13], directional interpolation in update stage is used wherein coefficients of interpolation filter adapt to the statistical property of image. In [14], allpass filter in the lifting stage is used to design orthogonal IIR (Infinite Impulse Response) filterbank. In [15], geometry of the image is used to design wavelet via lifting that leads to local and anisotropic filters. Authors have designed nonseparable filterbanks that are pixel-wise adaptive to the local image feature in [16]. In [5], authors designed predict stage by minimizing the energy in the wavelet subspace domain and designed update stage by minimizing the difference between reconstructed and original signal leaving wavelet signal.

However, all the methods discussed above design signal-matched wavelet using the full signal at the input. So far, to the best of our knowledge, no method has been proposed that designs signal or image-matched wavelets from compressive or partial measurements. This paper addresses this problem from the point of view of image reconstruction. Below are the salient contributions of this paper:

We propose design of image-matched separable wavelet in the lifting framework from compressively sensed image, that is also used to reconstruct the image. 2. 2.

In general, Gaussian or Bernoulli measurement matrices are used in compressive sensing application [2],[18]. We propose to use partial canonical identity matrix to sample data at sub-Nyquist rate such that sensing time is reduced considerably. 3. 3.

For the separable 2D wavelet transform, a new multi-level wavelet decomposition strategy is proposed that leads to improved reconstruction performance. We name this new wavelet decomposition strategy as multi-level L-Pyramid wavelet decomposition.

This paper is organized as follows: In section-II, we briefly present the theory on compressive sensing and lifting-based wavelet design. In section-III, we propose to use partial canonical identity (PCI) matrix to sub-sample images at sub-Nyquist rate. We show results on time complexity and performance of the conventionally used sensing matrices and PCI matrix to establish the use of latter in this work. In section-IV, we propose a new multi-level wavelet decomposition strategy for separable wavelet on images and compare the performance with the existing strategy. Section-V presents the proposed joint framework of matched wavelet estimation and image reconstruction in CS. We present experimental results in section-VI and conclusions in section-VII.

II Background

In this section, we briefly present the theory of compressive sensing and lifting framework of wavelets for the sake of self-completeness of the paper.

II-A Compressed Sensing

Classical compression method involves two steps: sensing and compression wherein, first, an analog data is sampled at or above the Nyquist-rate and then, it is compressed by an appropriate transform coding process. In general, natural signals are sparse or compressible in some transform domain. For examples, if a signal is smooth, it is compressible in Fourier domain and if it is piece-wise smooth, it is sparse in the wavelet domain. To understand this process, let us consider a signal $\mathbf{x}$ of dimension $N\times 1$ that has been sensed by a traditional sensing technique at or above the Nyquist rate. This signal is next transformed to sparse domain with the help of a sparsifying basis $\psi_{i},i=1,2,...,N$ and represented as:

[TABLE]

A signal s is $K-$ sparse if all but $K$ elements are zero, whereas a signal is compressible if its sorted coefficients obey the power low decay [19]

[TABLE]

where $s_{j}$ represent the sorted coefficients and $q$ represents decay power parameter. For large value of $q$, decay of coefficients is faster and correspondingly, signal is more compressible. In compression, some of the largest coefficients of the transformed signal are kept and all other coefficients are discarded. These coefficients along with their location information are sent to the receiver. Having the knowledge of the sparsifying basis and signal coefficients along with their positions in the original signal, signal is reconstructed back at the receiver end.

The above process consisting of first sensing the whole signal and then discarding many of its transform domain coefficients is inefficient. Compressive sampling or sensing [1, 2, 20] combines these two processes. Instead of sampling the signal at or above the Nyquist rate, signal’s linear projection on some measurement basis $\phi_{i}$ are obtained. If $\phi_{i}$ is the $i^{th}$ measurement basis, then $i^{th}$ observation of the projected signal is given by:

[TABLE]

where $M$ is the number of linear projections of the signal. In compact form, this can be written as:

[TABLE]

where $i^{th}$ measurement basis is stacked as a row of the matrix $\mathbf{\Phi}$ and $\mathbf{A}=\mathbf{\Phi\Psi}$. CS theory states that the original signal of length $N$ can be recovered with very high probability, if number of linear projections $M$ are taken such that [1]:

[TABLE]

where $K$ is the sparsity of the signal, C is some constant, and $M\ll N$ in general.

Equation (4) represents under-determined system of linear equations with $\mathbf{y}=\bm{\Phi}\mathbf{\hat{x}}$ having infinite many solutions $\mathbf{\hat{x}}$. But if the signal is sparse in some transform domain $\mathbf{\Psi}$, (4) can be solved using $l_{0}$ minimization as below:

[TABLE]

The above problem is NP-hard to solve. It has been shown in [21] that $l_{1}$ minimization

[TABLE]

provides the same solution as $l_{0}$ minimization. Here $||\mathbf{v}||_{1}$ denotes the $l^{1}$ norm or sum of the absolute values of the vector $\mathbf{v}$. $l^{1}$ minimization is known as Basis Pursuit (BP) in literature and can be solved by linear programming [22].

Compressive sensing is being used increasingly in image reconstruction, also called as compressive imaging (CI). For example, let us consider an image $\mathbf{X}$ of dimension $m\times n$, that is compressively sensed by a measurement matrix $\mathbf{\Phi}$. These measurements are given by

[TABLE]

where $vec(\mathbf{X})=\textbf{x}$ denotes the vector of length $N=mn$ of image $\mathbf{X}$ and the measured signal $\mathbf{y}$ is of dimension $M\times 1$, where $M$ is the number of compressive measurements. It has been observed that natural images, in general, are compressible in DCT (discrete cosine transform) [23] and wavelet domains [24]. Hence, DCT or wavelet can be applied as separable transforms on images and used as sparsifying basis $\mathbf{\Psi}$ in (4) in CS-based image reconstruction.

The measurement or the sensing basis $\mathbf{\Phi}$ can be chosen such that it satisfies Restricted Isometry Property (RIP) [25] and coherency property [2]. Some of the examples of the measurement matrices that satisfy these properties are random matrices with entries taken from i.i.d. Gaussian distribution [26], random matrices with entries taken from uniform Bernoulli distributions [20], and Fourier matrix [2]; although several other structured measurement matrices such as toeplitz and circulant matrices are also being used [27, 28, 29].

II-B Lifting Theory

Lifting is a technique for either factoring existing wavelet filters into a finite sequence of smaller filtering steps or constructing new customized wavelets from existing wavelets [30]. Lifting-based wavelet design is also known as second generation wavelet design. The design is modular, guarantees perfect reconstruction at every stage, and supports non-linear filters. A general lifting scheme consists of three steps: Split, Predict, and Update (Refer to Fig. 1).

Split: In the split step, input signal is split into two disjoint sets of samples, generally, even and odd indexed samples, labeled as $x_{e}[n]$ and $x_{o}[n]$, respectively. The original signal can be recovered perfectly by interlacing or combining the two sample streams. The corresponding filterbank is called as the Lazy Wavelet system [10] and is similar to the structure shown in Fig. 2 with analysis filters labeled as $H_{0}(z)=Z\{h_{0}[n]\}$, $H_{1}(z)=Z\{h_{1}[n]\}$ and synthesis filters as $F_{0}(z)=Z\{f_{0}[n]\}$, $F_{1}(z)=Z\{f_{1}[n]\}$.

Predict Step: In the predict step, also known as dual Lifting step, one of the two disjoint set of samples is predicted from the other set of samples. For example, in Fig. 1(a), we predict odd set of samples from the neighboring even samples by using the predictor $P\equiv T(z)$. Predict step is equivalent to applying a highpass filter on the input signal. Predict step modifies the highpass filter of the analysis end and lowpass filter of the synthesis end, without altering other filters, according to the following relations:

[TABLE]

[TABLE]

Update Step: In the update step, also known as primal lifting step, predicting samples of the predict step are updated with the predicted samples to provide the approximate coefficients of the signal. The signal is updated with $U\equiv S(z)$ (refer to Fig. 1). This step modifies the analysis lowpass filter and synthesis highpass filter according to the following relation:

[TABLE]

[TABLE]

Once all the filters are designed, Fig. 1 can be equivalently drawn as Fig. 2 or any existing wavelet system of Fig. 2 can be equivalently broken into lifting steps of Fig. 1. One of the major advantages of lifting scheme is that each stage (predict or update) is invertible. Hence, perfect reconstruction (PR) is guaranteed after every step.

III Compressive Sensing of Images using Proposed Sensing Matrix

In this section, first we establish the need for a different, other than conventional, sensing matrix. Next, we discuss the proposed matrix. And later, we present results to show the comparison of time complexity and reconstruction performance with the proposed matrix in CS based image reconstruction.

III-A Context

In today’s world, size of images are increasingly large and $N$ generally approaches to millions of samples. This large size imaging poses challenges for CS-based image reconstruction. First challenge is the huge size of measurement matrix $\Phi$ that poses problems with storage and computation. Other problems include design of imaging system with larger space bandwidth product (SBP) and difficult calibration requirements [31].

In an attempt to overcome the above challenges, single pixel camera hardware architecture has been proposed in [18]. It replaces the traditional camera architecture and captures the inner product between the scene under view and measurement basis. Thus, the camera captures one pixel at a time that is a linear combination of all pixel samples of the image. This process is repeated $M$ number of times with $M\ll N$. These are called the compressive measurements and are transmitted to the receiver where full sized image is reconstructed by employing the theory of CS-based reconstruction. For more information on single pixel camera, reader may refer to [18].

With the above explained architecture, a single pixel camera replaces the photon detector array of a traditional camera by a single photon detector; thus, reducing the size, cost, and complexity of the imaging architecture.

Although the above hardware architecture implements compressive imaging (CI) nicely, it suffers with some difficulties including sensor dynamic range, A/D quantization, and photon counting effects [18]. Also, this process is time consuming as one has to wait for $M$ samples that are captured sequentially. This is a serious problem in real-time applications, say, when one has to record a video using camera as the scene may change while capturing samples sequentially of the current scene under view. Also, since $M$ linear projections are captured instead of $N$ pixel samples, it “effectively” samples the image at Nyquist rate instead of actually sampling the image at the sub-Nyquist rate.

III-B Proposed Use of Partial Canonical Identity Sensing Matrix

We propose to use PCI sensing matrix that as per our knowledge, is the simplest sensing matrix proposed so far and “actually” senses the image at sub-Nyquist rate by capturing less number of pixels without sensing information about every pixel. This is explained as below.

Consider an image $\mathbf{X}$ of dimension $m\times n$. Instead of sampling all the $N(N=mn)$ pixels of the image using the sensor array of the traditional camera, we capture $M$ samples of the image using the proposed measurement matrix $\mathbf{\Phi}^{p}$, where $M\ll N$. The measurement matrix $\mathbf{\Phi}^{p}$ has the entries shown below:

[TABLE]

where $\Omega\subset\{1,2,...,N\}$ such that $|\Omega|=M$ and $|.|$ denotes the cardinality of the the set. We name the proposed sensing matrix as partial canonical identity (PCI) matrix because it consists of partially selected and permuted rows of the identity matrix.

The PCI matrix captures only $M$ samples of the actual image; thus, actually sub-samples the original image. This can be accomplished by using existing cameras by switching ON only $M$ sensors of the sensor array. This is unlike the single pixel camera where every captured pixel is the linear combination of the entire image pixel set. Also, in single pixel camera, one has to wait for $M$ units of time to sense $M$ number of samples, whereas all $M$ samples are sensed in one unit of time in the case of PCI matrix. Thus, PCI matrix reduces the sensing time by a factor of $M$ in comparison to a single pixel camera.

III-C Results using PCI Sensing Matrix

Consider a sub-sampled Lena image, shown in Fig. 3(a), with only 50% samples captured via PCI sensing matrix. The un-captured positions are filled with zeros. Fig. 3(b) shows the image reconstructed from this sub-sampled image using (7) with standard wavelet ‘db4’ as the sparsifying basis. From Fig. 3(b), we note that PCI matrix has good ability of reconstructing images from compressively sensed data and can be used in CS-based image reconstruction. Fig. 4 and 5 provide detailed results.

Fig. 4 compares the time taken in image reconstruction from the measured samples with sampling ratios varying from 10% to 90% using different measurement matrices. We compare the reconstruction time taken using the PCI matrix, random Gaussian matrix, and Bernoulli random matrix with $\pm 1$ as its entries. Random Gaussian matrix is preferred in wide range of applications because of the ease in theoretical analysis, while the Bernoulli matrix depicts physical implementation of single pixel camera. Image is reconstructed using equation (7), where we have used standard Daubechies orthogonal wavelet ‘db4’ as the sparsifying basis. We have used MATLAB solver spgl1 [32], [33] to solve (7) that implements Basis Pursuit (BP) [22].

Compressed sensing based reconstruction with Gaussian and Bernoulli sensing matrices is implemented using block compressed sensing [34]. This is to note that reconstruction with PCI matrix requires only the position information of the sampled pixels instead of the information of all entries of $M\times N$ sensing matrix that simplifies reconstruction with PCI matrix. We compare reconstruction results on three images: ‘Beads’, ‘Lena’, and ‘House’ as shown in Fig. 12(a), 12(c) and 12(h), respectively. Each image is of size $512\times 512$. We have chosen these images because they exhibit different spectral properties. For example, ‘Beads’ is rich in high frequencies, ‘Lena’ contains both low and high frequency contents, while ‘House’ is rich in low frequencies.

From Fig. 4, we note that the reconstruction time with Gaussian and Bernoulli sensing matrices is almost the same, whereas reconstruction time with the PCI matrix is extremely less. This huge reduction in time is owing to the implementation simplicity of PCI matrix. However, there is a trade-off between the reconstruction time and the accuracy. Fig. 5 compares reconstruction accuracy of images in terms of peak signal-to-noise ratio (PSNR) given by:

[TABLE]

where I is the reference image, $\hat{\textbf{I}}$ is the reconstructed image, and $m\times n$ is the size of the image. Results shown are averaged over 10 iterations. From Fig. 5, we observe that reconstruction accuracy with PCI matrix is approx. 2 dB lower than that with Bernoulli or Gaussian sensing matrices. Similar to this observation, in [35], authors have shown comparatively inferior CS-based reconstruction results with PCI matrix. In section-VI, we will show that this gap in reconstruction accuracy can be not only bridged, but rather enhanced to a large extent with matched wavelet. Since the use of PCI matrix expedites the complete pipeline, its use is overall useful, if a little bit inferior quality can be bridged.

IV Proposed Wavelet Decomposition Method on Images

In this Section, we propose a new optimized strategy of multi-level wavelet decomposition on images.

A separable wavelet transform is implemented on images by first applying 1-D wavelet transform along the columns and then along the rows of an image. This provides 1-level wavelet decomposition that consists of four components labeled as LL, LH, HL and HH, respectively. The same procedure is repeated on the LL part of the wavelet transform $k$-times to obtain $k$-level decomposition of an image (Fig. 6(a)). We call this decomposition as regular pyramid (R-Pyramid) wavelet decomposition.

In general, $k$-level wavelet decomposition of an image consists of following components: LL${}_{\mbox{{k}}}$, LH${}_{\mbox{{i}}}$, HL${}_{\mbox{{i}}}$ and HH${}_{\mbox{{i}}}$, where $i=1,2,...,k-1$. LH${}_{\mbox{i}}$, HL${}_{\mbox{i}}$ and HH${}_{\mbox{i}}$ components are obtained by applying wavelet transform on the columns and rows of LL${}_{\mbox{{i}-1}}$ component. LH${}_{\mbox{{i}}}$ is obtained by filtering LL${}_{\mbox{{i}-1}}$ column-wise using a lowpass filter and filtering it row-wise using a highpass filter. Thus, the current nomenclature of labeling subbands is: first character represents operation on columns and second character represents operation on rows, where operation implies highpass or lowpass filtering denoted by symbols ‘H’ and ‘L’, respectively.

In the conventional 2-D wavelet transform (Fig. 6(a)), wavelet decomposition is applied on LL${}_{\mbox{{i}}}$ part only to obtain the $(i+1)^{th}$ level coefficients. Since it is a separable transform, similar to 1-D wavelet transform wherein wavelet is applied repeatedly on lowpass filtered branches, we propose to apply wavelet in the lowpass filtered directions of LH${}_{\mbox{{i}-1}}$ and HL${}_{\mbox{{i}-1}}$ subbands in contrast to the conventional decomposition strategy wherein these subbands are left unaltered. Thus, the proposed 2nd level wavelet decomposition is as shown in Fig. 6(c).

Since we apply wavelet in only one direction of LH${}_{\mbox{{i}-1}}$ and HL${}_{\mbox{{i}-1}}$ subbands, we notate these subbands differently compared to the conventional scheme. We assign subscript with both ‘L’ and ‘H’ symbols of every subband to denote the no. of times wavelet has been applied in that direction. In order to understand this, let us first consider 1-level wavelet decomposition as shown in Fig. 6(b) that is similar to the conventional scheme shown in Fig. 6(a). However, the subbands are labeled as L${}_{\mbox{1}}$L${}_{\mbox{1}}$, L${}_{\mbox{1}}$H${}_{\mbox{1}}$, H${}_{\mbox{1}}$L${}_{\mbox{1}}$, and H${}_{\mbox{1}}$H${}_{\mbox{1}}$.

In the 2nd level wavelet decomposition, wavelet is applied in both directions of L${}_{\mbox{1}}$L${}_{\mbox{1}}$ subbands leading to L${}_{\mbox{2}}$L${}_{\mbox{2}}$, L${}_{\mbox{2}}$H${}_{\mbox{2}}$, H${}_{\mbox{2}}$L${}_{\mbox{2}}$, and H${}_{\mbox{2}}$H${}_{\mbox{2}}$ subbands. But in addition, wavelet is applied on the columns of L${}_{\mbox{1}}$H${}_{\mbox{1}}$ yielding two subbands L${}_{\mbox{2}}$H${}_{\mbox{1}}$ and H${}_{\mbox{2}}$H${}_{\mbox{1}}$. Also, wavelet is applied on the rows of H${}_{\mbox{1}}$L${}_{\mbox{1}}$ subband yielding two subbands H${}_{\mbox{1}}$L${}_{\mbox{2}}$ and H${}_{\mbox{1}}$H${}_{\mbox{2}}$. Applying similar strategy for the 3rd level decomposition, we obtain subbands as shown in Fig. 6(d). We name this wavelet decomposition as L-shaped pyramid (L-Pyramid) wavelet decomposition.

The efficacy of the proposed wavelet strategy is shown in CS-based image reconstruction using orthogonal Daubechies wavelet ‘db4’. Fig. 7 shows reconstruction accuracy in terms of PSNR (14) with sampling ratios ranging from 10% to 90% averaged over 10 iterations. We compare reconstruction accuracy with the existing wavelet decomposition strategy and with the proposed strategy on the same three images: ‘Beads’, ‘Lena’ and ‘House’. From Fig. 7, we note marked improvement with the proposed strategy over the existing strategy. Performance is particularly improved for image ‘House’ that is rich in low frequencies.

V Proposed Method of Designing Matched Wavelet from Compressively Sensed Images

As discussed earlier, wavelets are extensively used as sparsifying transform for image reconstruction in CS domain. Researchers often face challenge in choosing a wavelet because it is not known apriori as to which wavelet will provide better representation of the signal. In general, one uses compactly supported wavelets, either orthogonal Daubechies wavelets or biorthogonal 9/7 or 5/3 wavelets (note that first digit denotes the length of analysis lowpass filter and second digit denotes the length of analysis highpass filter). It is natural to think that a wavelet matched to a given signal will provide best representation of the given signal and hence, will provide better reconstruction in CS compared to any wavelet chosen arbitrarily. Although a number of methods exist for designing wavelet matched to a given signal [3, 4, 5, 6, 7, 8, 9], these methods work when full signal is available. On the other hand, in CS, only compressively sensed signal is available at the receiver. No method exists for estimation of matched wavelet from compressively sensed signal that can be utilized at the same time for efficient reconstruction of signal/image from compressively sensed signal/image. In this Section, we address this problem. We propose a joint framework for signal reconstruction in CS wherein we are estimating wavelet from the compressivley sensed image and use it for efficient image reconstruction at the same time.

Since the proposed work is on separable wavelets, we require to estimate matched wavelet for row and column directions separately. Thus, before proceeding with the matched wavelet design, we present the scanning mechanism of rows and columns data in images as used in this work.

V-A Scanning Mechanism for Row- and Column-wise Data

As stated earlier, we require to estimate matched wavelet for both the row and column directions. One easier strategy can be designing matched wavelet on row- or column-vectorized image and use the same wavelet, later, along both the columns and rows as a separable wavelet. Instead, we propose to design matched wavelet separately for the row- and column-directions. These directions can be scanned using the following two strategies:

• •

Strategy-1: We scan the image according to the scanning scheme as shown in Fig. 8, wherein rows or columns are stacked one after the other to obtain 1-D signal for both the directions. However, this will cause discontinuity in the 1-D signal at the transitions when one column ends and another starts and likewise, for the rows.

• •

Strategy-2: In order to avoid this discontinuity, other alternate strategy is to scan all even rows or columns in the reverse direction as shown in Fig. 9.

Since strategy-2 is robust to sudden transitions at the row- or column-endings, we use it in all our experiments.

V-B Proposed Methodology of Matched Wavelet Design

With the scanning strategy-2 proposed above, we convert a given image into two 1-D signals: one with column-wise scanning and another with row-wise scanning. Henceforth, in the present Section, we present matched wavelet design methodology using compressively sensed 1-D signal. We will use these designs in the next Section to note the performance over images as separable wavelet transforms.

The proposed methodology has three stage. In stage-1, we obtain coarse image estimate from compressively (partially) sensed samples using a standard wavelet. We call this a coarser estimate because the wavelet used is not matched to the given signal and hence, the original signal may not be that sparse over this wavelet compared to that with the matched wavelet. This will impact the reconstruction performance. In stage-2, we estimate matched analysis wavelet filter that provides sparser subband wavelet (detail) coefficients than those obtained from the standard wavelet in stage-1. Using these estimated filters and the coarser signal estimate of stage-1, we design all filters of the matched wavelet system. In stage-3, we reconstruct signal from measured sub-samples using the matched wavelet estimated in stage-2.

V-B1 Stage-1: Coarser Image Estimation

In this stage, we reconstruct a coarser estimate of the signal from compressively (partially) sensed measurement data y using Basis pursuit (BP) optimization method as below:

[TABLE]

where $\mathbf{\Psi}$ corresponds to any standard wavelet. We use biorthogonal 5/3 wavelet. The coarser approximation of the signal is obtained as $\tilde{\mathbf{x}}=\mathbf{\Psi}\mathbf{\tilde{s}}$. We solve the above optimization problem using MATLAB solver spgl1 [32, 33].

V-B2 Stage-2: Estimation of matched wavelet

We use the signal reconstructed in the previous step to design matched wavelet system by designing Predict and update stage filters of the lifting structure.

**(a) Predict Stage: ** For designing signal-matched analysis wavelet filter in the lifting framework, first, we consider Lazy wavelet with $H_{0}(z)=1$, $H_{1}(z)=z$, $F_{0}(z)=1$, and $F_{1}(z)=z^{-1}$, in Fig. 2. We start with this Lazy wavelet structure and proceed with designing the predict stage filter $T(z)$.

We apply coarser version of original signal $\tilde{\textbf{x}}$ estimated in Stage-1 as input to the Lazy wavelet filterbank and obtain even and odd sampled streams $\tilde{\textbf{x}}_{e}(n)$ and $\tilde{\textbf{x}}_{o}(n)$, respectively as shown in Fig. 1(a). We pass even indexed samples $\tilde{\textbf{x}}_{e}(n)$ through the predict stage filter and write the output of the lower subband signal, shown in Fig. 1(a), as below:

[TABLE]

where

[TABLE]

where ‘$\ast$’ is the convolution operator, $L_{t}$ is the length of the predict stage filter $t[n]$ with its Z-transform given by $T(z)=Z\{t[n]\}$. For good prediction, a sample (here odd indexed) should be predicted from its immediate past and immediate future neighbors that requires a careful choice on the predict stage filter provided by Theorem-1 as below.

Theorem 1**.**

The following structure of predict stage filter allows odd-indexed samples to be predicted from their nearest even-indexed samples, i.e., from their immediate past and immediate future samples:

[TABLE]

where $T(z)$ is an even-length filter. Even length $T(z)$ will ensure that equal number of past and future samples are used in prediction.

Proof.

On using (18) in (17), we obtain:

[TABLE]

This is to note that the $n^{th}$ sample of $p[n]$ predicts sample $\tilde{x}[2n+1]$. On expanding (19), we obtain

[TABLE]

Thus, the sample $\tilde{x}[2n+1]$ is predicted from its exact nearest even neighbors. For better clarity, Fig. 10 shows even neighboring samples that are being used to predict odd samples with a 2-tap and 4-tap predict filter $T(z)$. This proves the theorem. ∎

Subband signal $\tilde{\mathbf{d}}_{-1}$ can be considered as the noisy version of detail coefficients $\mathbf{d}_{-1}$ that could be obtained by passing the original signal x through the signal-matched analysis wavelet branch. This can be written as:

[TABLE]

where $\bm{\eta}_{1}$ is the corresponding error.

Further, we analyze this coarser approximate signal $\tilde{\textbf{x}}$ using the same biorthogonal 5/3 wavelet with which it had been reconstructed and obtain detail coefficients (subband coefficients of highpass filtered branch) $\hat{\mathbf{d}}_{-1}$. Again, this signal can be considered as the noisy version of detail coefficients $\mathbf{d}_{-1}$ and can be written as below:

[TABLE]

where $\bm{\eta}_{2}$ is the corresponding error.

From (21) and (22), we obtain:

[TABLE]

where $\bm{\eta}$ is the error that consists of two components: 1) because we are using approximate signal $\tilde{\mathbf{x}}$ instead of the original signal x and 2) because we are using biorthogonal 5/3 wavelet instead of the signal-matched wavelet.

On substituting for $\tilde{\mathbf{d}}_{-1}$ from (16), (17), and (18) in (23) and writing in the matrix form, we obtain

[TABLE]

where $\mathbf{A}$ is the convolution matrix consisting of even and odd indexed samples of $\tilde{\textbf{x}}$ and t denotes the vectorized form of predict stage filter $t[n]$ or $T(z)$. We solve for t in (24) using least squares method and substitute in (9) and (10) to update the analysis highpass and synthesis lowpass filters and obtain new filters $H_{1}^{new}(z)$ and $F_{0}^{new}(z)$, respectively. This ends the predict stage.

**(b) Update Stage: ** In the update stage, update polynomial $S(z)$ is required to be computed. In order to do this, we write the output of the upper subband signal using the lower subband signal, $\tilde{d}_{-1}[n]$ (refer Fig. 1) as below:

[TABLE]

where $s[n]$ is the time domain description of the update stage filter $S(z)$.

Similar to the predict stage filter, we require to choose $s[n]$ such that the elements of the upper branch are updated using nearest neighbors only. The corresponding structure for $s[n]$ is provided by Theorem-2 as below.

Theorem 2**.**

The following structure of the update stage filter allows the elements of the upper branch to be updated from nearest neighbors:

[TABLE]

where $L_{s}$ is the length of filter $S(z)$ or $s[n]$. Note that $S(z)$ is an even-length filter that ensures that sample update is done using equal number of past and future samples.

Proof.

For the sake of simplicity, let us consider two tap predict stage filter that provides the following detail coefficients:

[TABLE]

Once this signal is passed through the update stage filter in the update branch, we obtain:

[TABLE]

where ‘$\ast$’ is the convolution operator. With the choice of filter in (26), we obtain:

[TABLE]

On expanding (V-B2) and rearranging, we obtain

[TABLE]

The above signal updates approximate coefficients, i.e., $\tilde{x}[2n]$. It can be clearly noticed that coefficients $\tilde{x}[2n]$ are updated using the nearest neighbors from $\tilde{x}[2n+L_{s}]$ to $\tilde{x}[2n-L_{s}]$. For better clarity, Fig. 11 shows the neighboring samples that are being used to update even samples with a 2-tap predict filter $T(z)$ and a 2-tap update filter $S(z)$. ∎

This subband signal $a_{-1}[n]$ is passed through a 2-fold upsampler that provides:

[TABLE]

Next, signal $\tilde{x}_{1u}[n]$ is passed through the synthesis lowpass filter $f^{new}_{0}[n]$ that was updated in the predict stage mentioned earlier. This provides us the signal $\tilde{x}_{1}[n]$ (shown in Fig. 2) reconstructed from the upper subband only and is given by

[TABLE]

Assuming that the original signal of interest is rich in low frequency content, signal $\tilde{x}_{1}[n]$ reconstructed in the upper subband should be in close approximation to the input signal $\tilde{\textbf{x}}$. This allows us to solve for the update stage filter as below:

[TABLE]

It can be noted from (25), (31) and (32), that $\tilde{\textbf{x}}_{1}$ can be written in terms of update stage filter $s[n]$ obtained on solving (33) using least squares method. Correspondingly, analysis lowpass filter $H_{0}(z)$ and synthesis highpass filter $F_{1}(z)$ are updated to $H^{new}_{0}(z)$ and $F^{new}_{1}(z)$ using (11) and (12), respectively. This completes our design of matched wavelet.

This is to be noted that, since the lifting scheme is modular, more number of such predict and update stage filters can be estimated and appended in order to design higher order or larger length filters.

V-B3 Stage-3: Signal Reconstruction using Matched Wavelet

Once we have estimated matched wavelet, we employ (7) on measured subsampled measured signal y using CS with matched wavelet in W to recover the original signal x.

The above matched wavelet design for 1-D signals is applied on row-scanned and column-scanned signal of a given image and corresponding matched wavelet is designed for the row- and the column-space separately and hence, completes the design of image-matched separable wavelet system. Also, this is to note that the proposed method of designing matched wavelet from compressively sensed images is independent of the sensing matrix. The proposed matched wavelet can be designed from images compressively sensed using any sensing matrix.

VI Experiments and Results

In this section, we present CS-based reconstruction results of images using image-matched wavelets designed from compressively sensed images. We apply the proposed method on ten natural images shown in Fig. 12. Images with different spectral contents have been selected. For example, Balloon, House and cameraman are rich in low frequencies, Beads, Barbara and Mandrill are rich in high frequencies, and rest of the images have varied lower and higher frequency content.

VI-A Frequency Response and Filter Coefficients

We design image-matched wavelets for all the images shown in Fig. 12. Table-I shows the analysis filter designed with the image ‘Lena’ compressively sensed at a sampling ratio of 20% with PCI sensing matrix. Coefficients for filters designed are shown in Table-I for both column matched and row matched wavelets along with the coefficients of predict and update stage filters $T(z)$ and $S(z)$, respectively.

Fig. 13 shows normalized frequency responses of 5/3 length analysis filters designed for ‘Lena’ image. These frequency response have been shown for only column-matched wavelets for brevity purposes.

Discussion: From the above figure, we note that the filters designed with the proposed method also exhibit good linearity in the phase response. This is because the proposed method of matched-wavelet design is itself able to capture linear phase characteristics of the image, if present.

VI-B Experiment-1: Comparison of existing CS-based image reconstruction methodology with the proposed methodology

In the application of CS-based image reconstruction, the proposed methodology of this paper has three contributions:

Proposed use of PCI sensing matrix: that is computationally very inexpensive compared to the existing Gaussian matrix (Fig. 4), but provides approx. 2 dB lower performance compared to the existing Gaussian matrix (Fig. 5). 2. 2.

New wavelet decomposition strategy (Fig. 6): that provides better results in CS-based image reconstruction compared to the existing wavelet decomposition strategy (Fig. 7). 3. 3.

Design of image-matched wavelets: wherein wavelets have been designed from compressively sensed images and used for images reconstruction.

Based on the above observations, we would like to compare the performance of the proposed CS-based image reconstruction (with all three novelties: Proposed use of PCI sensing matrix, proposed wavelet decomposition, and proposed matched 5/3 length wavelet design) with the existing CS-based reconstruction (Gaussian sensing matrix, existing wavelet decomposition, and dB4 wavelet). Fig. 14 shows CS-based reconstruction results in terms of PSNR averaged over 10 iterations on three images ‘Beads’, ‘Lena’, and ‘House’.

Discussion: From Fig. 14 we observe that the proposed methodology performs consistently better at all sampling ratios with the following observations:

As the image changes from being rich in high frequency (‘Beads’) to mid-frequency (‘Lena’) to low frequency (‘House’) content, better and better reconstruction performance is observed. This is owing to the fact that the proposed matched wavelet design works best for signals rich in low frequency contents. 2. 2.

At the higher sampling ratio of 90%, performance gain over standard methodology is 3 dB with ‘Beads’, 4 dB with ‘Lena’, and 13 dB with ‘House’. Since at higher sampling ratios, most of the input image samples are available upfront, hence, matched wavelet design is optimum. This provides very good performance and huge improvement over the existing methodology, particularly, for images rich in low frequency content. 3. 3.

At the lower sampling ratio of 10%, performance gain over standard methodology is 6.5 dB with ‘Beads’, 9 dB with ‘Lena’, and 11 dB with ‘House’. In fact, standard methodology with standard wavelets almost fails in reconstructing images with any good quality at very low sampling ratios, while the proposed methodology still performs good. For visual clarity, Fig. 15 shows reconstructed images with the existing and the proposed methodology at 10% sampling ratio.

VI-C Experiment-2: Comparison of standard wavelets with matched wavelets in CS-based image reconstruction

In this subsection, we present CS-based image reconstruction results using the PCI sensing matrix and the proposed wavelet decomposition strategy. Table-II shows comparison on reconstruction accuracy in terms of PSNR averaged over 10 iterations obtained with standard wavelets and that obtained with image-matched wavelets designed from compressively sensed images with 5/3 length filters.

We have considered three standard wavelets: orthogonal Daubechies’ wavelets ‘db2’, ‘db4’, and biorthogonal standard 5/3 wavelet. While orthogonal wavelets are widely used in applications, biorthogonal wavelets take care of boundary effects and provide better compression results as compared to orthogonal wavelet [36], [37]. This is to note that our image-matched wavelets are biorthogonal by design.

Discussion: From Table-II, we observe much better reconstruction accuracy with matched wavelets compared to standard wavelets in terms of average PSNR over all the sampling ratio ranging from 90% to 10%. Orthogonal db4 wavelet performs better than other standard wavelets but its performance is still inferior to our matched wavelets. At the sampling ratio of 90%, matched wavelet provides an improvement upto 3 dB, while we observe an improvement upto 11.5dB at 10% sampling ratio.

Further, we observe that while standard wavelets fail almost completely at lower sampling ratio of 10% with a reconstruction PSNR of approximately 10dB only, matched wavelets are able to reconstruct images with a PSNR above 20dB.

For better visual clarity, we have also shown images reconstructed from compressively sensed images at 10% sampling ratio using the PCI matrix with matched 5/3 wavelet and standard db4 wavelet for all ten images in Fig. 16. From the figure, it can be clearly noticed that the existing wavelet db4 is not able to reconstruct full images whereas matched wavelets provide good reconstruction quality.

VII Conclusion

In this paper, we have proposed a joint framework wherein image-matched wavelets have been designed from compressively sensed images and later, used for reconstruction or recovery of the full image. We have also proposed to use a partial canonical identity sensing matrix for CS-based reconstruction of images that performs much faster compared to the existing Gaussian or Bernoulli matrices and hence, is suited for time-bound reconstruction based applications. Although there is a slight degradation in performance with the proposed sensing matrix but that is easily covered up by the matched wavelet design. We have also provided a new multi-level L-Pyramid wavelet decomposition strategy that works much more efficiently compared to the standard wavelet decomposition method. Overall, the proposed work with different sensing matrix, new wavelet decomposition strategy, and image-matched wavelets provide much better results in CS-based image reconstruction compared to the existing practices.