Beam hardening correction in polychromatic x-ray grating interferometry.

The beam hardening is one of the two causes of the fringe shift distortion in polychromatic X-ray grating interferometry. Based on the assumption of the uniform energy dependence, we developed a novel analytic approach to accurately retrieve the monochromatic attenuation function and fringe phase shift from the polychromatic measurement. This approach provides a useful tool for precise measurement of sample electron density distribution in X-ray grating interferometry.


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][2][3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18][19][20]. As is schematically shown in Fig. 1, this phase contrast imaging technique usually employs a Talbot-Lau interferometer [4,[20][21][22], which consists of an X-ray source, a source grating G 0 , a phase grating G 1 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][21][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 0 breaks up the source spot into an array of narrow sources to improve the spatial coherence of X-ray illumination. The phase grating G 1 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 g x, M g y; E) is given by the following expression: where M g = (R 1 + R 2 )/R 1 is the geometric magnification factor, R 1 is the G 0 -to-G 1 distance, and R 2 the G 1 -to-detector distance. In Eq. (1), each of the diffracted orders is labeled by an integer and represented by exp [i2πmx/p 1 ]. The amplitude of a given diffracted order is determined by the fringe-visibility coefficient C m (E), and the coherence degree γ(m). The amplitude is also modulated by the sample's attenuation function A 2 (x, y, E), which strongly depends on x-ray photon energy. Moreover, Eq. (1) indicates that each of the diffracted orders is encoded with a position-dependent fringe phase shift. For the m-th order, the fringe shift φ m (x, y; E) is given by [1][2][3][4]: where p 1 denotes the period of the phase grating, E is the x-ray photon energy, h the Planck's constant, c the speed of light, and r e the classical electron radius. In Eq. (2), ρ e, p (x, y) = ∫ Ray ρ e (x, y, s) ds denotes the sample's projected electron density, that is, the integration of sample's electron density ρ e (x, y, z) over the ray path z. The fringe shift is proportional to ∂ ρ e, p (x, y)/∂ x, which we call the density gradient thereafter. Equation (2) indicates that the fringe phase shift generated by a sample is inversely proportional to the square of photon energy. 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 g x, M g y; 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 φ 1 (x, y; E) in these cases. For π-grating setups, the m = ±2 orders dominate instead, so φ 2 (x, y; E) needs to be retrieved. In the monochromatic case, the fringe phase shifts φ 1 (x, y; E) or φ 2 (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][16][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: where D(E) denotes the normalized spectrum incorporating the source spectrum, the detector response and the phase grating substrate attenuation effect. One needs to measure the source spectrum using a spectrometer such as a CdTe x-ray spectrometer. The detector response, namely the quantum efficiency as a function of photon energy, can be calculated from the scintillator's chemical composition and the coating weight. Similarly, the attenuation of the grating substrate, as a function of photon energy, can be calculated from chemical composition, mass density and thickness of the grating substrate. In Eq. (3) we denote the fringe phase shift at the design energy E D by φ m,D (x, y). We have also made use of the relation φ m (x, y; E) = E 2 D E 2 φ m,D (x, y), provided the spectrum range does not include any absorption edge of the sample. In writing Eq. (3) note that the reduced coherence degree γ(m) associated with an X-ray tube is energy-independent. In experiments, one measures the polychromatic fringe phase-shift using either the phase stepping method or the Fourier analysis method [4-7, 13, 16, 17]. In either way, what one measured from the polychromatic fringe pattern is not the fringe shift φ m (x, y; E) at a given energy, rather one attains a distorted fringe shift that has a nonlinear relationship with the monochromatic fringe shift. In the following we denote the measured polychromatic fringe shift by φ m,Poly (x, y).
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 nonuniform 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][26][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][26][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][26][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.

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 2 (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 2 Poly (x, y) from the zero-th diffracted order of the fringe pattern [4]. Note that polychromatic attenuation A 2 Poly (x, y) is a spectral average is the spectrum. To extract the monochromatic attenuation function A 2 (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 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 2 (x, y; E), we solve Eq. (4) for p D (x, y) from the measured polychromatic attenuation A 2 Poly (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 and One can easily see that H(p(x, y)) is a strictly convex functional, since the function We want to retrieve the function p D (x, y), from given value A 2 Poly (x, y) = H(p D (x, y)) through iterative method: assuming p k (x, y) is acquired after the k-th iteration, and ∆(x, y) = A 2 Poly (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 2 Poly (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 Iter D (x, y). 1. Choose an initial starting value p 0 (x, y) = 0. 3. If the error ||δ(x, y)|| 2 is less than a preset threshold , we exit the iteration. Otherwise, set p k+1 (x, y) = p k (x, y) + δ(x, y) and go back to step 2. Here ||·|| 2 represent the standard 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)|| 2 would fall under 10 −13 after 15 iteration steps.
Once the sample attenuation function is extracted, we may incorporate this attenuation-filtration effect into a position-dependent effective 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 where the action of operator Arg {·} is to extract the phase angle of the expression in the bracket. This is the key equation for the retrieval of monochromatic fringe phase shifts. Owing to limited x-ray spatial coherence, for the setups with π/2-phase gratings, the dominant diffraction orders are m = ±1. Consequently, one only needs to measure just φ 1,Poly (x, y) from the fringe pattern, and then tries to retrieve the monochromatic fringe shift φ 1,D (x, y) based on Eq. (10). On the other hand, for the setups with π-phase gratings, the m = ±2 orders dominate, and one needs to measure φ 2,Poly (x, y) instead for retrieval of the monochromatic fringe shift φ 2,D (x, y). For a given effective spectrum, in order to solve for monochromatic fringe shift φ m,D (x, y) using Eq. (10), we need to know how to compute the fringe visibility coefficients C m (E). Fortunately we have already derived a general formula of C m (E) as follows [18]: In this formula ∆φ(E) is the phase shift step of the phase grating at energy E, and λ is x-ray wavelength. The factor (−1) 4kλR 2 /M g p 2 1 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 positiondependent 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 φ (q) m,D (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 φ < , where the operator ||·|| 2 indicates the L 2 -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). 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.

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 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 Iter D (x) by using the iterative algorithm of Eqs. (5)- (7). In this way we extracted the monochromatic attenuation function A 2 (x; E) = exp −p Iter D (x) · f (E)/ f (E D ) , 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 2 /M g (with M g = 1) was set to the 1 st Talbot distance R 2 /M g = p 2 1 /2λ D , where p 1 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 1 (E) given by Eq. (11), we found that the n-th position-dependent optical response coefficient for this setup is given by: 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 1 (n; x, y), we determined the monochromatic fringe shift φ Iter 1,Poly (x) at the design energy by using the iterative phase retrieval algorithm summarized in Eq. (13).  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 Theo D , 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 2 Poly (x)), which is derived from the zero-th diffracted order of the fringe pattern. The dashed red line, p Iter D (x), is the retrieved ray-sum at design energy E D using Algorithm 1. After 11 iterations, the relative error falls in 0.03%. 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 Poly (x) = − log(A 2 Poly (x)) measured from the zero-th order of the simulated fringe pattern. The blue trace plots the theoretical values of the monochromatic ray-sum p Theo D (x) at E D = 17 keV. As is shown in Fig. (4), the dashed red trace is the iteratively retrieved monochromatic ray-sum p Iter D (x) at the design energy E D = 17 keV. The retrieved monochromatic ray-sum p Iter D (x) agrees with the theoretical values p Theo D (x) within 0.03% through 11 iteration steps. Figure 5 shows the results for the fringe phase shifts in this simulation. In 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 φ Theo 1,D (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 (φ Theo 1,D (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 φ Iter 1,D (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.

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 5 E −3.5 + bZ 2 E −2 + 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 −n , 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 −n , thereby one can apply the Algorithm 1 to retrieve the monochromatic attenuation function A 2 (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][34][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 −n . Figure 6 shows the retrieving results. The solid blue line φ Theo 1,D (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 φ Theo 1,D (x). While the blue trace (φ Theo 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 −n , 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 (φ Theo 1,D (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 (φ Theo 1,D (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. 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.

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 −n 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 Iter D (x) to the position-dependent effective spectrums Θ(E; x). We then retrieved the monochromatic fringe phase shift φ Iter 1,D (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 φ Iter 1,D (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.

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 2 (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 2 (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][26][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 −n . 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.