Ghost Tomography

Ghost tomography using single-pixel detection extends the emerging field of ghost imaging to three dimensions, with the use of penetrating radiation. In this work, a series of spatially random x-ray intensity patterns is used to illuminate a specimen in various tomographic angular orientations with only the total transmitted intensity being recorded by a single-pixel camera (or bucket detector). The set of zero-dimensional intensity readings, combined with knowledge of the corresponding two-dimensional illuminating patterns and specimen orientations, is sufficient for three-dimensional reconstruction of the specimen. The experimental demonstration of ghost tomography is presented here using synchrotron hard x-rays. This result expands the scope of ghost imaging to encompass volumetric imaging (i.e., tomography), of optically opaque objects using penetrating radiation. For hard x-rays, ghost tomography has the potential to decouple image quality from dose rate as well as image resolution from detector performance.


INTRODUCTION
Ghost imaging (GI) first emerged in the domain of visible-light optics [1]. The term arose from Einstein's description of quantum entanglement as "spooky action at a distance," since initial realizations of the method utilized pairs of entangled photons. Classical implementations of GI have since been developed using pairs of correlated, coherent wave fields [2]. Very recently, GI has been achieved with atoms [3], electrons [4], and x-rays [5][6][7][8][9][10]. However, to date, none of the reported studies utilizing penetrating radiation have attempted to map the interior structure of a genuinely three-dimensional (3D) sample. GI clearly has the potential to achieve such tomographic reconstruction, constituting a natural extension of previously reported lower-dimensional ghost images. Here we report on the realization of ghost tomography (GT) using hard x-rays, whose penetrating power for optically opaque objects significantly extends both the applicability and utility of the technique.
Synthesizing images via the superposition of linearly independent intensity maps, random or otherwise, is the essence of GI [11][12][13][14]. These maps, when random, may be generated through quantum processes such as shot noise or through classical means such as spatially random masks. Nonrandom intensity maps may be generated using suitable deterministic masks. We restrict consideration to random illuminating intensity maps for the remainder of this paper on account of their ease of construction for x-ray fields. A key feature of GI is that the ensemble of superposed linearly independent illuminating intensity maps is formed by photons (or other imaging quanta) that never pass through the sample. A weak copy of the illuminating field, which may be obtained, e.g., using a beam splitter, does pass through the object but only the total number of transmitted quanta is measured by a single-pixel detector in a so-called "bucket signal." 3D GT of optically opaque objects has not been demonstrated in the literature; however, similar concepts have been presented. Direct (as opposed to computed) GT of optically transparent objects has been reported using optical coherence imaging [15]. Ghost topography (or 3D surface imaging by GI) has been developed in a remote sensing context using time-of-flight with a single-pixel camera [16]. 2D GT with terahertz radiation has been presented by Mohr et al. [17].
Since no imaging quanta that pass through the sample are ever registered by a position-sensitive detector, GI resolution is independent of the bucket detector. This is an important distinction between GT and computed tomography (CT): in CT, 3D volume resolution is limited by the pixel size of the detector. In CT, the pixel size of the detector and geometric magnification of the imaging system immediately suggest an appropriate discretization (i.e., voxel size) for the 3D volume; in GT, the correct 3D discretization must be found based on analysis of the ensemble of illuminating fields.
In two-dimensional (2D) GI applications, the parallelized intensity-intensity cross-correlation between the bucket and any one pixel of the random reference maps is used to compute the ghost image [11,12]. In what follows, we show that simply combining this method with tomography can be insufficient for 3D imaging, and we present new reconstruction schemes that give superior results. Our experimental proof-of-concept for GT significantly expands the scope of GI, empowering it to realize genuinely 3D reconstructions.

METHOD
A schematic of our experimental setup for x-ray GT is shown in Fig. 1. Illumination of a spatially random 1-mm-thick Ni foam with normally incident 26 keV hard x-rays from a synchrotron created spatially random intensity illumination patterns, such as that shown in Fig. 2A. An ensemble of such speckle patterns was obtained via transverse displacement of the foam over a 2D square grid with a step size of 400 μm. This transverse step size was chosen to be considerably larger than the width of the illuminating speckle intensity autocovariance (Fig. 2C), to reduce the degree of correlation between illuminating speckle images. The sample for our proof-of-concept x-ray GT experiment was an Al cylinder with diameter 5.60 mm, into which was drilled two cylindrical holes with respective diameters of 1.98 and 1.50 mm (Fig. 2B). This sample was secured to a rotation stage and illuminated with attenuated copies of the spatially random intensity maps, obtained by using a 220 Laue x-ray reflection from a (001) Si wafer beam splitter.
Approximately 2000 random-illumination intensity maps were used, forming a linearly independent mathematical basis [13] for the 2D ghost projection images. Noise-free simulations, assuming a similar experimental setup with 64 2 pixel illumination patterns, were conducted in Kingston et al. [18]. These showed image quality (using conventional GI) degraded as the number of illumination patterns was reduced, and that 1000 patterns per view-angle approached the lower limit of object resolvability. Regularization techniques, such as compressed sensing, were shown to produce significant image quality improvements; however, given that the data measured here would contain noise and other artifacts, we opted for 2000 measurements per view-angle.
GI spatial resolution [19] cannot be determined based on pixel size, so we have therefore used Fourier ring correlation (FRC)applied to the ensemble of illuminating speckle fields-to estimate the resolution of our imaging system as approximately 100 μm (Fig. 2E). FRC yields a best-case limit estimate for 2D spatial resolution. This is quite distinct from the point spread function (PSF) of the 2D imaging system (Fig. 2C) that is calculated as Fig. 1. Experimental setup for x-ray GT. Synchrotron x-rays from an undulator are passed through a spatially random mask (not shown). The resulting random 2D speckled beam is split into two copies by a crystal beam splitter working in a Laue diffraction condition. The diffracted beam, much weaker in intensity than the direct beam, is passed through the sample before being registered at the position-insensitive bucket detector. The direct beam, consisting of photons that never pass through the object, is measured over the position-sensitive detector. An ensemble of spatially random illuminating patterns is created by transversely displacing the mask. Note that only the spatially integrated signal (termed the "bucket signal") for each bucket-beam measurement is utilized in the x-ray GT. The process is repeated for a variety of angular orientations θ of the sample. FRC results from registered image subsets, used to determine 3D GI resolution (determined as the reciprocal distance at which correlation drops below 1 bit). The relevant (spck/bckt) resolution result of 100 μm, was used to select the 3D discretization for the tomographic reconstructions in Fig. 3. Image pairs include: spck/spck-speckle images compared at θ 0°and θ 68.750°; spck/blur-speckle image at θ 0°c ompared to blurred speckle image at θ 68.750°; bckt/bckt-bucket images compared at θ 0°and θ 68.750°; spck/bckt-speckle image at θ 0°compared to bucket image at θ 0°.

Research Article
the normalized autocovariance of the set of illuminating spatially random intensity fields [9,19]; this PSF estimates the resolution of conventional GI by cross-correlation [11,12]. In either case, the theoretically achievable resolution is limited below by twice the pixel size of the detector used to measure the illumination patterns. See Section 2A in Supplement 1 for further detail. These estimates of resolution allow us to select an appropriate discretization for our 3D reconstructed image.
For tomographic imaging, we repeated the 2000 randomillumination intensity maps for each of N 14 projection angles. Due to some instability of the beam-line vacuum, the x-ray beam dropped out at random time intervals during data acquisition. Unlike conventional CT, GI is insensitive to such random signal dropouts because it utilizes intensity-intensity correlations. Further, the object rotation angles were chosen using a quasi-random (or low-discrepancy) additive recurrence sequence of angles, θ, with step size equal to is the golden ratio. This equates to Δθ 111.25°and can be achieved equivalently with an angular step size for the object rotations of 180°− 111.25° 68.75°. Quasi-random sequences appear to be random locally but are highly ordered globally. Hence, at any time the experiment is ceased, the angle set acquired will be an approximately uniform sampling of 0, π radians. Each detected image was registered via an indirect detector, consisting of a scintillator screen, lens system and a 2560 × 2160 pixel pco.edge 5.5 (PCO AG, Germany) scientific CMOS-based camera with pixel pitch of 6.5 μm, and binned down to the resolution that was determined via FRC. Each object-free 2D reference-illumination beam was paired with a bucket-beam image containing the object (e.g., blue box in Fig. 2D corresponding to yellow box in Fig. 2A). The total signal in the blue-box region was summed to give the bucket signal B j,θ corresponding to the jth realization of the spatially random illuminating pattern, at sample rotation angle θ. The spatially random intensity pattern I j x, y illuminating this same region corresponds to the spatially resolved intensity map of the beam that did not pass through the object, where x, y are Cartesian coordinates in the detector plane.

ANALYSIS AND RESULTS
In 2D GI, the cross-correlation GI formula [11,12] may be used to estimate the 2D intensity transmission function T x, y; θ for a given fixed object rotation θ as the ensemble average (intensityintensity correlation): Here, B av,θ is the average bucket signal for a given θ, and M θ is the number of bucket measurements taken at each orientation. This ensemble average constitutes the superposition of linearly independent spatially random intensity maps I j x, y referred to above. Subsequent reconstruction of the 3D attenuation function using conventional tomography algorithms showed that the cross-correlation GI formula is inadequate for 3D imaging (see Figs. S8a and S8d of Supplement 1 and accompanying text in Section S4B). A posteriori information of the sample must be leveraged to produce a meaningful GT reconstruction. To achieve this, we employed iterative cross-correlation via the Landweber algorithm coupled with smoothness priors [18]. The relaxation parameter used was where J θ is the number of measurement pairs at angle θ, and σ 2 is the variance of the spatially random speckle patterns. Such a 2D reconstruction was performed for each of the 14 pseudo-random projection angles θ. Applying conventional tomographic reconstruction techniques to the resulting projection images produced a reasonable but very noisy tomogram (see Figs. S8b-c and S8e-f in Supplement 1 and accompanying text in Section S4B).
In the above two-step reconstruction scheme (ghost reconstruction followed by tomography), each projection image is reconstructed separately from the others. This is not the optimal approach, as projections at different angles are obviously related. A better result can be achieved by direct reconstruction, where one recovers the 3D volume directly from the bucket signals, thus using all measured information simultaneously; the intermediate step of recovering the 2D x-ray ghost projection images can be removed. A gradient descent (or the Landweber) algorithm for direct iterative tomographic reconstruction from bucket signals has very recently been developed in Section V of the simulation-based study of Kingston et al. [18]. A smoothness prior and enforced positivity in attenuation coefficient were included here to improve the result. Vertical and horizontal slices through the resulting x-ray GT reconstructions are shown in Figs. 3A and 3B, respectively. These may be compared to the conventional CT reconstructions obtained in the same set of experiments, as given in Figs. 3C and 3D. A semitransparent rendering of the 3D reconstructed ghost tomogram is given in Fig. 3E, with horizontal and vertical cutaway 3D renderings in Figs. 3F and 3G, respectively. The nontrivial preprocessing steps required to achieve the Research Article results in Fig. 3 are detailed in Supplement 1. This supplement also gives further GI tomographic slices in Fig. S9. The reconstructed sample densities, as obtained from the x-ray ghost tomograms, are quantitative. Using the XCOM (National Institute of Standards and Technology, NIST) database [20], the attenuation per unit density of Al at 26 keV is 1.65 cm 2 ∕g. The density of Al is 2.70 g∕cm 3 , giving an expected linear attenuation coefficient of 4.455 cm −1 . From the reconstructed x-ray ghost tomogram with 52 μm pixel dimension (i.e., binned ×8) the mean attenuation of the Al is measured as 4.80 cm −1 and corresponds to the attenuation of Al at 25.3 keV. This increase in attenuation is most likely due to inclusions of higher-Z metals to form the Al alloy, together with the difference in spectrum between the direct and diffracted beams.

DISCUSSION
The ability to achieve quantitative 3D imaging, in a GI geometry where none of the photons passing through the object are ever detected with a position-sensitive camera, is remarkable. A key observation is the previously mentioned impracticality of two-step GT achieved by simply combining 2D GI at each projection, with standard tomographic reconstruction concepts. Rather, we emphasize that direct GT was seen to be much more effective as iterative refinement occurs in a whole-of-data-set manner. We thereby reconstructed a 140 × 140 × 72 voxel ghost tomogram using approximately 26,000 bucket measurements spread over 14 sample-rotation orientations, equating to over 50 reconstructed voxels per bucket measurement. This efficiency was enabled by harnessing a posteriori assumptions (or enforcing priors). GT is particularly suited to such efficiencies, the further exploitation of which may aid in a long-term aim of reduced dose relative to conventional imaging.
The x-ray dose in conventional CT is the product of flux-perview-angle, and the number of viewing angles used. Dose reduction in conventional CT is necessarily achieved by some combination of: (i) reducing flux per viewing angle, decreasing tomogram signal-to-noise ratio; or (ii) reducing the number of viewing angles, typically resulting in artifacts in the tomogram since the reduced sampling conditions do not satisfy the incoherence requirements of compressed sensing.
An effective combination of compressed sensing and CT, known as compressive tomography (see e.g., [21][22][23]), uses illumination masks to obstruct patterns of pixels from the illuminating radiation. This method has the potential to reduce dose without introducing artifacts into the tomogram, as the illumination masks can be constructed to satisfy the incoherence requirements of compressed sensing. Compressive tomography can be achieved as a special case of GT, by using a pinhole mask translated to each unblocked pixel element per view-angle; thus, GT can be at least as effective at dose reduction, with the potential to be more so.
From a broader perspective, our demonstration of GT shows how the GI approach is naturally able to relax the constraints placed on image quality by dose rate, as well as on image resolution by detector performance. This is a fundamental departure from conventional imaging paradigms. GT affords the flexibility of independently varying a number of parameters, such as the illumination masks, exposure time, number of bucket readings, and number of object orientation angles. As a consequence, the resolution level can be optimized against dose rate in a manner that takes into account prior knowledge about the sample. For instance, illumination masks can be designed in a way to minimize dose to the sample (according to prior knowledge of it) while maintaining high resolution. This is not possible in direct imaging using a pixel array detector, which requires all pixels to be illuminated regardless of the object being imaged. Therefore, while it is important to compare the performance of ghost and conventional imaging-as done in this paper-it is crucial to recognize that GI is not just a different way of making images. Also, a GI or tomography system can be designed to be adaptive, in the sense that it can be optimized for the features of the object being imaged [24]. This may have great practical advantages when using ionizing radiation. For instance, it is not too far-fetched to imagine how GT with mask engineering could be used in future radiological practice. By using the available prior information, the dose could be spatially and angularly distributed to statistically match the object of interest (for instance the brain or the lungs) given that the size, shape, and density of these organs or body parts are well known.
We close this discussion with some speculations regarding possible future challenges and limitations of GT, the overcoming of which will progress the maturation of the technique. Limitations are not strictly technological per se; many can be overcome by further technique development. Significant work in GI remains to realize a genuine competitive advantage with regards to dose and photon budgets, relative to more conventional forms of imaging, or the discovery of at least one niche area in which the analyses GI provides are clearly superior. We expect a reasonably straightforward translation of GT to low-flux cone-beam laboratory x-ray sources through the use of computational GI, e.g., as performed in Zhang et al. [8]. A challenge with a true single-pixel bucket detector setup (as opposed to the accumulated bucket signal used here) will be that of alignment. One possibility is an alignment phantom and protocol that uses only signal intensity, i.e., maximum intensity equates to an aligned system; a second option would be an auxiliary alignment system. When considering a computational x-ray GI setup, the question arises: What magnifications could be reasonably achieved for tomographic ghost microscopy? In this case the mask used to structure the illumination can be preimaged with high resolution, e.g., with a transmission electron microscope. Assuming the mask has features up to the resolution of this prerecorded image, the limit to GI and GT resolution would be that of the translation accuracy of the mask, as well as the alignment with the assumed position of the single-pixel detector. We anticipate that hybrid systems, combining GI with conventional imaging approaches using a 2D position sensitive bucket detector, will be of future interest when considering optimization of dose and resolution. In addition to compressed sensing (e.g., [12,18]) and regularization in general (as exemplified here), we expect that artificial intelligence and deep learning are likely to play an important role in the future evolution of whole-of-data set reconstruction approaches for GT (as demonstrated by Shimobaba et al. [25] for GI).

CONCLUSION
We report the experimental demonstration of GT, obtained using hard x-rays. We demonstrate that GT is able to computationally measure the 3D internal distribution of a sample by a set of bucket readings of the total transmitted x-ray intensity from the sample. The task is accomplished by illuminating the sample with a known, varying set of 2D x-ray patterns for each rotation angle of the sample. We discussed our strategies for data Research Article acquisition and processing, showing how direct tomographic reconstruction from bucket readings is much more effective than the two-step approach of tomographic inversion following GI reconstruction of individual projections. These results outline how the flexibility of engineering a GT measurement marks a radical departure from the conventional tomographic imaging paradigm, being able to make optimal use of the available information to maximize tomogram quality and minimize the radiation dose used.