1. Introduction

In recent years, X-ray grating interferometry, a differential x-ray phase-contrast imaging technique, has been an active field of x-ray imaging research. This technique enables highly sensitive phase contrast imaging for soft tissues and low-Z materials, hence it has many potential applications in medical imaging and material science [1–20]. As is schematically shown in Fig. 1, this phase contrast imaging technique usually employs a Talbot-Lau interferometer [4, 20–22], which consists of an X-ray source, a source grating G₀, a phase grating G₁ and an imaging detector. The setup shown in Fig. 1 is also called Talbot-Lau interferometer in inverse geometry [23]. Note that in another Talbot-Lau interferometer setup one also employs an absorbing grating placed in front of the detector as the analyzer to extract the fringe shift generated by the sample [4, 20–22]. The beam hardening correction method described below in this work is applicable to all the Talbot-Lau interferometer setups. Serving as an aperture mask, the source grating G₀ breaks up the source spot into an array of narrow sources to improve the spatial coherence of X-ray illumination. The phase grating G₁ serves as a beam splitter to generate interference fringes with good intensity modulation [4, 20]. With presence of a sample, the information of the sample such as its x-ray attenuation and phase gradients are all encoded in the interference fringes. For one-dimensional phase grating interferometers with a monochromatic source of photon energy E, the Fourier expansion of the fringe intensity pattern I(M_gx, M_gy; E) is given by the following expression:

 

Fig. 1 Schematic of an x-ray phase grating interferometer with an X-ray source and source grating.

Download Full Size | PDF

The most important goal of the grating based differential phase contrast imaging is to accurately reconstruct the sample’s electron density distribution. To determine the electron density of the sample, the exact fringe phase shift of at a given energy is needed. To retrieve the fringe phase shifts ϕ_m(x, y; E) from the fringe pattern I(M_gx, M_gy; E), one can use either the phase stepping method or the Fourier analysis method [4–7, 13, 16, 17]. Once the fringe phase shift ϕ_m(x, y; E) is retrieved, one can simply recover the sample’s density gradients ∂ρ_e,p(x, y)/∂x by inverting Eq. (2). In practice, the selection of a specific diffracted order for phase retrieval depends on the grating setup. Owing to limited spatial coherence of the source assembly, the dominant diffraction orders are m = ±1 for π/2-phase grating setups. Consequently, one needs to retrieve ϕ₁(x, y; E) in these cases. For π-grating setups, the m = ±2 orders dominate instead, so ϕ₂(x, y; E) needs to be retrieved. In the monochromatic case, the fringe phase shifts ϕ₁(x, y; E) or ϕ₂(x, y; E) can be retrieved using either the phase stepping method or the Fourier analysis method. In the phase stepping method, an absorbing grating, which is of the same pitch as the fringe pattern, is placed in front of the detector. One scans the grating in steps to generate an intensity oscillation curve at each of the pixels. The oscillation phase gives the fringe phase shift [4, 20]. On the other hand, when the fringe is resolvable by the detector pixels, one may apply the Fourier analysis method for phase retrieval. In this method one crops the fringe’s Fourier spectrum around the diffraction peak to retrieve the fringe phase shift using inverse Fourier transform [18, 19, 21, 22]. Repeating the phase retrievals for a complete set of angular projections, one can reconstruct the quantitative images of sample electron density distribution [4, 15–17, 20].

In medical imaging and non-destructive testing, the differential phase contrast imaging technique is implemented mostly with polychromatic sources [4, 16, 17]. Under polychromatic illumination, the measured intensity fringe pattern is a sum of the intensity fringes of different photon energies, and each of them with a different fringe phase shift. Using Eq. (1), we can write this sum as:

The fringe shift distortion in polychromatic X-ray grating interferometry has two causes. First, the fringe phase shift and fringe visibility are all photon energy dependent, hence polychromatic fringe shift is a nonlinear function of the monochromatic fringe shifts of different energies [24]. The second cause of fringe shift distortion is the beam hardening, which refers to the non-uniform x-ray spectrum shifts caused by the non-uniform attenuation across the sample. The fringe shifts are further distorted by this uneven x-ray spectrum. In literature there are several published studies on how to correct the polychromatic effects on fringe phase shift. These studies conducted calibration experiments based on two concepts, namely, the effective energy and the ray-sum linearization [25–27]. The concept of effective energy hypothesizes that the polychromatic and monochromatic fringe phase shifts are proportional to each other and with the same sign. Specifically, it assumes that the polychromatic fringe phase shift ϕ_m,Poly(x, y) would be equivalent to a monochromatic fringe phase shift at the effective energy E_eff, that is, ϕ_m,Poly(x, y) = ϕ_m(x, y; E = E_eff) [25–27]. In this approach one needs to employ a phantom of reference material with known electron density and geometry to experimentally determine the effective energy value. Nevertheless, a sample may have a range of effective energies if the sample has uneven attenuation [26, 27]. This is because the effective energy tends to increase with increasing attenuation. Moreover, the effective energy concept is valid only for small fringe shifts of few degrees, otherwise the polychromatic shift is a non-linear function of monochromatic fringe phase shift. Furthermore, as we showed previously [24], sometimes these two fringe shifts may even have opposite signs, thereby no effective energy exists. As for the second concept, the beam hardening refers to that the x-ray exiting from more attenuated ray paths is with a higher fraction of high-energy photons than that with the original spectrum. Hence variable attenuation of the sample causes non-uniform spectrum shifts. In grating interferometry, the beam hardening causes significant errors in retrieving the density gradients of the sample, and consequently ruins the potential of quantitative imaging with the grating interferometry. In the literature, an experimental calibration procedure, namely the so-called ray-sum linearization technique, was reported for the beam hardening correction [25–27]. In this procedure one conducts calibration experiment with homogeneous sample of known attenuation coefficients, electron densities and thickness profile. This procedure starts to determine two effective energy values, one for attenuation and one for phase shift of the sample, from the polychromatic measurement. Assuming the effective energy is independent of the sample thickness, one determines the sample thickness from the measured polychromatic attenuation and sample phase shift. Comparing thus extracted thickness to the reference thickness profile, one derives the correction factors that converts the measured sample thickness to the true thickness profile. These correction factors are then used to correct the measured fringe phase shift values [26]. This technique is valid only for a homogeneous sample with known attenuation coefficients, electron density and geometry, although samples are mostly inhomogeneous. The burdensome procedure has to be repeated once any setup parameter of the grating interferometer is changed. In addition to these calibration approaches discussed, some researchers resorted to laborious wave-propagation simulations to analyze the relationship between the polychromatic and monochromatic fringe shifts [21, 22].

In order to reconstruct the electron density gradients of the sample, we recently developed an analytic approach that enables retrieving the monochromatic phase shift ϕ_m(x, y; E) from the polychromatic measurement [24]. We validated this approach for a number of interferometry setup geometries by using the wave-propagation based simulations [24]. However, our previous approach doesn’t address the X-ray beam hardening effects. In fact, in our previous approach we assumed the knowledge of the effective X-ray spectrum, which is assumed to be uniform across the sample. In this work, we address the uneven beam hardening problem by developing a new approach. In section 2, we present a novel analytic approach to incorporate the beam hardening correction in the fringe shift retrieval. In section 3, we apply these methods to retrieve the monochromatic fringe shift in the cases with uniform or approximately uniform energy dependence. Through several wave-propagation based simulations, we validate our phase retrieval formulas and demonstrate robustness of our method against noise. We discuss the limitations of our methods and conclude the paper in section 4.

2. Method

We assume that no absorption edge appears in the x-ray spectral range studied in this work. In order to account for the beam hardening effects, we need to recover the monochromatic attenuation function, which can be written as A²(x, y; E) = exp [−p(x, y; E), where p(x, y; E) denotes a ray-sum acquired at energy E. A ray-sum is an integral of the linear attenuation coefficient (LAC) µ(x, y, z; E) over the ray z, that is, p(x, y; E) = ∫_Ray µ(x, y, s; E) ds. In grating interferometry one can measure the polychromatic attenuation $A_{\text{Poly}}^2(x,y)$ from the zero-th diffracted order of the fringe pattern [4]. Note that polychromatic attenuation $A_{\text{Poly}}^2(x,y)$ is a spectral average of A²(x, y; E), that is, $A_{\text{Poly}}^2(x,y)={\displaystyle \int \mathfrak{D}(\mathrm{E})}\cdot A^2(x,y;\mathrm{E})\text{dE}$, where $\mathfrak{D}(\mathrm{E})$ is the spectrum. To extract the monochromatic attenuation function A²(x, y; E) from the polychromatic attenuation function, we observed that, one needs to know the energy functional form of the sample’s mass attenuation coefficient. Hence, in this work, we assume that the sample’s mass attenuation coefficient µ/ρ(x, y, z; E) scales with photon energy as µ/ρ(x, y, z; E) ∝ f (E), where f (E) is a known functional form of photon energy common to all substances in the sample. We call this assumption the uniform energy dependence condition, and f (E) the energy factor. Obviously any single-material sample with known mass attenuation coefficient satisfies this condition exactly by setting the energy factor f (E) = µ/ρ(E). In practice, inhomogeneous samples satisfy this condition only approximately under certain circumstances, as we will discuss in details below. Under this assumption of the uniform energy dependence, we can rewrite the attenuation function as $A^2(x,y;\mathrm{E})=\exp \left[-p_{\mathrm{D}}(x,y)\times \frac{f(\mathrm{E})}{f({\mathrm{E}}_{\mathrm{D}})}\right]$, where p_D(x, y) denotes the monochromatic ray-sum at the design energy, p_D(x, y) ≡ ∫_Ray µ(x, y, s; E_D) ds, and we have

To extract the monochromatic attenuation function A²(x, y; E), we solve Eq. (4) for p_D(x, y) from the measured polychromatic attenuation $A_{\text{Poly}}^2(x,y)$ on a pixel-by-pixel basis. For a given pixel, Eq. (4) constitutes a nonlinear equation of p_D(x, y). We developed an iterative algorithm to solve for p_D(x, y).

Consider the following functionals

One can easily see that H(p(x, y)) is a strictly convex functional, since the function f (x) = exp[−x] is strictly convex on ℝ, i.e. f″(x) > 0, for all x ∈ ℝ. We want to retrieve the function p_D(x, y), from given value $A_{\text{Poly}}^2(x,y)=H(p_{\mathrm{D}}(x,y))$ through iterative method: assuming p_k(x, y) is acquired after the k-th iteration, and $\mathrm{\Delta }(x,y)=A_{\text{Poly}}^2(x,y)-H(p_k(x,y))$ is the difference. We want to diminish the difference Δ(x, y) through the next step p_k+1(x, y) = p_k(x, y) + δ(x, y), i.e. we hope $H(p_{k+1}(x,y))=A_{\text{Poly}}^2(x,y)$. Then Δ(x, y) ≈ H(p_k+1(x, y))−H(p_k(x, y)) ≈ ∇H(p_k(x, y)) × δ(x, y), with the assumption of δ is small. This gives

So we can develop the algorithm as follows.

Algorithm 1 Algorithm to retrieve the monochromatic ray-sum $p_{\mathrm{D}}^{\text{Iter}}(x,y)$.

The algorithm is in fact a gradient decent algorithm [28]. Since H is strictly convex, a global minimum is guaranteed and the converging speed is fast. With our simulation data in following sections, we find the error ||δ(x, y)||₂ would fall under 10⁻¹³ after 15 iteration steps.

Once the sample attenuation function $A^2(x,y;\mathrm{E})=\exp \left[-p_{\mathrm{D}}^{\text{Iter}}(x,y)\cdot f(\mathrm{E})/f({\mathrm{E}}_{\mathrm{D}})\right]$ is extracted, we may incorporate this attenuation-filtration effect into a position-dependent eiffective spectrum as

The position-dependent effective spectrum Θ(E; x, y) is normalized as ∫Θ(E; x, y) dE = 1. Therefore, in terms of the effective spectrum Θ(E; x, y), we may rewrite the polychromatic intensity fringe pattern as

From above equation, the polychromatic fringe phase shift ϕ_m,Poly(x, y) is given by

In this formula Δϕ(E) is the phase shift step of the phase grating at energy E, and λ is x-ray wavelength. The factor ${\left(-1\right)}^{\lfloor 4k\lambda R_2/M_g{p}_1^2\rfloor }$ in this formula is defined with the floor-function ⌊x⌋, which is defined as the largest integer that is less or equal to x, so this factor takes values of 1 or −1, depending its exponent. Although Eq. (10) is the key basic equation for the retrieval of monochromatic fringe phase shifts, it is not a convenient formula to use for retrieval of ϕ_m,D(x, y). In a previous work [24], we developed an iterative approach to retrieve ϕ_m,D(x, y) from the polychromatic measurement provided the effective spectrum is uniform and independent of pixel position. In this work, we extended that approach to the cases with pixel position-dependent effective spectrum to account for the uneven beam hardening effects. To do so, we first compute the position-dependent coefficients Q_m(n; x, y), which is defined as:

We call Q_m(n; x, y) the n-th optical response coefficient of an interferometer setup. Note that the optical response coefficient is position dependent, since the effective spectrum is position-dependent due to the beam hardening effect. Following the reasoning and derivation in [24], we extended the previous iterative phase retrieval method to the cases with position-dependent X-ray spectra. Assume that ${\varphi }_{m,\mathrm{D}}^{(q)}(x,y)$ is an approximation of the monochromatic fringe shift after the q-th iteration, then the updated solution of ϕ_m,D(x, y) in the (q +1)-th iteration is given by

The iteration stops when ${\Vert {\varphi }_{m,\mathrm{D}}^{(q+1)}-{\varphi }_{m,\mathrm{D}}^{(q)}\Vert }_2<\mathit{ϵ}$, where the operator ||·||₂ indicates the L₂-norm and ϵ is a designated small number representing the error allowance. The smaller the ϵ is, the higher the accuracy one will achieve for the solution of ϕ_m,D(x, y).

3. Results

To validate the above formulas and algorithms for retrieving the monochromatic fringe shifts, we conducted numerical simulations of polychromatic X-ray interferometry with phantoms of breast tissues with known electron density distributions. For sake of clarity in presentation, we assumed that the thickness profiles of the phantoms vary only in x-direction. Specifically, Fig. 2 shows a 3D manifest of one of these phantoms. As is shown in the figure, the top layer (layer I) has an uneven thickness profile, but the profile varies only along the x-direction. Other layers (layers II, III and IV) all have even thickness profiles. Each layer is filled with a specific type of tissues or materials in the simulations. As is shown in Fig. 2, the x-ray projection thickness takes the maximum at the left plateau, then it decreases with a constant rate down the left slope, and it reaches the minimum at the flat valley. The phantom is left-right symmetric about the center. With this phantom design, the electron density gradient ∂ρ_e,p(x)/∂x is diminishing at the two plateaus and the valley, but takes constant magnitude with opposite signs along the two slopes. So does the monochromatic fringe shift, as which is proportional to ∂ρ_e,p(x)/∂x. By the phantom design, the phantom should generate fringe shifts of ±45° degree along its two slopes at this grating design energy. On the other hand, as x-ray traverses the phantom, different projection-thickness generates different X-ray spectrum shifts. Hence, for this phantom design, the beam hardening effects will manifest itself as the uneven magnitudes of the polychromatic fringe shift ϕ_m,Poly(x) over the two slopes. We will demonstrate how our retrieval method manages to recover accurately the monochromatic fringe shift ϕ_m,D(x) at the design energy. In the simulations, all the tested phantoms have the same design scheme but with different layer of materials. Numerical simulations of the wave propagation based on the Fresnel diffraction were conducted [29] for several interferometer setups. In subsection 3.1, we first consider the cases that satisfies the uniform energy dependence exactly. We test if the retrieval algorithms work well in the ideal cases of the uniform energy dependence. In subsection 3.2, we then consider the cases that satisfies the uniform energy dependence condition only approximately. These cases are more interesting from a practical point of view.

 

Fig. 2 3D manifest of the phantom employed in the simulations.

Download Full Size | PDF

3.1. Uniform energy dependence

As an ideal case of uniform energy dependence, we first consider a homogeneous phantom of glandular tissues with known mass attenuation coefficient µ/ρ_Gland(E). The shape of the phantom is as shown in Fig. 2 except it is made up of only one layer, the layer I. We calculated glandular mass attenuation coefficient according to µ/ρ_Gland(E) = ∑_iw_i(µ/ρ)_i(E), where (µ/ρ)_i(E) is the mass attenuation coefficient of the i-th constituent element in the tissue, and w_i is the elemental weight fraction. The elemental weight fractions are obtained from the published measurement data [30]. The experimentally determined elemental mass attenuation coefficients were tabulated in the literature [31], and these data were interpolated for all photon energies in the spectrum range. Substituting f (E) = µ/ρ_Gland(E) into Eq. (4), we retrieved the monochromatic ray-sum $p_{\mathrm{D}}^{\text{Iter}}(x)$ by using the iterative algorithm of Eqs. (5)–(7). In this way we extracted the monochromatic attenuation function $A^2(x;\mathrm{E})=\exp \left[-p_{\mathrm{D}}^{\text{Iter}}(x)\cdot f(\mathrm{E})/f({\mathrm{E}}_{\mathrm{D}})\right]$, and obtained the position-dependent effective spectrum Θ(E; x, y) by using Eq. (8).

The phase retrieval depends on the specifics of the interferometer setup employed in the phantom imaging. The simulated setup starts from the source selection. Since the coherence degree is independent of photon energy, thus the polychromatic fringe phase shift does not depend on the focal spot size and shape. Without loss of generality, we assumed a point-like focal spot in the simulations. In addition, we assumed a 26 kVp effective spectrum, as depicted in Fig. 3. The other interferometer setup parameters selected in the simulations are as follows. The phase grating is a π/2-grating of 5µm pitch at the design energy E_D = 17 keV. The reduced grating-to-detector distance R₂/M_g (with M_g = 1) was set to the 1 Talbot distance $R_2/M_g{p}_1^{2}2{\lambda }_{\mathrm{D}}$, where p₁ is the grating pitch and λ_D denotes x-ray wavelength at the design energy. In the π/2-grating setup, the dominant diffraction orders in the fringe pattern are m = ±1 orders. Using the fringe visibility coefficient C₁(E) given by Eq. (11), we found that the n-th position-dependent optical response coefficient for this setup is given by:

 

Fig. 3 26 kVp x-ray effective spectrum employed in the simulations.

Download Full Size | PDF

As for phase retrievals, the polychromatic fringe phase shift ϕ_1,Poly(x) was retrieved from the fringe pattern by using the Fourier analysis of the m = 1 diffracted peak [18, 19, 21, 22]. Based on the extracted ϕ_1,Poly(x) data and the computed optical response coefficients Q₁(n; x, y), we determined the monochromatic fringe shift ${\varphi }_{1,\text{Poly}}^{\text{Iter}}(x)$ at the design energy by using the iterative phase retrieval algorithm summarized in Eq. (13).

Figures 4 and 5 show the retrieving results for monochromatic attenuation and fringe phase shift. In Fig. (4), the dash-dot green profile depicts the polychromatic ray-sum $p_{\text{Ploly}}(x)=-\log \left(A_{\text{Poly}}^2(x)\right)$ measured from the zero-th order of the simulated fringe pattern. The blue trace plots the theoretical values of the monochromatic ray-sum $p_{\mathrm{D}}^{\text{Theo}}(x)$ at E_D = 17 keV. As is shown in Fig. (4), the dashed red trace is the iteratively retrieved monochromatic ray-sum $p_{\mathrm{D}}^{\text{Iter}}(x)$ at the design energy E_D = 17 keV. The retrieved monochromatic ray-sum $p_{\mathrm{D}}^{\text{Iter}}(x)$ agrees with the theoretical values $p_{\mathrm{D}}^{\text{Theo}}(x)$ within 0.03% through 11 iteration steps.

 

Fig. 4 Plots of the ray-sums: The design energy is set to E_D = 17 keV for π/2-grating interferometer with 26 kVp x-ray effective spectrum shown in Fig. 3. In this simulation, the phantom is made of homogeneous glandular tissue with the shape as shown in Fig. 2 except it has only one layer (layer I). The total height of the phantom is 40.5 mm. In the plot, the solid blue line, $p_{\mathrm{D}}^{\text{Theo}}$, represents the theoretical ray-sum at design energy E_D. The dash-dot green line, p_Poly(x), is the polychromatic ray-sum directly computed from $A_{\text{Poly}}^2(x)(p_{\text{Poly}}(x)=-\log (A_{\text{Poly}}^2(x))$, which is derived from the zero-th diffracted order of the fringe pattern. The dashed red line, $p_{\mathrm{D}}^{\text{Iter}}(x)$, is the retrieved ray-sum at design energy E_D using Algorithm 1. After 11 iterations, the relative error falls in 0.03%.

Download Full Size | PDF

 

Fig. 5 Plots of the fringe shifts: The solid blue line ${\varphi }_{1,\mathrm{D}}^{\text{Theo}}(x)$ represents the true theoretical values which are set to ±45° degrees at the two bumps. The dash-dot green line, ϕ_1,Poly, is the polychromatic fringe phase shifts retrieved from the simulated fringe pattern by using the Fourier analysis method. The dashed red line, ${\varphi }_{1,\mathrm{D}}^{\text{Iter}}(x)$, is the iterative solution using Eq. (13). After 29 iteration steps, the relative error falls in 0.11%.

Download Full Size | PDF

Figure 5 shows the results for the fringe phase shifts in this simulation. In Fig. 5, the solid blue trace plots the theoretical values of the monochromatic fringe shift ${\varphi }_{1,\mathrm{D}}^{\text{Theo}}(x)$ at E_D = 17 keV, while the dash-dot green trace plots the polychromatic fringe shift ϕ_1,Poly(x) measured from the simulated fringe pattern. Comparing the green and blue traces, there are gross discrepancies between ϕ_1,Poly(x) and ${\varphi }_{1,\mathrm{D}}^{\text{Theo}}(x)$. Moreover, note that the discrepancies are larger in magnitude at the plateau sides than that at the central valley side. In other words, the green trace (ϕ_1,Poly(x)) is not parallel to the blue trace $({\varphi }_{1,\mathrm{D}}^{\text{Theo}}(x))$. This is just a manifestation of the beam hardening effects, as the beam is harder at the plateau side than that at the valley side. In Fig. 5 the dashed red trace plots the retrieved monochromatic fringe shift ${\varphi }_{1,\mathrm{D}}^{\text{Iter}}(x)$ with relative error of 0.11% after 29 iteration steps. This excellent agreement validates out the monochromatic fringe shift retrieval algorithms of Eq. (13) and the monochromatic ray-sum retrieval Algorithm 1.

3.2. Approximately uniform energy dependence: single power energy dependence

Tissue mass attenuation coefficients generally do not exhibit uniform energy dependence. In fact, the elemental mass attenuation coefficients in diagnostic X-ray imaging can be well parameterized with multiple powers of photon energy, such as µ/ρ(E) = aZ⁵E^−3.5 +bZ²E⁻² +cZσ_KN(E), where Z is the elemental atomic number, σ_KN(E) denotes Klein-Nishina total cross-section, and a, b, c all are constants. Of course, there are many other but similar parameterizations reported in the literature [32]. Such a multiple power parameterization arises because the three contributing x-ray-matter interactions, namely the photoelectric absorption, the coherent and incoherent scattering have different characteristics in photon energy and atomic numbers [32]. As a sample is generally made up of different tissue or materials, which in turn contain different elements, hence one generally does not know the proportions of each of the energy powers for a sample to be imaged. However, the grating based x-ray phase contrast imaging currently is oftentimes carried out with low tube voltages up to about 50 kVp, and the spectral energy range is not very wide. We assumed that, in such a circumstance, the mass attenuation coefficients of the sample is given by µ/ρ(x, y, z; E) = c(x, y, z)E⁻ⁿ, where n is a number common to all the materials in a given sample. Here c(x, y, z) is an unknown function to account for variations in material composition over sample voxels. We call this assumption as the single power energy dependence. Under assumption, the energy factor in Eq. (5) becomes f (E) = E⁻ⁿ, thereby one can apply the Algorithm 1 to retrieve the monochromatic attenuation function A²(x, y; E). With low tube voltages up to about 50 kVp, we found that in general n is ranging from 2 to 3.5. With this uncertainty in the n-value, We should make sure that different n-values all result in reasonable accuracies in monochromatic fringe shift retrieval. Recently progress has been made in high-energy grating based imaging with photon energies as high as 100 keV [33–35]. It will be interesting to explore in future if the assumption of the single power energy dependence still works fine with such a high-energy grating technique.

To this purpose we reset the phantom to have multiple layers of different tissues. As is shown in Fig. 2, the phantom has 4 layers. The layers I, II, III, and IV represent 100% glandular tissue, 75% glandular and 25% adipose tissues, 50% glandular and 50% adipose tissues, and 25% glandular and 75% adipose tissues, respectively. The total height of the phantom is still 40.5 mm. The interferometer setup is the same as that described in secion 3.1. With the same 26 kVp beam and π/2-phase grating setup, we investigated the monochromatic fringe shift retrieval by using of the energy factors f (E) = E⁻ⁿ.

Figure 6 shows the retrieving results. The solid blue line ${\varphi }_{1,\mathrm{D}}^{\text{Theo}}(x)$ represents the true theoretical values which are set to ±45° degrees over the two slopes of the phantom. The dash-dot green line, ϕ_1,Poly(x), is the polychromatic fringe phase shifts retrieved from the simulated fringe pattern by using the Fourier analysis method. Comparing the green and blue traces, there are again gross discrepancies between ϕ_1,Poly(x) and ${\varphi }_{1,\mathrm{D}}^{\text{Theo}}(x)$. While the blue trace (ϕ_1,D (x)) is leveled in X-direction, the green trace (ϕ_1,Poly(x)) becomes slanted over the phantom’s slopes. This is just the beam hardening effects, as we discussed in details in section 3.1. The retrieved iterative solution using Eq. (13) with different position-dependent effective spectrums Θ(E; x) corresponding to f (E) = E⁻ⁿ, with n = 2, 2.5, 3, 3.5, are shown in the traces of different colors in the figure. There are two aspects to check when one compares these traces with the solid blue trace of $({\varphi }_{1,\mathrm{D}}^{\text{Theo}}(x))$. First, one checks the discrepancies between the retrieved values and the theoretical values of fringe shifts. We noted that the fringe shifts retrieved with n = 2, 2.5, 3 and 3.5 have errors 2.32%, 0.76%, 3.43% and 5.87% respectively, as is shown in Fig. 6. So the retrieval with n = 2.5 achieved the least error. Second, one checks if these traces are parallel to the blue trace $({\varphi }_{1,\mathrm{D}}^{\text{Theo}}(x))$ in each of their respective sections. In other words, it is to check the extent of beam hardening correction achieved through the phase retrieval. We noted that while the fringe shifts retrieved with n = 2 and 2.5 achieved almost perfect beam hardening correction, that with n = 3 and 3.5 over-compensate the beam hardening effect, as is shown in Fig. 6. The simulation with this phantom shows that the fringe shifts retrieved with n = 2.5 is the most accurate as compared to other n-values.

 

Fig. 6 Plots of the fringe shifts with different f (E) approximation: In this simulation, the phantom is reset to have multiple layers of different tissues with the shape shown in Fig. 2. There are total of 4 layers. The layers I, II, III, and IV represent 100% glandular tissue, 75% glandular and 25% adipose tissues, 50% glandular and 50% adipose tissues, and 25% glandular and 75% adipose tissues respectively. The total height of the phantom still 40.5 mm. The solid blue line ${\varphi }_{1,\mathrm{D}}^{\text{Theo}}(x)$ represents the true theoretical values which are set to ±45° degrees at the two bumps. The dash-dot green line, ϕ_1,Poly, is the polychromatic fringe phase shifts retrieved from the simulated fringe pattern by using the Fourier analysis method. The retrieved iterative solution using Eq. (13) with different position-dependent effective spectrums Θ(E; x) corresponding to f (E) = E⁻ⁿ, with n = 2, 2.5, 3, and 3.5 are shown in different colors in the figure.

Download Full Size | PDF

In another simulation, we replaced a part of layer IV of the previous phantom with a 3.5 mm layer of aluminum, a heavier material than glandular and adipose tissues, with otherwise the same interferometry setup. As is shown in Fig. 7, the errors of the fringe shifts retrieved with n = 2, 2.5, 3, 3.5 are 7.79%, 1.28%, 5.61%, and 12.14% respectively. The case of n = 2.5 is still the most accurate among the four exponent choices. When we take a closer look at Fig. 6 and Fig. 7, we realize a subtle difference between the two cases. In Fig. 6, the left bump in the yellow trace (n = 2.5) is slightly above the theoretical value. In Fig. 7, the yellow trace is slightly below the theoretical value. This indicates that with less penetrating samples one should use a larger n-value for beam hardening correction.

 

Fig. 7 Plots of the fringe shifts with different f (E) approximation: In this simulation, a part layer IV of the previous phantom is replaced with a 3.5 mm layer of aluminum. All other setups are the same. The solid blue line ${\varphi }_{1,\mathrm{D}}^{\text{Theo}}(x)$ represents the true theoretical values which are set to ±45° degrees at the two bumps. The dash-dot green line, ϕ_1,Poly, is the polychromatic fringe phase shifts retrieved from the simulated fringe pattern by using the Fourier analysis method. The retrieved iterative solution using Eq. (13) with different position-dependent effective spectrums Θ(E; x) corresponding to f (E) = E⁻ⁿ, with n = 2, 2.5, 3, and 3.5 are shown in different colors in the figure.

Download Full Size | PDF

3.3. Cases with noisy data

To evaluate the performance of our algorithm with noisy data, we added some Gaussian white noise (5%) to the fringe intensity image I and the reference grating only image. As is stated in the previous subsection, n = 2.5 in the energy factor f (E) = E⁻ⁿ is a good approximation. So we set f (E) = E^−2.5 in the monochromatic ray-sum retrieval algorithm and apply the retrieved monochromatic ray-sum $p_{\mathrm{D}}^{\text{Iter}}(x)$ to the position-dependent effective spectrums Θ(E; x). We then retrieved the monochromatic fringe phase shift ${\varphi }_{1,\mathrm{D}}^{\text{Iter}}(x)$. The results are shown in Fig. 8. The retrieved fringe phases attain good signal noise ratios. For example, the SNRs for ϕ_1,Poly(x) and ${\varphi }_{1,\mathrm{D}}^{\text{Iter}}(x)$ at the left bump are 32.8 and 52.7 respectively (Fig. 8). This result demonstrates that our phase retrieval method is robust. It can be seen the noise level towards the two ends is higher than the central. This is because the phantom is thicker towards the two ends. This being so, the relative noise level of the retrieved monochromatic attenuation function at the two ends is higher. This higher noise level at the two ends inevitably permeates through the position-dependent effective spectrum Θ(E; x, y). Consequently, the retrieved monochromatic fringe phase shift becomes noisier at the two ends.

 

Fig. 8 Performance of the iterative Eq. (13) with noisy data. The setup and phantom is the same as Fig. 6. We add some Gaussian white noise (5%) to fringe intensity image I and the reference grating only image. The solid blue line ${\varphi }_{1,\mathrm{D}}^{\text{Theo}}(x)$ represents the true theoretical values which are set to ±45° degrees. The dash-dot green line, ϕ_1,Poly, is the noisy polychromatic fringe phase shifts retrieved from the simulated fringe pattern by using the Fourier analysis method. The dashed red line, ${\varphi }_{1,\mathrm{D}}^{\text{Iter}}(x)$, is the iterative solution using Eq. (13) with f (E) = E^−2.5. The SNRs of ϕ_1,Poly(x) and ${\varphi }_{1,\mathrm{D}}^{\text{Iter}}(x)$ at the left bump are 32.8 and 52.7 respectively.

Download Full Size | PDF

4. Discussion and conclusions

In grating based X-ray phase contrast imaging, the most important task is to extract the monochromatic fringe phase shifts, which are critical to quantitative reconstruction of sample electron density distribution. In medical imaging and non-destructive testing, the grating setups are mostly implemented with polychromatic X-ray sources [4, 16, 17]. As an effort to develop an analytic approach to retrieve the monochromatic fringe phase shift from polychromatic measurement, we developed a simplified analytic approach [24]. That approach is applicable for samples with weak attenuation such as thin samples, or one has a prior knowledge of the spectral shift caused by sample attenuation. In practice, one likely does not have such prior knowledge. Moreover, samples usually have uneven thickness profiles. Uneven sample attenuation generates position-dependent x-ray spectrum shifts, which manifests itself as the beam hardening effects. The beam hardening causes gross distortion of fringe shifts. In this work, we address the beam hardening problem by developing a novel analytic approach. In this approach, we extract the monochromatic attenuation function A²(x, y; E) from the measured polychromatic attenuation by an iterative algorithm, provided that the uniform or single power energy dependence hold for the sample. We then compute the position-dependent effective spectrum Θ(E; x, y) by using A²(x, y; E). Finally, using the spectrum Θ(E; x, y) to iteratively solve the monochromatic fringe phase shift from the measured polychromatic fringe phase shift. Thus retrieved monochromatic fringe shift eventually gets rid of the beam hardening problem. This analytic approach is validated with numerical simulations with several imaging setups. Different from the laborious energy calibration reported in literature [25–27], which is applicable only to single-material samples, our approach is valid also for inhomogeneous multi-material samples.

The analytic approach developed in this work has several limitations. First, the assumption of the single-power energy dependence of sample attenuation may not hold for high tube voltage setups with very wide spectrum. Second, even when this assumption holds up, there is lack of a formula to determine the optimal exponent value in the energy factor f (E) = E⁻ⁿ. The optimal exponent value depends on the source spectrum, sample compositions and sample thickness profile. In general, the optimal n-values goes larger in the cases with less penetrating samples, or lower tube voltages, because in these cases the photoelectric absorption makes larger contribution to x-ray attenuation. Fortunately, as we show in section 3.2, any value within the range of 2 to 3.5 all result in reasonably good accuracies in the cases with low tube voltages.

In conclusion, assuming the uniform or single power energy dependence of sample attenuation, in this work we developed an analytic approach to retrieve the monochromatic attenuation function and fringe phase shift from the measured polychromatic attenuation. The retrieved monochromatic fringe shift gets rid of the beam hardening distortion, provides quantitative distribution of the projected electron density gradient of the sample. This analytic approach is validated with numerical simulations with several imaging setups. This work provides a useful tool in quantitative differential phase contrast imaging.