We present an experimental comparison between different iterative ghost imaging algorithms. Our experimental setup utilizes a spatial light modulator for generating known random light fields to illuminate a partially-transmissive object. We adapt the weighting factor used in the traditional ghost imaging algorithm to account for changes in the efficiency of the generated light field. We show that our normalized weighting algorithm can match the performance of differential ghost imaging.
© 2012 Optical Society of America
Classical ghost imaging (GI) [1–5] uses a series of random light patterns to illuminate an unknown object. For each pattern the reflected or transmitted light is measured using a single element detector. The series of single element measurements, combined with the known light patterns is used to deduce the object. In some systems the random light pattern is produced as a time varying laser speckle, and a beam splitter is used to illuminate both the unknown object and a reference camera, with which the pattern is recorded. Subsequently, the need for the beam splitter and camera has been removed by implementing a spatial light modulator (SLM) to produce a random, but known, pattern thereby reducing the number of components in the system necessary for GI experiments [6, 7]. This latter approach is known as computational GI and in terms of the experimental arrangement is closely related to the field of single pixel cameras .
In all approaches to GI an algorithm is employed to deduce the object using the series of measurements from the single element detector and either the recorded or computationally predicted random patterns. The algorithms employed fall into two categories, iterative ones that give a refined estimate of the object after every new light pattern and measurement, and inversion ones which infer an object based on the entire series of patterns and measurements.
Iterative algorithms use the measured signal to derive a weighting factor to the corresponding pattern that is then added to the iterative estimate of the object. In this paper we compare a number of these iterative algorithms within the context of computational GI. The algorithms we consider are traditional GI (TGI) and differential GI (DGI) . In a computational GI setup, TGI uses a weighting factor equal to the signal from the detector whereas DGI utilizes a weighting factor that depends on fluctuations in the measured signal and uses an additional detector to give a normalization. Beyond these two algorithms we introduce a variant of the TGI algorithm, normalized GI (NGI), which we show can match the performance of DGI.
Key to all these algorithms is that the changes in the measured signal should arise from the overlap of the known random pattern with the unknown object. Obviously other sources of signal change are possible; including fluctuations arising from changes in the source intensity and changes in the efficiency with which the pattern is imprinted. These later sources of noise scale with the signal level and hence become more significant when the signal is high.
2. Experimental setup
The experimental setup is shown in Fig. 1. Here a random light pattern is generated from a simulated superposition of plane waves using random numbers, which is then sent to an SLM to produce a synthesized speckle field. The SLM has 512 × 512 pixels in the window of size 3.584×3.584mm. We pass a collimated laser of wavelength λ = 632.8nm through a polarizing beam splitter and a half-wave plate, before illuminating the SLM window. The speckle field is generated by modulation of the SLM and the returning light field is then magnified by a simple telescope system consisting of 150mm and 450mm biconvex lenses. The object is located at the focus plane of the 450mm lens, which is also the image plane of the SLM window. A 50 : 50 beam splitter is placed before the object in order to split the speckle field into two beams; the object beam (I(xS)) and the reference beam (I(xR)). The object beam illuminates the object and is then collected by a bucket detector, thus providing an computational GI setup. The additional reference beam for monitoring the light differentiates our system from previous experimental computational GI configurations. Since we are generating a computer hologram that is then sent to the SLM to create the speckle field, we can therefore predict the light field at the reference arm, negating the demand for a CCD camera, and requiring only a second bucket detector. It should be noted that for TGI based on our computational GI setup, only the object bucket detector is needed. The additional bucket detector in the reference arm is only required for NGI and DGI. Light intensities detected by the object and reference bucket detectors are indicated by S and R respectively, and the speckle field is described by I(x, y). As we use a 50 : 50 beam splitter, it is understood that I(x, y) = 2I(xS, yS) = 2I(xR, yR).
3. Iterative ghost imaging algorithms
In all iterative GI techniques, the transmitting object located after the beam splitter, BS2, is reconstructed by correlating the speckle field intensity measured at S and R, then adding together each successive frame with a suitable weighting factor. The transmitted light power detected after the object can be expressed as
3.1. Traditional Ghost Imaging
In TGI, the reconstruction result of the object, O(x, y) is retrieved from the correlation between S and I(x, y). We define for each iteration, i, the contribution to the reconstruction to be 
3.2. Differential ghost imaging
Differential GI , first performed by Ferri et al, utilizes a second bucket detector to extract a reference signal which is used in the reconstruction to weight the speckle field based on the average transmission signal relative to the average reference signal. Similarly, each contribution to the reconstruction can be expressed asEq. (3) and Eq. (4) are both identical however the first term in brackets of Eq. (4) is now weighted according to the average value of S, which is normalized to the average value of R. As demonstrated in  the DGI algorithm improves by order of magnitude the SNR of the measurement with respect to TGI. Moreover, a key difference from TGI, it is no longer sensitive to other sources of noise. For example, fluctuations in the laser power or changes to the SLM efficiency will affect both the reference signal and the transmitted signal, and thus the contribution to the reconstruction will be weighted more appropriately.
3.3. Normalized ghost imaging
3.3.1. Normalized ghost imaging with two detectors
As seen in Eq. (4), larger values of S measured by the bucket detector results in a greater weight for that particular speckle field, therefore external noise sources can still affect the overall reconstruction. There exists another iterative algorithm which instead normalizes each individual measurement S, as well as the running average, according to the reference signal R, resulting in an arguably more intuitive approach for dealing with time varying noise sources. We call this approach normailized GI (NGI). The algorithm used to describe each contribution to the reconstruction in NGI is given byEqs. (4) and (5) we can summarize the difference between the two algorithms as
3.3.2. Normalized ghost imaging with a single detector
In a computational GI setup, we can show that the additional detector used to measure the reference signal in DGI and NGI can instead be estimated based on the known light field reflected from the SLM and the average measured signal S for an arbitrary number of previous iterations. Calculating R negates the requirement for an additional detector, whilst improving the performance of the reconstruction compared to TGI, thus single-detector NGI (SNGI) is identical to the TGI experimental setup, with only a modified algorithm.
3.4. Signal-to-noise ratio analysis
To make a quantitative comparison between the NGI and the existing algorithms, we adopt a similar approach as used by Ferri et al and investigate the theoretical contribution to the signal-to-noise ratio (SNR) for objects with varying transmission functions. In  the authors express the average quantity of Eq. (4) in terms of the object transmission fluctuation δT(x, y) = T(x, y) − T̄,Eq. (7) is obtained under the assumptions of uniform illumination (the average speckle beams are constant over their area) and perfect resolution (the speckle area is much smaller compared to features of the object). The corresponding signal of DGI can be defined as Eq. (6), we can express the signal of NGI as
The speckle patterns used in our experiment exhibit complex-Gaussian behaviour, such that the intensity is exponentially distributed (see Fig. 2), and the noise associated to the measurement of O(x, y) can be expressed asEq. (10) may be omitted. Again, under the assumptions of uniform illumination and perfect resolution, the noise of DGI can be expressed as
Finally, we show that the SNR contribution for NGI is9]. For comparison the SNR contribution for TGI was shown to be
4. Experiment results
We generated a series of random speckle patterns using an SLM by simulating the interference of many plane waves on a computer. The real and imaginary amplitude components and the wave vector k⃗ of each simulated plane wave is Gaussian distributed. Figure 2 shows a typical example of the speckle patterns generated on the SLM and the exponentially distributed intensity for many patterns, implying that the speckle hologram has complex-Gaussian statistics, thereby a good approximation for real speckle fields . A binary transmissive object, 5mm × 5mm in size, is located after a 3× magnification telescope in the image plane of the SLM. Since we know both the object and the random speckle field projected to the SLM, we are able to simulate the expected results for comparison with our experiment. Experimental and simulated reconstruction results after 10000 iterations are shown in Fig. 3. The simulated reconstruction is produced assuming no external noise sources. The partially transmissive object used is indicated in the bottom right of Fig. 3. It is clear that the DGI and NGI algorithms provide very similar results, as predicted from the theory, and both show improved background subtraction compared to TGI.
Compared with the traditional computational GI setup, the NGI algorithm requires a reference bucket detector. However, as discussed in section 3.3.2, the advantage of computational GI means that we can replace this bucket detector with a virtual reference detector generating a simulated R. Thus we can negate the requirement for the reference detector and return the system to a true single element camera, which we call single-detector NGI (SNGI). The two major factors that dominate the value of R are from the different speckle patterns displayed on the SLM and fluctuations of the incident laser power. We can computationally predict changes to the value of R due to the speckle pattern, whereas fluctuations of the laser power can be simulated by using a rolling average for a particular series of S measurements. The bottom row in Fig. 3 shows the experimental results for reconstructing the object using the SNGI algorithm. We observe similar results compared with DGI and NGI algorithms indicating an improved performance compared with the TGI algorithm for single element camera.
To demonstrate the effect of object transmission function on the performance of NGI compared with TGI and DGI algorithms we used a similar experimental approach to that in Ref. . By scanning a knife edge (located in the image plane of the SLM, as before) across the speckle field in well defined steps (for which ΔT = 1), we measured the SNR’s for the final object reconstruction obtained after 5000 random speckle iterations. The beam size used was 10 × 10mm and the speckle size at the plane of the object was found to be δs ∼ 90μm, providing around Ns ∼ 12500 speckles. The experimental results and theoretical predictions for the SNR’s of each iterative algorithm are shown in Fig. 4. Note that the y-axis has been normalized to the number of iterations. We observe close quantitative agreement between the theory and the measurements. The results indicate that for low transmissive objects, all algorithms reconstruct with similar SNR, while for more transmissive objects the DGI and NGI algorithms become more efficient in comparison to TGI due to the differential nature of the reconstruction. Furthermore, we observe that when using a single detector, SNGI is a more efficient algorithm for reconstructing objects of all transmissions compared to TGI. We observe that for increasing transmissive objects SNGI becomes less efficient than NGI, for which the reason is the subject of ongoing research. Similar to , we find a systematic discrepancy between the experimental results of TGI and the theoretical predictions.
5. Normalization in matrix inverse algorithms
5.1. Introduction to matrix inverse algorithms and compressive sensing
As an alternative to the iterative techniques discussed above, we can choose to record all the signals for a complete set of speckle patterns and then treat the image reconstruction as one of matrix inversion. The series of M speckle patterns, each containing N pixels can be represented by a M × N matrix. If the object is also represented as an N element column vector, then the vector containing the measured signals is a M element vector. This relationship is expressed as
In the case where the number of speckle patterns equals the number of pixels then the M × N matrix is square, such that its inverse can be calculated and the object vector determined. However when M < N and or N is large, the system is ill-conditioned and calculating the inverse of the matrix is not straightforward. Problems of this type are wide spread in physics and techniques for solving them have been developed. Within our system the appeal is to reconstruct the image of N pixels from M measurements where M < N. That this is possible is based on the fact that natural images are sparse and the reconstruction can be obtained by solving a convex optimization problem , which is a generalization of a linear least squares problem. In contrast to iterative methods, compressive GI (CGI) needs to take all measurements, represented here, in some compressible basis (in this case a discrete cosine transform which has been applied to each row of the M × N matrix). Solving the convex optimization problem requires minimizing the ℓ1 norm .
5.2. Normalized compressive ghost imaging
By normalizing the measured object signal relative to the reference signal as performed above, such that S′ ≡ S/R, we can apply the CGI technique  to reconstruct our object. Equation (17) can then be written for normalized CGI (NCGI) as
Performing both NCGI and CGI analyses using the same experimental data (acquired using the experimental setup in Fig. 1) we obtain the reconstruction in Fig. 5. We observe a clear improvement using the NCGI algorithm compared to the CGI algorithm, manifest as an increased SNR value. The efficiency with which NCGI can reconstruct sparse images over CGI is determined by the level of noise in the system. We find that when there is no system noise present, both reconstructions are essentially identical. Thus the main improvement in employing NCGI over CGI with the additional reference detector is the ability to protect the reconstruction from time varying noise sources.
In conclusion we have compared different iterative GI methods to reconstruct an object and studied a new GI algorithm, which we call normalized GI (NGI). The performance of the differential GI (DGI) and NGI algorithms show good quantitative agreement as predicted by the theoretical foundations that support them. Our results indicate that by normalizing the measured signal relative to a reference signal, a more appropriate weighting factor is applied to the ensemble average of the estimated object, compared to the traditional GI (TGI) algorithm. Our analysis of the measured SNR and the object transmission shows a significant improvement for more transmissive objects in comparison to TGI. Furthermore, we have shown it is possible to apply normalization to systems with a single detector, SNGI, by estimating the reference signal. We have also investigated normalization within a compressive matrix inversion method, showing similar results to an non-normalized algorithm but with enhanced noise suppression. We believe the NGI algorithm will be a useful resource for imaging where alternative techniques are required in the future.
At the time when this work was ready for publication the authors found through private communication with Alessandra Gatti and Fabio Ferri that they had similar findings, for whom we thank for useful discussion. MJP would like to thank the Royal Society and the Wolfson Foundation. The work of JHS was supported by the DARPA Information in a Photon (InPho) Program. We gratefully acknowledge the financial support from the UK EPSRC.
References and links
3. A. Gatti, E. Brambilla, M. Bache, and L. A. Lugiato, “Correlated imaging, quantum and classical,” Phys. Rev. A 70, 013802 (2004). [CrossRef]
5. F. Ferri, D. Magatti, A. Gatti, M. Bache, E. Brambilla, and L. A. Lugiato, “High-resolution ghost image and ghost diffraction experiments with thermal light,” Phys. Rev. Lett. 94, 183602 (2005).
6. J. H. Shapiro, “Computational ghost imaging,” Phys. Rev. A 78, 061802 (2008). [CrossRef]
7. Y. Bromberg, O. Katz, and Y. Silberberg, “Ghost imaging with a single detector,” Phys. Rev. A 79, 053840 (2009). [CrossRef]
8. M. Duarte, M. Davenport, D. Takhar, J. Laska, T. Sun, K. Kelly, and R. Baraniuk, “Single-pixel imaging via compressive sampling,” IEEE Signal Processing Magazine 25, 83–91 (2008). [CrossRef]
10. J. Goodman, Statistical Optics (Wiley, 2000).
11. S. Boyd and L. Vandenberghe, Convex Optimization (Cambridge University Press, 2004).
12. D. L. Donoho and Y. Tsaig, “Fast solution of l1-norm minimization problems when the solution may be sparse,” IEEE Trans. Inf. Theory 54, 4789–4812 (2008). [CrossRef]
13. O. Katz, Y. Bromberg, and Y. Silberberg, “Compressive ghost imaging,” Appl. Phys. Lett. 95(13), 131110 (2009). [CrossRef]