Optimization Methods of Compressively Sensed Image Reconstruction Based on Single-Pixel Imaging

: According to the theory of compressive sensing, a single-pixel imaging system was built in our laboratory, and imaging scenes are successfully reconstructed by single-pixel imaging, but the quality of reconstructed images in traditional methods cannot meet the demands of further engineering applications. In order to improve the imaging accuracy of our single-pixel camera, some optimization methods of key technologies in compressive sensing are proposed in this paper. First, in terms of sparse signal decomposition, based on traditional discrete wavelet transform and the characteristics of coe ﬃ cients distribution in wavelet domain, a constraint condition of the exponential decay is proposed and a corresponding constraint matrix is designed to optimize the original wavelet decomposition basis. Second, for the construction of deterministic binary sensing matrices in the single-pixel camera, on the basis of a Gram matrix, a new algorithm model and a new method of initializing a compressed sensing measurement matrix are proposed to optimize the traditional binary sensing matrices via mutual coherence minimization. The gradient projection-based algorithm is used to solve the new mathematical model and train deterministic binary sensing measurement matrices with better performance. Third, the proposed optimization methods are applied to our single-pixel imaging system for optimizing the existing imaging methods. Compared with the conventional methods of single-pixel imaging, the accuracy of image reconstruction and the quality of single-pixel imaging have been signiﬁcantly improved by our methods. The superior performance of our proposed methods has been fully tested and the e ﬀ ectiveness has also been demonstrated by numerical simulation experiments and practical imaging experiments.


Introduction
According to information theory, one can know that the traditional signal sampling and signal recovery must meet the Shannon-Nyquist theorem [1], which states that to recover the original signal without distortion from a sampled discrete signal, it is required that the signal sampling frequency must be at least twice the bandwidth of the original signal.However, with the rapid development of information technology, the efficient transmission and effective processing of massive amounts of data have put more and more increasing requirements on signals' sampling compression rate (SCR), processing speed and signals' recoverability [2][3][4].In 2006, compressive sensing (compressed sensing, CS) theory was proposed by Donoho et al. [1,[4][5][6].They invented an effective way of economically converting analog signals into digital forms to compress the original signals.Compressive sensing is based on the premise that the original signal can be recovered and be sampled with fewer measurements than Nyquist sampling, and a compressed representation of the original signal is directly acquired by a compressed sensing measurement matrix without the intermediate process of Nyquist sampling [3][4][5][6][7].The combination of signal sampling and signal compression is thus completed in one measurement process, which saves transmission and storage costs and reduces energy consumption and computational complexity [1,6,8,9].Finally, the original signal is reconstructed from the linear measurement data through corresponding reconstruction algorithms [1,6].Thus, compressive sensing greatly reduces the data amount of signal processing and breaks the sampling frequency limitation of the Nyquist sampling theorem.
Because of the excellent performance of CS in the field of signal processing and image processing, compressive sensing has been widely applied in many fields since it was first proposed [7].Some typical application research areas include radar imaging [10][11][12], images' recovery and wireless image sensor networks [13,14], applications in classification problems [15,16], biomedical signal processing [8,[17][18][19], such as nuclear magnetic resonance imaging (MRI) [20,21], electrocardiographic (ECG) signal processing [22,23] and single-pixel imaging [24][25][26][27][28].For the single-pixel imaging, in 2008, Rice University invented a single-pixel imaging system (single-pixel camera) based on CS theory [25], which pioneered the application field of CS-based single-pixel imaging.In their research, a photodiode was used as the optical signal detector to measure the encoded information of an imaging scene, and they achieved CS-based single-pixel imaging for the first time.After Rice University, other researchers have contributed some works based on single-pixel imaging [26][27][28][29].In [26], Sun et al., improved the quality of single-pixel imaging using digital microscanning.In [27], degradation factors added into the algorithm were used to avoid image blurring and achieve single exposure super-resolution compressive imaging.In [28], the authors used a dynamic single-pixel CS system to acquire super-resolution images.In [29], the performance of an imaging approach using two detectors simultaneously was evaluated for real-time video.
Although a lot of works have contributed to optimize single-pixel imaging in image reconstruction algorithms and some important technologies, there still exist some problems in CS-based single-pixel imaging.On the one hand, for the image sparsifying basis, in order to be able to reconstruct images using compressive sensing, the imaging scene needs to be represented as a sparse combination of elements.Hence, if the scene is not sparse, it must first undergo some sparsifying transform, and discrete wavelet transform (DWT) is usually used as the signal decomposition basis [30,31].However, the sparse representation performance of the traditional wavelet transform is not good enough, and reconstruction of higher precision images with a lower sampling compression rate (SCR) requires high enough sparsity of the original signals.On the other hand, in CS-based single-pixel imaging system, a digital mirror device (DMD) and a spatial light modulator (SLM) are usually used to encode and modulate spatial optical signals.However, due to the hardware limitations of DMD or SLM, the sensing matrix of CS must only be deterministic binary matrix in single-pixel imaging system, so deterministic binary matrices need to be designed and constructed [32][33][34][35].In order to facilitate the hardware implementation, the traditional binary sensing matrices such as 0-1 binary sparse random matrix (with the elements of +1 and 0) and part Hadamard random binary matrix (with the elements of +1 and −1) are usually fixed as the deterministic binary matrices in single-pixel imaging system.However, the structure and performance of traditional binary sensing matrices are usually not optimal, which extremely limits the improvement space of expanding compressed sensing measurement matrix and further enhancing imaging quality in single-pixel camera.
Therefore, on the basis of traditional methods, it is very necessary to further optimize the signals sparse decomposition and the binary sensing matrix in compressed sensing, and to that end related work is reported in this paper.First, based on the traditional DWT, the original wavelet decomposition basis was optimized by our proposed constraint condition, which enhancing its sparse representation performance.Then, we propose a new algorithm model and a new initialization method to construct more optimized binary sensing matrices than traditional binary sensing matrices by gradient projection (GP) strategy.Finally, the above proposed key technologies were applied to our single-pixel imaging system and effectively improved the image reconstruction accuracy and imaging quality.The remaining contents of this paper are organized as follows: Section 2 introduces the basic theoretical background and analyzes some related works as the foundation of our contributions; in Section 3, our single-pixel imaging system is built, we design a new optimization method based on the traditional wavelet decomposition basis and propose a new compressed sensing measurement matrix optimization algorithm; in Section 4, we chose appropriate parameters that perform relatively well by the simulation experiments, and test the performance and demonstrate the effectiveness of our methods by practical imaging experiments using our single-pixel camera; then our proposed methods are analyzed and the experimental results discussed in Section 5. Finally, Section 6 summarizes our work and this paper is concluded.

Related Works and Theory Formulation
In this section, we focus on related theory works and the technology background that are the basis of our proposed methods.We present these contents that are closely related to our work on the three parts: compressive sensing, wavelet decomposition analysis and sensing matrix.

Compressive Sensing
Compressive sensing-based signal processing includes two major processes: signal measurement and signal reconstruction.First, the original signal is denoted as x, x = [x 1 , x 2 , . . ., x n ] T , x ∈ C n , [] T is the transpose operator.Considering the signal x is measured in a noisy environment, the measured value is expressed as the vector y, y = [y 1 , y 2 , . . ., y m ] T , y ∈ C m , m is the number of measurements for sampling the original signal.So: where Φ is the sensing matrix [34], Φ ∈ C m×n (m n), reducing the number of measurements compared with the vector x. e is the Gaussian white noise, e ∈ C m .The dimensions of y are far less than the dimensions of x, so the original signal is compressed by reducing the signal's dimensions without the intermediate stage of Nyquist sampling.Therefore, after compressed sensing on the original signal x, it is easier to transmit and store the data y than x.y is the data that needs to be processed for restoring the original scene image x through the computer in CS-based single-pixel imaging system.Since the high-dimensional data x is reconstructed from the low-dimensional data y, the underdetermined equations have to be solved.There are more unknown variables than equations in the process of compressed sensing, thus making it underdetermined (n m), and it is very difficult to solve the underdetermined equations.In order to restore the original signal x by solving the underdetermined equations, it is usually required that the original signal is compressible (sparse) or can be represented sparsely [3,5,36,37].
If x is sparse, the sparsity level of is usually expressed using the "L0 norm" x 0 .x 0 denotes the number of non-zero elements in.So x 0 is added as the regularization penalty term in Equation ( 1) via L0 norm minimization, Equation ( 1) is updated as: where ε is a non-negative upper bound that denotes the noise level to limit the error between the reconstructed signal and the original signal.If x is non-sparse, it is usually necessary to use an orthogonal decomposition basis or sparsifying dictionary for representing the original signal x sparsely [3] as follows: where s is the sparse coefficients that represent x by the sparsifying decomposition basis So in order to acquire the sparse representation s, Equation ( 2) is transformed as: where ŝ is the estimation of s, Θ = ΦΨ is called the equivalent dictionary [38], reducing the number of the basic decomposition dictionary Ψ. Equation ( 4) is the L0 norm minimization in signal sparse decomposition [12], so we can get the estimation x of the original signal from x = Ψŝ through CS-based reconstruction algorithms.
For CS-based reconstruction algorithms, the commonly used algorithms mainly include greedy pursuit algorithms [39][40][41][42], minimum norm optimization algorithms [43][44][45] and iterative threshold algorithms [46,47].Among these algorithms, orthogonal matching pursuit (OMP) algorithm is one of the greedy pursuit algorithms and has stable performance in reconstructing signals because of its simplicity and efficiency [41,42], so the OMP algorithm is usually used quite effectively to get an approximated solution in compressed sensing.

Wavelet Analysis
We use a wavelet transform as the decomposition basis to represent the original signal sparsely in this paper.The wavelet transform can perform the localized analysis on time (space) frequency by scale shifting and time (space) shifting of a mother wavelet.Thus, the information in signal's time domain (space domain) and signal's frequency domain can be analyzed simultaneously at multiple scales.Wavelet transform performs the property of time (space) subdivision at high frequency and frequency subdivision at low frequency [48].In continuous time domain t, the signal x(t) can be expressed by the wavelet transform as: where t denotes time, W a,τ (t) is a mother wavelet, a and τ are the parameters of scale shifting and time shifting respectively [49].Approximate coefficients and detailed coefficients in wavelet domain are generated by a low-pass wavelet filter and a high-pass wavelet filter, respectively.For continuous wavelets, scale a, time t, and time shifting τ are all continuous.If the wavelet transform is used to process images or digital signals, Equation (5) should be discretized as discrete wavelet transform (DWT), and our work of optimizing wavelet decomposition basis in this paper is exactly based on DWT.

Compressed Sensing Measurement Matrix
Compressed sensing measurement matrix is an important intermediate bridge between x and y in CS-based signal processing to reduce the measurements number of the vector x, and it is directly related to the efficiency of signal sampling and the quality of signal reconstruction.It has been proved that sensing matrix should satisfy some necessary conditions to accurately reconstruct the original signal with high probability [4,5,36,37].Among these conditions, the most important one is that sensing matrix Φ should satisfy the restricted isometry property (RIP) [50] as follows: where δ K is the restricted isometry constant (RIC), δ K ∈ (0, 1), K is the sparsity level of x with K non-zero elements.RIP condition is a necessary condition to ensure that the original signal can be recovered with high probability.However, it is very complicated to verify whether the sensing matrix satisfies the RIP and also difficult to compute δ K from the known sensing matrix Φ [51].Thus, finding a simple and efficient alternative condition that can replace RIP is very important to construct a sensing matrix of good performance.In [37], it has been proved that a sensing matrix would meet RIP with large probability if the sensing matrix Φ and sparsifying basis Ψ (or orthogonal dictionary) were as incoherent (orthogonal) as possible, and the recovered signal would have higher reconstruction precision [52].That is to say, if Φ and Ψ are incoherent or uncorrelated, sampling data from the original signal can add more new information that has not already represented by the known basis Ψ.In other words, in the sampling process of CS, the incoherence between Φ and Ψ can ensure that the new information that has not be represented by the known decomposition matrix Ψ can be sampled.Some classic random compressed sensing measurement matrices and binary sensing matrices satisfy the RIP condition well and can guarantee the CS-based signals reconstruction with high probability [50,[53][54][55][56][57].These commonly used classical compressed sensing matrices include Gaussian random matrix, Bernoulli random matrix, Fourier random matrix, 0-1 sparse binary matrix, part Hadamard binary matrix and so on.Among them, 0-1 sparse binary matrix and part Hadamard binary matrix are usually applied and are fixed as the conventional deterministic binary sensing matrices in single-pixel imaging system.
According to the theoretical analysis in [58,59], the mutual incoherence between sensing matrix and sparsifying basis means that the corresponding Gram matrix is as close as possible to an identity matrix, and the Gram matrix is defined as: where Θ H is the conjugate transpose operator of Θ, Θ = ΦΨ.Thus, minimizing the mutual coherence is equivalent to shrinking the absolute off-diagonal elements of the corresponding Gram matrix.The more closely Gram matrix approximates the identity matrix, the less the coherence between Φ and Ψ is.In order to obtain the compressed sensing matrix of better performance, the following optimization problem is denoted as: where I is an identity matrix of size n × n.Therefore, to construct an optimized sensing matrix, Equation (8) has to be solved.

Methods
In this section, we introduce our proposed methods that improve the quality of CS-based image reconstruction in our single-pixel imaging system, optimizing the original wavelet decomposition basis and the traditional binary sensing matrices.

Imaging Method in Our Single-Pixel Camera
We have built a CS-based single-pixel imaging system (single-pixel camera) in our laboratory and reconstructed two-dimensional imaging scenes by the methods of CS-based signal processing.Figure 1a,b respectively show the schematic diagram and a physical view of our single-pixel imaging system.
As shown in Figure 1, our single-pixel imaging system consists of: the imaging lens 1, the Digital Mirror Device (DMD, V-7001, ViAlUX, Chemnitz, Germany, operating at 22.7 kHz/s), the mirror, the light collection lens 2, the Silicon Photomultiplier (SPM) Detector (C13369, HAMAMASTSU, Shizuoka, Japan, 320 nm to 900 nm, photoelectric sensitivity at 1.0 × 10 9 V/W), the Analog Digital Converter (ADC, PCI50641, TDEC, Chengdu, China, sampling frequency at 50 MHz) and the computer.The lens of our single-pixel camera are designed for distant scenes, and the imaging objects are distant buildings.The SPM detector is very sensitive to optical signals with detection ability of photon level, which can detect valid signals quickly and accurately and it is beneficial to improve the signal-to-noise ratio of signal detection.
Our single-pixel camera's specific imaging process is: firstly, the lens 1 collects the light of the distant scene (equivalent to the original signal x) and projects the light onto DMD, which means that the image plane of the imaging scene falls exactly on the DMD surface; then the computer encodes the row masks of sensing matrix Φ to the DMD scanning sequence, and DMD performs spatial light modulation on the image plane, which means that DMD picks the necessary pixel light and filters out the unnecessary pixel light (equivalent to the signal sensing measurement process of Φ in CS) as shown in Figure 2  It should be noted that DMD is a spatial light modulator (SLM) with digital micromirror arrays and facilitates the hardware implementation and engineering application of the sensing matrix in the CS-based imaging system.The DMD sequence indicates the elements distribution status in binary sensing matrix.Due to the hardware limitation, as shown in Figure 2, every micromirror unit of DMD arrays only has binary states with the two flipping angles of +12 degrees and −12 degrees, and the two flipping angles are corresponding to the "1" and "0" (or "−1") binary states, respectively, in the sensing matrix, so only deterministic binary sensing matrices can be used in the single-pixel camera, such as the 0-1 binary sparse matrix with the elements of only "0" and "1" and the part Hadamard matrix with the elements of only "+1" and "−1".In the 0-1 binary sensing matrix, the element "0" means that the corresponding pixel micromirror of DMD is off (−12 degrees) removing the unnecessary pixel light, and the element "1" indicates the corresponding pixel of DMD is on (+12 degrees) sampling the necessary pixel light.The sampled signal value Y is detected by the photodetector, as shown in Figure 1a.However, for the part Hadamard matrix, it is very hard for DMD to physically realize the "−1" state of the Hadamard matrix directly.Thus, in order to realize the distribution of the Hadamard matrix in DMD, our approach is: before coding the Hadamard matrix, first we set all the effective micromirror pixel units of DMD to on ("+1"), and the photodetector acquires the corresponding measurement value as Y1.Then the elements "−1" and "+1" of the Hadamard matrix code the pixel micromirror units of DMD to off and on respectively (at this time, the element "−1" is equivalent to the element "0" of 0-1 binary matrix, which sets pixel micromirror unit to off), detecting the corresponding measurement value as Y2.The "+1" states in Y2 can correctly and truly represent the "+1" elements of the Hadamard matrix, while the "−1" states in Y2 represent the off states (elements "0") of micromirror units and the photodetector can't acquire the corresponding signals of "−1".There are the same and common "+1" states of micromirror units in Y1 and Y2.Thus, in the Hadamard matrix, for the "+1" elements, 2Y2-Y1 means 2 × +1 − 1 × +1 = 1 × +1 that is exactly the "+1" elements of the Hadamard matrix, while for the "−1" elements, 2Y2-Y1 means 2 × 0 − 1 × +1 = −1 × +1 = −1 that can actually achieve the representation of "−1" elements of the Hadamard matrix.Consequently, for Hadamard matrix, the measurement value is obtained as Y = 2Y2 − Y1.
We measure m times in this same way as Y and obtain the measurement vector y, y = [y 1 , y 2 , . . ., y m ] T .In order to eliminate the interference of system noise and dark noise as much as possible, our adopted method is to set all the effective micromirror pixel units of DMD to "0", which means that all the micromirror units are turned off and SPM detects the background noise as Y0.Thus, the effective measurement value at each measurement is Y-Y0 that has most of the ambient noise and dark noise removed.In addition, to avoid the influence of ambient light noise, our whole imaging system was blackened and sealed with a cover.In our single-pixel imaging system, the projection method of column-by-column scanning is adopted to acquire the imaging information, and the imaging resolution of image reconstruction is 256 × 256.Because of the stable and efficient operation of OMP algorithm [41,42] in reconstructing small data (column signals of length 256 are small data), we choose OMP algorithm as the imaging algorithm to reconstruct imaging scenes in our single-pixel camera.
During the process of SPM signal detection, due to the fact the light intensity fluctuates over time, the detected electric signal intensity at each imaging session may also fluctuate with the optical signal.In order to as much as possible eliminate the dynamic deviation of grayscale at each imaging and make the grayscale distribution of reconstructed images more balanced and more stable, we normalize the grayscale range of every imaging picture and equalize the reconstructed images into the grayscale range of 0 to 255.The grayscale equalization equation is expressed as: the final reconstructed images are obtained by Equation (9).

DWT-Based Decomposition Optimization
The sparsity of signals is directly related to the signal sampling compression rate (SCR) and the signal reconstruction quality.This wouldn't mean that there must be K non-zero elements in sparse signals and all the remaining elements are 0 [61].If all the coefficients of this signal were arranged from large to small in the sparse domain and the arranged coefficients were approximately an exponential decay distribution of approaching to zero, the signal would be also called an approximately sparse signal or a compressible signal [61].DWT is usually used to compress signals, the coefficients distribution of non-sparse signals in DWT domain is very close to the exponential decay distribution, as shown in Figure 3. Figure 3a is a random non-sparse column signal of length 256, and the non-sparse column signal is transformed to the sparse coefficients distribution by DWT (the wavelet decomposition basis), as shown in the blue curve of Figure 3b.One can see that the wavelet transform can make the non-sparse signal become more sparse from Figure 3a to Figure 3b.As shown in the blue and red curves respectively of Figure 3b, the wavelet coefficients are gradually decreasing and are very close to the exponential decay distribution.Because the coefficients in DWT domain show a sparse distribution of approximately exponential decay, we use DWT to sparsify the non-sparse signal.In order to make wavelet coefficients more sparse and closer to a true exponential decay distribution, we design a diagonal constraint matrix with its diagonal elements exponentially decreasing as the constraint condition, as shown in Figure 5.This constraint matrix is denoted as W and the diagonal elements of W are (p) i−1 in which p is a parameter that can adjust the constraint degree, 0 < p < 1, i is the i-th diagonal element and n is the length of the original signal.Therefore, the constraint matrix can be expressed by the following equation as: In Equation (10), the first term is 1, and all the remaining elements are positive and less than 1 with the exponential decay distribution, so this constraint matrix can perform a constraint role of forcing the decaying wavelet coefficients to distribute more exponentially and more sparsely.
This proposed constraint condition is added into DWT, so Equation ( 4) is updated as: The estimation ŝ of wavelet coefficients is obtained by solving the problem (11) with quadratically constrained linear program (QCLP), and the estimation of the reconstructed original signal is gotten by Equation ( 12) as follows: x = ΨW −1 ŝ.
Because our imaging scenes are non-sparse optical signals, this method of proposing the constraint matrix can theoretically make our imaging scene be represented more sparsely and improve the imaging quality, and the effectiveness of this method will be demonstrated by experiments in Section 4.

Binary Sensing Matrix Optimization
According to the introduction in Section 2, it is required that the elements of sensing matrix must be binary in single-pixel camera, and the traditional binary sensing matrices usually used in single-pixel camera are 0-1 binary sparse matrix and part Hadamard matrix.Thus, in this subsection, our proposed sensing matrix optimization method is based on the traditional random binary sensing matrices.We focus on seeking a more optimized binary sensing matrix that can minimize the mutual coherence between sensing matrix and sparsifying basis and further improve the accuracy of signal reconstruction.
Since there may be both positive and negative elements in binary measurement matrices, such as the Hadamard matrix, we divide Φ into the positive part U and the negative part V, so Φ is denoted as: where u i,j = φ i,j + and v i,j = −φ i,j + , i = 1, 2, . . ., m and j = 1, 2, . . ., n. φ i,j + is the positive-part operator, such as φ i,j + = max 0, φ i,j .u i,j , v i,j and φ i,j are the elements in the matrices U, V and Φ respectively.
From the model (8), one can know that, due to the coupling of sensing matrix and sparsifying matrix, it is very complicated to directly solve Equation (8) and also difficult to guarantee the construction accuracy of the sensing matrix.In order to construct a more optimized sensing matrix, we propose a new method of solving Equation ( 8) by adding a constraint term of approximate equivalent dictionary (matrix) Θ ≈ ΦΨ, for ensuring higher accuracy in constructing sensing matrix.Therefore, the following optimization problem is further obtained by updating Equation ( 8): where ξ is the upper bound of the approximate sensing matrix term and Φ is the estimated solution of Φ. Equation ( 14) is a convex constrained optimization problem.We convert the constrained optimization problem to an unconstrained optimization problem by lagrangian method and propose a new mathematical model based on Equation ( 14), as expressed in Equation ( 15): so that the constraint becomes a penalty.For a proper choice of τ, the two problems ( 14) and ( 15) are equivalent, where τ is the weighting factor chosen empirically for balancing the two items of Equation (15).F( Φ) is the loss function (or objective function) that needs to be solved via seeking the minimum, so Equation ( 15) is an unconstrained convex optimization problem.
The term Θ − ΦΨ , but also is the accuracy guarantee of Φ by solving the penalty term Θ − ΦΨ 2 F .For the new model (15) with the specific two terms, the descent direction of F( Φ) is not only the gradient direction but also the direction of gradient projection, which better ensures that our method is able to approximate the optimal solution of Φ.Thus, the iterative gradient projection (GP) based strategy is used to seek the minimum value of the loss function F( Φ).The gradient ∇F( Φ) of F( Φ) is expressed as: The sensing matrix Φ is updated by continuous iterations, and Φ (k) is the intermediate estimated solution of Φ during the iterations.That is, from the k-th iteration to the k+1-th iteration, Φ (k) is updated to Φ (k+1) .The following equation is defined with a positive scalar parameter α (k) , α (k) > 0, as: At each iteration of Φ (k) , the algorithm searches along the negative gradient direction −∇F(Φ (k) ), projecting onto the non-negative orthant and taking a backtracking line search until sufficient decrease is sought out in F(Φ (k) ).
Then, the second scalar parameter λ (k) is defined, Using λ (k) and α (k) reduces the possibility that F( Φ) may increase at some time during iterations and enhances the efficiency to search the global optimal solution of Φ in our algorithm.

Algorithm Operation
The previous subsection is the basic idea of our proposed algorithm, the specific running process and details of our algorithm is described in this subsection for solving the loss Equation (15).The gradient projection based specific steps of our algorithm are presented as follows.
In step 1, our proposed initialization method of Φ (0) and Θ (0) will be described in details in Section 3.5.In the above algorithm, the sensing matrix Φ is trained and optimized by making full use of the known information of sparsifying decomposition matrix Ψ.It is worth noting that, after updating Φ (k) at each iteration, every column of Φ (k) is regularized to facilitate algorithm operation and make experimental results more balanced, as expressed in Equation ( 21).Because 0-1 binary random matrix and part Hadamard binary matrix are used in our single-pixel imaging system, Φ (k) needs to be set to binary matrices after regularizing Φ (k) at each iteration, as denoted in Equations ( 22) and (23).After the iteration stops, the last optimized binary sensing matrix Φ (k) is fixed as the deterministic binary compressed sensing matrix that is applied to the single-pixel camera.From the above theoretical analysis, one can know that the new optimized sensing matrix Φ obtained by our proposed method can be self-adaptive to the known wavelet decomposition matrix Ψ, and would have lower mutual coherence with Ψ and higher construction accuracy in theory.

New Initialization
A wise initialization of sensing matrix could be of great worth, which directly determines the iteration efficiency of our algorithm and affects the line backtracking search direction of Φ (k) .Thus, we proposed a new method of initializing Φ based on the characteristics of the loss Equation (15).Firstly, Θ is initialized according to the pseudo inverse Θ H ⊥ of Θ H at the assumption that Gram matrix is equal to an identity matrix, ), the equation is followed as: and: So the initial Θ (0) is expressed as: Thus, the initial sensing matrix Φ (0) is gotten from Θ (0) = Φ (0) Ψ, and Φ (0) is denoted as: The Φ of Equation ( 30) is the random binary sensing matrix that needs to be initialized randomly at first.For example, in our single-pixel imaging system, Φ can be initialized as 0-1 binary random matrix and part Hadamard random matrix.
The above method of initializing Φ takes full advantages of the two approximate terms in the loss Equation ( 15) assuming the two ideal situations of Θ = ΦΨ and Θ H Θ = I and simultaneously combines the characteristics of the two terms, which conforms to the properties of solving the objective Equation (15).Consequently, Θ (0) obtained by our initialization and the fixed Ψ are input to the step 1 of our algorithm as known variables for training more optimized sensing matrix Φ.After applying our initialization, the risk of local minima and instability the solution of Φ may fall into is ameliorated by the fact that the initialization of Φ (0) by Equation ( 30) has already been in the neighborhood of a desirable attraction basin of the global optimum or a local optimum.Thus, compared with other initialization methods such as the direct random binary initialization on Φ, our initialization method can train more optimized sensing matrix, theoretically approximating the global optimum that is almost equivalent to the ideal solution of Θ = ΦΨ and (ΦΨ) T ΦΨ = I simultaneously.In other words, this new initialization method facilitates the loss function F(Φ) to approximate the global optimal solution with higher efficiency and greater probability.

Experiments
In this section, we describe experiments conducted to test the performance of our optimized wavelet decomposition basis and optimized binary sensing matrix for compressively sensed image reconstruction, compared with the traditional wavelet decomposition basis and the conventional binary measurement matrices that are usually used in single-pixel imaging respectively.Finally our proposed methods were applied to improve the traditional imaging methods of our single-pixel camera by conducting practical imaging experiments.
The imaging resolution of our single-pixel camera is 256 × 256.In order to better balance the sampling speed and the image reconstruction accuracy, we set the usually-used sampling rate to 0.25 in our imaging system with m = 64 and n = 256.Because our single-pixel camera used the imaging method of column-by-column scanning via DMD, we measured the each column of images or imaging scenes separately for 64 times (64 rows of the sensing matrix) by the binary sensing matrix of size 64 × 256, and then reconstructed images' columns of size 256 × 1 one by one, finally added the 256 reconstructed columns back to the full image of size 256 × 256.Our DMD worked at 22.7 kHz/s, so the imaging time of column-by-column was 1×64×256 22700 ≈ 0.72 seconds.The imaging time of scanning was very short, which could reconstruct a full 2D image quickly and well.
Since our single-pixel camera experiment system was designed to take images of distant scenes such as buildings, the imaging field of view is some limited.Before the practical imaging experiments, in order to determine the proper empirical parameters in our algorithm, we conducted the simulation experiments in reconstructing images that are similar to our real imaging scenes.The images used in our simulation experiments are images Building and House that are similar to the actual imaging scenes, as shown in Figure 6, and the two images can well simulate the real imaging scenes.In the experiments, 0-1 binary sparse sensing matrix and part Hadamard sensing matrix are used, and the image reconstruction algorithm is OMP algorithm in which the iteration number is m, as mentioned in Sections 2.1 and 3.1.The 0-1 sparse binary sensing matrix was produced by the function "SparseRandomMtx(m, n, m/2)" of Matlab, in which m denotes the row number, n denotes the column number and m/2 denotes the number of "1" in every column, so the elements "0" and "1" are equally weighted with the same 50% respectively in the 0-1 binary matrix.The Hadamard matrix sequence was generated by the function "hadamard(n)" of Matlab, and the Hadamard matrix H a is a square matrix and  In the simulation experiments, we used the peak signal to noise ratio (PSNR) and the structural similarity (SSIM) to evaluate the accuracy and quality of image reconstruction, PSNR and SSIM are defined as: where x is the original image, and x is the reconstructed image.H and W are the height and the width of the image respectively.µ x and µ x are respectively the mean values of x and x. σ 2 x , σ 2 x and σ x x are the variances and covariance of x and x respectively.c 1 and c 2 are the constants, L is the range of pixel greyscale, k 1 = 0.01, k 2 = 0.03, L = 255.

Setting in DWT Optimization
In our proposed optimization method of wavelet decomposition basis, one can know that the size of p in Equation ( 10) directly determines the restraint degree of the restraint matrix on wavelet coefficients.The distribution of wavelet coefficients is regulated by adjusting the value of p.If p is too large, the affect of constraint matrix constraining the wavelet coefficients is weakened, and the sparsity of wavelet coefficients cannot be effectively improved, which is not good for improving the quality of image reconstruction.However, if p is too small, the constraint degree of constraint matrix in DWT domain is higher, which would increase the computational difficulty of our algorithm operation and cause information loss, this is also not helpful to reconstruct images.Therefore, choosing a proper p value that performs well is very important for the optimization of wavelet decomposition basis, and directly affects the results of image reconstruction.
In order to choose a relatively optimal p value, when m = 64 and n = 256, images House and Building were reconstructed separately with different p values under the 0-1 binary sensing matrix and the part Hadamard sensing matrix, then we determined an appropriate parameter p according to the simulation experimental results.We obtained the PSNR results of images reconstruction with different p values by our method optimizing DWT under 0-1 binary matrix and part Hadamard matrix, respectively.Since the 0-1 binary matrix and the part Hadamard matrix were generated randomly, in order to avoid experiments' randomness and make experimental results more stable and reliable, the averages of five experiments' results were calculated and the average experimental results are shown in Figure 7.As shown in Figure 7, whether for the 0-1 binary sensing matrix or the part Hadamard sensing matrix, the PSNR of reconstructed images increases first and then decreases as p increases, and the optimal range of p value is about [0.995, 0.997].Thus, after giving full consideration to the balance of images Building and House and the balance of two types of binary sensing matrices, p value is set to 0.996 in our actual imaging algorithm.
As shown in Figure 7, after adding the constraint condition of exponential decay to the wavelet decomposition basis, the PSNRs of reconstructed images are effectively improved.For both of the 0-1 binary sensing matrix and the part Hadamard sensing matrix, the PSNRs of reconstructed image Building are improved by about 1.5 dB, and the PSNRs of reconstructed image House are improved by about 2 dB.When p is set to 0.996, the optimized DWT improves PSNRs more significantly by about 2 to 2.5 dB.The experiments of Figure 7 also implicitly demonstrate the effectiveness of our method in improving the quality of compressively sensed images reconstruction.

Experiments for DWT Optimization
According to the good improvement results of above simulation experiments in reconstructing images Building and House, parameter p was set to 0.996 and our method of optimizing DWT was further applied to our single-pixel imaging system for practical imaging experiments.0-1 binary sensing matrix and part Hadamard binary sensing matrix are used as the masks to drive and encode the DMD sequences in our imaging system, as shown in Figure 1a.We used the single-pixel camera to take images of two scenes of the distant buildings at the sampling rate of 0.25, and reconstructed images of the two imaging scenes are shown in Figures 8a-d and 9a-d.The close-up regions marked with a red solid rectangle were enlarged, as shown in Figures 8a and 9a.
According to the practical imaging results of Figures 8a-d and 9a-d, we can give the qualitative evaluation from the visual quality.Our wavelet optimization method effectively improves the imaging quality.For both of the 0-1 binary sensing matrix and the part Hadamard sensing matrix, the reconstructed images by our method have better visual quality with clearer image details, lower noise and better contrast, and the window contours are more obvious and sharper in Figure 8c,d and Figure 9c,d.
Because there are not original images as the evaluation benchmarks for the two imaging scenes, in order to quantitatively analyze and evaluate the imaging quality, we define the signal to noise ratio (SNR) to quantitatively calculate the imaging accuracy, as expressed in Equation (33): where I f is the average grayscale intensity of the feature (here calculated from the pixel grayscale highlighted by the solid green rectangle in Figures 8a and 9a), I b is the average grayscale intensity of the background (calculated from the pixel grayscale highlighted by the solid blue rectangle in Figures 8a and 9a), σ f and σ b are the standard deviations of the grayscale intensities in the feature and the background respectively.The imaging results of scene 1 before and after using the proposed optimization methods at the sampling rate of 0.25 ("00" denotes the traditional imaging method before optimization; "01" denotes the imaging method after only optimizing wavelet decomposition basis; "10" denotes the imaging method after only optimizing binary sensing matrices; "11" denotes the imaging method after optimizing both wavelet decomposition basis and binary sensing matrices simultaneously).
Figure 9.The imaging results of scene 2 before and after using the proposed optimization methods at the sampling rate of 0.25 ("00" denotes the traditional imaging method before optimization; "01" denotes the imaging method after only optimizing wavelet decomposition basis; "10" denotes the imaging method after only optimizing binary sensing matrices; "11" denotes the imaging method after optimizing both wavelet decomposition basis and binary sensing matrices simultaneously).
According to the evaluation index of Equation ( 33), we obtained the statistic SNRs table of Figures 8a-d and 9a-d, as shown in Table 1.In order to express each imaging method conveniently and concisely, we used the "00" "01" "10" "11" to denote the four different imaging methods.From Table 1, one can know that, whether for the 0-1 binary matrix or the part Hadamard matrix, the SNRs of reconstructed images for both two imaging scenes have been very evidently improved after our method optimizing wavelet decomposition basis, as shown by the bold in Table 1.Therefore, the data results of Table 1 quantitatively demonstrate the superiority and effectiveness of our proposed DWT optimization method.
Table 1.SNRs (dB) of the imaging results of our single-pixel camera before and after optimizing wavelet decomposition basis ("00" denotes the traditional imaging method before optimization; "01" denotes the new imaging method after optimizing wavelet decomposition basis).

Sensing Matrix
Imaging

Experiments for Binary Sensing Matrix Optimization
In order to test the performance of our proposed algorithm of optimizing traditional binary sensing matrix and get the empirical values of parameters in Section 3.4, we also conducted simulation experiments of reconstructing images Building and House by using our optimized binary measurement matrices.After testing some values, we found empirically that these parameters perform relatively well that ξ = 0.01, α min = 1 × 10 −3 , α max = 1 × 10 3 , and τ = 0.35 in the loss Equation (15).DWT matrix was used as the known decomposition basis in our experiments.Table 2 is the experiments' statistic results of reconstructing images Building and House before and after optimizing the binary sensing matrices with m = 64 and n = 256.Table 2. PSNRs (dB)/SSIMs of reconstructed images Building and House before and after optimizing binary sensing matrices ("00" denotes the traditional method before optimization; "10" denotes the new method after optimizing binary sensing matrices).

Sensing Matrix
Building House From Table 2, one can know that, after optimizing the traditional 0-1 binary sensing matrix and part Hadamard binary sensing matrix by using our proposed algorithm, the optimized binary measurement matrices enhance the PSNRs of reconstructed images very significantly by about 2 to 3 dB and enhance the SSIMs by about 0.1, and the reconstructed images have higher precision.Then, the new optimized binary measurement matrices were fixed as the deterministic binary sensing matrices that were further applied to our single-pixel camera for practical imaging experiments.The experimental results are shown in Figures 8e,f and 9e,f and the corresponding SNRs are given in Table 3.
Table 3. SNRs (dB) of the imaging results of our single-pixel camera before and after optimizing binary sensing matrices ("00" denotes the traditional imaging method before optimization; "10" denotes the new imaging method after optimizing binary measurement matrices).

Sensing Matrix
Imaging  8a,b and 9a,b of traditional methods, after using the new optimized binary measurement matrices in our imaging system, the imaging quality of two scenes has been evidently improved with less visual artifacts and higher contrast, the reconstructed images have less distortion and noise is evidently suppressed.In Figures 8e,f and 9e,f, the contours become more clear and the edges get sharper, and the details of the locally enlarged regions are also richer.At the same time, the imaging SNRs have been greatly improved as shown in Table 3.These typical experiments have demonstrated well that our method of optimizing the traditional binary sensing matrices can further enhance the accuracy of reconstructed images and favorably improve the imaging quality of our single-pixel camera.

Compressive Experiments
In order to test the comprehensive performance of our methods in optimizing wavelet decomposition basis and binary sensing matrices, both the two optimization methods of wavelet decomposition basis and binary measurement matrices were simultaneously applied to our single-pixel camera at the sampling rate of 0.25, the imaging results are shown in Figures 8g,h and 9g,h and the corresponding SNRs are given in Table 4.
Table 4. SNRs (dB) of the imaging results by optimizing binary sensing matrices and wavelet decomposition basis ("00" denotes the traditional imaging method before optimization; "01" denotes the imaging method after only optimizing wavelet decomposition basis alone; "10" denotes the imaging method after only optimizing binary sensing matrices; "11" denotes the imaging method after optimizing both wavelet decomposition basis and sensing measurement matrices simultaneously).It can be seen from Figures 8g,h and 9g,h that, after simultaneously optimizing the wavelet decomposition basis and binary sensing matrices by our proposed methods, the imaging visual quality of Figures 8g,h and 9g,h is significantly higher than that of Figures 8a,b and 9a,b, and the proposed methods can produce clearer and sharper results and recover the smaller details of the imaging scenes, and the reconstructed images are also closer to the real imaging scenes.Furthermore, compared with Figures 8a,b and 9a,b, our methods have favorably suppressed the column visual artifacts that were produced by the imaging method of column-by-column scanning and our methods play the important roles in reducing the noise of reconstructed images.From Table 4, one can know that the comprehensive performance of our two optimization methods is superior than either one optimization method alone.In addition, the two imaging scenes were reconstructed at six different sampling rates and the statistic SNRs results of the two reconstructed imaging scenes were calculated, and the curves of experimental results are shown in Figure 10.

Sensing Matrix
One can see that from Figure 10, for different compression sampling rates (SCR), our proposed optimization methods can effectively improve the accuracy and quality of images reconstruction.Compared to the traditional imaging method ("00") before optimization, whether the wavelet decomposition basis and the binary sensing matrices are optimized separately or simultaneously, the imaging SNRs and the accuracy of reconstructed images have been enhanced pretty obviously.Furthermore, our proposed methods of optimizing both the wavelet decomposition basis and the binary sensing matrices simultaneously ("11") achieve the best performance compared with the traditional method or either one optimization method alone.These experimental results have fully demonstrated the effectiveness of our proposed methods in optimizing the traditional wavelet decomposition basis and the conventional binary sensing matrices.
Figure 10.The evolution relationship between the imaging SNRs of reconstructed scenes and different compression sampling rates (SCR) for (a) 0-1 binary sensing matrix and (b) part Hadamard binary sensing matrix ("00" denotes the traditional imaging method before optimization; "01" denotes the imaging method after only optimizing wavelet decomposition basis; "10" denotes the imaging method after only optimizing binary sensing matrices; "11" denotes the imaging method after optimizing both wavelet decomposition basis and binary sensing matrices simultaneously).

Analysis and Discussion
These above experiments have well tested the superior performance of our proposed methods in optimizing traditional wavelet decomposition basis and constructing new binary measurement matrices based on the conventional binary sensing matrix.According to the experimental results, one can know that our methods have better performance in extracting information from the original signals and restoring the original signals from the sampled information.Our proposed methods of optimizing the traditional wavelet decomposition basis and the conventional binary sensing matrices can favorably suppress the noise of reconstructed images better compared with before optimization, as shown in Figures 8 and 9.In addition, the following four points are worth being analyzed and discussed:

•
Firstly, because we applied the method of column-by-column scanning to reconstruct the image, there are some column visual artifacts in the reconstructed image, which is the most obvious in Figure 8a Secondly, when the Hadamard matrix is applied in our single-pixel camera, although each signal measurement value is Y = 2Y2 − Y1, it is required that Y1 is measured only once at the beginning when all the effective micromirror pixel units of DMD are set to on.This means that the measurement times of Hadamard matrix with the added measurement Y1 is only once more than that of 0-1 binary matrix.DMD works at 22.7 kHz/s, which makes a full 2D image be reconstructed quickly and well for about 0.72 s.Thus, the calculation method of Y = 2Y2 − Y1 will not extra increase more additional noise.

•
Thirdly, in the process of signal sampling, the amount of information extracted from the original signal increases with the rise of sampling compression rate (SCR), which makes it easier to recover the original signal from more sampled information.This is why the improvement degree of imaging SNR in our single-pixel camera tends to decrease slowly as sampling rate increases, as shown in Figure 10.Thus, we have to consider the trade-off between SCR and recovery accuracy in our imaging experiments, and 0.25 is an appropriate middle sampling rate that can give full consideration to balancing the imaging speed and the imaging accuracy.

•
Fourth, one can generally know that, in image reconstruction, the higher SNR is, the more difficult it is to further improve SNR.In the comprehensive experiments, we think adding the constraint matrix could affect the coherence between sensing matrix and wavelet decomposition basis to some extent, because the sensing matrix and the wavelet decomposition basis may be coupled from each other.Consequently, the simultaneous optimization ("11") of both wavelet decomposition basis and sensing matrix is a little superior than either one optimization method alone ("01" or "10").However, compared with the traditional imaging method ("00") before optimization, the imaging SNRs have significantly been enhanced by our methods of whether separately ("01" or "10") or simultaneously ("11").

Conclusions
Based on the basic theory of compressive sensing, a single-pixel imaging system was built in our laboratory.We proposed the constraint condition of exponential decay to optimize the traditional wavelet decomposition basis, which effectively improved the sparse representation performance of discrete wavelet transform.At the same time, in order to facilitate DMD hardware implementation and engineering application of binary sensing matrix in our single-pixel imaging system, we proposed a new optimization algorithm based on the new mathematical model and the new initialization method to construct more optimized binary measurement matrices of better performance.Our methods not only effectively optimize the signal sparse decomposition and binary sensing matrices and improve the imaging methods of our single-pixel camera, but also greatly enhance image reconstruction accuracy and imaging quality.Our research may enlighten widely exploration in this direction of the field and have certain application value in expanding wider range researches.In the future works, we will make further efforts in the researches of scanning projection mode of DMD, signal sparsifying transform and signal reconstruction algorithms in compressive sensing and single-pixel imaging.
; then the sampled pixel lights are reflected by the mirror and collected by the lens 2 into the photodetector (SPM); the SPM converts optical signals into electric signals, and then ADC converts analog electric signals to digital signals (the digital signals are equivalent to the measured signals y); finally our computer reconstructs the imaging scene (equivalent to x) from the digital signals y through CS reconstruction algorithms, which means that the imaging scene is recovered from the small sampled data by single-pixel imaging.

Figure 1 .
Figure 1.(a) Schematic diagram, (b) experimental set-up of our single-pixel imaging system.

Figure 2 .
Figure 2. DMD modulates the spatial light and and the flipping micromirrors on DMD arrays (this picture has been modified according to our content on the basis of the original one, Figure 3 of [60]).

Figure 3 .
Figure 3. (a) The original random non-sparse column signal, (b) the coefficients distribution in DWT domain of the original column signal.Based on the above analysis, we propose a new constraint condition of the positive exponential decay distribution of which the first term is 1, and this new constraint condition is applied to constrain the existing wavelet coefficients.The curve of this constraint function and the corresponding discrete histogram of this exponential decay are shown in Figure 4, in which the red solid curve is an exponential constraint distribution, and the blue columns are the constraint factors of wavelet coefficients of the corresponding sequence position.The constraint factors of blue columns are all positive real number less than 1 except the first term being 1.In other words, the wavelet coefficients are multiplied by the discrete exponential decay distribution with the first term being 1, and this constraint condition can make the wavelet coefficients distribution get more sparse and closer to the actual and real exponential decay.

Figure 4 .
Figure 4.The exponential decay function and the corresponding discrete histogram.

2 F 2 F 2 F. 2 F
is based on the approximate solution of Θ = ΦΨ and makes it easier to solve Φ by balancing Θ − ΦΨ and ΦΨ T ΦΨ − I Our proposed new model (15) not only can ensure smaller mutual coherence between sensing matrix and decomposition matrix by the cost ΦΨ T ΦΨ − I H a T H a = nI, where [n n]=size(H a ) and I = eye(n,n), and the part Hadamard matrix was constructed by extracting the corresponding rows or columns from Hadamard square matrix.The algorithms operated in the environment of MATLAB 2017b on a standard PC with a 2.20 GHz Intel Core i5 CPU and 4 GB memory, running on the Windows 10 operating system.

Figure 7 .
Figure 7. Relationships between PSNR and p in reconstructing images Building and House for (a) 0-1 binary sensing matrix and (b) part Hadamard sensing matrix by our method optimizing wavelet decomposition basis.

Figure 8 .
Figure8.The imaging results of scene 1 before and after using the proposed optimization methods at the sampling rate of 0.25 ("00" denotes the traditional imaging method before optimization; "01" denotes the imaging method after only optimizing wavelet decomposition basis; "10" denotes the imaging method after only optimizing binary sensing matrices; "11" denotes the imaging method after optimizing both wavelet decomposition basis and binary sensing matrices simultaneously).
,b and Figure 9a,b before optimization.Our optimization methods have effectively reduced the generation of column artifacts and reconstructed the images with more balanced structure and clearer texture, as shown in Figures 8c-h and 9c-h.•