We use coded aperture x-ray scatter imaging to interrogate scattering targets with a pencil beam. Observations from a single x-ray exposure of a flat-panel scintillation detector are used to simultaneously determine the along-beam positions and momentum transfer profiles of two crystalline powders (NaCl and Al). The system operates with a 3 cm range resolution and a momentum transfer resolution of 0.1 nm−1. These results demonstrate that a single snapshot can be used to estimate scattering properties along an x-ray beam, and serve as a foundation for volumetric imaging of scattering objects.
© 2012 OSA
We propose and demonstrate coded aperture x-ray scatter imaging (CAXSI). We use coded apertures to encode multiplex sampling to improve the efficiency of scatter tomography. Coded apertures have long been used to improve sensor throughput . Major coded aperture applications include uniformly redundant arrays for lensless imaging  and Hadamard transform spectroscopy . Coded apertures may also be viewed as “light field” encoders that enable radiance measurement using irradiance detectors . Building on studies of “reference structures” for compressive tomographic imaging [5–7], our group developed coded aperture snapshot spectral imaging (CASSI) [8, 9]. In 2009, we proposed an extension of this approach to compressive x-ray tomography . In each of these examples, coded apertures are used to alleviate space-time-spectral trade-offs and enable snapshot acquisition of data conventionally recorded sequentially. As with light field imagers, the coded aperture may be regarded as a reference structure that encodes radiance to remove range/angle ambiguity for scatter detected on irradiance sensors.
X-ray scatter imaging has shown promise for a wide variety of applications, including detection of abnormal structures in biological tissue [11–13], measurements of surface structure , and detection of explosives and other controlled substances [15–18]. Reference  gives an overview of x-ray scatter imaging for explosives detection and shows reconstructions of buried landmines using Compton backscatter imaging, as well as reconstructions of various plastics (nylon, PMMA, PE, PTFE, PVC) using coherent scatter computed tomography (CSCT). CSCT  has been applied to bone mineral density measurements  and detection of urinary stones . In addition, reference  demonstrates a fan beam energy-dispersive CSCT system which can detect various plastics in an aluminum case.
CSCT uses a series of images recorded at multiple angles to estimate an object’s coherent scatter properties. Another approach to scatter tomography is energy-dispersive x-ray diffraction tomography (EXDT) , which scans an object voxel-by-voxel using collimators and provides an effectively isomorphic mapping between the object voxels and the measurements. EXDT was demonstrated with an x-ray tube  and with a synchrotron source . It has been used to probe polymer and bone surfaces , to reduce the false alarm rate of luggage scanners in airline security , and to probe mineral content in thick cement samples .
The goal of this paper is to demonstrate that the momentum transfer in a scattering object may be estimated at each point along a 1D “pencil” beam by acquiring a single irradiance image. Our experimental system is depicted in Fig. 1. We use a 2D irradiance detector array perpendicular to the beam and a coded aperture between the object and detector to modulate the scattered radiation. Use of the coded aperture allows the irradiance detector to act more like a radiance detector by introducing an angular sensitivity to each pixel. To obtain a volumetric scatter image, the pencil beam could be scanned over a 3D object with estimation performed for each transverse position. Other multiplexing strategies, such as using fan beam or cone beam geometries, are expected to reduce the total x-ray exposure needed for tomographic reconstructions.
In the following section we describe a measurement model for the pencil beam coded aperture system. We present our experimental setup in Sec. 3, reconstruction techniques in Sec. 4, and discuss our experimental results in Sec. 5.
2. Pencil beam system model
In the pencil beam system depicted in Fig. 1 an x-ray source is filtered by a pinhole to produce a thin beam propagating along the z axis. The scattering object is placed between the primary and secondary apertures so that it is penetrated by the beam. The primary beam is stopped by the secondary mask to prevent it from flooding the detector image. Scattered x-rays diverge from the main beam to strike the aperture, where they are either absorbed or transmitted to the detector plane. Each pixel in the detector array receives scattered power from multiple points along the beam, and the structure of this multiplexing is controlled by the aperture code.
Upon scattering, an x-ray changes its wavevector by q = k′ − k, where k is the incident wavevector and k′ is the scattered wavevector. We consider elastic scattering, for which |k| = |k′|. This condition holds for Rayleigh scattering (including coherent scattering) and approximates Compton scattering at low energies or at low scattering angles. Elastic scattering follows Bragg’s law q = 2ν sin(θ/2), where q = |q|, ν is the x-ray frequency, and θ is the scattering angle as shown in Fig. 1.
Let the “scattering density” f(z, q) be the ratio of scattered power to incident power at position z along the beam and with momentum transfer q. This model only depends on the magnitude of the momentum transfer and not its direction, so it applies to liquids, fine powders (as in this experiment), and amorphous compounds. In the absence of a coded aperture, scattering at angle θ from position z produces an irradiance at radius ρ on the detector proportional to 1/(z2 +ρ2). The coded aperture is modeled by the transmission function t(ρ, ϕ) ∈ [0, 1] in the plane z = d, where (ρ, ϕ) are the polar radius and angle relative to the beam. The total irradiance at the detector point (ρ, ϕ) is
In contrast with previous studies of coded aperture imaging that emphasize codes that are orthongonal under translation , coded aperture range imaging requires codes that are orthogonal under magnification. Equation (1) is a scale transformation of the aperture code t(ρ, ϕ) as a function of object position z. The projected image of the aperture code is magnified by z/(z − d). One can disambiguate images corresponding to different values of z by applying aperture codes which are orthogonal under changes in scale (i. e. magnification). Harmonic codes (e.g. t(ρ) = cos(ρ)) have this property. To also disambiguate q the code must also vary as a function of ϕ. As a simple easily manufactured example, we choose the square grid
Equation (1) is the continuous measurement model for this system. This model is discretized by expanding the scattering density over compact voxel functions in the coordinates z and q. For this purpose we use the function rect(x) which is equal to unity for |x| < 1/2 and zero everywhere else. The voxels are chosen with sampling rates Δz in z and Δq in q and have centers (zj, qj). The discrete model for the scattering density isEqs. (3) and (4) are used in Eq. (1), the discrete model is expressed as the linear system Eq. (5)) can be used with numerical methods to estimate the object vector f, given measurements g and functional models for P(ν) and t(ρ, ϕ). We have computed the forward matrix H for our experimental system and used an algorithm derived from Ref.  to reconstruct the underlying object vector f for each object configuration (the procedure is described in detail in Sec. 4). The reconstructed objects are presented in Sec. 5, but first we describe the pencil beam CAXSI experiment in the following section.
3. Experimental methods
In order to build the pencil beam CAXSI experiment shown in Fig. 1, a standard diagnostic x-ray system, which has been described in detail previously [25, 26], was modified to include an optical bench, a collimator, a coded aperture, and a sample stage at adjustable positions in the beam. The x-ray source used was a General Electric (GE) model MX100 that has a tungsten target with a 12° anode angle. A 0.6 mm focal spot size was chosen to minimize focal spot blurring. The source produced Brehmsstrahlung radiation in the energy range 20–116 keV and also characteristic lines from Tungsten’s Kα transition doublet at 58.0–59.3 keV and from the Kβ transitions at 66.7–67.7 keV. The x-ray beam had an inherent filtration equivalent to 1.1 mm thick aluminum at 80 kVp. The normal acquisition mode was set at 116 kVp, 500 mAs .
Since broadband illumination is expected to degrade specificity of material classification , spectral shaping is critical to CAXSI. Toward that end, the beam was shaped by a 0.1 mm thick tungsten filter which served as a band-pass between approximately 30 keV and tungsten’s K-edge at 69.5 keV. The expected source spectrum was modeled using the semi-empirical x-ray spectrum modeling program XSPECT . XSPECT produced a model for the mean spectral number density N (ν) of photons illuminating the object. This model is plotted with a normalized maximum value in Fig. 2, and was used to calculate the power spectral density P(ν).
Scattered x-rays were collected with a stationary amorphous silicon indirect cesium iodide (CsI) flat panel detector (Paxscan, 4030 CB series, Varian Medical Systems) designed to perform with extended dynamic range. The detector had a pixel size of 194 μm and a matrix size of 2048 × 1536. The source-to-image distance was 201 cm. The acquisition setup was finely calibrated to maximize the signal-to-noise ratio of the acquired images. Specifically, the detector was gain-calibrated at the expected photon flux. Furthermore, it was offset calibrated before each acquisition with 16 dark frames to correct for structured noise. Post-calibrated images were acquired using the image acquisition and processing software ViVA (Varian Medical Systems).
A pencil beam was achieved using a primary aperture at a distance of about 130 cm from the source. The primary aperture consisted of a hole 2 mm in diameter drilled into a 6 mm thick lead sheet. Taking into account the focal spot size of 0.6 mm, the beam divergence half-angle is estimated to have been about 0.06°, which approximately satisfies the parallel ray condition for an ideal pencil beam. The secondary aperture was placed 180 cm from the source and consisted of another 6 mm thick lead sheet oriented parallel to the detector. The aperture, with its x-ray projection shown in Fig. 3, had a square center piece was designed to block the primary beam and the horizontal bar structures served as supports. The aperture code was designed to implement Eq. (2) with u = 9.9 cm−1 so that the period was 0.64 cm.
Two crystalline powders, sodium chloride (NaCl) and aluminum (Al), were chosen as scattering targets for their strong coherent-scatter cross sections and applicability to powder detection. These samples were placed in separate Nalgene vials of 1 cm diameter that were determined to have negligible contributions to the scatter images. The vials were placed so that the beam penetrated each sample over the full diameter of its vial. Exposures were taken with each of the two samples alone in the beam with their centers 60.2 cm from the detector, and a third exposure was taken with both NaCl and Al placed at 60.2 cm and 52.6 cm from the detector, respectively. This last exposure was meant to test the system’s ability to image multiple samples at different ranges in one snapshot. The detector images for each of these cases are shown normalized in Figs. 4(a), 4(b), and 4(c), where one can clearly see Debye rings modulated by the coded aperture.
In addition to the diffraction images, a background frame was acquired with no samples in the beam in order to measure any additional radiation from secondary sources in the system. The diffraction images along with the background images were used to reconstruct the scattering profiles of each test object as a function of position and momentum transfer, and these were compared with individual reference profiles measured by a Panalytical X’Pert PRO x-ray diffractometer with known sample positions. In the next section we describe the algorithm used to reconstruct the test objects, and in Sec. 5 we discuss the reconstruction results.
4. Object reconstruction
The data collected by the system discussed above consisted of superpositions of the scattered radiation from different test objects in the beam. In this section, we describe how the spatial and momentum transfer distribution of each object was estimated from the measurements.
Given the discrete measurement model in Eq. (5), each diffraction image is represented by a vector g. Our images also contained a noisy background image with mean μb so that the expected value at each pixel is given by g = Hf+μb. Treating the x-ray detection as a statistical process, our actual measurements y are approximated by the Poisson process
We propose to first estimate μb from b using a Poisson image denoising algorithm, and use the resulting estimate μ̂b of μb to reconstruct f. In particular, we estimate μb using a maximum penalized likelihood estimation method discussed in , according to which:24]. A pseudocode of this iterative reconstruction method is provided below:
- Initialize f(0) = HTy.
- The final estimate is given by f = f(i−1).
In our experimental setup, the detector array consisted of 2048 × 1536 square pixels. In order to reduce computational complexity, we perform “polar downsampling” of the measured diffraction images. Aside from modulation by the aperture, the diffraction patterns consist of concentric rings which can be effectively represented over bins in the polar coordinates (ρ, ϕ) relative to the beam position. In practice, few polar bins, relative to the size of the detector, are sufficient to reliably capture the information content in the diffraction images. In our experiments, the images were partitioned into 233 uniform radius bins between ρ = 2.5 cm and ρ = 11.5 cm. The polar angle was similarly segmented over its entire range into 120 bins. As a result of this strategy, we were able to reduce sampling from 2048 × 1536 to 233 × 120, which afforded significant savings in the computation of H.
We applied the GML algorithm described above to our diffraction data (Fig. 4) in order to estimate f for each test object, and we compare these reconstructions to reference data in the next section.
5. Results and discussion
The forward matrix H for the pencil beam system was calculated by sampling the object space with voxels of width 0.33 cm in z and 0.027 nm−1 in q. Here and in the following we replace numerical values of q with q/4π to agree with recent literature. From the single-frame diffraction images shown in Figs. 4(a) and 4(b), each test object’s coefficient vector f representing the scattering density f (z,q) was estimated using the methods described in the last section. Figures 5(a) and 5(b) show the estimated spatial distribution f (z) = ∫ dq f (z,q) of the scattering density for NaCl and Al placed separately in the beam at z = −60.2 cm. Although the beam penetrated only 1 cm of each sample, the spatial reconstructions both have FWHM equal to 3 cm. The peak of the spatial profile occurs at z0 = −59.3 cm for both samples, corresponding to a position error of 1.5%. This demonstrates the along-beam ranging capability of the coded aperture system.
Figures 5(c) and 5(d) show the estimated momentum transfer profiles f (z0, q) at the spatial peak z0 = −59.3 cm for NaCl and Al, respectively. The dominant peaks in both profiles were accurately reconstructed, and there is evidence of some of the smaller peaks. The dominant peak for NaCl was reconstructed at q = 1.767 nm−1 (0.2% error) with a FWHM of 0.06 nm−1 (3.6%). The dominant peak for Al was reconstructed at q = 2.149 nm−1 (0.4% error) and a FWHM of 0.07 nm−1 (3.3%). These results show that the pencil beam coded aperture system can be used to estimate the scattering structure of target samples without a priori position information.
To test the system’s ability to distinguish different objects in the beam within a single snapshot, one vial of NaCl and one vial of Al were placed along the beam at z = −60.2 cm and z = −52.6 cm, respectively. In this configuration the pencil beam passed first through the NaCl sample and then through the Al sample, producing the diffraction image shown in Fig. 4(c). As before, the coefficients f were reconstructed and produced the spatial scattering profile f (z) shown in Fig. 6(a). The spatial distribution shows a peak for NaCl at z = −59.3 cm (1.5% error) and another for Al at z = −52.0 cm (1.1% error). Momentum transfer profiles for these two locations are shown in Figs. 6(b) and 6(c). The dominant peak for NaCl was reconstructed at q = 1.79 nm−1 with FWHM equal to 4.7%. The dominant peak for Al was estimated to lie at q = 2.137 nm−1 (0.3% error) with FWHM equal to 4.7%. The reconstructed momentum transfer profiles are consistent, whether the objects are measured separately or placed together in the beam.
Finally, the combined NaCl and Al diffraction pattern (Fig. 4(c)) is plotted in Fig. 7(a) over the polar coordinates (ρ, ϕ). Figure 7(b) shows the modeled diffraction pattern Hf̂ based on the corresponding object estimate f̂. The RMS error between these diffraction patterns is approximately 10% of the peak signal value in Hf̂, indicating agreement between the two images and providing combined validity to the measurement model, the experimental process, and the reconstruction algorithm. We suspect that shot noise and dark current noise are the primary contributors to the differences between the two images.
There are several key limitations to the system demonstrated here. Since we have used an irradiance detector, the spectral characteristics of the source play a role in achieving good spatial and momentum transfer resolution. We expect that a superposition of several sharp spectral features will improve the accuracy and resolution of the reconstruction. Also more careful design of the aperture code, including finer features, will likely yield improved spatial and momentum transfer resolution.
We expect that the accuracy of the estimated momentum transfer should be useful in performing material classification at points along the beam. Our results demonstrate the use of coded apertures and an irradiance detector to simultaneously range and estimate the momentum transfer from two coherent-scatter samples in a single snapshot. Future work will include construction of fan beam and cone beam coded aperture systems and investigation of novel spatial and spectral coding techniques for the source and scattered radiation, including energy-resolved detection. Looking forward, we regard the pencil beam system as a solid foundation for understanding more sophisticated CAXSI experiments.
This work was supported by the Department of Homeland Security, Science and Technology Directorate Explosives Division through contract HSHQDC-11-C-0083.
References and links
1. D. J. Brady, Optical Imaging and Spectroscopy (Wiley-OSA, 2009). [CrossRef]
3. M. Harwit and N. J. A. Sloane, Hadamard Transform Optics (Academic Press, 1979).
4. A. Veeraraghavan, R. Raskar, A. Agrawal, A. Mohan, and J. Tumblin, “Dappled photography: Mask enhanced cameras for heterodyned light fields and coded aperture refocusing,” ACM Transactions on Graphics 26, 69-1–69-12 (2007).
5. D. J. Brady, N. P. Pitsianis, and X. Sun, “Reference structure tomography,” J. Opt. Soc. Am. A 21, 1140–1147, (2004). [CrossRef]
10. K. Choi and D. J. Brady, “Coded aperture computed tomography,” in “Adaptive Coded Aperture Imaging, Non-Imaging, and Unconventional Imaging Sensor Systems ,” SPIE 7468, 74680B-1–74680B-10, (2009).
12. J.-P. Schlomka, A. Harding, U. van Stevendaal, M. Grass, and G. L. Harding, “Coherent scatter computed tomography: a novel medical imaging technique,” SPIE 5030, 256–265, (2003). [CrossRef]
13. M. T. M. Davidson, D. L. Batchelar, S. Velupillai, J. D. Denstedt, and I. A. Cunningham, “Laboratory coherent-scatter analysis of intact urinary stones with crystalline composition: a tomographic approach,” Phys. Med. Biol. 50, 3907 (2005). [CrossRef] [PubMed]
14. R. J. Cernik, K. H. Khor, and C. Hansson, “X-ray colour imaging,” Journal of the Royal Society Interface 5, 477–481 (2008). [CrossRef]
15. G. Harding and B. Schreiber, “Coherent x-ray scatter imaging and its applications in biomedical science and industry,” Radiat. Phys. Chem. 56, 229–245, (1999). [CrossRef]
16. G. Harding, “X-ray scatter tomography for explosives detection,” Radiat. Phys. Chem. 71, 869–881 (2004). [CrossRef]
17. R. W. Madden, J. Mahdavieh, R. C. Smith, and R. Subramanian, “An explosives detection system for airline security using coherent x-ray scattering technology,” SPIE 7079, 707915-1–707915-11, (2008).
18. C. Crespy, P. Duvauchelle, V. Kaftandjian, F. Soulez, and P. Ponard, “Energy dispersive x-ray diffraction to identify explosive substances: Spectra analysis procedure optimization,” Nucl. Instrum. Methods Phys. Res. A 623, 1050 – 1060, (2010). [CrossRef]
20. J. Delfs and J.-P. Schlomka, “Energy-dispersive coherent scatter computed tomography,” Appl. Phys. Lett. 88, 243506 (2006). [CrossRef]
21. G. Harding, M. Newton, and J. Kosanetzky, “Energy-dispersive x-ray diffraction tomography,” Phys. Med. Biol. 35, 33 (1990). [CrossRef]
22. C. Hall, P. Barnes, J. Cockcroft, S. Colston, D. Husermann, S. Jacques, A. Jupe, and M. Kunz, “Synchrotron energy-dispersive x-ray diffraction tomography,” Nucl. Instrum. Methods Phys. Res. B 140, 253 – 257 (1998). [CrossRef]
24. W. H. Richardson, “Bayesian-based iterative method of image restoration,” J. Opt. Soc. Am. 62, 55–59, (1972). [CrossRef]
25. A. Chawla and E. Samei, “Geometrical repeatability and motion blur analysis of a new multi-projection x-ray imaging system,” IEEE Nuclear Science Symposium Conference Record 5, 3170 –3173, (2006). [CrossRef]
26. A. Chawla, S. Boyce, L. Washington, H. McAdams, and E. Samei, “Design and development of a new multi-projection x-ray system for chest imaging,” IEEE Trans. Nucl. Sci. 56, 36–45, (2009). [CrossRef]
28. S. R. Beath and I. A. Cunningham, “Pseudomonoenergetic x-ray diffraction measurements using balanced filters for coherent-scatter computed tomography,” Med. Phys. 36, 1839–1847, (2009). [CrossRef] [PubMed]
29. C. Dodge and M. Flynn, “Advanced integral method for the simulation of diagnostic x-ray spectra,” Med. Phys. 33, 1983 (2006).
30. E. Kolaczyk and R. Nowak, “Multiscale likelihood analysis and complexity penalized estimation,” The Annals of Statistics 32, 500–527, (2004). [CrossRef]