Gerchberg-Saxton-like ghost imaging

Correlation is widely used to reconstruct the object image in ghost imaging (GI). But it only offers a linear proportion of the signal-to-noise ratios (SNR) to the number of measurements. We develop a Gerchberg-Saxton-like technique for GI image reconstruction in this manuscript. The proposed technique takes the advantage of the integral property of the Fourier transform, and treat the captured data as constraints for image reconstruction. We numerically and experimentally demonstrate the technique, and observe a nonlinear growth of the SNR value with respect to the number of measurements in the simulation. The proposed technique provides a different perspective of image reconstruction of GI, and will be beneficial to further explore its potential. © 2015 Optical Society of America OCIS codes: (110.1758) Computational imaging; (110.3010) Image reconstruction techniques; (100.5070) Phase retrieval. References and links 1. T. B. Pittman, Y. H. Shih, D. V. Strekalov, and A. V. Sergienko, “Optical imaging by means of two-photon quantum entanglement,” Phys. Rev. A 52, 3429 (1995). 2. A. Valencia, G. Scarcelli, M. DAngelo, and Y. Shih, “Two-photon ghost imaging with thermal light,” Phys. Rev. Lett 94, 063601 (2005). 3. W. Gong, P. Zhang, X. Shen, and S. Han, “Ghost ’pinhole’ imaging in Fraunhofer region,” Appl. Phys. Lett 95, 071110 (2009). 4. K. W. C. Chan, M. N. O’Sullivan, and R. W. Boyd, “High-order thermal ghost imaging,” Opt. Lett 34, 213343 (2009). 5. X. F. Liu, X. H. Chen, X. R. Yao, W. K .Yu, G. J. Zhai, and L. A. Wu, “Lensless ghost imaging with sunlight,” Opt. Lett39, 002314 (2014). 6. J. H. Shapiro, “Computational ghost imaging,” Phys. Rev. A 78, 061802R (2008). 7. Y. Bromberg, O. Katz, and Y. Silberberg, “Ghost imaging with a single detector,” Phys. Rev. A 78, 053840 (2009). 8. B. Sun, M. P. Edgar, R. Bowman, L. E. Vittert, S. Welsh, A. Bowman, and M. J. Padgett, “3D computational imaging with single-pixel detectors,” Science 340, 844–847 (2013). 9. C. Zhao, W. Gong, M. Che, E. Li, H. Wang, W. Xu, and S. Han, “Ghost imaging lidar via sparsity constraints,” Appl. Phys. Lett101, 141123 (2012). 10. P. Clemente, V. Duŕan, V. Torres-Company, E. Tajahuerce, and J. Lancis “Optical encryption based on computational ghost imaging,” Opt. Lett 35, 2391-2393 (2010). 11. F. Ferri, D. Magatti, L. A Lugiato, and A. Gatti, “Differential ghost imaging,” Phys. Rev. Lett 104, 253603 (2010). 12. O. Katz , Y. Bromberg, and Y. Silberberg, “Compressive ghost imaging,” Appl. Phys. Lett 95, 131110 (2009). 13. W. Wang, Y. P. Wang, Y. Wu, X. X. Yang, and Y. Wu, “Optical eigenmode imaging with a sparse constraint,” Opt. Lett39, 092614 (2014). 14. X. M. Hu, J. L. Suo, T. Yue, L. H. Bian, and Q. H. Dai, “Patch-primitive driven compressive ghost imaging,” Opt. Express23, 011092 (2015). #249703 Received 8 Sep 2015; revised 15 Oct 2015; accepted 15 Oct 2015; published 21 Oct 2015 © 2015 OSA 2 Nov 2015 | Vol. 23, No. 22 | DOI:10.1364/OE.23.028416 | OPTICS EXPRESS 28416 15. W. L. Gong and S. S. Han, “High-resolution far-field ghost imaging via sparsity constraint,” Sci. Rep 5, 9280 (2015). 16. W. Wang, Y. P. Wang, J. H. Li, X. X. Yang, and Y. Wu, “Iterative ghost imaging,” Opt. Lett 39, 175150 (2014). 17. X. R. Yao, W. K. Yu, X. F. Liu, L. Z. Li, L. A. Wu, and G. J. Zhai, “Iterative denoising of ghost imaging,” Opt. Express22, 024268 (2014). 18. R. J. Marks II,Handbook of Fourier Analysis and Its Applications (Oxford University Press, 2009). 19. R. Gerchberg and W. Saxton, “Phase determination for image and diffraction plane pictures in the electron microscope,” Optik34, 275-284 (1971). 20. G. Zheng, R. Horstmeyer, and C. Yang, “Wide-field, high-resolution Fourier Ptychographic microscopy,” Nat. Photonics7, 739-745 (2013). 21. S. Y. Dong, P. Nanda, R. Shiradkar, K. K. Guo, and G. A. Zheng, “High-resolution fluorescence imaging via pattern-illuminated Fourier ptychography,” Opt. Express. 22, 020856 (2014). 22. Y. S. G. Nashed, D. J. Vine, T. Peterka, J. J. Deng, R. Ross, and C. Jacobsen, “Parallel ptychographic reconstruction,” Opt. Express22, 32082-32097 (2014). 23. S. M. M. Khamoushi, Y. Nosrati, and S. H. Tavassoli, “Sinusoidal ghost imaging,” Opt. Lett 40, 003452 (2015). 24. M. Guizar-Sicairos and J. R. Fienup, “Phase retrieval with transverse translation diversity: a nonlinear optimization approach,” Opt. Express 16, 7264-7278 (2008). 25. L. Tian, X. Li, K. Ramchandran, and L. Waller, “Multiplexed coded illumination for Fourier Ptychography with an LED array microscope,” Biomed. Opt. Express 5 , 2376 2389 (2014). 26. F. Zhang, I. Peterson, J. Vila-Comamala, A. Diaz, F. Berenguer, R. Bean, B. Chen, A. Menzel, I. K. Robinson, and J. M. Rodenburg, “Translation position determination in ptychographic coherent diffraction imaging,” Opt. Express21, 13592-13606 (2013). 27. A. Alpers, G. T. Herman, H. F. Poulsen, and S. Schmidt, “Phase retrieval for superposed signals from multiple binary objects,” J. Opt. Soc. Am. A 27, 1927-1937 (2010).


Introduction
Ghost imaging (GI) is a novel imaging technique that allows the imaging through harsh environments [1][2][3][4][5][6][7][8], with great potential in applications in remote sensing [9] and optical encryption [10], etc. In GI, the camera does not need to "see" the object to be imaged directly. Alternatively, the image of the object is reconstructed by correlating two beams. One of the beams illuminates the object and its energy was recorded by a single pixel detector with no spatial resolution (i.e., a bucket detector). The other beam is captured using a multipixel camera. Actually, this camera can be removed in computational GI as the fields that produces the illumination patterns can be computer-programmed [6,7]. However, this correlation-based reconstruction technique only offers a signal-to-noise ratio (SNR) proportional to the number of measurements [11]. Recently, a lot of efforts have been taken to increase the SNR for a given number of measurements. Three of the most important methods are: (1) Differential ghost imaging (D-GI) [11]. This technique offers a considerable improvement of the SNR. But its performance is object-dependent. (2) Compressive ghost imaging. This method is based on the compressive sensing technique [12][13][14][15]. It only works for sparse objects, and the optimization procedure is time-consuming. (3) Iterative ghost imaging. This method employs the image transmission matrix to recover the image iteratively using GI or DGI [16,17]. It has been demonstrated that the iterative GI has a better performance than the DGI in terms of SNR.
In this article, we demonstrate a Gerchberg-Saxton-like technique for GI image reconstruction by taking the integral property of the Fourier transform into account. According to this property [18], we can regard the data recorded by the bucket detector as the Fourier transform of the object image evaluated at the origin (k = 0). Note that one needs to apply a sequence of two-dimensional speckles to illuminate the object. Each of these speckle patterns randomly select certain spectral components of the object and shift them to the origin in the Fourier space. One can thus use this as a priori information in reconstructing the object image. Figure 1 shows the schematic setup we used in simulation and experiments. The unknown transmission object, T (x), is illuminated by a sequence of speckles, I m (x), where the subscript integer (m = 1 . . . M) denotes the m th illumination. We define the amplitude immediately behind the object, I tm (x) = I m (x)T (x), as the target image. Then, the signal collected by the bucket detector can be written as

Principle of the proposed technique and simulation results
The sequence [S 1 , S 2 , . . . , S M ] is the data that is experimentally captured in the computational GI system [6,7]. One of the most important tasks of GI is to reconstruct the dimensional object, Traditionally, the reconstruction is obtained by correlating the data [S m ] with the intensity I m (x) of speckle patterns as where · ≡ 1 M ∑ m denotes an ensemble average over M measurements, δ S m = S m − S m , δ I m = I m − I m . According to the correlation theory, the recovered image O(x) is approximately equal to the convolution between T (x) and the correlation function γ(x) of the speckle fields. The resolution has the scale of the speckle size δ x which is decided by the width of γ(x). Therefore, the object plane is divided into N speckle resolution units each with the area of δ x 2 [11]. One of the biggest issues with this method is that the SNR value of the reconstructed image is statistically proportional to the ratio of the number of patterns M to the number of speckles N speckle . In the case of M = N speckle , the SNR value is only about 1.
Here we propose an alternative technique to reconstruct O(x) from the experimentally recorded data sequence [S m ], (m = 1, . . . , M), rather than solving Eq. (2). In the proposed technique, we treat the image reconstruction problem as the estimation of O(x) with the knowledge of the speckle patterns I m (x), so that the error between I m (x)O(x)dx and the experimental data S m is minimized. According to the integral property of the Fourier transform [18], we have where . This inspires us to develop a Gerchberg-Saxton-like algorithm with the series of constraints [S m ] sequentially imposed on the Fourier spectra of the corresponding target image I tm (x). The initial guess of the object is chosen as the image reconstructed using conventional GI as in Eq. (2). Assuming the last guess of the object is written as O(x), the iteration routine of the proposed algorithm is as follows: Step 1. Multiplication of O(x) with I m (x) produces an estimation of the target image I tm (x) = I m (x)O(x). Perform the Fourier transform of I tm (x), and apply the constraint Eq. (3) to the resulting spectrum.
Step 2. Inverse Fourier transform the modified Fourier spectrum back to the spatial domain and produce an updated target image, I updated tm (x). Then the current estimation of the object can be written as: Step 3. Replace O by O updated , and repeat the above two steps for all the M measurements again and again until the iteration converges. The convergence can be determined by whether the mean-square error of two consecutively solutions is smaller than a pre-defined threshold value (which is 10 −10 in our study). Now we perform numerical simulations for an object with the size of 64 * 64 pixels to demonstrate the proposed technique. The reconstructed images using the traditional correlation technique and the proposed technique are shown in the top and bottom rows, respectively, in Fig. 2. It clearly suggests that the noise decreases as M increases no matter which methods are used. However, the noise in the images recovered using the proposed technique is significantly suppressed comparing to GI. To quantify the error between the recovered and the ground-truth image (shown in the inset), we define the SNR as [17] whereT is the mean of T (x). The SNR value of the image reconstructed by solving Eq. (2) is increased from 1.05 to 1.84 as M increases from 2000 to 5000. However, this value can be increased from 1.65 to as high as 66.58 by using the proposed method in our study, meaning that the noise has been significantly reduced, as shown in the bottom row of Fig. 2.  It is well known that the SNR value of the recovered image is linearly proportional to M when using correlation techniques. However, Figure 2 suggests that this linear relationship is broken by using the proposed method. To be more specific, we plot the SNR value as the function of M in Fig. 3. One can clearly see the nonlinear growth of the SNR value of the image reconstructed using the proposed method in contrast to the linear behavior for GI. We fit the curve with an exponential function, and found the R 2 factor equal to 0.995.
We also observed that the image quality can be close to the ground-truth image when M > N speckle . This can be understood as follows: The effect of the speckle illumination is essentially to shift randomly-selected spectral components of the object O(x) to the origin in the Fourier domain. As a consequence, these randomly selected spectral components can be captured by the bucket detector, and encoded in S m . This means that the sequence S m , (m = 1, . . . , M), acquired by the bucket detector contains the information about the spectrum of the object in a statistics manner. For an object of N × N pixels in size, its spectrum also has N speckle = N 2 pixels. Statistically speaking, with N speckle measurements, there is a large possibility that one can obtain sufficient information about the spectrum of the object. The more number of speckle patterns one may use, the more information about the object spectrum one can collect, resulting in a higher SNR value for the reconstructed image, as demonstrated by the red curve in Fig. 3.

Experimental demonstration and Discussion
We demonstrate the proposed method by using the experimental setup as shown in Fig. 1. The light emitted from an incoherent light source (Philips, 200 W) was collimated to illuminate the digital mirror device (DMD, DLP Discovery-4100, TI). Then, a sequence of binary random patterns was generated by using a computer program and loaded to the DMD. A projector was used to project the binary patterns onto the object plane. Finally, a bucket detector (Thorlabs Det100 silicon photodiode) was used to record the energy of the target image. The detector recorded a sequence of dimensionless energy values S m , m = 1, . . . , M, while changing the illumination pattern. The frame rate of the DMD was set to be 100 Hz in our experiment while the acquisition frequency of the was 10000 Hz. This means that the detector took 100 measurements for each speckle pattern so that the noise can be averaged out. To avoid the influence of the elapse caused by either the DMD or the detector, the final value of S m was set to the average value of the intermediate 90 measurements. The object used in the experiment was the 'ghost' image shown in Fig. 4(a). GI reconstructs the transmission image of the object using the experimental setup. In our experiment, we generated 4096 binary patterns and illuminated the object. The single-pixel detector then recorded a sequence of 4096 energy values S m , (m = 1, . . . 4096). This in principle allows us to reconstruct an image of the size 64 × 64 with an SNR value equal to 1 using the standard GI. Experimentally, we found that the SNR value of the image [see Fig. 4(b)] recovered by solving Eq. (2) equal to 0.913, which is close to the theoretical value 1. In contrast, using the proposed method, the SNR value of the reconstructed image shown in Fig. 4(d) is 2.994, which is 2 times larger than that of Fig. 4(b). This is expectable and in consistence with the simulation results in Fig. 2. We observed that the SNR value of Fig. 4(d) is even larger than that of Fig. 4(c), which was reconstructed using DGI [11] . It is expected that the SNR value can be further increased by measuring the actual speckle patterns on the object plane, or increasing the number of patterns to be more than the critical number N speckle . It was demonstrated in the simulation that the SNR of the reconstructed image reaches 27.74 for the 'ghost' image when M is set to 5000.
The convergence of the proposed method is shown in Fig. 5. The number of iteration is denoted as N ite . When it reaches 50 − 100, the SNR converges to the optimum value. The fast-convergence property is in consistent with other Gerchberg-Saxton algorithms [19][20][21]. Moreover, the computational time can be further decreased by using a graphical processing unit (GPU) [20,22] .
The proposed algorithm has no requirement on the illumination pattern. Thus, any arbitrary illumination can be generated as a priori information in the method. Structured light, such as sinusoidal patterns [23], may provide better results in an ideal system. In the article, random patterns are chosen as the illumination sources because the optical system is simple and the acquisition time is fast using the DMD. Moreover, random patterns are nature in harsh environments.

Conclusion
In conclusion, we have numerically and experimentally demonstrated an efficient method for GI image reconstruction. The proposed method treats the recorded data sequence S m in a way significantly different from the perspective of the widely-used correlation technique: It takes S m as the constraints to be imposed on the origin in the Fourier domain. The proposed method is similar to the Gerchberg-Saxton algorithm [19] and ptychography [20,21,[24][25][26][27], except that the constraints here are sequentially applied to a single point during iteration. Our method does not require any change to the standard GI or computational GI acquisition system, but recover the object images with better quality (exponential growth of the SNR value along with the increase of M). This will be benefitial to further explore the potential of GI.