Image quality of compressive single-pixel imaging using different Hadamard orderings

Pedro G. Vaz; Daniela Amaral; L. F. Requicha Ferreira; Miguel Morgado; Miguel Morgado; João Cardoso

doi:10.1364/OE.387612

1. Introduction

Over the past several years, conventional cameras have relied on the Shannon-Nyquist sampling theorem. This well known theorem stipulates that the sampling frequency should be, at least, twice the highest frequency contained in the signal in order to avoid aliasing. Traditionally, higher image resolution implies an increase in the number of sensor elements and, consequently, an increased cost and complexity [1].

Recently, compressive sensing (CS) [2] has emerged as an alternative sampling approach that enables signal reconstruction, under certain conditions, using a limited sub-set of samples. CS theory is built upon two major fundamental conditions: sparsity and incoherence [3]. The first one is related to the nature of the signal, i.e., it must be sparse when expressed in a specific orthonormal basis ($\Psi$), for example the discrete cosine transform (DCT) or the wavelet transform. This characteristic implies that only a few coefficients will contain the majority of the signal information which, fortunately, is the case of most natural signals [4]. According to CS theory, the signal can be represented as:

(1)$$x = \Psi\theta \; ,$$

where $x \in \mathbb {R}^{n \times 1}$ is the vectorized signal and $\theta \in \mathbb {R}^{n \times 1}$ is the sparse vector that contains the projections of $x$ in the basis $\Psi \in \mathbb {R}^{n \times n}$. A two dimensional signal (image) must be reshaped for a column vector by stacking its columns or rows.

Incoherence is associated with the way the signal is sampled. The sensing matrix ($\Phi$), used to sample the scene, must have low coherence with $\Psi$. This is achieved using a matrix with independent identically distributed entries, e.g. $\pm 1$ binary entries [4,5]. One of the advantages of CS is that the matrix $\Phi$ is independent of the type or nature of the signal of interest, resulting in a sampling process that can be mathematically described as:

(2)$$y = \Phi x = \Phi \Psi \theta = A \theta \; ,$$

where $y \in \mathbb {R}^{m \times 1}$ is a column vector containing the projections of the sampled data on the sensing basis ($\Phi \in \mathbb {R}^{m \times n}$) and $A = \Phi \Psi$ is an $m \times n$ matrix denominated reconstruction matrix.

When these two conditions are satisfied, the signal can be reconstructed even when the sampled data ($y$) has much lower dimension than the original data ($x$) ($m \ll n$) [4]. In a CS approach, $\Phi$ is used to sample the scene while the reconstruction process is done using the matrix $A$ [6].

One application of CS is the single-pixel imaging (SPI) [7,8] which combines an innovative and simple hardware architecture with the mathematical theory behind CS. This technique uses structured illumination to sample the scene, while a single photodetector collects the light reflected by the target. Each illumination pattern corresponds to one vector of the sensing basis ($\Phi$), while the photodetector readout corresponds to the vector $y$. This setup allows the acquisition of an image with the projection of a number of patterns lower than the number of its pixels.

Since its discovery, several SPI systems have been developed [9]. Specifically, Studer et al. [10] developed a compressive fluorescence microscope with the possibility to execute hyperspectral imaging up to 128 channels. Moreover, Edgar et al. presented a single-pixel camera able to record, simultaneously, visible and infra-red real-time video [11]. In addition, several authors [12–14] demonstrated that a SPI system is able to detect the image even through scattering media. Finally, X-ray imaging, with a bucket detector, was also developed without the necessity of a highly coherent X-ray source [15,16], or as a way to reduce the radiation dose in medical imaging [17].

Considering that CS relies on the acquisition of a restricted number of measurements, the patterns used for structured illumination, or selective detection, must be carefully chosen. A good set of patterns will allow a better reconstruction of the image using less samples, decreasing the amount of time required to sample the target, and reducing the amount of memory necessary to store both the projection patterns and the sampled information.

Several works were recently conducted to study the influence of patterns selection on image reconstruction quality, especially using the Hadamard transform as sensing basis. Regarding compressive sensing, Yu [18] proposed a Hadamard ordering, denominated Cake-cutting, where the projection patterns are sorted by increasing number of connected component. Finaly, Yu and Liu [19] presented a complex ordering method termed Origami ordering. They have compared this method with the Russian Doll ordering [20], in an experimental setup, with better results.

This paper describes a study on the influence of the pattern projection order on the image reconstruction quality of low resolution planar images. In applications such as microscopy [21] and x-ray imaging [16], targets with shallow depth of field, as the ones used in this work, are commonly imaged. The experiment was performed both on simulated data and real data, acquired using a structured illumination single-pixel camera (SPC). This type of acquisition is used in applications with time constraints and when the target should be excited prior to acquisition.

The selected sensing basis corresponds to the Hadamard transform, which is composed by a set of Walsh functions. This sensing matrix was used to assess two simulated and two experimental planar targets. Five different sensing basis orderings (Natural, Walsh, Cake-cutting, High Frequency and Random) were tested using compression ratios (CR) from 1 to 0.1. Excluding the High Frequency order, these orders were selected because they represent standard arrangements of the Hadamard matrix that can be found in the literature and are already implemented in most programming languages. Finally, two image reconstruction algorithms, total variation minimization by augmented Lagrangian and alternating direction algorithms (TVAL3) and Nesterov’s algorithm (NESTA), were used in the reconstruction process.

2. Methods

2.1 Sensing matrix

The sensing basis is a key element in the performance of a CS system. The Hadamard transform matrix, known to be well-conditioned and suitable for applications like SPI [22,23], corresponds to a square matrix composed only by $\pm 1$ values. This matrix is formed by the set of Walsh functions that consist of rectangular pulses containing only values of +1 and -1. Together, they represent an orthonormal basis that can be interpreted as the Hadamard matrix and both its rows and columns form an orthogonal basis [20].

This work is focused on this type of sensing matrix, since it is commonly used due to easy hardware implementation [12,22,24]. However, other matrices, like Fourier [25], wavelet [26], Radon [27] or randomly distributed speckles [28] can also be used.

The Natural order Hadamard transform matrix can be computed as:

(3)$$H_{2^k} = \begin{bmatrix} H_{2^{k-1}} & H_{2^{k-1}}\\ H_{2^{k-1}} & -H_{2^{k-1}} \end{bmatrix} = H_2 \otimes H_{2^{k-1}} \; ,$$

where $H_1 = 1$, $2^k$ is the order of the Hadamard matrix, with $k$ a non-negative integer greater than or equal to 2, and $\otimes$ represents the Kronecker product. Due to its nature, this matrix is easily computed using recursive methods. Figure 1(a) shows an example of a Natural order Hadamard matrix for $k=4$.

Fig. 1. Different orderings used during this work. This schematic contains the 16^th order Hadamard matrices. White squares represent the value 1 while black squares represent -1. The numbers next to the matrices lines correspond to the position of the line in the natural order.

Download Full Size | PDF

In spite of both Natural ordered and Walsh ordered Hadamard transform matrices containing the same information, the latter corresponds to a different disposition of the Walsh functions. In this case, the Hadamard matrix is rearranged in such a way that its Walsh functions are ordered by an increasing number of zero crossings. Figure 1(b) represents a Hadamard matrix with Walsh arrangement. Zhuoran et al. [29] developed a simple algorithm to reorder the Natural into the Walsh order.

The Cake-cutting order Hadamard matrix was determined using the method described in [18]. For this arrangement, the Hadamard matrix rows are organized by increasing number of connected components appearing in the projected patterns. A connected component corresponds to a continuous group of pixels with the same value. In order to determine this ordering, the Walsh functions were reshaped to the desired projection pattern. Then, the Matlab function bwconncomp was used to determine the connected components. Figure 1(c) shows the 16^th order Hadamard matrix with Cake-cutting ordering. Finally, the High Frequency ordered Hadamard matrix is obtained by flipping, vertically, the Walsh ordered Hadamard matrix. Figure 1(d) presents the Hamadard matrix with High frequency arrangement for $k=4$.

2.2 Reconstruction algorithms

Total variation minimization by augmented Lagrangian and alternating direction algorithms (TVAL3) [30] and Nesterov’s algorithm (NESTA) [31] are the algorithms adopted in this work to reconstruct the image. This choice is due to their execution speed and good image reconstruction (low root mean square error) with only a few projections (down to 10$\%$), compared with other algorithms also used to reconstruct images associated to CS [32,33].

2.2.1 TVAL3

Unlike other reconstruction algorithms, total variation methods consider the signal’s gradient to be sparse [30]. This is equivalent to say that the transform $\Psi$ is the inverse gradient transform. So, the goal is to find the sparsest gradient that satisfies the constraint imposed by Eq. (2).

The optimization problem, as used in this work, can be defined as:

(4)$$\min_{\omega_{i},x}\Sigma_{i}\Vert \omega_{i}\Vert_2 , \; \textrm{subject to} \; \Phi x=y \; \textrm{and} \; D_{i}x=\omega_{i} \; \forall i,$$

where $\omega _i = D_{i}x \in \mathbb {R}^{2 \times 1}$ is the discrete gradient of $x$ at position $i$, horizontally and vertically, and $\Vert \dots \Vert _2$ is the $l_2$ norm. In spite of $x$ being one dimensional, its spatial characteristics must be preserved in order to determine the horizontal and vertical gradients.

The TVAL3 solver applies the augmented Lagrangian multiplier method to the objective function in the form:

(5)$$\begin{aligned}\min \mathcal{L}_{A}(\omega_{i}, x) &= \min_{\omega_i,x} \sum_{i} (\Vert \omega_{i}\Vert_2 - \nu_{i}^{T}(D_{i}x - \omega_{i} ) + \\ & + \frac{\beta_{i}}{2} \Vert D_{i}x - \omega_{i} \Vert_{2}^{2}) - \lambda^{T}(\Phi x-y) + \frac{\mu}{2} \Vert \Phi x-y \Vert_{2}^{2}\; ,\end{aligned}$$

where $\nu _i \in \mathbb {R}^{2 \times 1}$ and $\lambda \in \mathbb {R}^{m \times 1}$ are Lagrange multipliers and $\beta _i \in \mathbb {R}$ and $\mu \in \mathbb {R}$ are regularization parameters associated with penalty terms of each constraint.

In this work, the algorithm implementation provided by Li et al. [30] was used with the following set of parameters as starting conditions: $\nu _i = \texttt {zeros(2,1)} \; \forall \; i$, $\lambda = \texttt {zeros(m,1)}$, $\beta _i = 2^5 \; \forall \; i$ and $\mu = 2^8$. Moreover, other solver parameters have been used as stated in [34].

2.2.2 NESTA

The Nesterov’s algorithm is frequently used to minimize smooth convex functions that are differentiable everywhere [35]. This method solves the $l_1$ norm minimization problem under a quadratic constraint:

(6)$$\min_{x} \Vert x \Vert_{1} \; \textrm{subject to} \; \Vert y- \Phi x \Vert_{2}\leq\epsilon ,$$

where $\epsilon$ is a parameter that accounts for the uncertainty of the measurements $y$. Since $\Vert x \Vert _1$ is not a smooth function [31], this objective function should be reformulated as:

(7)$$\min_{x \in \mathcal{Q}_{p}} f_{\mu}(x) = \min \Big( \max_{\theta \in \mathcal{Q}_{d}}\;<\;\theta, x\;>\;{-} \frac{\mu}{2} \Vert \theta \Vert^{2}_2 \Big) ,$$

where $\mathcal {Q}_{d} = \{ \theta : \Vert \theta \Vert _{\infty }\leq 1 \}$ and $\mathcal {Q}_{p} = \{ x:\Vert y- \Phi x \Vert _{2}\leq \epsilon \}$ and $\mu$ corresponds to the smoothness parameter. Higher values of $\mu$ result in a faster algorithm convergence while lower values for $\mu$ provide better accuracy [31].

The algorithm implementation provided by J. Bobin and S. Becker [36] was used in this work with the following set of parameters as starting conditions: $\epsilon = 4 \times 10^{-2}$, $\mu = 0.02$ and the stopping criteria $\nabla f_\mu (x) \leq 10^{-6}$. The basis $\Psi$ was selected as the DCT. Moreover, other solver parameters have been used as stated in [31].

In both reconstruction algorithms, permutation of the columns of the sensing matrix was not allowed. In a simulation environment, this permutation consists of mixing the pixels of the target image before the sampling procedure while in a real application it is achieved by projecting on the target a randomized version of each Walsh function.

2.3 Simulated data

A simulation tool was developed to test the sensing matrices using controlled conditions. This tool performs both simulation of target sampling and image reconstruction. Two different images from image processing datasets were used (Fig. 2). Both images are in greyscale with a resolution of 128 $\times$ 128 pixels.

Fig. 2. Images analysed during the simulated and experimental test.

Download Full Size | PDF

The image sampling procedure was performed according to the Eq. (2) where $x$ is a vectorized version of the image and $\Phi$ the desired sensing matrix. Furthermore, noise was added to the projections to simulate the noise present in the acquisition process. Normally distributed noise was added, using Matlab function randn, according to the following equation:

(8)$$y_s = y + c\times\overline{|y|}\times \sigma ,$$

where $y_s$ are the simulated signal projections, $c$ a constant with values of 0.1, 0.5 and 1, $\overline {|y|}$ is the mean of the absolute value of all projections and $\sigma$ an independent and normally distributed random variable with zero mean and unit standard deviation.

The simulated projections are then yielded to the two reconstruction algorithms. The number of projections ($m$) supplied to the algorithm was changed to achieve compression ratios (CR) of 0.1, 0.2, 0.3, 0.5, 0.8 and 1. The CR can be defined as the number of projections used in the reconstruction divided by the number of pixels of the image ($m/n$). This implies that a smaller value of CR results in higher data compression. Five different simulated acquisitions were performed for each set of parameters (CR, noise level and reconstruction algorithm). For each simulated acquisition, a different set of normally distributed random variable was also generated, resulting in variations of the reconstructed images. The results for the evaluation metrics are presented as mean $\pm$ standard deviation in order to illustrate these variations.

2.4 Experimental data

An experimental setup was used to perform the same test using a SPC prototype. The experimental SPC is composed by three main components: the pattern projection system, the amplified photodiode and the data acquisition system. All these components were mounted in an optical table to improve mechanical stability.

Figure 3 shows the experimental setup assembly. The selected pattern projection system was based on a digital light projector (DLP LightCrafter ^TM 4500 - Texas Instruments) that includes a digital micromirror device (DMD) as the intensity modulation element. A photodiode receiver module (Edmund Optics Si (300-1000nm), 10.0mm, photodiode receiver module #57-627) was coupled with two absorptive neutral density filters (0.6 OD and 0.9 OD - Edmund Optics #55-222) and used as sensor element. The neutral density filters were used prevent the photodiode from saturating. A data acquisition system (DAS - National Instruments ^TM DAQ X USB-6361) was used to digitize the photodiode output voltage. Finally, a focusing lens (+40 mm) was placed in front of the DLP optical engine to focus the projected pattern into a target placed at short distance. Both the DAS and the DLP were controlled using a Python 3.7^TM program.

Fig. 3. Experimental setup. DAS - Data acquisition system; DLP - Digital light projector; PD - Photodiode; FL - Focusing lens; T - Target (cardboard).

Download Full Size | PDF

All the experimental acquisitions were obtained with a sampling frequency of 20kS/s, a pattern exposure time of 83 microseconds and illumination by white light. This exposure time was selected because it corresponds to a multiple of the DLP working frequency (60 Hz). The experimental tests were done using reproductions of the images in Fig. 2 with dimensions of 5 cm $\times$ 5 cm.

A sensing matrix of order 4096 ($k=12$) was used to produce the illumination patterns. Each pattern corresponds to a reshape of each line of the sensing matrix. Since these Walsh functions contain positive and negative numbers, but the DLP can only reproduce two states (On/Off), both the positive part and negative part of each pattern were projected on the target [22,37]. This procedure results in the determination of two coefficients (one for the positive part and other for the negative part) for each Walsh function, i.e., a total of 8192 coefficients.

The data collection only occurred after the stabilization of the photodiode output, resulting in around 1300 samples for each pattern (65 ms). The mean value of these samples was considered to be the coefficient for the projected pattern. The resulting coefficient for each Walsh function ($y_e$) was determined as:

(9)$$y_e = \overline{y_+} - \overline{y_-}\;,$$

where $y_e$ is the coefficient, $\overline {y_+}$ is the mean value obtained for the positive part of each Walsh function and $\overline {y_-}$ is the mean value obtained for the negative part of the Walsh function. The complete set of projections was then used to reconstruct the target images with TVAL3 and NESTA using the same parameters of Section 2.3.

2.5 Evaluation metrics

During simulation, the reconstructed images ($\hat {I}$) were evaluated using two different metrics apart from visual inspection. The image contrast was computed as [38]:

(10)$$K = \dfrac{\sigma_{\hat{I}}}{\langle \hat{I} \rangle} \: ,$$

where $\sigma _{\hat {I}}$ is the standard deviation of the reconstructed image and $\langle \hat {I} \rangle$ its mean.

The second metric was the structural similarity index (SSIM) [39]. This metric is focused on three image parameters: luminance ($l$), contrast ($c$) and structure ($s$). Each of these parameters is evaluated using the following equations:

(11)$$\begin{aligned} l(\hat{I},I) & = \frac{2\mu_{\hat{I}} \mu_I + C_1}{\mu_{\hat{I}}^2 + \mu_I^2+ C_1} \: ,\\ c(\hat{I},I) & = \frac{2\sigma_{\hat{I}} \sigma_I + C_2}{\sigma_{\hat{I}}^2 + \sigma_I^2 + C_2} \: ,\\ s(\hat{I},I) & = \frac{\sigma_{\hat{I},I} + C_3}{\sigma_{\hat{I}} \sigma_I+C_3} \: ,\\ \textrm{SSIM}(\hat{I},I) & = l(\hat{I},I)^\alpha \times c(\hat{I},I)^\beta \times s(\hat{I},I)^\gamma \: , \end{aligned}$$

where $\mu$ stands for the image mean value, $\sigma$ is the image standard deviation, $\sigma _{\hat {I},I}$ represents the covariance between $\hat {I}$ and $I$ and $C_1$, $C_2$ and $C_3$ are constants used to prevent indeterminate expressions.

This metric was determined using Matlab function ssim with $C_1 = 2.55^2$, $C_2 = 7.65^2$, $C3 = C2/2$ and $\alpha = \beta = \gamma = 1$. The peak signal-to-noise ratio (PSNR), a common metric in the image processing field, was not used in this work because there is evidence that, in some cases, it does not produce results coherent with the visual perception [39].

3. Results

The results hereby presented include those of both the simulation and experimental setup. All the dataset, including reconstructed images and values of contrast and SSIM can be found in [40].

3.1 Simulated data

The main numerical results regarding the simulation data can be found in Figs. 4 and 5 while the reconstructed images are presented in Figs. 6 and 7. Only the reconstructed images with compression ratios 0.1, 0.3 and 0.5, TVAL3 reconstruction algorithm and noise level $c=0.5$ are presented.

Fig. 4. Numerical results of image reconstruction using TVAL3 and a noise level of $c = 0.5$. The left side plots represent the contrast of each reconstructed image while the right side represents SSIM, both as a function of compression ratios for cameraman (top) and boat (bottom) targets.

Download Full Size | PDF

Fig. 5. Numerical results of image reconstruction using NESTA and a noise level of $c = 0.5$. The left side plots represent the contrast of each reconstructed image while the right side represents SSIM, both as a function of compression ratios for cameraman (top) and boat (bottom) targets.

Download Full Size | PDF

Fig. 6. Image reconstruction of cameraman using TVAL3 algorithm and a noise level of $c = 0.5$.

Download Full Size | PDF

Fig. 7. Image reconstruction of boat using TVAL3 algorithm and a noise level of $c = 0.5$.

Download Full Size | PDF

Regarding the TVAL3 algorithm, it can be seen that SSIM increases with the increase of the CR, for all the orders. The Walsh and the Cake-cutting orders show better results, followed, in most cases, by the Random order. Walsh order presents, in general, higher contrast values while Cake-cutting order presents a significant better SSIM, specially for lower CRs. For instance, the results for the boat image, reconstructed with the Walsh and Cake-cutting orders and CR=0.1 (Figs. 7(j) and 7(m)), achieve a SSIM of 0.26 and 0.38, respectively. For the cameraman case, the values are closer, but the Cake-cutting order still do better than the Walsh order. In fact, for the SSIM with CR = 0.1, the error bars of Walsh and Cake-cutting orders intersect (Fig. 4(b))

The Walsh order consistently achieves higher values for contrast with the exception of the boat image for CR = 0.1. In this case, the Natural order achieves a value of 0.35, which is much larger than the other cases. This effect occurs because, when a low CR is used, the natural order reconstructs the image with a strong aliasing effect, producing several copies of the real image. These copies result in an image with two white and one black stripes (Fig. 7(a)) leading to a high value of contrast. As the CR increases, the contrast decreases as well as the aliasing effect, rising again when the CR approaches 1 (Fig. 4(c)). The same effect is visible in the cameraman case (Fig. 4(a)). The differences in the SSIM between both orders are attenuated with the increase of CR. For the boat image, Walsh and Cake-cutting orders achieve almost the same SSIM when considering a CR=0.5 (Fig. 4(b)).

Regarding NESTA, the Walsh and Cake-cutting orders achieve, again, the best results in terms of contrast and SSIM. As in the TVAL3 case, the Walsh order shows higher contrast values for low CR while the Cake-cutting shows a better SSIM. As the CR increases, the values of SSIM achieved by Walsh order match the ones obtained by Cake-cutting order. For example, for the cameraman image and CR=0.3, the values of SSIM for the Walsh and Cake-cutting orders are inside each other error bars (Fig. 5(d)).

Unlike the previous algorithm, the third best order is not clear in this case. If we consider SSIM, the Natural order occupies the third position, followed by the random order. By comparing the boat image reconstructed using NESTA for a CR of 0.5 (dataset [40] - Figs. 34(a) and 34(d)) we can state that the Random order presents the correct image structure while the Natural order shows superimposed copies of the original image. The value of SSIM contradicts this conclusion with the Natural order obtaining a value of 0.47 while the Random order achieves a value of 0.37. Therefore, a visual analysis should never be discarded.

Finally, the TVAL3 algorithm showed, globally, higher values of SSIM when comparing with NESTA. Due to length restrictions, we just present the images from the TVAL3 algorithm, since this was the algorithm that obtained better numerical results.

The structural differences between the reconstructed images are very clear both on cameraman (Fig. 6) and boat (Fig. 7). The Natural order presents superimposed copies of the image, the Random order shows a watercolor like pattern, the High Frequency order highlights the figure edges, and both Walsh and Cake-cutting orders can reconstruct the image with higher fidelity.

The reconstruction quality of the Walsh and Cake-cutting orders, even in the presence of noise and for a CR of 0.3, allows to identify background details, thus confirming their good performance. For lower CRs, the Cake-cutting order shows the best image quality. This fact is coherent with the better SSIM obtained for this order when CR=0.3. The background details in Fig. 6(n), are visible, although with a small SNR. This does not occur for the other orders. The same conclusion can be reached by analysing Fig. 7, where the only order that recovered the structural boat shape with CR=0.1 was the Cake-cutting.

Moreover, it is also clear that the images reconstructed with Walsh order for CR=1 presents a vertical blurring due to the orientation of the Hadamard matrices. The blurring would have been horizontally if the matrices had been transposed before projection. On the contrary, Cake-cutting order does not present a principal blurring orientation because the projected matrices are similar to its transposed versions. This effect can be regarded as an advantage of the Cake-cutting order over the Walsh order.

From the analysis of Figs. 6 and 7 it can be concluded that the Walsh order reproduces both the cameraman and boat images with a better contrast than the Cake-cutting order. This difference is particularly evident in a direct comparison between Fig. 6(k) and Fig. 6(n).

Regarding the High Frequency order, its performance is better on the boat image. This is expected, since the boat image has many high frequencies, and can be confirmed by visual inspection and numerical analysis.

3.2 Experimental data

Figures 8 and 9 show the reconstructed images for all the orders, using CR from 0.1 to 0.5 and TVAL3 algorithm, obtained with the experimental setup described in section 2.4. These images have been reconstructed with a resolution of 64 $\times$ 64 pixels.

Fig. 8. Image reconstruction of cameraman using TVAL3 with data from the bench experiment.

Download Full Size | PDF

Fig. 9. Image reconstruction of boat using TVAL3 with data from the bench experiment.

Download Full Size | PDF

Comparing these results with the simulation, a global decrease in the quality of the reconstructed images can be observed. This is due to the presence of unaccounted noise sources in the experimental apparatus. Among those noise sources, the large noise equivalent power of the photodiode ($10^{-6}$ V/Hz$^{1/2}$) and the projection system focus that slightly blurs the pattern on the test image must be highlighted. This is, in fact, one of the limitations of SPI using structured illumination instead of selective light collection [8].

Despite the decrease of SNR, the effects observed in the simulations are also present in the experimental data. The Cake-cutting and Walsh orders obtained a good performance and allow for reliable image reconstruction even when the CR is equal to 0.3. In addition, the Natural order causes image artifacts overlay; the Random order evidences a watercolor pattern; and, the High Frequency order highlights the image boundaries.

By comparing the Cake-cutting order with the Walsh order, the experimental results follow the same trend as those from the simulated data. The Cake-cutting order achieved a better performance, at low CRs, than the Walsh order. In the cameraman target, the structure of the image foreground is better reconstructed with the Cake-cutting order (Fig. 8(m)) compared with the Walsh order (Fig. 8(j)), especially the zones of the head and shoulder. Moreover, the Walsh order reconstruction presents the vertical blurring, similar to an aliasing artifact, already presented in the experimental results, that is not visible in the Cake-cutting reconstruction. The same conclusion can be made for the boat case, where the boat outline was reconstructed when using the Cake-cutting order (Fig. 9(m)) but it was not when using the Walsh order (Fig. 9(j)).

By comparing the two algorithms, it can be concluded that, there is a significant difference between the quality of the reconstructed image using NESTA or TVAL3 (dataset [40] - Figs. 49-56). In addition, a boundary conditions artifact can be seen on the images reconstructed with Cake-cutting order and NESTA (dataset [40] - Figs. 54(m) and Figs. 56(m)) at the top left of each image. This artifact leads to a lower contrast image when comparing with the Walsh order. The artifact disappears for reconstruction without compression. This situating could be corrected by removing the image boundaries. However it will lead to the reduction of the image resolution. The lack of a digital reference image invalidates the use of SSIM.

4. Conclusions

Five different Hadamard basis ordering were tested during this work. Both simulation and experimental results have demonstrated that the Cake-cutting order is the best choice to reconstruct natural images using CR of 0.1. For this compression ratio, the Walsh order presented a vertical blurring that Cake-cutting ordering avoided. For CR of 0.3 and 0.5, both the Cake-cutting and Walsh orders presented similar results. However, the Walsh order presented images with better contrast.

It is also important to mention the performance of the High Frequency order. From the simulation results, one can see that this order highlights the objects in the foreground making easy to detect their contours. This is, somehow, similar to a high-pass spatial filtering. The knowledge that emerges from this study opens the possibility to develop an adaptive method to select Walsh functions that maximize the image quality using a priory information about the analyzed scene.

Regarding the Natural order, all the images presented a superimposed artifact that strongly reduces the image perception. This effect occurs because, during our experiment, the permutation of the columns of the sensing matrix was not allowed. When column permutation is allowed, the influence of the Hadamard matrix order is strongly reduced because the projected patterns become similar to a speckle pattern, independently of the order. Also, the use of column permutation introduces a new variable that should be studied in a future work.

The influence of the basis ordering on the quality of the reconstructed image was demonstrated using qualitative and quantitative methods. A good image quality reconstruction, using an experimental setting with structured illumination, was achieved with a CR of 0.3. This value was referenced as a threshold when using compressive sensing in low resolution images, which is the case of this study [18]. The coefficients obtained with SPI for low resolution images present a low degree of sparsity, which causes this mathematical tool to fail [37]. In contrast, better CR can be achieved, especially when reconstructing images with larger resolutions [41] or instrumentation with higher SNR is used [42].

Funding

Fundação para a Ciência e a Tecnologia (PTDC/EMD-TLM/30295/2017); European Regional Development Fund; Competitiveness and Internationalization Operational Programme (PT-COMPETE 2020).

Disclosures

The authors declare no conflicts of interest.

References

1. E. Candes, T. Tao, J. Romberg, and R. Baraniuk, “Compressed sensing makes every pixel count,” in What’s Happening in the Mathematical Sciences, vol. 7D. Mackenzied, ed. (American Mathematical Society, 2009), chap. 7, pp. 114–127.

2. Y. C. Eldar and G. Kutyniok, eds., Compressed Sensing (Cambridge University, 2009).

3. R. G. Baraniuk, “Compressive sensing [lecture notes],” IEEE Signal Process. Mag. 24(4), 118–121 (2007). [CrossRef]

4. E. Candes and M. Wakin, “An introduction to compressive sampling,” IEEE Signal Process. Mag. 25(2), 21–30 (2008). [CrossRef]

5. T. Hong and Z. Zhu, “Online learning sensing matrix and sparsifying dictionary simultaneously for compressive sensing,” Signal Process. 153, 188–196 (2018). [CrossRef]

6. I. Rish and G. Grabarnik, Sparse modeling: theory, algorithms, and applications (CRC, 2014).

7. D. Takhar, J. N. Laska, M. B. Wakin, M. F. Duarte, D. Baron, S. Sarvotham, K. F. Kelly, and R. G. Baraniuk, “A new compressive imaging camera architecture using optical-domain compression,” in Computational Imaging IV, C. A. Bouman, E. L. Miller, and I. Pollak, eds. (SPIE, 2006).

8. M. F. Duarte, M. A. Davenport, D. Takhar, J. N. Laska, T. Sun, K. F. Kelly, and R. G. Baraniuk, “Single-pixel imaging via compressive sampling,” IEEE Signal Process. Mag. 25(2), 83–91 (2008). [CrossRef]

9. M. Rani, S. B. Dhok, and R. B. Deshmukh, “A systematic review of compressive sensing: Concepts, implementations and applications,” IEEE Access 6, 4875–4894 (2018). [CrossRef]

10. V. Studer, J. Bobin, M. Chahid, H. S. Mousavi, E. Candes, and M. Dahan, “Compressive fluorescence microscopy for biological and hyperspectral imaging,” Proc. Natl. Acad. Sci. 109(26), E1679–E1687 (2012). [CrossRef]

11. M. P. Edgar, G. M. Gibson, R. W. Bowman, B. Sun, N. Radwell, K. J. Mitchell, S. S. Welsh, and M. J. Padgett, “Simultaneous real-time visible and infrared video with single-pixel detectors,” Sci. Rep. 5(1), 10669 (2015). [CrossRef]

12. E. Tajahuerce, V. Durán, P. Clemente, E. Irles, F. Soldevila, P. Andrés, and J. Lancis, “Image transmission through dynamic scattering media by single-pixel photodetection,” Opt. Express 22(14), 16945–16955 (2014). [CrossRef]

13. L. Martínez-León, P. Clemente, Y. Mori, V. Climent, J. Lancis, and E. Tajahuerce, “Single-pixel digital holography with phase-encoded illumination,” Opt. Express 25(5), 4975 (2017). [CrossRef]

14. R. Dutta, S. Manzanera, A. Gambín-Regadera, E. Irles, E. Tajahuerce, J. Lancis, and P. Artal, “Single-pixel imaging of the retina through scattering media,” Biomed. Opt. Express 10(8), 4159–4167 (2019). [CrossRef]

15. H. Yu, R. Lu, S. Han, H. Xie, G. Du, T. Xiao, and D. Zhu, “Fourier-transform ghost imaging with hard x rays,” Phys. Rev. Lett. 117(11), 113901 (2016). [CrossRef]

16. Y. Klein, A. Schori, I. P. Dolbnya, K. Sawhney, and S. Shwartz, “X-ray computational ghost imaging with single-pixel detector,” Opt. Express 27(3), 3284–3293 (2019). [CrossRef]

17. D. Pelliccia, A. Rack, M. Scheel, V. Cantelli, and D. M. Paganin, “Experimental x-ray ghost imaging,” Phys. Rev. Lett. 117(11), 113902 (2016). [CrossRef]

18. W.-K. Yu, “Super sub-nyquist single-pixel imaging by means of cake-cutting hadamard basis sort,” Sensors 19(19), 4122 (2019). [CrossRef]

19. W.-K. Yu and Y.-M. Liu, “Single-pixel imaging with origami pattern construction,” Sensors 19(23), 5135 (2019). [CrossRef]

20. M.-J. Sun, L.-T. Meng, M. P. Edgar, M. J. Padgett, and N. Radwell, “A russian dolls ordering of the hadamard basis for compressive single-pixel imaging,” Sci. Rep. 7(1), 3464 (2017). [CrossRef]

21. A. Rodríguez, P. Clemente, E. Tajahuerce, and J. Lancis, “Dual-mode optical microscope based on single-pixel imaging,” Opt. Lasers Eng. 82, 87–94 (2016). [CrossRef]

22. Z. Zhang, X. Wang, G. Zheng, and J. Zhong, “Hadamard single-pixel imaging versus fourier single-pixel imaging,” Opt. Express 25(16), 19619 (2017). [CrossRef]

23. Y. Tsaig and D. L. Donoho, “Extensions of compressed sensing,” Signal Process. 86(3), 549–571 (2006). [CrossRef]

24. A. Hedayat and W. D. Wallis, “Hadamard matrices and their applications,” Ann. Statist. 6(6), 1184–1238 (1978). [CrossRef]

25. Z. Zhang, X. Ma, and J. Zhong, “Single-pixel imaging by means of fourier spectrum acquisition,” Nat. Commun. 6(1), 6225 (2015). [CrossRef]

26. K. M. Czajkowski, A. Pastuszczak, and R. Kotyński, “Single-pixel imaging with morlet wavelet correlated random patterns,” Sci. Rep. 8(1), 466 (2018). [CrossRef]

27. S. Dongfeng, H. Jian, M. Wenwen, Y. Kaixin, S. Baoqing, W. Yingjian, Y. Kee, X. Chenbo, L. Dong, and Z. Wenyue, “Radon single-pixel imaging with projective sampling,” Opt. Express 27(10), 14594–14609 (2019). [CrossRef]

28. Y. Bromberg, O. Katz, and Y. Silberberg, “Ghost imaging with a single detector,” Phys. Rev. A 79(5), 053840 (2009). [CrossRef]

29. C. Zhuoran, Z. Honglin, J. Min, W. Gang, and S. Jingshi, “An improved hadamard measurement matrix based on walsh code for compressive sensing,” in 2013 9th International Conference on Information, Communications & Signal Processing, (IEEE, 2013).

30. C. Li, W. Yin, H. Jiang, and Y. Zhang, “An efficient augmented lagrangian method with applications to total variation minimization,” Comput. Optim. Appl. 56(3), 507–530 (2013). [CrossRef]

31. S. Becker, J. Bobin, and E. J. Candès, “NESTA: A fast and accurate first-order method for sparse recovery,” SIAM J. Imaging Sci. 4(1), 1–39 (2011). [CrossRef]

32. Y. Kang, Y.-P. Yao, Z.-H. Kang, L. Ma, and T.-Y. Zhang, “Performance analysis of compressive ghost imaging based on different signal reconstruction techniques,” J. Opt. Soc. Am. A 32(6), 1063 (2015). [CrossRef]

33. L. Bian, J. Suo, Q. Dai, and F. Chen, “Experimental comparison of single-pixel imaging algorithms,” J. Opt. Soc. Am. A 35(1), 78–87 (2018). [CrossRef]

34. C. Li, W. Yin, and Y. Zhang, “User’s guide for tval3: Tv minimization by augmented lagrangianand alternating direction algorithms,” https://www.caam.rice.edu/~optimization/L1/TVAL3/v.beta/User_Guide_beta2.4.pdf. Date of access: 12/11/2018.

35. Y. Nesterov, “Smooth minimization of non-smooth functions,” Math. Program. 103(1), 127–152 (2005). [CrossRef]

36. J. Bobin and S. Becker, “Nesta - a fast and accurate first-order method for sparse recovery,” https://statweb.stanford.edu/~candes/software/nesta/. Date of access: 21/12/2018.

37. M.-J. Sun, M. P. Edgar, G. M. Gibson, B. Sun, N. Radwell, R. Lamb, and M. J. Padgett, “Single-pixel three-dimensional imaging with time-based depth resolution,” Nat. Commun. 7(1), 12010 (2016). [CrossRef]

38. P. G. Vaz, A. Humeau-Heurtier, E. Figueiras, C. Correia, and J. Cardoso, “Laser speckle imaging to monitor microvascular blood flow: A review,” IEEE Rev. Biomed. Eng. 9, 106–120 (2016). [CrossRef]

39. Z. Wang, A. Bovik, H. Sheikh, and E. Simoncelli, “Image quality assessment: From error visibility to structural similarity,” IEEE Trans. on Image Process. 13(4), 600–612 (2004). [CrossRef]

40. P. G. Vaz and J. Cardoso, Supplementary data of single pixel imaging experiment, Mendeley Data, http://dx.doi.org/10.17632/cymynmkxmj.6 (2020).

41. S. J. Olivas, Y. Rachlin, L. Gu, B. Gardiner, R. Dawson, J.-P. Laine, and J. E. Ford, “Characterization of a compressive imaging system using laboratory and natural light scenes,” Appl. Opt. 52(19), 4515–4526 (2013). [CrossRef]

42. C. Zhou, T. Tian, C. Gao, W. Gong, and L. Song, “Multi-resolution progressive computational ghost imaging,” J. Opt. 21(5), 055702 (2019). [CrossRef]

Image quality of compressive single-pixel imaging using different Hadamard orderings

Abstract

1. Introduction

2. Methods

2.1 Sensing matrix

2.2 Reconstruction algorithms

2.2.1 TVAL3

2.2.2 NESTA

2.3 Simulated data

2.4 Experimental data

2.5 Evaluation metrics

3. Results

3.1 Simulated data

3.2 Experimental data

4. Conclusions

Funding

Disclosures

References

Cited By

Figures (9)

Equations (11)

Optics Express