Forward model for propagation-based x-ray phase contrast imaging in parallel- and cone-beam geometry

Elisabeth R. Shanblatt; Yongjin Sung; Rajiv Gupta; Brandon J. Nelson; Shuai Leng; William S. Graves; Cynthia H. McCollough

doi:10.1364/OE.27.004504

1. Introduction

Propagation-based X-ray phase-contrast imaging (PB-XPCI) is a promising technique that has shown exciting results in pre-clinical studies [1–3]. This technique involves illuminating a sample with coherent X-rays and allowing the scattered wave to propagate in free space at various distances, and then using a forward model for wave propagation (e.g. the transport-of-intensity equation (TIE) [4–6]) to recover the phase of the object’s exit surface wave. PB-XPCI has found a wide range of applications across various length scales, particularly for computed tomography (CT). Applications include clinical-scale large field-of-view (FOV) studies of joint and bone structure [7,8], breast imaging [9], brain imaging [10,11], and cardiac applications [12]. PB-XPCI is also used for high resolution nanotomography [13–15]. Numerous systems have been studied with phase contrast nanotomography, including materials [16,17], and biological systems such as lungs [18], bone ultrastructure [19,20], neurological systems [21], and the organism Caenorhabditis elegans [22].

Forward models of PB-XPCI are important for system prototyping and design, as well as for iterative tomographic reconstruction. Current forward models suffer from limitations on the type of objects that can be modeled, and can be computationally expensive. Non-paraxial forward models such as the Mie solution [23–25] and full-wave Rytov approach [26–29] provide very high accuracy for parallel-beam geometry, but are not able to simultaneously incorporate multiple scattering, arbitrarily large objects, and cone-beam effects.

A simple method for modeling a three-dimensional object is to treat it as a two-dimensional projection of its complex refractive index along the propagation axis. We refer to this as the projected-object approximation (rather than the projection approximation as it is sometimes called) to distinguish from when a projection is assumed along the entire ray path, which is a routine assumption in many tomographic reconstruction methods. The phase and amplitude changes to a wave caused by an object treated as a projection can be expressed as [30],

Δ φ (x, y) = - \frac{2 π}{λ} \int_{0}^{T} δ (x, y, z) d z

a (x, y) / a_{0} (x, y) = exp [- \frac{2 π}{λ} \int_{0}^{T} β (x, y, z) d z]

where the change in phase

Δ φ

, induced by the object over its thickness

T

, is simply the integral of the real part of the spatially varying refractive index

δ (x, y, z)

, scaled by the wavelength

λ

. The quantity

a (x, y) / a_{0} (x, y)

is the ratio of the amplitude at the entrance plane of the object to the amplitude at the exit plane of the object, and

β (x, y, z)

is the imaginary part of the object’s spatially varying refractive index. The full refractive index is expressed then as

n (x, y, z) = [1 - δ (x, y, z)] + i β (x, y, z)

.

While the projected-object approximation described in Eqs. (1) and (2) is useful for many imaging systems, it begins to break down once diffraction effects within the sample become important. To account for intra-object scattering and diffraction in optical and certain X-ray applications, the beam propagation method (BPM) [31,32] has been used [33,34], which utilizes a multi-slice approach and does not require that the object for a given imaging system be valid under the projected-object approximation [35–37]. BPM has been used in X-ray ptychography [36,38] and other applications [39,40], where the parallel-beam assumption is a good approximation. However, current implementations of BPM do not account for cone-beam effects within the sample. Cone-beam effects are important in many applications where imaging is performed with a diverging (or converging) beam [14], and is of crucial importance for nanotomography [15,41–44]. Cone-beam geometry is also important in medical X-ray radiography and CT [45] due to requirements of wide beam coverage and short scan times.

Furthermore, iterative methods are important for both tomographic reconstructions and phase retrieval applications [46,47]. An accurate and fast forward model is crucial for developing iterative reconstruction suitable for coherent phase-contrast tomography. Thus, our motivation here is two-fold: firstly, an iterative algorithm can only be as good as its forward model, and as short computation times are essential for high-throughput and clinical feasibility, we choose to forgo high-angle accuracy in favor of a paraxial model for the flexibility to incorporate arbitrary objects and cone-beam effects. Secondly, an accurate forward model can be used for imaging system design and prototyping, for instance in determining tradeoffs between phase sensitivity, FOV, and resolution.

Here, we combine BPM with cone-beam effects by rescaling the diffracted wave at multiple planes within the object, rather than treating the object under the projected-object approximation and only rescaling once. Hence, we call our method the cone-beam beam propagation method, or CB-BPM. We make this distinction from the traditional approach by referring to it as parallel-beam BPM, or PB-BPM.

2. Forward models used for propagation-based X-ray phase-contrast imaging

The appropriate choice of forward model for a scattering process depends on the scattering potential, type of illumination, and competing needs of accuracy and low computational expense. This section discusses some of the forward models in present use for PB-XPCI and the types of objects for which they are most appropriate.

Many current implementations of PB-XPCI rely on the assumption that the object can be treated under the projected-object approximation (Eqs. (1) and (2)). This assumption is valid for thin and weakly scattering objects, and for other modalities such as incoherent X-ray absorption imaging used in clinical CT. However, for coherent phase imaging, and particularly for imaging thick objects, diffraction within the object can be significant. Figure 1 demonstrates object and beam geometries used in current PB-XPCI forward models.

Fig. 1 Scattering geometries in propagation-based X-ray phase contrast imaging. (a) Parallel-beam single-scattering geometry, often used in the Mie solution. (b) Parallel-beam beam propagation method (PB-BPM), which does not account for cone beam effects but does allow multiple scattering. (c) Cone-beam model that only allows a single scattering event, as the object is treated under the projected-object approximation. (d) The model proposed here, CB-BPM, uses cone-beam geometry and allows for a thick, multi-scattering object to be considered.

Download Full Size | PDF

2.1 Mie solution

The Mie solution (sometimes referred to as the Lorenz-Mie solution) is an exact solution of Maxwell’s equations in the form of an infinite series, making it an ideal candidate against which to test other models. While it was originally demonstrated for a plane wave scattering from a single homogeneous sphere [23], solutions have recently been demonstrated for a wide range of particle numbers, sizes, shapes, and illumination types. One approach, known as generalized Lorenz-Mie theory, allows exact solutions for scattering by non-plane wave illumination such as Gaussian beams [25,48,49]. Another approach, T-Matrix formulation [24,50,51], can account for multiple scattering caused by aggregates of particles. While these approaches are very powerful and highly general, their convergence is not well-established for near-field solutions of multiple scattering [52].

For certain types of scattering assemblies, the scattered field can be found by solving by the wave equation with appropriate boundary conditions, which is particularly important for convergence and accuracy of near-field solutions [53,54]. However, this type of solution is only computationally feasible for very limited types of systems. In addition to single spheres, solutions can be found for one or more infinite cylinders [53].

While the Mie solution’s accuracy makes it an excellent gold standard against which to compare other methods, its long computation time is problematic. The computation time scales as a ratio of the object’s diameter to the wavelength [26], which is a problem particularly for X-ray imaging of human-sized objects, where this ratio can exceed $10^{10}$ . It is also impractical to model large, complex biological structures as particle aggregates.

2.2 First-order solutions

The full-wave Rytov (FWR) approach uses the scalar wave theory for the propagation of X-rays across the entire object’s volume as well as the free space between the object and detector [26,27]. This contrasts with other existing approaches using the projected-object approximation for the object and applying a wave approach only for the free-space propagation. FWR can incorporate a phantom defined as non-uniform rational B-splines (NURBS), which enabled XPCI simulation using the realistic XCAT phantom without numerical artifacts due to discretized voxels [28]. FWR can be applied to cone-beam geometry, but it is computationally expensive. FWR does not account for multiple scattering, which may be important in some applications.

While not yet widely used in the X-ray regime, much recent work has been done on nonlinear first-order methods for reconstruction in optical diffraction tomography, which is a related problem to what is treated here, but typically uses visible light. Some of these approaches include inversion methods [55–57], which can account for multiple scattering and reflections. Recently, applications of machine learning [34,58] have been demonstrated for tomographic reconstruction, as have other methods for phase retrieval [59,60].

2.3 Beam propagation method

The multi-slice PB-BPM has been used for a variety of applications in optics. Although it is capable of accounting for a wide variety of objects, it does not incorporate geometric cone-beam effects, which are important for cone-beam CT and other diverging-beam applications where geometric magnification is non-negligible.

BPM does not require that an entire object be valid under the projected-object approximation. Rather, each individual slice is approximated as a projection, with propagation occurring in between each slice. This allows for multiple scattering events and diffraction within the object, so long as the slices are chosen to be sufficiently thin such that the refractive index does not vary significantly within an individual slice along the propagation direction. For real-space imaging, a convenient choice of slice thickness is to set the resolution to be isotropic, or an integer number of the transverse resolution. While an interpolation between slices can potentially be performed, the fundamental resolution limit is set by the camera’s pixel size and the system magnification.

In the case of diffraction imaging, where the resolution is typically much higher and the FOV typically much smaller, the minimum diffraction-limited resolution is set by the illuminating wavelength, detector size, and propagation distance of the system. It therefore may be more convenient to set the axial resolution to a larger number. BPM relies heavily on the fast Fourier transform (FFT) algorithm, and therefore the computation time generally depends on the number of points in the array, as well as the number of slices. Physical parameters such as wavelength typically do not change the computation time of this method, although they will affect the distance over which the projected-object approximation is valid.

One drawback of BPM is that in most implementations, the paraxial approximation is assumed to hold for the wave. This means that high-angle scattering events and wave curvature are not incorporated into the model. In the context of PB-XPCI, many implementations of the transport of intensity equation (TIE) already assume paraxial scattering [4,61,62]. Then, the application of CB-BPM requires no further approximations. The following section discusses the paraxial approximation in further detail.

3. Cone-beam beam propagation method

3.1 The paraxial approximation

While high angle wave information is sacrificed in the BPM approach, numerous steps in the PB-XPCI process routinely make the paraxial approximation. Many implementations of the TIE rely on the paraxial approximation in order to perform phase retrieval [5], as does the Fresnel scaling theorem, which is used to model cone-beam effects [30]. The paraxial approximation can be stated as [63],

\begin{array}{l} | λ f_{x} | ≪ 1 \\ | λ f_{y} | ≪ 1 \end{array}

where

λ

is the wavelength and

f_{x}

and

f_{y}

are the spatial frequencies in the

x

and

y

dimensions respectively. In an image-space (as opposed to Fourier space) near-field imaging approach such as PB-XPCI, the spatial frequencies are approximately the same at the sample and detector planes, and are typically set by the pixel size and overall dimension of the detector, as this determines how the spatial frequencies of the wave are sampled. The highest possible product on the left-hand sides of Eq. (3) will be given by the largest possible spatial frequency the imaging system can sample, which is determined by the inverse of the detector pixel size. Photon energies as low as 1 keV still have wavelengths of the order 1-2 nm, so even an imaging system with very small detector pixel size (of the order 1 μm) will produce a product several orders of magnitude below unity. Thus, an image-space near-field imaging system satisfies the paraxial approximation, Eq. (3).

Other imaging techniques utilizing X-rays, such as coherent diffractive imaging (CDI) use a very tightly focused beam to allow data to be collected in Fourier space [64]. CDI is capable of sampling very high spatial frequencies compared to the wavenumber $1 / λ$ , and therefore may not be valid under the paraxial approximation.

3.2 Fresnel scaling theorem and PB-BPM

The Fresnel scaling theorem is a convenient way to incorporate magnification into a computational wave simulation. It maps the intensity $I_{C B}$ obtained by cone- beam propagation in terms of the intensity that would be obtained by parallel- beam propagation, $I_{P B}$ ,with a different propagation distance and pixel size [30],

I_{C B} (x, y, z_{2}) = \frac{1}{M^{2}} I_{P B} (\frac{x}{M}, \frac{y}{M}; z_{e f f}) .

In this geometry, the magnification is defined as

M = (z_{1} + z_{2}) / z_{1}

and

z_{e f f} = z_{2} / M

, where

z_{1}

is the source-object distance,

z_{2}

is the object-detector distance, and

z_{e f f}

is the effective propagation distance. The factor

1 / M^{2}

conserves energy.

Traditional PB-BPM works as an iterative process in the following way,

\begin{array}{l} U_{1} (x, y) = U_{0} (x, y) T_{1} (x, y) \\ U_{2} (x, y) = P_{Δ z} {U_{1} (x, y)} T_{2} (x, y) \\ U_{3} (x, y) = P_{Δ z} {U_{2} (x, y)} T_{3} (x, y) \\ ... \\ U_{k + 1} (x, y) = P_{Δ z} {U_{k} (x, y)} T_{k + 1} (x, y)} \end{array}

where

U_{k} (x, y)

is the wave at the

k^{t h}

plane,

T_{k} (x, y)

is the transmission function at the

k^{t h}

plane, and the operator

P_{Δ z}

is the propagator (for instance, angular spectrum of plane waves propagation) over the slice thickness

Δ z

[63]. The angular spectrum propagator, which is used here, can be expressed in operator notation as [30]

P_{Δ z} = F^{- 1} exp (i Δ z \sqrt{k^{2} - f_{x}^{2} - f_{y}^{2}}) F

where

k = 1 / λ

with

λ

being the wavelength,

f_{x}

and

f_{y}

are the spatial frequencies in x and y respectively,

F

is the Fourier transform operator, and

F^{- 1}

is the inverse Fourier transform operator.

An advantage of BPM is that it treats waves at any given plane as two dimensional, greatly simplifying the computation. In the limit where the amount of diffraction over the entire sample thickness is negligible, Eq. (5) reduces to simply multiplying all of the transmission functions together, and the result is equivalent to the projected-object approximation in Eqs. (1) and (2).

CB-BPM is implemented by applying the Fresnel scaling theorem (Eq. (4)) to the beam propagation method (Eq. (5)) at each slice. The method is described as follows.

3.3 Steps for applying CB-BPM

A primary advance of CB-BPM—one which allows it to account for multiple scattering and arbitrarily thick objects —is the application of the Fresnel scaling theorem within the object, rather than only at the exit surface wave. This is accomplished by adding an interpolation and distance re-scaling step at every plane in the object.

The procedure for going from plane $k$ to plane $k + 1$ is as follows. First, the effective slice thickness $Δ z$ is calculated at the $k + 1$ plane based on the change in magnification between adjacent slices (at the very first slice, the previous magnification $M_{k}$ is simply 1),

Δ z_{k + 1} = Δ z \frac{M_{k}}{M_{k + 1}} .

Then, the propagated wave at the $k + 1$ plane is found by applying the propagation operator to the wave at the $k^{t h}$ plane $U_{k} (x, y)$ ,

U_{P, k + 1} (x, y) = P_{Δ z_{_{k + 1}}} [U_{k} (x, y)]

where the subscript

P, k + 1

denotes that the wave has been propagated to the

k + 1

plane.

So far, $x$ and $y$ have been treated as general transverse coordinates. At this point, it becomes important to consider the scaling of the coordinates at each plane. The coordinates $x$ and $y$ will now be considered to be the specific coordinates at the detector plane.

The exit wave at the $k + 1$ plane is found by interpolating the wave it onto the appropriate grid based on the change in magnification. Here, the interpolation operation to transform the coordinates from plane $k$ to plane $k + 1$ for a wave that has already been propagated to the $k + 1$ plane is defined as

I_{k \to k + 1} [U_{P, k + 1} (\frac{x}{M_{k}}, \frac{y}{M_{k}})] \equiv U_{P, k + 1} (\frac{x}{M_{k}} \to \frac{y}{M_{k + 1}}, \frac{y}{M_{k}} \to \frac{y}{M_{k + 1}})

which takes the scaling from that of plane

k

to that of plane

k + 1

. While the wave has been propagated to plane

k + 1

, it needs to be scaled to plane

k + 1

as well. It is assumed that steps are taken to ensure that energy is conserved properly during the interpolation process. The final step is to multiply by the object’s (appropriately interpolated) transmission function at the

k + 1

plane. Putting it all together, the procedure can be described as follows: first, obtain the wave at the first plane by multiplying the first wave’s transmission function

T_{1} (x, y)

by the initial input wave

U_{0} (x, y)

, and interpolating this product to have the correct magnification,

U_{1} (\frac{x}{M_{1}}, \frac{y}{M_{1}}) = I_{1 \to M_{1}} [U_{0} (x, y) T_{1} (x, y)] .

Next, this wave is propagated to the second plane,

U_{P, 2} (\frac{x}{M_{1}}, \frac{y}{M_{1}}) = P_{Δ z_{_{2}}} [U_{1} (\frac{x}{M_{1}}, \frac{y}{M_{1}})]

Then, using the quantity calculated in Eq. (11), the wave is interpolated and multiplied by the appropriately scaled transmission function at the second plane to get the final wave at the second plane,

U_{2} (\frac{x}{M_{2}}, \frac{y}{M_{2}}) = I_{M_{2} \to M_{1}} [U_{P, 2} (\frac{x}{M_{1}}, \frac{y}{M_{1}})] T_{2} (\frac{x}{M_{2}}, \frac{y}{M_{2}}) .

The subscript on wave $U_{_{_{P, 2}}}$ denotes that the wave has been propagated to plane 2, but not yet scaled. Once the wave has been interpolated and multiplied by the transmission function, the operations at plane 2 are complete, and the wave $U_{2}$ has a subscript of simply 2.

The process is then repeated to go from the second to third plane and so on until the end of the object is reached. The general formula can then be expressed as,

U_{k + 1} (\frac{x}{M_{k + 1}}, \frac{y}{M_{k + 1}}) = I_{M_{k} \to M_{k + 1}} [U_{_{P, k}} (\frac{x}{M_{k}}, \frac{y}{M_{k}})] T_{k + 1} (\frac{x}{M_{k + 1}}, \frac{y}{M_{k + 1}}) .

Section 5 shows examples of simulations this method compared to other forward models.

4. Computational simulation methods

This section describes simulations carried out to test and validate CB-BPM, as well as illustrative simulations that show the important of multiple scattering, and the ability of CB-BPM to handle arbitrary objects. All simulations were performed in MATLAB.

Current non-projected object forward models do not incorporate cone-beam effects without exorbitant computational expense, and current implementations of the Mie solution for multiple scattering use parallel beam geometry. Therefore, it was necessary to perform accuracy comparisons using the parallel-beam version of CB-BPM, which reduces to the traditional PB-BPM for X-ray wavelengths.

The Mie scattering implementations used are those described by Schäfer et al. [53] because software is available in MATLAB, and the solutions were developed to be accurate in the near-field.

When possible, the physical parameters used in these simulations were matched to our experimental system under development [29,65]. However, the physical parameters had to be changed in certain cases in order to keep computation times with the Mie solution reasonable. The computation time for a Mie solution increases with the ratio of the object diameter to the wavelength. For the wavelength and length scales of interest for our imaging system design [29], this ratio is quite large ( $\approx 10^{7}$ or more), and it was therefore untenable to do calculations matching these parameters with the Mie solution. Therefore, the photon energy, object size, grid size, pixel size, and propagation distances had to be reduced in some cases. However, the refractive index was kept constant and assumed to be that of water at 30 keV [66] in order to keep the simulation as close as possible to the design parameters in the system under development [29]. All input waves were treated as monochromatic plane waves.

A guiding consideration in carrying out the parameter rescaling was the Fresnel number. The Fresnel number is defined as,

F_{N} = \frac{a^{2}}{λ z}

such that

a

is the aperture size or characteristic length scale of the object,

λ

is the wavelength, and

z

is the propagation distance [63]. The near-field is where

F_{N} > 1

, and the far-field is where

F_{N} ≪ 1

[63]. Many implementations of PB-XPCI that have large fields of view are performed in the near-field. In the near-field, there is little to moderate diffraction, while in the far-field, the diffraction is significant. However, it is important that the Fresnel number be small enough so that some diffraction effects can be observed. To maintain a Fresnel number that ensured the measured field would still be in the near-field (

F_{N} > 1

), the other simulation parameters, such as the pixel size and propagation distance, were scaled down as well. Furthermore, the Schäfer implementation for cylinders [53] only calculates the scattered field at a set distance based on the number of pixels used in the simulation and the physical size of the pixels.

While the Fresnel number is straight forward to compute for a case of a sphere or aperture, choosing a characteristic length scale is less simple for a large heterogeneous object. Therefore, in the case of the XCAT phantom the characteristic object size was considered to be about 16 pixels, although the phantom itself takes up nearly the entire FOV. When multiple object sizes or propagation distances were used, the Fresnel number is listed in Table 1 as a range.

Table 1. Parameters used in simulations

View Table

The reported relative differences between the Mie (or FWR) solution and the BPM are the mean differences in intensity along a lineout, namely, the mean of $(I_{2} - I_{1}) / I_{1}$ at each pixel along the lineout, where $I_{1}$ is chosen to be the Mie solution, or the FWR solution if the comparison is between BPM and FWR. This calculation was performed after the resulting intensities were normalized, and in the case of rotationally symmetric objects, after the intensity pattern was radially averaged to smooth aliasing artifacts that result from using a discrete propagation grid.

4.1 Single homogeneous sphere

First, a simulation was performed on a single homogeneous sphere and PB-BPM was compared to the Mie solution and the full-wave Rytov (FWR) solution. While this system falls within the single and weak scattering approximation, it is nonetheless important to verify PB-BPM on a simple system to compare computation times and verify that the paraxial approximation is appropriate.

The photon energy was reduced to 3 keV due to computation limitations of the Mie solution. The refractive index was set to be that of water at 30 keV. Other simulation parameters are listed in Table 1.

In order to fairly compare BPM to the FWR solution, it is important to use precisely the same scattering potential. Therefore, in this implementation, the scattering potential input to PB-BPM was set to be as close as possible to that of the FWR solution; a Hann-filtered sphere, in order to ensure that the two used the same object. This was achieved by running the FWR simulation with a propagation distance of 10⁻¹⁶ m, and using the complex exit wave as the input to BPM. This necessitated using the projected-object approximation, but in this particular case of a single weakly-scattering object, the projected-object approximation is valid.

A Hann-filtered sphere was used for the scattering potential here in order to suppress ringing artifacts caused by discretized voxels. This provides a fairer comparison to the Mie solution. Because the Mie solution uses an analytic expression, it does not suffer from the same discontinuity artifacts, such as spatial frequency aliasing or ringing at the boundary of the sphere, as the FWR and BPM approaches do.

The BPM intensities were radially averaged for clear presentation and to mitigate aliasing and discretization artifacts. In principle, the intensity pattern should be perfectly rotationally symmetric, since the objects and illumination are rotationally symmetric. However, because of propagation artifacts caused by utilizing a discrete propagation grid, the patterns exhibit minor distortions. Radially averaging produces a better representation of the overall pattern.

4.2 Single cylinder

Here, and for the multiple cylinders described in the following section (Sec. 4.3), the Schäfer implementation [53] was used to calculate the Mie solution for an infinite cylinder or cylinders [67,68]. These simulations are computationally intensive, particularly for high photon energies, so the photon energy had to be reduced. To maintain a Fresnel number that kept the solution in the near-field, the pixel and cylinder size had to be scaled as well. Furthermore, the geometry of the Schäfer implementation is such that the propagation distance is controlled by the pixel size and number of points, so the comparisons were performed for short distances.

For objects that are represented by discrete binary functions, each interface contains high frequency content that in turn causes ringing, making it difficult to directly compare to methods that use smooth or filtered functions to represent objects. To mitigate this artifact, the object was filtered at every axial slice with a Gaussian function (width = 4), and any smoothing effects beyond the physical boundary of the cylinder were masked out. Because the cylinder is infinite, the intensity pattern only varies in one dimension, and therefore the lineouts shown contain all the information about the intensity pattern.

4.3 Multiple cylinders

Next, multiple infinite cylinders separated by various distances were tested. This allows for a solution involving multiple scattering events to be calculated.

The validity of PB-BPM was tested against this case for two different cylinder geometries. One involves all cylinders being aligned with the only offset being along the beam propagation direction, such that the farther ones are ‘shadowed’ by the near ones. In parallel-beam geometry under the projected-object approximation, such a configuration would look like a single cylinder (with a different refractive index). The second case investigated was with the cylinders displaced along the beam propagation direction and perpendicular to beam propagation direction. The computational expense for this simulation is very high, so the photon energy had to be decreased to 0.3 keV. The refractive index was maintained at that of water for 30 keV, as discussed previously.

4.4 Multiple spheres

Next, the intensities were compared for scattering with and without the projected object approximation. The object was a system of two spheres aligned along the axis of propagation. Water and calcium spheres were each tested at 30 keV for a variety of separations, and for parallel- and cone-beam geometry. The mean differences were calculated comparing treating the object with and without the projected-object approximation. The resulting intensities were radially averaged as discussed in Sec. 4.1.

4.5 XCAT breast phantom

To demonstrate the use of CB-BPM on a clinically relevant phantom, the model was tested with the XCAT breast phantom [69–71], developed by researchers at Duke University, the University of Massachusetts, and the University of California Davis, to model coherent propagation through approximately 5 cm of breast tissue. The detailed phantom includes different refractive indices for adipose tissue, mammary tissue, skin, blood vessels, and ligaments [72]. The pixel and grid sizes were specified by the available XCAT phantom.

5. Results

5.1 Homogeneous sphere

Figure 2 shows good agreement between the Mie, FWR, and PB-BPM approaches for a single homogeneous sphere. In the example shown in Figs. 2(a)-2(c), the computation time for PB-BPM was 9 ms, for FWR it was 5 ms, and for the Mie solution it was 11.8 hours. Compared to the Mie solution, for a modest loss in accuracy (<1%), PB-BPM offers a computation speed increase of >10⁶.

Fig. 2 Results of different forward models for a homogeneous sphere of water. Example images from (a) beam propagation method, (b) full-wave Rytov, and (c) Mie solution for a 0.5 mm sphere of water propagated a distance of 2 m. The scale bar shown is shared among (a)-(c). (d) Comparisons of relative error between the three methods as a function of propagation distance. (e) Example lineout from the intensities shown in (a)-(c). All simulations used a 0.5 mm diameter sphere. The photon energy was reduced to 3 keV due to computation limitations of the Mie solution. The object for PB-BPM was set to be equal to the scattering potential used for the FWR method in order to ensure comparison of nearly-identical objects. The intensities were radially averaged.

Download Full Size | PDF

5.2 Scattering from a single cylinder

In the example shown in Fig. 3(b), the relative error is below $2 \times 10^{- 5}$ . The computation time for PB-BPM was 0.17 s, while for the Mie solution shown in Fig. 3(b) it was 8.3 hours. In the case of the largest cylinder diameter, 0.3 mm, the computation time of the Mie solution was 12.3 hours. PB-BPM is >10⁵ times faster, and conveniently, does not depend on any of the physical parameters in the simulation, but only on the grid size necessary for the calculation.

Fig. 3 Results of PB-BPM and Mie solution for an infinite cylinder. (a) Comparison of relative difference versus the cylinder diameter between PB-BPM and the Mie solution. (b) Lineout of scattering from a single cylinder. In this case, the diameter of the cylinder was 0.2 mm and the propagation distance was 0.256 mm, as set by the Schäfer implementation’s geometry. Here, the illuminating photon energy was 3 keV. The cylinder object treated with PB-BPM shown in (b) was Gaussian filtered to avoid ringing artifacts due to discretization.

Download Full Size | PDF

5.3 Multiple scattering from 3 cylinders

The examples shown in Figs. (4) and (5) illustrate that PB-BPM produces scattered intensities with mean differences of the order 10⁻⁵ compared to the Mie solution for the multiple scattering events investigated. The computation time for PB-BPM was 1.4 seconds, and the computation time for the Mie solution was 344 seconds.

Fig. 4 Simulation of three axially aligned cylinders. (a) Illustration of cylinder geometry and propagation direction. (b) Lineout of the scattered intensity. Here, the pixel size was 0.25μm, the cylinder diameters were 2 μm, and the photon energy was 0.3 keV. The cylinders were displaced −8 μm, 0 μm, and 12 μm from the center, respectively. Here, the mean difference was <10⁻⁵.

Download Full Size | PDF

Fig. 5 Simulation of three cylinders offset in two dimensions. (a) Schematic of the three cylinders and propagation direction. (b) Lineout of the intensity for the Mie solution compared to PB-BPM. In this simulation, the pixel size was 0.125 μm, the diameter was 2 μm, and the photon energy was 0.3 keV. The cylinders displaced in the (transverse, axial) coordinates from the center by (−2, −4), (0, 0), and (3, 6) μm, respectively. The mean difference between the intensities was <10⁻⁵.

Download Full Size | PDF

The multiple-cylinder solution is the most computationally demanding of all, and therefore it was necessary to use very low photon energies, small pixel sizes, and small diameter cylinders. Furthermore, the geometrical constraints of this implementation, namely the pre-set propagation distance, necessitated using an unphysically small detector pixel size and a long wavelength. Importantly, this change of parameters allows us to examine non-trivial diffraction of the cylinders, shown by the visible oscillations in in Figs. 4(b) and 5(b).

5.4 Simulations with cone-beam effects and multiple scattering within the object

The multiple scattering aspect of the BPM approach becomes particularly important for thick and strongly scattering objects. An illustrative case is that of two spheres separated along the axis of propagation. As expected, when the spheres are very close to each other, the projected-object approximation is reasonably good and the difference between treating the object with or without intra-sample scattering yields similar results. However, once the distance between the spheres increases—even modestly—the difference becomes significant.

Figure 6 demonstrates that two spheres—separated even a small amount—produces non-negligible multiple scattering events, and the intensity at the detector is significantly different if treated under the projected-object approximation. The difference is even more pronounced in cone-beam geometry compared to parallel beam geometry, as well as for more strongly scattering material, for instance calcium rather than water. Not accounting for multiple scattering, as is the case in Fig. 6(d), produces propagated intensity that simply looks like what would result from a single sphere of a higher refractive index. Meanwhile, Fig. 6(c) demonstrates a more complex intensity pattern, with interference caused by the second sphere.

Fig. 6 Comparison of intensity between the projected-object approximation and multiple scattering events. (a) Schematic of the simulation geometry. (b) Mean difference caused by the projected-object approximation for two axially separated spheres as a function of separation distance. Errors are shown for calcium and water spheres, and for cone beam and parallel beam propagation. (c) Example of two Ca spheres separated by 10 cm treated with CB-BPM and (d) treated under the projected-object approximation followed by the Fresnel scaling theorem. These simulations used 0.5 mm diameter spheres, and propagation distance of 1.5 m. In the cone beam simulation, the distance between the source and the object was 3.5 m. The illumination was 30 keV. The images were radially averaged to smooth propagation artifacts. The color scale is shared between (c) and (d).

Download Full Size | PDF

These errors do not arise as readily for incoherent imaging, such as traditional CT, where the Fresnel fringes caused by intra-object diffraction do not coherently interfere due to the relatively low spatial coherence of the beam, and therefore wave and refraction effects are not readily observed. It is important to note that beam spatial and temporal coherence plays a significant role in the visibility of diffraction effects in all types of imaging.

5.5 Simulations with XCAT breast phantom

Figure 7 shows the results of the XCAT breast phantom treated under different geometries and approximations, demonstrating that CB-BPM is capable of handling arbitrary objects. It is however important to ensure that the detector pixel size is small enough to properly sample Fresnel fringes from the highest frequency features in the object. The grainy appearances in Fig. 7(c), 7(f), and 7(i) are likely due to the fine structures such as ducts and blood vessels being undersampled. The relatively minor amount of diffraction, particularly in the cone beam case, suggests that in order to get more phase sensitivity for this type of object, a combination of smaller detector pixels and/or a longer propagation distance should be used. Still, the demonstration of CB-BPM with these XCAT phantoms shows a promising capability to sensitively model soft tissues in cone-beam geometry.

Fig. 7 Simulation with XCAT breast phantom, showing results from PB-BPM (top row), CB-BPM (middle row), and projected-object approximation followed by parallel beam propagation (bottom row). (a)-(c) PB-BPM for propagation distances of 0.1, 5, and 8 m, respectively. (d)-(f) show CB-BPM with a source-object (z₁) distance of 3.5 m, and object-detector distances (z₂) of 0.1, 1.5, and 4.5 m, respectively. (g)-(i) show the projected-object approximation, followed by parallel beam propagation of the wave over the remaining distance to the detector. Color scales are shared among all images. The photon energy used was 30 keV.

Download Full Size | PDF

6. Discussion

Many imaging applications such as conventional medical X-ray CT use multiple projections to tomographically reconstruct a 3D image [73]. However, for phase imaging where diffraction and cone-beam effects become important, simply performing filtered backprojection (FBP) may not be sufficient [74]. Indeed, a forward and/or backward model that accounts for diffraction is essential for coherent phase imaging. Multi-slice parallel beam reconstruction has been demonstrated in the optical regime [37], and here the forward model has been explored in silico for X-rays, along with novel cone-beam multi-slice propagation.

Other propagation methods, such as the wave propagation method (WPM), are accurate for non-paraxial waves but require increased computational complexity [33,37]. PB-BPM and CB-BPM utilize fast Fourier transforms, making them a good choice for iterative reconstruction, where the forward model has to be computed many times. Efficient proof-of-principle computational prototyping is also essential for imaging system design.

CB-BPM loses accuracy for high-angle, non-paraxial scattering, and can quickly become computationally intensive for large numbers of slices within the object. It is therefore less appropriate for longer-wavelength imaging of large, strongly scattering objects. Furthermore, it has been demonstrated here well within the Fresnel regime, and is not intended to propagate over distances that would be described by small Fresnel numbers that fall within the Fraunhofer regime. It is also important to note that CB-BPM does not take reflections or incoherent scattering events into account.

7. Conclusion

In this work, the use of PB-BPM has been validated for X-ray wavelengths, and a cone-beam multi-slice beam propagation method, called CB-BPM, has been described and demonstrated. This technique can model arbitrary sized objects and account for diffraction and magnification within that object. This allows for increased flexibility in a forward model that is capable of treating a variety of objects, only requiring that the paraxial approximation is valid and that the slices are sufficiently thin. Furthermore, the approximations made by this forward model are already within the scope of assumptions made for PB-XPCI when phase retrieval is achieved with the transport-of-intensity equation.

The computation time of CB-BPM depends only on the size of the grids, for instance the number of pixels and how many slices the object is divided into. Future work of interest will be validating this model experimentally and incorporating it into a reconstruction framework, and combining it with machine learning approaches. Interesting extensions also include implementations with asymmetric grid or pixel sizes or in regimes with partial spatial coherence.

Funding

Team Science Award from Mayo Clinic and Arizona State University, U.S. Department of Defense (W81XWH-13-2-0067, CDMRP Air Force Contract FA8650-17-C-9113), and USAMRAA (W81XWH-15-C-0052, W81XWH-17-C-0068).

Acknowledgment

We would like to thank Dr. Joseph Lo from Duke University for providing the XCAT phantom.

Disclosures

Dr. McCollough receives grant funding from Siemens Healthcare for work unrelated to this project. No other authors report a potential conflict of interest.

References

1. A. Bravin, P. Coan, and P. Suortti, “X-ray phase-contrast imaging: from pre-clinical applications towards clinics,” Phys. Med. Biol. 58(1), R1–R35 (2013). [CrossRef] [PubMed]

2. T. E. Gureyev, S. C. Mayo, D. E. Myers, Y. I. Nesterets, D. M. Paganin, A. Pogany, A. W. Stevenson, and S. W. Wilkins, “Refracting Röntgen’s rays: propagation-based x-ray phase contrast for biomedical imaging,” J. Appl. Phys. 105(10), 102005 (2009). [CrossRef]

3. M. Endrizzi, “X-ray phase-contrast imaging,” Nucl. Instrum. Methods Phys. Res. A 878, 88–98 (2018). [CrossRef]

4. T. E. Gureyev and K. A. Nugent, “Phase retrieval with the transport-of-intensity equation. II. Orthogonal series solution for nonuniform illumination,” J. Opt. Soc. Am. A 13(8), 1670–1682 (1996). [CrossRef]

5. D. Paganin, S. C. Mayo, T. E. Gureyev, P. R. Miller, and S. W. Wilkins, “Simultaneous phase and amplitude extraction from a single defocused image of a homogeneous object,” J. Microsc. 206(Pt 1), 33–40 (2002). [CrossRef] [PubMed]

6. J. P. Guigay, M. Langer, R. Boistel, and P. Cloetens, “Mixed transfer function and transport of intensity approach for phase retrieval in the Fresnel region,” Opt. Lett. 32(12), 1617–1619 (2007). [CrossRef] [PubMed]

7. A. Horng, E. Brun, A. Mittone, S. Gasilov, L. Weber, T. Geith, S. Adam-Neumair, S. D. Auweter, A. Bravin, M. F. Reiser, and P. Coan, “Cartilage and soft tissue imaging using X-rays: propagation-based phase-contrast computed tomography of the human knee in comparison with clinical imaging techniques and histology,” Invest. Radiol. 49(9), 627–634 (2014). [CrossRef] [PubMed]

8. T. Geith, E. Brun, A. Mittone, S. Gasilov, L. Weber, S. Adam-Neumair, A. Bravin, M. Reiser, P. Coan, and A. Horng, “Quantitative Assessment of Degenerative Cartilage and Subchondral Bony Lesions in a Preserved Cadaveric Knee: Propagation-Based Phase-Contrast CT Versus Conventional MRI and CT,” AJR Am. J. Roentgenol. 210(6), 1317–1322 (2018). [CrossRef] [PubMed]

9. S. T. Taba, T. E. Gureyev, M. Alakhras, S. Lewis, D. Lockie, and P. C. Brennan, “X-ray phase-contrast technology in breast imaging: principles, options, and clinical application,” AJR Am. J. Roentgenol. 211(1), 133–145 (2018). [CrossRef] [PubMed]

10. L. C. P. Croton, K. S. Morgan, D. M. Paganin, L. T. Kerr, M. J. Wallace, K. J. Crossley, S. L. Miller, N. Yagi, K. Uesugi, S. B. Hooper, and M. J. Kitchen, “In situ phase contrast X-ray brain CT,” Sci. Rep. 8(1), 11412 (2018). [CrossRef] [PubMed]

11. F. Pfeiffer, O. Bunk, C. David, M. Bech, G. Le Duc, A. Bravin, and P. Cloetens, “High-resolution brain tumor visualization using three-dimensional x-ray phase contrast tomography,” Phys. Med. Biol. 52(23), 6923–6930 (2007). [CrossRef] [PubMed]

12. G. Shinohara, K. Morita, M. Hoshino, Y. Ko, T. Tsukube, Y. Kaneko, H. Morishita, Y. Oshima, H. Matsuhisa, R. Iwaki, M. Takahashi, T. Matsuyama, K. Hashimoto, and N. Yagi, “Three dimensional visualization of human cardiac conduction tissue in whole heart specimens by high-resolution phase-contrast CT imaging using synchrotron radiation,” World J. Pediatr. Congenit. Heart Surg. 7(6), 700–705 (2016). [CrossRef] [PubMed]

13. P. Cloetens, W. Ludwig, J. Baruchel, D. Van Dyck, J. Van Landuyt, J. P. Guigay, and M. Schlenker, “Holotomography: quantitative phase tomography with micrometer resolution using hard synchrotron radiation x rays,” Appl. Phys. Lett. 75(19), 2912–2914 (1999). [CrossRef]

14. S. C. Mayo, P. R. Miller, S. W. Wilkins, T. J. Davis, D. Gao, T. E. Gureyev, D. Paganin, D. J. Parry, A. Pogany, and A. W. Stevenson, “Quantitative X-ray projection microscopy: phase-contrast and multi-spectral imaging,” J. Microsc. 207(Pt 2), 79–96 (2002). [CrossRef] [PubMed]

15. R. Mokso, P. Cloetens, E. Maire, W. Ludwig, and J. Y. Buffière, “Nanoscale zoom tomography with hard x rays using Kirkpatrick-Baez optics,” Appl. Phys. Lett. 90(14), 144104 (2007). [CrossRef]

16. E. Maire and P. J. Withers, “Quantitative X-ray tomography,” Int. Mater. Rev. 59(1), 1–43 (2014). [CrossRef]

17. A. Beerlink, S. Thutupalli, M. Mell, M. Bartels, P. Cloetens, S. Herminghaus, and T. Salditt, “X-ray propagation imaging of a lipid bilayer in solution,” Soft Matter 8(17), 4595–4601 (2012). [CrossRef]

18. M. Krenkel, M. Töpperwien, C. Dullin, F. Alves, and T. Salditt, “Propagation-based phase-contrast tomography for high-resolution lung imaging with laboratory sources,” AIP Adv. 6(3), 035007 (2016). [CrossRef]

19. M. Langer, A. Pacureanu, H. Suhonen, Q. Grimal, P. Cloetens, and F. Peyrin, “X-ray phase nanotomography resolves the 3D human bone ultrastructure,” PLoS One 7(8), e35691 (2012). [CrossRef] [PubMed]

20. M. Langer and F. Peyrin, “3D X-ray ultra-microscopy of bone tissue,” Osteoporos. Int. 27(2), 441–455 (2016). [CrossRef] [PubMed]

21. A. Khimchenko, C. Bikis, A. Pacureanu, S. E. Hieber, P. Thalmann, H. Deyhle, G. Schweighauser, J. Hench, S. Frank, M. Müller-Gerbl, G. Schulz, P. Cloetens, and B. Müller, “Hard x-ray nanoholotomography: large-scale, label-free, 3D neuroimaging beyond optical limit,” Adv. Sci. (Weinh.) 5(6), 1700694 (2018). [CrossRef] [PubMed]

22. C. Olendrowitz, M. Bartels, M. Krenkel, A. Beerlink, R. Mokso, M. Sprung, and T. Salditt, “Phase-contrast x-ray imaging and tomography of the nematode Caenorhabditis elegans,” Phys. Med. Biol. 57(16), 5309–5323 (2012). [CrossRef] [PubMed]

23. M. Born and E. Wolf, Principles of Optics, 7th ed. (Cambridge University, 1999).

24. M. I. Mishchenko, L. D. Travis, and D. W. Mackowski, “T-matrix computations of light scattering nonspherical particles: a review,” J. Quant. Spectrosc. Radiat. Transf. 55(5), 535–575 (1996). [CrossRef]

25. G. Gouesbet, J. A. Lock, and G. Grehan, “Generalized Lorenz–Mie theories and description of electromagnetic arbitrary shaped beams: localized approximations and localized beam models, a review,” J. Quant. Spectrosc. Radiat. Transf. 112(1), 1–27 (2011). [CrossRef]

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

27. Y. Sung, C. J. R. Sheppard, G. Barbastathis, M. Ando, and R. Gupta, “Full-wave approach for x-ray phase imaging,” Opt. Express 21(15), 17547–17557 (2013). [CrossRef] [PubMed]

28. Y. Sung, W. P. Segars, A. Pan, M. Ando, C. J. R. Sheppard, and R. Gupta, “Realistic wave-optics simulation of X-ray phase-contrast imaging at a human scale,” Sci. Rep. 5(1), 12011 (2015). [CrossRef] [PubMed]

29. Y. Sung, R. Gupta, B. Nelson, S. Leng, C. H. McCollough, and W. S. Graves, “Phase-contrast imaging with a compact x-ray light source: system design,” J. Med. Imaging (Bellingham) 4(4), 043503 (2017). [CrossRef] [PubMed]

30. D. M. Paganin, Coherent X-Ray Optics (Oxford University, 2006).

31. J. A. Fleck Jr., J. R. Morris, and M. D. Feit, “Time-dependent propagation of high energy laser beams through the atmosphere,” Appl. Phys. (Berl.) 10(2), 129–160 (1976). [CrossRef]

32. J. M. Cowley and A. F. Moodie, “The scattering of electrons by atoms and crystals. I. A new theoretical approach,” Acta Crystallogr. 10(10), 609–619 (1957). [CrossRef]

33. U. S. Kamilov, I. N. Papadopoulos, M. H. Shoreh, A. Goy, C. Vonesch, M. Unser, and D. Psaltis, “Optical tomographic image reconstruction based on beam propagation and sparse regularization,” IEEE Trans. Comput. Imaging 2(1), 59–70 (2016). [CrossRef]

34. U. S. Kamilov, I. N. Papadopoulos, M. H. Shoreh, A. Goy, C. Vonesch, M. Unser, and D. Psaltis, “A learning approach to optical tomography,” Optica 2(6), 517–522 (2015). [CrossRef]

35. W. G. Tam, “Split-step Fourier-transform analysis for laser-pulse propagation in particulate media,” J. Opt. Soc. Am. A 72(10), 1311–1316 (1982). [CrossRef]

36. A. M. Maiden, M. J. Humphry, and J. M. Rodenburg, “Ptychographic transmission microscopy in three dimensions using a multi-slice approach,” J. Opt. Soc. Am. A 29(8), 1606–1614 (2012). [CrossRef] [PubMed]

37. X. Ma, W. Xiao, and F. Pan, “Optical tomographic reconstruction based on multi-slice wave propagation method,” Opt. Express 25(19), 22595–22607 (2017). [CrossRef] [PubMed]

38. T. M. Godden, R. Suman, M. J. Humphry, J. M. Rodenburg, and A. M. Maiden, “Ptychographic microscope for three-dimensional imaging,” Opt. Express 22(10), 12513–12523 (2014). [CrossRef] [PubMed]

39. K. Li, M. Wojcik, and C. Jacobsen, “Multislice does it all-calculating the performance of nanofocusing X-ray optics,” Opt. Express 25(3), 1831–1846 (2017). [CrossRef] [PubMed]

40. J. Lim, A. Goy, M. H. Shoreh, M. Unser, and D. Psaltis, “Learning tomography assessed using Mie theory,” Phys. Rev. Appl. 9(3), 34027 (2018). [CrossRef]

41. E. Maire and P. J. Withers, “Quantitative x-ray tomography,” Int. Mater. Rev. 59(1), 1–43 (2014). [CrossRef]

42. M. Krenkel, A. Markus, M. Bartels, C. Dullin, F. Alves, and T. Salditt, “Phase-contrast zoom tomography reveals precise locations of macrophages in mouse lungs,” Sci. Rep. 5(1), 9973 (2015). [CrossRef] [PubMed]

43. M. Guizar-Sicairos, J. J. Boon, K. Mader, A. Diaz, A. Menzel, and O. Bunk, “Quantitative interior x-ray nanotomography by a hybrid imaging technique,” Optica 2(3), 259–266 (2015). [CrossRef]

44. M. Bartels, M. Krenkel, P. Cloetens, W. Möbius, and T. Salditt, “Myelinated mouse nerves studied by X-ray phase contrast zoom tomography,” J. Struct. Biol. 192(3), 561–568 (2015). [CrossRef] [PubMed]

45. S. Steckmann, M. Knaup, and M. Kachelriess, “Algorithm for hyperfast cone-beam spiral backprojection,” Comput. Methods Programs Biomed. 98(3), 253–260 (2010). [CrossRef] [PubMed]

46. M. Beister, D. Kolditz, and W. A. Kalender, “Iterative reconstruction methods in X-ray CT,” Phys. Med. 28(2), 94–108 (2012). [CrossRef] [PubMed]

47. M. Nilchian, Z. Wang, T. Thuering, M. Unser, and M. Stampanoni, “Spline based iterative phase retrieval algorithm for X-ray differential phase contrast radiography,” Opt. Express 23(8), 10631–10642 (2015). [CrossRef] [PubMed]

48. J. A. Lock and G. Gouesbet, “Generalized Lorenz–Mie theory and applications,” J. Quant. Spectrosc. Radiat. Transf. 110(11), 800–807 (2009). [CrossRef]

49. N. J. Moore and M. A. Alonso, “Closed form formula for Mie scattering of nonparaxial analogues of Gaussian beams,” Opt. Express 16(8), 5926–5933 (2008). [CrossRef] [PubMed]

50. P. C. Waterman, “Matrix formulation of electromagnetic scattering,” Proc. IEEE 53(8), 805–812 (1965). [CrossRef]

51. M. I. Mishchenko, “Light scattering by randomly oriented axially symmetric particles,” J. Opt. Soc. Am. A 8(6), 871–882 (1991). [CrossRef]

52. F. M. Kahnert, “Numerical methods in electromagnetic scattering theory,” J. Quant. Spectrosc. Radiat. Transf. 79–80(1), 775–824 (2003). [CrossRef]

53. J. Schäfer, S. C. Lee, and A. Kienle, “Calculation of the near fields for the scattering of electromagnetic waves by multiple infinite cylinders at perpendicular incidence,” J. Quant. Spectrosc. Radiat. Transf. 113(16), 2113–2123 (2012). [CrossRef]

54. S.-C. Lee, “Dependent scattering of an obliquely incident plane wave by a collection of parallel cylinders,” J. Appl. Phys. 68(10), 4952–4957 (1990). [CrossRef]

55. G. A. Tsihrintzis and A. J. Devaney, “Higher order (nonlinear) diffraction tomography: inversion of the Rytov series,” IEEE Trans. Inf. Theory 46(5), 1748–1761 (2000). [CrossRef]

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

57. E. Soubies, T.-A. Pham, and M. Unser, “Efficient inversion of multiple-scattering model for optical diffraction tomography,” Opt. Express 25(18), 21786–21800 (2017). [CrossRef] [PubMed]

58. M. Lee, S. Shin, and Y. Park, “Reconstructions of refractive index tomograms via a discrete algebraic reconstruction technique,” Opt. Express 25(22), 27415–27430 (2017). [CrossRef] [PubMed]

59. C. T. Putkunz, M. A. Pfeifer, A. G. Peele, G. J. Williams, H. M. Quiney, B. Abbey, K. A. Nugent, and I. McNulty, “Fresnel coherent diffraction tomography,” Opt. Express 18(11), 11746–11753 (2010). [CrossRef] [PubMed]

60. M. Chen, L. Tian, and L. Waller, “3D differential phase contrast microscopy,” Biomed. Opt. Express 7(10), 3940–3950 (2016). [CrossRef] [PubMed]

61. M. R. Teague, “Deterministic phase retrieval: a Green’s function solution,” J. Opt. Soc. Am. 73(11), 1434–1441 (1983). [CrossRef]

62. 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]

63. J. W. Goodman, Introduction to Fourier Optics, 3rd ed. (Roberts & Company Publishers, 2005).

64. J. Miao, T. Ishikawa, I. K. Robinson, and M. M. Murnane, “Beyond crystallography: diffractive imaging using coherent x-ray light sources,” Science 348(6234), 530–535 (2015). [CrossRef] [PubMed]

65. W. S. Graves, J. Bessuille, P. Brown, S. Carbajo, V. Dolgashev, K.-H. Hong, E. Ihloff, B. Khaykovich, H. Lin, K. Murari, E. A. Nanni, G. Resta, S. Tantawi, L. E. Zapata, F. X. Kärtner, and D. E. Moncton, “Compact x-ray source based on burst-mode inverse Compton scattering at 100 kHz,” Phys. Rev. Spec. Top. Accel. Beams 17(12), 120701 (2014). [CrossRef]

66. B. L. Henke, E. M. Gullikson, and J. C. Davis, “X-ray interactions: photoabsorption, scattering, transmission, and reflection at E=50-30000 eV, Z=1-92,” At. Data Nucl. Data Tables 54(2), 181–342 (1993). [CrossRef]

67. K. S. Morgan, K. K. W. Siu, and D. M. Paganin, “The projection approximation and edge contrast for x-ray propagation-based phase contrast imaging of a cylindrical edge,” Opt. Express 18(10), 9865–9878 (2010). [CrossRef] [PubMed]

68. K. S. Morgan, K. K. W. Siu, and D. M. Paganin, “The projection approximation versus an exact solution for x-ray phase contrast imaging, with a plane wave scattered by a dielectric cylinder,” Opt. Commun. 283(23), 4601–4608 (2010). [CrossRef]

69. C. M. Li, W. P. Segars, G. D. Tourassi, J. M. Boone, and J. T. Dobbins 3rd, “Methodology for generating a 3D computerized breast phantom from empirical data,” Med. Phys. 36(7), 3122–3131 (2009). [CrossRef] [PubMed]

70. D. W. Erickson, J. R. Wells, G. M. Sturgeon, E. Samei, J. T. Dobbins 3rd, W. P. Segars, and J. Y. Lo, “Population of 224 realistic human subject-based computational breast phantoms,” Med. Phys. 43(1), 23–32 (2016). [CrossRef] [PubMed]

71. G. M. Sturgeon, N. Kiarashi, J. Y. Lo, E. Samei, and W. P. Segars, “Finite-element modeling of compression and gravity on a population of breast phantoms for multimodality imaging simulation,” Med. Phys. 43(5), 2207–2217 (2016). [CrossRef] [PubMed]

72. D. R. White, R. V. Griffith, and I. J. Wilson, “ICRU Report 46: photon, electron, proton, and neutron interaction data for body tissues,” J. ICRUos24(1), (1992).

73. J. Hsieh, Computed Tomography Principles, Design, Artifacts, and Recent Advances, 3rd ed. (SPIE Press, 2015).

74. G. Gbur and E. Wolf, “Relation between computed tomography and diffraction tomography,” J. Opt. Soc. Am. A 18(9), 2132–2137 (2001). [CrossRef] [PubMed]

Object type	Photon energy (keV)	Pixel size (μm)	Grid size	Propagation distance	Refractive index material and photon energy	Object size	Fresnel number
Homogeneous sphere	3	16	64x64	0.5-2 m	Water at 30 keV	0.5 mm	302 - 1.2x10³
Single cylinder	3	8	64x64	0.256 mm	Water at 30 keV	0.05-0.3 mm	2.4x10^{4 -} 8.5x10⁵
Three aligned cylinders	0.3	0.25	128x128	16 μm	Water at 30 keV	2 μm	60
Three offset cylinders	0.3	0.125	128x128	8 μm	Water at 30 keV	2 μm	121
Two aligned spheres	30	3.45	256x256	1.5 m	Water, Ca at 30 keV	0.5 mm	4.0x10³
XCAT phantom	30	78.5	1000x1000	0.1-4.5 m	Various tissues at 30 keV	1.25 mm	8.4x10³ - 3.8x10⁵

Forward model for propagation-based x-ray phase contrast imaging in parallel- and cone-beam geometry

Abstract

1. Introduction

2. Forward models used for propagation-based X-ray phase-contrast imaging

2.1 Mie solution

2.2 First-order solutions

2.3 Beam propagation method

3. Cone-beam beam propagation method

3.1 The paraxial approximation

3.2 Fresnel scaling theorem and PB-BPM

3.3 Steps for applying CB-BPM

4. Computational simulation methods

4.1 Single homogeneous sphere

4.2 Single cylinder

4.3 Multiple cylinders

4.4 Multiple spheres

4.5 XCAT breast phantom

5. Results

5.1 Homogeneous sphere

5.2 Scattering from a single cylinder

5.3 Multiple scattering from 3 cylinders

5.4 Simulations with cone-beam effects and multiple scattering within the object

5.5 Simulations with XCAT breast phantom

6. Discussion

7. Conclusion

Funding

Acknowledgment

Disclosures

References

Cited By

Figures (7)

Tables (1)

Equations (14)

Optics Express