1. Introduction

X-ray grating interferometry is a differential X-ray phase-contrast imaging technique that has attracted intensive research efforts in recent years. This is because this technique may provide highly sensitive means for imaging soft tissues and low-Z materials, hence it has many potential applications in medical imaging and material science [1–15]. X-ray grating interferometry is usually implemented as a Talbot-Lau interferometer [4, 15–17], which consists of an X-ray source, a source grating G₀, a phase grating G₁ and an imaging detector, as is schematically shown in Fig. 1. Briefly speaking, the source grating G₀ serves as an aperture mask to break the source spot into an array of narrow sources to enhance the fringes visibility through improved spatial coherence of X-ray illumination. The phase grating G₁ serves as a beam splitter to generate interference fringes, which is recorded by the imaging detector. In addition, one may place an absorbing grating in front of the detector, and this grating serves as a fringe analyzer, which will be used when the fringe period is too small to be resolved by the detector pixels [4,15].

 

Fig. 1 Schematic of an x-ray phase grating interferometer with a microfocus source.

Download Full Size | PDF

In grating interferometry, sample information is encoded in the perturbed interference fringes. Consider first the case of monochromatic X-ray. Without a sample the G₁ grating generates periodic interference fringe pattern I_E, which can be expressed as [18,19]:

One of the most important goals of the grating based phase contrast imaging is to provide quantitative maps {ρ_e(x, y, z)} of a sample’s electron densities through tomographic acquisitions. To fulfill this goal in X-ray interferometry, one needs to determine the fringe phase shifts in each of the angular projections. In practice, for π/2-phase gratings the dominant diffraction orders are the m = 1 order and its complex conjugate m = −1, and for π-phase gratings the m = ±2 orders dominate, owing to limited spatial coherence of the source assembly. To retrieve the fringe phase shifts ϕ₁(x, y; E) or ϕ₂(x, y; E) from the distorted fringe pattern I_E, one can use either the phase stepping method or the Fourier analysis method. In the phase stepping method, an absorbing grating, which has the same period as the fringe pattern, is placed in front of the detector. This grating scans the fringe in steps to generate phase stepping curve to extract the fringe phase shifts [4, 15]. On the other hand, when the fringe is resolvable by the detector, one can apply the Fourier analysis method for phase retrieval. In this method one crops the fringe’s Fourier spectrum around a diffraction peak m/M_g p₁ to extract the fringe phase shifts [13,14,16,17]. Once the fringe phase shift ϕ_m is retrieved, the gradients ∂ρ_e,p/∂x of the sample’s projected electron densities can be simply recovered from ϕ_m by simply inverting Eq. (3). Repeating this retrieval procedure for all angular projections, one will be able to compute the quantitative volumetric distribution of sample’s electron densities by using the tomographic reconstruction method [6–8,20]. Hence the retrieval of fringe phase shifts is a task crucial to the grating based phase contrast imaging.

However, phase contrast imaging techniques are mostly implemented with polychromatic sources [4, 20–24]. Under polychromatic x-ray exposure, the intensity fringe pattern I_Poly is a sum of the fringe patterns {I_E} of different photon energies. We denote the fringe phase shift retrieved from the polychromatic intensity fringe I_Poly by ϕ_m,Poly(x, y). As such, one important question arises: how to retrieve the monochromatic fringe shifts from the measured polychromatic fringe shifts, as Eq. (3) holds only for the monochromatic case. To answer this question, one faces two challenges. First, as is shown Eq. (2), each monochromatic fringe shift ϕ_m is the phase angle of a complex vector representing the m-th diffraction order. The phase and amplitude of the complex vector are all energy-dependent. The polychromatic fringe shift is the phase angle of spectrum-weighted sum of those complex vectors of different energies. Hence, the first challenge is how to reconstruct monochromatic fringe shifts from the polychromatic fringe shifts generated by a beam with known effective X-ray spectrum. [21, 22]. The second challenge is how to deal with the beam hardening problem associated with uneven sample attenuation. The sample preferentially filters out the low-energy photons through sample attenuation, which make the effective spectrum shifted for higher energies. This is the well-known beam hardening effects of sample attenuation. But uneven sample attenuation makes the spectral shifts position-dependent. Consequently, uneven sample attenuation renders effective spectrum unknown for the monochromatic fringe shift retrieval. In this work we only address the first challenge, and leave the beam hardening problem for future research.

In order to address the first challenge, an energy calibration approach based on the effective energy concept was proposed [21, 22]. This approach hypothesizes that the polychromatic and monochromatic fringe phase shifts are proportional to each other and they have the same signs. Under this assumption polychromatic fringe phase shift ϕ_m,Poly(x, y) would be equivalent to a monochromatic fringe phase shift at an effective energy, which is denoted by E_eff [20–24]. Under this assumption one needs to conduct the energy calibration experiments to determine the effective energy that satisfies ϕ_m(x, y; E_eff) = ϕ_m,Poly(x, y). In the energy calibration experiments one employs phantoms of reference materials with known shape [20, 23, 24], and one must repeat the laborious calibration process whenever the interferometer setup or X-ray spectrum are modified. Moreover, the energy calibration approach is valid only for small fringe phase shifts of few degrees, as for larger fringe shifts the relationship between ϕ_m(x, y; E) and ϕ_m,Poly(x, y) is nonlinear and the concept of effective energy becomes invalid [21, 22]. For the cases with large fringe phase shifts researchers have to resort to laborious numerical simulations to analyze the relationship between the polychromatic and monochromatic fringe shifts [21,22].

In this work, we develop a novel analytic method to retrieve the monochromatic fringe phase shifts, consequently the sample gradients ∂ρ_e,p(x, y)/∂x, directly from the measured polychromatic fringe phase shifts ϕ_m,Poly(x, y), assuming the effective x-ray spectrum is known and available throughout this paper. In section 2, starting from the Fourier expansion of the distorted fringe pattern I_Poly, we present two analytic methods for computing the gradients ∂ρ_e,p(x, y)/∂x from the polychromatic fringe phase shift ϕ_m,Poly(x, y). In section 3 we apply these methods to derive the phase retrieval formulas for interferometers with π-phase gratings and π/2-phase gratings respectively. These analytic formulas are validated by the numerical simulation performed. We discuss the limitations of our method and conclude the paper in section 4.

2. Method

In order to recover the monochromatic fringe phase shifts ϕ_m(x, y; E) from the polychromatic fringe phase shifts ϕ_m,Poly(x, y), we developed an analytic approach as follows. Consider a grating interferometer with a polychromatic X-ray source and an energy-integration detector. We assume a known X-ray spectrum, which incorporates the source spectrum, the detector response and spectrum shift owing to X-ray attenuation of the sample. In this work we assume that the sample attenuation generates the same spectral shift for all detector pixels, so the attenuation effect is accounted for by the effective spectrum, which is denoted by $\mathfrak{D}\left(\mathrm{E}\right)$and is normalized such that ${\displaystyle \int \mathfrak{D}\left(\mathrm{E}\right)\text{dE}=1}$. In following discussion, we always assume that the effective spectrum is given. Our task is to study how to reconstruct monochromatic fringe shifts from the polychromatic fringe shifts generated with a given x-ray spectrum.

In this polychromatic case the detected fringe intensity I_Poly(M_gx, M_gy) is the spectrum-weighted sum of monochromatic fringe patterns of different energies, that is, $I_{\text{Poly}}\left(M_{g}x,M_{g}y\right)={\displaystyle \int \mathfrak{D}\left(\mathrm{E}\right)}\cdot I_E\left(M_{g}x,M_{g}y\right)\text{dE}$. Substituting I_E(M_gx, M_gy) of Eq. (2) into I_Poly(M_gx, M_gy), and including the sample attenuation in the effective spectrum $\mathfrak{D}\left(\mathrm{E}\right)$, we found that

2.1. First order approximations

To solve for ϕ_m,D(x, y) from Eq. (7), we keep only the linear terms in ϕ_m,D(x, y), and ignore the contributions of the higher powers of ϕ_m,D(x, y) in Eq. (7). We call this approximation the 1^st order approximation, and denote this approximate solution by ${\varphi }_{m,\mathrm{D}}^{\left(1\right)}\left(x,y\right)$. This would be a good approximation in the cases with |ϕ_m,D| ≪ 1. Under this first approximation, we found:

2.2. An iterative method for higher order approximations

In practice sometimes fringe phase shifts can be as large as several tens degree. In these cases, the linear solution (10) becomes inaccurate, and one should include higher-degree powers of ϕ_m,D(x, y) in the search for solution from Eq. (7). We developed an iterative method to solve Eq. (7) for ϕ_m,D(x, y) as follows. Assume ${\varphi }_{m,\mathrm{D}}^{\left(q\right)}\left(x,y\right)$, q ≥0, is the solution after the q-th iteration, then the updated solution of ϕ_m,D(x, y) in the (q + 1)-th iteration is given by

3. Results

As we mentioned earlier, in practice the dominant diffraction orders in the fringe patterns are the m = 1 order and its complex conjugate, for π/2-phase gratings, and the m = 2 order as well as its complex conjugate, for π-phase gratings. Hence one only needs to solve Eq. (7) for ϕ_1,D(x, y), the fringe phase shift of the m = 1 order at the design energy, from the polychromatic fringe shift ϕ_1,Poly(x, y) for interferometers with π/2-phase gratings, and solve Eq. (7) for ϕ_2,D(x, y) from ϕ_2,Poly(x, y) for interferometers with π-phase gratings.

3.1. For π-grating interferometers

Assume that for achieving high fringe visibility one sets the grating-to-detector distance the j-th fractional Talbot distance at the design energy E_D. In other words, one sets $R_2/M_g=j\cdot \left(p_1^2/8{\lambda }_{\mathrm{D}}\right)$, where λ_D = c·h/E_D is the X-ray wavelength at the design energy, and j = 1, 3, 5, ⋯, which denotes the order of the Talbot distances. Note that the phase shift Δϕ of the grating varies with photon energy as Δϕ(E) = (E_D/E) Δϕ(E_D), provided in the energy range there is no absorption edges. As such, substituting the fringe visibility coefficient C₂(E) of Eq. (6) into Eq. (8), we found the optical response coefficients Q₂(n) of this π-grating setup as follows:

To validate the above formulas for extracting the monochromatic fringe shifts with a π-phase grating interferometer, we conducted numerical simulations of polychromatic X-ray interferometry with phantoms of known electron distributions. From the known electron densities of phantoms we calculated the theoretical values of the monochromatic fringe shifts ${\varphi }_{2,\mathrm{D}}^{\text{Theory}}\left(x,y\right)$ at the design energy. On the other hand, we determined the polychromatic fringe phase shifts ϕ_2,Poly(x, y) from the fringe patterns generated by the Fresnel diffraction simulations [25, 26]. We applied the derived formulas of Eq. (13) or Eq. (14) to retrieve ${\varphi }_{2,\mathrm{D}}^{\left(1\right)}\left(x,y\right)$ and ${\varphi }_{2,\mathrm{D}}^{\left(q\right)}\left(x,y\right)$from the polychromatic fringe phase shifts ϕ_2,Poly(x, y). The determined ${\varphi }_{2,\mathrm{D}}^{\left(1\right)}$ and ${\varphi }_{2,\mathrm{D}}^{\left(q\right)}$ will be compared to the theoretical values of ${\varphi }_{2,\mathrm{D}}^{\text{Theory}}\left(x,y\right)$. Specifically, Fig. 2 shows the shape-profile of a one-dimensional phantom of glandular tissues. Each of the numerical one-dimensional phantoms consists of three leveled line sections and a pair of scant line sections of opposite orientations, as is shown in Fig. 2. Different phantoms differ from each other in the slopes of the scant line sections. Under the assumption of large source-to-grating distance, the pair of scant line sections generate the projected electron density-gradients of equal magnitude but opposite signs, while the leveled line sections in the phantom generate zero-gradients. Using Eq. (3), we computed the theoretical values ${\varphi }_{2,\mathrm{D}}^{\text{Theory}}\left(x,y\right)$ of the monochromatic fringe shifts at the design energy. On the other hand, we retrieved the polychromatic fringe shifts ϕ_2,Poly(x, y) from the fringe patterns generated by simulated Fresnel diffraction. Note that the values of ϕ_2,Poly(x, y) depends on the phantom, the effective X-ray spectrum and the interferometer setup employed. But ϕ_2,Poly(x, y) does not dependent on the focal spot size and shape, as the coherence degree is independent of photon energy, as is shown in Eq. (4). In the simulations, we assumed a 35 kVp x-ray effective spectrum as is shown in Fig. 3.

 

Fig. 2 Profile of phantom thickness.

Download Full Size | PDF

 

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

Download Full Size | PDF

To validate the effectiveness of Eqs. (13) and (14), we list the changes of the coefficients $\frac{1}{n!}\frac{Q_2\left(n\right)}{Q_2\left(1\right)}$ with respect to n in Table 1. It can be seen that, with the assumed 35 kVp x-ray effective spectrum, the coefficients $\frac{1}{n!}\frac{Q_2\left(n\right)}{Q_2\left(1\right)}$ decrease rapidly with n for the first and third Talbot distances. The convergence of Eq. (14) can be guaranteed when ϕ_2,D(x, y)≤ 90°.

Table 1. Decay of the coefficients $\frac{1}{n!}\frac{Q_m(n)}{Q_m(1)}$ with increasing n, in the iterative Eqs. (14) and (17). Here the effective spectrum is assumed 35 kVp shown in Fig. 3, and the Talbot distance was set to the first and third order (j = 1, 3) Talbot distances respectively.

View Table

Without loss of generality in the study of ϕ_2,Poly(x, y), we assumed a point-like focal spot in the simulation. The interferometer setup parameters selected in the simulations are as follows. The period of the π-grating was set to 5μm with design energies varying from 15 keV to 28 keV. The reduced grating-to-detector distance R₂/M_g (with M_g = 1) was set to the 1^st and 3^rd Talbot distances at the design energy. The polychromatic fringe phase shifts ϕ_2,Poly(x, y) was retrieved from the simulated fringe pattern by using the Fourier analysis method, as reported in [13, 14, 16,17]. Based on the simulated ϕ_2,Poly(x, y) data, we determine the monochromatic fringe shift ${\varphi }_{2,\mathrm{D}}^{\left(1\right)}\left(x,y\right)$ at the design energy by using Eq. (13) as the 1^st order approximation, or using the iterative method of Eq. (14) for higher order approximations ${\varphi }_{2,\mathrm{D}}^{\left(q\right)}\left(x,y\right)$. We then checked the accuracies of thus determined solutions ${\varphi }_{2,\mathrm{D}}^{\left(1\right)}\left(x,y\right)$ and ${\varphi }_{2,\mathrm{D}}^{\left(q\right)}\left(x,y\right)$against the theoretical values ${\varphi }_{2,\mathrm{D}}^{\text{Theory}}\left(x,y\right)$ of the monochromatic fringe shifts.

Figure 4 plots the simulation results for a π-grating interferometer setup, in which the grating design energy E_D was set to 18 keV and the grating-detector distance was set to the 1^st fractional Talbot distance, i.e., $R_2/M_g=p_1^2/8{\lambda }_{\mathrm{D}}$, where λ_D is the X-ray wavelength at the design energy. The theoretical values ${\varphi }_{2,\mathrm{D}}^{\text{Theory}}(x,y)$ of the monochromatic fringe shifts for this phantom is given by the blue profile, which traces zeros for the three leveled line sections, and 5 degrees for one slope and -5 degrees for the other. The green profile represents the values of the polychromatic fringe shifts generated with the effective x-ray spectrum given in Fig. 3. Figure 4 shows that the polychromatic fringe shift ϕ_2,Poly(x, y) on the two slopes differs from ${\varphi }_{2,\mathrm{D}}^{\text{Theory}}(x,y)$ by relative errors of 11% or less. The cyan trace depicts the values of the monochromatic fringe shift ${\varphi }_{2,\mathrm{D}}^{(1)}(x,y)$ computed with the 1^st order approximation as is summarized by Eq. (13). Thus determined values of ${\varphi }_{2,\mathrm{D}}^{(1)}(x,y)$ differ from ${\varphi }_{2,\mathrm{D}}^{\text{Theory}}(x,y)$ by only 1.1%. The red trace in Fig. 4 depicts the profile of ϕ_2,D(x, y) computed with the iterative method of Eq. (14) after 8 iterations. The red profile shows that thus determined ${\varphi }_{2,\mathrm{D}}^{(q)}(x,y)$ values are almost identical to ${\varphi }_{2,\mathrm{D}}^{\text{Theory}}(x,y)$ with relative errors of 0.9% or less. Summarizing the results shown in Fig. 4, we see that when the magnitudes of polychromatic fringe shifts ϕ_2,Poly(x, y) are only few degrees, the 1^st approximation (Eq. (13)) is good enough to determines the monochromatic fringe shift ϕ_2,D(x, y) at the design energy from the measured polychromatic fringe shifts ϕ_2,Poly(x, y).

 

Fig. 4 Plots of the simulation results: The design energy is set to 18 keV with π-grating interferometer. In this simulation, the true theoretical values of ${\varphi }_{2,\mathrm{D}}^{\text{Theory}}(x,y)$ (the solid blue line) are set to a relatively small value of 5 degrees. The dash-dot green line, ϕ_2,Poly(x, y), is the polychromatic fringe phase shifts retrieved from the simulated fringe pattern by using the Fourier analysis method. The dotted cyan line is the 1^st order approximate solution ${\varphi }_{2,\mathrm{D}}^{(1)}(x,y)$ using Eq. (13). The dashed red line is the higher order approximate solution ${\varphi }_{2,\mathrm{D}}^{(q)}(x,y)$ using Eq. (14).

Download Full Size | PDF

However, when polychromatic fringe shifts ϕ_2,Poly(x, y) are large, the 1^st approximation of Eq. (13) is not good enough to determine monochromatic fringe shifts ϕ_2,D(x, y) at the design energy. Figure 5 shows the results with the same interferometer settings but with a different phantom that generates polychromatic fringe shifts ϕ_2,Poly(x, y) as large as ±39.5 degrees, as is shown by the green profile in Fig. 5. The theoretical values ${\varphi }_{2,\mathrm{D}}^{\text{Theory}}(x,y)$ of the monochromatic fringe shifts for this second phantom is given by the blue profile in Fig. 5, which traces zeros for the three leveled line sections, and ±45 degrees for the two slopes of the phantom. Figure 5 shows that in this simulation the polychromatic fringe shift ϕ_2,Poly(x, y) on the two slopes differs from ${\varphi }_{2,\mathrm{D}}^{\text{Theory}}(x,y)$ by relative errors of 12.1% or less. The cyan trace depicts the values of the monochromatic fringe shift ${\varphi }_{2,\mathrm{D}}^{(1)}(x,y)$ computed with the 1^st order approximation of Eq. (13). Figure 5 shows that ${\varphi }_{2,\mathrm{D}}^{(1)}(x,y)$ overshoots the theoretical values by errors as large as 19.1%. This circumstance indicates that the 1^st approximation of Eq. (13) is not accurate enough, and one should use the iterative solution of Eq. (14). The red trace in Fig. 5 depicts the profile of the iterative solution ${\varphi }_{2,\mathrm{D}}^{(q)}(x,y)$ after 28 iterations in Eq. (14). It shows that thus determined ${\varphi }_{2,\mathrm{D}}^{(q)}(x,y)$ values are almost identical to ${\varphi }_{2,\mathrm{D}}^{\text{Theory}}(x,y)$ with relative errors of 0.1% or less. This results demonstrates clearly the advantage of the iterative method of Eq. (14) as compared to the 1^st approximation method of Eq. (13).

 

Fig. 5 Plots of the simulation results: The design energy is set to 18 keV with π-grating interferometer. In this simulation, the true theoretical values of ${\varphi }_{2,\mathrm{D}}^{\text{Theory}}(x,y)$ (the solid blue line) are set to a large value of 45 degress. The dash-dot green line, ϕ_2,Poly(x, y), is the polychromatic fringe phase shifts retrieved from the simulated fringe pattern by using the Fourier analysis method. The dotted cyan line is the 1^st order approximate solution ${\varphi }_{2,\mathrm{D}}^{(1)}(x,y)$ using Eq. (13). The dashed red line is the higher order approximate solution ${\varphi }_{2,\mathrm{D}}^{(q)}(x,y)$ using Eq. (14).

Download Full Size | PDF

To evaluate the performance of Eqs. (13) and (14) with noisy data, we added Gaussian white noise (20%) to the fringe intensity image I_poly and the grating only image I_g. We then retrieved the polychromatic fringe phase shifts ϕ_2,Poly(x, y) from the simulated noisy fringe pattern using the Fourier analysis method. The 1^st order approximation solution ${\varphi }_{2,\mathrm{D}}^{(1)}(x,y)$ as well as the higher order solution ${\varphi }_{2,\mathrm{D}}^{(q)}(x,y)$, were computed using Eqs. (13) and (14). The retrieved fringe phases attain good signal noise ratios. For example, for the left bump as is shown in Fig. 6, ϕ_2,Poly(x, y), ${\varphi }_{2,\mathrm{D}}^{(1)}(x,y)$, and ${\varphi }_{2,\mathrm{D}}^{(q)}(x,y)$ have SNRs of 64.3, 45.7, and 63.8 respectively. It shows that our methods are robust against noise.

 

Fig. 6 Plots of the simulation results: The design energy is set to 18 keV with π-grating interferometer. In this simulation, the true theoretical values of ${\varphi }_{2,\mathrm{D}}^{\text{Theory}}(x,y)$ (the solid blue line) are set to 45 degress. Gaussian white noise (20%) was added to the fringe intensity image I_poly and the grating only image I_g. In the plot, The dash-dot green line, ϕ_2,Poly(x, y), is the polychromatic fringe phase shifts retrieved from the simulated noisy fringe pattern by using the Fourier analysis method. The dotted cyan line is the 1^st order approximate solution ${\varphi }_{2,\mathrm{D}}^{(1)}(x,y)$ using Eq. (13). The dashed red line is the higher order approximate solution ${\varphi }_{2,\mathrm{D}}^{(q)}(x,y)$ using Eq. (14). The signal noise ratios (SNRs) at the left bump for ϕ_2,Poly(x, y), ${\varphi }_{2,\mathrm{D}}^{(1)}(x,y)$, and ${\varphi }_{2,\mathrm{D}}^{(q)}(x,y)$ are 64.3, 45.7, and 63.8 respectively.

Download Full Size | PDF

We have also tested these two methods for the interferometer settings with higher order fractional Talbot distances such as the third or fifth fractional Talbot distances. We observed similar results as that presented above in details. But as we will show below, π/2-phase grating interferometers may exhibit a different picture.

3.2. For π/2-grating interferometers

Different from π-phase grating interferometers, the dominant diffraction order of a π/2-phase grating interferometer is the m = 1 order and its complex conjugate. Setting the grating-to-detector distance to one of the fractional Talbot distance at the design energy E_D, one has $R_2/M_g=j\cdot \left(p_1^2/2{\lambda }_{\mathrm{D}}\right)$, where j = 1, 3, 5,⋯, which denotes the order of the Talbot distances. Substituting the fringe visibility coefficient C₁(E) of Eq. (6) into Eq. (8), we found that the optical response coefficients Q₁(n) for this π/2-phase grating setup are given by:

The changes of the coefficients $\frac{1}{n!}\frac{Q_1(n)}{Q_1(1)}$ with respect to n are shown in Table 1. It can be seen that the coefficients $\frac{1}{n!}\frac{Q_1(n)}{Q_1(1)}$ decrease rapidly with n. So Eqs. (16) and (17) is suitable in seeking the solution of the gradient of the electron density at design energy.

At the first glance the solution Eqs. (16)–(17) look similar to the solution Eqs. (13)–(14) for π-grating setups. However, there is one important difference between them. For π-grating setups, the ratio Q₂(0)/Q₂(1) in Eqs. (13)–(14) is always positive, which means that the polychromatic and monochromatic fringe shifts will always have the same sign (for fringe shift value less than π/2). But, for π/2-grating setups, the correction factor Q₁(0)/Q₁(1) in Eqs. (16)–(17) can be either positive or negative, depending on the Talbot distance order, the spectrum and the design energy of the setup. For example, if the 1^st order fractional Talbot distance is selected (j = 1 in Eq. (16)), then all the integrands in Eq. (16) is positive, so the ratio Q₁(0)/Q₁(1) holds positive regardless of the spectrum and grating design energy. But if the 3^rd (j = 3), or 5^th (j = 5) order Talbot distance is selected, the ratio Q₁(0)/Q₁(1) may turn to negative. Hence, for π/2-grating setups, the polychromatic and monochromatic fringe shifts can differ not only in magnitudes, but also in the signs. Hence Eqs. (16)–(17) should be used to recover the true monochromatic fringe phase shifts for π/2-grating setups.

The following simulation results illustrated this unusual feature of π/2-grating setups. Figure 7 shows the simulation results for a π/2-grating interferometer setup, in which the grating design energy E_D = 26.5 keV and the reduced grating-to-detector distance was set to the 3^rd order fractional Talbot distance, i.e., $R_2/M_g=3\times p_1^2/2{\lambda }_{\mathrm{D}}$. Usually, the polychromatic fringe shifts are of the same sign as the theoretical values ${\varphi }_{1,\mathrm{D}}^{\text{Theory}}(x,y)$. But along the left bump in Fig. 7, the polychromatic fringe shift ϕ_1,Poly (green profile) register a value of −2.85 degree, while the theoretical value of the monochromatic fringe shift ${\varphi }_{1,\mathrm{D}}^{\text{Theory}}(x,y)$ (blue profile) is 5 degree. Along the right bump in Fig. 7, the green and blue profiles show that ϕ_1,Poly = 2.85 degree and ${\varphi }_{1,\mathrm{D}}^{\text{Theory}}(x,y)=-5$ degree. Under this peticuliar circumstance, if one blindly recovers the gradients ∂ρ_e,p (x, y)/∂x directly from the polychromatic fringe shift ϕ_1,Poly, then originally positive projected electron density gradients could be mistaken as negative gradients. Fortunately, Eq. (16) or Eq. (17) can be used to recover the true monochromatic fringe phase shifts. In this example, the iterative solution ${\varphi }_{1,\mathrm{D}}^{(q)}(x,y)$, i.e., the monochromatic fringe shifts computed with the iterative method, recovers the theoretical values ${\varphi }_{1,\mathrm{D}}^{\text{Theory}}(x,y)$, as is shown by the overlapped red and blue profiles in Fig. 7.

 

Fig. 7 Plots of the simulation results: The design energy is set to 26.5 keV with π/2-grating interferometer. In this simulation, the true theoretical values of ${\varphi }_{1,\mathrm{D}}^{\text{Theory}}(x,y)$ (the solid blue line) are set to 5 degrees. The dash-dot green line, ϕ_1,Poly(x, y), is the polychromatic fringe phase shifts retrieved from the simulated fringe pattern by using the Fourier analysis method. The dotted cyan line is the 1^st order approximate solution ${\varphi }_{1,\mathrm{D}}^{(1)}(x,y)$ using Eq. (16). The dashed red line is the higher order approximate solution ${\varphi }_{1,\mathrm{D}}^{(q)}(x,y)$ using Eq. (17).

Download Full Size | PDF

Figure 8 shows the simulation results for a π/2-grating interferometer setup, which was set to the 1^st fractional Talbot distance, i.e., $R_2/M_g=p_1^2/2{\lambda }_{\mathrm{D}}$, and the grating design energy is also set to E_D = 18 keV. The green, blue, cyan and red profiles corresponding to the same quantities as in Figs. 4 and 5 for π-grating interferometers. Figure 8 shows that the polychromatic fringe shifts ϕ_1,Poly(x, y) take values between ±4.4 degrees. The green profile (ϕ_1,Poly) deviates from the blue profile of the theoretical values $\left({\varphi }_{1,\mathrm{D}}^{\text{Theory}}\right)$ by upto 11.9%. Figure 8 also shows that the computed values of the monochromatic fringe shift, by using either the 1^st order approximation of Eq. (16), or the iterative solution of Eq. (17), are almost identical to the theoretical values $\left({\varphi }_{1,\mathrm{D}}^{\text{Theory}}\right)$ with relative errors of 0.23% and 0.04% or less respectively.

 

Fig. 8 Plots of the simulation results: The design energy is set to 18 keV with π/2-grating interferometer. In this simulation, the true theoretical values of ${\varphi }_{1,\mathrm{D}}^{\text{Theory}}\left(x,y\right)$ (the solid blue line) are set to 5 degrees. The dash-dot green line, ϕ_2,Poly(x, y), is the polychromatic fringe phase shifts retrieved from the simulated fringe pattern by using the Fourier analysis method. The dotted cyan line is the 1^st order approximate solution ${\varphi }_{1,\mathrm{D}}^{\left(1\right)}\left(x,y\right)$ using Eq. (16). The dashed red line is the higher order approximate solution ${\varphi }_{1,\mathrm{D}}^{\left(q\right)}\left(x,y\right)$ using Eq. (17).

Download Full Size | PDF

4. Discussion and conclusions

As is well known, one of the most important goals of X-ray phase contrast imaging is to provide quantitative maps of the sample’s electron densities. To fulfill this goal in the grating based interferometry, one needs to determine the fringe phase shift fringes in angular projections. In the monochromatic case the gradient of sample’s projected electron densities, or the sample density gradients for short, can be calculated from the measured fringe phase shift, since they are simply proportional to each other, and the proportional constant is unambiguously defined, as is shown in Eq. (3). However, in the polychromatic case, as is shown in previous sections, the relationship between the sample’s density gradient and fringe shift becomes much more complicated. In the literature, researchers proposed the energy calibration approach based on the effective energy concept [21–23], as we discussed in section 1. This approach is valid only if the polychromatic and monochromatic fringe phase shifts are of few degrees and both have the same sign. In fact, if the polychromatic and monochromatic fringe phase shifts have opposite signs, no effective energy can be defined. In this work we proposed a general analytic approach, assuming the effective spectrum is known. Specifically, we present the 1^st order approximation method and the iterative method in our analytic approach. The derived formulas, which are summarized in Eqs. (13)–(14) for π-grating interferometer setups, and in Eqs. (16)–(17) for π/2-grating setups, enable us to compute the monochromatic fringe phase shifts and the sample density gradients from polychromatic fringe phase shifts as large as few tens of degrees. This approach is validated with numerical simulations of several grating setups.

The nonlinear relationship between the monochromatic and polychromatic fringe phase shifts, as is shown in previous section, has its roots in the dependence of the fringe visibility on photon energy. For example, for π/2-grating interferometer setups, the fringe visibility coefficient C₁(E) ∝ sin (πE_D/(2E)) × sin (jπE_D/(2E)), as is shown in Eq. (15). For the setups with high order Talbot distances, the product sin (πE_D/(2E)) × sin (jπE_D/(2E)) oscillates between positive and negative as photon energy varies. This causes the intensity modulation varies with energy not only in its magnitude, but also in polarity. Under polychromatic x-ray illumination, this C₁(E) oscillation with energies has two consequences. First, it significantly reduces the fringe visibility because of the cancellation of intensity modulation among the intensities generated by photons of different energies. Such a phenomenon of the fringe visibility loss was observed in experiments [27]. Second, this C₁(E) oscillation with energies give rise to the negative Q₁(0)/Q₁(1) ratios, which consequently makes the polychromatic and monochromatic fringe shifts to have opposite signs. This peculiar feature vividly demonstrates the crucial role that x-ray spectrum plays in setting up the optical properties of a grating interferometer.

Several limitations of this work should be mentioned. First, we found that the performance of the iterative methods of Eq. (14) and Eq. (17) will degraded if the polychromatic fringe phase shifts are as large as 80 to 90 degrees. This is because more higher order terms should be considered in the iterations, and the solution updates get stagnated or blow out due to numerical round-of error. Fortunately, it is rare to encounter such large fringe shifts. Second, in this work we assumed the knowledge of the effective X-ray spectrum, which takes into account also the spectrum shift caused by X-ray attenuation of the sample. This assumption implies that the sample attenuation is uniform, and one should have a prior knowledge of its functional relation with photon energy. But in practice one encounters uneven sample attenuation, which causes spectrum shifts that are varying over detector pixels. Such uneven beam hardening effects have caused gross inaccuracies in the energy calibration approach [22, 24]. Our analytic approach cannot deal with this kind of beam hardening problem, and more research is needed.

In conclusion, assuming the knowledge of the effective X-ray spectrum, we proposed in this work a general analytic approach that enables one to compute the monochromatic fringe phase shifts and the sample density gradients from polychromatic fringe phase shifts as large as few tens of degrees. This approach, which includes the first approximation method and the iterative method, explores the non-linear relationship between the polychromatic and monochromatic fringe phase shifts. This approach is validated with numerical simulations of several grating interferometry setups. Hence, the analytic approach developed in this work provides a useful tool in quantitative imaging based on the polychromatic grating interferometry.