Single pixel imaging at high pixel resolutions

The usually reported pixel resolution of single pixel imaging (SPI) varies between $32 \times 32$ and $256 \times 256$ pixels falling far below imaging standards with classical methods. Low resolution results from the trade-off between the acceptable compression ratio, the limited DMD modulation frequency, and reasonable reconstruction time, and has not improved significantly during the decade of intensive research on SPI. In this paper we show that image measurement at the full resolution of the DMD, which lasts only a fraction of a second, is possible for sparse images or in a situation when the field of view is limited but is a priori unknown. We propose the sampling and reconstruction strategies that enable us to reconstruct sparse images at the resolution of $1024 \times 768$ within the time of $0.3~$s. Non-sparse images are reconstructed with less details. The compression ratio is on the order of $0.4 \%$ which corresponds to an acquisition frequency of $7~$Hz. Sampling is differential, binary, and non-adaptive, and includes information on multiple partitioning of the image which later allows us to determine the actual field of view. Reconstruction is based on the differential Fourier domain regularized inversion (D-FDRI). The proposed SPI framework is an alternative to both adaptive SPI, which is challenging to implement in real time, and to classical compressive sensing image recovery methods, which are very slow at high resolutions.


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]. Digital micromirror devices (DMD) are the most frequently used spatial light modulators in optical SPI set-ups and in this paper we will take into consideration their typical technical specifications such as resolution, modulation frequency, binary operation, and bandwidth. Modulation frequency of modern DMDs is on the order of twenty kilohertz. This is not a lot, as at least several thousand exposures are needed per single image measurement. A non-compressive sequential measurement at a resolution of 1024 × 768 would involve projecting patterns that take 77GB at the bitrate of 17.7Gb/s. As a result, either the image acquisition time must be very long or the resolution is reduced to the commonly reported range between 32×32 and 256×256 falling far below imaging standards with classical methods. Additionally, for compressive imaging the time needed for digital image reconstruction may be substantial at higher resolutions. Most real-time reconstruction SPI 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][9].
There exist alternatives to the DMD technology, which we will not consider in this paper. Far higher frame-rates are possible with structured illumination with LEDs arrays, although such setups for ghost imaging and SPI have been demonstrated only at low resolutions [10][11][12][13]. Modulation with rotary elements with fixed patterns is a high-speed cost-efficient alternative to using dynamic modulators in THz [14][15][16]. The modulation speed may be also increased by combining various light modulation techniques together or by using arrayed light sources with a fast modulation rate [17,18]. High resolution images may be obtained by data fusion techniques when SPI is combined with high resolution images acquired with classical cameras [19,20]. Block compressive imaging with use of multiple detectors or with a focal plane-array, as well as parallel detection with a pushframe camera also effectively increase the sampling frequency per pixel [21][22][23][24].
To reach the resolutions above 128 × 128 one has to accept strong compression (i.e. a low compression ratio) which usually results in a poor image quality. One way of retaining a good image quality is to apply adaptive sampling [25][26][27][28] with a sequence of sampling patterns selected dynamically during the measurement. Having in mind the large bitrate at which subsequent patterns should be calculated and sent to the DMD, adaptive sampling is difficult to implement in real time and currently has a rather theoretical significance.
In this paper we make an attempt to construct an SPI framework optimal in terms of using DMD modulation at high resolutions. We will consider non-adaptive high resolution binary sampling and a low compression ratio. Similarly as in our recent work [29] we will also assume that all sampling patterns contain approximately but not exactly half pixels in the on-state and that the measurement is differential. This provides an increased signal entropy and an improved signal-to-noise-ratio (SNR). We will also commence the image reconstruction with the differential Fourier domain (D-FDRI) regularization method proposed in [29]. The novel elements proposed in this paper are the sampling patterns based on multiple image maps and the second stage of image reconstruction that makes use of these maps. Every map defines a distinct partitioning of all image pixels into nonoverlapping sectors. In a noise-free scenario, the initial image reconstruction is guaranteed to provide the correct mean values for each sector of every map. The purpose of the second stage of image reconstruction is to determine the actual field of view based on the locations of empty sectors of every map. Non-empty sectors are corrected accordingly. Overall, the initial image reconstruction provides a low-quality image which will be the final result for dense images sampled at a low compression ratio. The second stage of the algorithm improves the reconstructed image, potentially up to the resolution of the DMD, when the image is sparse.

Objectives
Our aim is to introduce a sampling and image reconstruction framework for SPI that would satisfy a number of practically driven requirements. We want to make efficient use of the full DMD's spatial and frequency bandwidths for image sampling. For modern DMDs this implies binary sampling at the resolution on the order of 1024 × 768 and at the modulation frequency of 22.5 kHz. This corresponds to the transmission bandwidth of 17.7 Gb/s. At such a high bandwidth, it is difficult to implement adaptive sampling. Therefore, our choice is to use binary non-adaptive sampling at a full DMD resolution. To make the measurement practically feasible, the acquisition time should not exceed a fraction of a second. This implies an extremely strong compression -for instance the compression ratio of 0.5% corresponds to the acquisition rate of 5.6 Hz, which is still reasonable.
We want to use differential sampling. This is a simple technique to improve the SNR and to make the measurement independent of a constant or slowly varying background bias signal. At the same time, we do not want to compromise the available DMD bandwidth. We will use the differences between measurements with subsequent binary sampling patterns to reconstruct the image. This means that effectively we increase the number of measurements by one, rather than by a factor of 1/2, 1/3 or 1/ with ∈ N that is often accepted with more straightforward differential methods [30][31][32].
We want to be able to determine the field of view composed of possibly non adjacent non-empty areas of the image easily. With this aim, we propose a novel differential sampling scheme with binary sampling patterns that will serve a double role. First, the sampling should give the usual information about the spatial contents of the measured image. On top of this it should give some guarantees about the possibility of identifying empty regions of the image. We will consider multiple partitionings of the image surface with the help of auxiliary maps, and then construct the actual binary sampling patterns using these maps. By an image map we understand any partitioning of all image pixels into distinct subsets. These subsets may consist of both isolated as well as of neighbouring pixels. Sampling functions deduced from a single map should enable us to find the mean values of the measured image within every subset of the pixels in that map. Moreover, we want to keep approximately half of the pixels in the on-state in every sampling pattern.
Finally, we want to have a fast image reconstruction algorithm, capable of reconstructing the image at a quality dependent on image sparsity or on the field of view. When objects cover the whole surface of the image, the compression ratio is low (on the order of 0.5% or less) and one can at most expect a low-quality reconstruction. When the field of view is limited, the quality may be improved. The challenge is that the first situation requires low spatial frequency sampling, while the later, a high frequency sampling. In effect, neither random sampling with white noise, nor low-frequency e.g. Fourier or DCT sampling perform well. We propose to use image maps based on somehow arbitrarily chosen realizations of spatially correlated Gaussian noise to get a desired trade-off between low and high frequency sampling, and to construct the maps. The reconstruction algorithm begins with a non-negative differential Fourier-domain regularized solutions (D-FDRI) [29]. The advantages of the D-FDRI are that its implementation requires just one single matrix-vector multiplication, it is applicable with high resolution binary differential sampling, and it provides a regularized solution to the inverse problem. Following, the result is iteratively improved by assigning the value of zero to the determined areas of the image maps, and by scaling other pixels to repair the mean values within other regions. This procedure is iterated for all maps.
2.2. Background on the differential Fourier Domain Regularized Inversion image reconstruction method (D-FDRI) [29,33] Compressive measurement y 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 as where the rows of the measurement matrix M contain the patterns displayed on the DMD during the measurement and · denotes 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 [34]. 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, the sampling is binary and makes use of all the DMD pixels. Therefore, M is a binary matrix, = 1024 · 768, and is on the order of 3 · 10 3 . A linear reconstruction followed by the truncation of negative values (here denoted with a function) takes the formx = (x 0 ), wherex 0 = P · y.
Depending on the choice of the measurement matrix M and the reconstruction matrix P, one can obtain various SPI schemes. A frequent choice is to take M and P as consisting of selected rows of a linear transform matrix and selected columns of its inverse eg. of a direct and inverse Fourier, DCT or Walsh-Hadamard transforms. In a more general approach, P may be the Moore-Penrose pseudoinverse + of M (i.e. P = M + ). The pseudoinverse approach allows for using an arbitrary form of sampling patterns but without regularization may result in a poor noise robustness. Even more generally, P may be any matrix that is a generalized inverse of M, (i.e. P = M ). The generalized inverse is not unique and this approach allows for including some optimization or regularization in the choice of P. A generalized inverse is used in the FDRI method [5] to regularize the inverse problem in the Fourier domain. Finally, a regularization may be combined with a differential measurement and then P = (D · M) · D, where D is the 1D finite-difference operator (discrete gradient , = , −1 − , , where is the Kronecker delta). This is what we have proposed in [29] as the differential Fourier domain regularized inversion (D-FDRI) method, which is the starting point for the present work. We note that it isx 0 and notx that fulfills the measurement equation (1) in the absence of noise. Still in practice, in terms of the peak signal to noise ratio (PSNR), the nonnegative imagẽ x is a better approximation to the object x. For D-FDRI, the reconstruction matrix P is calculated from M as [29], where F is the 2D Fourier transform, * is the complex conjugate transpose, and the filterˆis a diagonal matrix [5] where and are used to tune the properties of the regularization and , are the spatial frequencies. In the present paper we assume that = 0.5 and = 10 −7 . D-FDRI reconstruction may be combined with arbitrary measurement matrices, for instance with binary nonorthogonal patterns. D-FDRI is inherently differential, so the reconstruction is blind to the mean value of the noise n d . It could be also directly applied with complementary sampling schemes if the SNR needs to be further improved [29]. Calculation of the reconstruction matrix P is computationally and memory intense but is done once. Afterwards, the image is reconstructed from the measurement y with use of Eq. (2) in single matrix-vector multiplication.

Image maps
Our measurement matrix M consisting of rows with binary sampling patterns is calculated based on a set of image maps m 1 , m 2 , ...m . Every map consists of the same number of pixels as the DMD and is defined by labeling the pixels with integer numbers between 1 and , i.e. m ∈ {1, ... } , with ∈ {1, } . In this paper the number of maps is = 100, and = 31 but other odd values of are also possible. Our maps are arbitrarily defined with the only assumption that they should be dominated with low spatial frequency contents but still contain some high spatial frequency information as well. A more extensive study of possible map compositions is beyond the scope of this paper but certainly some further optimization of the maps for a given kind of images is possible. Here the maps are obtained by generating spatially correlated Gaussian complex-valued noise and by assigning the label ∈ {1, } to each pixel based on the uniform quantization of the phase level of a noise realization. Some samples of the maps are illustrated in Fig.1(a). The resolution of these maps is 1024 × 768 and each of them defines an image partitioning into = 31 pixel subsets. Then, each map is translated into a sequence of binary patterns which are stacked together to form the measurement matrix M. Fig. 1(b) illustrates some of the patterns created using the map from Fig. 1(a).
Translation of maps into binary patterns is graphically explained in Fig. 2. Translation is performed with the help of subsequent rows of an auxiliary binary matrix A which play the role of lookup tables. To create a single sampling pattern we select one row from A and replace the  The binary sampling pattern obtained from the map in subplot (a) using the row of matrix A marked with a blue dashed rectangle. Pattern encoding is repeated for subsequent rows of A and for subsequent maps. The first two rows of A are used only for the first map and are omitted for other maps to remove linear dependence from the set of sampling patterns. Overall this gives a sequence of = ( + 1) + ( − 1) · ( − 1) = 3002 binary differential 1024 × 768 sampling patterns that form the rows of the measurement matrix M.
-th region in the map (with ∈ {1, .. }) with the binary value taken from the -th column of the lookup table. Then we repeat the same using subsequent rows from A. Matrix A is shown in Fig. 2(b) in a graphical form with the values of 0 and 1 indicated with yellow and blue . A contains + 1 rows and columns so a sequence of + 1 sampling patterns could be obtained from a single map. However, only the first map is used to create + 1 binary patterns, and each of the following maps is translated into − 1 patterns (by omitting the first two rows of A). This assures that the differential measurement includes nonredundant information about the mean intensity of every region of every map and that the sampling patterns are linearly independent. Otherwise, for instance the information about the mean value of the entire image could be found independently using the sampling patterns obtained from every map. Overall, the number of binary sampling patterns equals = ( + 1) + ( − 1) · ( − 1), the number of measurements effectively used for image reconstruction equals − 1 (one measurement is lost due to the differential treatment of data), and the number of sectors in all the maps is equal to · (which exceeds ). Matrix A is created in the same way as in [29] but is used differently. A is obtained by brute force numerical search under the condition that rank of the matrix obtained by subtracting subsequent rows of A is equal to , i.e.
(D · A) = , and that each of the rows of A contains ± 1 ones and zeros. As a result, the proposed differential sampling probes independently every region of the map, it is binary, and it consists of patterns with approximately half of the pixels in the on and off states. These properties of sampling patterns are advantageous in terms of signal entropy and sensitivity to truncation noise [29].
Let us introduce a denotation for the set of indexes pointing to image pixels that belong to -th sector of -th map where ∈ {1, .. } enumerates the pixels. A useful property of the proposed sampling is that in a noise-free scenario an imagex 0 reconstructed from the compressive differential measurement using formula (2) will have the correct mean values calculated within every region of every map, i.e.
This is true becausex 0 satisfies the noise free measurement equation, and the differential sampling with proposed patterns based on maps is equivalent to non-differential sampling with patterns consisting of all individual sectors extracted from the maps. For nonnegative images x, finding regions with mean value ofx (approximately) equal to zero enables us to mark them as empty and to eliminate them from image reconstruction.

Reconstruction algorithm
The proposed algorithm consists of two stages. First, an initial approximate reconstructionx is calculated using equation (2) with the reconstruction matrix P given by Eqs. (3) and (4). For dense images this will be also the final result. In the second stage, empty image regions are identified from the mean values calculated over every sector of every map for whichx 0 [v ( , ) ] ≤ and pixels from these sectors are set to zero. If this change affects some sector of another map, the remaining pixels from that sector are scaled to restore the proper mean value. The algorithm is detailed below.
The numerical cost of this algorithm ( · ) is proportional to the number of pixels , so the proposed method may be used with high resolution images. At the same time, the compression ratio / has to be low, otherwise the matrix P with · elements becomes too large to be calculated or even stored in computer memory. Later in this paper we will use the denotation MD-FDRI for the proposed modified map-based D-FDRI method. ⊲ Image reconstruction algorithm 2:x 0 ← P · y 3:x ← (x 0 ) ⊲ Initial D-FDRI reconstruction 4: for = 1 to do ⊲ Loop over maps 5: for = 1 to do ⊲ Parallel loop over sectors of map m i 6: ⊲ If yes, set it to zero 8: ⊲ Otherwise, correct the sector's mean 10: end if 11: end for 12: end for 13: returnx ⊲ Return the reconstructed image 14: end function

Discussion
MD-FDRI gives a high quality reconstruction when the image is spatially sparse. For gray-scale dense images it works reasonably well but other SPI methods may give better results. The optimal choice of the parameters an and of the geometric shapes of map regions is object dependent and is beyond the scope of this paper. It should somehow reflect the complexity of the geometric shapes of the sparse regions of images. The most desired situation is when as many region maps as possible entirely overlap with sparse parts of the measured image. In fact, a comparison of results obtained with different values of and in the Sect. S5 of the Supplementary Materials indicates that the reconstruction quality does not vary in a simple or regular way with and .
= 31 was the largest value for which we were able to calculate the lookup table A using brute force optimization and this value works well for various kinds of sparse images.
It may be interesting to make it clear why the lookup table A is at all needed to encode the sampling functions instead of directly using separate image regions as distinct sampling functions. Numerically, the two approaches would give similar results. However, with the MD-FDRI lookup-table encoding, the optical signal measured with a photodiode is higher by a factor of approximately /2=31/2=15 which is substantial. This is because approximately half of the DMD pixels take the "on" or "off" positions for each sampling pattern, while the regions of image maps contain only 1/31 of pixels. A similar result could be obtained with a Hadamard matrix used as the lookup table. However unlike the proposed lookup table, a Hadamard matrix contains a single row consisting of ones, and other rows with half of the elements containing ones. This leads to approximately twice larger signal for one sampling pattern than for the others, which unnecessarily raises the demand for the bit-resolution of the A/D converter used in the measurement. Our lookup table only consists of rows with approximately the same numbers of ones and of zeros. On top of this, using the lookup table A assures that information about the map regions is found from a differential measurement, and if the detector adds systematically some constant bias to the detection signal, this bias does not affect the reconstruction result, while at the same time, the mean value of the image is still measured correctly.
In this paper we pay special attention to the kind of sparsity obtained by limiting the field of view of the observed image. This is practically interesting because by using a diaphragm in an optical set-up one could improve the image quality within the limited FOV without any other modifications to the set-up, sampling functions or to the reconstruction algorithm which would still work at the original resolution ( = 1024 × 768). A full resolution dense image can not be measured accurately at the compression ratio of 0.4% but one could scan the image with a diaphragm in a sequence of measurements to obtain a high quality measurement.
More detailed information on the generation of image maps, on constructing the measurement matrix from the image maps, on the numerical form of the lookup tables with different sizes and on removing linear dependence from the sampling patterns can be found in the Supplementary Materials.

Numerical results
Later in this paper we will use the denotation MD-FDRI for the modified map-based D- In Fig. 3 we present exemplary reconstructions of several test images obtained with the proposed MD-FDRI method in a simulated SPI measurement. The complete image set used in the simulation contains 50 non-sparse test images. We further consider both sparse and non-sparse images. The sparse images are obtained from the same image set by limiting the visible field of view, i.e. by replacing values of all pixels in the image, apart from a selected image area, with zeroes. To compare images with different field of view, we introduce a parameter (0 ≤ ≤ 1), which denotes the proportion of the image area which contains relevant visual information. The placement of the non-empty regions in the image may be arbitrary. We note that to some extent, the significance of parameter is similar to a measure of image sparsity, however they are not strictly equivalent, as the non-empty parts of the image may still be compressible. The reconstructions of both sparse and non-sparse images obtained with the MD-FDRI method show good quality, considering the extreme compression ratio used in the SPI sampling. Moreover, by reducing the field of view to a small region of the image, the quality of reconstruction of this region is strongly enhanced, as shown in Fig. 3(d-f).
Another example of the performance of MD-FDRI is shown in Visualization 1 and selected frames from the visualization are presented in Fig. 4. A synthetic high resolution image filled with a large number of objects is measured through a simulated diaphragm with a varied aperture. A similar diaphragm could find application in SPI microscopy, where the aperture size and location could be moved to the locations of interest which will be then measured at the full resolution of the DMD. The trade off between resolution and field of view is something common to various imaging techniques involving an information channel with a limited bandwidth [35,36], and a similar trade-off is known in classical microscopy as well. We see that although the same non-adaptive sampling patterns are used for all frames, the resolution of reconstructed images vary, and approaches the resolution of the DMD for small aperture sizes.
The enhancement of the reconstruction quality for sparse images is obtained in the second stage of the MD-FDRI reconstruction. MD-FDRI consists of two stages. The first stage produces an initial reconstruction with a single matrix-vector product using Eq. (2), similarly to our previous work [29]. In the second stage, the algorithm iteratively clears empty sectors of the image and corrects the mean value in each non-empty region of every image map. Therefore, the second stage of MD-FDRI yields the better reconstruction quality, the more empty sectors the image actually contains. Fig. 5    These two images as well as every frame in Visualisation 1 are reconstructed each from a separate differential compressive measurement with exactly the same = 3002 binary sampling patterns. For the modulation frequency of 22.5kHz, the measurement time is 0.14s. The reconstruction time of a single frame is on the order of 0.3s. Sampling is nonadaptive and the algorithm does not need to know the location of the diaphragm. This example shows a realistic way to obtain SPI at the full resolution of the DMD with image measurement and reconstruction times on the order of a fraction of a second. second stage of the reconstruction algorithm still allows us to remove the reconstruction artifacts from the empty sectors of the image. However, some improvement to the reconstruction quality of the non-empty regions is possible only if the SNR of the measurement is relatively high. The robustness of the MD-FDRI method to additive measurement noise is further illustrated in Fig. 6.
Finally, in Fig. 7 and Fig. 8 we compare the proposed map-based SPI sampling and MD-FDRI reconstruction routine to two other commonly used sampling and reconstruction scenarios. The first one comprises Walsh-Hadamard (WH) sampling patterns and the inverse transform used for image reconstruction. In the second scenario, binarized DCT sampling patterns are used instead and the reconstruction is obtained using FDRI method, similarly to our previous work [5].    In each case, the sampling matrix consists of patterns at full resolution of the DMD and the compression ratio is 0.4%. The reconstructed images are compared in terms of two metrics: the peak signal-to-noise ratio (PSNR) and the structural similarity index (SSIM). While DCT+FDRI seems to offer the best reconstruction quality for SPI imaging of non-sparse images, in the case of images with high sparsity, MD-FDRI produces much richer reconstructions with significantly more detail than the other considered methods. For = 0.01, average PSNR obtained with MD-FDRI is by over 3 higher than for DCT+FDRI and by over 5 higher than for WH. In general, the reconstructions obtained with MD-FDRI are usually of higher quality that the ones obtained with WH for most values of . We note, that the two metrics: PSNR and SSIM are not always consistent in evaluation of the reconstructions, especially of sparse images. For instance, for some values of the reconstructions obtained with DCT+FDRI have the highest PSNR of all considered methods and the lowest SSIM at the same time. Also MD-FDRI outperforms other methods for much broader range of in terms of SSIM than in terms of PSNR. This is because, these two metrics differently weight the reconstruction errors occurring in the empty regions of the image. The empty sectors are on average reconstructed with much smaller absolute value of the errors than the non-empty ones, but the relative error is still large. According to the PSNR criterion, the reconstruction errors are averaged over the whole image, therefore mostly the accuracy of reconstructing the non-empty sectors of the image contributes to this metric. On the other hand, for SSIM criterion the errors are averaged locally within a Gaussian window around each pixel of the image and the final metric is the mean value of the local ones. Therefore, the accuracy of reconstructing both empty and non-empty areas of the image are equally important, when this criterion is used. The advantage of MD-FDRI over other reconstruction methods in terms of SSIM visible in Fig. 8 reflects the unique ability of this method to locate the empty sectors in the recovered images and clear them from reconstruction artifacts.

Optical results
We have validated the performance of the proposed MD-FDRI algorithm by implementing it in a classical SPI set-up with a DMD used at its full resolution. We were able to obtain high resolution SPI imaging of sparse objects. To this end we have prepared a translucent test pattern consisting of a set of holes in nontransparent metal foil with diameter in the range from 0.15mm to 0.41mm. The object is observed in transmission. The high contrast and brightness of the spots assure a high SNR of the measurement. Our experimental setup is depicted in Fig. 9. The object is illuminated by unpolarized light from a collimated LED source. Light transmitted through the object is collected by achromatic doublet and its image is formed on the DMD modulator (Vialux V-7001 XGA with DLP7000 chip). Compared to imaging at a lower resolution with the DMD pixels grouped into larger rectangular superpixels, here imaging of the object onto the DMD requires pixel-level precision. The DMD modulates the incident beam into two reflected beams. We will use only a single channel here but in fact the D-FDRI allows one to process the complementary channels using just a single reconstruction matrix P [29]. Then, the channels could be used for parallel detection through different spectral or polarization filters [29]. The modulator is operated at 22.7 kHz and at a spatial resolution of 1024 × 768. We use a VIS-NIR amplified photodiode (Thorlabs PDA100A2) as the light detector. The signals are digitized with a digital oscilloscope (PicoScope 5000) at a sampling rate of 1/(256 ns) and streamed through a USB bus for real-time processing. Figure 10 illustrates the performance of MD-FDRI method applied to the high-resolution optical SPI imaging. Figure 10(a) depicts the ground truth image of the sparse test object captured with a 1278x1020 camera. The difference in the brightness of the subsequent spots comes from the  Fig. 10 (b) we present the same image measured using binarized DCT sampling patterns at the resolution 256 × 256 and compression ratio of 3%, which is one of the commonly used SPI state-of-the-art sampling protocols. While in both cases the time needed for image acquisition is comparable, the reconstruction obtained with MD-FDRI provides significantly higher quality and accuracy in recovering image details.

Conclusions
In this paper we have demonstrated that using a classical DMD-based single-pixel imaging (SPI) optical set-up with non-adaptive binary sampling one may be capable of capturing and reconstructing sparse images at a high resolution within a fraction of a second. This is important as after over a decade of research on SPI, the reported set-ups hardly ever exceeded the resolution of 256 × 256 falling far behind current resolution standards.
To this point we have proposed a novel differential binary high-resolution (n=1024 × 768) highly compressive sampling scheme combined with an image reconstruction algorithm with a numerical cost of ( · ). Here is the number of the DMD pixels and is equal to the resolution of the reconstructed images, and ≈ 0.4% · is the number of sampling functions, the acquisition of a single image takes 0.14s, and a typical image reconstruction time is ≈ 0.3s on a Intel i9-10900X 3.70GHz CPU. The same algorithm is capable of reconstructing non-sparse images with reduced quality due to the strong compression. This work paves the way for practical high resolution SPI applications without more sophisticated modulation or detection elements than just a DMD and a photodiode. Additionally, we point a way to vary the SPI imaging quality using a simple diaphragm. Besides the usefulness of this property for instance for SPI microscopy, we would like to point that the trade off between resolution and field of view is something common to various imaging techniques involving an information channel with a limited bandwidth. A similar trade-off is known in classical microscopy as well. By limiting the field of view and decreasing the information content within the measured image it is possible to improve SPI imaging quality up to the full resolution of the DMD without any further adjustments to the measurement set-up or SPI algorithm. Data Availability Statement. Source data and source code will be provided by the authors at a reasonable request.
Disclosures. The authors declare no conflicts of interest.

S1. Introduction
This document contains supplementary materials to our Optics Express paper entitled "Single pixel imaging at high pixel resolutions" which we will further call the main paper. Supplementary materials provide more detailed information on the generation of image maps (see Section S2), on constructing the measurement matrix from the image maps and removing linear dependence from the sampling patterns (see Section S3), and on the resolution difference between the Fourier domain regularized inversion (FDRI) and image map-based Fourier domain regularized inversion (MD-FDRI which is introduced in the main paper) methods (see Section S4). In the last section S5 an atlas of the lookup tables used for encoding maps with different number of image regions is given and sample SPI results obtained with varying parameter values are compared.

S2. Calculation of image maps and of the measurement matrix
In the main paper, the binary patterns included in the rows of the measurement matrix M play a double role. First, the sampling provides the usual information about the spatial contents of the measured image. On top of this it should give some guarantees about the possibility of identifying empty regions of the image. With this aim, sampling patterns are generated with the help of auxiliary image maps m 1 , m 2 , ...m whose role is discussed in Sect. 2.3 of the main paper. Every map defines a partitioning of all image pixels into groups . The first map is used to create + 1 sampling patterns, and every other map is used to create − 1 patterns. The particular form of the maps used by us is rather arbitrary. It could be further optimized but image maps very different from those proposed in the main paper could also work fine with the MD-FDRI algorithm. This is why in the main paper we have not paid too much focus on the way we generate the maps. In short, the maps are obtained by generating spatially correlated Gaussian complex-valued noise and by assigning the label ∈ {1, } based on the uniform quantization of the phase level of a noise realization. The precise procedure is defined by the Python program shown in Listing 1. Sample image maps generated by this program are shown in Fig. S1 1 import numpy as np 2 def create_image_map(dim=(768, 1024), m=31):   The measurement matrix M is composed of rows that represent binary sampling patterns. These sampling patterns are obtained by processing the list of image maps. Translation of image maps into binary patterns is performed with the help of subsequent rows of an auxiliary binary matrix A which play the role of lookup tables. This procedure is described in Sect. 2.3 and in Fig. 1. of the main paper. The following Python function may be used to build a list of image maps and to calculate the measurement matrix: Create the binary measurement matrix M with rows representing sampling patterns

S3. The lookup table -a differential binary code for encoding binary patterns
Here we overview the properties of the lookup table A that we use to translate the image maps into sampling patterns. A is a binary + 1 × matrix. The main property that A should fulfill is that D · A is a full-rank matrix, i.e.
(D · A) = . Given an image map, subsequent rows of A are used to construct a sequence of binary sampling patterns. The columns of A are linked to regions of pixels of that image map. Since D · A is a full-rank matrix, the differential measurement with the so created sampling patterns is equivalent to a measurement with patterns defined by separate pixel regions of the map. Yet, the proposed approach is better in terms of light efficiency and it is differential so a constant detector bias signal (e.g. from the dark current) is automatically disregarded in the measurement. By keeping similar surfaces of pixel regions within a map and by choosing matrices A with a similar number of ones and zeros in every row we assure that the SPI detection signal does not experience large variations. This comes in contrast to SPI schemes in which the mean value of image is measured with all DMD pixels in one position while the other patterns consist of approximately half of pixels in one position, a situation which leads to unnecessarily high requirements for the bit-depth of the DAQ.
We note that any two matrices A with permuted columns may be considered equivalent because the columns are arbitrarily assigned to pixel regions of an image map. We also note that when two or more image maps are used in a sequence to produce the corresponding sampling patterns, one could easily end up with measuring redundant information. One origin for this redundancy is that all the regions of any two maps cover exactly the same area (namely all image pixels) and for instance the image intensity could be measured with sampling patterns built from any map. The second is less obvious and results from the differential operator acting between the last and first sampling patterns from two subsequent maps. In effect to avoid redundancy only the first image map should be used to build + 1 sampling patterns and for each following map only − 1 patterns need to be created (per map). For this purpose we neglect the first two rows of A for all image maps but the first one. For this to work it is not sufficient to know that (D · A) = but also it is necessary to make sure that a block matrix composed of the last row of A as one block, and A with first two rows taken away as the second block, after applying the difference operator and adding a row consisting of ones and minus ones (representing an equation for the image intensities measured with two image maps) is a full rank matrix. This condition may be easier to formulate with a piece of Python code than with an equation:   ((1,m)),-np.ones ((1,m))) ) )) 21 tst=np.linalg.matrix_rank(b)==m-1 # is it ok to disregard first two rows of A for subsequent image maps?

S4. Resolution test
Nonadaptive compressive imaging at a very low compression ratio (like the one considered here / ≈ 0.4%) is challenging when the images are neither sparse in the spatial nor in the Fourier domains. The most common approach in optical SPI is to use sampling with a low spatial frequency subset of some basis functions (such as DCT, Fourier, Walsh-Hadamard basis etc.). In fact, even for sparse high resolution images their spatial spectra have usually the highest amplitudes at low spatial frequencies. On the other hand, retaining a high resolution is impossible without high spatial frequency information.
In Fig. 5 of the main paper we have shown MD-FDRI imaging results of a Siemens resolution test (with and without noise, and at the full or limited field of view). Here we compare these results with DCT-FDRI imaging (FDRI with binarized DCT sampling functions not based on image maps, and with a single step FDRI image reconstruction). The DCT-FDRI results are presented in Fig. S2(b,c) which may be compared to Fig. 5(c,e) of the main paper. The DCT-FDRI performs better for an image with a full field of view both in a noise-free and noisy situations. This is confirmed by the values of PSNR and SSIM criteria as well as by visual examination. However, the center of the resolution test appears blurred, also when the field of view is limited. Therefore DCT-FDRI cannot be called a high-resolution method. On the other hand, for all the  Fig. 5(a,c,e) of the main paper which shows similar results obtained with MD-FDRI obtained with the same number of sampling patterns ( = 3002) and at the same resolution ( = 1024 × 768). a) ground truth; b) noise-free measurement; c) noisy measurement situations with a limited field of view, the MD-FDRI performs better in terms of PSNR and SSIM. Also a visual examination allows to appreciate the resolution obtained within the central part of the Siemens star (See the right columns of Fig. 5(c,e) of the main paper).

S5. Atlas of binary lookup tables
The optimal choice of the parameters and is highly object dependent. For gray-scale dense images, other SPI methods may give better results than MD-FDRI. For sparse images or images with a limited FOV, optimal and and kind of image maps depend on the image complexity and kind of features included in an image. In Figs. S3, S4 we show sample image reconstructions with and without noise for different values of an taken in such a way that the total number of sampling patterns is approximately constant. The results from this comparison indicate that the reconstruction quality does not vary in a simple or regular way with and . In the end of this section we provide a list of binary lookup tables of different sizes (for odd values of up to = 31). We note that these matrices are not unique. From our experience, usually the larger is , the better for the SPI imaging quality, but the results vary depending on the composition of sparse images. The lookup table used by us in this work is the largest we have calculated (with brute force binary search that becomes practically impossible to continue for larger ).                 Disclosures. The authors declare no conflicts of interest.