Polychromatic X-ray effects on fringe phase shifts in grating interferometry

: In order to quantitatively determine the projected electron densities of a sample, one needs to extract the monochromatic fringe phase shifts from the polychromatic fringe phase shifts measured in the grating interferometry with incoherent X-ray sources. In this work the authors propose a novel analytic approach that allows to directly compute the monochromatic fringe shifts from the polychromatic fringe shifts. This approach is validated with numerical simulations of several grating interferometry setups. This work provides a useful tool in quantitative imaging for biomedical and material science applications.


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][2][3][4][5][6][7][8][9][10][11][12][13][14][15]. X-ray grating interferometry is usually implemented as a Talbot-Lau interferometer [4,[15][16][17], which consists of an X-ray source, a source grating G 0 , a phase grating G 1 and an imaging detector, as is schematically shown in Fig. 1. Briefly speaking, the source grating G 0 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 1 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].
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 1 grating generates periodic interference fringe pattern I E , which can be expressed as [18,19]: In this expansion, the diffraction orders are labeled by integer m, and C m (E) denotes the Fourier coefficient associated with the harmonics exp i2πmx/p 1 , and C m (E) is a function of the photon energy E. Note that C m (E) determines the modulation of the m-th order fringes. For this reason we will call C m (E) the fringe visibility coefficient of the m-th order fringes throughout this paper. In Eq. (1) γ(m) denotes the reduced coherence degree associated with the m-th diffraction order. The coherence degree is determined by the source and source-grating G 0 assembly, and it also affects the visibility of the interference fringes [4,[15][16][17]. In addition, in Eq. (1) I in is the entrance intensity at G 1 plane, and M g = (R 1 + R 2 )/R 1 denotes the geometric magnification factor, where R 1 and R 2 are the G 0 -to-G 1 and G 1 -to-detector distances respectively. In the presence of a sample, the intensity fringe is distorted and the sample's information is encoded into the fringe pattern. Specifically, the intensity fringe pattern I E is modified to [4,16,17]: It indicates that each of the diffraction orders of the fringe gets a position-dependent phase shift, which we call the fringe phase shift throughout the paper. In Eq. (2) we denote the fringe phase shift of the m-th diffraction order by φ m (x, y; E), which is related to the sample's electron densities by the following expression: where p 1 denotes the period of the phase grating G 1 , E is the photon energy, h the Planck's constant, c the speed of light, and r e the classical electron radius. In Eq. (3) ρ e, p (x, y) . = Ray ρ e (x, y, s) ds denotes the sample's projected electron density, that is, the integration of the sample's electron density ρ e (x, y, z) over the ray path. Equation (3) shows that the fringe phase shifts φ m (x, y; E) are proportional to the projected electron density gradients ∂ρ e, p (x, y)/∂x of the sample, and inversely proportional to the squared photon energy E 2 . Note that this energy dependence of fringe shifts holds provided no absorption edge appears in the photon-energy range studied. In addition, in Eq. (2) A 2 denotes the sample's attenuation, which is energy dependent as well.
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 φ 1 (x, y; E) or φ 2 (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 1 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][7][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][21][22][23][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][21][22][23][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.

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 D(E) and is normalized such that D(E) 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 g x, M g y) is the spectrum-weighted sum of monochromatic fringe patterns of different energies, that is, , and including the sample attenuation in the effective spectrum D(E), we found that where I in denotes the entrance detected energy fluence. Note that for incoherent polychromatic X-ray sources such as X-ray tubes, the reduced coherence degree γ(m) is independent of Xray photon energy. In Eq. (4) we write φ m (x, y; E) in terms of (E 2 D /E 2 )φ m,D (x, y), where E D is the design energy of the phase grating G 1 , and φ m,D (x, y) denotes the fringe phase shift at energy E D . This reflects that the monochromatic fringe phase shift generated by a sample is inversely proportional to the square of photon energy, providing the photon energy is away from the absorption edges of the sample elements. Under polychromatic X-ray illumination, one may retrieve the polychromatic fringe phase shift φ m,Poly (x, y) from the fringe pattern by using either the phase stepping method or the Fourier analysis method [4,[13][14][15][16][17]. We realized that, regardless of the method used, the retrieved polychromatic fringe phase shift φ m,Poly (x, y) is given by where the action of operator Arg{·} is to extract the phase angle of the expression in the bracket. For a given spectrum D(E), if one knows the fringe visibility coefficients C m (E), one should in principle be able to recover φ m,D (x, y) and thereby the sample's projected electron density gradients ∂ρ e, p (x, y)/∂x. We have derived the general formula for C m (E) in our previous works: [18, 19]: In this equation, Δφ is the phase-shift step of the phase grating. Since the floor-function x is defined as the largest integer that is less or equal to x, the factor (−1) 4kλ R 2 /M g p 2 1 swings between 1 and −1 as its exponent changes. Substituting Eq. (6) into Eq. (5) and completing the integral over photon energies, one could compute the fringe phase shift φ m,Poly (x, y) if we know φ m (x, y; E). In practice, substituting Eq. (6) into Eq. (5), we want to invert Eq. (5) to recover φ m,D (x, y) from the measured polychromatic fringe phase shift φ m,Poly (x, y) and consequently the gradients ∂ρ e, p (x, y)/∂x. Using Taylor expansion of exp , and completing the integrals over photon energies in Eq. (5), we found that: where the position-independent coefficients Q m (n) is defined as follows: These coefficients characterize the optical responses of the interferometer to a given X-ray spectrum. We will call {Q m (n)} the n-th optical response coefficient of a interferometer setup. Especially, the 0-th optical response coefficient Q m (0) represents the visibility of the m-th diffraction order of an interferometer setup. Equation (7) shows that in general the polychromatic fringe phase shift φ m,Poly (x, y) has a non-linear relationship with the monochromatic fringe phase shift φ m,D (x, y).

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 φ (1) m,D (x, y). This would be a good approximation in the cases with φ m,D 1. Under this first approximation, we found: If the measured polychromatic fringe shifts φ m,Poly (x, y) are few degrees only, then Eq. (9) shows that the monochromatic fringe phase shift at the design energy That is, when φ m,Poly 1 and φ m,D 1, φ m,D (x, y) and φ m,Poly (x, y) have a linear relationship, and φ m,D (x, y) is equal to φ m,Poly (x, y) multiplied by a position-independent correction factor Q m (0)/Q m (1). In the literature, this linear relationship between φ m,D (x, y) and φ m,Poly (x, y) has prompted the concept of the effective energy, which is defined as the photon energy that satisfying φ m (x, y; E = E eff ) = φ m,Poly (x, y), as we mentioned in section 1 [21,22,24]. Since the monochromatic fringe phase shift φ m (x, y; E) inversely proportional to the square of photon energy, hence Eq. (9) implies thereby E eff = √ Q m (0)/Q m (1) · E D , provided the correction factor Q m (0)/Q m (1) > 0. Then an important problem is how to compute this correction factor. In the literature, the correction factor Q m (0)/Q m (1) has been determined through numerical simulations of X-ray wave propagation [21], or through energy calibration experiments with reference phantoms of known shape and composition [20,23,24]. In contrast to these previous works, we will compute Q m (0)/Q m (1) directly by substituting the fringe visibility coefficients C m (E) of Eq. (6), into Eq. (8). Once (1) is computed for a given interferometer setup, the gradients ∂ρ e, p (x, y)/∂x of sample's projected electron densities can be found from the measured polychromatic fringe phase shift φ m,Poly (x, y) as follows:

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 that φ (q) m,D (x, y), 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 Using Eq. (11) and starting from q = 0, we get Poly (x, y) at the 1 st approximation. In the (q+1)-th iteration the contributions from the terms up to (2q+1)-th powers of φ m,D (x, y) are included in Eq. (11). The iteration stops when φ where > 0 is a designated small number representing the error allowance and ||·|| 2 indicates the l 2 norm. The smaller the is, the more accurate the solution is. The resulting φ (q+1) m,D (x, y) from the last iteration is deemed to the solution of Eq. (7).

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.

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 · p 2 1 /8λ D , 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 2 (E) of Eq. (6) into Eq. (8), we found the optical response coefficients Q 2 (n) of this π-grating setup as follows: n = 1, 2, 3, · · · . Equipped with Eq. (9) and Eq. (12), we found following relationship between the monochromatic and polychromatic fringe shifts in the 1 st approximation: Equation (13) shows that, for a given π-phase grating interferometer, the correction factor Q 2 (0)/Q 2 (1) depends on the effective x-ray spectrum D(E), the design energy E D , and the selected order j of the fractional Talbot distances. Similarly, the iterative solution Eq. (11) becomes: 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 φ  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. To validate the effectiveness of Eqs. (13) and (14), we list the changes of the coefficients 1 n! Q 2 (n) Q 2 (1) with respect to n in Table 1. It can be seen that, with the assumed 35 kVp x-ray effective spectrum, the coefficients 1 n! Q 2 (n) Q 2 (1) 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 • .
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 2 /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 φ (1) 2,D (x, y) 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 φ (q) 2,D (x, y). We then checked the 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.    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 2 1 /8λ D , where λ D is the X-ray wavelength at the design energy. The theoretical values φ Theory 2,D (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 φ  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). 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 φ Theory 2,D (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 φ    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 us- ing the Fourier analysis method. The 1 st order approximation solution φ (1) 2,D (x, y), as well as the higher order solution φ (q) 2,D (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), φ (1) 2,D (x, y), and φ 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.

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-todetector distance to one of the fractional Talbot distance at the design energy E D , one has R 2 /M g = j · p 2 1 /2λ D , where j = 1, 3, 5, · · · , which denotes the order of the Talbot distances. Substituting the fringe visibility coefficient C 1 (E) of Eq. (6) into Eq. (8), we found that the optical response coefficients Q 1 (n) for this π/2-phase grating setup are given by: Using Eq. (9) and Eq. (15), we derived the following relationship between the monochromatic and polychromatic fringe shifts in the 1 st approximation: Poly (x, y) , Obviously the iterative solution Eq. (11) becomes: The changes of the coefficients 1 n! Q 1 (n) Q 1 (1) with respect to n are shown in Table 1. It can be seen that the coefficients 1 n! 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 2 (0)/Q 2 (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 1 (0)/Q 1 (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 1 (0)/Q 1 (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 1 (0)/Q 1 (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 × p 2 1 /2λ D . Usually, the polychromatic fringe shifts are of the same sign as the theoretical values φ  Fig. 7. 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 2 1 /2λ 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 (φ Theory 1,D ) 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 (φ Theory 1,D ) with relative errors of 0.23% and 0.04% or less respectively.

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 π/2grating 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 1 (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 1 (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 1 (E) oscillation with energies give rise to the negative Q 1 (0)/Q 1 (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.