Iterative denoising of ghost imaging

Xu-Ri Yao; Wen-Kai Yu; Xue-Feng Liu; Long-Zhen Li; Ming-Fei Li; Ling-An Wu; Guang-Jie Zhai

doi:10.1364/OE.22.024268

1. Introduction

Due to its features of nonlocality and simple experimental configuration, ghost imaging (GI) has attracted wide attention in the last decade. It was first observed by using entangled photon pairs [1], then later it was shown both theoretically and experimentally that thermal light can also be used to realize GI [2, 3, 4, 5, 6, 7, 8, 9]. This promise displays great potential with respect to conventional imaging techniques because it allows imaging in scenarios where using a detector array is restricted or difficult. Several GI configurations have been proposed and demonstrated for imaging in harsh environments [10, 11, 12, 13], optical coherence tomography [14] and optical encryption [15, 16, 17, 18]; however, the poor visibility and signal-to-noise ratio (SNR) of the images have restricted practical applications. To address these problems, various protocols have been suggested and demonstrated, among which differential ghost imaging (DGI) [19] and normalized ghost imaging (NGI) [20] show remarkable enhancement; the new time-correspondence schemes [21, 22, 23] not only improve in enhancing the SNR but also reduce the computing time. Compressed sensing combined with GI [24, 25, 26] can greatly improve the image quality for a given number of measurements though at the cost of more computation time. Another work that appeared while our paper was under review also proposed a method to improve the SNR using an iterative algorithm, and gave a numerical simulation [27].

In this paper, we present a different approach that we call iterative denoising of ghost imaging (IDGI) to pinpoint then remove the noise in GI. We first discuss the cause and effect of the noise, then describe our method to obtain an optimized approximation to this noise. A threshold is utilized in the iteration to reduce the effects of too large speckle areas. It is shown that with this method the image quality can be improved greatly, even more so than by DGI. This represents another step forward towards real practical applications.

2. Iterative ghost imaging algorithms

2.1. Noise in ghost imaging

The optical setup for lensless GI with thermal light (for convenience a pseudothermal light source is often used [3]) is shown schematically in Fig. 1. In GI, two detectors are used, a bucket detector which only measures the total light intensity transmitted/reflected by an object, and a spatially resolving reference detector such as a charge-coupled device (CCD) camera in another beam reflected/transmitted by from a beamsplitter placed before the object. Neither detector alone can image the object, but an image can be retrieved by intensity correlation measurements between the two detectors. In conventional GI, the background-subtracted correlation function is [28]

Δ G^{(2)} (x) = \frac{1}{M} \sum_{i = 1}^{M} I_{B}^{(i)} I_{R}^{(i)} (x) - \frac{1}{M} \sum_{i = 1}^{M} I_{B}^{(i)} \frac{1}{M} \sum_{i = 1}^{M} I_{R}^{(i)} (x) .

Here

I_{R}^{(i)} (x)

is the intensity distribution on the pixel array detector for the ith measurement, M denotes the total number of the measurements, and

I_{B}^{(i)} = \sum_{x} I_{R}^{(i)} (x) O (x)

is the total intensity collected by the bucket detector placed behind the object in the ith measurement, where O(x) is the transmission function of the object. It can be shown that

\begin{array}{l} \frac{1}{M} \sum_{i = 1}^{M} I_{B}^{(i)} I_{R}^{(i)} (x) & = \frac{1}{M} \sum_{i = 1}^{M} \sum_{x^{'}}^{N} O (x^{'}) (I_{R}^{(i)} (x^{'}) - \bar{I} (x^{'})) (I_{R}^{(i)} (x) - \bar{I} (x)) + \bar{I} (x) \sum_{x^{'}}^{N} \bar{I} (x^{'}) O (x^{'}) \\ = σ^{2} (x, x) O (x) + \sum_{x^{'} \neq x} σ^{2} (x, x^{'}) O (x^{'}) + \bar{I} (x) \sum_{x^{'}}^{N} \bar{I} (x^{'}) O (x^{'}) \end{array}

and

\frac{1}{M} \sum_{i = 1}^{M} I_{B}^{(i)} \frac{1}{M} \sum_{i = 1}^{M} I_{R}^{(i)} (x) = \bar{I} (x) \sum_{x^{'}}^{N} \bar{I} (x^{'}) O (x^{'}),

where σ²(x, x′) is the covariance between the intensities at x and x′,

\bar{I} (x) = \frac{1}{M} \sum_{i = 1}^{M} I_{R}^{(i)} (x)

is the mean of the intensity at pixel x, and N denotes the total number of pixels. Assuming that the intensity at each pixel is independent of that at every other pixel, i.e., σ²(x, x′) = 0 (x ≠ x′), and neglecting other noise sources for the time being [28], then for a pseudothermal light source we should have

Δ G^{(2)} (x) = σ^{2} O (x)

where σ²(x, x) has been abbreviated to σ² for short.

Fig. 1 Schematic of ghost imaging. BS: beamsplitter; CCD1: CCD camera acting as a bucket detector; CCD2: CCD reference detector.

Download Full Size | PDF

Unfortunately, in actual situations, various factors such as the transverse coherence length of the light source and the finite number of measurements prevent the covariance σ²(x, x′) from being zero, thus giving rise to additional noise. This term, which represents the actual noise, must therefore be included and added to Eq. (4), but to avoid confusion with the variance term σ²(x, x) we use n(x, x′) to denote σ²(x, x′):

Δ G^{(2)} (x) = σ^{2} O (x) + \sum_{x^{'} \neq x} n (x, x^{'}) O (x^{'}) .

Though the value of n(x, x′) is usually very small, the total noise term

\sum_{x^{'} \neq x} n (x, x^{'}) O (x^{'})

may be non-negligible and can drown the reconstructed image if the object has a large transmission ratio or low contrast. Conventional GI is far from practical application since an enormous number of measurements have to be taken to overcome the noise [29].

2.2. Iterative ghost imaging algorithms

To help understand the purpose of our iterative method, we start with the definition of DGI, which was proposed and demonstrated by Ferri et al [19]:

\begin{matrix} DGI (x) = \frac{1}{M} \sum_{i = 1}^{M} I_{B}^{(i)} I_{R}^{(i)} (x) - \frac{\bar{I_{B}}}{I_{R}} \frac{1}{M} \sum_{i = 1}^{M} (I_{R}^{(i)} (x) \sum_{x}^{N} I_{R}^{(i)} (x)) \\ = Δ G^{(2)} (x) - \frac{\bar{I_{B}}}{I_{R}} (\frac{1}{M} \sum_{i = 1}^{M} (I_{R}^{(i)} (x) \sum_{x}^{N} I_{R}^{(i)} (x)) - I_{R} \frac{1}{M} \sum_{i = 1}^{M} I_{R}^{(i)} (x)) \\ = σ^{2} O (x) + \sum_{x^{'} \neq x} n (x, x^{'}) O (x^{'}) - T \sum_{x^{'}} n (x, x^{'}) \end{matrix}

Here $\bar{I_{B}} = \frac{1}{M} \sum_{i = 1}^{M} I_{B}^{(i)}$ , $I_{R} = \frac{1}{M} \sum_{i = 1}^{M} \sum_{x}^{N} I_{R}^{(i)} (x)$ , and $T = \frac{1}{N} \sum_{x} O (x)$ denotes the transmission ratio of the object. Comparing Eq. (6) with Eq. (5) of conventional GI, we see that there is a variable quantity $T \sum_{x^{'}} n (x, x^{'})$ which is approximately equal to the noise $\sum_{x^{'} \neq x} n (x, x^{'}) O (x^{'})$ and is subtracted at each pixel. Thus from the viewpoint of denoising, DGI can also be regarded as a process of noise subtraction.

As seen in Eq. (6), the noise term $T \sum_{x^{'}} n (x, x^{'})$ depends on the transmission ratio T of the object, while the real noise $\sum_{x^{'} \neq x} n (x, x^{'}) O (x^{'})$ comes from the object’s transmission function O(x). Differential GI shows dramatic improvement when the object is nearly transparent (T ≈ 1) or opaque (T ≈ 0) because in these circumstances the assumed value of the noise is close to the real noise. However, without prior knowledge of the object, DGI may not be so successful, whereas with our IDGI algorithm we use the initial image retrieved within the context of GI instead of the transmission ratio of the object T. This will give a more accurate estimate of the noise, and a better image IDGI(x) can be obtained, as follows:

IDGI (x) = σ^{2} O (x) + \sum_{x^{'} \neq x} n (x, x^{'}) O (x^{'}) - \sum_{x^{'} \neq x} n (x, x^{'}) O (x^{'}),

where O′ (x) is the image retrieved by some conventional algorithm. This improved image can then be used to give a better estimate of the noise, which is again used to obtain a better image, and so on. This iterative process can be repeated k times until the SNR reaches a maximum limit to give the final IDGI image

\begin{matrix} {IDGI}^{(k + 1)} (x) = σ^{2} O (x) + \sum_{x^{'} \neq x} n (x, x^{'}) O (x^{'}) - \sum_{x^{'} \neq x} n^{'} (x, x^{'}) {IDGI}^{(k)} (x^{'}), \\ n^{'} (x, x^{'}) = n (x, x^{'}), n (x, x^{'}) \leq t; n^{'} (x, x^{'}) = 0, n (x, x^{'}) > t . \end{matrix}

Numerical simulation and actual experiments were performed to demonstrate the efficacy of our IDGI algorithm. The initial value of IDGI used was the image reconstructed by DGI. As O(x) varies over the range [0, 1], the image used as the initial value was approximately normalized to [0, 1] by subtracting the mean of the intensities from those pixels where no light was transmitted by the object, and then divided by the maximum pixel intensity value. A certain threshold t, which will be discussed later, was used to set the zero value in some positions in n′ (x, x′). Generally, the SNR of IDGI was found to reach its maximum with an iteration number of just k = 3.

To illustrate the basic concept underlying IDGI, we consider the effects of the error in image estimation. For convenience, we take a binary object as shown in Fig. 2(a), which is totally opaque except for a square in the center. We assume that the average size of the speckles of the pseudothermal light on the object plane is larger than the size of one pixel of the CCD camera. For a pixel x₀ located on the reference camera corresponding to some spot within the opaque area of the object, the distribution of its covariance σ²(x₀, x′) with neighboring pixels exhibits a Gaussian-like distribution (Fig. 2(b)). However, as seen in Eq. (5), the noise at x₀ is $\sum_{x^{'} \neq x} n (x_{0}, x^{'}) O (x^{'})$ (n(x₀, x′) = σ²(x₀, x′)), so the fluctuations of the σ²(x₀, x′) that are transmitted through the transparent square is the main factor that contributes to the noise at that pixel (Fig. 2(c) and Fig. 2(d)). As IDGI(x′) is not completely identical with the object O(x′) and may not be zero at the pixels neighboring x₀, then the relatively large values of σ²(x₀, x′) at these pixels will result in considerable difference between the estimated noise $\sum_{x^{'} \neq x_{0}} n (x_{0}, x^{'}) IDGI (x^{'})$ and the real noise $\sum_{x^{'} \neq x_{0}} n (x_{0}, x^{'}) O (x^{'})$ . Thus we use a threshold t to set the large values in this area adjoining x₀ to be zero. In this way, we reduce the errors resulting from the blurring due to too large speckles in the transverse field, and obtain a more accurate estimated noise. On the other hand, when the pixel x′₀ is located at the edge of or inside the transparent area, as shown in Fig. 2(b1), the situation is reversed. From Figs. 2(c1) and 2(d), it can be seen that in conventional GI the large values within the area neighboring x′₀ contribute to the noise at x′₀. If all these large values are set to zero, the estimated noise $\sum_{x^{'} \neq x_{0}} n (x_{0}, x^{'}) IDGI (x^{'})$ will no longer be accurate in the transparent area, and as a result, the retrieved image will have blurry edges. Therefore, an appropriate threshold is needed to eliminate this. For an object with a continuously varying gray scale, large speckles will affect the areas with low or high transmittance differently, just as with a binary object, so an appropriate threshold is needed as well.

Fig. 2 (a) The binary object which is opaque except for an transparent square at the center. (b) The distribution of σ²(x₀, x′) in which x₀ is placed far from the transparent area of the object. (c) Image obtained by multiplying (a) with (b). (b1) The distribution of σ²(x′₀, x′) in which x′₀ is placed at the edge of the transparent rectangle. (c1) Image obtained by multiplying (a) with (b1). (d) The conventional GI image, the value at x₀ is the sum of the values in the square in (c), while the value at x′₀ is the sum of the values in the square in (c1).

Download Full Size | PDF

A numerical simulation has been performed to show the relation between the SNR_image and the threshold t, where SNR_image is defined as [22]:

{SNR}_{image} = \frac{\sum_{x}^{N} {[O (x) - T]}^{2}}{\sum_{x}^{N} {[R (x) - O (x)]}^{2}},

in which R(x) is the retrieved image. A digital object “NSSC” of size 128×128 pixels with a continuously varying gray scale is shown in Fig. 3(a). As the threshold t varies over the range [min{n(x, x′)}, max{n(x, x′)}], to simplify the analysis, we use the normalized threshold

t^{'} = \frac{t - \min {n (x, x^{'})}}{\max {n (x, x^{'})} - \min {n (x, x^{'})}}

. From Fig. 3(b), it can be seen that the image retrieved with a low threshold has a very smooth background but blurry edges, while the image in Fig. 3(c), retrieved with a high threshold, has a fluctuating background but sharp edges. The image retrieved with the optimal threshold t′ = 0.5 is shown in Fig. 3(d). These results are consistent with our analysis above. To see the change in value of the SNR more clearly, we plot SNR_image versus t in Fig. 3(e). We can see that it reaches a maximum when the threshold t′ = 0.5, i.e., t is approximately equal to the mean value of n(x, x′), thus this mean value is used as the threshold in all the simulations and experiments below.

Fig. 3 Dependence of retrieved image on threshold. (a) The object. (b) Image retrieved with t′ = 0.1. (c) Image retrieved with t′ = 0.9. (d) Image retrieved with t′ = 0.5. (e) SNR_image vs. normalized threshold t′.

Download Full Size | PDF

Figure 4 shows the numerical simulation results for the same digital object “NSSC” (Fig. 4(a)). The speckle patterns with a negative exponential probability distribution were produced artificially and the average speckle size was about 3 times the size of a pixel. The number of measurements was 5 × 10⁴, and white Gaussian noise was randomly superimposed upon the bucket detector signal to give an SNR of SNR_signal = 21.49 dB, where

{SNR}_{signal} = 20 log (\frac{σ (signal)}{σ (noise)}),

in which σ is the standard deviation. The images reconstructed by conventional GI and DGI are shown in Figs. 4(b) and 4(c), respectively, and the image retrieved by IDGI after 3 iterations is shown in Fig. 4(d). To provide a quantitative comparison of the image quality obtained by these different methods, all the images were normalized to [0, 1]. We find that the SNR_image of GI is 1.32 and that of DGI is 2.11, while for IDGI it is 4.52, a dramatic enhancement.

Fig. 4 Simulated results for GI, DGI and IDGI, averaged over 50,000 frames. (a) The object. (b) GI image, with SNR_image = 1.32. (c) DGI image, with SNR_image = 2.11. (d) IDGI image from k = 3 iterations, with SNR_image = 4.52.

Download Full Size | PDF

3. Experiment

An experiment was performed with a real light source in the setup as shown in Fig. 1, but with a digital object [22, 23]. The pseudothermal source was formed by passing a He-Ne laser beam of wavelength 632.8 nm through a slowly rotating ground glass disk. The light beam was divided by a 50:50 beam splitter(BS) into two spatially correlated beams, the object and reference beams, and then detected by the two CCD cameras, CCD1 and CCD2. The object and CCD2 were symmetrically placed with respect to BS, namely, Z₁ = Z₂ = 190 mm, and CCD1 just behind the object played the role of a bucket detector. The digital mask of Lena, which is widely used as a standard in imaging processing, of 64 × 64 pixel size is shown in Fig. 5(a). The number of exposure frames was 1.8 × 10⁴. The results of different intensity correlation methods are shown in Figs. 5(b), 5(c), and 5(d), respectively. We can see that IDGI with an iteration number of 3 produced the maximum SNR_image of 6.93.

Fig. 5 Experimental results for GI, DGI and IDGI. (a) Digital mask. (b) GI image, with SNR_image = 1.39. (c) DGI image, with SNR_image = 2.45. (d) IDGI image, with SNR_image = 6.93.

Download Full Size | PDF

4. Conclusion

In conclusion, we have analyzed the influence of noise in GI, and presented a theoretical analysis together with an experimental demonstration of a new technique which we call iterative ghost imaging. Through alternately estimating the noise and updating the retrieved image, we can gradually approach the real noise values, and thus finally obtain an optimized image. An appropriate threshold is applied in the iteration to eliminate the detrimental influence of coherence. Our results indicate that IDGI gives a much better performance than DGI even with a few number of iterations, although of course the process of iteration will require somewhat more computation time. This new protocol should be very effective in practical applications of intensity correlation imaging.

Acknowledgments

This work was supported by the National Major Scientific Instruments Development Project of China (Grant No. 2013YQ030595), the Hi-Tech Research and Development Program of China (Grant No. 2013AA122902), and the National Basic Research Program of China (Grant No. 2010CB922904).

References and links

1. T. B. Pittman, Y. H. Shih, D. V. Strekalov, and A. V. Sergienko, “Observation of two-photon ’ghost’ interference and diffraction,” Phys. Rev. A 52, R3429 (1995). [CrossRef]

2. R. S. Bennink, S. J. Bentley, and R. W. Boyd, “’Two-photon” coincidence imaging with a classical source,” Phys. Rev. Lett. 89, 113601 (2002). [CrossRef]

3. A. Gatti, E. Brambilla, M. Bache, and L. A. Lugiato, “Ghost imaging with thermal light: comparing entanglement and classical correlation,” Phys. Rev. Lett. 93, 093602 (2004). [CrossRef] [PubMed]

4. J. Cheng and S. S. Han, “Incoherent Coincidence Imaging and Its Applicability in X-ray Diffraction,” Phys. Rev. Lett. 92, 093903 (2004). [CrossRef] [PubMed]

5. D. Z. Cao, J. Xiong, S. H. Zhang, L. F. Lin, L. Gao, and K. G. Wang, “Enhancing visibility and resolution in Nth-order intensity correlation of thermal light,” Appl. Phys. Lett. 92, 201102 (2008). [CrossRef]

6. Y. J. Cai and S. Y. Zhu, “Ghost imaging with incoherent and partially coherent light radiation,” Phys. Rev. E 71, 056607 (2005). [CrossRef]

7. D. Zhang, Y. H. Zhai, L. A. Wu, and X. H. Chen, “Correlated two-photon imaging with true thermal light,” Opt. Lett. 30, 2354 (2005). [CrossRef] [PubMed]

8. X. H. Chen, Q. Liu, K. H. Luo, and L. A. Wu, “Lensless ghost imaging with true thermal light,” Opt. Lett. 34, 695 (2009). [CrossRef] [PubMed]

9. K. W. C. Chan, M. N. O’Sullivan, and R. W. Boyd, “High-order thermal ghost imaging,” Opt. Lett. 34, 3343 (2009). [CrossRef] [PubMed]

10. J. Cheng, “Ghost imaging through turbulent atmosphere,” Opt. Express 17, 7916–7921 (2009). [CrossRef] [PubMed]

11. R. E. Meyers, K. S. Deacon, and Y. Shih, “Turbulence-free ghost imaging,” Appl. Phys. Lett. 98, 111115 (2011). [CrossRef]

12. N. D. Hardy and J. H. Shapiro, “Reflective ghost imaging through turbulence,” Phys. Rev. A 84, 063824 (2011). [CrossRef]

13. X. F. Liu, M. F. Li, X. R. Yao, W. K. Yu, G. J. Zhai, and L. A. Wu, “High-visibility ghost imaging from artificially generated non-Gaussian intensity fuctuations,” AIP Advances 3, 052121 (2013). [CrossRef]

14. X. F. Liu, X. R. Yao, X. H. Chen, L. A. Wu, and G. J. Zhai, “Thermal light optical coherence tomography for transmissive objects,,” J. Opt. Soc. Am. A 29, 1922 (2012). [CrossRef]

15. P. Clemente, V. Durán, V. Torres-Company, E. Tajahuerce, and J. Lancis, “Optical encryption based on computational ghost imaging,” Opt. Lett. 35, 2391 (2010). [CrossRef] [PubMed]

16. M. Tanha, R. Kherdmand, and S. Ahmadi-Kandjani, “Gray-scale and color optical encryption based on computational ghost imaging,” Appl. Phys. Lett. 101, 101108 (2012). [CrossRef]

17. S. Li, X. R. Yao, W. K. Yu, L. A. Wu, and G. J. Zhai, “High-speed secure key distribution over an optical network based on computational correlation imaging,” Opt. Lett. 38, 2144 (2013). [CrossRef] [PubMed]

18. W. K. Yu, S. Li, X. R. Yao, X. F. Liu, L. A. Wu, and G. J. Zhai, “Protocol based on compressed sensing for high-speed authentication and cryptographic key distribution over a multiparty optical network,” Appl. Opt. 52, 7882 (2013). [CrossRef]

19. F. Ferri, D. Magatti, L. A. Lugiato, and A. Gatti, “Differential ghost imaging,” Phys. Rev. Lett. 104, 253603 (2010). [CrossRef] [PubMed]

20. B. Sun, S. S. Welsh, P. Edgar, J. H. Shapiro, and M. J. Padgett, “Normalized ghost imaging,” Opt. Express 20, 16892 (2012). [CrossRef]

21. K. H. Luo, B. Q. Huang, W. M. Zheng, and L. A. Wu, ”Nonlocal Imaging by Conditional Averaging of Random Reference Measurements,” Chin. Phys. Lett. 29, 074216 (2012). [CrossRef]

22. M. F. Li, Y. R. Zhang, K. H. Luo, L. A. Wu, and H. Fan, “Time-correspondence differential ghost imaging,” Phys. Rev. A 87, 033813 (2013). [CrossRef]

23. M. F. Li, Y. R. Zhang, X. F. Liu, X. R. Yao, K. H. Luo, H. Fan, and L. A. Wu, “A double-threshold technique for fast time-correspondence imaging,” Appl. Phys. Lett. 103, 211119 (2013). [CrossRef]

24. O. Katz, Y. Bromberg, and Y. Silberberg, “Compressive ghost imaging,” Appl. Phys. Lett. 95, 131110 (2009). [CrossRef]

25. J. Du, W. l. Gong, and S. S. Han, “The influence of sparsity property of images on ghost imaging with thermal light,” Opt. Lett. 37, 1067 (2012). [CrossRef] [PubMed]

26. W. K. Yu, M. F. Li, X. R. Yao, X. F. Liu, L. A. Wu, and G. J. Zhai, “Adaptive compressive ghost imaging based on wavelet trees and sparse representation,” Opt. Express 22, 7133 (2014). [CrossRef] [PubMed]

27. W. Wang, Y. P. Wang, J. Li, X. Yang, and Y. Wu, “Iterative ghost imaging,” Opt. Lett. 39, 5150 (2014). [CrossRef] [PubMed]

28. P. Zerom, Z. Shi, M. N. OSullivan, K. W. C. Chan, M. Krogstad, J. H. Shapiro, and R. W. Boyd, “Thermal ghost imaging with averaged speckle patterns,” Phys. Rev. A 86, 063817 (2012). [CrossRef]

29. B. I. Erkmen and J. H. Shapiro, “Signal-to-noise ratio of Gaussian-state ghost imaging,” Phys. Rev. A 79, 023833 (2009). [CrossRef]

Iterative denoising of ghost imaging

Abstract

1. Introduction

2. Iterative ghost imaging algorithms

2.1. Noise in ghost imaging

2.2. Iterative ghost imaging algorithms

3. Experiment

4. Conclusion

Acknowledgments

References and links

Cited By

Figures (5)

Equations (10)

Optics Express