X-ray phase-contrast imaging with three 2D gratings.

X-ray imaging is of paramount importance for clinical and preclinical imaging but it is fundamentally restricted by the attenuation-based contrast mechanism, which has remained essentially the same since Roentgen's discovery a century ago. Recently, based on the Talbot effect, groundbreaking work was reported using 1D gratings for X-ray phase-contrast imaging with a hospital-grade X-ray tube instead of a synchrotron or microfocused source. In this paper, we report an extension using 2D gratings that reduces the imaging time and increases the accuracy and robustness of phase retrieval compared to current grating-based phase-contrast techniques. Feasibility is demonstrated via numerical simulation.

1. INTRODUCTION

Since its introduction in 1973, X-ray computed tomography (CT) has revolutionized medical imaging and become a cornerstone of modern radiology. Improved resolution and reduced dose are two critical requirements in biomedical applications and remain the focus of CT research and development. The rapid development of small animal models of various human diseases has also generated the need for preclinical imaging. Spurred by the refinements of CCD cameras and microfocus X-ray tubes in the 1990s, a number of micro-CT systems were constructed, reaching image resolutions between 10 and 100μm. These scanners, while producing high-resolution images, do not provide high-contrast or low-dose imaging of in vivo processes in small animal models. Most importantly, many normal and diseased tissues such as cancers display poor image contrast on CT/micro-CT images.

Current X-ray CT techniques generate contrast from differences in attenuation. Hence, weakly absorbing tissues display poor contrast [1]. The X-ray absorption coefficient is roughly proportional to the fourth power of the atomic number Z, apart from the jumps at absorption edges [2]. Therefore, soft tissues are difficult to image in terms of only the magnitude of the X-ray beam. However, X-ray phase-contrast imaging uses the diffraction properties of X-rays and reveals significant differences indistinguishable in an attenuation coefficient distribution. The wave propagation of X-rays is characterized macroscopically by the complex refraction index of the underlying medium. The refraction index n of soft biological tissues is approximately 1, namely, n = 1 − δ + iβ where β and δ quantify the intensity attenuation and phase shift, respectively. At X-ray energies of 15–25 keV, δ is about three orders of magnitude greater than β, depending on the X-ray wavelength and electron density. Thus, phase-contrast imaging enjoys a sensitivity to soft-tissue 1000 times greater than that with attenuation-based imaging [3]. Moreover, because the cross-section of the X-ray phase shift varies slower than the attenuation counterpart at higher energies, X-ray phase-contrast imaging can be performed at higher energies, further reducing radiation dose. More importantly, the phase-contrast mechanism is potentially a new design criterion for contrast agents, including those that can be functionalized for molecular imaging.

The use of phase shift as an X-ray contrast mechanism has generated considerable interest over the past decade [1]. X-ray phase-contrast imaging can be implemented using interferometry or diffractometry approaches which measure the first order derivative of the phase shift, or inline holography approaches which measure the X-ray intensity. Because traditional interferometry and diffractometry require monochromatic X-rays and high-precision crystals, they are restricted to synchrotron radiation. The inline method is feasible with a polychromatic X-ray tube, and thus more practical.

Gabor originally proposed the inline holography approach in 1948, for which he was awarded the Nobel Prize. In 2000, Takeda et al. demonstrated how phase-contrast imaging captures structures in carcinoma with 30 μm resolution, and observed that the X-ray energy for phase-contrast imaging can be increased from 17.7 keV to 30 keV. Wilkins et al. developed inline holography techniques with a polychromatic microfocus X-ray source [4]. Wu et al. generalized the theory and improved the techniques along this direction [5, 6]. In these kinds of systems, the phase retrieval is done by either free-space propagation or grating-based shearing interferometer. The advantage of using the interferometer is both sensitivity in the phase shift measurement and insensitivity to mechanical drift or vibration. However, such systems are limited by the spatial coherence requirement. While synchrotrons are expensive, microfocus (< 100 μm) X-ray tubes require a low-duty cycle due to limited spot size and resultant anode heating.

Recently, grating-based phase-contrast X-ray interferometry has been developed with very promising results [7-11]. The working principle is nearly the same as that of Talbot interferometry [12-14] in the visible light regime. Such a system uses a conventional (noncoherent) X-ray tube and gratings to overcome the well-known problem of spatial and temporal incoherence [11, 15]. Specifically, their approach uses three 1D gratings to produce coherent wavelets from a hospital X-ray tube, construct interference patterns at an appropriate distance, and detect phase changes from wave-fronts distorted by an object. Thus, quantitative phase-contrast imaging becomes possible in the case of centimeter-sized objects. This grating-based interferometric approach has advantages over other X-ray phase-contrast imaging techniques because it allows the cone-beam imaging geometry, large field of view, and efficient use of polychromatic low-brilliance sources [16]. Since the total phase shift is an integral over the X-ray path length that is nearly straight, it generates projections of the refractive index that can be inverted using X-ray CT techniques [8]. For example, X-ray phase-contrast tomography was first performed using crystal X-ray interferometry in combination with phase-shifting interferometry [1, 17].

The developments outlined above represent both an outstanding opportunity and major challenges in the field for practical widespread use. In this paper, we develop the theory, methods, and top-level system design so that a conventional hospital-grade X-ray tube may be used in a practical way for phase-contrast imaging of small animals. Our general hypothesis is that 2D grating-based phase-contrast imaging techniques can be developed to produce more informative projective, and tomographic images of biomedical interest than current techniques.

The rest of the paper is organized as follows. In the second section, we will propose a 2D grating-based phase-contrast imaging system and formulate the forward imaging model. In the third section, we will describe our numerical simulations and report representative results. Finally, we conclude by discussing relevant issues.

2. SYSTEM AND METHOD

2.1.Top-level design of a 2D grating-based system

Our proposed 2D grating-based system architecture requires three gratings: the source grating G0, the phase grating G1, and analyzer grating G2, which are all aligned along the principal optical axis. The source grating G0 is placed right behind an X-ray source, acts as a set of apertures to form a 2D array of virtual sources. The object under imaging is placed in front of the phase grating G1. The phase grating G1 is placed at a distance L from the source grating G0. The phase grating G1 induces a periodic spatial modulation in either amplitude or phase of the resultant wave front after X-ray penetrates the phase object. The phase and analyzer gratings G1 and G2 are of the same periodicity p and p in the x - and y- directions, respectively, and must be separated at a Talbot distance d The distance d can also be one of the fractional Talbot distances. However, a smaller d is preferable in practice. We choose d to be the first Talbot distance Z in (8) in this study, though there are alternative choices at a smaller fractional Talbot distances.

In this setup, each virtual source from the source grating G0 is incoherent with respect to the others, and small enough to create a fringe pattern with good visibility in the image plane after the analyzer grating G2. As the source mask G0 can contain a large number of individual apertures, each creating a sufficiently coherent virtual point source, standard X-ray generators with source sizes of more than a square millimetre can be used efficiently [11]. This key observation implies that hospital X-ray tubes can be used for phase-contrast imaging [11, 15, 16, 18]. To ensure that each virtual source contributes constructively to the fringe patterns at appropriate distances downstream, by a geometrical argument as in [11, 15, 16, 18], the periodicity of the source grating G0 is required to be The longitudinal or lateral periodicity of the fringes does not depend on the wavelength of the radiation [19]. Hence the setup is achromatic, allowing the use of a broadband source. These fringe patterns can be detected using a CCD camera. The phase-stepping approach [20-22] can be adapted to retrieve the gradients of the object phase shift, from which it is possible to reconstruct the phase profile of the object.

The gratings are usually made of gold and/or silicon. The thickness of gold pillars of several tens of microns is considered sufficient to block X-rays effectively. Roughly speaking, 10 μm is needed at 20 keV, 30 μm at 30 keV, and 60 μm at 40 keV [23]. Furthermore, the pitch of the grating must be about the spatial coherent length of primary X-rays, which is on the order of microns. Such an X-ray amplitude grating can be fabricated by lithography and in the case of the G0 and G2 gratings, gold electroplating, followed by coating in epoxy and bonding to a frame for mechanical stability [23].

2.2.Theory of 2D Talbot interferometry

In this section, we present the mathematical theory for the setup in the previous section.

2.2.1.Talbot effect under plane wave illumination

Talbot effect of 2D grating was first analyzed in [24]. The treatment was refined further in the following years. We consider the case of a 2D grating G1 illuminated coherently by a unit-amplitude plane-wave X-rays exp (ikz) of wavelength λ where k = 2π/λ is the wave number, though the case of Gaussian beams can also be investigated [25]. The optical axis is parallel to the z-axis. The 2D grating G1 on the xy-plane at z = 0 is of periods p and p in the x- and y, respectively. The complex transmittance function of G1 is given by the following Fourier series: We assume that gratings are infinitely extended in this preliminary study. Under a paraxial approximation, the wave field at a distance z behind the grating G1 is, by the Fresnel-Kirchhoff diffraction formula [25-27], where The intensity of the wave field is given by where Let the periods p and p be of a rational ratio in the following sense: for some positive numbers M, N, and p, where M and N are integers. Let Then we have It follows that for any integer s ∈ ℤ. Hence, the wave field is longitudinally periodic in the z-direction.

This is the Talbot effect discovered in [12]. The same complex amplitude of the wave field as that of the complex transmission function of the grating is generated at z = sZ, which is referred to as the Talbot images or self-images in [25, 28]. At other distances such as The wave field can be expressed as a superposition of laterally shifted replicas of the exit wave field right after the grating [26]. It is not possible to find reassembled images at every fractional Talbot distance D. It depends on the fraction v/u and the shape of the grating groove [28].

Some examples are as follows. For u = 2 and v = 1, for u = 4 and v = 1, which is quite different from the 1D grating case [26, 29].

2.2.2.Phase stepping technique

When a phase object that causes a phase shift Φ(x, y) is placed in front of the grating G1, under a stationary phase approximation, the wave field at a distance z behind the grating G1 is Because we obtain where Note that a higher-order expansion is necessary when the phase shift is not small [25, 30, 31]. The intensity of the wave field is then given by Let the complex transmittance function of G2 be with the same periods as those of the phase grating G1, a moiré pattern after the analyzer grating G2 is given by [22, 23] where χ = k p/K, k = 0, 1,…, K − 1 and χ = k p/k, k = 0, 1,…, k − 1 are the scanned displacements of G2 in the x - and y-directions, respectively. The phase gradients can be reconstructed by the inverse discrete Fourier transform (DFT) as follows. By the inverse DFT, if the grating coefficients a and b vanish rapidly when |m| and |n| are large compared to K and K. The coefficients ρ(z) are of real values at the Talbot distance Z hence The lower-order terms such as m = ± 1 and n = ± 1 are typically used because their contributions are dominant in the measured fringes. Higher-order harmonics in (18) may cause error in the determination of the phase gradients. The error can be reduced by sophisticated fringe analyze methods with sufficiently large K and K [8, 21–23].

3. NUMERICAL EXPERIMENTS

We report numerical simulation results in this section. The grating transmission functions are chosen to be of the Ronchi ruling with the Fourier coefficients in (2) and (19) given by The first experiment is to simulate the Talbot effect in Section 2.2.1. An incident coherent plane wave of energy of 40 keV is used. The grating periods are p = 2μm and p = 3 μm. We set M = 2 and N = 3 according to (7). The Talbot distance in (8) is Z = 232cm. The result is shown in Figure 2.

Figure 2

Intensity distribution for an incident coherent plane wave of energy of 40 keV behind the grating G1 of periods p = 2 μm and p = 3 μm The grating transmission function is assumed to be a Ronchi grating. The Talbot distance is Z = 232cm. (a) The distribution of the grating transmission function. (b) The distribution of the wave intensity at the Talbot distance. The spatial unit is μm.

The second experiment is to simulate the phase stepping measurement in Section 2.2.2 for phase objects. The X-ray energy is 40 keV and has a wavelength λ = 0.03nm. The parameters of the gratings G1 and G2 are p = 4μm and p = 4μm. The parameters are chosen according to those in [11]. The grating transmission function is assumed to be a Ronchi grating as in the previous experiment. The first Talbot distance is Z = 1, 652 cm. When the grating periodicities are smaller, the Talbot distances will be smaller, too. The phase object is given by the following function: where S and S are the numbers of basis functions. A and A are the weights for the basis functions. a and a are the scale factors of the basis functions. b and b are the centers of the basis functions. The basis function ψ(t) is as follows. Let β(t) be a nonnegative function with support in [− 1, 1] and let For a > 0 and b ∈ R, let It follows that Each φ is of the support [b − a, b+a] . a and b are chosen such that the supports of the φ are disjoint. A < 10−3 are chosen such that the ϕ and ϕ in (17) are within the principle value interval [−π, π] to avoid the phase wrapping issue in finding the phase angle by (22). This magnitude 10−3 of A is chosen to simulate weak phase objects. The a, b and A are chosen in the same manner. We have used in our experiments the following basis functions: We have run the simulation for both basis functions, different grating periodicities, and various stepping lengths. Consistent results are obtained. Representative results are shown in the following figures. Images are all sampled at the step 1000λ = 0.03 μm in [−p/2, p/2] × [−p/2, p/2].

In Figure 3, the wave intensity distributions right behind the grating G1 are z = 0 μm, and at the first Talbot distance are z = 1, 652cm by (14), and wave intensity at the first Talbot distance computed by the first-order approximation by (18) is presented. The result demonstrates that the first-order formula in (15) is a good approximation for weak phase objects in computing the downstream wave field. The parameters for the phase object in (24) are S = 5 and S = 5 The derivatives are presented in Figure 4.

Figure 3

Wave intensity distributions after a phase object. (a) Right behind the grating G1; (b) At the first Talbot distance; (c) At the first Talbot distance by the first-order approximation. The spatial unit is μm.

Figure 4

Distributions of phase gradients.The phase object is specified in (24) with S = 5 and S = 5. (a) The true x-gradient; (c) The measured x-gradient by phase stepping measurement. (b) and (d) are for the y-gradient. The spatial unit is μm.

In Figure 4, the phase stepping technique is used to measure the derivatives of the phase objects. The same parameters as in Figure 3 are used. The phase steps are χ = 17 and χ = 17. The original distributions of gradients Φ and Φ and those reconstructed gradient distributions from (22) for m = 1, n = 0 and m = 0, n = 1 are presented in Figure 4. The errors in phase measurement compared to their exact values are smaller than the floating error at the orders 10−17 and 10−9% for the maximal absolute error and the maximal relative error, respectively.

4. DISCUSSIONS AND CONCLUSION

Talbot interferometry techniques can be applied for X-ray phase-contrast imaging using a polychromatic hospital grade source coupled with an amplitude grating to form a set of equidistant source wavelets. As a result, we can avoid using a synchrotron-radiation source or a microfocus X-ray tube. This may eventually lead to various grating-based X-ray phase-contrast imaging systems. Further research topics include better grating design, more accurate phase retrieval, optimized prototyping, and identification of the optimal X-ray energy balancing image resolution, contrast and dose, as well as in vivo testing in different applications such as mouse studies. We have demonstrated this technique at the first-order Talbot distance Z. It is reported that the same technique works at the fractional Talbot distances in (11) with appropriate approximation for wave intensities under measurement.

Using the 2D grating-based phase-contrast imaging system, phase-contrast tomography is feasible. Tomography requires the phase deformation profile Φ(x, y, θ) which is the projection of the complex refractive index δ: where (X, Y, z) denotes a 3D position in an object to be reconstructed, and θ the angle between the xy and XY axes. To recover the projection in the form (30), we may calculate Φ(x, y, θ) by integrating its gradient along the x-direction or y-direction, or by a filtering technique as in [15, 18]. Then it is straightforward to reconstruct an underlying distribution of the refractive index δ using available CT methods. When the source emits a cone beam or beam shapes other than plane waves, G1 and G2 should have different periodicities and patterns. As there are theories on the Talbot effect available for spherical wave or Gaussian beams [25], the results in this manuscript can be extended to those cases. For example, in [32], the capillary plates with hexagonal patterns have been applied to phase-contrast imaging.

The results in the above are based on the paraxial approximation of the Fresnel integral for plane waves. The Talbot effect, based on a nonparaxial formulation, as well as spherical and Gaussian beams has been studied with similar results and can be applied for phase-contrast imaging with gratings [25]. The grating-based X-ray phase-contrast tomography can also be achieved using more complicated methods based on the Helmholtz equation, a fundamental governing equation [27]. In that case, the forward model would become more accurate at a much higher computational cost. Accordingly, we may have to use an iterative approach to reconstruct phase-contrast images from sufficiently many projections. This is a typical inverse coefficient problem of partial differential equations. The Born approximation can be used to find approximate solutions of the Helmholtz equation [33], with higher-order Born approximations used to improve the solution accuracy.

The main contribution of this paper is the use of 2D (chessboard) gratings so that a continuous operating hospital X-ray tube can be converted into a 2D array of coherent source wavelets, and truly 2D phase-contrast projections optimally formed. The use of 2D gratings instead of the state of the art 1D gratings has the major advantage of achieving nearly isotropic phase detection with fewer alignment problems. It is theoretically impossible to obtain both the x- and y-derivatives of a 2D phase object from 1D phase-stepping measurement along only one direction. We have simulated both 1D and 2D phase-stepping measurements and demonstrated this point. Two-axis scan with one 1D grating requires not only translations but also rotations, which means extra mechanical hardware. On the other hand, the 2D grating interferometry technique can reduce the mechanical rotations as compared to the 1D grating phase-stepping measurement. Hence, the 2D grating interferometry technique is more efficient than the 1D case. Up to date, no group has reported such a 2D-grating-based x-ray phase-contrast imaging system, which we are developing for biomedical applications. During the revision of this work, we noted the work [34], which was submitted on 2 August 2007, after our initial submission of this paper on 14 June 2007. The principal difference between our work and theirs is that our approach is interferometry-based while their “coded-aperture approach is not an interferometric one, and it is therefore substantially different from the grating interferometer method” [34].

In conclusion, we have described a 2D grating-based system and theory for superior phase retrieval and reported our numerical simulation results showing its feasibility. We are actively constructing the system for testing our algorithm with real data.