Robust Entangled-Photon Ghost Imaging with Compressive Sensing

This work experimentally demonstrates that the imaging quality of quantum ghost imaging (GI) with entangled photons can be significantly improved by properly handling the errors caused by the imperfection of optical devices. We also consider compressive GI to reduce the number of measurements and thereby the data acquisition time. The image reconstruction is formulated as a sparse total least square problem which is solved with an iterative algorithm. Our experiments show that, compared with existing methods, the new method can achieve a significant performance gain in terms of mean square error and peak signal–noise ratio.


Introduction
Ghost imaging (GI) has raised increasing interest recently due to its wide applications ranging from biological sciences to security protocols [1,2]. In an entangled-photon GI system, the object reconstruction is based on two correlated optical beams, i.e., the object beam and the reference beam. The object beam emits light through an object, which is monitored by a bucket detector. The reference beam does not interact with the object and the light is monitored by a spatial-resolution detector. The first entangled-photon GI experiment was demonstrated by using entangled photons generated by spontaneous parametric down conversion (SPDC) [3][4][5][6][7]. Quantum imaging with entangled photons suffers from low-efficiency due to the low-flux of entangled photons, and it is also time-consuming. To solve this problem, compressive sensing (CS) was introduced into quantum GI [8][9][10], which greatly reduces the number of measurements and thereby the acquisition time [10][11][12]. Some investigations have been conducted to improve the quality of imaging in terms of peak signal-noise ratio (PSNR) and the mean square error (MSE). To improve the quality of reconstruction and reduce the number of samples in the GI system, many CS methods have been employed, which include Orthogonal Matching Pursuit (OMP) [13], Gradient Projection for Sparse Reconstruction (GPSR), etc. Compared to traditional quantum GI, the use of these algorithms could provide higher PSNR and lower MSE.
GI was originally performed using entangled-photon pairs [14], and then was realized with thermal light [15]. In thermal GI, a laser beam is used to illuminate an object, and a light is collected by a single-pixel bucket with no spatial resolution. By combining CS and GI, the spatial resolution of recovered images can beat the diffraction limit of GI [16]. Thermal GI has potential in practical applications [17], such as X-ray tomography [18], astronomy [19] and single-pixel imaging [20][21][22]. Compared with the thermal-light GI [23], quantum GI can obtain higher-visibility and imaging quality [24]. Quantum imaging can break the resolution limit of Rayleigh diffraction [25,26] to achieve super-resolution and strong anti-interference capability. However, the performance of the GI system can suffer from errors caused by the imperfection of optical devices and measurement. In particular, the spatial light modulator (SLM) which modulates the entangled photons in the GI system causes modulation errors. The errors caused by the imperfection of the optics devices were not considered in the literature.
In this work, we experimentally show that, by properly handling the errors caused by the imperfection of optical devices in an entangled-photon GI system, significant performance improvement can be achieved. Compressive GI is considered to reduce the data acquisition time. The image reconstruction is formulated as a sparse total least square (STLS) problem, that is solved using an iterative algorithm. Both simulation and experimental results are provided to demonstrate that a significant performance gain can be achieved by the proposed method in terms of PSNR and MSE, compared with existing compressive sensing-based methods.

Robust Ghost Imaging Based on STLS
We use a quantum compressive GI system where SPDC is used to generate photons, and entangled photons are used as the light source. The quantum GI setup [3][4][5][6][7] is shown in Figure 1. Figure 1. The experimental schematic of STLS quantum GI: HWP, half-wave plate; BBO, β-barium borate crystal; BS, beam splitter; Random patterns placed on a spatial light modulator (SLM); SPCM, single photon counting modules; C.C, coincidence measurement between SPCM1 and SPCM2.
SPDC is caused by random vacuum fluctuations, and the generation of entangled-photon pairs is random. The conversion efficiency of this process is extremely low. The bi-photon state can be represented as: where | is the bi-photon state, ( , ) is the optical-field of the pump, ( ) and ( ) are the creation operators of signal light and reference light respectively. and are the positions of two optical paths at BS, respectively. and are the positions of photons at the cross-section of SPCM1 and SPCM2, and |0,0 is the vacuum state [27]. The coincidence count signal for the detectors SPCM1 and SPCM2 can be expressed by the fourth-order correlation function of optical field intensity as follows: where ( ), ( ), ( ), and ( ) are the positive and negative frequency part of the optical field operator at coordinate (x, y) [19]. The optical fields on the detection plane, which appear in the second line of (2), are given by: SPDC is caused by random vacuum fluctuations, and the generation of entangled-photon pairs is random. The conversion efficiency of this process is extremely low. The bi-photon state can be represented as: where |ψ is the bi-photon state, ϕ(x s , x i ) is the optical-field of the pump,â + s (x s ) andâ + i (x i ) are the creation operators of signal light and reference light respectively. x s and x i are the positions of two optical paths at BS, respectively. x and y are the positions of photons at the cross-section of SPCM1 and SPCM2, and |0, 0 is the vacuum state [27]. The coincidence count signal for the detectors SPCM1 and SPCM2 can be expressed by the fourth-order correlation function of optical field intensity as follows: whereÊ + s (x),Ê + i (y),Ê − s (x), andÊ − i (y) are the positive and negative frequency part of the optical field operator at coordinate (x, y) [19]. The optical fields on the detection plane, which appear in the second line of (2), are given by: The free-space propagation function of the object arm is: where h(x, α) is the free-space propagation function from α to x, and the free-space propagation function of the reference arm is expressed as follows: where the L(β) = exp iπ/(λ f )β 2 is the optical transfer function of the lens. Here, f = 100 mm is the focal length of the lens. Based on the Fresnel approximation, the two-photon amplitude can be expressed as: where the free-space propagation function can be written as h(x, is the transmission function of the object, ∆e represents measurement error, and A i (ρ) is a random pattern loaded onto the SLM. The bi-photon state generated by SPDC can be approximated by [28,29]. When the experiment satisfies the thin lens equation [30], and where (d 1 + d 4 )/d 2 = 1.5 is the theoretical magnification factor of our imaging system, where f is the focal length of the lens. Accordingly, two-photon amplitude is given by: Substituting (1), (3), (4) into (2), we have: Then, the integrated coincidence signal becomes: where m ∈ 1, 2, . . . , M with M being the total number of measurements, and n ∈ 1, 2, ..., N with N being the number of pixels corresponding to the object. In (10), |A m (−α n ) + ∆A| 2 can be rewritten as: where ∆S mn = 2 × A m (−α n ) × ∆A + |∆A| 2 accounts for the errors caused by the imperfection of optical devices, and A mn = |A m (−α n )| 2 is the (m,n)th element of the sensing matrix A M×N . Equation (10) can be rewritten in matrix form as: Sensors 2019, 19, 192 4 of 11 where C M×1 is the measurement vector, and A M×N is the sensing matrix, M is the number of samples. The N-dimensional unknown signal T N×1 is a vector constructed by T n = |T(α n )| 2 . ∆E is a measurement error vector constructed by ∆e. In our method, pseudo-random measurement matrices A M×N (M < N, N = 64 × 64) are selected as the sensing matrices. Given the initial parameters, we can obtain a complete sequence of numbers. The general form of the pseudo-random sequence can be written as Then, we rearrange the sequence into a matrix with size M × N. This ill-posed [31][32][33] problem can be solved by exploiting the sparsity of the signal if A M×N satisfies the restricted isometry property (RIP) [34,35]. Compared with TLS and other CS algorithms, STLS integrates the advantages of compressive sensing and TLS. On one hand, STLS can reconstruct images with only a small number of samples by exploiting the sparsity of the signal, i.e., STLS can deal with ill-posed problems while TLS cannot. On the other hand, STLS is able to handle the errors caused by both the imperfection of optical devices and measurement, thereby achieving robust GI and high imaging quality, which makes it superior to conventional CS algorithms as the conventional compressed sensing algorithm does not consider errors in the dictionary matrix.
We use symlet wavelet in discrete wavelet transform (DWT), to transform the image into a sparse domain, then reconstruct the image with STLS algorithms in the sparse domain. We assume that the wavelet transform matrix is W N×N , then W N×N × T N×1 = θ N×1 , where T N×1 is an unknown signal vector, and θ N×1 is the representation of T N×1 in the sparse domain. Accordingly, Equation (12) can be rewritten as: Then, we can obtain the flowchart of the STLS algorithm as shown in Figure 2, and the two steps are executed iteratively until an ideal solution is obtained.  Due to the errors existing in the model shown in Equation (8), the optimization problem can be formulated as: . .
the cost function includes two parts: the error term ‖[Δ , Δν]‖ and the regularization term ‖ ‖ where γ is a parameter to control the sparsity of the solution. Clearly, when is equal to zero, the problem is reduced to the TLS. The optimization problem is non-convex. We set = 300 and =  Due to the errors existing in the model shown in Equation (8), the optimization problem can be formulated as: the cost function includes two parts: the error term [∆E, ∆ν] 2 F and the regularization term λ θ 1 1 where γ is a parameter to control the sparsity of the solution. Clearly, when γ is equal to zero, the problem is reduced to the TLS. The optimization problem is non-convex. We set γ = 300 and δ = 1 × 10 −10 . We use a coordinate descent method to solve the problem. The method fixes a parameter between ∆ν and θ, while optimizing the other one, until a stop criterion is satisfied. The flowchart of the algorithm is illustrated in Figure 2. Once obtaining θ, we use T N×1 = W −1 N×N × θ N×1 to transform θ to T, and then use T to obtain the reconstructed image.

Numerical Simulation Results
In order to examine the effectiveness of the proposed scheme with STLS, we compare the performance of the scheme with that of the scheme with OMP, GPSR, Method proposed in [19] and TLS algorithms. The method proposed in [19] is direct CS method. We use MSE and PSNR to evaluate the reconstruction quality. For an L × W image, the MSE and PSNR are defined as: where C i,j represents the true image and C i,j denotes the reconstructed image, lg is the base-10 logarithm function, and C max is the maximum pixel value of the image. In our simulations, L = 64 and W = 64. Generally, the larger the PSNR, the better the quality of the reconstructed image. All the simulations are run using MATLAB 2014a in a computer with configuration: Intel(R) Core (TM) i7-7700 CPU 3.6 GHz and 16 GB memory.
In the simulations, we assume that the image scene includes 64 × 64 pixels, and there are nine objects in the scene as shown in Figure 3a, where each object is represented by a pixel and the interval between two adjacent objects is a pixel. This is used to evaluate the performance of the algorithm in the case of small objects. We added the disturbance ∆S and error ∆E into the image and assume that ∆S is Gaussian distributed with mean 0 and variance 0.1, and error ∆E in (12) is also Gaussian distributed with mean 0 and variance 0.5 × C max . The maximum intensity of the original image is 255.
The sampling number for OMP, GPSR, Method proposed in [19] and STLS is 300, and the sampling number for TLS is 4500 because TLS requires M > N. It can be seen in Figure 3 that OMP and TLS exhibit the worst performances; GPSR and Method proposed in [19] work slightly better. We can also see that STLS outperforms other algorithms significantly. Simulations show that the PSNRs of OMP, GPSR, Method proposed in [19], TLS and STLS are 34.9166 dB, 37.4358 dB, 37.7993 dB, 33.2462 dB and 38.4405 dB, respectively, and the MSEs of OMP, GPSR, Method proposed in [19], TLS and STLS are 20.9612, 11.7356, 10.7932, 30.7939 and 9.3118, respectively. STLS achieves much better PSNR and MSE than OMP, GPSR, Method proposed in [19] and TLS schemes. The convergence of the STLS algorithm and the typical execution times are shown in Figure 4 and Table 1, respectively. Figure 4 illustrates the MSE of the reconstruction image with different numbers of iterations. Normally, the algorithm converges within 20 iterations. between two adjacent objects is a pixel. This is used to evaluate the performance of the algorithm in the case of small objects. We added the disturbance Δ and error Δ into the image and assume that Δ is Gaussian distributed with mean 0 and variance 0.1, and error Δ in (12) is also Gaussian distributed with mean 0 and variance 0.5 × . The maximum intensity of the original image is 255.  The sampling number for OMP, GPSR, Method proposed in [19] and STLS is 300, and the sampling number for TLS is 4500 because TLS requires M > N. It can be seen in Figure 3 that OMP and TLS exhibit the worst performances; GPSR and Method proposed in [19] work slightly better. We can also see that STLS outperforms other algorithms significantly. Simulations show that the   [19] and TLS schemes. The convergence of the STLS algorithm and the typical execution times are shown in Figure 4 and Table 1, respectively. Figure 4 illustrates the MSE of the reconstruction image with different numbers of iterations. Normally, the algorithm converges within 20 iterations.

Experimental Results and Discussions
Our quantum GI experimental system is shown in Figure 5. A continuous-wave laser with 460 nm wavelength was used to pump a BBO crystal, which was cut according to type-II collinear SPDC. The power of the pump laser is 300 mW, and the central wavelength and the bandwidth of the filter after BBO are 920 nm and 10 nm, respectively. The photon detection efficiency of single-photon detectors (SPCM-AQRH-FC, from Excelitas Technologies) at 920 nm is about 35%. In order to maximize the SPDC efficiency, we need to polarize the pump laser by a half-wave plate (HWP). Because only a small number of photons can be converted into entangled-photon pairs, we have to filter out the unconverted pump light by placing a filter behind BBO. A beam splitter (BS) was used to divide the entangled-photon pairs into the signal and reference arms. The entangled photons go through the object and are collected by the photon counting module (SPCM) in the object arm. The reference entangled photons are modulated by the SLM in the reference arm and are collected by the other SPCM. The SLM used in the experiment is HOLOEYE HES 6001-NIR with phase and amplitude type, and its resolution is 1920 × 1080 with pixel size 8 × 8 µm 2 . To demonstrate the performance of STLS, we use double-slit as the object with 64 × 64 pixels shown in Figure 6. The singles counts we collected from SPCM in the object arm and reference arm are 5.05 × 10 counts/s and 5.10 × 10 counts/s. In the experiment, the center-to-center distance of the double-slit is 1200 μm (22 pixels) and the pixel pitch is 54 μm. Figure 6 shows the reconstructed images by OMP, GPSR, Method proposed in [19] and STLS with sampling numbers of 500, 1000 and 1500, respectively. When the sampling number is 500, it is hard to identify the double-slit for OMP and GPSR. When the sampling number is 1000, we can see a blurry double slit in the images by OMP, GPSR and Method proposed in [19] but a clear one in the image by STLS, which demonstrates the advantage of the STLS scheme. To demonstrate the performance of STLS, we use double-slit as the object with 64 × 64 pixels shown in Figure 6. The singles counts we collected from SPCM in the object arm and reference arm are 5.05 × 10 4 counts/s and 5.10 × 10 4 counts/s. In the experiment, the center-to-center distance of the double-slit is 1200 µm (22 pixels) and the pixel pitch is 54 µm. Figure 6 shows the reconstructed images by OMP, GPSR, Method proposed in [19] and STLS with sampling numbers of 500, 1000 and 1500, respectively. When the sampling number is 500, it is hard to identify the double-slit for OMP and GPSR. When the sampling number is 1000, we can see a blurry double slit in the images by OMP, GPSR and Method proposed in [19] but a clear one in the image by STLS, which demonstrates the advantage of the STLS scheme.
PSNR and MSE were also used to evaluate the quality of the CS image reconstruction. Figures 7  and 8 show the MSE and PSNR of three schemes, where we can see that the MSE of STLS is much smaller than the other two schemes, and the PSNR of STLS is much higher than the other two schemes. We can see that the STLS scheme achieves better reconstruction results than others with the same number of measurements. As the number of samples increases, better recovery results can be obtained. The performance of our quantum ghost imaging scheme can be optimized by improving the coincidence rate and the beam quality of the pump laser. In the experiment, the object has 64 × 64 = 4096 pixels. We can obtain a very high-quality image with 1500 measurements and a clear image with 1000 measurements.
Sensors 2019, 19, 192 8 of 11 are 5.05 × 10 counts/s and 5.10 × 10 counts/s. In the experiment, the center-to-center distance of the double-slit is 1200 μm (22 pixels) and the pixel pitch is 54 μm. Figure 6 shows the reconstructed images by OMP, GPSR, Method proposed in [19] and STLS with sampling numbers of 500, 1000 and 1500, respectively. When the sampling number is 500, it is hard to identify the double-slit for OMP and GPSR. When the sampling number is 1000, we can see a blurry double slit in the images by OMP, GPSR and Method proposed in [19] but a clear one in the image by STLS, which demonstrates the advantage of the STLS scheme. PSNR and MSE were also used to evaluate the quality of the CS image reconstruction. Figure 7 and Figure 8 show the MSE and PSNR of three schemes, where we can see that the MSE of STLS is much smaller than the other two schemes, and the PSNR of STLS is much higher than the other two schemes. We can see that the STLS scheme achieves better reconstruction results than others with the same number of measurements. As the number of samples increases, better recovery results can be obtained. The performance of our quantum ghost imaging scheme can be optimized by improving the coincidence rate and the beam quality of the pump laser. In the experiment, the object has 64 × 64 = 4096 pixels. We can obtain a very high-quality image with 1500 measurements and a clear image with 1000 measurements.  PSNR and MSE were also used to evaluate the quality of the CS image reconstruction. Figure 7 and Figure 8 show the MSE and PSNR of three schemes, where we can see that the MSE of STLS is much smaller than the other two schemes, and the PSNR of STLS is much higher than the other two schemes. We can see that the STLS scheme achieves better reconstruction results than others with the same number of measurements. As the number of samples increases, better recovery results can be obtained. The performance of our quantum ghost imaging scheme can be optimized by improving the coincidence rate and the beam quality of the pump laser. In the experiment, the object has 64 × 64 = 4096 pixels. We can obtain a very high-quality image with 1500 measurements and a clear image with 1000 measurements.   [19] and STLS with different sampling numbers.  [19] and STLS with different sampling numbers.

Conclusions
In this work, we have investigated robust GI with a relatively small number of measurements to deal with the imperfection of optical devices, measurement error and noise in the optical path. The proposed method uses the pseudo-random matrix as the measurement matrix, and signals are transformed into a sparse domain by discrete wavelet transform. Then, the reconstruction is formulated as a sparse total least square problem which is solved iteratively. Both simulation and experimental results have been provided to show the superiority of the proposed STLS scheme, which can achieve significantly better PSNR and MSE than the system with other reconstruction algorithms such as TLS, OMP, GPSR and Method proposed in [19]. This work demonstrates the significant potential of handling the imperfection of optical devices in improving the quality of reconstruction. Figure 8. PSNR of the reconstructed quantum ghost imaging with OMP, GPSR, Method proposed in [19] and STLS with different sampling numbers.

Conclusions
In this work, we have investigated robust GI with a relatively small number of measurements to deal with the imperfection of optical devices, measurement error and noise in the optical path. The proposed method uses the pseudo-random matrix as the measurement matrix, and signals are transformed into a sparse domain by discrete wavelet transform. Then, the reconstruction is formulated as a sparse total least square problem which is solved iteratively. Both simulation and experimental results have been provided to show the superiority of the proposed STLS scheme, which can achieve significantly better PSNR and MSE than the system with other reconstruction algorithms such as TLS, OMP, GPSR and Method proposed in [19]. This work demonstrates the significant potential of handling the imperfection of optical devices in improving the quality of reconstruction.