Single Photon Compressive Imaging Based on Digital Grayscale Modulation Method

In single-pixel imaging or computational ghost imaging, the measurement matrix has a great impact on the performance of the imaging system, because it involves modulation of the optical signal and image reconstruction. The measurement matrix reported in the existing literatures is first binarized and then loaded onto the digital micro-mirror device (DMD) for optical modulation, that is, each pixel can only be modulated into on-off states. In this paper, we propose a digital grayscale modulation method for more efficient compressive sampling. On the basis of this, we demonstrate a single photon compressive imaging system. A control and counting circuit, based on field-programmable gate array (FPGA), is developed to control DMD to conduct digital grayscale modulation and count single-photon pulse output from the photomultiplier tube (PMT) simultaneously. The experimental results show that the imaging reconstruction quality can be improved by increasing the sparsity ratio properly and compressive sampling ratio (SR) of these gray-scale matrices. However, when the compressive SR and sparsity ratio are increased appropriately to a certain value, the reconstruction quality is usually saturated, and the imaging reconstruction quality of the digital grayscale modulation is better than that of binary modulation.

The measurement matrix selected in the single photon compressive imaging system can obtain the relevant information of the image and determine the quality of image reconstruction [18,19], which influence the performance of the imaging system. To study the better matrix performance of the measurement matrix in the single photon compressive imaging system, some researchers have designed some new measurement matrices in recent years. In 2015, Tiwari et al. [20] designed the sparse sensing matrix to obtain high resolution medical images with least incurred computational cost. In 2017, based on Bernoulli and Gaussian random matrices, Nouasria et al. [21] constructed a new measurement matrix for compressed sensing by enhancing the orthogonality of the measurement matrix column. In 2018, Ding et al. [22] proposed a global perception matrix design method, which improved the performance of image compression in the multi-dictionary CS system. The measurement matrix constructed by the above researchers based on compressed sensing theory improves the quality of image reconstruction. But in the single photon compressive imaging system, the measurement matrix loaded on digital micromirror device (DMD) is required to be a binary matrix, and the measurement matrix composed of 0 and 1 on the DMD contains less information, and if it is applied to gray image reconstruction, most of the information will be lost and it is difficult to ensure the accurate recovery of image information. The signal-to-noise ratio (SNR) is improved due to the large number of modulations of the gray matrix. Therefore, the modulation effect of the gray matrix is better than that of the binary matrix.
This paper proposes a digital grayscale modulation method, and constructs three different random measurement matrices with gray level through matrix transformation of different schemes. These matrices contain quaternary Gaussian random matrix, octal Gaussian random matrix, and hexadecimal Gaussian random matrix. Meanwhile, we construct a single photon compressive imaging system. In this system, the gray image is displayed by spatial light modulation of DMD and high precision timing control of field-programmable gate array (FPGA). Finally, we quantify the performance of the single photon compressive imaging system under different measurement matrices.

2.
Principle and realization of experimental system The single photon compressive imaging system is constructed as shown in Fig. 1 Fig. 1 Schematic diagram of a single photon compressive imaging system, including a parallel single photon source (it consists of light emitting diode (LED) lamp, collimator, attenuator, and diaphragm), imaging target, imaging lens, DMD, focusing lens, and PMT.
The system adopts the idea of compressed sensing for compressed sampling and image reconstruction. The modulation process of each frame matrix for the optical image on the DMD can be regarded as the inner product between the measurement matrix on the DMD and the original image. The inner product is the number of photons detected in the direction of reflection of +12° during each measurement time. The measurement process is carried out under the synchronous control module based on the FPGA. In the experiment, the measurement matrix is loaded sequentially by the FPGA to the DMD for modulating the image, and synchronously records the number of photons output by a photon counting photomultiplier tube (Hamamatsu Photonics H10682-110 PMT) with an effective photosensitive area of Φ8 mm. And the number of photons is proportional to the intensity of light. Finally, the image is reconstructed by the compressed sensing algorithm according to the measurement matrix loaded to the DMD and the photon sequence y is received by the synchronization control module.
Let the image x be a one-dimensional signal of length N and the sparsity K (i. e., containing K non-zero values). The mathematical model of the above compression measurement process can be expressed as follows: In the formula, e is the measurement noise, and Φ is the measurement matrix of the M N  dimension, and M is much smaller than N. We can reconstruct the image x from the obtained M-dimensional measured values. The measured value y is an M-dimensional vector, which reduces the measurement object from N-dimensional to M-dimensional.
where x is the original image. In the theory of compressed sensing, there are many algorithms for solving (3). The greedy algorithm directly based on the L 0 norm solution [such as matching pursuit (MP) and orthogonal matching pursuit (OMP)] converts L 0 into L 1 relaxation algorithm of L 1 problem [such as base tracking (basis pursuit, BP)] and total variation algorithm [such as total variation augmented lagrangian alternating direction algorithm, TVAL3)] [23,24]. In this paper, combining the performance of various algorithms, the TVAL3 algorithm is used as the reconstruction algorithm of the target image, and it has lower complexity and faster reconstruction.

Implementation of DMD digital gray matrix modulation
It has been proved theoretically and practically that matrices satisfying the restricted isometry property (RIP) can be used as the compressed sensing measurement matrix. Bernoulli and Gaussian matrices are proposed to satisfy restricted isometry property with high probability and have the advantages of a simple structure and a high reconstruction efficiency [25][26][27]. At present, the Gaussian random measurement matrix is generally transformed into 0-1 binary matrix and applied to the DMD, which can modulate the simple binary image. But when applied to gray-level images, the binary matrix carries less information and increases the computational complexity in the observation process and has a great impact on the recovery accuracy. If the Gaussian random measurement matrix is converted into a corresponding gray matrix, the same imaging effect can be achieved correspondingly, and the reconstruction accuracy of the gray image is also improved.

Construction of the digital gray-scale matrix
Firstly, this scheme generates a matrix Φ of M N  , and makes every element between 0 and 1 in matrix Φ obey the Gaussian distribution with the mean value of 0 and the variance of 1/M independently.
The measurement matrix has strong randomness. It can be proved that when the number of measurements of the Gaussian random matrix is lg( / ) M cK N K  , the RIP condition is satisfied with a maximum probability [28]. Gaussian random matrix is widely used in CS, mainly because it is not related to most orthogonal basis or orthogonal dictionary. At the same time, because of its simple structure and it is easy to be constructed, the image can be accurately reconstructed from a small amount of measured data.
Then, each element in the Gaussian matrix constructed is converted into a measurement matrix whose matrix element values are independently distributed between 0 and 3, 0 and 7, and 0, and 15 according to corresponding coefficients of a certain proportion. Finally, the values in the matrix are rounded and converted into the corresponding matrix elements of Gaussian random matrices in quaternary, octal, and hexadecimal, respectively.

Implementation of DMD grayscale modulation
In this paper, the method of DMD grayscale modulation is applied to the single photon compressive imaging system for the modulation of gray matrix, and the binary pulse width modulation (PWM) method is used to realize the full digital control of the image gray level by the DMD. PWM uses binary and division of time pulse width modulation technology to process binary weighted image data. The length of the pulse width indicates the reflection time of the position of the micromirror to determine the magnitude of the gray value at that location, so as to control the time of a single micromirror switch to form the gray modulation of each pixel of the image [29,30]. The method of grayscale display based on the DMD adopted in this paper is the same as accessory light modulator package (ALP) in principle [31]. By this method, the input image signal of the DMD can be converted into gray image with rich levels. Binary PWM modulation can be expressed as follows: where N is the gray level of the pixel, T PWM is the weight of each bit pixel, and T LSB is the display time of the lowest bit (LSB). As shown in Fig. 2, it is an example of a pulse sequence controlled by PWM for 4-bit binary image data with 10 gray values. It has 2 4 gray levels, i.e., 0-15 gray values, and the binary number of 10 gray values is 1010. The process of grayscale modulation is equivalent to dividing the 4-bit image into four binary images, in which gray represents the time when the micromirror is on and white represents the time when the micromirror is off. From the lowest significant bit (LSB) to the highest significant bit (HSB) modulation, a frame is divided into four separate bit time slices. The length of time slices is proportional to the binary weight of the bit. When the HSB is allocated to 8/15 frames, HSB-1 is allocated to 4/15 frames, HSB-2 is allocated to 2/15 frames, and LSB is allocated to 1/15 frames. The binary "1" of the PWM control signal will drive the DMD to the "on" position, while "0" will make the DMD move to the "off" position. In HSB time slice and HSB-2 time slice, the DMD is in "on" state, while in HSB-1 time slice and LSB time slice, the DMD is in "off" state. After the integration of this display sequence, the modulation effect equivalent to gray value 10 is produced, so four binary images are superimposed on the display screen to complete the display of this 4-bit image. 8-bit and 16-bit images can also be implemented, which mainly depends on the reset time of the lens and the field time of each frame of image. If the gray level of an image is N, the display of each grayscale image in the DMD can be divided into n＝log 2 N number of time slices of data bits. The image is converted into n binary black and white images, and f (x, y) represents a grayscale image, then: where ( , ) i b x y is the binary image displayed by the DMD within the ith number of time slices of data bits. Therefore, adjusting the DMD display frame rate is equal to adjusting the stabilization time of the micromirrors. In summary, the process of adjusting the display frequency of the DMD is the process of modulating the grayscale of the DMD, that is, controlling the PWM.

Control module and timing design based on FPGA
The scheme adopts random measurement matrices with the gray level. We design a high-precision control and counting module based on the FPGA, which can control the DMD working in the external synchronization mode to load random measurement matrices sequentially and record the number of photons of each frame synchronously, that is, the measured value synchronously.
The control and counting module based on the FPGA is implemented on the Cyclone IV chip of Altera DE2-115 development board. The main working sequence of synchronous control pulse is shown in Fig. 3. The sampling frequency pulse is generated by the FPGA, and the sampling frequency pulse signal period is consistent with the time of the last pulse control signal. The number of the gated square wave signal is consistent with the random measurement matrix, and the corresponding control pulse rising edge in each square wave is aligned with the sampling frequency pulse signal generated by dividing frequency. The number is equal to the preset sampling and the interval between each pulse is proportional distribution. When the number of gated square-wave control pulses is 3 in Fig. 3, the control signal pulse time interval is T1:T2:T3=4:2:1. When a rising edge of this synchronization pulse signal is received, a random binary matrix is loaded into the DMD, and the measured data can be obtained by the PMT. When the next synchronization pulse signal is coming, the measured data and random binary matrix are sent to personal computer (PC) for image reconstruction by CS algorithm. The numerical octal random matrix based on the transformation of the Gaussian matrix is shown in Fig. 4, and the values of some elements without numbers are 0. In this scheme, the measurement matrix is loaded on the DMD with a size of 1 024 768  . The random measurement matrix is composed of M N  matrix elements, and the elements of each matrix are no longer "0" or "1", but also have other numbers. Therefore, the DMD micromirror arrays corresponding to a matrix element are grouped together as one pixel. In this scheme, if the value of the representative measurement matrix coefficient is higher, the display time on the corresponding DMD pixel position is longer, and the coefficient value is proportional to the display time. Taking the octal Gaussian random matrix of Fig. 4 as an example, the process of independent representation of the corresponding numerical matrix is shown in Fig. 5. First, the octal numerical correspondence in the matrix is transformed into a binary representation. Then, the corresponding high-to-low-bit elements are taken out sequentially, and three measurement matrices are separately obtained. Under the synchronous control pulse signal of the FPGA, three measurement matrices are displayed from high to low sequentially, and the display time is determined by the synchronous control pulse signal.

Detection limit and experimental performance verification
We assume that the image is K-sparse, and the average photon count rate reflected by the micromirror is R ph cps. If the measurement time is t for the imaging target, the measured photon count is The pulse count output by the detector is the sum of the signal photon number W ph , the background photon count W b , and the dark count W d . As the illumination on the object is at the single photon level, the measurement noise in our imaging system is dominated by Poisson shot noise, and the total noise component can be analyzed as follows: We assume that the background photon count rate is R m and the dark count rate is R n . The SNR of single photon imaging system will be In the photon counting mode, if we define the detection limit as the light level where the SNR equals to 1 and the measurement time is one second, the photon counting rate at the detection limit can be approximately expressed as Through the actual system experiment, the photon counting rate range is 10 2 counts/s-10 3 counts/s. To compare the reconstruction performance of different matrices more accurately, the photon counting rate of the detector is used to represent the intensity of light. The photon counting rate at the light intensity level is approximately 300 counts/s. The image pixels in the experimental system are 6464, and only 13 756 detected photons are contained in the measured data, corresponding to 3.36 photons per pixel.
In the process of compressed sensing, the performance of the measurement matrix is closely related to the sparsity ratio of the matrix and the compressive sampling ratio (SR). The performance of the measurement matrix can be evaluated by the reconstruction quality of image. For a quantitative comparison of the image quality, we introduce three commonly used full-reference evaluation indexes, namely, the mean square error (MSE), peak signal to noise ratio (PSNR), and mean structure similarity index (MSSIM). In order to evaluate the resolution of the reconstructed image, the customized mask pattern is used as the imaging target of the experimental system. The image reconstruction algorithm is implemented by the TVAL3 algorithm.

Influence of sparsity ratio on imaging quality
Because this scheme adopts random measurement matrix with the gray level, the influence of sparsity ratio of the measurement matrix on matrix performance is studied through simulation experiments. Taking the construction of quaternary Gaussian random matrix as an example, the sparsity ratio of the matrix is defined as the proportion of non-zero elements of the matrix. In order to obtain multiple sets of data to prove the universality of the results, reconstruction experiments at different compressive SRs are set up. The compressive SRs of the three groups are set to 0.25, 0.50, and 0.75, respectively, and the matrix sparsity ratio of each group increases continuously.
The experimental results are shown in Fig. 6. In order to obtain the simulation results of the measurement matrix under different sparsity ratios, Figs. 6(a), 6(b), and 6(c) correspond to the indexes MSE, PSNR, and MSSIM, respectively, and each data point has an average of 500 independent simulation experiments. We can see that the quality of gray image reconstruction will increase with an increase in compressive SR, and the quality of the reconstructed gray image is closely related to the matrix sparsity ratio. When the sparsity ratio of matrix is small, the quality of the reconstructed image will increase with an increase in the matrix sparsity ratio. But for a larger sparsity ratio (greater than 0.15), the reconstruction quality stays almost unchanged. Therefore, there is a tradeoff between the matrix sparsity ratio and the compressive SR. When the matrix sparsity ratio is greater than 0.15, the SR will be increased as much as possible to improve the quality of gray image reconstruction. However, in the actual single photon compressive imaging experiment, while reducing the sparsity ratio of the matrix, sufficient luminous flux should be guaranteed to obtain a higher SNR.
In order to validate the above conclusions, we use custom mask of patterns as imaging targets. The compressive SRs of three groups of experiments shown in Figs. 7(a), 7(b), and 7(c) are set to 0.25, 0.50, and 0.75, and the sparsity ratio of each group is increased in turn. Three groups of experiments are set to 0.05, 0.10, and 0.15, respectively, and Fig. 7 shows the experimental results of the quaternary Gaussian random matrix under these conditions. The first group of experiments shown in Figs. 7(a1)  7(a3) represents the reconstruction results of increasing the sparsity ratio of the matrix sequentially when the compressive SR is 0.25. The second group of experiments shown in Figs. 7(b1)  7(b3) represents the reconstruction result of increasing the sparsity ratio of the matrix sequentially when the compressive SR is 0.50. The third group of experiments shown in Figs. 7(c1)  7(c3) represents the reconstruction result of increasing the sparsity ratio of the matrix sequentially when the compressive SR is 0.75. As can be seen from Fig. 7, the quality of the reconstructed image is improved with an increase in matrix sparsity ratio.

Influence of compressive sampling ratio on imaging quality
In order to analyze the influence of the number of measurements on the quality of image reconstruction, we set different compressive SRs to verify the corresponding experiments results.
Gaussian binary random matrix, quaternary Gaussian random matrix, octal Gaussian random matrix, and hexadecimal Gaussian random matrix are used as measurement matrices for compressed sampling and image reconstruction. At the same time, the orthogonal matrix of partial Fourier is used for the experiment by contrast.
The imaging results corresponding to the above five random matrices are shown in Fig. 8. The pattern of the flower is used as the imaging target of the experimental system corresponding to the Gaussian binary random matrix. Taking the Chinese character pattern of "light" as the imaging target, it corresponds to the quaternary Gaussian random matrix, and the United States Air Force (USAF) resolution test board is selected as the imaging target of the experimental system corresponding to the octal Gaussian random matrix. The "aircraft" pattern is taken as the imaging target of the experimental system corresponding to the hexadecimal Gaussian random matrix. The pattern of the "triangle combination" is used as the imaging target of the experimental system corresponding to the hexadecimal Fourier matrix. Among them, Figs. 7(a1)  7(a4), 7(b1)  7(b4), 7(c1)  7(c4), 7(d1)  7(d4), and 7(e1)  7(e4) correspond to the above five experiments of image reconstruction. In a single group of the above five experiments, there are four independent groups with corresponding compressive SRs of 0.10, 0.25, 0.50, and 0.75, respectively. The sparsity ratio of the selected measurement matrix is 0.15, and the imaging resolution is 64 64  . As can be seen from Fig. 8, based on the compressed sensing theory, it can recover the original signal from a small number of measurements, and with an increase in the number of measurements, the quality of image reconstruction becomes better and better. In order to compare the quality of the reconstructed image of the imaging object in the experimental system more intuitively, the Gaussian gray matrix and the Fourier gray matrix are evaluated by the MSE, PSNR, and MSSIM. The reconstructed image indexes change with the compressive SRs as shown in Fig. 9. From the MSE, PSNR, and MSSIM evaluation indices corresponding to the reconstructed gray images, it can be seen that when the compressive SR is less than 0.40, the reconstruction performance of the three gray-scale matrices is relatively poor, while when the compressive SR is greater than 0.40, the performance of the random matrices is improved with an increase in the compressive SR. In addition, when the compressive SR is close to 0.10 and 0.75, the reconstruction performance of the three gray-scale matrices has little difference.
Compared with the Gaussian gray matrix, the Fourier gray matrix has better performance. And with an increase in the gray level of gray matrix, the reconstruction performance of these two kinds of gray matrices is much better than that of binary matrices.

Comparison of imaging performance of different measurement matrices
The performances of different types of measurement matrices are analyzed and verified by simulation experiments. In order to obtain the evaluation indexes of reconstructed images under different measurement matrices, the simulation results are averaged for 1 000 independent experimental results to reduce the impact of randomness on experimental results. The compressive SR is 0.50, and the sparsity ratio is 0.15. Table 1 presents the reconstruction results of the MSE, PSNR, and MSSIM of the reconstructed image under the same experimental conditions. From the experimental results in Table 1, we can see that the hexadecimal Gaussian random matrix has the best performance, and the MSSIM is 0.900 3, followed by the octal Gaussian random matrix.

Conclusions
In this paper, a digital grayscale modulation method is proposed to improve the efficiency of compressive sampling and a single photon compressive imaging system is established. With this system, we realize the reconstruction of the gray images by digital grayscale modulation of the DMD based on gray-scale matrices, which are improved based on the Gaussian random matrix. The control and photon counting module are designed based on the FPGA independently. This paper analyzes and constructs new gray-scale measurement matrices. The quality of gray image reconstruction is compared among three kinds of gray matrices and binary matrix, and the gray matrix has better performance than binary matrix. In addition, the quality of image reconstruction can be improved by reducing the sparsity ratio of the matrix appropriately and increasing the compressive SR. In order to obtain better image reconstruction quality under conditions of limited acquisition and the sparsity ratio of the matrix, the sparsity ratio is preferably 0.15, and the compressive SR is preferably 0.50. The experimental results show that the image reconstruction performance of the gray matrix with a higher gray level is better.