Expand this Topic clickable element to expand a topic
Skip to content
Optica Publishing Group

Detailed analysis of the interference patterns measured in lab-based X-ray dual-phase grating interferometry through wave propagation simulation

Open Access Open Access

Abstract

In this work, we analyze the interference patterns measured in lab-based dual-phase grating interferometry and for the first time explain the spatial dependencies of the measured interference patterns and the large visibility deviations between the theoretical prediction and the experimental results. To achieve this, a simulator based on wave propagation is developed. This work proves that the experimental results can be simulated with high accuracy by including the effective grating thickness profile induced by the cone-beam geometry, the measured detector response function and a non-ideal grating shape. With the comprehensive understanding of dual-phase grating interferometry, this provides the foundations for a more efficient and accurate algorithm to retrieve sample’s structure information, and the realistic simulator is a useful tool for optimizing the set-up.

© 2023 Optica Publishing Group under the terms of the Optica Open Access Publishing Agreement

1. Introduction

X-ray Talbot grating interferometry (XT-GI) is a powerful multi-modality imaging technique [15], and widely implemented in medical and material imaging. In order to realize Talbot grating interferometry, a phase grating $G_{1}$ is used to introduce a spatially periodic phase modulation on the incident wavefront. According to the Talbot effect [6], periodic interference patterns will alternatively appear and disappear along the propagation axis. Then, a detector is fixed at one of the fractional Talbot distances to record the fringe. However, since the fringe is normally too small to be directly resolved by the detector, an analyser absorption grating $G_{2}$ is required immediately in front of the detector to help resolve the pattern by forming a Moiré fringe [1] or transforming local fringe position into intensity oscillation with the assist of the phase stepping approach [3]. Moreover, if the X-ray source cannot provide enough transverse spatial coherence, a source grating $G_{0}$ is utilized to decouple spatial resolution from spatial coherence [4]. In other cases, $G_{0}$ is not necessary, for example, if X-rays are generated by a micro-focus tube or synchrotron radiation is used [7,8]. Three types of complementary information of the sample can be retrieved simultaneously from Talbot grating interferometry: absorption contrast, differential phase contrast (DPC) and dark-field signal [5]. Absorption and differential phase contrasts result from the line integrals along the ray paths of the imaginary and real parts of the sample’s refractive indexes, respectively. The dark-field signal originates from the decoherence of the wavefront due to small angle scattering when X-rays propagating through the sample, which provides access to sub-micro structure information beyond the direct resolution of the interferometry [9].

However, the analyser absorption grating $G_{2}$ in Talbot grating interferometry blocks almost half of information-carrying photons, which greatly reduces the dose efficiency, and it is also challenging to fabricate high aspect ratio and large-field absorption gratings. Therefore, resolving the fringe without using an absorption grating provides obvious advantages. One solution is to make use of the magnification factor of the cone beam, so-called inverse geometry grating interferometry [8,10], whose fringe is magnified enough to be directly resolved. However, the achievable resolution of this approach is greatly influenced by the source size due to penumbra. Another promising solution is to place two phase gratings close together (see Fig. 1). Caused by the small difference in geometrical magnifications of two gratings, the generated Moiré fringe whose period is hundreds of times larger than the grating pitch can be easily resolved by the detector. This set-up is commonly called dual-phase grating interferometry (DPGI) [11,12]. Compared to Talbot grating interferometry, another benefit offered by DPGI is that tuning the correlation length doesn’t influence the sample magnification, which provides a solution to find the relation between the sample’s dark-field signal and the correlation length without image registration, named tunable dark-field imaging [12]. This relation is the key to quantitatively access sub-micron structure information of the sample [9]. Following these pioneer works, an analytical model for DPGI has been formulated in Wigner distribution and a fringe-formation mechanism is provided [13,14]. Due to its importance to this paper, a brief overview of their theoretical work will be given. As the schematic of the set-up is shown in Fig. 1, each phase grating is considered as a beam splitter that splits the incoming wavefront into a series of tilted planar waves according to Fourier coefficients of the grating transmission function, as described by the angular spectrum theory [15]. After being modulated by two gratings successively, the final wavefront on the detector plane is formed by a sum of all tilted planar waves, some of which are diffracted by two gratings whose diffraction orders are the same but opposite in sign, generating signals with much lower spatial frequencies than others. These low-frequency fringes are called beat patterns in DPGI, and empirically, only the beat patterns for the first three diffraction orders are large enough to be resolved due to the finite pixel size and fast decrease in magnitude for the signals of higher frequencies. Other unresolvable oscillations only contribute to the background. With $l$ indexing the diffraction orders, the final normalized oscillation intensity along the transverse direction ($x$) can be formulated as [14]:

$$\frac{I_{\textrm{res}}(x)}{I_{0}} = 1 + \sum_{l=1}^{3} V_l \cdot \cos\left(2\pi l \frac{x}{p_{\textrm{fr}}}\right),$$
with $I_{0}$ normalization constant. The visibility coefficient $V_l$ is related to Fourier coefficients of the grating transmission function, set-up geometry, source size, effective spectrum, and detector response function [14,16]. If the two gratings have the same period $p_{\textrm {g}}$, the first-order fringe period ($l = 1$) is given by:
$$p_{\textrm{fr}} = \frac{R_{\textrm{s}} + R_{\textrm{g}} + R_{\textrm{d}}}{R_{\textrm{g}}} \cdot p_{\textrm{g}},$$
with $R_{\textrm {s}}$ the distance between source and $G_1$, $R_{\textrm {g}}$ the inter-grating distance and $R_{\textrm {d}}$ the distance between $G_2$ and detector (see Fig. 1).

 figure: Fig. 1.

Fig. 1. Schematic of dual-phase grating interferometry. The X-rays are generated by a micro-focus tube, propagate through two adjacent phase gratings, and are recorded by an integrating detector.

Download Full Size | PDF

Although several simulation frameworks have been developed to model fringes formed in DPGI [13,14,17] and the theoretical predictions for the global oscillation intensity given by Eq. (1) are in good agreement with their idealized simulation results [14], these models fail to explain some properties of typical measurements. Firstly, these models expect a homogeneous fringe [14,17], while the experimentally recorded fringe as in Fig. 2(a) presents a clear spatial dependency. This interference pattern is measured with 2.24 mm inter-grating distance, and other measurement details can be found in subsection 3.1. To better display the intensity oscillation on the x-axis, a one-dimensional fringe profile (see Fig. 2(b)) is extracted from Fig. 2(a) by averaging 50 pixels along the y-axis for a better statistics, and the averaged region is denoted by a gold rectangle. It can be seen from Fig. 2(b) that the interference pattern is spatially symmetrical. In the central part of the pattern, there are intermediate peaks between every two main peaks, indicating that the interference pattern in DPGI does contain signals of higher frequencies according to Eq. (1). However, these higher frequency signals gradually disappear from the center to the sides. Secondly, the analytical model cannot accurately predict the fringe contrast of the measured result [14]. To quantify the fringe contrast, the visibility $V$ is used, defined by the maximum ($I_{\textrm {max}}$) and minimum ($I_{\textrm {min}}$) intensity as:

$$V = \frac{I_{\textrm{max}} - I_{\textrm{min}}}{I_{\textrm{max}} + I_{\textrm{min}}}.$$

The visibility for each pixel in Fig. 2(b) was extracted by searching the maximum and minimum intensities within one period around this pixel according to Eq. (3) and plotted in Fig. 2(c). The period is determined by locating the first harmonics of the Fourier transform of intensity oscillation in Fig. 2(b). From Fig. 2(c) we can see that the visibility increases from the center and reaches the maximum on both sides of the center. Then, the visibility constantly decreases when further away from the center. This spatial dependency of the measured visibility is also mentioned by Bopp et al. [17]. However, both analytical model [14] and simulation frameworks [14,17] result in a uniform visibility distribution. Moreover, the maximum visibility shown in Fig. 2(c) is below $0.15$, while the analytical model [14], as implemented based on our experimental settings, predicts a constant visibility of $0.287$ along the x-axis.

 figure: Fig. 2.

Fig. 2. Overview of the measured interference pattern. (a) shows a recorded image when the inter-grating distance is at 2.24 mm, and the gold square region is averaged along the y-axis to yield a line plot of the interference pattern in (b). The visibility $V$ for each pixel in (b) is extracted and plotted in (c).

Download Full Size | PDF

It is crucial to correctly understand what causes the deviations between the analytical prediction and the experimental results, which is a prerequisite for accurately extracting the sample’s three-modality information or optimizing the set-up. Therefore, this paper aims at filling the gap between the theoretical and measured results by implementing a wave-propagation simulation.

2. Simulation methods

In this section, we will first introduce the assumptions we made to interpret spatial-depending interference patterns, followed by elaborating the methods used in our DPGI simulator to justify this hypothesis, including modelling the effective thickness profile of the grating illuminated by a spherical wave, implementing the spherical wave propagation in DPGI, handling polychromatic source and extended source size, and modelling detector response.

The strong spatial dependency of the fringe (see Fig. 2) is hypothesized to result from the grating thickness profile varying with the X-ray incident angle in a cone-beam geometry, which is supported by two pieces of evidence. Firstly, the phase gratings used have a very high aspect ratio (around $56$), making the effective thickness profile very sensitive to the incident angle. As shown in Fig. 3, a rather small incident angle $\alpha$ will lead to a significant change in the effective thickness profile $t_1(x)$, since the non-normal incidence x-rays often cross the material-air boundary when propagating through a grating with larger aspect ratio. Secondly, according to the results of Shashev et al. [18], the visibility of the interference pattern generated by a single phase grating at the 1st fractional Talbot distance varies with the X-ray’s incident angle, and the maximum visibility is not at normal incidence but with a certain inclination. Similarly, the visibility in DPGI symmetrically peaks at both sides of the center as if a certain incident angle is met (see Fig. 2(c)). To prove this hypothesis, a wave propagation based simulation is implemented, and the physical principles have been developed as follows.

 figure: Fig. 3.

Fig. 3. Schematic for a ray trace simulation to generate an effective thickness profile of the grating with rectangular bars in a spherical wave. The effective thickness profile $t_1(x)$ is calculated by accumulating the thickness density map $g_1(x,z)$ along the ray paths $S(z;x)$. Due to the high aspect ratio of the grating, the effective thickness profile $t_1(x)$ is very sensitive to even a small change in the incident angle $\alpha$. The range of the incident angle for our set-up is from $-1.9^{\circ }$ to $1.9^{\circ }$.

Download Full Size | PDF

2.1 Effective thickness profile of grating for a spherical wave

As shown in Fig. 3, under the paraxial approximation, a discretized monochromatic spherical wavefront is generated at the entrance of $G_1$ and can be described in one dimension [15]:

$$u_1(x) = \frac{\exp[ik(x^2/(2R_{\textrm{s}}))]}{R_{\textrm{s}}},$$
with $k = \frac {\lambda }{2\pi }$ the wave number and $x$ the transverse position. The incident angle at each $x$ can be calculated according to the geometry of the system: $\alpha _1(x) = \frac {x}{R_{\textrm {s}}}$, or by the phase of the wavefront $\varphi _1(x)$ [15]:
$$\alpha_1(x) = \frac{\lambda}{2\pi}\frac{\partial}{\partial x} \varphi_1(x).$$

For a grating made of silicon, a thickness density map of $G_1$ is created as $g_1(x,z)$ where the air-region: $g_1(x,z) = 0$ and the silicon-region: $g_1(x,z) = \frac {d}{N}$ with the designed grating thickness $d$ and the number of discretized slices in the z-axis $N$. With the ray path $S(z;x)$ determined by the incident angle for each $x$ at the wavefront $u_1(x)$, the effective thickness profile $t_1(x)$ is calculated by summing up the thickness density map along the paths:

$$t_1(x) = \sum_{i=0}^{N-1} g_1(S(i;x),i).$$

Since the thin grating can be considered as a sufficiently weak scattering object introducing negligible disturbance on the ray path, which satisfies the precondition of projection approximation [19], the wavefront at the exit surface $u'_{1}(x)$ is approximated by:

$$u'_{1}(x) = u_1(x) \cdot e^{{-}ik\delta_{\textrm{Si}}t_1(x)} \cdot e^{{-}k\beta_{\textrm{Si}}t_1(x)},$$
with the refractive index defined as: $n = 1 - \delta + i\beta$.

2.2 Spherical wave propagation in DPGI

To find the final interference pattern at the detector plane, the wavefront behind the first grating $u'_{1}(x)$ is propagated through $G_2$ and then reaches the detector. However, it is challenging to numerically propagate a spherical wave, since properly discretizing a spherical wavefront normally requires impractically dense sampling due to its quadratic phase terms [20], as a large field of view is a design specification for this work. This problem has already been addressed by previous works [1922], where the basic concept is to separate rapid oscillating terms to relieve the sampling requirement. In this work, we adopt the method from Munro [21], which has been implemented in multi-slice propagation in inhomogeneous space. Firstly, a divergent beam is converted into a planar wave by applying a coordinate transformation to subtract its quadratic phase term. As a result, the wavefront $u'_{1}(x)$ is transformed from the global coordinate $(x)$ into the primed coordinate system $(x')$ through:

$$v_1(x') = \frac{u'_1(x) \cdot R_{\textrm{s}}}{\exp[ikx^2/(2R_{\textrm{s}})]}.$$

Then, the wavefront $v_1(x')$ is propagated from the exit of $G_1$ to the entrance of $G_2$ in free space and we choose to use the Fresnel transfer function approach [15]. To simulate accurate results with negligible aliasing we still need to properly sample the aperture function $v_1(x')$ and the propagator. According to the criteria given in [15], the minimum sampling rate for our simulation is found at 250 µm−1 and the sampling ratio $Q$ is set at $1.1$. The sampling ratio $Q$ determines the size of zero padding on both sides of the aperture function $v_1(x')$ as $\frac {Q-1}{2}$. Then, the wavefront at the entrance surface of $G_2$ is formulated as:

$$v_2(x') = \mathfrak{F}^{{-}1}\{\mathfrak{F}\{ZP(v_1(x'))\} \cdot \exp[{-}j \pi \lambda f_{x}^{2} \cdot (R_{\textrm{g}} / M_1)]\},$$
where $\mathfrak {F}$ is Fourier transform, $ZP$ represents zero padding, $f_x$ is the spatial frequency in the x-axis, and the propagation distance from $G_1$ to $G_2$ in the primed coordinate system is reduced to $R_{\textrm {g}} / M_1$ by a magnification factor $M_1$ given by:
$$M_1 = \frac{R_{\textrm{s}}+R_{\textrm{g}}}{R_{\textrm{s}}},$$

The next step is to transform $v_2(x')$ back into the global coordinate system to retrieve the complete wavefront before entering $G_2$:

$$u_2(x) = \frac{1}{R_{\textrm{s}}+R_{\textrm{g}}} \cdot v_2\left(\frac{x}{M_1}\right) \cdot \exp\{ikx^2/[2(R_{\textrm{s}}+R_{\textrm{g}})]\}.$$

Compared to $u'_{1}(x)$, the x coordinate in $u_2(x)$ is magnified by a factor $M_1$ after the propagation, and the amplitude is dropped by $\frac {1}{M_1}$. Then, the projection approximation and Fresnel scaling algorithm are repeated to calculate the wavefront at the detector plane. It is important to mention that the X-ray’s incident angle on $G_2$ is approximated by geometry instead of the phase of the wavefront as $\alpha _2(x) = \frac {x}{R_{\textrm {s}}+R_{\textrm {g}}}$ since the disturbances on ray paths caused by $G_1$ are small enough (around tens of micro-degree) to be ignored. Moreover, the sampling rate for the transmission function of $G_2$ is lower than $G_1$ due to the transversely magnified coordinate.

2.3 Polychromatic and extended source

The algorithm described above is only able to simulate the wavefront at the detector plane for a monochromatic point source. However, X-rays generated by a tube source normally contain a continuous spectrum. With the proof that the fringe generated by a single source with a polychromatic spectrum is the sum of intensities for different frequencies provided by Cowley [23], the interference pattern from the polychromatic source is then simulated by weighted summing up intensity pattern for each photon energy, and the weights are given by the normalized effective intensity spectrum.

Furthermore, an extended source can be represented by a group of incoherent point sources which reproduce the intensity distribution [24], and the final interference pattern is a sum of individual intensities generated by each point source [23]. However, simulating interference patterns for every point source is computationally heavily. According to the theoretical results obtained by Yan et al. [25], the effect of an extended source in DPGI can be modelled by applying a convolution kernel, which greatly increases computing efficiency. When rewriting the formula which describes the relationship between an off-center line source and the total fringe phase shift in DPGI [25], we can derive the lateral displacement relationship between a point source and the corresponding interference pattern if two gratings are identical:

$$I(x+\Delta_{\textrm{s}}) = 1 + \sum_{l=1}^{\infty} V_l \cdot \cos\left[2\pi l \frac{1}{p_{\textrm{fr}}}(x+\Delta_{\textrm{s}})\right],$$
with $\Delta _{\textrm {s}}$ as the lateral shift of a point source. Combining both the superposition and space-invariant properties, DPGI can be proven to be a spatial invariant imaging system [15], and the lateral shift of a point source leads to a displacement of the fringe with the same distance but in opposite direction according to Eq. (12). Therefore, we can model the effect of the extended source by a convolution between the interference pattern simulated by a point source and the intensity distribution of the X-ray source.

2.4 Detector response

Only the intensity of the wavefront on the detector is measurable, and correctly modelling the detector response function is important in DPGI since the fringe is directly resolved without analyser grating. The simplest model of the response function is to re-bin the simulated fringe to the pixel size, with an assumption of shape edge. However, due to the charge sharing effect, readout sampling or electronics, the edge response of a real detector is not so sharp [26], instead, the intensity will be shared by a few neighbouring pixels resulting in Gaussian or exponential line spread function (LSF) which is the first difference of an edge response [27]. Once the LSF is determined, a more realistic fringe can be simulated through a convolution between the LSF and the interference pattern simulated with a ideal sharp edge detector.

3. Experiments and results

3.1 Experimental settings

The experiments were carried out on a transmission-type microfocus X-ray tube (FXT-160.51, FEINFOCUS GmbH, Germany) with 3 µm tungsten coated CVD diamond transmission anode, operated at 40 kV and 200 mA under the microfocus mode. The phase gratings (GRATXRAY, Switzerland) made of silicon were identical with 1 µm period and $0.5$ duty cycle, and the thickness of grating is 28.17 µm introducing a $\pi$ phase shift at the design energy at 22 keV. The fringes were recorded by an sCMOS camera (Gsense, PHOTONIC SCIENCE, UK) coupled to a 100 µm CsI:TI scintillator. The detector had an active input area of $67 \times 67 \;\mathrm {\mu }^{2}$ and its effective pixel size was $16.4 \times 16.4 \;\mathrm {\mu }^{2}$. The first phase grating $G_1$ was placed 500 mm away from the X-ray tube and the detector was located 500 mm from $G_1$, resulting in a source-to-detector distance of 1000 mm. The second grating $G_2$ was placed downstream $G_1$ and the inter-grating distance was tuned by a linear stage (CLS-5282, SmarAct GmbH, Germany) from 2.24 mm to 9.24 mm in steps of 1 mm. For each inter-grating distance, the fringe was recorded with an exposure time of 25 s.

3.2 Simulation settings

All parameters set in the simulation were set identical to in the experiment. An initial spectrum emitted from the tube was modelled by a Monte Carlo simulation [28]. Then this spectrum was filtered by several uniform attenuation materials, including air (1000 mm), a carbon fibre input window of the detector (1 mm), silicon substrates of two gratings ($2 \times 0.222\;\textrm {mm}$), and a diamond window at the exit surface of X-ray tube (0.26 mm). Finally, the average deposited energy in the scintillator for each photon energy is incorporated to simulate the effective spectrum. Then the effective spectrum was normalized and presented in Fig. 4. The spectrum started at 8 keV since lower energy photons were significantly absorbed, mainly by the silicon substrates, and their contributions were very limited. The sampling rate for the spectrum is 4 keV−1. Based on the measurement of the JIMA pattern, the FWHM of the source size was around 5 µm. Therefore, a Gaussian kernel with 5 µm FWHM was used to represent the effect of the extended source in the simulation. In the simulation, a 25.42 mm aperture size was applied at the $G_1$ plane, facing the center of the source. Due to the two-fold magnification, the effective aperture size at the detector plane was 50.84 mm.

 figure: Fig. 4.

Fig. 4. Normalized effective intensity spectrum employed in the simulation for tungsten target at 40 kV, including uniform attenuation materials and quantum efficiency of the detector.

Download Full Size | PDF

3.3 Comparison between experiment and simulation

The interference pattern recorded at the 2.24 mm inter-grating distance is shown in Fig. 2(a). To better reveal the properties of the fringe, the intensity oscillation in real space is transformed into the visibility coefficients in Fourier space according to Eq. (1). For each pixel in Fig. 2(a), one period of intensity oscillation around the pixel in the transverse direction is segmented and a Windowed Fourier transform (WFT) [29] is applied to retrieve the first three orders’ visibility coefficients. This means one recorded image can be converted into three visibility coefficient images whose values are uniform along the vertical direction due to the accurate alignment except for statistical errors. Therefore, the measured spatial dependencies of visibility coefficients for the first three orders are extracted and shown as blue solid lines in Fig. 5(a), (b), and (c), respectively, calculated by averaging the visibility coefficient images along the vertical direction for better statistics. The corresponding error bars (one standard deviation) are presented as shaded regions. To validate the assumption that the spatial dependency of the interference pattern results from the divergent beam, a simulation is performed with angular-dependent effective thickness profiles (see Fig. 3), a perfectly sharp-edge detector, and ideal rectangular grating bars. Then, the simulated fringe is transformed into the first three orders’ visibility coefficients by applying the WFT and the results are plotted in Fig. 5 as black dotted lines (angular-dependent effective thickness profile: ADTP). In addition, the x-axis is transformed into the incident angle according to the geometry of the set-up. Comparing the results between the experiment (Exp) and simulation (ADTP), except for the third order, the simulator can correctly model the trends of how the visibility coefficients change with the incident angles. The visibility of the first order is low when the incident angle is around zero and gradually increases and peaks at $+/-0.5^\circ$, while the maximum magnitude of the second order appears when normal incidence and its visibility continue to drop as the incident angle increases. Furthermore, the corresponding theoretical baselines predicted by the analytical model [14] are plotted as red solid lines (Theory). Since the analytical model ignores the angular-dependent effective thickness profiles, it is constant at all incident angles. Moreover, as the rectangular shape of the grating is presupposed in the analytical model, the theoretically predicted visibilities match very well with the simulation results only if the incident angle is zero, giving a rectangular effective thickness profile. However, visibility coefficients for all diffraction orders are generally overestimated in the simulation, and one significant difference is that the second order suppresses the first and third orders more effectively in the simulation than measurement with a much higher magnitude over a wider range. Although the simulation results clearly indicate that the spatial dependency of the interference pattern comes from the changing effective thickness profile of the grating in a divergent beam, there is still a significant difference between the experimental and the simulated interference patterns.

 figure: Fig. 5.

Fig. 5. Comparison of the visibility coefficients for the first three orders between experimental, simulated and theoretical prediction results. The experimental results measured at 2.24 mm inter-grating distance are plotted as blue solid lines (Exp) in figures (a), (b), and (c) for the first three orders, respectively. The corresponding error bars for one standard deviation are represented by shaded regions. Moreover, the theoretical baselines for each order are plotted as red solid lines. The simulation with angular-dependent effective thickness profiles (ADTP) of the gratings, an ideal sharp-edge detector and perfect rectangular grating bars are plotted as black dotted lines. Once the real detector response function (RDR) is included (see Fig. 6), the simulator can produce more accurate results plotted as black dash-dot lines (ADTP and RDR). When the non-ideal grating shape (NIGS) is included in the simulation (see Fig. 7), not only are the spatial dependencies of the measured visibility coefficients correctly predicted by the simulator but the magnitudes of these coefficients are well reproduced with tolerable deviations, plotted as black dashed lines (ADTP, RDR and NIGS). These errors most probably come from the non-periodic and irregular shape of the gratings which cannot be accurately modelled.

Download Full Size | PDF

To further improve the simulator, the edge response of the detector was measured and its line spread function (LSF) was included. A sharp edge of a 0.09 mm Tantalum foil attaching to the detector was imaged a half meter away from the source, which greatly restrained the blurring effect from the extended source size. Then, the edge response of the detector was recorded and its first difference (LSF) was extracted and fitted by two symmetrical exponential functions (Laplace distribution), as shown in Fig. 6. This fitted LSF was convolved with the previously simulated results to include the effect of the detector unsharpness (real detector response, RDR), and the improved simulation results are presented in Fig. 5 as black dash-dot lines (ADTP and RDR). Compared to the simulation results that only consider the angular dependency of the effective thickness, when the effect of the detector unsharpness is also included the visibility coefficients for all diffraction orders are generally reduced and the higher orders’ signals suffer from more severe damping effects since the LSF of the detector acts as a low-pass filter. Though the improved simulator produces more accurate results with more similar visibility coefficients to the measured ones, there is still a margin for improvement of the simulation accuracy.

 figure: Fig. 6.

Fig. 6. Measured edge response for the detector (a) and the measured line spread function (LSF) is the first difference of the edge response represented by the red dots in (b). The fitted LSF from two symmetrical exponential functions is plotted as the blue line in (b).

Download Full Size | PDF

According to the cross-sectional image from scanning electron microscopy (SEM) of the grating used (see Fig. 7(a)), presupposing the shape of the grating as rectangular is not a proper approximation. Instead, a trapezoidal shape is more similar to the real grating. Moreover, in Fig. 7(b) a zoomed-in image indicates that the silicon and air regions are not equivalent, with less area occupied by silicon. To include the influence of a more realistic shape of the grating, a non-ideal grating shape is simulated and one period of it is drawn in Fig. 7(c), with a trapezoidal shape and unequal silicon and air areas. Then, the corresponding effective thickness profiles of two gratings are determined in the ray trace simulation according to Eq. (6). Once the non-ideal grating shape (NIGS) is also included, plotted as black dashed lines (ADTP, RDR and NIGS) in Fig. 5, the simulator can model the experimental results very well with limited deviations. It is clear from the above results that, although the experimental conditions are approximated as much as possible, they do not perfectly match. As an example, in Fig. 7(b) it can be observed that the real grating bars are not perfectly periodic, while the exact shape of Fig. 7(a) is very hard to model.

 figure: Fig. 7.

Fig. 7. (a) Cross-sectional image from scanning electron microscopy (SEM) of the grating used in the experiments. A zoomed-in region from (a) indicated by a black rectangle is shown in (b). Based on the measured morphology of the grating, a non-ideal grating is simulated and the shape for one period is presented in (c).

Download Full Size | PDF

In DPGI, one of the advantages is that the interference pattern can be conveniently tuned by changing the inter-grating distance, allowing for tunable dark-field imaging. Therefore, the experimental and simulated results are compared at three different inter-grating distances. In Fig. 8(a), (c), and (e), the comparisons for the first two orders’ visibility coefficients in Fourier space between the measurement (solid lines) and the simulation (dash lines) at 2.24 mm, 4.24 mm and 6.24 mm inter-grating distances, respectively. The third order is ignored here since their contributions are limited as indicated in Fig. 5. In general, the simulation can accurately reproduce the experimental results at different inter-grating distances across the full field of view. It can be seen from Fig. 8(a), (c), and (e) that among all three different inter-grating distances, the visibility coefficients for the first order present similar tendencies where a valley always appears in the A1 region and two peaks are located at the two sides of the A1, while the visibility coefficients for the second order show different tendencies with the inter-grating distance increased. Three regions in Fourier space for different incident angles are selected and the corresponding interference patterns in real space are compared between the measurement (solid lines) and simulation (dash lines), shown in Fig. 8(b), (d), and (f). In the A1 region, the magnitudes of the two visibility coefficients are small and comparable, therefore, not only do the corresponding real space fringes present low contrast but the observable intermediate peaks appear. In the A2 region, the 1st-order visibility coefficients reach the maximum and show much larger magnitudes than the second order, leading to high-contrast and single-frequency patterns in real space. As the incident angle increases further, in the A3 region, both visibility coefficients drop but the first orders are dominant, as a result, the interference patterns present lower contrast and have no intermediate peaks. It is not surprising that once the simulator can model the experimental results well in Fourier space, the interference patterns in real space can be correctly simulated.

 figure: Fig. 8.

Fig. 8. (a), (c) and (e) compare the experimental (solid line) and simulated (dash line) results in Fourier space represented by the visibility coefficient for the first two orders at three inter-grating distances ($2.24$, $4.24$ and 6.24 mm). Three regions according to the incident angle are selected and the corresponding real space fringes are compared in (b), (d) and (f).

Download Full Size | PDF

Finally, the macroscopic properties of the interference patterns are compared at eight different inter-grating distances from $2.24$ to 9.24 mm. Figure 9(a) compares the periods between the theoretical prediction, simulation and measurement. The theoretical periods are calculated according to Eq. (2), and the periods for simulation and experiment are extracted from the first harmonics of the Fourier transform of one-dimensional interference patterns. We can see from Fig. 9(a) that the extracted periods are almost identical with the limited standard deviation averaged among all inter-grating distances of $0.07$ pixels. Figure 9(b) compares the averaged visibility for the first two orders over the full field of view between the experiment and simulation at different inter-grating distances. Although the peaking position and maximum visibility for both orders are modelled accurately, with the inter-grating distance increased, the averaged visibility in measurement drops slower than the simulated results after the peaking position. This mismatch might result from two factors. Firstly, the irregular and non-uniform shape of the grating cannot be exactly reproduced in the simulation (see Fig. 7), and we believe this is the main reason. Secondly, the identical detector response as in the real experiment is also very hard to model, especially considering the pixelated response function.

 figure: Fig. 9.

Fig. 9. The macroscopic properties of the interference patterns are compared. (a) Comparison for the periods among theoretical prediction, measurement and simulation at eight different inter-grating distances from $2.24$ to 9.24 mm. (b) compares the averaged visibility for the first two orders over the full field of view between the experiment and simulation at different inter-grating distances.

Download Full Size | PDF

4. Discussion and conclusions

As mentioned in the previous sections, a recently developed technique, dual-phase grating interferometry (DPGI) has the potential to improve the performance of Talbot grating interferometry. Directly resolving the fringe without requiring an analyser absorption grating provides an advantage in dose efficiency and eliminates the difficulties in grating fabrication [11,12]. Furthermore, changing the correlation length in DPGI does not influence the sample’s magnification factor, which provides a more convenient way to implement tunable dark-field imaging to investigate the sample’s sub-micro structure information [12].

Although a rigorous mathematical model has been derived to understand DPGI [13,14], this model cannot explain the spatially depending interference pattern and the huge visibility difference between the theory and our experimental results. This paper provides a DPGI simulator based on wave propagation to successfully diminish the gap between the theoretical prediction and experimental results. Firstly, since the simulation can reproduce the spatial dependency with very high accuracy, the spatial-depending feature is explained by the local effective thickness profile changed with the incident angle. Due to the large aspect ratio of the used grating, the local effective thickness distribution is greatly influenced by the cone-beam geometry. Meanwhile, in the near-field Fresnel propagation, the major contributions to a certain point in the propagated wavefront come from a limited region around this point at the aperture plane [15,30]. Consequently, the locally recorded interference pattern is directly determined by the local effective thickness profile according to Eq. (1). Secondly, the simulation with a sharp edge detector and perfect rectangular gratings which are also assumed in the theoretical model [14] cannot precisely model the measurement results. Once the simulator includes the measured camera response function and non-ideal grating shape, it can reproduce the experimental results with high accuracy. These results indicate that the camera response function and precise grating shape play a very important role in DPGI, which should always be considered when we implement this technique.

This simulator also reveals several potential developments in the future. Firstly, the small angle scattering information or dark-field signal of the sample can be represented by the visibility reduction in Talbot grating interferometry with the assistance of the phase stepping approach [5]. By stepping one of the gratings one or more periods, for each pixel the visibility (normalised amplitude) is extracted from the phase stepping curve after approximating the intensity oscillation as a sinusoidal function [3,5]. A similar phase stepping approach is also used in DPGI for decoupling the image resolution from the period of the fringe [12]. However, the interference patterns in DPGI consist of multiple harmonics (different diffraction orders) which present distinct spatial dependencies (see Fig. 8). Therefore, when extracting the dark-field signal in DPGI, any phase stepping approach that ignores higher order harmonics would be inaccurate. Sufficient sampling is therefore needed to retrieve the magnitude of the visibility coefficients. With a comprehensive understanding and accurate simulation of the spatial dependencies of the visibility coefficients, a better algorithm can be developed to efficiently and accurately extract the dark-field signal from the sample in DPGI. Furthermore, this simulator can be used to optimize the set-up performance further. For example, Fig. 8 shows that the highest visibility coefficients are achieved in the A2 region for different inter-grating distances for current set-up configurations and the higher visibility is always appreciated in the grating interferometry for better signal-to-noise ratio. Therefore, this simulator can be used to locate the high visibility regions for different set-up configurations and find the best configuration leading to an optimum visibility. Moreover, the maximum visibility coefficients result from a non-rectangular effective thickness profile, so this implies a bended grating with this optimum shape would be able to generate uniform high contrast fringes in DPGI.

In conclusion, this paper presents a realistic DPGI simulator which successfully models experimental results and interprets the differences between the theoretical prediction and measurement. This simulator paves the way to develop a better algorithm to retrieve the sample’s structure information and optimize the set-up.

Funding

Interreg Vlaanderen-Nederland (Smart*Light); European Regional Development Fund (Smart*Light); Provincie Oost-Vlaanderen (Smart*Light); Fonds Wetenschappelijk Onderzoek (3179I12018, 3G036518, 3G010820); the Swiss LOS Lottery Fund of the Kanton of Aargau, CH.

Acknowledgments

We acknowledge our technicians Sander Vanheule, Yen Decappelle and Iván Josipovic from the Radiation Physics research group at Ghent University for their technical support.

Disclosures

The authors declare no conflicts of interest related to this article.

Data availability

Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.

References

1. A. Momose, S. Kawamoto, I. Koyama, Y. Hamaishi, K. Takai, and Y. Suzuki, “Demonstration of x-ray talbot interferometry,” Jpn. J. Appl. Phys. 42(Part 2, No. 7B), L866–L868 (2003). [CrossRef]  

2. A. Momose, “Recent advances in x-ray phase imaging,” Jpn. J. Appl. Phys. 44(9A), 6355–6367 (2005). [CrossRef]  

3. T. Weitkamp, A. Diaz, C. David, F. Pfeiffer, M. Stampanoni, P. Cloetens, and E. Ziegler, “X-ray phase imaging with a grating interferometer,” Opt. Express 13(16), 6296–6304 (2005). [CrossRef]  

4. F. Pfeiffer, T. Weitkamp, O. Bunk, and C. David, “Phase retrieval and differential phase-contrast imaging with low-brilliance x-ray sources,” Nat. Phys. 2(4), 258–261 (2006). [CrossRef]  

5. F. Pfeiffer, M. Bech, O. Bunk, P. Kraft, E. F. Eikenberry, C. Brönnimann, C. Grünzweig, and C. David, “Hard-x-ray dark-field imaging using a grating interferometer,” Nat. Mater. 7(2), 134–137 (2008). [CrossRef]  

6. H. F. Talbot, “Lxxvi. facts relating to optical science. no. iv,” The London, Edinburgh, Dublin Philos. Mag. J. Sci. 9(56), 401–407 (1836). [CrossRef]  

7. M. Engelhardt, J. Baumann, M. Schuster, C. Kottler, F. Pfeiffer, O. Bunk, and C. David, “High-resolution differential phase contrast imaging using a magnifying projection geometry with a microfocus x-ray source,” Appl. Phys. Lett. 90(22), 224101 (2007). [CrossRef]  

8. A. Momose, H. Kuwabara, and W. Yashiro, “X-ray phase imaging using lau effect,” Appl. Phys. Express 4(6), 066603 (2011). [CrossRef]  

9. M. Strobl, “General solution for quantitative dark-field contrast imaging with grating interferometers,” Sci. Rep. 4(1), 7243 (2015). [CrossRef]  

10. T. Donath, M. Chabior, F. Pfeiffer, O. Bunk, E. Reznikova, J. Mohr, E. Hempel, S. Popescu, M. Hoheisel, M. Schuster, J. Baumann, and C. David, “Inverse geometry for grating-based x-ray phase-contrast imaging,” J. Appl. Phys. 106(5), 054703 (2009). [CrossRef]  

11. H. Miao, A. Panna, A. A. Gomella, E. E. Bennett, S. Znati, L. Chen, and H. Wen, “A universal moiré effect and application in x-ray phase-contrast imaging,” Nat. Phys. 12(9), 830–834 (2016). [CrossRef]  

12. M. Kagias, Z. Wang, K. Jefimovs, and M. Stampanoni, “Dual phase grating interferometer for tunable dark-field sensitivity,” Appl. Phys. Lett. 110(1), 014105 (2017). [CrossRef]  

13. A. Yan, X. Wu, and H. Liu, “Quantitative theory of x-ray interferometers based on dual phase grating: fringe period and visibility,” Opt. Express 26(18), 23142–23155 (2018). [CrossRef]  

14. A. Yan, X. Wu, and H. Liu, “Predicting fringe visibility in dual-phase grating interferometry with polychromatic x-ray sources,” J. X-Ray Sci. Technol. 28(6), 1055–1067 (2020). [CrossRef]  

15. J. W. Goodman, Introduction to Fourier Optics. Goodman (McGraw-Hill, 1968).

16. A. Yan, X. Wu, and H. Liu, “A general theory of interference fringes in x-ray phase grating imaging,” Med. Phys. 42(6Part1), 3036–3047 (2015). [CrossRef]  

17. J. Bopp, V. Ludwig, M. Seifert, G. Pelzer, A. Maier, G. Anton, and C. Riess, “Simulation study on x-ray phase contrast imaging with dual-phase gratings,” Int. J. Comput. Assist. Radiol. Surg. 14(1), 3–10 (2019). [CrossRef]  

18. Y. Shashev, A. Kupsch, A. Lange, R. Britzke, G. Bruno, B. R. Müller, and M. P. Hentschel, “Talbot-lau interferometry with a non-binary phase grating for non-destructive testing,” in Talbot-Lau interferometry with a non-binary phase grating for non-destructive testing, (2016), pp. Tu_3_G_2–1.

19. D. Paganin, Coherent X-ray optics, 6 (Oxford University Press on Demand, 2006).

20. O. Chubar and R. Celestre, “Memory and cpu efficient computation of the fresnel free-space propagator in fourier optics simulations,” Opt. Express 27(20), 28750–28759 (2019). [CrossRef]  

21. P. R. Munro, “Rigorous multi-slice wave optical simulation of x-ray propagation in inhomogeneous space,” J. Opt. Soc. Am. A 36(7), 1197–1208 (2019). [CrossRef]  

22. M. Mansuripur, “Certain computational aspects of vector diffraction problems,” J. Opt. Soc. Am. A 6(6), 786–805 (1989). [CrossRef]  

23. J. M. Cowley, Diffraction physics (Elsevier, 1995).

24. M. Engelhardt, C. Kottler, O. Bunk, C. David, C. Schroer, J. Baumann, M. Schuster, and F. Pfeiffer, “The fractional talbot effect in differential x-ray phase-contrast imaging for extended and polychromatic x-ray sources,” J. Microsc. 232(1), 145–157 (2008). [CrossRef]  

25. A. Yan, X. Wu, and H. Liu, “Clarification on generalized lau condition for x-ray interferometers based on dual phase gratings,” Opt. Express 27(16), 22727–22736 (2019). [CrossRef]  

26. M. Boone, “New imaging modalities in high resolution x-ray tomography,” Ph.D. thesis, Ghent University (2013).

27. S. Smith, Digital signal processing: a practical guide for engineers and scientists (Elsevier, 2013).

28. J. Dhaene, “Development and application of a highly accurate polychromatic x-ray microtomography simulator,” Ph.D. thesis, Ghent University (2017).

29. A. Torre, Linear ray and wave optics in phase space: bridging ray and wave optics via the Wigner phase-space picture (Elsevier, 2005).

30. K. M. Mechlem, “Towards low-dose spectral phase-contrast x-ray imaging,” Ph.D. thesis, Technische Universität München (2020).

Data availability

Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.

Cited By

Optica participates in Crossref's Cited-By Linking service. Citing articles from Optica Publishing Group journals and other participating publishers are listed here.

Alert me when this article is cited.


Figures (9)

Fig. 1.
Fig. 1. Schematic of dual-phase grating interferometry. The X-rays are generated by a micro-focus tube, propagate through two adjacent phase gratings, and are recorded by an integrating detector.
Fig. 2.
Fig. 2. Overview of the measured interference pattern. (a) shows a recorded image when the inter-grating distance is at 2.24 mm, and the gold square region is averaged along the y-axis to yield a line plot of the interference pattern in (b). The visibility $V$ for each pixel in (b) is extracted and plotted in (c).
Fig. 3.
Fig. 3. Schematic for a ray trace simulation to generate an effective thickness profile of the grating with rectangular bars in a spherical wave. The effective thickness profile $t_1(x)$ is calculated by accumulating the thickness density map $g_1(x,z)$ along the ray paths $S(z;x)$. Due to the high aspect ratio of the grating, the effective thickness profile $t_1(x)$ is very sensitive to even a small change in the incident angle $\alpha$. The range of the incident angle for our set-up is from $-1.9^{\circ }$ to $1.9^{\circ }$.
Fig. 4.
Fig. 4. Normalized effective intensity spectrum employed in the simulation for tungsten target at 40 kV, including uniform attenuation materials and quantum efficiency of the detector.
Fig. 5.
Fig. 5. Comparison of the visibility coefficients for the first three orders between experimental, simulated and theoretical prediction results. The experimental results measured at 2.24 mm inter-grating distance are plotted as blue solid lines (Exp) in figures (a), (b), and (c) for the first three orders, respectively. The corresponding error bars for one standard deviation are represented by shaded regions. Moreover, the theoretical baselines for each order are plotted as red solid lines. The simulation with angular-dependent effective thickness profiles (ADTP) of the gratings, an ideal sharp-edge detector and perfect rectangular grating bars are plotted as black dotted lines. Once the real detector response function (RDR) is included (see Fig. 6), the simulator can produce more accurate results plotted as black dash-dot lines (ADTP and RDR). When the non-ideal grating shape (NIGS) is included in the simulation (see Fig. 7), not only are the spatial dependencies of the measured visibility coefficients correctly predicted by the simulator but the magnitudes of these coefficients are well reproduced with tolerable deviations, plotted as black dashed lines (ADTP, RDR and NIGS). These errors most probably come from the non-periodic and irregular shape of the gratings which cannot be accurately modelled.
Fig. 6.
Fig. 6. Measured edge response for the detector (a) and the measured line spread function (LSF) is the first difference of the edge response represented by the red dots in (b). The fitted LSF from two symmetrical exponential functions is plotted as the blue line in (b).
Fig. 7.
Fig. 7. (a) Cross-sectional image from scanning electron microscopy (SEM) of the grating used in the experiments. A zoomed-in region from (a) indicated by a black rectangle is shown in (b). Based on the measured morphology of the grating, a non-ideal grating is simulated and the shape for one period is presented in (c).
Fig. 8.
Fig. 8. (a), (c) and (e) compare the experimental (solid line) and simulated (dash line) results in Fourier space represented by the visibility coefficient for the first two orders at three inter-grating distances ($2.24$, $4.24$ and 6.24 mm). Three regions according to the incident angle are selected and the corresponding real space fringes are compared in (b), (d) and (f).
Fig. 9.
Fig. 9. The macroscopic properties of the interference patterns are compared. (a) Comparison for the periods among theoretical prediction, measurement and simulation at eight different inter-grating distances from $2.24$ to 9.24 mm. (b) compares the averaged visibility for the first two orders over the full field of view between the experiment and simulation at different inter-grating distances.

Equations (12)

Equations on this page are rendered with MathJax. Learn more.

I res ( x ) I 0 = 1 + l = 1 3 V l cos ( 2 π l x p fr ) ,
p fr = R s + R g + R d R g p g ,
V = I max I min I max + I min .
u 1 ( x ) = exp [ i k ( x 2 / ( 2 R s ) ) ] R s ,
α 1 ( x ) = λ 2 π x φ 1 ( x ) .
t 1 ( x ) = i = 0 N 1 g 1 ( S ( i ; x ) , i ) .
u 1 ( x ) = u 1 ( x ) e i k δ Si t 1 ( x ) e k β Si t 1 ( x ) ,
v 1 ( x ) = u 1 ( x ) R s exp [ i k x 2 / ( 2 R s ) ] .
v 2 ( x ) = F 1 { F { Z P ( v 1 ( x ) ) } exp [ j π λ f x 2 ( R g / M 1 ) ] } ,
M 1 = R s + R g R s ,
u 2 ( x ) = 1 R s + R g v 2 ( x M 1 ) exp { i k x 2 / [ 2 ( R s + R g ) ] } .
I ( x + Δ s ) = 1 + l = 1 V l cos [ 2 π l 1 p fr ( x + Δ s ) ] ,
Select as filters


Select Topics Cancel
© Copyright 2024 | Optica Publishing Group. All rights reserved, including rights for text and data mining and training of artificial technologies or similar technologies.