1. Introduction

The idea of compressive imaging (CI) has been used when a large amount of data need to be acquired [1], while most of the collected data will be discarded [2]. The method is also useful when high-resolution sensor arrays are not available due to limitations of manufacturing technology, for example, in LIDAR 3D imaging [3, 4] or terahertz imaging [5]. CI is also suitable for applications where the data acquisition device is expensive, such as low-light-level imaging (L³-imaging) [6–8]. Using CI for L³-imaging takes advantage of better measurement signal-to-noise ratio (SNR). In such a system, the measurements are linear combinations of multiple pixels, consequently the measurement power is much larger than that of a single pixel.

Nevertheless, inherent challenges in CI make it difficult for implementation in many fields for real-world applications. One difficulty is its sensitivity to system calibration [9]. System calibration error can also be amplified by a system’s data processing. Another difficulty is the size of the signal space [9]. This is a challenge for computational memory allocation and computation power. Besides these difficulties, CI system’s speed also limits its applications, especially for real-time imaging. Several groups, including us, have investigated the use of GPU or FPGA to accelerate the data processing [8,10–14]. However, few have considered accelerating the data acquisition. In this paper, we develop a methodology to compute the optimized binary sensing matrices to improve the speed in the system’s measurement collection.

In compressive imaging, sensing matrix is always a crucial factor. Initially, random sensing matrices, such as matrices with elements selected from Bernoulli or Gaussian distribution, or matrices with row vectors randomly selected from a Fourier transformation matrix, have been studied. A CI system using these matrices can successfully reconstruct sparse signals, as they satisfy the restricted isometry property (RIP) [15, 16]. However, as noise increases, the performances degrades fast [17, 18]. Even without noise, given a sparse representation dictionary, searching an optimal sensing matrix is an open question for CI. Therefore, several research groups work on sensing matrix design for different applications [19,20]. A major trend focuses on the mutual coherence of the multiplication of the sensing matrix and the sparse representation dictionary [21–24]. When this coherence is small, a wide set of signals can be reconstructed successfully using the sensing matrix. Different from using mutual coherence, by analyzing the oracle estimator MSE as Gaussian noise present in system measurements, Chen et al. [25, 26] find that a good sensing matrix ought to be close to a Parseval tight frame. Based on the work in literature [21,26], Yang et al. [27] further optimize the sensing matrix combining with a matrix decomposition, thus to have the optimal reconstruction and the optimal mutual coherence performance simultaneously. Different from the works discussed above, in which the sparse representation dictionary is predefined, Duarte-Carvajalino and Sapiro [22] design the dictionary and the sensing matrix jointed from a set of training images to obtain better reconstructions. Besides the mutual coherence, several other groups use information-theoretic framework to design a sensing matrix [28–32].

Note that the design methods discussed above are for gray value sensing matrices, which has their limitation in CI system implementation. In many CI systems investigated in literature, a DMD is the spatial light modulator for measurement collection. A DMD consists of thousands of micromirrors. These mirrors flip between two angles at a predefined frequency. During one measurement collection time T₀, if the mirrors use half of the time to reflect light to output, then the object light is modulated by a factor 0.5. Therefore, light modulation is based on time division in DMD. Due to this modulation methodology, for a fixed DMD working frequency f, collecting a measurement using a binary projection needs one period 1/f. However, for gray values, such as values which require 8 bit binary numbers for representation, DMD needs 256 periods to realize the light modulation. Therefore, the measurement needs much longer time(≫1/f). This slows down a CI system’s measurement collection and imaging speed. Besides this major disadvantage, using gray value projections, CI systems require more complicated system calibration process, thus potentially have larger systematical error compared to using binary projections. Because of these issues, we study binary sensing matrix design in this work.

In CI, we cannot design a sensing matrix without considering a reconstruction method. In the classical definition of CI, an object is reconstructed by solving an optimization problem leveraging object sparsity. Besides such nonlinear methods, in our previous works we also investigated the linear reconstruction Wiener operator [18,33]. The operator minimizes reconstruction mean square error (MSE). It can be written as $\mathrm{W}={\mathrm{R}}_{\overrightarrow{x}}{\mathrm{F}}^{\mathrm{T}}{({\text{FR}}_{\overrightarrow{x}}{\mathrm{F}}^{\mathrm{T}}+{\sigma }^2\mathrm{I})}^{-1}$, where ${\mathrm{R}}_{\overrightarrow{x}}$ is the object autocorrelation matrix, F is the sensing matrix, and σ² is the noise variance. This operator does not use object sparsity for reconstruction. Its performance is slightly worse than nonlinear methods [34, 35] for sparse signals when noise is low. However, the performance difference is small [18]. When random features are collected for objects of size (80×80) with dynamic range [0,255] and sparsity K = 400 in the wavelet domain, the reconstruction RMSE (root mean square error) values are 9.21 and 6.54 using Wiener operator and a nonlinear Nowak method [36], respectively [18]. The noise standard deviation is 0.21. Same level RMSE difference was observed using PCA (Principal Component Analysis) features. Although Wiener operator does not use object sparsity, it is defined based on an object’s second order statistical parameter, the object’s autocorrelation matrix ${\mathrm{R}}_{\overrightarrow{x}}$. This statistical parameter is not used in the nonlinear reconstruction methods of CI. Additionally, because noise is directly used to define Wiener operator, the operator is more tolerant to noise. Thus, when noise is high or object light level is low, to obtain a better reconstruction, object sparsity is not the only object prior which can be used for reconstruction. Wiener operator presents even better reconstruction in terms of smaller RMSE. Besides these advantages, more importantly to the work discussed in this paper, the computational speed using the Wiener operator is much faster than nonlinear iterative methods, especially with GPU implementation [8,11,12]. Therefore we use the Wiener operator for sensing matrix design.

Note that, Wiener-type operators recently have been used in geophysics, and medical imaging specifically in diffuse optical tomographic imaging [37–39]. In a near-surface geophysical problem, observed data is emulated using a linear layered-earth model by representing the data as a multiplication between a matrix and a model vector [37]. To solve the geophysical problem is to understand how well the calculated data match with the observed data, and then select those well predicted data. A Wiener-type matrix, data-resolution matrix is used to do the selection. Then another Wiener-type matrix, model-resolution matrix is used to evaluate how well the layered-earth model works with the selected data. In diffuse optical tomographic imaging, the forward problem is different from CI or geophysical problem [38,39]. However, the natural logarithm of the amplitude of experimental data can still be written as a function of biological tissue’s absorption coefficients. Using the Jacobian of this function and a regularization parameter, data-resolution matrix and model-resolution matrix for this problem can be defined. Then based on data-resolution matrix, independent measurements are chosen to reduce data collection time without compromising reconstructed image quality much [38]. Model-resolution matrix is used with basic pursuit deconvolution algorithm to improve reconstructed image quality [39].

With the focuses on binary sensing matrices and Wiener operator for reconstruction, the paper is organized as follows. Before studying binary projection design, we first use maximizing measurement SNR and minimizing reconstruction MSE as the criteria to search the optimal gray value sensing matrixes in Section 2. Then based on the solutions, we design binary matrices in Section 3. Experimental results are presented in Section 4 to demonstrate the idea. Conclusions are drawn in Section 5.

2. Sensing matrix design using measurement SNR and reconstruction MSE

2.1. Sensing matrix design based on maximizing measurement SNR

We first define the mathematic model for CI. In compressive imaging, system measurements can be represented as

In Eq. (2), the object autocorrelation matrix ${\mathrm{R}}_{\overrightarrow{x}}=E\{{\Vert \overrightarrow{x}{\overrightarrow{x}}^{\mathrm{T}}\Vert }^2\}$ can be decomposed as QDQ^T, where Q is an orthonormal matrix, and D is diagonal. The i_th column of Q, ${\overrightarrow{q}}_i$, is the eigenvector of matrix R. The diagonal elements of D are sorted in descending order, such as d₁ ≥ d₂ ≥ ⋯ ≥ d_N. d_i is the eigenvalue of R corresponding to ${\overrightarrow{q}}_i$. Note for natural images, all d_i are positive.

To obtain the optimal projection in terms of maximizing SNR, we solve the following optimization problem,

Using the decomposition of ${\mathrm{R}}_{\overrightarrow{x}}$, we rewrite Tr{FR_xF^T} as Tr{(FQ)D(FQ)^T}. If we redefine FQ as H, then problem P₁ becomes

The matrix H can be further written as H = [H₁ H₂], where H₁ and H₂ have sizes (M×M) and (M×(N−M)).

One special case for P₁ or P_1a is ${\mathrm{H}}_1{\mathrm{H}}_1^{\mathrm{T}}=\mathrm{I}$ and H₂ = 0. The cost function for this case is ${\displaystyle \sum _{i=1}^M{d}_i}$. If we define ${\mathrm{Q}}_{\text{PCA}}=[{\overrightarrow{q}}_1{\overrightarrow{q}}_2\cdots {\overrightarrow{q}}_M]$, then this special case can be written as $\mathrm{F}={\mathrm{H}}_1{\mathrm{Q}}_{\text{PCA}}^{\mathrm{T}}$. Note if H₁ = I, then $\mathrm{F}={\mathrm{Q}}_{\text{PCA}}^{\mathrm{T}}$ is the PCA projection. We modify F slightly as

2.2. Sensing matrix design based on minimizing the reconstruction MSE

In compressive imaging, the final goal is to obtain high quality object reconstructions. Instead of improving the intermediate measurement SNR, in this subsection we study the matrix design to minimize the reconstruction error.

As discussed in Sec. 1, we use Wiener operator for object reconstruction. The Wiener operator is defined as

Using W, the object is reconstructed as ${\overrightarrow{x}}_{\text{est}}=\mathrm{W}\overrightarrow{y}$. Then the reconstruction MSE becomes [40]

The first term in Eq. (7) is independent of F. To minimize the reconstruction error MSE_CS is equivalent to maximizing the second term. Thus, we define the design problem as

Due to the term ${\left({\text{FR}}_{\overrightarrow{x}}{\mathrm{F}}^{\mathrm{T}}+{\sigma }^2\mathrm{I}\right)}^{-1}$, solving P₂ with a close form solution, such as using KKT theory [41], is hard. Note that, for $\mathrm{F}={\text{AQ}}_{\text{PCA}}^{\mathrm{T}},\text{Tr}\left\{{\text{FR}}_{\overrightarrow{x}}^2{\mathrm{F}}^{\mathrm{T}}{\left({\text{FR}}_{\overrightarrow{x}}{\mathrm{F}}^{\mathrm{T}}+{\sigma }^2\mathrm{I}\right)}^{-1}\right\}$ is ${\displaystyle \sum _{j=1}^M{d}_i^2/(d_i+{\sigma }^2)}$. Same as in Section 2.1, if we assume $\mathrm{F}={\mathrm{Q}}_{\text{PCA}}^{\mathrm{T}}$ is slightly modified as

Note the multiple optimal solutions to P₂ provide more options for CI system implementation. A practical issue in CI is the high feature measurement dynamic range. To keep a good reconstruction performance, CI system measurement SNR needs to be large. Then, for projections such as PCA and Hadamard, the first feature collected is the summation of all object pixel values. This feature is significantly larger than others. Thus detectors having large dynamic range are required. However, as what we find in section 2.1 and 2.2, any matrix $\mathrm{F}={\text{AQ}}_{\text{PCA}}^{\mathrm{T}}$ presents the optimized measurement SNR or reconstruction MSE. Therefore, we can choose one among all of these matrices, which has the most relaxed dynamic range requirement to system detectors.

3. Optimized binary sensing matrices

In Sec. 2.2, we conclude that, ${\mathrm{F}}_{\text{opt}}={\text{AQ}}_{\text{PCA}}^{\mathrm{T}}$ presents the minimal MSE for any orthonomal matrix A. Note the row vectors of F_opt are M orthonomal vectors in a space ℚ, which is spanned by $\{{\overrightarrow{q}}_1,{\overrightarrow{q}}_2,\cdots ,{\overrightarrow{q}}_M\}$. A matrix F presents a small MSE, if its distance from the space ℚ, defined as ${\Vert \mathrm{F}-{\text{AQ}}_{\text{PCA}}^{\mathrm{T}}\Vert }_F$, is small. Hence we can consider searching the optimal binary projection to minimize MSE as to search a binary matrix $\mathit{sgn}\left\{{\text{AQ}}_{\text{PCA}}^{\mathrm{T}}\right\}$ close to ℚ. sgn() is the signum function, which has value {1,0, −1} for positive, zero, or negative input, respectively. Note that, if all elements in $\mathit{sgn}\left\{{\text{AQ}}_{\text{PCA}}^{\mathrm{T}}\right\}$ are assumed nonzero, then the L₂ norm of each sensing vector is N. To make both of the binary and the non-binary vectors having same L₂ norms, we multiply $\sqrt{N}$ to ${\text{AQ}}_{\text{PCA}}^{\mathrm{T}}$. Then designing binary matrices becomes solving the following problem,

Using these definitions, the Lagrangian for P₃ can be written as

P₃ is an equality-constrained optimization problem. To solve P₃, we use the Newton-Raphson method [41], which is specifically designed for nonlinear equality-constrained optimization problem. The algorithm can be represented using the following equations,

In Eq. (13) and (14), ${\nabla }_{\text{AA}}^2ℒ$ is the Hessian matrix of the Lagrangian, J is the Jacobi of the constraint for P₃, and ∇ε is the gradient vector of the cost function. The dimensions of ${\nabla }_{\text{AA}}^2ℒ$, J, and ∇ε are (M²×M²), $\left(\frac{M(M+1)}{2}\times M^2\right)$, and (M²×1), respectively. Different from Eq. (12), the variables in Eq. (13) and (14) are vectors $\overrightarrow{a}$ and $\overrightarrow{\lambda }$, which are obtained by lexicographically ordering the elements of A and the upper triangle part of λ. The vector $\overrightarrow{a}$ is defined as $\overrightarrow{a}={\left[A_{11}\cdots A_{1M}\cdots A_{M1}\cdots A_{MM}\right]}^{\mathrm{T}}$. $\overrightarrow{\lambda }$ is $\overrightarrow{\lambda }={\left[{\lambda }_{11}\cdots {\lambda }_{1M}{\lambda }_{22}\cdots {\lambda }_{2M}\cdots {\lambda }_{(M-1)M}{\lambda }_{MM}\right]}^{\mathrm{T}}$. Similarly, $\overrightarrow{g}$ is $\overrightarrow{g}={\left[{\mathrm{g}}_{11}\cdots {\mathrm{g}}_{1M}{\mathrm{g}}_{22}\cdots {\mathrm{g}}_{2M}\cdots {\mathrm{g}}_{(M-1)M}{\mathrm{g}}_{MM}\right]}^{\mathrm{T}}$. In this section, the definitions of ${\nabla }_{\text{AA}}^2ℒ$, J, ∇ε are presented directly. More details about the calculation of these parameters can be found in the Appendix section.

The elements of the Hessian matrix ${\nabla }_{\text{AA}}^2ℒ({\overrightarrow{a}}^{(\nu )},{\overrightarrow{\lambda }}^{(v)})$ are defined as

Function δ(i) is the delta function. For the Jacobi matrix of G, note that G is symmetric. Therefore only the upper triangle part of G is used to find J. We have

The third parameter in the Newton-Raphson method is $\nabla \epsilon (\overrightarrow{a})$ with

To ensure the step direction $\mathrm{\Delta }{\overrightarrow{a}}^{(\nu )}$ is a descent direction for $\nabla \epsilon ({\overrightarrow{a}}^{(v)})+J{\left({\overrightarrow{a}}^{(\nu )}\right)}^{\mathrm{T}}{\overrightarrow{\lambda }}^{(\nu )}$, the Hessian matrix ${\nabla }_{\text{AA}}^2ℒ$ is modified to become positive definite [41, 42]. In addition, the stepsize α in Eq. (14) is set to be small, and the number of the iterations ν is set to be large, so that ∇ε will converge to a cost value without oscillation.

4. Numerical results for fast compressive L³-imaging

As discussed in section 1, collecting features or linear combinations of object pixels will benefit system measurement SNR. Therefore, CI is well suitable for low-light-level imaging (L³-imaging) [6,7]. Figure 1 is a block-based system diagram for L³-imaging. In such a system, objects are defined into blocks. Then, for each block, multiple compressive measurements are collected using a detector. The advantage of using block-based CI (BCI) instead of CI with a single detector is to obtain high resolution object reconstructions [33].

 

Fig. 1 A system diagram for compressive L³-imaging.

Download Full Size | PDF

Figure 2(a) and 2(b) present the two kinds of objects used for the experimental study [43,44]. They are low-light-level images and are rather dark, so for illustrations only we enhanced their brightness in the figures. The sizes of the two kinds of objects are (198×128) and (1056×1920). Both of them have minimum 0 and maximum 76.5, although the mean value is 45.8 for the first kind, and 33.5 for the second. The block size is set to be (32×32). M = 8 and 16 features are collected for each block. Noise is assumed to be Gaussian noise with standard deviation σ in the range [5,150]. With such noise level, the SNR range is [−13.02dB,19.24dB]. To make this SNR range having more physical meaning, we use Hamamatsu SI PIN photodiode S1223 as an example to calculate the input light power. The detector S1223 has NEP (noise equivalent power) 9.4 × 10⁻¹⁵W/Hz^1/2. If the system bandwidth is assumed as 30MHz, then when SNR is 0dB, the input light power or the noise power is 51.486pW. In our experimental study, the lowest SNR is −13.02dB. At this SNR level, the input light power is 2.568pW. In other words, for such level of light power, directly using PIN photodiode S1223 for measurement collection, the input signal will be overwhelmed by noise. However, using a BCI system as shown in Fig. 1, we can still make the measurements, and get an object reconstructed.

 

Fig. 2 Object examples of size (a) (196 × 128) [43,44], and (b) (1056 × 1920).

Download Full Size | PDF

In the first experiment, we studied seven kinds of projections, such as three Hadamard projections, two DCT projections, PCA projection ${\mathrm{Q}}_{\text{PCA}}^{\mathrm{T}}$, and the projection ${\text{AQ}}_{\text{PCA}}^{\mathrm{T}}$. In the first Hadamard projection, H_freq, we sort the sensing vectors based on the variation of elements’ signs. For example, the first vector has all one values. The element’s sign varies 0 time. In the second vector, the first half elements are +1, and then the rest are −1. The sign varies once. In the third vector, the element values change from +1 to −1, and then back to +1. The sign varies twice. Note that, using this sorting method, the system collects the DC value, the low frequency components, and the high frequency components of the object x sequentially. In the second Hadamard projection, H_mean, we use a training sample average to sort the sensing vectors. The sensing vectors are in the same order as the Hadamard coefficients of the sample average. In the third Hadamard projection, H_snr, we sort the vectors using measurement SNR [7]. The two DCT for this experiment are the regular DCT, F_DCT, and the DCT with sensing vectors sorted using measurement SNR, ${\mathrm{F}}_{{\text{DCT}}_{\text{snr}}}$. Table 1 presents the cost function ε(F) = Tr{FR_xF^T} for problem P₁ using these projections. It can be observed, the features obtained using PCA and ${\text{AQ}}_{\text{PCA}}^{\mathrm{T}}$ have the largest ε(F) as we discussed in Section 2.1.

Table 1. The cost function ε(F) = Tr{FR_xF^T} in P₁ when M = 32 features are obtained using PCA, DCT, and Hadamard matrices. (^*The values listed in the table are ε(F)/10⁶.)

View Table | View all tables in this article

In Fig. 3(a) and 3(b), we present the normalized reconstruction root mean square error (RMSE) vs. noise σ for the seven kinds of projections. The RMSE is defined as $\Vert {\overrightarrow{x}}_{\text{est}}-\overrightarrow{x}\Vert /\Vert \overrightarrow{x}\Vert$. As expected, for both kinds of objects, ${\mathrm{Q}}_{\text{PCA}}^{\mathrm{T}}$ and ${\text{AQ}}_{\text{PCA}}^{\mathrm{T}}$ present the same reconstruction error performance at different noise levels. These errors are smaller than the errors obtained using other kinds of projections. This is consistent with the conclusion obtained in section 2.2, that the projection ${\text{AQ}}_{\text{PCA}}^{\mathrm{T}}$ for any orthonomal matrix A is optimal in terms of minimizing RMSE. Note that, comparing Fig. 3(a) and 3(b), we find the RMSE values for the large size object are smaller. The reason is that, for a large size object, pixel values in object blocks vary slowly, or blocks are sparse in frequency domain. Therefore, using a few feature measurements can reconstruct the block with small error.

 

Fig. 3 Reconstruction RMSE for objects of size (a) (196×128) and (b) (1056×1920), when M = 32 Hadamard, DCT, and PCA features are collected for (32×32) object blocks.

Download Full Size | PDF

In the second experiment, we use the iterative algorithm discussed in section 3 to search for the optimal binary projection for the object shown in Fig. 2(b). Note that, the solution to this algorithm is a local minimizer. To find a global solution, 10M² initializations with 4000 iterations in each trail are used, where M is the number of features, and M² is the number of unknown variables in problem P₃. For example, if M = 8, the number of initializations used to search the optimal matrix A of size (8×8) is 640. In addition, to maximize the object signal in a feature measurement, we normalize the maximum absolute value of each sensing vector to 1.

In Fig. 4(a), the 8 Hadamard vectors in H_snr, which are selected based on measurement SNR [7], are presented. Figure 4(b) presents the optimization result, F_rmse, for M = 8. It can be observed, except the all-one vector, the solution for P₃ is not directly related to Hadamard projection.

 

Fig. 4 (a) The Hadamard vectors sorted using measurement SNR, (b) the optimal binary vector set or the optimal solution for problem P₃, (c) the suboptimal solution for P₃ which has the minimal feature measurement dynamic range, when M is 8 for 32 × 32 blocks.

Download Full Size | PDF

Besides the optimal solution, we also compare the 640 local minimizers obtained by the 10M² trials. The RMSE values obtained using these local solutions are calculated. We find 217 local minimizers including the optimal one presents very similar RMSE. The RMSE difference between these local minimizers and the optimal solution is less than 10⁻⁴. Then among these sub-optimizers, the solution having the smallest feature dynamic range, labeled as F_dyn, is presented in Fig. 4(c). The all-one vector is not among the result. Table 2 summarizes the minima, the maxima, and the dynamic ranges of features obtained using H_mean, H_snr, ${\mathrm{F}}_{{\text{DCT}}_{\text{snr}}}$, PCA, F_rmse, and F_dyn. It is clear that, F_dyn needs about 30% less dynamic range than the other five projections.

Table 2. The dynamic ranges of M = 8 features obtained using different sensing matrices. (^*The values listed in the table are the actual values divided by 10⁴.)

View Table | View all tables in this article

Figure 5(a) presents the reconstruction errors using H_mean, H_snr, ${\mathrm{F}}_{{\text{DCT}}_{\text{snr}}}$, ${\mathrm{Q}}_{\text{PCA}}^{\mathrm{T}}$, Frmse, and F_dyn vs. noise σ. It can be observed that F_rmse presents even smaller RMSE than PCA, although it is binary. The reason is, after normalizing each vector having maximum 1, a system using F_rmse can collect more object light into a measurement than using ${\mathrm{Q}}_{\text{PCA}}^{\mathrm{T}}$. Thus system measurement SNR and reconstruction RMSE are improved. From the figure we can also observe that F_dyn presents similar RMSE values as F_rmse. As a comparison, we also present the reconstruction RMSE for the standard Hadamard and random binary projections in Fig. 5(b). It is clear that, both projections have larger reconstruction errors. To investigate the reason that the optimized binary projections work well, we calculate the Fourier transform of each projection vector. Figure 6 shows the normalized FFT of the vectors in the standard Hadamard, H_snr, and F_dyn matrices. The FFT for F_rmse is similar to F_dyn. From the figure, it is clear that, F_dyn will collect object information at more frequencies. H_snr will collect at fewer frequencies, while the standard Hadamard collects at the fewest. Thus the reconstruction using F_dyn has the smallest RMSE.

 

Fig. 5 When M = 8 features are collected for 32×32 blocks, the normalized RMSE using (a) the sorted Hadamard, DCT, PCA, F_rmse, and F_dyn matrices, and (b) the standard Hadamard, random binary, and F_rmse matrices.

Download Full Size | PDF

 

Fig. 6 When M = 8 features are collected for 32×32 blocks, the normalized Fourier transforms of the vectors in (a) the standard Hadamard matrix, (b) H_snr, and (c) F_dyn.

Download Full Size | PDF

The experiment is repeated for M = 16. Sensing matrices F_rmse and F_dyn are shown in Fig. 7(a) and 7(b). The reconstruction RMSE for H_mean, H_snr, ${\mathrm{F}}_{{\text{DCT}}_{\text{snr}}}$, ${\mathrm{Q}}_{\text{PCA}}^{\mathrm{T}}$, F_rmse, and F_dyn are plotted as functions of noise in Fig. 7(c). Once again, the optimal and the suboptimal sensing matrices present very similar reconstruction performance. We also calculate the feature dynamic range values. We find the feature dynamic range is reduced 28% from 7.47 × 10⁴ for H_snr to 5.39×10⁴ for F_dyn. Figure 7(d) and 7(e) presents the normalized FFT of the vectors in Hsnr and F_dyn matrices. Once again, we observe that F_dyn will collect more object information at more frequencies. Thus the system using this matrix has better reconstruction performance.

 

Fig. 7 (a) The optimal binary vector set or the optimal solution for problem P₃, (b) the suboptimal solution for P₃ which has the minimal feature measurement dynamic range, (c) the normalized RMSE for sorted Hadamard, DCT, PCA matrices, F_rmse, and F_dyn, and the normalized Fourier transforms of the vectors in (d) H_snr and (e) F_dyn, when M is 16 for 32 × 32 blocks.

Download Full Size | PDF

Figure 8 presents 4 reconstructions using different imaging systems and projections. Figure 8(a) is the average of 16 conventional imaging measurements. (b), (c), and (d) are the reconstructions obtained using the sensing matrix H_mean, H_snr, and F_dyn, where M = 16 measurements are collected, and the noise level is σ = 148.67. To make the reconstruction more visually pleasant, we increase the intensity 3 times. It is clear, all reconstructions using CI have smaller RMSE compared with the conventional result. The RMSE values for the conventional imaging, H_mean, H_snr, and F_dyn are 0.7312, 0.0733, 0.0645 and 0.0635, respectively. Note that, the reconstructions using H_snr, and F_dyn are very similar to each other, and are close to the original. To have a better observation, Fig. 9 shows two magnified detail parts of both reconstructions. Figure 9(a) and 9(b) are for H_snr. Figure 9(c) and 9(d) are for F_dyn. It can be observed that, the reconstruction using F_dyn has slightly smoother reconstruction with less block effect compared with the reconstruction using H_snr.

 

Fig. 8 The object reconstruction using (a) conventional imaging, (b) H_mean, (c) H_snr, and (d) F_dyn, when M = 16 features are collected for 32 × 32 blocks and σ is 148.67.

Download Full Size | PDF

 

Fig. 9 The object reconstructions (a)&(b) using H_snr, and (c)&(d) using F_dyn when M = 16 features are collected for 32 × 32 blocks and σ is 148.67.

Download Full Size | PDF

In the next experiment, we add a bias to each projection to make it nonnegative. Then the maximum of the obtained projection is normalize to 1. For the PCA and DCT projections, the vector values are quantized to 128 levels. As discussed in the beginning of the paper, CI systems using gray value projections, require longer time to make measurements. In the time to collect a measurement with a 8-bit projection, a system can make 128 measurements using binary vectors. Therefore, binary projections benefit a CI system imaging speed. Figure 10 presents the reconstruction RMSE with these projections. We get the same conclusions as in Fig. 7(c). The RMSE values for F_rmse and F_dyn are almost same. They are smaller than the RMSE values using other projections. Among the rest four projections, H_snr has the smallest reconstruction at almost all noise levels. Table 3 summarizes the maximal feature values or the dynamic ranges use all six projections. Once again, we can see F_dyn has the smallest dynamic range.

 

Fig. 10 Normalized RMSE for different sets of nonnegative sensing vectors when M = 16 features are collected for 32 × 32 block.

Download Full Size | PDF

Table 3. The dynamic ranges of M = 16 features obtained using different nonnegative matrices. (^*The values listed in the table are the actual values divided by 10⁴.)

View Table | View all tables in this article

In the last set of experiments, we studied the performance of F_rmse and F_dyn with sparse objects. The object blocks in Fig. 2(b) are transformed into the wavelet domain using the 4-level wavelet transformation matrix. Except the largest 64 coefficients in each block, the rest are set to zero. Using this sparse object, we repeat the experiment as in Fig. 10. Besides the six matrices we discussed in last several experiments, wavelet and curvelet projections are also studied. The wavelet projection matrix is for 4 level transformation. Both kinds of projection vectors are sorted using measurement SNR. These vectors are also added with a bias to make them nonnegative, and then quantized to 128 discrete levels. Figure 11 presents the reconstruction RMSE vs. noise σ. We observe similar trend as in Fig. 10, except that the RMSE values are smaller because the object sparsity benefits the reconstruction. For wavelet and curvelet transforms, the reconstruction error are larger than other projections at most noise levels. The reason is that, both transformation vectors are spatial localized. They are suitable to extract object’s local information and make the transformation sparse. However, this property is not helpful in compressive low-light-level imaging, in which more light from an object are required to surpass detector noise. Between wavelet and curvelet transforms, the performance of the latter is worse. This is caused by the non-orthogonality among the transformation vectors.

 

Fig. 11 Normalized RMSE for different sets of nonnegative sensing vectors when M = 16 features are collected for 32 × 32 blocks. The object blocks are sparse in the wavelet domain. The sparsity is K = 64.

Download Full Size | PDF

In the last experiment, we make the object further sparse with K = 4 for each block. The sparse object is shown in Fig. 12(a). Then we use a nonlinear method, fast linearized Bregman iterative algorithm, to reconstruct the object. In Fig. 12(b), 16 random binary vectors are used to make the feature measurements. In Fig. 12(c), F_dyn is used. The noise level in the experiment is σ = 148.67. We can observe that the reconstruction using F_dyn is visually better. The RMSE for Fig. 12(b) is 0.0623, while the RMSE for 12(c) is 0.0268. The designed binary sensing matrix works well for the sparse object. Because it is more structured in Fourier domain, it works better than random projections. Additionally, we also present the reconstruction using Wiener operator in Fig. 12(d) when F_dyn is used. The performance is even better than using the linearized Bregman iterative algorithm. As we discussed in the beginning of the paper, although Wiener operator does not use sparsity for reconstruction, it incorporates object second order statistical parameter, object autocorrelation matrix ${\mathrm{R}}_{\overrightarrow{x}}$. It is more suitable when system measurements have much noise. The RMSE value obtained using Wiener operator is 0.0131 for the sparse object.

 

Fig. 12 (a) The original object which is sparse in the wavelet domain. The sparsity for each 32 × 32 block is K = 4. The object reconstructions using (b) the fast linearized Bregman iterative method with a binary random matrix, (c) the fast linearized Bregman iterative method with F_dyn, and (d) Wiener operator with F_dyn, when M = 16 features are collected for each block and σ is 148.67.

Download Full Size | PDF

Table 4 summarizes the measurement collection time and the object reconstruction time using different kinds of sensing matrices and linear/nonlinear reconstruction methods. We assume TI Discovery 4100 development tool set is used to make measurements. Then the DMD working frequency is set as 32kHz. With this working frequency, for 16 binary features, we need 0.5ms to make the measurements. For 128-level gray sensing matrix, the time becomes 64ms. Using the linear reconstruction method, Wiener operator, 12ms is required to obtain the reconstruction. Using the nonlinear method, fast linearized Bregman iterative algorithm, the reconstruction time is more than 100s. Thus, to make the system imaging speed faster, binary sensing matrix and Wiener operator should be used.

Table 4. The (measurement collection time, reconstruction time) for different kinds of sensing matrices and linear/nonlinear methods. 16 feature are collected for each block. The DMD working frequency is set as 32kHz.

View Table | View all tables in this article

5. Conclusions

In this paper, we study to design binary sensing matrices for fast compressive L³-imaging. Two optimization problems, maximizing measurement SNR and minimizing RMSE, are solved first. We find the solutions to both problems are matrices ${\text{AQ}}_{\text{PCA}}^{\mathrm{T}}$, where A is orthonomal. Based on the observation, a binary matrix design problem is defined as to search a matrix $\mathit{sgn}({\text{AQ}}_{\text{PCA}}^{\mathrm{T}})$. Then from the experimental study, we find the optimal solution presents similar reconstruction as PCA, or even better than PCA when the maximum absolute values of projection vectors are normalized to 1. More interestingly, many suboptimal binary matrices are found to have almost the same reconstructions as the optimal. Among them, we can choose the matrix presenting the least measurement dynamic range to relax the requirement to CI sensors.

We focus on Gaussian noise in this work. When the light from object is extremely low, the noise becomes Poisson. For this case, the compressive imaging idea will still help the system’s reconstruction performance, because system measurement SNR will be improved using BCI. However, the linear Wiener operator and the nonlinear reconstruction algorithms used for CI need to be modified for Poisson noise. Besides Gaussian noise, in a physical system there are other kinds of detector noise, long-term drift error, and systematic errors from stray-light. Different kinds of detectors or imaging elements have different noise. For photomultipliers, photoconductive and photovoltaic detectors, shot noise, generation-recombination noise or g-r noise, thermal noise, and 1/f noise are the four main factors. Shot noise matches Poisson noise probability density distribution. When the signal frequency is not very high, g-r noise can be modeled as white noise. Thermal noise generally is white Gaussian noise. 1/f noise can be ignored in most cases when a signal frequency is higher than 1kHz. For different kinds of noise, using BCI, we can always have larger detector signal. When the signal power is larger, Poisson noise becomes Gaussian noise. For long-term drift, except calibrating the imaging system, using the optimized binary sensing matrix will shorten the measurement collection time. Thus it will help the system to deal with such kind of error. Stray-light will cause detector collecting undesired background or scattered light. For active imaging, modulating the illumination source and then preprocessing the detected signal using a bandpass filter will suppress stray-light in a system.

We study fast compressive L³-imaging in this paper. Thus, we use binary sensing matrices to make measurements and linear Wiener operator for object reconstruction. However, for applications in which imaging speed is not restricted and object light is not low or noise is not high, gray-value sensing matrices and nonlinear reconstruction methods have their advantages. To further improve the reconstruction performance of fast compressive L³-imaging, in the future we will investigate linear methods including object sparsity as an prior.