## Abstract

We present a closed-form image reconstruction method for single-pixel imaging based on the generalized inverse of the measurement matrix. Its numerical cost scales proportionally with the number of measured samples. Regularization of the inverse problem is obtained by minimizing the norms of the convolution between the reconstructed image and a set of spatial filters. The final reconstruction formula can be expressed in terms of matrix pseudoinverse. At high compression, this approach is an interesting alternative to the methods of compressive sensing based on l1-norm optimization, which are too slow for real-time applications. For instance, we demonstrate experimental single-pixel detection with real-time reconstruction obtained in parallel with measurement at a frame rate of 11 Hz for highly compressive measurements with a resolution of 256 × 256. To this end, we preselect the sampling functions to match the average spectrum obtained with an image database. The sampling functions are selected from the Walsh-Hadamard basis, from the discrete cosine basis, or from a subset of Morlet wavelets convolved with white noise. We show that by incorporating the quadratic criterion into the closed-form reconstruction formula, we can use binary rather than continuous sampling and reach similar reconstruction quality as is obtained by minimizing the total variation. This makes it possible to use cosine- or Morlet-based sampling with digital micromirror devices without advanced binarization methods.

© 2018 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

## 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 *ℓ*^{1}-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 *ℓ*^{2} 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:

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):

The proposed criterion is defined as the squared *ℓ*^{2} norm of the (circular) convolution of *x* and a filter *h*:

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:

*C*is a circulant positive semi-definite square matrix. Its Fourier transform $\widehat{C}=F\cdot C\cdot {F}^{\ast}$ is diagonal and consists of non-negative elements ${\widehat{C}}_{i,j}={\delta}_{i,j}|{\widehat{h}}_{i}{|}^{2}$ , where $\widehat{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 $|{\widehat{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 $\left|{\widehat{h}}_{i}\right|>0$. Under these assumptions, the solution of Eq. (2) may be obtained with Lagrange’s multipliers method, with

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

with where ${C}^{-1}=F*\cdot {\widehat{C}}^{-1}\cdot F$ and ${\widehat{C}}_{i,j}^{-1}={\delta}_{i,j}|{\widehat{h}}_{i}{|}^{-2}$. The result from Eq. (7) may be rewritten with the help of the Moore-Penrose matrix pseudoinverse*A*

^{+}=

*A*

^{*}· (A · A

^{*})

^{−1}. By introducing a new diagonal matrix $\widehat{\Gamma}={\widehat{C}}^{-1/2}$ we obtain where ${\widehat{\Gamma}}_{i,j}={\delta}_{i,j}|{\widehat{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*

^{2}·

*n*multiplications. On the other hand, matrix inversion by Gaussian elimination requires only

*k*

^{3}multiplications. However, for large matrices, the cost of (7) becomes dominated by matrix multiplication (with

*k*

^{2}·

*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*

^{−1}, 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*

_{0}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 *ℓ*^{1} 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:

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}{\alpha}^{(p)}{C}^{(p)}$ with ${\widehat{C}}_{i,j}^{(p)}={\delta}_{i,j}|{\widehat{h}}_{i}^{(p)}{|}^{2}$. This yields the definitions of ${\widehat{C}}^{-1}$ and $\widehat{\Gamma}$ necessary for Eqs. (7) and (8),

We further assume that the composite criterion (9) is constructed with the following filters: discrete gradient filters *h*^{(1)} = *∂ _{x}*,

*h*

^{(2)}=

*∂*; a penalty filter for high spatial frequencies ${\widehat{h}}^{(3)}=\sqrt{{\omega}_{x}^{2}+{\omega}_{y}^{2}}$; and an identity filter

_{y}*h*

^{(4)}=

*δ*(

*x, y*). Gradient filters, whose transfer functions are ${\widehat{h}}^{(1)}=\mathrm{sin}({\omega}_{x})$ and ${\widehat{h}}^{(2)}=\mathrm{sin}({\omega}_{y})$, are responsible for minimizing image nonuniformity. The high-frequency penalty filter imparts a preference for a typical 1/|

*ω*| image spectrum. Finally,

*h*

^{(4)}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 $\widehat{\Gamma}$,

*ω*) is the spatial frequency corresponding to index

_{x}, ω_{y}*i*and

*ω*

_{x}_{/}

*∈ (−*

_{y}*π*, +

*π*). A small constant of 10

^{−5}is assigned to

*ε*. Sample spectra $\widehat{\Gamma}({\omega}_{x},{\omega}_{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.

To summarize this section, we reconstruct the image *x*_{0} from the measurement *y* using Eq. (6), where the matrix *P* is precalculated in advance with Eq. (7) or (8), and $\widehat{\Gamma}(\mu )$ 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

*κ*and

*N*assure that the wavelet function

*g*is normalized ‖

*g*‖ = 1 and has zero mean $\overline{g}=0$. Parameters

*σ, n*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.

_{p}, θ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}}{\sigma}$ in Eq. (12). Then, we calculated *n _{p}* and

*σ*using the following equation:

*σ*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.

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).

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.

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.

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.

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.

## 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.

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.

## 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 ℓ_{1} -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]