Real-time single-pixel video imaging with Fourier domain regularization

Krzysztof M. Czajkowski; Anna Pastuszczak; Rafał Kotyński

doi:10.1364/OE.26.020009

1. Introduction

Indirect image measurement techniques called single-pixel imaging and computational ghost imaging [1,2] contribute to many novel ideas in optics. Prospect applications of these measurement methods include spectral imaging [3–6], polarimetric imaging [7–9], 3D imaging [10–12], around-the-corner imaging [13], imaging through scattering media [14], imaging at low light levels [15], spectroscopy [16], and pattern recognition [17]. The inherently indirect principle of the measurement in single-pixel detectors makes it necessary to reconstruct the image numerically after the measurement. This is computationally demanding, and for compressive measurements, the reconstruction is ambiguous and typically involves optimization methods. The framework of compressive sensing [18–22] offers methods that provide very good reconstruction quality and are based on constrained ℓ¹-norm minimization of a sparse image representation. However, single-step reconstruction methods with regularization [22,23] used to solve the ill-posed inverse problem are significantly faster, at least for moderately sized measurements.

One of the greatest challenges for single-pixel cameras is real-time imaging at a video rate [24–28]. The video frame rate is limited by the method of sampling itself, requiring a spatial light modulator (SLM) to be switched to different states as many times as the number of captured image samples. For compressive measurement, this number is usually a small fraction of the pixel resolution of the image. State-of-the-art SLMs offer sampling rates above 20 kHz. The modulation speed of structured illumination in computational ghost imaging setups may be further increased by combining various light modulation techniques together or by using arrayed light sources with a fast modulation rate [29–31]. However, iterative optimization algorithms are usually not suffuciently fast to allow real-time image reconstruction at a similar rate. Therefore, effucient single-step reconstruction methods are necessary for single-pixel video imaging, providing reasonable quality of the recovered images.

In this paper we focus on the closed-form solutions to the inverse problem, which may be cast to the calculation of the Moore-Penrose pseudoinverse. The product of the pseudoinverse of the measurement matrix and the vector containing measurements captured by a single-pixel camera provides an estimate for the measured image, and has been used for single-pixel detection in conjunction with various nonorthogonal sampling protocols [32–34]. On the other hand, when the measurement patterns represent a subset of functions of an orthogonal or unitary transform, the pseudoinverse-based solution is equivalent to that obtained with a respective inverse transform from the incomplete data. Whereas for an overdetermined linear system, the pseudoinverse may be used to find the minimal mean-square-error solution, in the case of an underdetermined system it yields the solution with the minimal Euclidean norm. This kind of regularization is not necessarily what one needs in single-pixel detection. We refine this result by including minimization of a quadratic criterion into the solution. The criterion considered in this study combines the minimization of the image gradient with a penalty function for the high spatial frequency components of the image. It is expressed as an ℓ² norm of the convolution between the image and the filter proposed in this work. The final result is also a closed-form pseudoinverse-based method. It shares the advantage of the direct pseudoinverse based reconstruction in terms of its linear numerical cost as a function of the dimension of the measurement. It makes sense to precalculate the reconstruction matrix before the measurements take place and to store it in memory. This approach is feasible for a reasonably small dimension of the measurement vector, typically on the order of several thousands of elements. For a resolution of 256 × 256, 2000 elements correspond to a compression ratio of 3%. With a 22 kHz digital micromirror device (DMD), this allows for compressive single-pixel detection at a rate of over 11 Hz with image reconstruction obtained on the fly.

2. Image reconstruction algorithm

Single-step image reconstruction and a linear cost of the reconstruction algorithm as a function of the size of the compressive measurement are two factors desired for real-time single-pixel imaging. In our case, the reconstructed image is obtained as the solution to a minimization of a quadratic criterion with a linear constraint. The mathematical formulation of the reconstruction algorithm is presented in this section. In the text we call it the Fourier domain regularized inversion (FDRI) method. The usual single-pixel detection measurement model is assumed:

M \cdot x = y,

where x is the measured image consisting of n real-valued pixels, y is the result of measurement that includes k elements, and M is a rectangular k × n measurement matrix. Rows of M include the sampling patterns, and k < n for a compressive measurement. The image is reconstructed from the measurement y through the minimization of a quadratic criterion constrained by Eq. (1):

x_{0} = \underset{x}{\arg \min} E (x) subject to M \cdot x = y .

The proposed criterion is defined as the squared ℓ² norm of the (circular) convolution of x and a filter h:

E (x) = {‖ x * h ‖}^{2} .

Criterion E(x) measures the overall intensity of the result of a freely selected shift-invariant linear filtering applied to the image. It may be connected to various kinds of image properties. In fact, convolving the image and a filter is used for a broad range of purposes, including edge extraction, differentiation, averaging, noise removal, sampling, pattern recognition, or feature analysis. The same criterion (3) may be expressed in quadratic form:

E (x) = x^{T} \cdot C \cdot x = {\hat{x}}^{T} \cdot \hat{C} \cdot \hat{x},

where C is a circulant positive semi-definite square matrix. Its Fourier transform

\hat{C} = F \cdot C \cdot F^{*}

is diagonal and consists of non-negative elements

{\hat{C}}_{i, j} = δ_{i, j} | {\hat{h}}_{i} |^{2}

, where

\hat{h}

is the transfer function of the filter h. We use the caret to denote the Fourier transformed vectors and matrices; F is the respective 2D Fourier transform matrix. Moreover, we assume that

| {\hat{h}}_{i} |^{2}

is centrally symmetric (this is always true when the filter h is real-valued). Then, C is a symmetric real-valued matrix. Furthermore, we assume that C is invertible, which requires

| {\hat{h}}_{i} | > 0

. Under these assumptions, the solution of Eq. (2) may be obtained with Lagrange’s multipliers method, with

ℒ (x, λ) = \frac{1}{2} x^{T} \cdot C \cdot x - λ * \cdot (M \cdot x - y) .

Setting the derivatives of $\partial ℒ / \partial x_{i}$ to zero, together with Eq. (1), yields the following solution to the optimization problem (2):

x_{0} = P \cdot y,

with

P = C^{- 1} \cdot M * \cdot {(M \cdot C^{- 1} \cdot M *)}^{- 1}

where

C^{- 1} = F * \cdot {\hat{C}}^{- 1} \cdot F

and

{\hat{C}}_{i, j}^{- 1} = δ_{i, j} | {\hat{h}}_{i} |^{- 2}

. The result from Eq. (7) may be rewritten with the help of the Moore-Penrose matrix pseudoinverse A⁺ = A^* · (A · A^*)⁻¹. By introducing a new diagonal matrix

\hat{Γ} = {\hat{C}}^{- 1 / 2}

we obtain

P - F * \cdot \hat{Γ} \cdot F \cdot {(M \cdot F * \cdot \hat{Γ} \cdot F)}^{+},

where

{\hat{Γ}}_{i, j} = δ_{i, j} | {\hat{h}}_{i} |^{- 1}

. Either of the two formulas (7) or (8) may be used equivalently to calculate matrix P. The pseudoinverse is usually calculated with a singular value decomposition requires k² · n multiplications. On the other hand, matrix inversion by Gaussian elimination requires only k³ multiplications. However, for large matrices, the cost of (7) becomes dominated by matrix multiplication (with k² · n multiplications needed to calculate M · M^∗) and is similar to that of (8). Still, we have used formula (7) because in practice it was almost twice as fast as (8) for the calculations presented in this paper. This result is obviously affected by many factors, including the efficiency of parallel processing or memory bandwidth. We have also neglected here the cost of multiplication by F or C⁻¹, assuming the use of the fast Fourier transform algorithm. In any case, whichever of the formulas (7) or (8) is chosen, the calculation of P is computationally demanding. Fortunately, after P has been precalculated and stored in memory, x₀ may be reconstructed with Eq. (6) in only k · n multiplications, which is interesting for real-time image reconstruction. Because P is a large matrix, we store it in memory in a single-precision floating-point representation. A fixed-point representation probably could be used as well, contributing to further savings in memory and computation time.

The proposed reconstruction algorithm is non-adaptive. It is also not based on the ℓ¹ norm used in compressive sensing, so it does not tend to find a sparse representation of the image. Therefore, in a compressive measurement it is crucially important to use sampling, which by itself is an image compression method, rather than to rely on the incoherence [35] between the sampling and compression bases. For instance, the sampling functions may be chosen to have a similar power spectrum to the typical power spectrum of real-world images with dominant low-frequency information. Then, optimizing the criterion (3) may be seen as a way to regularize the result. Depending on the choice of the filter h, various features of the reconstructed image may be enhanced or suppressed, as long as this does not contradict the measurement (1).

In this study, we considered a simple generalization of the criterion (3) by defining a similar composite criterion. The modified criterion adds together the contributions from the norms obtained with different filters:

E (x) = \sum_{p} α^{(p)} {‖ x * h^{(p)} ‖}^{2},

where the coefficients α⁽^p⁾ >= 0 are responsible for weighting, in a trade-off sense, the convex criteria obtained with every filter h⁽^p⁾. With the modified criterion, we must adjust the definition of matrix C, which now becomes

C = \sum_{p} α^{(p)} C^{(p)}

with

{\hat{C}}_{i, j}^{(p)} = δ_{i, j} | {\hat{h}}_{i}^{(p)} |^{2}

. This yields the definitions of

{\hat{C}}^{- 1}

and

\hat{Γ}

necessary for Eqs. (7) and (8),

{\hat{Γ}}_{i, j} = \frac{δ_{i, j}}{\sqrt{\sum_{p} α^{(p)} \cdot {| {\hat{h}}_{i}^{(p)} |}^{2}}} and {\hat{C}}_{i, j}^{- 1} = {\hat{Γ}}_{i, j}^{2},

We further assume that the composite criterion (9) is constructed with the following filters: discrete gradient filters h⁽¹⁾ = ∂_x, h⁽²⁾ = ∂_y; a penalty filter for high spatial frequencies ${\hat{h}}^{(3)} = \sqrt{ω_{x}^{2} + ω_{y}^{2}}$ ; and an identity filter h⁽⁴⁾ = δ(x, y). Gradient filters, whose transfer functions are ${\hat{h}}^{(1)} = \sin (ω_{x})$ and ${\hat{h}}^{(2)} = \sin (ω_{y})$ , are responsible for minimizing image nonuniformity. The high-frequency penalty filter imparts a preference for a typical 1/|ω| image spectrum. Finally, h⁽⁴⁾ is needed to remove singularities from Eq. (10).

In the numerical part of this work, we show the image recovery results as a function of a single parameter μ ∈ [0, 1] (instead of α⁽^p⁾). With our four filters h⁽^p⁾, p = 1, ..., 4, we have the following form of $\hat{Γ}$ ,

{\hat{Γ}}_{i, j} = \frac{δ_{i . j}}{\sqrt{(1 - μ) (\sin^{2} (ω_{x}) + \sin^{2} (ω_{y})) + μ^{2} \frac{ω_{x}^{2} + ω_{y}^{2}}{2 π^{2}} + ε}},

where (ω_x, ω_y) is the spatial frequency corresponding to index i and ω_x_/_y ∈ (−π, +π). A small constant of 10⁻⁵ is assigned to ε. Sample spectra

\hat{Γ} (ω_{x}, ω_{y})

are plotted in Fig. 1 for various values of the parameter μ. We note that the proposed method converges to the pseudoinverse (P ≈ M⁺), when ε ≫ 1.

Fig. 1 Filter $\hat{Γ} (ω_{x}, ω_{y})$ defined in Eq. (11) for various values of the parameter μ (logarithmic scale).

Download Full Size | PDF

To summarize this section, we reconstruct the image x₀ from the measurement y using Eq. (6), where the matrix P is precalculated in advance with Eq. (7) or (8), and $\hat{Γ} (μ)$ is defined as (11). Straightforward pseudoinverse-based reconstruction, on the other hand, is based on Eq. (6) with the matrix P = M⁺.

3. Sampling protocols

The FDRI reconstruction method introduced in the previous section may be used with either orthogonal or nonorthogonal sampling. In this study, we considered three compressive sampling protocols, namely discrete cosine transform (DCT)-based sampling, Walsh-Hadamard sampling, and Morlet wavelet-based sampling, which all consist of small (3%-6%) subsets of the complete bases sets. We note that orthogonality is usually not preserved by binarization, which is needed for displaying the functions on DMD modulators. For instance, the DCT basis or Fourier basis binarized for use with DMD displays are no longer orthogonal. Sub-gridding, multiple exposure, or sophisticated binarization algorithms [27] have been used to approximate the continuous functions with binary patterns more accurately. However, as we show here, it may be sufficient to apply binarization directly to the continuous functions and to improve the reconstruction method instead. Other useful nonorthogonal bases interesting for single-pixel detection are based on nonorthogonal wavelet bases or on speckle patterns. As for the orthogonal and inherently binary Walsh-Hadamard basis, the proposed reconstruction method may improve the reconstruction quality at a low compression ratio k/n, compared to the direct reconstruction with the inverse transform or with a straightforward application of pseudoinverse (x ≈ M⁺ · y).

A novel kind of nonorthogonal sampling for single-pixel imaging has been recently proposed, with sampling functions obtained as a convolution between the real part of Morlet wavelets and realizations of white Gaussian noise [34]. We now briefly review its principles. The Morlet-based sampling proposed in [34] is equivalent to a uniform random sampling of the representation of the image in the feature space constructed with Morlet wavelets. A feature space is built out of vectors, whose elements correspond to specific features of images. Simple features may be associated with the spatial and frequency contents of an image. Here, the feature space is constructed using Morlet wavelets, and its part is preselected on the basis of the analysis of an image database. Morlet wavelets, also known as Gabor wavelets, consist of products of Gaussian functions with exponential functions. In two dimensions, a Morlet wavelet is equal to

g_{σ, n_{p,} θ} (x, y) = N e^{- \frac{x^{2} + y^{2}}{2 σ^{2}}} (e^{i (π n_{p} / 2 σ) (x \cos (θ) + y \sin (θ))} - κ),

where the constants κ and N assure that the wavelet function g is normalized ‖g‖ = 1 and has zero mean

\bar{g} = 0

. Parameters σ, n_p, θ are related to the size of the Gaussian envelope, number of periods within the envelope, and the orientation of modulation. The Morlet wavelet-based sampling functions considered in this study are obtained as convolutions of Morlet wavelets from Eq. (12) with realizations of white Gaussian noise. Similar to the DCT basis, the binarized Morlet-based sampling functions considered here are obtained from the continuous ones with a direct thresholding at the mean values of the respective continuous functions.

Morlet wavelets have well-established unique Fourier domain properties that determine their suitability for applications that require joint image description in the image and Fourier spaces. In fact, visual information is often analyzed in either or both of these two spaces. The Fourier transform of a Morlet wavelet is also a Morlet wavelet. Moreover, two Morlet wavelets constituting a Fourier function pair equalize the uncertainty principle [36]. Therefore, Morlet wavelets are a preferable choice for creating an image representation maintaining the best possible resolution tradeoff between image and Fourier spaces. In spite of the importance of Morlet wavelets in quantum mechanics with similarity to quantum coherent states, their role in human and artificial vision, and usefulness for multi-scale analysis shared with other wavelet representations, they also have drawbacks that limit their applications. In particular, they are in general not orthogonal; the expansion into Morlet wavelets is not unique and is computationally costlier than finding the Fourier representation, the kernel of which in fact consists of a particular orthogonal subset of Morlet wavelets.

We selected the sampling functions proportionally to the frequencies of their occurrences in a compressive decomposition of an image database. For the DCT and Walsh-Hadamard bases, a set of 49 images was decomposed with the respective two transforms. Then, the basis functions with the maximal average contribution to the expansions were selected. The same approach could not be used to select nonorthogonal functions obtained from Morlet wavelets. In this case, by decomposing a set of 49 images in the discrete Fourier transform (DFT) basis, we found the optimal range of values for the ratio $\frac{n_{p}}{σ}$ in Eq. (12). Then, we calculated n_p and σ using the following equation:

σ = σ_{m i n} + (σ_{m a x} - σ_{m i n}) \cdot \frac{ω_{m a x} - ω}{ω_{m a x} - ω_{m i n}},

where

ω = \frac{n_{p}}{2 σ}

. The largest values of σ are obtained at small spatial frequencies ω, and the smallest at the highest spatial frequencies.

The average spatial spectra of the sampling functions are compared in Fig. 2. These spectra are dominated by low-frequency contents, and binarization extends them towards higher spatial harmonics. We note that without weighting the sampling functions, these spectra do not need to match an average spectrum of the images, and they clearly differ from each other.

Fig. 2 Fourier spectra of sampling functions used in various sampling schemes $(\sum_{i} | {\hat{m}}_{i} | / ‖ {\hat{m}}_{i} ‖)$ , where m_i are the rows of the measurement matrix M, and the caret denotes the 2D Fourier transform. B denotes binarized sampling functions.

Download Full Size | PDF

Making a just comparison between different sampling schemes is beyond the scope of this paper. Because the proposed reconstruction method could be applied with any kind of sampling function, it is likely that the best sampling protocol for a specific application may be obtained by numerical optimization or machine learning, rather than by choosing one of the three methods considered here. In any case, it is crucial that the sampling functions should provide a reasonably good compression base for the measured images. The optimal choice of the sampling functions depends on the image formation model. For single-pixel imaging of typical real-world images, the sampling functions should be concentrated at low spatial frequencies. Although we proposed a simple and efficient methodology for selecting small subsets of the three bases, we are aware that there possibly exist other, better ways for doing so, which might potentially lead to better image reconstruction results.

4. Numerical results

In this section we summarize the simulations of computational imaging with the proposed FDRI reconstruction algorithm. In Fig. 3 we compare the FDRI reconstructions obtained with different values of the parameter μ introduced in Eq. (11) with a reconstruction based on the pseudoinverse of the measurement matrix. For μ = 0, FDRI minimizes the discrete gradient norm, whereas for μ ≈ 1 it tends to suppress the high spatial harmonics of the image. The presented test image was sampled with binary DCT functions at a compression ratio of 3%. In the case of the pseudoinverse-based reconstruction, the high level of undersampling and binarization of the sampling functions results in undesired artifacts and loss of image quality. With the FDRI method, the quality of the reconstructed image is visibly improved for all values of μ. Further, in this work, we assumed that μ = 0.5, which assures a balanced contribution of the filters to the composite criterion defined in Eq. (9).

Fig. 3 Single-pixel imaging at a compression ratio of 3%. Pseudoinverse and FDRI image reconstructions (top). The dependence of the reconstruction quality obtained with FDRI on the value of μ (bottom). Sampling functions consist of a subset of the binarized DCT basis. PSNR reaches maximal values for μ ∈ (0.5, 0.7).

Download Full Size | PDF

In Fig. 4, we compare the performance of the FDRI method when either continuous or binary sampling protocols are used. For continuous sampling functions composed of Morlet wavelets convolved with white noise or of DCT basis functions, reconstructions of comparably good quality are obtained with both pseudoinverse and FDRI recovery methods. The binarization of the sampling patterns modifies both their spatial and spectral properties. The quality of pseudoinverse image reconstructions is strongly affected by the binarization of the sampling functions. On the other hand, the quality of images reconstructed with the FDRI method remains similar for continuous and for binarized sampling functions. The FDRI method adjusts the Fourier spectrum of the solutions, and it is capable of diminishing the error introduced by using binary, instead of grayscale, sampling patterns.

Fig. 4 Image reconstructions obtained with the proposed FDRI (top row) and with the pseudoinverse (bottom row) methods for DCT and Morlet wavelet-based sampling schemes. When the sampling functions are continuous, FDRI and pseudoinverse reconstructions have a similar quality. On the other hand, for binary sampling, the FDRI method is much better than the pseudoinverse. The compression ratio is always equal to 3%.

Download Full Size | PDF

We also compared the image reconstructions obtained with the FDRI method and with the NESTA [37] package, which utilizes total variation (TV) minimization to reconstruct the image. TV optimization is one of the most common reconstruction methods used in single-pixel imaging. The results in Fig. 5 indicate that the FDRI method allows reconstruction of compressively measured images with a quality similar to that achieved with NESTA. In the presented example, for sampling with 3% binary DCT patterns, the peak signal-to-noise ratio (PSNR) for these two methods differs by only 0.2 dB in favor of NESTA. At the same time, the FDRI method is considerably faster. Because it requires only a single matrix by vector multiplication, the image reconstruction time is proportional to the number of sampling functions as well as to the number of pixels at the assumed resolution. For an image consisting of 256 × 256 pixels and at a compression ratio of 3%, image reconstruction takes 0.1 s on a medium-class PC (in the optical experiment we further decreased the computation time by using single-precision arithmetic in a C++ program). Reconstruction conducted with NESTA in the same conditions requires 7 s of computations, which precludes the use of methods based on TV minimization for real-time single-pixel imaging.

Fig. 5 Numerical simulations of compressive imaging, presenting a comparison between various reconstruction methods for DCT sampling functions at a 3% compression ratio (top). The PSNR values and reconstruction times as a function of the number of sampling functions (bottom). The proposed method (FDRI) achieves performance similar to that of total variation regularization (NESTA), while the calculation time is decreased by a factor of ~ 70.

Download Full Size | PDF

As a summary, a set of 49 test images were reconstructed with all three methods, namely FDRI, pseudoinverse, and NESTA, and which were analyzed in terms of PSNR for several different sampling protocols. The results are shown in Fig. 6. In most cases, the use of NESTA results in only marginally higher PSNR than those obtained with the FDRI method. When binarized sampling patterns are used, FDRI offers considerable improvement in reconstruction quality in comparison with pseudoinverse. For the binary Morlet sampling protocol, the FDRI method also outperforms the recoveries obtained with TV minimization. On the other hand, for continuous sampling protocols, PSNR values obtained with pseudoinverse and FDRI methods are similar, whereas a minor increase can still be obtained with NESTA. The highest qualities of the reconstructed images are obtained for continuous sampling methods. However, when only binary sampling functions are considered, the use of binarized DCT functions together with the FDRI or NESTA reconstruction methods offer the highest PSNR values. Only slightly worse results are obtained with Hadamard sampling.

Fig. 6 Comparison of peak signal-to-noise ratio (averaged over 49 test images) obtained for various sampling and reconstruction methods. The improvement gained by using the FDRI reconstruction method is especially visible for the binary sampling methods denoted with B.

Download Full Size | PDF

In real experimental conditions the measurement is always affected by noise. Figure 7 shows how the reconstruction quality (PSNR) depends on the level of additive noise for the FDRI, TV (NESTA) and pseudoinverse reconstructions. We assumed that the noise is uncorrelated and Gaussian. The compression ratio is equal to 3 %. Up to a certain level, all the methods are robust to the presence of noise. With further increases in noise intensity, the values of PSNR decrease monotonically in a similar way for the three methods.

Fig. 7 Robustness of the FDRI, pseudoinverse, and TV (NESTA) reconstruction methods to additive Gaussian noise for the three kinds of binary sampling at a compression ratio of 3%. The PSNR averaged over 49 test images is plotted as a function of σ_y/μ_y, where σ_y is the standard deviation of the Gaussian noise and μ_y is the mean value of the measurement vector y.

Download Full Size | PDF

5. Experimental results

We verified our results experimentally with an optical single-pixel camera setup, schematically illustrated in Fig. 8. The projection of a three-dimensional scene was sampled with binary patterns displayed by a DMD with 1024 × 768 pixel resolution at a sampling rate of 22 kHz. The signal was then measured using a differential photodetection technique [38,39], with two photodiodes simultaneously collecting light reflected from the areas of the DMD set to the on and off states, respectively. This approach, as compared to typical single-pixel cameras equipped with only one photodetector, significantly increases the experimental SNR by compensating the effects of background illumination, variations in light intensity over time, and signal oscillations induced by electronic equipment. The signals were sampled with a 15-bit digital oscilloscope at a rate of 7.8 MS/s, and the measurements were streamed on-the-fly to a PC, where the reconstruction of the video was conducted in parallel. Together with a cost-efficient image recovery method, in which only a single product of the precalculated matrix with a vector is needed per image reconstruction, the camera is capable of continuous sampling and reconstructing video images in real time, with a frame rate dependent only on the number of sampling patterns used in a compressive measurement. For instance, with resolution of 256 × 256 and 3% compression ratio, the video frame rate exceeds 11.3 Hz.

Fig. 8 Schematic view of the single-pixel camera. The DMD spatial light modulator samples the image with a set of binary patterns. A differential photodetection technique, with two photodetectors gathering the complementary optical signals reflected from the DMD in two directions, was used to increase the experimental SNR.

Download Full Size | PDF

The sampling functions used in the experiment were chosen in the same manner as discussed in Section 3; however, only binary functions were considered. Although grayscale patterns may also be achieved with DMDs by time-division multiplexing, only binary patterns are displayed at full speed. The resolution of the sampling functions was set to 256 × 256, while the patterns were resized to fill the whole area of the DMD display. In each set of sampling functions, a single white pattern (with all pixels set to the on state) was included for the purpose of synchronization between the DMD and the data acquisition unit. In Fig. 9 we present the reconstructions of a single representative video frame obtained with our camera using three different types of sampling functions at compression ratios of 3% and 6%. The reconstructions calculated using the proposed FDRI method are compared with recoveries based on the pseudoinverse of the measurement matrix. The improvement of the reconstruction quality with the FDRI method is especially robust for highly compressive measurements, when the Fourier domain regularization allows for smoothing the undesired artifacts resulting from the lossy compression accompanying the measurement. At the same time, the sharpness of the edges and information about small details in the image, such as the resolution chart, are preserved with high accuracy. Finally, in Fig. 10 we illustrate a selection of video frames recorded by our single-pixel camera using 3% of DCT basis functions for sampling the image at a frame rate of 11.3 Hz. Video frames were reconstructed on the fly using the FDRI algorithm with parameter μ = 0.5. Full video is included in Visualization 1.

Fig. 9 Comparison of reconstructions obtained with the FDRI and pseudoinverse methods for experimental data measured with an optical single-pixel camera with different types of sampling protocols at a resolution of 256 × 256 and compression ratios of 6% (top 6 images) and 3% (bottom 6 images).

Download Full Size | PDF

Fig. 10 Several representative video frames recorded with our optical single-pixel camera at a rate of 11.3 Hz and reconstructed on the fly using the proposed FDRI algorithm. The sampling patterns consist of 3% binarized DCT basis functions with a resolution 256 × 256. For full video see Visualization 1.

Download Full Size | PDF

6. Conclusions

We propose a fast, linear, closed-form reconstruction method for single-pixel real-time video imaging, experimentally demonstrating its feasibility to reconstruct good-quality video with 256 × 256 resolution at a frame rate of 11 Hz. Image reconstruction takes place on the fly in parallel with data acquisition. The reconstruction method is based on the multiplication of the generalized inverse of the rectangular measurement matrix by the measurement vector. This matrix is obtained by constrained minimization of a quadratic criterion. The criterion is defined using a combination of spatial filters corresponding to the minimization of discrete gradients of an image and a penalty filter for high spatial frequencies. The proposed method shares a major advantage with the direct pseudoinverse-based reconstruction. Namely, the recovery of an image requires only a single-step calculation with linear numerical cost as a function of the number of measurements. Therefore, the reconstruction is extremely fast and may be conducted on the fly with a single-pixel video camera equipped with a state-of-the-art DMD. For instance, in our experiment the image was sampled by a DMD at a rate of 22 kHz with binary patterns obtained from 3% of the discrete cosine transform basis. However, the reconstruction matrix must be calculated beforehand and stored in memory, which may be challenging for imaging at high resolution. Fortunately, it is possible to reduce the memory requirements by choosing a well-designed sampling protocol, which allows for highly compressive measurements. In our case, this is achieved by selecting the sampling functions proportionally to the frequency of their occurrence in the compressive decomposition of an image database.

The criterion included in the definition of the reconstruction procedure provides a better regularization of the inversion of the measurement matrix than the direct pseudoinverse. Therefore, it leads to a notably higher quality of the reconstructed images, with reduced noise and artifacts. We proved that the proposed reconstruction method is comparable with solvers on the basis of optimization of the total variation, such as NESTA, in terms of the PSNR of the recovered images, while the reconstruction time in the same conditions is reduced multiple times with our method. Moreover, the proposed reconstruction method remains valid for nonorthogonal sampling functions. Nonorthogonal functions are obtained, for instance, by binarizing continuous functions, such as DCT or Fourier basis, or Morlet wavelets convolved with white noise. Binarized sampling functions may be easily displayed on DMDs. However, they have different spectral properties from the original continuous functions. This causes some quality degradation in images reconstructed with pseudoinverse. We show that the proposed method is well-suited to mitigate this effect.

Funding

Narodowe Centrum Nauki (grant No. UMO-2014/15/B/ST7/03107).

References and links

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

2. J. H. Shapiro, “Computational ghost imaging,” Phys. Rev. A 78, 061802 (2008). [CrossRef]

3. L. Bian, J. Suo, G. Situ, Z. Li, J. Fan, F. Chen, and Q. Dai, “Multispectral imaging using a single bucket detector,” Sci. Rep. 6, 24752 (2016). [CrossRef] [PubMed]

4. Z. Li, J. Suo, X. Hu, C. Deng, J. Fan, and Q. Dai, “Effucient single-pixel multispectral imaging via non-mechanical spatio-spectral modulation,” Sci. Rep. 7, 41435 (2016). [CrossRef]

5. J. Senlin, H. Wangwei, W. Yunlong, H. Kaicheng, S. Qiushuai, Y. Cuifeng, L. Dongqi, Ye Qing, Z. Wenyuan, and T. Jianguo, “Hyperspectral imaging using the single-pixel Fourier transform technique,” Sci. Rep. 7, 45209 (2017). [CrossRef]

6. Q. Pian, R. Yao, N. Sinsuebphon, and X. Intes, “Compressive hyperspectral time-resolved wide-field fluorescence lifetime imaging,” Nat. Photon. 11, 411–414 (2017). [CrossRef]

7. V. Duran, P. Clemente, M. Fernandez-Alonso, E. Tajahuerce, and J. Lancis, “Single-pixel polarimetric imaging,” Opt. Lett. 37(5), 824–826 (2012). [CrossRef] [PubMed]

8. F. Soldevila, E. Irles, V. Durán, P. Clemente, M. Fernández-Alonso, E. Tajahuerce, and J. Lancis, “Single-pixel polarimetric imaging spectrometer by compressive sensing,” Appl. Phys. B 113(3), 551–558 (2013). [CrossRef]

9. J. Fade, E. Perrotin, and J. Bobin, “Polarizer-free two-pixel polarimetric camera by compressive sensing,” Appl. Opt. 57(7), B102–B113 (2018). [CrossRef] [PubMed]

10. B. Sun, M. P. Edgar, R. Bowman, L. E. Vittert, S. Welsh, A. Bowman, and M. J. Padgett, “3D computational imaging with single-pixel detectors,” Science. 340(6134), 844–847 (2013). [CrossRef] [PubMed]

11. M. Sun, M. Edgar, G. Gibson, N. R. B. Sun, R. Lamb, and M. Padget, “Single-pixel three-dimentional imaging with time-based depth resolution,” Nat. Commun. 7, 12010 (2016). [CrossRef]

12. L. Li, W. Xiao, and W. Jian, “Three-dimensional imaging reconstruction algorithm of gated-viewing laser imaging with compressive sensing,” Appl. Opt. 53(33), 7992–7997 (2014). [CrossRef]

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

14. V. Durán, F. Soldevila, E. Irles, P. Clemente, E. Tajahuerce, P. Andrés, and J. Lancis, “Compressive imaging in scattering media,” Opt. Express 23(11), 14424–14433 (2015). [CrossRef] [PubMed]

15. X. Liu, J. Shi, X. Wu, and G. Zeng, “Fast first-photon ghost imaging,” Sci. Rep. 8, 5012 (2018). [CrossRef] [PubMed]

16. D. J. Starling, I. Storer, and G. A. Howland, “Compressive sensing spectroscopy with a single pixel camera,” Appl. Opt. 55(19), 5198–5202 (2016). [CrossRef] [PubMed]

17. D. Pastor, A. Pastuszczak, M. Mikolajczyk, and R. Kotynski, “Compressive phase-only filtering at extreme compression rates,” Opt. Commun. 383, 446–452 (2017). [CrossRef]

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

19. J. Romberg, “Imaging via compressive sampling,” IEEE Signal Process. Mag. 25(2), 14–20 (2008). [CrossRef]

20. D. L. Donoho and Y. Tsaig, “Fast solution of ℓ₁ -norm minimization problems when the solution may be sparse,” IEEE Trans. Inf. Theory 54(11), 4789–4812 (2008). [CrossRef]

21. Y. C. Eldar, Sampling theory, Beyond Bandlimted Systems(Cambridge University, 2015).

22. 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]

23. I. Dokmanic, M. Kolundzija, and M. Vetterli, “Beyond Moore-Penrose: Sparse pseudoinverse,” in “2013 IEEE International Conference on Acoustics, Speech and Signal Processing,” (IEEE, 2013), pp. 6526–6530.

24. 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, 10669 (2015). [CrossRef] [PubMed]

25. G. M. Gibson, B. Sun, M. P. Edgar, D. B. Phillips, N. Hempler, G. T. Maker, G. P. A. Malcolm, and M. J. Padgett, “Real-time imaging of methane gas leaks using a single-pixel camera,” Opt. Express 25(4), 2998–3005 (2017). [CrossRef] [PubMed]

26. N. Huynh, E. Zhang, M. Betcke, S. Arridge, P. Beard, and B. Cox, “Single-pixel optical camera for video rate ultrasonic imaging,” Optica 3(1), 26–29 (2016). [CrossRef]

27. Z. Zhang, X. Wang, G. Zheng, and J. Zhong, “Fast Fourier single-pixel imaging via binary illumination,” Sci. Rep. 7, 12029 (2017). [CrossRef] [PubMed]

28. C. F. Higham, R. Murray-Smith, M. J. Padgett, and M. P. Edgar,“Deep learning for real-time single-pixel video,” Sci. Rep. 8, 2369 (2018). [CrossRef] [PubMed]

29. W. Yuwang, L. Yang, S. Jinli, S. Guohai, Q. Chang, and D. Qionghai, “High Speed Computational Ghost Imaging via Spatial Sweeping,” Sci. Rep. 7, 45325 (2017). [CrossRef]

30. Z.-H. Xu, W. Chen, J. Penuelas, M. Padgett, and M.-J. Sun, “1000 fps computational ghost imaging using led-based structured illumination,” Opt. Express 26(3), 2427–2434 (2018). [CrossRef] [PubMed]

31. C. Liu, J. Chen, J. Liu, and X. Han, “High frame-rate computational ghost imaging system using an optical fiber phased array and a low-pixel APD array,” Opt. Express 26(8), 10048–10064 (2018). [CrossRef] [PubMed]

32. C. Zhang, S. Guo, J. Cao, J. Guan, and F. Gao, “Object reconstitution using pseudo-inverse for ghost imaging,” Opt. Express 22(24), 30063–30073 (2014). [CrossRef]

33. W. Gong, “High-resolution pseudo-inverse ghost imaging,” Photonics Res. 3(5), 234 (2015). [CrossRef]

34. K. M. Czajkowski, A. Pastuszczak, and R. Kotynski, “Single-pixel imaging with Morlet wavelet correlated random patterns,” Sci. Rep. 8, 466 (2018). [CrossRef] [PubMed]

35. E. Candes and J. Romberg, “Sparsity and incoherence in compressive sampling,” Inverse Probl. 23, 969 (2007). [CrossRef]

36. J. G. Daugman, “Uncertainty relation for resolution in space, spatial frequency, and orientation optimized by two-dimensional visual cortical filters,” J. Opt. Soc. Am. A 2(7), 1160–1169 (1985). [CrossRef] [PubMed]

37. S. Becker, J. Bobin, and E. J. Candes, “NESTA: a fast and accurate first-order method for sparse recovery,” SIAM J. Imaging Sci. 4(1), 1–39 (2009). [CrossRef]

38. Y. Wen-Kai, L. Xue-Feng, Y. Xu-Ri, W. Chao, Z. Yun, and Z. Guang-Jie, “Complementary compressive imaging for the telescopic system,” Sci. Rep. 4, 5834 (2014).

39. F. Soldevila, P. Clemente, E. Tajahuerce, N. Uribe-Patarroyo, P. Andrés, and J. Lancis, “Computational imaging with a balanced detector,” Sci. Rep. 6, 29181 (2016). [CrossRef] [PubMed]

Real-time single-pixel video imaging with Fourier domain regularization

Abstract

1. Introduction

2. Image reconstruction algorithm

3. Sampling protocols

4. Numerical results

5. Experimental results

6. Conclusions

Funding

References and links

Supplementary Material (1)

Cited By

Figures (10)

Equations (13)

Optics Express