Differential real-time single-pixel imaging with Fourier domain regularization -- applications to VIS-IR imaging and polarization imaging

The speed and quality of single-pixel imaging (SPI) are fundamentally limited by image modulation frequency and by the levels of optical noise and compression noise. In an approach to come close to these limits, we introduce a SPI technique, which is inherently differential, and comprises a novel way of measuring the zeroth spatial frequency of images and makes use of varied thresholding of sampling patterns. With the proposed sampling, the entropy of the detection signal is increased in comparison to standard SPI protocols. Image reconstruction is obtained with a single matrix-vector product so the cost of the reconstruction method scales proportionally with the number of measured samples. A differential operator is included in the reconstruction and following the method is based on finding the generalized inversion of the modified measurement matrix with regularization in the Fourier domain. We demonstrate $256 \times 256$ SPI at up to $17~$Hz at visible and near-infrared wavelength ranges using two polarization or spectral channels. A low bit-resolution data acquisition device with alternating-current-coupling can be used in the measurement indicating that the proposed method combines improved noise robustness with a differential removal of the direct current component of the signal.


Introduction
Indirect image measurement techniques called single-pixel imaging (SPI), since their introduction over a decade ago [1,2], have led to a considerable amount of novel ideas about image measurement at various wavelength ranges, spectral imaging, imaging through scattering media, 3D imaging etc. [3,4].
The speed of single-pixel imaging is limited by the frequency of spatial modulation and by the time needed for digital image reconstruction. Most real-time reconstruction approaches rely on a single-step image reconstruction using a fast transform, e.g. the Fourier (FFT), Walsh-Hadamard (FWHT), or discrete cosine (DCT) transforms, or on the evaluation of a matrix-vector product [5,6]. There is also growing interest in using neural networks for image reconstruction or for removing artifacts caused by compression [7,8].
Modulation speed obtained with modern digital micromirror devices (DMD), such as the one used by us, is on the order of tens of kilohertz. This is not a lot, as many thousand exposures are needed per single image measurement. Far higher frame-rates are achievable with structured illumination with LEDs arrays, although such set-ups for ghost imaging and SPI have been demonstrated only at low resolutions [9][10][11][12]. Modulation with rotary elements with fixed patterns is a high-speed cost-efficient alternative to using dynamic modulators in THz [13,14].
Block compressive imaging also effectively increases the sampling frequency per pixel with use of multiple detectors or a focal plane-array [15][16][17].
Differential photodetection techniques allow to eliminate some of the optical noise, and to effectively increase the SNR of the detection signal. Differential imaging and normalized imaging have been used in computational ghost imaging prior to SPI [18,19]. In SPI, a reference detector may be used either to measure the unmodulated reference signal and to normalize the detection signal in the presence of source intensity fluctuations, or to measure in parallel a differently modulated signal. The latter possibility is often based on the fact that DMD micromirrors may take exactly one of the two positions (let us call them "0" or "1"). Intensities of two reflected beams are proportional to the result of image sampling with complementary (binary negated) patterns [20,21]. Their difference provides a signal with improved SNR and with removed bias. This is possible when the collimated incident beam comes from a perpendicular direction onto the DMD, for instance in a telescopic set-up [20]. Otherwise the set-up must be carefully calibrated to retain a similar level of signal from both detectors, which becomes challenging for 3D moving and rotating objects. Still, this kind of balanced photodetection has been used by many groups [5,22,23]. A more straightforward way to obtain differential detection is to subtract the measurements from two subsequent measurements [24]. This is a way not only to remove a slowly changing bias from the signal but also to encode both positive and negative values of real-valued sampling functions [25][26][27][28]. A more elaborated encoding allows to represent complex-valued sampling functions with three [29] instead of four [30] real-valued non-negative patterns. This approach has its origins in interferometry and also allows to remove a constant-bias at the same time. Encoding a complex pattern with a larger number of non-negative intensity signals has been also proposed [31]. A simplex encoding may also be used to represent an arbitrary bunch of real-valued sampling functions using non-negative patterns whose number is increased by one. Once again the constant bias of the measurement is removed, and by varying the simplex dimension one may select an optimal trade-off between optical noise and compression noise, depending on the experimental conditions [32]. The same method may be easily combined with complementary detection to improve the SNR even further [32]. Encoding real or complex valued functions onto non-negative intensity patterns usually does not produce binary patterns yet, and a further binarization step is necessary. For instance error diffusion is commonly used with success thanks to the difference in DMD resolution and the actual imaging resolution.
In this paper we introduce a differential method for measuring the DC component of images (the zeroth spatial frequency) applicable to various SPI architectures with binary modulation. In a classical approach, the DC component is measured using a single pattern filled with all pixels in the state "1", possibly followed by another pattern with all pixels in the state "0" [20,21]. The proposed measurement method makes use of a small number of binary patterns with approximately (but not exactly) half of DMD pixels in the "1" and "0" states to encode the DC component. This allows to limit the required operating range of the detector and to improve the noise robustness of the measurement. Unlike in [31,32], the sampling is readily binary. At the same time it retains the advantageous of other differential methods such as the immunity of the measurement to some constant or slowly varying bias of the detection signal. On top of this, we propose a complete fast differential SPI framework with binary sampling and a single step image reconstruction. The set of sampling patterns includes several patterns required to capture the DC component of the image, and other patterns taken from some standard basis but binarized with variable thresholding. The purpose of varying the binarization level is to make the best use of the available detector measurement range. Image reconstruction is based on the Fourier domain regularized inversion (FDRI) [33] with a minor modification that includes the differential treatment of the measurement signal. We will use the short name D-FDRI for the modified method. In a simple comparison made under experimental conditions, the proposed reconstruction method gives a similar image quality as optimization with a compressive sensing method [34]. Finally, we validate D-FDRI experimentally. Concurrent imaging in the infrared and visible spectral bands, likewise polarisation imaging are fields where SPI is an interesting alternative to classical optical approaches [35][36][37][38][39][40][41][42][43]. Making use of rather standard single pixel camera set-ups we demonstrate the capability of conducting a real time continuous measurement and reconstruction at a resolution of 256 × 256 in two independent channels, using moderate computational resources.

Differential SPI with Fourier-domain regularized inversion
Compressive measurement of an image x (with pixel values arranged in a vector) in the presence of additive signal noise n s and detector noise n d may be expressed with a matrix-vector multiplication where the rows of the measurement matrix M contain the patterns displayed on the DMD during the measurement, and · is the dot product. In Eq. (1) noise has been decomposed into a signal independent part n d primarily attributed to the detector dark current and the signal dependent part n s from sources such as background illumination [44]. The measurement is compressive when the dimension of y which will be denoted as is smaller than the number of pixels in the image x. In our case = 256 2 and the compression ratio = / ≥ 2%. Consider applying discrete difference operator of the order p on y, Matrices D have the size of ( − ) × and are defined as where is the Kronecker delta symbol. Applying the gradient or Laplace operators on y makes y ′ invariant to the level of constant bias present in the detection signal represented in our model by the mean value ofn d . In the case of Laplace operator, a bias changing linearly with time is also removed. In practice, using y ′ for image reconstruction allows to operate the data acquisition device (DAQ) with alternatingcurrent (AC) coupling at a reduced voltage range with the bit-resolution better adapted to cover the signal range.
The generalized inverse of the matrix product D · M multiplied by y ′ gives an approximation for image x,x = (D · M) · (D · M) · x = P g · y, where P g = (D ·M) ·D and denotes the generalized matrix inverse,which is not unique. The most commonly used generalized matrix inverse operator is the Moore-Penrose (pseudo)inverse. Pseudoinverse is known to minimize the mean square error (MSE) in linear regression and is sometimes used in single-pixel imaging. However, in a compressive measurement, it may be more advantageous to use a different kind of matrix generalized inverse. For instance based on Fourier-domain regularization (FDRI) [5], the solution to Eq. (1) may be expressed as where + denotes the pseudoinverse, and F is the 2D Fourier transform.ˆ is a diagonal matrix where and are used to tune the properties of the regularization [5] and , are the spatial frequencies. Taking into account the Laplace or gradient operators from Eq. (2) which effectively modify the measurement matrix, the following matrix may be used to reconstruct the image in a differential measurement, Calculation of matrix P g is computationally costly, and in practice only possible at a small compression ratio = / . On the other hand, the measurement matrix M as well as P g can be precalculated before the measurement. We use the singular-value decomposition (SVD) to calculate the pseudoinverse in Eq. (6), and the FFT algorithm to calculate the 2D Fourier transforms. Then the image reconstruction requires us to evaluate a single matrix-vector product x = P g · y. Additionally, any negative values in the reconstructed image are replaced with zeros. The cost of matrix-vector multiplication is ( · ). For = 256 × 256-pixel images sampled at a compression ratio of = 2%, a real-time (17 Hz) reconstruction does not require using GPUaccelerated computation, even when two independent channels are present. In fact, we have developed a program in Python based on integer and floating-point single-precision arithmetic capable of receiving, interpreting and averaging data from the DAQ and reconstructing images in two bands in real time on an Intel i5-4590 CPU. Also at a higher compression ratio, e.g. when = 6% with a lower level of compression noise, the operating speed is limited by the DMD frame-rate and not by digital processing.
In comparison to reconstructions with inverse transforms, for which there exist fast algorithms such as FFT, DCT or FWHT, the proposed D-FDRI method is slower and more memory demanding. Its advantages are that it works well with arbitrary sampling matrices M, also with nonorthogonal and binary patterns. The computational cost of FFT is ( · log 2 ). On the other hand, that of the matrix-vector multiplication is ( · ) or ( · 2 ). Therefore in practice the method is applicable at small compression ratios . The reconstruction properties, including noise robustness, may be further optimized by the choice of and , and the reconstructed image is invariant to the presence of a constant or linearly changing in time bias in the detection signal.

Differential measurement of the DC-component with binary sampling functions
The usual method of measuring the zeroth spatial frequency (the DC-component) of the image is to display a pattern consisting of all pixels equal to "1" on the DMD, possibly followed by a pattern consisting of all pixels equal to "0". We have replaced a direct measurement with a measurement distributed over several sampling functions. We wanted them to be binary, with approximately but not exactly half of the pixels in states "0" and "1". On top of this, information about the DC component of the image should be preserved after applying the finite difference operator on the measurement with these patterns.
To this end, we group the DMD pixels into an odd number of classes . The class number is randomly assigned to every pixel. To encode the measurement of the zeroth spatial frequency of the image we define + binary sampling patterns. Each of these patterns consists of either ( + 1)/2 or ( − 1)/2 pixel classes in state "1". We note that such patterns contain approximately half of the pixels in states "0" and "1". To find a possible mapping between the pixel classes and the patterns we define a small binary matrix A of size ( + ) × . Each row of the matrix defines a single sampling pattern, and each column represents one class of pixels. We performed a numerical search for such binary matrices with the condition that every row of the matrix contains ( ± 1)/2 equal values (zeros or ones) and that (D + · A ) = . This assures that a row consisting of ones is linearly dependent with the rows of D + · A , implying that information about the DC component is contained in the measurement also after applying the Laplace or gradient operator D + . Among many possible solutions, we have chosen such matrices A that minimize the dispersion of coefficients needed to reconstruct the DC component from the measurement with D + · A . As an example, the following matrices satisfy these conditions for = 3 and = 7 and for the order of the finite difference operator = 1 and = 2, After finding A and mapping its columns to groups of DMD pixels, + additional binary sampling patterns are created. Figure 1 graphically illustrates how the binary sampling patterns are constructed with the help of matrix A . In Fig. 1, three colors indicate the three pixel classes. Every DMD pixel is randomly assigned to exactly one of the classes, and pixels within each class take the same value equal to the corresponding element of matrix A .   We use a subset of binarized DCT basis corresponding to the lowest spatial frequencies. These components are selected in a way shown in Fig. 2(a). The binarization level is selected randomly in the range [( − 1)/2 , ( + 1)/2 ]. By the binarization level we understand the percentage of pixels in state "0" after binarization. A sample DCT pattern binarized at various levels is shown in Fig. 3, while Fig. 4 illustrates the effect of randomizing the binarization level on the number of pixels in state "1" for the proposed sampling in comparison with a classical sampling protocol. SPI measurement of a uniform image would produce a detection signal y with the same distribution as that shown in Fig. 4. An image dominated with low-frequency contents is likely to give a similar signal as well, and the binarization range gives a rough indication of the operating range of the detector. Fig. 2(c) shows that binarization of the real-valued sampling patterns enhances the high spatial frequency contents of their power spectral density (PSD). On the other hand, Fig. 2(c)(d) indicate that the PSD is little affected by the finite difference operators. The effective measurement matrix D · M consists of linear combinations of patterns in M, and in terms of PSD, the differential sampling is nearly equivalent to the original one consisting of binarized DCT patterns.

Tuning the method
The properties of the D-FDRI framework presented in this section depend on the choice of several parameters, namely , , and .
The value of determines the number of additional samples that probe the DC component of an image and the thresholding range used in the binarization of the sampling functions. We used = 7 most of the time, which implies that the binarization level of sampling functions is selected randomly in the range [43%, 57%], although values of between 3 and 11 also work well. As a rule of thumb, the thresholding level should be close to 50%, but varying it, increases the entropy of the measured signal, and improves noise robustness of the method.
The order of the differential operator ∈ {1, 2} is not as important as we expected. We did not observe any practical advantage of using = 2 over = 1 in the presence of noise under experimental conditions. When = 2 the calculation of pseudoinverse in Eq. (6) usually resulted in obtaining some very small eigenvalues in the SVD. Still, in practice we reached a similar imaging quality with both operators. As a remark, Eq. (6) should be evaluated using double-precision floating point arithmetic, even if P g is later stored with a lower precision for the reconstruction stage.
In Fig. 5 we demonstrate how the parameters and influence the performance of the D-FDRI method without noise and in the presence of measurement noise. To measure the quality of the reconstructed images, we use the peak signal-to-noise ratio (PSNR) metric defined as: where max(x) is the maximal pixel brightness of the original image. PSNR is affected both by the measurement noise and compression noise. Without measurement noise PSNR primarily depends on the compression CR. However at a certain level of experimental noise, PSNR can no longer be increased by increasing CR. This situation is illustrated in Fig. 5, which shows the drop in PSNR due to the measurement noise for various values of and . Noise robustness improves with the increase of + . In fact, PSNR explicitly depends on the mean value of the image, which is found from the measurement with + patterns. The more patterns are used to determine the DC component, the lower is its measurement uncertainty.  The significance of parameters and in Eq. (5) has been discussed in Ref. [5]. The two parameters require fine-tuning. While ≈ 0.5 works reasonably well in most situations [5], the value of may significantly affect noise robustness of the method. It is reasonable to take << 1 and then to test by trial and error the value of and respective magnitude of the elements in matrix P g . In the case |P g | contains elements with a large magnitude, image reconstruction with the formulax = P g · y would enhance the measurement noise.

Complementary sampling
SPI optical set-ups with complementary sampling make use of the two beams reflected from the DMD. A single measurement gives two measurement vectors, y 1 = M · x and y 2 = (1 − M) · x (where we have neglected the noise terms, and 1 stands for a × matrix filled with ones). The usual approach is to take the difference y = (y 1 −y 2 )/2 as the result of a differential measurement. This allows for increasing the SNR and removing an equal bias signal from the two detectors.
Since in D-FDRI the differential operator is included in the definition of P g in Eq. (6), the reconstructed image does not change when some constant value is added to y P g · y = P g · (y + c).
This makes it possible to acquire data from the detectors with AC coupling. On the other hand, it allows to construct two independent data channels with the two measurements y 1 and −y 2 . The same reconstruction matrix P g may be used in the two channels. The final formula for image reconstruction in the two channels is thereforẽ where = 1, 2 is the channel number, and ( ) = if > 0, and 0 otherwise. In the present work, images reconstructed in the two channels are represented with different colors.

Imaging quality, noise robustness, and signal entropy
Entropy of the signals measured by the SPI system is calculated as where ( ) is the probability of an outcome of a single measurement ( = 1, 2, ..., ) taking value . The total number of possible outcomes of a single measurement equals 2 , where is the number of bits used in the analogue-to-digital (A/D) conversion of the signal. In the information theory, entropy is a common measure of the average information carried by the possible outcomes of a stochastic variable. The less probable a specific outcome is, the more information is conveyed in its occurrence. While entropy is limited by the number of possible states of the random event ≤ , its maximal value is attained in case of events with uniform probability distribution ( ). It may be intuitively rephrased that the outcome of an experiment is most uncertain, and therefore most informative, when all possible results are equally probable. Entropy is also an essential measure for data compression, as it represents the absolute limit of how well data may be compressed without loss. In Fig. 6(a) we compare our proposed differential SPI sampling protocol with the standard non-differential one in terms of the average entropy per measurement / . The respective probability distributions ( ) for both sampling methods are presented in Fig. 6(b). We estimate the probability distributions ( ) experimentally, based on the histograms of simulated measurements calculated for over 180 thousand images. The images are taken from the Caltech 256 Image Dataset [45], which contains 30 thousand images organized into 256 distinct categories. Additional data augmentation (rescaling, random cropping, flipping, rotation) was used to further increase the amount of data. For each of the images, we calculate the SPI signal y obtained by sampling the image with patterns stored in matrix M, according to Eq. (1) excluding the noise. Then, each signal y is normalized. For the proposed sampling, we assume that the mean value of the measurementsȳ = 0 (i.e. only the AC component of the signal is measured) and the signal fits in range −1 ≤ ≤ 1 for = 1, 2, .., . For the classical sampling, all measurements are non-negative and, after normalization, they fit into the range 0 ≤ ≤ 1, where the maximal value 1 is practically obtained only for the measurement with the pattern which has all pixels in the "1" state. Finally, the signals are discretised using between 8 and 16 bits as , which corresponds to the range of the bit-resolution of the DAQ. The estimates of ( ) are calculated for each value of separately as the average probability of measuring value over all of the sampling patterns and the whole image dataset. While the standard SPI protocol results in the concentration of the measured values around the center of the DAQ range, the proposed differential method yields much broader spread of the measured signal. This broadening is achieved thanks to the following properties of the proposed method: (1) the elimination of the constant bias from the measurements by introducing the differential measurement protocol, (2) decomposition of the DC measurement into several patterns with transparency close to (but not equal) 50%, and (3) the application of randomized thresholding to the DCT functions. The resulting entropy of the proposed sampling is close to its upper limit and takes an average value / ≈ 0.96 , while for the standard sampling protocol / ≈ 0.4 on average.
For the purpose of numerical evaluation of the proposed method, in the present and the following simulations we have taken = 7, = 2, = 0.5 and = 10 −8 . The size of the scene is = 256 2 , and the number of DCT sampling functions is equal to = 1974 + + , which corresponds to the compression ratio ≈ 3%. DCT basis corresponding to the lowest spatial frequencies are selected, and they are thresholded randomly to contain between 43% and 57% pixels in the "1" and "0" states. As a reference, we take the non-differential FDRI [5] with the same basis functions discretized at 50%, and with = 0.5. This way, the comparison is focused on the significance of including the differential operator and on the varied thresholding in D-FDRI. In Figs. 7 and 8 we analyze the robustness of the proposed method to noise coming from two mutually independent sources: the statistical detector noise n d , modeled as additive white Gaussian noise, and the discretization noise applied to the signal at the stage of A/D conversion. For this reason, we have performed simulated measurements and reconstructions with a set of  Fig. 8. Sensitivity of the proposed sampling method (red) to the A/D discretization. Standard non-differential SPI sampling with the same set of DCT patterns (blue) for comparison. Average PSNR (a) and SSIM (b) values obtained from simulated SPI measurements and reconstructions with 5120 images. Before the image reconstruction, the SPI signals are converted into signals with b significant bits, i.e. only 2 distinct values are measured by the detector. We use b between 8 and 16, which corresponds to the range of the bit-resolution of our experimental equipment. We assume no presence of other types of experimental noise. 5120 images selected randomly from the Caltech 256 database. For evaluation of the quality of the reconstructed images, we use two metrics, namely: the PSNR defined in Eq. (9) and the mean similarity index (SSIM).
For calculating SSIM, we use the built-in Matlab function ssim with default parameters, which implements the mean SSIM metric from [46].
In contrast to the standard SPI protocol, D-FDRI is significantly less sensitive to both types of noise present in the SPI cameras. Especially its robustness to the discretization noise allows to lower the cost of the SPI set-up by taking advantage of a cheaper, lower-resolution DAQ with almost no impact on the quality of the recorded images. We have found that in the case of using an 8-bit DAQ, the average decrease of PSNR due to the discretization is ∼ 0.6 dB for the proposed method, as compared to almost 2.5 dB for the standard non-differential sampling.

Experimental polarization imaging and VIS-IR imaging
In this section we validate the proposed D-FDRI SPI framework experimentally. The schematics of our optical SPI set-ups for polarization imaging and imaging in the visible and near infrared wavelength bands are shown in Fig. 9(a) and (b), respectively. The modulator splits the incoming optical beam into two complementarily modulated beams that are measured independently with two bucket detectors. In the case of polarization imaging, linear orthogonally oriented polarizers are put in front of the two photodiodes, and after image reconstruction, the polarization information is indicated with pseudocoloring. In the presented situation, the scene is illuminated with unpolarized light which passes through a linear polarizer and a moving birefringent object. In the case of VIS-IR imaging, two different detectors are used for measuring the intensities of the beams reflected from the DMD. A similar set-up could allow for a complementary differential measurement, however here the two channels are used to measure independently two polarization states or two wavelength bands, respectively, while the sampling and reconstruction method assure the differential properties of the measurement.
The DMD (Vialux V-7001 XGA with DLP7000 chip) is operated at 22.7 kHz and at a spatial resolution of 256 × 256 (with 3 × 4 pixels blocks assigned equal values). The infrared channel makes use of reflective optics. The detector (Vigo System, PVI-4TE-1x1mm) is a four-stage thermoelectrically cooled IR photovoltaic detector based on HgCdTe heterostructures optimised   Fig. 10. Flow-chart explaining the processing of acquired data. Data is acquired continuously at a sampling of 256ns. The data gathered during the exposure of a single sampling pattern (marked in color) are selected and averaged to obtain a single measurement . Averaging is performed independently in each of the two channels. Then, monochromatic images from the two channels are reconstructed independently and integrated by color multiplexing to produce a single pseudocolor frame.
for the range 1-5.5 m. The detector is equipped with integrated transimpedance amplifier (DC 100k V/A). For the visible range, we use Thorlabs PDA100A2 amplified silicon photodiodes. The signals are digitized with a Picoscope 5000 DAQ with a sampling of 256ns and streamed through a USB bus for real-time processing. A multiprocess Python program reads, synchronizes, and averages the data, and reconstructs images from two independent channels in real time. The flowchart in Fig. 10 shows the data processing stages in the two channels and image integration into a single pseudocolor frame. For polarization imaging, the two polarization channels provide the red and blue components of the red-green-blue color representation. For visible-infrared imaging, the visible channel is represented as a gray-scale image, and is replaced with the infrared channel represented with an orange image at those locations where the infrared channel exceeds a threshold value that corresponds to a temperature of approximately 100 • C. Sample image reconstructions are shown in Fig. 11. The following parameters were used  The compression ratios between 2% to 6% allow for imaging at in between 17 Hz and 6 Hz. For still objects, image quality naturally increases with the compression ratio / which is Table 1. Performance of D-FDRI reconstruction method in comparison with the compressive sensing (CS) techniques in terms of the average reconstruction time and PSNR obtained for compression ratios = 2%, 3%, or 6%. In both cases, the same measurement matrices have been used for each value of , and the measurements have been obtained using a single-channel visible SPI for the still robotic toy scene presented in Fig. 12(a). CS reconstructions have been calculated via solving the basis pursuit denoise problem with the total variation regularisation using NESTA solver [34] inversely proportional to the framerate. However, for non-stationary objects the value of PSNR is a trade-off between compression ratio and acquisition time. This trade-off is characterized in Fig. 12 showing a rotating Lego robot. Fig. 12(a) and (b) show the SPI measurements of the robot at the compression ratio equal to 100%, 6%, 3% and 2%, respectively. Rotating wheels give strong artifacts shown in Fig. 12(c), and the PSNR drops the most rapidly with angular velocity when the compression ratio is the largest (See Fig 12(d)). These results are also shown in Visualizations 4-6, made for a varying angular speed of robot wheels. Finally, in Table. 1 we compare the proposed D-FDRI reconstruction method with the classical compressive sensing (CS) approach in terms of both image quality and reconstruction time. In both cases, the same differential binarized DCT measurement matrices are used with = 7 and = 1. For the purpose of CS reconstruction, we solve the basis pursuit denoise problem with the total variation regularisation using NESTA solver [34]. As the differential measurement matrix D p k · M is not semi-orthogonal (i.e. (D p k · M) · (D p k · M) ≠ I, where I is an identity matrix), we reformulate the CS optimization problem into the following equivalent form: where D p k · M = U · · V T is the SVD decomposition of matrix D p k · M. The parameter in Eq. (13) controls the trade-off between sparsity and fidelity of the reconstruction. Denotations v 2 and v TV stand for the ℓ 2 -norm and total-variation norm of a vector v, respectively. The reconstruction times presented in Table. 1 do not include either the time necessary to compute SVD for CS reconstruction or to compute matrix P g for D-FDRI, as they both may be computed only once for each measurement matrix and stored. The reconstruction times have been obtained on an Intel i5-4590 CPU without any GPU acceleration. For the given sampling, D-FDRI is faster than TV-norm minimization by over two orders of magnitude. As a matter of fact, being a single-step reconstruction method based on calculating a single matrix-vector product per image, D-FDRI is fast enough to reconstruct video images with resolution 256 × 256 in real-time on a desktop CPU using over 20 kHz DMD for image sampling. At the same time, both methods produce reconstructions of comparable quality, while D-FDRI yields on average 0.23 dB higher PSNR owing to better robustness to noise. This result does not prove that the proposed method would always give a higher PSNR than the methods of compressive sensing but it indicates that both approaches typically provide comparable results.

Conclusions
We propose a fast, differential, single-pixel imaging framework and demonstrate real-time SPI at the resolution of 256 × 256 at frame-rates between 6 and 17 Hz in two independent channels assigned to different spectral bands or polarizations. The D-FDRI reconstruction method is based on a matrix-vector multiplication (one multiplication per channel per frame) and resembles the solution proposed in Ref. [5] but is modified to account for the differential treatment of the measured data with either a Laplace or gradient discrete operators. It gives a similar or even better image quality as could be obtained with compressive sensing,but is significantly faster than iterative optimization techniques. Measurement of the zeroth spatial frequency is decomposed into several measurements with binary patterns with help of an auxiliary binary numerically optimized matrix. D-FDRI could use various bases of patterns, and we focus on binarized DCT basis with varied thresholding. The entropy of the detection signal is increased in comparison to standard SPI protocols enabling to reduce the bit-resolution requirements of the DAQ or otherwise improving noise robustness of the proposed SPI method. A trade-off behaviour between imaging quality and compression ratio is demonstrated as a function of the frame-rate.

Disclosures. The authors declare no conflicts of interest.
Data Availability Statement. Source data will be provided by the authors at a reasonable request. Source code for calculating measurement and reconstruction matrices will be made public at [47] after article publication. The FDRI part of the code is based on [48]. We have used the Caltech256 image database [45] in this work.