Full-wave approach for x-ray phase imaging

Yongjin Sung; Colin J. R. Sheppard; George Barbastathis; Masami Ando; Rajiv Gupta

doi:10.1364/OE.21.017547

1. Introduction

Since their original discovery in 1895, x-rays have been extensively used in medical imaging, inspection, and non-destructive evaluation [1, 2]. With the aid of tomographic methods and contrast agents, x-ray imaging has been extended to the diagnosis of a multitude of diseases in various organ systems. In conventional radiography, the measured intensity variation is primarily a record of the attenuation of an x-ray beam due to photoelectric absorption and incoherent scattering [3]. Tomographic imaging extends this by ascribing a linear attenuation coefficient to each unit voxel in the imaged object. In contrast, x-ray phase imaging (XPI) [4–12] records intensity variations that encode coherent scattering or refraction of light due to an object. The information acquired in XPI is attributed to the electron density of the object [13] and is complementary to the atomic number available in conventional radiography or tomography. In soft-tissues of different types, electron density differences predominate. XPI, therefore, can image soft tissues with much higher contrast resolution than traditional radiography [14–17].

Experimental techniques for XPI can be divided into three categories: analyzer-based imaging [4, 5], interferometry [6–8], and propagation-based methods [9–12]. Irrespective of the experimental setup, analysis of XPI systems and development of reconstruction algorithms require an accurate forward model that explains the image formation process. In XPI, the image is formed largely by three processes: (i) x-ray/object interaction, (ii) propagation of the beam in the free space, and (iii) interaction of the beam with other imaging elements such as an analyzer crystal or grating. In propagation-based XPI, which is of main interest in this paper, the third process plays no role and will not be considered. For x-ray/object interaction, a 2-D transmittance function typically represents the effect of an object on the incident x-ray beam profile. In this approach, the magnitude and argument of the complex transmittance function are obtained from integrating the attenuation coefficient and optical path length of volume elements (i.e., voxels), respectively, along the optical axis [18, 19]. The so-called projection approach is reasonable when the object is thin and the angular deviation of the rays can be ignored. In thick or large objects, ray tracing may be more appropriate [20]; however, the computational cost of ray tracing can be excessive for a source of finite aperture size. In both the projection and ray tracing approaches, the choice of an object plane, for which the transmittance function is calculated, is arbitrary and light diffraction within the object is ignored. For free-space propagation of an x-ray beam, the Fresnel kernel approach has been widely adopted under the assumption that the overall angular deviation of the x-ray beam is small and that there is no beam divergence such as that in fan or cone-beam geometry [18]. More rigorously, the Fresnel-Kirchhoff integral has been also adopted, which includes such a beam divergence effect [19, 20]. In this case, however, the effect of a finite source size is modeled as a convolution operator and the projection approximation usually needs to be adopted for the x-ray/object interaction.

In this study, we combine the x-ray-object interaction and free-space propagation into an integrated wave optics framework. For this, we spatially and angularly decompose x-ray beams emitted from a finite source aperture into multiple plane waves. For each plane wave input, we calculate the light field at the detector plane by solving the wave equation under the first Rytov approximation [21, 22]. The total light field is calculated as the coherent sum of the fields for all plane wave inputs, from which the intensity at the detector plane can be simply obtained. Our model does not rely on the paraxial approximation and is suitable for simulating x-ray phase imaging of a large object in a practical cone-beam setup. By using the Rytov approximation, we obtain results for thick specimens in more robust fashion than the models adopting the projection or first Born approximations. This should be emphasized in the x-ray imaging regime, where the refractive index is very close to one [23].

2. Full-wave approach for x-ray phase imaging

2.1 Scalar diffraction theory

Figure 1 shows a schematic diagram of the geometry considered in this study. An x-ray source, represented by a planar aperture with a known distribution of intensity, is located at z = -R₁. The object to be imaged is located at z = 0, and the detector that records the intensity of the total field (straight-through plus scattered from the object) is placed at z = R₂. Temporal coherence of an x-ray source is negligibly small. We model this by assuming a planar aperture filled with incoherent point sources. The complex refractive index of each voxel in an object encodes its interaction with an x-ray beam and it can be written in the following form [13]:

n = 1 - δ + i β .

Here, the parameters δ and β are responsible for the phase delay and absorption, respectively, of the x-ray beam caused by the object. In the x-ray cone-beam geometry of Fig. 1, the incident field on a specimen can be decomposed into multiple spherical waves, each emanating from a point source in the aperture. The incoherent superposition of these point sources is equivalent to all the x-rays emitted from the source aperture.

Fig. 1 Schematic diagram of the imaging geometry.

Download Full Size | PDF

From Weyl’s formula [24], the complex amplitude of the spherical wave emitted from a point $(x_{0}, y_{0})$ can be written on the detector as:

U^{(i)} (x, y; x_{0}, y_{0}) = \int \int_{- \infty}^{\infty} \exp [- i k (p_{1} x_{0} + q_{1} y_{0})] u^{(i)} (x, y; p_{1}, q_{1}) d p_{1} d q_{1} .

Here, the variable λ is the wavelength and

k = 2 π / λ

is the wave number of the x-ray beam in the background medium. The function

u^{(i)} (x, y; p_{1}, q_{1})

represents the complex amplitude of the plane wave propagating in the direction

(p_{1}, q_{1}, m_{1})

:

u^{(i)} (x, y; p_{1}, q_{1}) = (i / λ m_{1}) \exp {i k [p_{1} x + q_{1} y + m_{1} (R_{2} + R_{1})]},

where the dependent variable

m_{1}

is defined by

m_{1} = \sqrt{1 - p_{1}^{2} - q_{1}^{2}}

. Note that the term

\exp [- i k (p_{1} x_{0} + q_{1} y_{0})]

in Eq. (2) represents the phase shift due to the displacement of the point source from the center of the aperture. For an incident plane wave, the scattered field from the object can be written in a simple form with the first Born or Rytov approximation [25]. In the Born approximation, the strength of scattered field is assumed to be much smaller than that of incident field. In the Rytov approximation, the gradient of phase of scattered field is assumed to be very small. The Born approximation is valid when the optical thickness of the object is small, while the Rytov approximation places no restriction on the size of the object but only requires that the gradient of refractive index at the boundaries of the object be small [25]. In a recent study, we provided the accuracy of the Rytov approximation in x-ray phase imaging [23]. The scattered field for each individual component of spherical wave illumination in Eq. (2) may be obtained by summing up the scattered fields for all the constituent plane waves. The total field recorded on the detector can be written as

U (x, y; x_{0}, y_{0}) = \int \int_{- \infty}^{\infty} \exp [- i k (p_{1} x_{0} + q_{1} y_{0})] u (x, y; p_{1}, q_{1}) d p_{1} d q_{1},

where the function

u (x, y; p_{1}, q_{1})

is

u^{(i)} (x, y; p_{1}, q_{1}) [1 + {\bar{u}}^{(s)} (x, y; p_{1}, q_{1})]

according to the first Born approximation [24] and

u^{(i)} (x, y; p_{1}, q_{1}) \exp [{\bar{u}}^{(s)} (x, y; p_{1}, q_{1})]

according to the first Rytov approximation [26, 27]. The dimensionless function

{\bar{u}}^{(s)} (x, y; p_{1}, q_{1})

is related to the complex scattering potential

f (X, Y, Z)

of the specimen as

{\bar{u}}^{(s)} (x, y; p_{1}, q_{1}) = \iint {[i 4 π (W + m_{1} / λ)]}^{- 1} \tilde{f} (U, V, W) \exp [i 2 π (U x + V y + W R_{2})] d U d V,

where the scattering potential

f (X, Y, Z)

is a function of the complex refractive index defined as

f (X, Y, Z) = k^{2} (1 - {(n (X, Y, Z) / n_{0})}^{2})

, and the function

\tilde{f} (U, V, W)

is the 3-D Fourier transform of

f (X, Y, Z)

. Note that the variables

(U, V, W)

are related by

W = \sqrt{{(1 / λ)}^{2} - {(U + p_{1} / λ)}^{2} - {(V + q_{1} / λ)}^{2}} - m_{1} / λ

.

The intensity recorded on the detector can be written as:

I (x, y) = \iint I_{0} (x_{0}, y_{0}) {| U (x, y; x_{0}, y_{0}) |}^{2} d x_{0} d y_{0} .

Substituting the total field

U (x, y; x_{0}, y_{0})

from Eq. (4), and using the convolution theorem, Eq. (6) can be rewritten as:

I (x, y) = \iint \iint {\tilde{I}}_{0} (λ^{- 1} (p_{1} - {\bar{p}}_{1}), λ^{- 1} (q_{1} - {\bar{q}}_{1})) u (x, y; p_{1}, q_{1}) u^{*} (x, y; {\bar{p}}_{1}, {\bar{q}}_{1}) d p_{1} d q_{1} d {\bar{p}}_{1} d {\bar{q}}_{1} .

Equation (7) is a master equation that can be applied to the apertures of an arbitrary size and shape.

2.2 Coherent limit

First, let us consider the coherent limit of Eq. (7), where we assume that an infinitesimally small x-ray source is located on the optical axis and emits a simple cone-shaped x-ray beam. Therefore, the Fourier transform of the source strength S₀ can be written as

{\tilde{I}}_{0} (u, v) = S_{0} .

From Eq. (7), the intensity recorded on the detector can be simplified to:

I (x, y) = S_{0} {| \iint u (x, y; p_{1}, q_{1}) d p_{1} d q_{1} |}^{2} .

Note that the function

u (x, y; p_{1}, q_{1})

contains a phase term that oscillates very fast as

p_{1}

and

q_{1}

are varied. The step size for the evaluation of this integral scales down with the wavelength; therefore, the computational cost for direct evaluation of Eq. (9) can be prohibitively large, especially in x-ray imaging. Instead, this integral can be evaluated using the method of stationary phase [28], which is valid when

k r \to \infty

. In a typical x-ray geometry,

k r

is usually as large as

10^{11}

; therefore, the assumption

k r \to \infty

is quite valid. Using the method of stationary phase, the following asymptotic formula can be obtained:

\iint u (x, y; p_{1}, q_{1}) d p_{1} d q_{1} ~ (1 / r) \exp (i k r) \bar{u} (x, y; x / r, y / r),

where

\bar{u} (x, y; x / r, y / r)

is

1 + {\bar{u}}^{(s)} (x, y; x / r, y / r)

according to the first Born approximation and

\exp [{\bar{u}}^{(s)} (x, y; x / r, y / r)]

according to the first Rytov approximation. Here, the variable

r

is defined as

r = \sqrt{x^{2} + y^{2} + {(R_{1} + R_{2})}^{2}}

. After replacing the integral in Eq. (9) with Eq. (10), the intensity of the scattered field in the coherent limit can be written as:

First Born approximation : I (x, y) = (I_{0} / r^{2}) {| 1 + {\bar{u}}^{(s)} (x, y; x / r, y / r) |}^{2},

First Rytov approximation : I (x, y) = (I_{0} / r^{2}) \exp [2 Re {{\bar{u}}^{(s)} (x, y; x / r, y / r)}] .

2.3 Arbitrary source aperture

Returning to the master equation Eq. (7), let us consider the case of a finite aperture with arbitrary shape. As in the coherent limit case, direct evaluation of the integral in Eq. (7) is practically impossible because of the rapidly varying phase terms. For simplicity, we change the variables in Eq. (7) as

λ^{- 1} (p_{1} - {\bar{p}}_{1}) = P, λ^{- 1} (q_{1} - {\bar{q}}_{1}) = Q, λ^{- 1} (m_{1} - {\bar{m}}_{1}) = M .

\begin{array}{l} I (x, y) = \iint \iint {\tilde{I}}_{0} (P, Q) \exp {i 2 π [P x + Q y + M (R_{2} + R_{1})]} \\ \times m_{1}^{- 1} {(m_{1} - λ M)}^{- 1} \bar{u} (x, y; p_{1}, q_{1}) {\bar{u}}^{*} (x, y; p_{1} - λ P, q_{1} - λ Q) d p_{1} d q_{1} d P d Q . \end{array}

In a typical XPI arrangement, the size of the source aperture

I_{0} (x_{0}, y_{0})

is on the order of 100 μm; therefore, the cut-off spatial frequency of

{\tilde{I}}_{0} (P, Q)

is around 10⁴ m⁻¹. The variables

p_{1}

,

q_{1}

,

{\bar{p}}_{1}

, and

{\bar{q}}_{1}

are connected to the divergence angle of x-ray beams, which is on the order of 0.1 m⁻¹. The wavelength

λ

is less than 10⁻¹¹ m. Since the variables

λ P

and

λ Q

are about 10⁶ times smaller than

p_{1}

and

q_{1}

, the following approximations can be made safely:

\begin{array}{l} {(m_{1} - λ M)}^{- 1} ≅ m_{1}^{- 1}, \bar{u} (x, y; p_{1} - λ P, q_{1} - λ Q) ≅ \bar{u} (x, y; p_{1}, q_{1}), \\ M = λ^{- 1} (m_{1} - {\bar{m}}_{1}) ≅ - (p_{1} / m_{1}) P - (q_{1} / m_{1}) Q + ϑ (λ P^{2} / m_{1}) . \end{array}

The variables determining the phase, e.g.,

M

, usually need more caution when the approximations as above are made. It is because a small difference in the phase can make a big change in the calculated intensity as in the case of destructive interference. However, in the integrand of Eq. (14) the phase error due to the ignored term

ϑ (λ P^{2} / m_{1})

is on the order of 1 mrad, which can be safely ignored. Introducing these approximations, we can simplify Eq. (14) to:

I (x, y) = \iint I_{0} (x - x_{p}, y - y_{p}) (1 / m_{1}^{2}) {| \bar{u} (x, y; p_{1}, q_{1}) |}^{2} d p_{1} d q_{1},

where

x_{p} = (p_{1} / m_{1}) (R_{2} + R_{1})

,

y_{p} = (q_{1} / m_{1}) (R_{2} + R_{1})

, and

m_{1} = \sqrt{1 - p_{1}^{2} - q_{1}^{2}}

.

Changing the variables of integration from $(p_{1}, q_{1})$ to $(x_{p}, y_{p})$ , we can obtain the following formulae for the intensity distribution in the case of an arbitrary source aperture:

First Born approximation:

I (x, y) = \iint (1 / r_{p}^{2}) I_{0} (x - x_{p}, y - y_{p}) {| 1 + {\bar{u}}^{(s)} (x, y; x_{p} / r_{p}, y_{p} / r_{p}) |}^{2} d x_{p} d y_{p},

First Rytov approximation:

I (x, y) = \iint (1 / r_{p}^{2}) I_{0} (x - x_{p}, y - y_{p}) \exp [2 Re {{\bar{u}}^{(s)} (x, y; x_{p} / r_{p}, y_{p} / r_{p})}] d x_{p} d y_{p},

where

r_{p} = \sqrt{x_{p}^{2} + y_{p}^{2} + {(R_{1} + R_{2})}^{2}}

.

We note that when the focal spot size of the source gets smaller, as in a micro-focus source, the magnitude of the term $ϑ (λ P^{2} / m_{1})$ increases to the order of 0.1 rad and the validity of the approximation for $M$ may be undermined. However, in the limiting case of an infinitesimally small source, which corresponds to substituting $I_{0} (x, y) = I_{0} δ (x, y)$ in Eqs. (17) and (18), we obtain the formulas for the coherent limit in Eqs. (11) and (12), respectively. Considering that the latter two equations were obtained with the stationary phase approximation, which can be safely justified in XPI, we may assume that the approximations made for Eqs. (17) and (18) are as valid for a small source as for a large one.

3. Results and Discussion

3.1 Near-field x-ray images of multiple size spheres

In this section, we calculate the near-edge fringe patterns of homogeneous spheres (beads) to investigate the effect of object size on the contrast of signal and required detector resolution in propagation-based XPI. We assume that the source is infinitesimally small (i.e., a spatially-coherent source) and a monochromatic beam of 30 keV. The source-to-object distance $R_{1}$ and the object-to-detector distance $R_{2}$ are fixed at 1 m, while the radius of sphere is varied from 10 μm to 1 mm. The corresponding Fresnel number $F = a^{2} / λ R_{2}$ , which is based on the radius $a$ of the sphere and the object-to-detector distance $R_{2}$ , ranges from 2.4 to 2.4 × 10⁴. For the complex refractive index of the sphere, we use δ = 3.516 × 10⁻⁷ and β = 1.233 × 10⁻¹⁰ referring to the value for water at the given source energy [30]. Under these conditions, the error due to the Rytov approximation is negligible [23].

Figures 2 and Fig 3 compare x-ray diffraction patterns of homogeneous beads of different radii from 10 μm to 1 mm. For the smallest bead, the intensity image clearly shows the diffracting nature of an x-ray beam. This effect is localized to the boundary of bead and is less visible when the size of bead is increased. For the bead of 500 μm radius, under the given test conditions the energy of diffracted beam is concentrated over a 7 μm wide rim around the edge of the beam. In practice, this sharp edge response will be blurred due to a large focal spot size, a broad spectrum of the source, and filtering by the modulation transfer function of the detector. When there is a large absorption background, this edge response may not be visible above the background due the limited dynamic range of the detector. Image noise will further blunt this edge response. In propagation-based imaging, the easiest way to enhance this edge response is to increase the object-to-detector distance. This strategy, however, comes at the cost of reduced beam flux, increased detector area due to geometric magnification, and overall larger system footprint.

Fig. 2 Near-field x-ray diffraction patterns of homogeneous, spherical water beads of different radii: (a) 10 μm, (b) 100 μm, and (c) 1 mm. The images were numerically generated assuming a cone-beam x-ray set-up with an infinitesimally small focal spot and 30keV monochromatic x-ray beam. The source-to-object distance R₁ and the object-to-detector distance R₂ are each 1 m.

Download Full Size | PDF

Fig. 3 Near-field x-ray diffraction patterns of homogeneous water beads of different radii. (a) Intensity profile across the center of bead for the bead radii: 10, 20, 50, and 100 μm; (b) near-edge portion of the intensity profile for two bead radii (30keV monochromatic x-ray beam): 500 μm and 1 mm. The detector coordinate (horizontal axis) was scaled so that the edges of the bead are located at ± 1.

Download Full Size | PDF

3.2 Simulation with a discretized object model

An XPI simulator requires the complex refractive index of an object as an input. In the previous section, we used an analytic expression for the sphere [29] with published values of refractive index [30]. Alternatively, a more complex object model can be generated using computer-aided design (CAD) software [20]; or it can be acquired with crystal analyzed-based computed tomography (CT) [31] or dual energy x-ray CT [32]. The object model generated in either way is typically discretized into pixels or voxels with sharp edges. Such discretization can generate artifacts in simulated phase images such as diffraction patterns from the sharp voxel edges [20]. A simple way to mitigate this problem is to decrease the pixel size, thus increasing the size of an object model and, consequently, the computation time that scales as O(n³) where n is the number of voxels. Instead, the artifact can be effectively suppressed by applying a smoothing filter to the object model. To demonstrate this, we generated a discretized object model by modifying a 3-D Shepp-Logan phantom [33] and calculated phase images of the phantom with and without a smoothing filter. Here we model a 30 keV point x-ray source located far upstream so that the incident beam can be assumed to be a plane wave. Both the object model and detector plane were discretized into 100 μm pixels.

Figure 4(a) shows real and imaginary parts of a horizontal cross-section of the 3-D phantom. The refractive index values of the phantom were determined referring to those of water for the given energy 30 keV [30]. Figures 4(b) and 4(c) show the image calculated in the detector plane at 2 m from the center of the phantom. Without a smoothing filter, the calculated image is contaminated with ringing artifact due to the discretization of the object model. The following smoothing filter was applied to the object model before processing it with our phase propagation system:

h (x, y, z) = {(2 π σ^{2})}^{- 3 / 2} \exp [- (x^{2} + y^{2} + z^{2}) / 2 σ^{2}],

Here the FWHM of the Gaussian profile along each axis is about 2.35σ. Figures 4(d) and 4(e) show the image calculated from the object model after applying the Gaussian smoothing kernel of Eq. (19) (FWHM = 200 μm). As clearly shown, most ringing artifacts are significantly suppressed by the Gaussian filter. The figure correctly displays both the absorption (reduced beam flux within the phantom) and phase signatures (enhanced edges) of the phantom. We note that as the FWHM of the applied Gaussian kernel increases, the image calculated at the image plane will be more smoothened.

Fig. 4 Simulation with a modified 3-D Shepp-Logan phantom: (a) real (δ) and imaginary (β) parts of a horizontal cross-section of the refractive index phantom shown in 3-D in the image below the cross-sectional image. (The 3-D image was obtained by modifying the code from https://sites.google.com/site/hispeedpackets/Home/shepplogan.) (b, d) simulated x-ray images without (b) and with (d) applying Gaussian smoothing to the phantom. The FWHM of the applied Gaussian kernel is about 200 μm along each axis. (c) and (e) are the zoom-in views of (b) and (d), respectively.

Download Full Size | PDF

Next, we demonstrate our method using the refractive index map acquired with x-ray dark field imaging (XDFI), a crystal analyzer-based method. The experimental setup comprises an x-ray source, an asymmetrically-cut Bragg-type monochromator-collimator (MC), a Laue-case angle analyzer (LAA) and a CCD camera (X-FDI 1.00:1, Photonic Science Ltd.). For the light source, we use the beam line BL14C on a 2.5 GeV storage ring (KEK Photon Factory, Tsukuba, Japan), and for LAA we use a Si crystal of 150 μm thickness. The specimen is placed between the MC and LAA. More detailed explanation of this method can be found elsewhere [34]. We prepared an iliac artery specimen fixed with formalin and placed it in a tube filled with formalin. XDFI can provide the absorption (β) and phase (δ) parts of the complex refractive index map of the specimen (Fig. 5(a)). Example cross-sections are shown in Fig. 5(b). For each cross-section, 600 projection images were processed using a tomographic reconstruction algorithm [35]. We note that due to dehydration after dissection, the tissue specimen shrank creating an empty space, which is replaced by formalin when the specimen is placed in a plastic tube for imaging. These regions are shown in dark grey (d) in Fig. 5(b) (ii). In the absorption image (i) only the wall of plastic tube (a) is vaguely visible, while in the phase image (ii) an artery wall (c) and the tissue shrinkage regions (d) are clearly visible. The plastic rod (b), which was inserted for comparison, is shown to be highly uniform inside, which proves the fidelity of this method. Figure 6(i) shows a simple projection of the phase map, which can be connected under the projection approximation to the phase profile of distorted wave front after the sample. Figures 6(ii)-(iv) show the x-ray phase images simulated with the model proposed in this study for various focal spot size of the source: (ii) 20 μm, (iii) 100 μm, and (iv) 500 μm. The source-to-object and object-to-detector distance were fixed at 0.2 m. For the 20 μm focal spot, which may be produced with a micro-focus source [36], the edges of plastic tube, plastic rod, and artery wall are clearly visible. The space formed by tissue shrinkage is filled with formalin whose refractive index is very different from the surrounding tissue. Therefore, small structures of the shrunk tissue generate strong diffraction, which is shown as speckle patterns (d) dispersed in Fig. 6(ii). For the 100 μm focal spot, which may be produced with a mammographic x-ray tube [37] or a carbon nanotube field emission x-ray source [38], the edges of plastic tube, plastic tube, and artery wall are still visible although these edges are shown blurred due to penumbral smoothing. The speckle pattern due to the shrunk tissue parts is also clearly distinguished from the other background region. For the 500 μm focal spot corresponding to typical x-ray tubes used in absorption-based imaging, all these features (a-d) are significantly blunted.

Fig. 5 Refractive index map of an iliac artery sample acquired with crystal analyzer-based CT: (a) 3-D rendered image of the data cube; (b) Sample cross-sections of the complex refractive index map: (i) absorption (β), and (ii) phase (δ) part. a: plastic tube; b: plastic rod (inserted for comparison); c: three-layer artery wall; and d: tissue shrinkage.

Download Full Size | PDF

Fig. 6 X-ray phase imaging simulation using the complex refractive index map in Fig. 5: (i) Simple projection of the phase map δ; (ii)-(iv): x-ray phase images simulated with the proposed model for various focal spot size of the source (20 μm, 100 μm, and 500 μm). A cone-shaped x-ray beam illuminates the sample from the side of the data cube shown in Fig. 5. The source-to-object and object-to-detector distance were fixed at 0.2 m.

Download Full Size | PDF

4. Conclusion

In this study, we have developed a forward model based on scalar diffraction theory for x-ray phase contrast imaging. Our model, which uses the first Rytov approximation, is robust and can handle optically thick objects. It can also be used with a source of any size, shape, and energy spectrum. By not relying on the projection or paraxial approximations, the model can be applied to large objects or objects located off-axis in a cone-beam x-ray setup. By unifying the steps for x-ray-object interaction and free-space propagation, the model can be conveniently adopted for an inverse problem as well as for an x-ray phase imaging simulator.

Acknowledgments

This work was supported by the Department of Homeland Security’s Science and Technology Directorate through contract HSHQDC-11-C-0083, the National Research Foundation of Singapore through the Singapore-MIT Alliance for Research and Technology Centre, and the DARPA AXiS program (Grant No. N66001-11-4204, P.R. No. 1300217190). The experiment was performed under the approval of the PF User Association (PF-UA) at KEK under No. 2008S2-002, 2011G-672 for use of the Photon Factory. The authors thank Dr. Synho Do and Dr. Omid Khalilzadeh for helpful discussion.

References and links

1. A. Stanton, “Wilhelm Conrad Röntgen on a new kind of rays: translation of a paper read before the Würzburg Physical and Medical Society, 1895,” Nature 53, 274–276 (1896).

2. E. D. Pisano, M. J. Yaffe, and C. M. Kuzmiak, Digital Mammography (Lippincott Williams & Wilkins, 2004).

3. B. Henke, E. Gullikson, and J. C. Davis, “X-Ray Interactions: Photoabsorption, Scattering, Transmission, and Reflection at E = 50-30,000 eV, Z = 1-92,” At. Data Nucl. Data Tables 54(2), 181–342 (1993). [CrossRef]

4. A. Momose, T. Takeda, and Y. Itai, “Blood Vessels: Depiction at Phase-Contrast X-ray Imaging without Contrast Agents in the Mouse and Rat-Feasibility Study 1,” Radiology 217(2), 593–596 (2000). [PubMed]

5. M. Ando, A. Maksimenko, H. Sugiyama, W. Pattanasiriwisawa, K. Hyodo, and C. Uyama, “Simple x-ray dark-and bright-field imaging using achromatic Laue optics,” Jpn. J. Appl. Phys. 41(Part 2, No. 9A/B), L1016–L1018 (2002). [CrossRef]

6. U. Bonse and F. Beckmann, “Multiple-beam X-ray interferometry for phase-contrast microtomography,” J. Synchrotron Radiat. 8(1), 1–5 (2001). [CrossRef] [PubMed]

7. T. Weitkamp, C. David, O. Bunk, J. Bruder, P. Cloetens, and F. Pfeiffer, “X-ray phase radiography and tomography of soft tissue using grating interferometry,” Eur. J. Radiol. 68(3Suppl), S13–S17 (2008). [CrossRef] [PubMed]

8. F. Pfeiffer, M. Bech, O. Bunk, P. Kraft, E. F. Eikenberry, Ch. Brönnimann, C. Grünzweig, and C. David, “Hard-X-ray dark-field imaging using a grating interferometer,” Nat. Mater. 7(2), 134–137 (2008). [CrossRef] [PubMed]

9. S. W. Wilkins, T. E. Gureyev, D. Gao, A. Pogany, and A. Stevenson, “Phase-contrast imaging using polychromatic hard X-rays,” Nature 384(6607), 335–338 (1996). [CrossRef]

10. P. Spanne, C. Raven, I. Snigireva, and A. Snigirev, “In-line holography and phase-contrast microtomography with high energy x-rays,” Phys. Med. Biol. 44(3), 741–749 (1999). [CrossRef] [PubMed]

11. Y. S. Kashyap, P. S. Yadav, T. Roy, P. S. Sarkar, M. Shukla, and A. Sinha, “Laboratory-based X-ray phase-contrast imaging technique for material and medical science applications,” Appl. Radiat. Isot. 66(8), 1083–1090 (2008). [CrossRef] [PubMed]

12. Z. Zaprazny, D. Korytar, V. Ac, P. Konopka, and J. Bielecki, “Phase contrast imaging of lightweight objects using microfocus X-ray source and high resolution CCD camera,” JINST 7(03), C03005 (2012). [CrossRef]

13. D. M. Paganin, Coherent X-ray Optics (Oxford University Press, New York, 2006).

14. M. Ando, K. Yamasaki, F. Toyofuku, H. Sugiyama, C. Ohbayashi, G. Li, L. Pan, X. Jiang, W. Pattanasiriwisawa, D. Shimao, E. Hashimoto, T. Kimura, M. Tsuneyoshi, E. Ueno, K. Tokumori, A. Maksimenko, Y. Higashida, and M. Hirano, “Attempt at visualizing breast cancer with x-ray dark field imaging,” Jpn. J. Appl. Phys. 44(17), L528–L531 (2005). [CrossRef]

15. R. M. Aspden and D. W. L. Hukins, “Collagen organization in articular cartilage, determined by X-ray diffraction, and its relationship to tissue function,” Proc. R. Soc. Lond. B Biol. Sci. 212(1188), 299–304 (1981). [CrossRef] [PubMed]

16. R. A. Lewis, “Medical phase contrast x-ray imaging: current status and future prospects,” Phys. Med. Biol. 49(16), 3573–3583 (2004). [CrossRef] [PubMed]

17. D. Shimao, H. Sugiyama, T. Kunisada, and M. Ando, “Articular cartilage depicted at optimized angular position of Laue angular analyzer by X-ray dark-field imaging,” Appl. Radiat. Isot. 64(8), 868–874 (2006). [CrossRef] [PubMed]

18. A. Pogany, D. Gao, and S. W. Wilkins, “Contrast and resolution in imaging with a microfocus x-ray source,” Rev. Sci. Instrum. 68(7), 2774 (1997). [CrossRef]

19. X. Wu and H. Liu, “A general theoretical formalism for X-ray phase contrast imaging,” J. XRay Sci. Technol. 11(1), 33–42 (2003). [PubMed]

20. A. Peterzol, J. Berthier, P. Duvauchelle, C. Ferrero, and D. Babot, “X-ray phase contrast image simulation,” Nucl. Instrum. Methods Phys. Res., Sect. B 254, 307–318 (2007).

21. L. A. Chernov and R. A. Silverman, Wave Propagation in a Random Medium (McGraw-Hill, 1960)

22. V. I. Tatarski, Wave Propagation in a Turbulent Medium (McGraw-Hill, 1961).

23. Y. Sung and G. Barbastathis, “Rytov approximation for x-ray phase imaging,” Opt. Express 21(3), 2674–2682 (2013). [CrossRef] [PubMed]

24. E. Wolf, “Three-dimensional structure determination of semi-transparent objects from holographic data,” Opt. Commun. 1(4), 153–156 (1969). [CrossRef]

25. A. C. Kak and M. Slaney, Principles of Computerized Tomographic Imaging (SIAM, 1988).

26. Y. Sung, W. Choi, C. Fang-Yen, K. Badizadegan, R. R. Dasari, and M. S. Feld, “Optical diffraction tomography for high resolution live cell imaging,” Opt. Express 17(1), 266–277 (2009). [CrossRef] [PubMed]

27. A. J. Devaney, “Inverse-scattering theory within the Rytov approximation,” Opt. Lett. 6(8), 374–376 (1981). [CrossRef] [PubMed]

28. M. Born, E. Wolf, and A. B. Bhatia, Principles of Optics: Electromagnetic Theory of Propagation, Interference and Diffraction of Light (Cambridge Univ Press, 1999).

29. E. M. Stein and G. L. Weiss, Introduction to Fourier Analysis on Euclidean Spaces (Princeton University Press, Princeton, 1971).

30. B. Henke, E. Gullikson, and J. C. X. Davis, “X-Ray Interactions: Photoabsorption, scattering, transmission, and reflection at E = 50-30,000 eV, Z = 1-92,” At. Data Nucl. Data Tables 54(2), 181–342 (1993). [CrossRef]

31. N. Sunaguchi, T. Yuasa, Q. Huo, S. Ichihara, and M. Ando, “Refraction-contrast tomosynthesis imaging using dark-field imaging optics,” Appl. Phys. Lett. 99(10), 103704 (2011). [CrossRef]

32. R. E. Alvarez and A. Macovski, “Energy-selective reconstructions in x-ray computerized tomography,” Phys. Med. Biol. 21(5), 733–744 (1976). [CrossRef] [PubMed]

33. M. Schabel, “3D Shepp-Logan phantom,” MATLAB Central File Exchange (2006).

34. M. Ando, N. Sunaguchi, Y. Wu, S. Do, Y. Sung, A. Louissaint, T. Yuasa, S. Ichihara, and R. Gupta, “Crystal Analyser-based X-ray Phase Contrast Imaging in the Dark Field: Implementation and Evaluation using Excised Tissue Specimens,” Invest. Radiol. under review.

35. N. Sunaguchi, T. Yuasa, Q. Huo, and M. Ando, “Convolution reconstruction algorithm for refraction-contrast computed tomography using a Laue-case analyzer for dark-field imaging,” Opt. Lett. 36(3), 391–393 (2011). [CrossRef] [PubMed]

36. L. Poletto, M. Caldon, G. Tondello, and A. Megighian, “A system for high-resolution x-ray phase-contrast imaging and tomography of biological specimens,” Proc. SPIE 7078, 70781P, 70781P-10 (2008). [CrossRef]

37. T. Tanaka, C. Honda, S. Matsuo, K. Noma, H. Oohara, N. Nitta, S. Ota, K. Tsuchiya, Y. Sakashita, A. Yamada, M. Yamasaki, A. Furukawa, M. Takahashi, and K. Murata, “The first trial of phase contrast imaging for digital full-field mammography using a practical molybdenum x-ray tube,” Invest. Radiol. 40(7), 385–396 (2005). [CrossRef] [PubMed]

38. G. Cao, Y. Z. Lee, R. Peng, Z. Liu, R. Rajaram, X. Calderón-Colon, L. An, P. Wang, T. Phan, S. Sultana, D. S. Lalush, J. P. Lu, and O. Zhou, “A dynamic micro-CT scanner based on a carbon nanotube field emission x-ray source,” Phys. Med. Biol. 54(8), 2323–2340 (2009). [CrossRef] [PubMed]

Full-wave approach for x-ray phase imaging

Abstract

1. Introduction

2. Full-wave approach for x-ray phase imaging

2.1 Scalar diffraction theory

2.2 Coherent limit

2.3 Arbitrary source aperture

3. Results and Discussion

3.1 Near-field x-ray images of multiple size spheres

3.2 Simulation with a discretized object model

4. Conclusion

Acknowledgments

References and links

Cited By

Figures (6)

Equations (19)

Optics Express