Rytov approximation for x-ray phase imaging

Yongjin Sung; George Barbastathis

doi:10.1364/OE.21.002674

1. Introduction

X-ray phase imaging (XPI) refers to the techniques that measure the amount of wavefront distortion induced by the object under investigation [1]. Interest in XPI keeps increasing since it can provide information on the electron density of materials, which was not available in conventional radiography images. XPI techniques can be largely divided into three categories: analyzer-based imaging, grating interferometry, and propagation-based methods [2]. In analyzer-based imaging [3–5] and grating interferometry [6–8] the angle of refracted rays or the gradient of the distorted wavefront are measured with an analyzer crystal and a pair of gratings, respectively. In propagation-based methods [9–11], the detector records one or more intensity images downstream the optical axis; the variation along the axis is then connected to the wavefront profile of the original beam using an appropriate propagation model. In the simplest approach for XPI simulation, the 2-D transmittance function of a sample is obtained using line integrals of the sample’s 3-D refractive index map along rays, and the intensity image at the detector is calculated by convolution of the transmittance function with the Fresnel kernel [1]. However, the projection approach is valid only with small objects, where one can assume small angle scattering; thus, a more rigorous approach including the diffraction within the sample is required for large objects.

As with other regimes of electromagnetic radiation, the intensity and wavefront profile of an x-ray beam after an object can be described by the wave equation [1]. Since the rigorous solution to the scalar wave equation can be quite involved, approximations such as those attributed to Born [12] and Rytov [13] have been adopted to provide approximate but explicit solutions. The validity of the approximations depends on a number of parameters, such as the geometry of the object, especially its thickness; and the refractive index contrast between the object and the surrounding medium or between regions of the object itself that have variable index. Regarding the physical dimension of the object, the validity condition for the Rytov approximation has caused considerable controversy [14, 15]. Applying the method of renormalization group (RG) to homogeneous slabs, where explicit solutions for the scattered field are available, Kirkinis [16] recently confirmed that the Rytov approximation is more robust than the Born approximation for homogeneous slabs as the slab thickness increases. However, similar analyses have not been carried out, to our knowledge, for more general cases. Instead, the validity of the approximation is based on implicit arguments without comparison to explicit solutions for the scattered field [17].

In this paper, we investigate the validity and accuracy of the Rytov approximation in x-ray phase imaging of a homogeneous sphere with an arbitrary complex refractive index and radius, a case in which explicit solutions for the scattered field are again available. We show that the validity of the Rytov approximation can be determined by a single parameter V defined in terms of the complex refractive index of the sphere and the Fresnel number. The parameter V relieves us of the need to calculate the exact Mie solution for very large spheres, when the number of terms required in the Mie series expansion becomes prohibitive. We show that the error can be calculated from the value of V, without having to explicitly calculate the Mie series.

By comparing with the exact Mie solution (either explicitly for small spheres or using the principle of similarity for large ones), we explain how to predict the accuracy of the Rytov approximation in specific cases of complex refractive index and radius of the sphere. We also present the maximum size of a water sphere for which the Rytov approximation is valid within 1% accuracy. This may provide a practical criterion to determine the validity of the Rytov approximation in x-ray phase imaging of large objects such as luggage or a human patient. We consider the cases of intensity and phase measurement separately.

2. Scattered field calculation under the first Rytov approximation

Let a plane wave or a collimated electromagnetic beam be incident on a homogeneous sphere. Figure 1 shows the schematic diagram and definition of the spherical coordinates.

Fig. 1 Schematic diagram of the geometry.

Download Full Size | PDF

The complex amplitude of the field after its interaction with the object can be described by the scalar wave equation [18].

(\nabla^{2} + k {(\vec{r})}^{2}) Ψ (\vec{r}) = 0,

where

k (\vec{r}) = k_{0} n (\vec{r})

,

k_{0} = 2 π / λ_{0}

,

λ_{0}

is the wavelength of the incident beam, and

n (\vec{r})

is the complex refractive index of the object.

Suppose the complex amplitude of the field can be written in the following form [13, 19]:

Ψ (\vec{r}) = Ψ_{0} (\vec{r}) \exp {ϕ_{s} (\vec{r})},

where

Ψ_{0} (\vec{r})

is the complex amplitude of the incident field, and

ϕ_{s} (\vec{r})

is the scattered complex phase. The real part of

ϕ_{s}

represents distortion of the wavefront, while its imaginary part represents attenuation of the incident light.

Substituting the expression for $Ψ$ , the first-order perturbation solution of Eq. (1) can be expressed as an integral equation

(Ψ_{0} ϕ_{s}) (\vec{r}) = \int G (\vec{r} - {\vec{r}}^{'}) Ψ_{0} ({\vec{r}}^{'}) [{(\nabla ϕ_{s} ({\vec{r}}^{'}))}^{2} - Q ({\vec{r}}^{'})] d^{3} r^{'},

where

G (\vec{r}) = \exp (i k | \vec{r} |) / (4 π | \vec{r} |)

,

Q (\vec{r}) = k_{0}^{2} (1 - n {(\vec{r})}^{2})

is the scattering potential of the object. Equation (3) is implicit in the sense that the unknown function

ϕ_{s} (\vec{r})

is included within the integral on the right hand side. The first Rytov approximation greatly simplifies the solution by assuming [13, 19]

| {(\nabla ϕ_{s})}^{2} | < < | Q | .

Then, Eq. (3) can be simplified to

ϕ_{s} (\vec{r}) = - \frac{1}{Ψ_{0} (\vec{r})} \int G (\vec{r} - {\vec{r}}^{'}) Ψ_{0} ({\vec{r}}^{'}) Q ({\vec{r}}^{'}) d^{3} r^{'} .

After some manipulation,

φ_{s}

and Q can be simply related in the spatial frequency space as follows [13, 20]:

{\tilde{Φ}}_{s} (U, V; z) = \frac{1}{i 4 π w} \exp (i 2 π W z) \tilde{Q} (U, V, W),

where

{\tilde{Φ}}_{s}

is the 2-D Fourier transform of

ϕ_{s}

, and

\tilde{Q}

is the 3-D Fourier transform of Q. The variables w and W are defined as

w = \sqrt{{(1 / λ)}^{2} - U^{2} - V^{2}}

and

W = w - 1 / λ

, respectively. Under the projection and paraxial approximations, Eq. (6) can be written as

{\tilde{Φ}}_{s} (U, V; z) = \frac{λ}{i 4 π} \exp [- i π λ z (U^{2} + V^{2})] \tilde{Q} (U, V, 0) .

3. Explicit validity condition of the first Rytov approximation for a homogeneous sphere

For a homogeneous sphere with radius R and complex refractive index n, $\tilde{Q}$ can be written as [21]

\tilde{Q} (\vec{k}) = Q_{0} \tilde{S} (\vec{k}),

\tilde{S} (\vec{k}) = {\begin{matrix} 4 π R^{3} / 3 & for | \vec{k} | = 0 \\ \frac{4 π R}{{(2 π | \vec{k} |)}^{2}} [\frac{\sin (2 π | \vec{k} | R)}{2 π | \vec{k} | R} - \cos (2 π | \vec{k} | R)] & otherwise \end{matrix},

where

Q_{0} = k_{0}^{2} (1 - n^{2})

.

Taking the inverse Fourier transform of Eq. (9) [21], and substituting in (7) and then (2), we find

ϕ_{s} (\vec{r}) = \frac{λ R}{i 2 π} Q_{0} \int_{0}^{\infty} \frac{\exp (- i π λ z ρ^{2})}{ρ} J_{0} (2 π ρ r) [\frac{\sin (2 π ρ R)}{2 π ρ R} - \cos (2 π ρ R)] d ρ,

where

r = \sqrt{x^{2} + y^{2}}

, and

J_{α} (x)

is the Bessel function of the first kind and α is the order of the Bessel function. From this, we can explicitly calculate the neglected term from the Rytov approximation as

{(\nabla ϕ_{s})}^{2} = - k_{0}^{2} {(1 - n^{2})}^{2} {\int_{0}^{\infty} \exp (- i t^{2} / F) J_{1} (\frac{r}{R} t) [\frac{\sin (t)}{t} - \cos (t)] d t}^{2},

where

F \equiv 4 π R^{2} / (λ z)

is the Fresnel number.

Substituting Eq. (11) for ${(\nabla ϕ_{s})}^{2}$ in Eq. (4), the validity condition of the first Rytov approximation may be expressed for a homogeneous sphere as

V \equiv | 1 - n^{2} | A (F) < < 1,

where the function A(F) is defined as follows:

A (F) \equiv \max_{β \in [0, \infty)} | {\int_{0}^{\infty} \exp (- i t^{2} / F) J_{1} (β t) [\frac{\sin (t)}{t} - \cos (t)] d t}^{2} | .

We obtain A(F) by numerically integrating Eq. (13) for different values of β and searching for the maximum. From the projection approximation and the known profile of the sphere, we should expect

| \nabla ϕ_{s} |

to have a maximum near the edge of the sphere; for that reason, we may well search only the region near β ; 1. Figure 2 shows A(F) calculated for seven different values of F, indicating that for large enough values of F, A(F) is well approximated by a linear fit (R² = 0.9993), as

A (F) \approx 0.632 F^{0.493} \approx 0.552 \sqrt{F} (F > 10^{7}) .

We found that the square-root approximation works very well for F ≳ 10⁷ (95% confidence interval, [0.5 0.5001]). For smaller values of F, the accuracy of (14) as an approximation to (13) worsens.

Fig. 2 F vs. A(F) plotted in the logarithmic scale of base 10. The data points are fitted to log₁₀(y) = a + blog₁₀(x) (R² = 0.9999), for which a = −0.199 (95% confidence interval, [-0.238 −0.160]), and b = 0.493 (95% confidence interval, [0.487 0.499]).

Download Full Size | PDF

Note that the scattered complex phase in Eq. (4) cannot be known a priori; therefore, the validity of the Rytov approximation cannot be checked until the scattered field is calculated. On the other hand, the parameter V in Eq. (12) is a function of known parameters, and thus it can be used to predict the validity of the approximation without the need for a full calculation of the scattered field. More importantly, two spheres of a different size may provide a same V value when the other parameters such as the wavelength λ and the distance z are properly chosen; the two cases can be simulated with the same accuracy using the Rytov approximation. This is analogous to the principle of similarity in fluid mechanics, where the Reynolds number similarly serves as a scaling parameter. As expected, the parameter V depends on the size of the object being imaged, although in the case of a sphere the size also determines the radius of curvature “imparted” on the scattered wavefront.

4. Error estimation of the Rytov approximation: comparison with Mie solution

Next, we consider the validity of the Rytov approximation in x-ray imaging of large objects such as luggage or a human patient, for which the ratio R⁄λ is on the order of 10⁹ and, hence, so is the F parameter. Mie theory provides an exact formula for the complex vector field scattered from a homogeneous sphere [18, 22]; therefore, it can serve as a gold standard with which the accuracy of the Rytov approximation can be checked. Since the Mie solution is an infinite series of vector spherical harmonics, in practice higher order terms are truncated to provide a solution within a given accuracy [23]. The number of terms required to achieve certain accuracy increases linearly with R⁄λ, where R is the radius of the sphere and λ is the wavelength of the incident beam. However, in case the ratio R⁄λ is extremely large as in our case, alternative approaches need be sought, e.g., adopting an asymptotic formula for the Mie series, transforming the series into an equivalent but a more rapidly converging sum, etc [24]. Here, we suggest an alternative method, scaling down the sphere with V in Eq. (12) as a scaling parameter and predict the accuracy of the Rytov approximation based on the model system.

Suppose a sphere of water with radius 5 cm is illuminated with an incident beam of 60 keV, and the scattered field from the sphere is recorded at 1 m distance from the center of the sphere. The complex refractive index of water at 60 keV can be expressed as n = 1-δ + iβ, where δ = 6.844 × 10⁻⁸ and β = 3.379 × 10⁻¹¹ [25]. The value of V in this case is 0.0027. Note that the full Mie series cannot be directly evaluated because of the large R⁄λ ratio (2.42 × 10⁹). Instead, we consider a model system with the same V value, but with R⁄λ = 10⁴. We put the measurement plane at z = 4πR to guarantee the scalar field assumption in Eq. (1) holds. For these choices of parameters, the value of F is 10⁴. Now, we find the wavelength of the incident beam to guarantee the same V value. Since most materials are highly dispersive in x-rays, we need to solve $| 1 - n {(λ)}^{2} | = V / A (F)$ . For water, we obtain λ = 0.256 nm, and the other parameters are accordingly determined as R = 2.56 μm, z = 32.2 μm. In a nutshell, the accuracy of the Rytov approximation for the model system (λ = 0.256 nm, R = 2.56 μm, and z = 32.2 μm) is the same as that for the original system (λ = 0.0207 nm, R = 5 cm, and z = 1 m). Since the radius of sphere in the model system is small enough, we can check the validity of the Rytov approximation using the Mie solution.

The intensity of the light recorded with a detector is given by the component of the Poynting vector on the direction perpendicular to the detector. When the detector plane is perpendicular to the optical axis [26],

I = \frac{1}{2} Re {\cos (θ) [E_{θ}^{t} H_{ϕ}^{t *} - E_{ϕ}^{t} H_{θ}^{t *}] - \sin (θ) [E_{ϕ}^{t} H_{r}^{t *} - E_{r}^{t} H_{ϕ}^{t *}]},

where the superscript t indicates the total (sum of the incident and scattered) field, and the superscript * represents the complex conjugate. The variables θ and ϕ represent the polar and azimuthal coordinates and as subscripts they represent the corresponding components of a vector field, respectively. On the other hand, the phase profile after the sphere can be obtained from the argument of the polar-angle component of the electric field when the profile is calculated at a large distance that the radial component can be ignored. The error of the Rytov approximation in predicting the intensity and phase profiles after a homogeneous sphere may be defined as

E r r_{I} = \frac{\int | I_{M i e} - I_{R y t o v} | d r}{\int | I_{M i e} | d r},

E r r_{φ} = \frac{\int | φ_{M i e} - φ_{R y t o v} | d r}{\int | φ_{M i e} | d r},

where I and φ are the intensity and phase profiles after the sphere, respectively.

Figure 3 plots Err_I and Err_φ for different values of V. For less than 1% accuracy in intensity profile, one needs V < 0.0049, while for the same accuracy in phase profile, one needs V < 0.030. The validity condition is less stringent in case of estimating the phase profile than the intensity profile. However, the error quickly increases for a large value of V.

Fig. 3 Error of the first-order Rytov approximation in predicting intensity and phase profiles after a homogeneous sphere.

Download Full Size | PDF

In the example we considered above, imaging 5 cm radius sphere at 1 m distance using the 60 keV incident beam, we obtained V = 0.0027, and this corresponds to intensity error of 0.38% and phase error of 0.18%. Figure 4 compares the results obtained with the Rytov approximation and the Mie series for the model system (λ = 0.256 nm, R = 2.56 μm, and z = 32.2 μm). Figures 4(c) and 4(d) are the cross-sections of Figs. 4(a) and 4(b), respectively. Figures 4(e,f) plot the difference in the intensity and profiles, respectively, provided by the Rytov approximation and the Mie series. Note that the center part of the profile in Fig. 4(c) is smaller than one due to the material absorption of the sphere, while the oscillation near the edges is the phase signature. Figure 4(e) shows that the Rytov approximation provides a more accurate profile in the center part than the edge regions.

Fig. 4 (a) Intensity and (b) phase images obtained by the first-order Rytov approximation. (c) and (d) are cross-section profiles across the center of (a) and (b), respectively. (e) Difference between (c) and the corresponding profile obtained by the exact Mie solution. (f) Difference between (d) and the corresponding profile obtained by the exact Mie solution.

Download Full Size | PDF

Figure 5 compares the cross-section profiles calculated by the Rytov approximation and the exact Mie solution near the edge of the sphere, the dotted region in Fig. 4(c). The two curves accurately match in the outer part of the sphere, which corresponds to the detector coordinate larger than 2.5 μm. However, in the inner part, the profile given by the Rytov approximation is shifted to the inside compared the Mie solution, which is responsible for the large error near the edge regions in Fig. 4(a). The reason of this shift is not clear.

Fig. 5 Comparison of the near-edge profiles calculated by the Rytov approximation and the exact Mie solution.

Download Full Size | PDF

Suppose we fix the energy of the incident beam at 60 keV, and measure the intensity and phase at 1 m distance from the center of the sphere. Using Eq. (12) and Fig. 3, we can calculate the maximum radius of the sphere R_max, which renders the Rytov approximation to be valid within 1% accuracy. For the water sphere, we obtain R_max = 9.16 cm for the intensity measurement and R_max = 58.0 cm for the phase measurement. These values are reasonable for medical applications or luggage inspection.

5. Conclusion

In this study, we investigated the validity and accuracy of the Rytov approximation in predicting the intensity and phase profiles of the diffracted beam after a homogeneous sphere. We showed that the validity condition of the Rytov approximation can be expressed in terms of a single parameter V, which is defined in terms of the complex refractive index of sphere and the Fresnel number. In comparison with the exact Mie solution, we calculated the accuracy of the Rytov approximation for different values of V. Using the principle of similarity, we estimated the maximum size of water sphere that can be accurately simulated with the first-order Rytov approximation.

Acknowledgments

This research was funded by the Department of Homeland Security’s Science and Technology Directorate (Contract No. 6924804) and the National Research Foundation of Singapore through the Singapore-MIT Alliance for Research and Technology Centre.

References and links

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

2. P. C. Diemoz, A. Bravin, and P. Coan, “Theoretical comparison of three X-ray phase-contrast imaging techniques: propagation-based imaging, analyzer-based imaging and grating interferometry,” Opt. Express 20(3), 2789–2805 (2012). [CrossRef] [PubMed]

3. E. Förster, K. Goetz, and P. Zaumseil, “Double crystal diffractometry for the characterization of targets for laser fusion experiments,” Kristall und Technik 15(8), 937–945 (1980). [CrossRef]

4. D. Chapman, W. Thomlinson, R. E. Johnston, D. Washburn, E. Pisano, N. Gmür, Z. Zhong, R. Menk, F. Arfelli, and D. Sayers, “Diffraction enhanced x-ray imaging,” Phys. Med. Biol. 42(11), 2015–2025 (1997). [CrossRef] [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. A. Momose, W. Yashiro, Y. Takeda, Y. Suzuki, and T. Hattori, “Phase tomography by X-ray Talbot interferometry for biological imaging,” Jpn. J. Appl. Phys. 45(6A), 5254–5262 (2006). [CrossRef]

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

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

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

10. X. Wu, H. Liu, and A. Yan, “Optimization of X-ray phase-contrast imaging based on in-line holography,” Nucl. Instrum. Methods Phys. Res. B 234(4), 563–572 (2005). [CrossRef]

11. A. Snigirev, I. Snigireva, V. Kohn, S. Kuznetsov, and I. Schelokov, “On the possibilities of x‐ray phase contrast microimaging by coherent high‐energy synchrotron radiation,” Rev. Sci. Instrum. 66(12), 5486–5492 (1995). [CrossRef]

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

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

14. J. B. Keller, “Accuracy and Validity of the Born and Rytov Approximations,” J. Opt. Soc. Am. 59, 1003–1004 (1969).

15. F. Lin and M. Fiddy, “The Born-Rytov controversy: I. Comparing analytical and approximate expressions for the one-dimensional deterministic case,” J. Opt. Soc. Am. A 9(7), 1102–1110 (1992). [CrossRef]

16. E. Kirkinis, “Renormalization group interpretation of the Born and Rytov approximations,” J. Opt. Soc. Am. A 25(10), 2499–2508 (2008). [CrossRef] [PubMed]

17. M. Slaney and A. Kak, Principles of Computerized Tomographic imaging (SIAM, 1988).

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

19. M. Slaney, A. C. Kak, and L. E. Larsen, “Limitations of imaging with first-order diffraction tomography,” IEEE Trans. Microw. Theory Tech. 32(8), 860–874 (1984). [CrossRef]

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

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

22. C. F. Boliren and D. R. Huffman, Absorption and Scattering of Light by Small Particles (J Wiley & Sons, 1983).

23. A. A. Neves and D. Pisignano, “Effect of finite terms on the truncation error of Mie series,” Opt. Lett. 37(12), 2418–2420 (2012). [CrossRef] [PubMed]

24. W. G. Melbourne, Radio Occultations Using Earth Satellites: a Wave Theory Treatment (Wiley-Interscience, 2005).

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

26. F. Slimani, G. Grehan, G. Gouesbet, and D. Allano, “Near-field Lorenz-Mie theory and its application to microholography,” Appl. Opt. 23(22), 4140–4148 (1984). [CrossRef] [PubMed]

Rytov approximation for x-ray phase imaging

Abstract

1. Introduction

2. Scattered field calculation under the first Rytov approximation

3. Explicit validity condition of the first Rytov approximation for a homogeneous sphere

4. Error estimation of the Rytov approximation: comparison with Mie solution

5. Conclusion

Acknowledgments

References and links

Cited By

Figures (5)

Equations (17)

Optics Express