# X-Ray Phase-Contrast Imaging with Three 2D Gratings

Authors: Ming Jiang; Christopher Lee Wyatt; Ge Wang
Journal: International Journal of Biomedical Imaging
 (2008)

---

## Abstract

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.


---

## Body

## 1. Introduction

Since its introduction in 1973, X-ray computed tomography (CT) has revolutionize 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 and100μ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 distanceL 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 px and py 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 ZT 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 line 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
(1)p0,x=Ldpx,p0,y=Ldpy. 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 px and py in
the x- and y-directions, respectively. The complex
transmittance function of G1 is given by the following Fourier series:
(2)T1(x,y)=∑m,nam,nexp(i2πmxpx)exp(i2πnypy). 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],
(3)E(x,y,z)=exp(ikz)∑m,nβm,n(z)exp(i2πmxpx+i2πnypy),
where
(4)βm,n(z)=am,nexp(−iπλm2zpx2−iπλn2zpy2). The intensity of the wave field is given by
(5)I(x,y,z)=|E(x,y,z)|2=∑m,nρm,n(z)exp(i2πmxpx+i2πnypy), where (6)ρm,n(z)=∑m',n'βm+m',n+n'(z)βm',n'(z)¯ Let the periods px and py be of a rational ratio in the following
sense:
(7)px=Mp,py=Np, for some positive numbers M,N,
and p,
where M and N are integers. Let
(8)ZT=2M2N2p2λ.
Then we have
(9)βm,n(z)=am,nexp(−2πiN2m2zZT−2πiM2n2zZT).
It follows that
(10)βm,n(z+sZT)=βm,n(z),
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=sZT,
which is referred to as the Talbot images or self-images in [25, 28]. At other distances such as
(11)Du,v=vuZT,u=1,2,…,v=1,2,…,
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 Du,v.
It depends on the fraction v/u and the shape of the grating groove [28].Some examples are as follows. Foru=2 and v=1,(12)E(x,y,z|z=Du,v)=exp(ikz)T1(x+MNp2,y+MNp2);
for u=4 and v=1,(13)E(x,y,z|z=Du,v)=exp(ikz)2i[T1(x,y)−T1(x+MNp2,y+MNp2)],
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 (14)ES(x,y,z)≈exp(ikz)∑m,nβm,n(z)exp(2πimxpx+2πinypy+iΦ(x−mλzpx,y−nλzpyy)). Because
(15)Φ(x−mλzpx,y−nλzpy)≈Φ(x,y)−mλzpx∂Φ∂x(x,y)−nλzpy∂Φ∂x(x,y), we obtain
(16)ES(x,y,z)≈exp(ikz+iΦ(x,y))∑m,nβm,n(z)×exp(2πimpxϕx(x,y)+2πinpyϕy(x,y)), where
(17)ϕx(x,y)=x−iλz2π∂Φ∂x(x,y),ϕy(x,y)=y−iλz2π∂Φ∂y(x,y). 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
(18)IS(x,y,z)=∑m,nρm,n(z)exp(2πimpxϕx(x,y)+2πinpyϕy(x,y)). Let the complex transmittance function of G2 be
(19)T2(x,y)=∑m,nbm,nexp(2πimxpx)exp(2πinypy),
with the same periods as those of the phase grating G1, a moiré pattern
after the analyzer grating G2 is given by [22, 23]
(20)Mkx,ky(x,y,z)=IS(x,y,z)×T2(x+χx,y+χy)=∑m,nbm,nρm,n(z)×exp(2πimpx(ϕx(x,y)+χx))×exp(2πinpy(ϕy(x,y)+χy)),
where χx=kxpx/Kx,kx=0,1,…,Kx−1  and χy=kypy/Ky,ky=0,1,…,Ky−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,
(21)ηm,n(x,y)=bm,nρm,n(z)exp(2πimpxϕx(x,y)+2πinpyϕy(x,y))=1KxKy[∑kx=0Kx−1∑ky=0Ky−1Mkx,ky(x,y,z)exp(2πikxmKx+2πikynKy)] if the grating coefficients am,n and bm,n vanish rapidly when |m| and |n| are large compared to Kx and Ky.
The coefficients ρm,n(z) are of real values at the Talbot distance ZT, hence (22)2πmpxϕx(x,y)+2πnpyϕy(x,y)=arg[ηm,n(x,y)]. 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 Kx and Ky [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
(23)sincmπ2sincnπ2. 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 px=2μm and py=3μm.
We set M=2 and N=3 according to (7). The Talbot distance
in (8) is ZT=232cm. The result is shown in Figure 2.Figure 1
Schematics of our 2D grating-based phase-contrast
image system. See the text for detailed
configuration. Note that the scales of the grating features are exaggerated for illustration.Intensity distribution for an incident coherent plane wave of energy of 40 keV behind the grating G1 of periodspx=2μm and py=3μm.
The grating transmission function is assumed to be a Ronchi grating. The Talbot
distance is ZT=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.
(a)(b)The second experiment is to simulate the
phase stepping measurement in Section2.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 px=4μm and py=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 ZT=1,652cm.
When the grating periodicities are smaller, the Talbot distances will be
smaller, too. The phase object is given by the following function:
(24)Φ(x,y)=Φx(x)Φy(y), where
(25)Φx(x)=∑i=1SxAx,iψax,i,bx,i(x),Φy(y)=∑j=1SyAy,jψay,j,by,j(y).Sx and Sy are the numbers of basis functions. Ax,i and Ay,j are the weights for the basis functions. ax,i and ay,j are the scale factors of the basis functions. bx,i and by,j are the centers of the basis functions. The
basis function ψa,b(t) is as follows. Let β(t) be a nonnegative function with support in [−1,1] and let
(26)α(t)=∫−∞tβ(s)ds. For a>0 and b∈R,
let
(27)φa,b(t)=β(x−ba),ψa,b(t)=aα(x−ba).
It follows that
(28)∂Φ∂x(x,y)=Φy(y)∑i=1SxAx,iφax,i,bx,i(x),∂Φ∂y(x,y)=Φx(x)∑j=1SyAy,jφay,j,by,j(y). Each φa,b is of the support [b−a,b+a]⋅ax,i     and bx,i are chosen such that the supports of the φax,i,bx,i are disjoint. Ax,i<10−3 are chosen such that the ϕx and ϕy 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 Ax,i is chosen to simulate weak phase objects. The ay,j,by,j,
and Ay,j are chosen in the same manner. We have used in
our experiments the following basis functions:
(29)β(t)={1−|x|,if|x|<1,0,if|x|<1,β(t)={1,if|x|<1,0,if|x|<1. 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 [−px/2,px/2]×[−py/2,py/2].In Figure3, 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 Sx=5 and Sy=5.
The derivatives are presented in Figure 4.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.
(a)(b)(c)Distributions of phase gradients.The phase object is specified in (24) with Sx=5 and Sy=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.
(a)(b)(c)(d)In Figure4, 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 χx=17 and χy=17.
The original distributions of gradients Φx and Φy 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 distanceZT.
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 δ:
(30)Φ(x,y,θ)=2πλ∫δ(X,Y,z)dz,
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 thex- 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.


## 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 distanceL 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 px and py 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 ZT 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 line 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
(1)p0,x=Ldpx,p0,y=Ldpy. 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 px and py in
the x- and y-directions, respectively. The complex
transmittance function of G1 is given by the following Fourier series:
(2)T1(x,y)=∑m,nam,nexp(i2πmxpx)exp(i2πnypy). 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],
(3)E(x,y,z)=exp(ikz)∑m,nβm,n(z)exp(i2πmxpx+i2πnypy),
where
(4)βm,n(z)=am,nexp(−iπλm2zpx2−iπλn2zpy2). The intensity of the wave field is given by
(5)I(x,y,z)=|E(x,y,z)|2=∑m,nρm,n(z)exp(i2πmxpx+i2πnypy), where (6)ρm,n(z)=∑m',n'βm+m',n+n'(z)βm',n'(z)¯ Let the periods px and py be of a rational ratio in the following
sense:
(7)px=Mp,py=Np, for some positive numbers M,N,
and p,
where M and N are integers. Let
(8)ZT=2M2N2p2λ.
Then we have
(9)βm,n(z)=am,nexp(−2πiN2m2zZT−2πiM2n2zZT).
It follows that
(10)βm,n(z+sZT)=βm,n(z),
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=sZT,
which is referred to as the Talbot images or self-images in [25, 28]. At other distances such as
(11)Du,v=vuZT,u=1,2,…,v=1,2,…,
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 Du,v.
It depends on the fraction v/u and the shape of the grating groove [28].Some examples are as follows. Foru=2 and v=1,(12)E(x,y,z|z=Du,v)=exp(ikz)T1(x+MNp2,y+MNp2);
for u=4 and v=1,(13)E(x,y,z|z=Du,v)=exp(ikz)2i[T1(x,y)−T1(x+MNp2,y+MNp2)],
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 (14)ES(x,y,z)≈exp(ikz)∑m,nβm,n(z)exp(2πimxpx+2πinypy+iΦ(x−mλzpx,y−nλzpyy)). Because
(15)Φ(x−mλzpx,y−nλzpy)≈Φ(x,y)−mλzpx∂Φ∂x(x,y)−nλzpy∂Φ∂x(x,y), we obtain
(16)ES(x,y,z)≈exp(ikz+iΦ(x,y))∑m,nβm,n(z)×exp(2πimpxϕx(x,y)+2πinpyϕy(x,y)), where
(17)ϕx(x,y)=x−iλz2π∂Φ∂x(x,y),ϕy(x,y)=y−iλz2π∂Φ∂y(x,y). 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
(18)IS(x,y,z)=∑m,nρm,n(z)exp(2πimpxϕx(x,y)+2πinpyϕy(x,y)). Let the complex transmittance function of G2 be
(19)T2(x,y)=∑m,nbm,nexp(2πimxpx)exp(2πinypy),
with the same periods as those of the phase grating G1, a moiré pattern
after the analyzer grating G2 is given by [22, 23]
(20)Mkx,ky(x,y,z)=IS(x,y,z)×T2(x+χx,y+χy)=∑m,nbm,nρm,n(z)×exp(2πimpx(ϕx(x,y)+χx))×exp(2πinpy(ϕy(x,y)+χy)),
where χx=kxpx/Kx,kx=0,1,…,Kx−1  and χy=kypy/Ky,ky=0,1,…,Ky−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,
(21)ηm,n(x,y)=bm,nρm,n(z)exp(2πimpxϕx(x,y)+2πinpyϕy(x,y))=1KxKy[∑kx=0Kx−1∑ky=0Ky−1Mkx,ky(x,y,z)exp(2πikxmKx+2πikynKy)] if the grating coefficients am,n and bm,n vanish rapidly when |m| and |n| are large compared to Kx and Ky.
The coefficients ρm,n(z) are of real values at the Talbot distance ZT, hence (22)2πmpxϕx(x,y)+2πnpyϕy(x,y)=arg[ηm,n(x,y)]. 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 Kx and Ky [8, 21–23].


## 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 px and py in
the x- and y-directions, respectively. The complex
transmittance function of G1 is given by the following Fourier series:
(2)T1(x,y)=∑m,nam,nexp(i2πmxpx)exp(i2πnypy). 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],
(3)E(x,y,z)=exp(ikz)∑m,nβm,n(z)exp(i2πmxpx+i2πnypy),
where
(4)βm,n(z)=am,nexp(−iπλm2zpx2−iπλn2zpy2). The intensity of the wave field is given by
(5)I(x,y,z)=|E(x,y,z)|2=∑m,nρm,n(z)exp(i2πmxpx+i2πnypy), where (6)ρm,n(z)=∑m',n'βm+m',n+n'(z)βm',n'(z)¯ Let the periods px and py be of a rational ratio in the following
sense:
(7)px=Mp,py=Np, for some positive numbers M,N,
and p,
where M and N are integers. Let
(8)ZT=2M2N2p2λ.
Then we have
(9)βm,n(z)=am,nexp(−2πiN2m2zZT−2πiM2n2zZT).
It follows that
(10)βm,n(z+sZT)=βm,n(z),
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=sZT,
which is referred to as the Talbot images or self-images in [25, 28]. At other distances such as
(11)Du,v=vuZT,u=1,2,…,v=1,2,…,
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 Du,v.
It depends on the fraction v/u and the shape of the grating groove [28].Some examples are as follows. Foru=2 and v=1,(12)E(x,y,z|z=Du,v)=exp(ikz)T1(x+MNp2,y+MNp2);
for u=4 and v=1,(13)E(x,y,z|z=Du,v)=exp(ikz)2i[T1(x,y)−T1(x+MNp2,y+MNp2)],
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 (14)ES(x,y,z)≈exp(ikz)∑m,nβm,n(z)exp(2πimxpx+2πinypy+iΦ(x−mλzpx,y−nλzpyy)). Because
(15)Φ(x−mλzpx,y−nλzpy)≈Φ(x,y)−mλzpx∂Φ∂x(x,y)−nλzpy∂Φ∂x(x,y), we obtain
(16)ES(x,y,z)≈exp(ikz+iΦ(x,y))∑m,nβm,n(z)×exp(2πimpxϕx(x,y)+2πinpyϕy(x,y)), where
(17)ϕx(x,y)=x−iλz2π∂Φ∂x(x,y),ϕy(x,y)=y−iλz2π∂Φ∂y(x,y). 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
(18)IS(x,y,z)=∑m,nρm,n(z)exp(2πimpxϕx(x,y)+2πinpyϕy(x,y)). Let the complex transmittance function of G2 be
(19)T2(x,y)=∑m,nbm,nexp(2πimxpx)exp(2πinypy),
with the same periods as those of the phase grating G1, a moiré pattern
after the analyzer grating G2 is given by [22, 23]
(20)Mkx,ky(x,y,z)=IS(x,y,z)×T2(x+χx,y+χy)=∑m,nbm,nρm,n(z)×exp(2πimpx(ϕx(x,y)+χx))×exp(2πinpy(ϕy(x,y)+χy)),
where χx=kxpx/Kx,kx=0,1,…,Kx−1  and χy=kypy/Ky,ky=0,1,…,Ky−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,
(21)ηm,n(x,y)=bm,nρm,n(z)exp(2πimpxϕx(x,y)+2πinpyϕy(x,y))=1KxKy[∑kx=0Kx−1∑ky=0Ky−1Mkx,ky(x,y,z)exp(2πikxmKx+2πikynKy)] if the grating coefficients am,n and bm,n vanish rapidly when |m| and |n| are large compared to Kx and Ky.
The coefficients ρm,n(z) are of real values at the Talbot distance ZT, hence (22)2πmpxϕx(x,y)+2πnpyϕy(x,y)=arg[ηm,n(x,y)]. 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 Kx and Ky [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
(23)sincmπ2sincnπ2. 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 px=2μm and py=3μm.
We set M=2 and N=3 according to (7). The Talbot distance
in (8) is ZT=232cm. The result is shown in Figure 2.Figure 1
Schematics of our 2D grating-based phase-contrast
image system. See the text for detailed
configuration. Note that the scales of the grating features are exaggerated for illustration.Intensity distribution for an incident coherent plane wave of energy of 40 keV behind the grating G1 of periodspx=2μm and py=3μm.
The grating transmission function is assumed to be a Ronchi grating. The Talbot
distance is ZT=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.
(a)(b)The second experiment is to simulate the
phase stepping measurement in Section2.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 px=4μm and py=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 ZT=1,652cm.
When the grating periodicities are smaller, the Talbot distances will be
smaller, too. The phase object is given by the following function:
(24)Φ(x,y)=Φx(x)Φy(y), where
(25)Φx(x)=∑i=1SxAx,iψax,i,bx,i(x),Φy(y)=∑j=1SyAy,jψay,j,by,j(y).Sx and Sy are the numbers of basis functions. Ax,i and Ay,j are the weights for the basis functions. ax,i and ay,j are the scale factors of the basis functions. bx,i and by,j are the centers of the basis functions. The
basis function ψa,b(t) is as follows. Let β(t) be a nonnegative function with support in [−1,1] and let
(26)α(t)=∫−∞tβ(s)ds. For a>0 and b∈R,
let
(27)φa,b(t)=β(x−ba),ψa,b(t)=aα(x−ba).
It follows that
(28)∂Φ∂x(x,y)=Φy(y)∑i=1SxAx,iφax,i,bx,i(x),∂Φ∂y(x,y)=Φx(x)∑j=1SyAy,jφay,j,by,j(y). Each φa,b is of the support [b−a,b+a]⋅ax,i     and bx,i are chosen such that the supports of the φax,i,bx,i are disjoint. Ax,i<10−3 are chosen such that the ϕx and ϕy 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 Ax,i is chosen to simulate weak phase objects. The ay,j,by,j,
and Ay,j are chosen in the same manner. We have used in
our experiments the following basis functions:
(29)β(t)={1−|x|,if|x|<1,0,if|x|<1,β(t)={1,if|x|<1,0,if|x|<1. 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 [−px/2,px/2]×[−py/2,py/2].In Figure3, 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 Sx=5 and Sy=5.
The derivatives are presented in Figure 4.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.
(a)(b)(c)Distributions of phase gradients.The phase object is specified in (24) with Sx=5 and Sy=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.
(a)(b)(c)(d)In Figure4, 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 χx=17 and χy=17.
The original distributions of gradients Φx and Φy 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 distanceZT.
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 δ:
(30)Φ(x,y,θ)=2πλ∫δ(X,Y,z)dz,
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 thex- 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.



---

*Source: 827152-2008-03-24.xml*