Moiré artifacts reduction in Talbot-Lau X-ray phase contrast imaging using a three-step iterative approach

Siwei Tao; Yueshu Xu; Ling Bai; Zonghan Tian; Xiang Hao; Cuifang Kuang; Cuifang Kuang; Xu Liu; Xu Liu; Xu Liu

doi:10.1364/OE.466277

1. Introduction

In recent years, X-ray imaging has played a significant role in medical diagnostics and industrial inspection since it can provide internal structures of tested objects [1–3]. In general, X-ray absorption image is quite useful to distinguish strong absorbing objects. However, as for materials made of low electron density, such as tissues made of carbon, hydrogen and oxygen, the contrast of absorption is often low and internal structures are nearly indistinguishable. X-ray phase contrast imaging has been put forward to solve this problem [4], which is based on the fact that the parameter relates to phase shift is three orders of magnitude higher than that relates to absorption [5]. Talbot-Lau X-ray phase contrast imaging is one of the most promising methods to build an X-ray phase contrast imaging system at a relative low cost as it relaxes the restriction of highly coherent X-ray source [6, 7], which means even a conventional low brilliance X-ray tube can be utilized in the imaging system [8]. Three images including transmission, differential phase contrast (DPC) and dark-field (DF) can be retrieved simultaneously via this method [9]. It is reported that this method has already been conducted in medical diagnostics and non-destructive industrial imaging nowadays [3, 10], along with some commercialized system arises [11].

The common image retrieval method in Talbot-Lau X-ray phase contrast imaging can be divided into two types: single-shot and phase-stepping [8, 12–14]. A typical method of the single-shot type is the spatial harmonic method based on Fourier transforms [13]. The single-shot method can retrieve three images from two single raw image (the sample and the background). This fast method accelerates the imaging speed and lowers the dose at the cost of poor resolution and signal-to-noise ratio (SNR) compared to the phase-stepping method [14]. In that case, the phase-stepping method is usually performed in order to reconstruct images with high resolution and SNR. With the development of micro-nano processing industry, the period of the grating used in Talbot-Lau X-ray imaging system is often on the order of microns [15]. Precious motorized translation stages are needed to perform accurate phase-stepping as stepping errors will cause artifacts in the reconstructed three images. Hauke found that the dependency of artifacts due to the phase-stepping errors is mathematically derived by a Taylor series expansion and verified it via simulation [16]. Since a conventional X-ray tube is usually utilized in the Talbot-Lau X-ray phase contrast imaging system, the dose fluctuations is another factor that cannot be ignored along with the phase-stepping errors. In particular, continuous multiple exposures are inevitable during the phase-stepping procedure. The continuous exposures of X-ray tube will lead to dose fluctuations between each raw image and cause artifacts in the reconstructed transmission, DPC and DF images. Based on this problem, Hashimoto derived the effects of the phase-stepping errors and the dose fluctuations on DPC images [17]. However, the influence of the dose fluctuation on transmission and DF images is not derived and discussed in this passage.

In order to cut the expense on the high accuracy motorized translation stages and X-ray tube with extremely high stability, research towards the moiré artifact reduction via software approach has been carried out. Okada put forward an algorithm based on the least square method to obtain the phase distribution and scanning phase shift simultaneously in phase shift interferometer at first in 1991 [18]. Wang and Han came up with an advanced iterative algorithm (AIA) where inaccurate initial values of phase shift can be used for iterations to reduce artifacts [19]. Vargas proposed a fast moiré artifacts reduction method called the Generalization of the Principal Component Analysis algorithm (GPCA) which saved the processing time of the image reconstruction [20]. Kaeppler put forward an improved reconstruction algorithm of phase-stepping data with the optimization of the total variance of the image [21]. Hauke came up with an algorithm which utilized two different cost functions to optimize the DPC and DF image separately for moiré artifact suppression [22]. However, the algorithm above only considers the effects caused by phase-stepping errors. Hashimoto proposed an algorithm for image reconstructions of phase stepping data with stepping errors and dose fluctuations [17]. The artifacts are reduced and the precision of the reconstructed images is guaranteed but at the cost of more steps. The deep learning techniques have also been introduced in moiré artifact reduction recently. Chen used the CNN network to remove the moiré artifact in the reconstructed images [23], and a large number of data should be used as the train datasets to ensure the performance of the network.

In this work, the effects of the stepping errors and the dose fluctuations on the transmission, DPC and DF images are theoretically derived and systematically summarized. Moreover, a novel three-step least square iterative algorithm for the image reconstruction in Talbot-Lau interferometry with stepping errors and dose fluctuations is proposed. This algorithm is demonstrated to be effective in both numerical simulation and experiment. This work may be beneficial for the commercialization of X-ray phase contrast imaging towards clinic use since it relaxes the limitation of high precision motorized stages and high stability of X-ray tube.

2. Methods

The setup of a typical Talbot-Lau X-ray phase contrast imaging system is shown in Fig. 1. It consists of an X-ray tube, an object, a detector and three gratings (named G0, G1 and G2). Normally, the G1 works as a beam splitter. If the incident X-rays meet the necessary spatial coherence of the system, the self-image of the G1 will occur at periodic distance downstream of the G1, which is well-known as the Talbot effect. However, the spot size of the conventional X-ray tube usually cannot meet the necessary spatial coherence. Therefore, an absorption grating called G0 is needed to modulate the X-ray tube source into an array of individually coherent, but mutually incoherent line sources. Each line source acts as an individual X-ray source to illuminate the G1. The periodic distributions produced by each line source and the G1 will be incoherently superimposed. By reasonably setting the distances between the G0 to the G1 and the G1 to the detector, the fringe visibility will be greatly enhanced. As the period of the G1 used in this system is usually a few microns and the pixel size of a common X-ray detector is tens even hundreds of microns. An absorption grating G2 is placed after the G1 to generate large period fringes, which can be captured by the detector. In order to generate moiré fringes with large fringe period, the period of the G2 is designed the same as the period of the periodic distributions generated by the G1 at the position of the detector. And the G2 is always placed a little tilted relative to the G1.

Fig. 1. Schematic of a typical Talbot-Lau X-ray phase contrast imaging system.

Download Full Size | PDF

When an object is placed into the beam path, the X-ray photon distributions received by the detector will change due to the absorption, diffraction and scattering of X-ray when it interacts with the sample. However, these three signals are mixed and stacked together in the raw image, which can be expressed as:

(1)$${\rm I}({x,y} )= {I_0}({x,y} )+ M({x,y} )\cdot \cos ({\varphi ({x,y} )} ),$$

here, ${I_0}({x,y} )$ represents the intensity distributions, $M({x,y} )$ the modulation distributions and $\varphi ({x,y} )$ the phase distributions [19]. In order to solve these variables to separate absorption, diffraction and scattering signals, phase-stepping is performed. It is clear that at least three steps are needed for solving three parameters in the conventional phase-stepping procedure [7]. At the same time, another three steps of the background should be performed to retrieve transmission, DPC and DF of the tested sample. The X-ray distributions measured by the detector when performing phase-stepping can be expressed as:

(2)$$I({x,y,n} )= {I_0}({x,y} )+ M({x,y} )\cdot \cos \left( {\Phi ({x,y} )+ \frac{{2\pi n}}{N}} \right),$$

where N is the number of steps used in phase-stepping, ${{2\pi n} / N}$ denotes the phase steps with the n ranging from 0 to $N - 1$. We define ${L_1}({{\rm I}({x,y} )} )$ and ${L_2}({{\rm I}({x,y} )} )$ as:

(3)$${L_1}({{{\mathbf {I}}}({x,y} )} )= \frac{1}{N}\sum\limits_{n = 0}^{N - 1} {I({x,y,n} )} = {I_0}({x,y} ),$$

(4)$${L_2}({{\mathbf {{I}}}({x,y} )} )= \frac{2}{N}\sum\limits_{n = 0}^{N - 1} {I({x,y,n} ){e^{ - \frac{{2\pi n}}{N}i}}} = M({x,y} ){e^{i\Phi ({x,y} )}},$$

here, ${\mathbf {{I}}}({x,y} )$ denotes $I({x,y,0} )$ to $I({x,y,N - 1} )$. In that case, the intensity, visibility and phase distributions can be extracted from

(5)$${I_0}({x,y} )= {L_1}({{\mathbf {{I}}}({x,y} )} ),$$

(6)$$V({x,y} )= \frac{{|{{L_2}({{\mathbf {{I}}}({x,y} )} )} |}}{{{L_1}({{\mathbf {\textrm {I}}}({x,y} )} )}},$$

(7)$$\Phi ({x,y} )= \arg ({{L_2}({{\mathbf {{I}}}({x,y} )} )} ).$$

However, due to the low accuracy of the motorized translation stage, a variable $\delta (n )$ representing phase-stepping error at each step should be considered in the calculation. Also, when conducting continuous exposures, the temperature of the X-ray tube will continue to rise. The stability of X-ray tube cannot be ensured due to the increase of temperature, which results in the dose fluctuations during exposures. For example, when the filament voltage is set constant during the exposure, the resistance of the filament circuit increases with the temperature. This will cause the reduction in filament current and tube current. In that case, the X-ray intensity will attenuate during the exposure. The dose fluctuations usually vary in different X-ray tubes and the physical model of dose fluctuations is complicated since it relates to various factors. It is appropriate to use the Taylor series expansion (TSE) to fit the dose fluctuations. As the intensity and modulation amplitude usually do not change between frame to frame, the light distribution at the detector for each step can be expressed as

(8)$$I({x,y,n} )= \sum\limits_{i = 0}^k {{a_i}{n^i}} \left[ {{I_0}({x,y} )+ M({x,y} )\cdot \cos \left( {\Phi ({x,y} )+ \frac{{2\pi n}}{N} + \delta (n )} \right)} \right],$$

in which ${a_i}$ is the coefficient which is related to dose fluctuations. In order to evaluate the effects caused by the phase-stepping errors and the dose fluctuation, the TSE is also used for analysis. The errors caused by the above two factors is

(9)$$\Delta {e_{Z\_factor}}({x,y} )= \sum\limits_{n = 0}^{N - 1} {\frac{{\partial Z({x,y} )}}{{\partial \delta (n )}}\Delta \delta (n )} + \sum\limits_{n = 0}^{N - 1} {\frac{{\partial Z({x,y} )}}{{\partial \mu (n )}}} \Delta \mu (n ),$$

here, $Z({x,y} )\in ({{I_0}({x,y} ),V({x,y} ),\Phi ({x,y} )} )$, $\mu (n )= \sum\limits_{i = 0}^k {{a_i}{n^i}}$. The detailed calculations of the impacts on the reconstructed images caused by the above two factors are illustrated in Supplement 1. The results are shown as follows:

(10)$$\Delta {e_{{I_0}\_step}}({x,y} )={-} \frac{{M({x,y} )}}{N}\sum\limits_{n = 0}^{N - 1} {\sin \left( {\Phi ({x,y} )+ \frac{{2\pi n}}{N}} \right)} \Delta \delta (n ),$$

(11)$$\Delta {e_{{I_0}\_dose}}({x,y} )= \frac{1}{N}\sum\limits_{n = 0}^{N - 1} {\left( {{I_0}({x,y} )+ M({x,y} )\cdot \cos \left( {\Phi ({x,y} )+ \frac{{2\pi n}}{N}} \right)} \right)} \Delta \mu (n ),$$

(12)$$\Delta {e_{\Phi \_step}}({x,y} )= \frac{1}{N}\sum\limits_{n = 0}^{N - 1} {\left( {1 - \cos \left( {2\left( {\Phi ({x,y} )+ \frac{{2\pi n}}{N}} \right)} \right)} \right)} \Delta \delta (n ),$$

(13)$$\Delta {e_{\Phi \_dose}}({x,y} )={-} \frac{1}{N}\sum\limits_{n = 0}^{N - 1} {\left( {\frac{{2{I_0}({x,y} )}}{{M({x,y} )}}\sin \left( {\Phi ({x,y} )+ \frac{{2\pi n}}{N}} \right) - \sin \left( {2\left( {\Phi ({x,y} )+ \frac{{2\pi n}}{N}} \right)} \right)} \right)} \Delta \mu (n ),$$

(14)$$\Delta {e_{V\_step}}({x,y} )= \frac{1}{N}\sum\limits_{n = 0}^{N - 1} {\left( {{{\left( {\frac{{M({x,y} )}}{{{I_0}({x,y} )}}} \right)}^2}\sin \left( {\Phi ({x,y} )+ \frac{{2\pi n}}{N}} \right) - \frac{{M({x,y} )}}{{{I_0}({x,y} )}}\sin \left( {2\left( {\Phi ({x,y} )+ \frac{{2\pi n}}{N}} \right)} \right)} \right)} \Delta \delta (n ),$$

(15)$$\Delta {e_{V\_dose}}({x,y} )= \frac{1}{N}\sum\limits_{n = 0}^{N - 1} {\left( {\left( {2 - {{\left( {\frac{{M({x,y} )}}{{{I_0}({x,y} )}}} \right)}^2}} \right)\cos \left( {\Phi ({x,y} )+ \frac{{2\pi n}}{N}} \right) + \frac{{M({x,y} )}}{{{I_0}({x,y} )}}\cos \left( {2\left( {\Phi ({x,y} )+ \frac{{2\pi n}}{N}} \right)} \right)} \right)} \Delta \mu (n ).$$

2.1 Phase distributions calculation

Three steps of the least square method are conducted and iterative process are introduced for accurate calculation of each parameter in Eq. (8) to reduce the moiré artifacts in the reconstructed images. We first start with the phase distributions calculation of the tested sample. Here, we define ${A_1}({x,y} )= M({x,y} )\cos ({\Phi ({x,y} )} )$, ${A_2}({x,y} )={-} M({x,y} )\cos ({\Phi ({x,y} )} )$ and Eq. (8) can be written as

(16)$$I({x,y,n} )= \sum\limits_{i = 0}^k {{a_i}{n^i}} \left[ {{I_0}({x,y} )+ {A_1}({x,y} )\cdot \cos \left( {\frac{{2\pi n}}{N} + \delta (n )} \right) + {A_2}({x,y} )\cdot \sin \left( {\frac{{2\pi n}}{N} + \delta (n )} \right)} \right].$$

For the sake of subsequent iteration, the loss function is set as the root mean square error (RMSE), which is

(17)$$RMSE = \sqrt {\frac{1}{P}\frac{1}{Q}\sum\limits_{x = 1}^P {\sum\limits_{y = 1}^Q {{{({{I^{exp }}({x,y,n} )- {I^{recon}}({x,y,n} )} )}^2}} } } ,$$

where $I_n^{exp }({x,y} )$ is the X-ray photon distributions measured in the experiment, $I_n^{r\textrm{econ}}({x,y} )$ is the reconstructed X-ray photon distributions by calculation. The least square method is adopted for calculation of phase distributions. The solution to unknowns in Eq. (16) is

(18)$$\begin{array}{l} \left[ {\begin{array}{*{20}{c}} {{I_0}({x,y} )}\\ {{A_1}({x,y} )}\\ {{A_2}({x,y} )} \end{array}} \right] = \frac{1}{{\sum\limits_{i = 0}^k {{a_i}{n^i}} }} \cdot Matrix_1^{ - 1} \times \left[ {\begin{array}{*{20}{c}} {\sum\limits_{n = 0}^{N - 1} {{I^{real}}({x,y,n} )} }\\ {\sum\limits_{n = 0}^{N - 1} {{I^{real}}({x,y,n} )\cos \left( {\frac{{2\pi n}}{N} + \delta (n )} \right)} }\\ {\sum\limits_{n = 0}^{N - 1} {{I^{real}}({x,y,n} )\sin \left( {\frac{{2\pi n}}{N} + \delta (n )} \right)} } \end{array}} \right].\\ Matri{x_1} = {\left[ {\begin{array}{*{20}{c}} {\sum\limits_{n = 0}^{N - 1} 1 }&{\sum\limits_{n = 0}^{N - 1} {\cos \left( {\frac{{2\pi n}}{N} + \delta (n )} \right)} }&{\sum\limits_{n = 0}^{N - 1} {\sin \left( {\frac{{2\pi n}}{N} + \delta (n )} \right)} }\\ {\sum\limits_{n = 0}^{N - 1} {\cos \left( {\frac{{2\pi n}}{N} + \delta (n )} \right)} }&{\sum\limits_{n = 0}^{N - 1} {{{\cos }^2}\left( {\frac{{2\pi n}}{N} + \delta (n )} \right)} }&{\sum\limits_{n = 0}^{N - 1} {\cos \left( {\frac{{2\pi n}}{N} + {\delta_n}} \right)\sin \left( {\frac{{2\pi n}}{N} + \delta (n )} \right)} }\\ {\sum\limits_{n = 0}^{N - 1} {\sin \left( {\frac{{2\pi n}}{N} + \delta (n )} \right)} }&{\sum\limits_{n = 0}^{N - 1} {\cos \left( {\frac{{2\pi n}}{N} + \delta (n )} \right)\sin \left( {\frac{{2\pi n}}{N} + \delta (n )} \right)} }&{\sum\limits_{n = 0}^{N - 1} {{{\sin }^2}\left( {\frac{{2\pi n}}{N} + \delta (n )} \right)} } \end{array}} \right]^{ - 1}} \end{array}$$

From this calculation, the phase, intensity and modulation distributions can be obtained simultaneously. The phase distributions are retrieved as

(19)$$\Phi ({x,y} )= {\tan ^{ - 1}}\left( { - \frac{{{A_2}({x,y} )}}{{{A_1}({x,y} )}}} \right).$$

And the modulation distributions are determined as

(20)$$M({x,y} )= \sqrt {{{({{A_1}({x,y} )} )}^2} + {{({{A_2}({x,y} )} )}^2}} .$$

It is obvious that when the phase-stepping errors are set as 0 and dose fluctuations are ignored, the formula for calculating phase distributions is the same as the traditional phase-stepping retrieval method.

2.2 Phase-stepping errors calculation

After obtaining the phase distributions, the phase-stepping errors can be calculated in similar way as above. It is also necessary to define ${B_1}(n )= M({x,y} )\cos ({{{2\pi n} / N} + \delta (n )} )$, ${B_2}(n )={-} M({x,y} )\sin ({{{2\pi n} / N} + \delta (n )} )$. In that case, Eq. (8) can also be written as:

(21)$$I({x,y,n} )= \sum\limits_{i = 0}^k {{a_i}{n^i}} [{{I_0}(n )+ {B_1}(n )\cdot \cos ({\Phi ({x,y} )} )+ {B_2}(n )\cdot \sin ({\Phi ({x,y} )} )} ].$$

The loss function is also chosen as the RMSE. With the performance of the least square method again, phase-stepping errors can be obtained. It is assumed that the image size of raw data is $P \times Q$. The unknowns in Eq. (8) are calculated via

(22)$$\begin{array}{l} \left[ {\begin{array}{*{20}{c}} {{I_0}(n )}\\ {{B_1}(n )}\\ {{B_2}(n )} \end{array}} \right] = \frac{1}{{\sum\limits_{i = 0}^k {{a_i}{n^i}} }} \cdot Matrix_2^{ - 1} \times \left[ {\begin{array}{*{20}{c}} {\sum\limits_{x = 1}^P {\sum\limits_{y = 1}^Q {{I^{real}}({x,y,n} )} } }\\ {\sum\limits_{x = 1}^P {\sum\limits_{y = 1}^Q {{I^{real}}({x,y,n} )} \cos ({\Phi ({x,y} )} )} }\\ {\sum\limits_{x = 1}^P {\sum\limits_{y = 1}^Q {{I^{real}}({x,y,n} )} \sin ({\Phi ({x,y} )} )} } \end{array}} \right].\\ Matri{x_2} = {\left[ {\begin{array}{*{20}{c}} {\sum\limits_{x = 1}^P {\sum\limits_{y = 1}^Q 1 } }&{\sum\limits_{x = 1}^P {\sum\limits_{y = 1}^Q {\cos ({\Phi ({x,y} )} )} } }&{\sum\limits_{x = 1}^P {\sum\limits_{y = 1}^Q {\sin ({\Phi ({x,y} )} )} } }\\ {\sum\limits_{x = 1}^P {\sum\limits_{y = 1}^Q {\cos ({\Phi ({x,y} )} )} } }&{\sum\limits_{x = 1}^P {\sum\limits_{y = 1}^Q {{{\cos }^2}({\Phi ({x,y} )} )} } }&{\sum\limits_{x = 1}^P {\sum\limits_{y = 1}^Q {\cos ({\Phi ({x,y} )} )\sin ({\Phi ({x,y} )} )} } }\\ {\sum\limits_{x = 1}^P {\sum\limits_{y = 1}^Q {\sin ({\Phi ({x,y} )} )} } }&{\sum\limits_{x = 1}^P {\sum\limits_{y = 1}^Q {\cos ({\Phi ({x,y} )} )\sin ({\Phi ({x,y} )} )} } }&{\sum\limits_{x = 1}^P {\sum\limits_{y = 1}^Q {{{\sin }^2}({\Phi ({x,y} )} )} } } \end{array}} \right]^{ - 1}} \end{array}$$

After the calculation, the phase-stepping errors can be retrieved via

(23)$$\delta (n )= {\tan ^{ - 1}}\left( { - \frac{{{B_2}(n )}}{{{B_1}(n )}}} \right).$$

2.3 Dose fluctuation coefficients calculation

As the initial values of intensity, phase, modulation distributions and phase-stepping errors are calculated, the coefficients of dose fluctuations can be calculated via the third step of the least square method. Here, we define the reconstructed X-ray photon intensity of a single pixel without dose fluctuations as:

(24)$${I^{wodf}}({x,y,n} )= {I_0}({x,y} )+ M({x,y} )\cdot \cos \left( {\Phi ({x,y} )+ \frac{{2\pi n}}{N} + \delta (n )} \right).$$

The calculation formula for the dose fluctuation coefficients of each pixel is given as below:

(25)$$\begin{array}{l} \left[ {\begin{array}{*{20}{c}} {{a_0}({x,y} )}\\ {{a_1}({x,y} )}\\ \ldots \\ {{a_k}({x,y} )} \end{array}} \right] = Matrix_3^{ - 1} \times \left[ {\begin{array}{*{20}{c}} {\sum\limits_{n = 0}^{N - 1} {{I^{real}}({x,y,n} )} }\\ {\sum\limits_{n = 0}^{N - 1} {{I^{real}}({x,y,n} )\cdot n} }\\ \ldots \\ {\sum\limits_{n = 0}^{N - 1} {{I^{real}}({x,y,n} )\cdot {n^k}} } \end{array}} \right].\\ Matri{x_3} = \left[ {\begin{array}{*{20}{c}} {\sum\limits_{n = 0}^{N - 1} {{n^k} \cdot {I^{wodf}}({x,y,n} )} }& \ldots &{\sum\limits_{n = 0}^{N - 1} {n \cdot {I^{wodf}}({x,y,n} )} }&{\sum\limits_{n = 0}^{N - 1} {{I^{wodf}}({x,y,n} )} }\\ {\sum\limits_{n = 0}^{N - 1} {{n^{k + 1}} \cdot {I^{wodf}}({x,y,n} )} }& \ldots &{\sum\limits_{n = 0}^{N - 1} {{n^2} \cdot {I^{wodf}}({x,y,n} )} }&{\sum\limits_{n = 0}^{N - 1} {n \cdot {I^{wodf}}({x,y,n} )} }\\ \ldots & \ldots & \ldots & \ldots \\ {\sum\limits_{n = 0}^{N - 1} {{n^{2k}} \cdot {I^{wodf}}({x,y,n} )} }& \ldots &{\sum\limits_{n = 0}^{N - 1} {{n^{k + 1}} \cdot {I^{wodf}}({x,y,n} )} }&{\sum\limits_{n = 0}^{N - 1} {{n^k} \cdot {I^{wodf}}({x,y,n} )} } \end{array}} \right] \end{array}$$

In order to minimize the random errors of the dose fluctuation coefficients, the final dose fluctuation coefficients are the average of the calculated dose fluctuation coefficients for each pixel.

2.4 Iteration

The initial values of coefficients ${a_i}$ should be first calculated using polynomial fitting with the time and the mean value of each raw image. Then, the phase distributional and phase-stepping errors can be calculated sequentially. The optimization iterations start after the acquisition of these initial values. A single optimization iteration cycle includes the calculation of phase distributions, phase-stepping errors and dose fluctuation coefficients. The least square method is applied to each calculation procedure. This is the reason why this approach is called the three-step least square iterative algorithm. The convergence criteria of this iteration are given as:

(26)$$\left\{ {\begin{array}{l} {|{({{\delta^{(j)}}(n )- {\delta^{(j)}}(1 )} )- ({{\delta^{(j - 1)}}(n )- {\delta^{(j - 1)}}(1 )} )} |< {\varepsilon_1}}\\ {\left|{\sum\limits_{i = 0}^k {{a_i}^{(j)}{n^i}} - \sum\limits_{i = 0}^k {{a_i}^{(j - 1)}{n^i}} } \right|< {\varepsilon_2}} \end{array}} \right.,$$

where $(j)$ indicates the jth iteration. ${\varepsilon _1}$ and ${\varepsilon _2}$ are the given accuracies (for example, ${\varepsilon _1} = 1 \times {10^{ - 5}}, {\varepsilon _2} = 1 \times {10^{ - 5}}$). The iteration will continue unless the convergence criteria are met. The given accuracies are set according to the requirements for the precision of the reconstruction results. Higher accuracies result in precious results but longer calculation time and vice versa.

According to the principle above, the whole algorithm can be expressed as Algorithm 1:

oe-30-20-35096-i001

This method should be applied to both raw images with the objects and the background. Finally, the transmission $T\textrm{ra}$, differential phase contrast $Dpc$ and dark-field $Dfc$ images can be retrieved via the formulas below:

(27)$$Tra = \frac{{I_0^w(x,y)}}{{I_0^{wo}(x,y)}},$$

(28)$$Dpc = {\Phi ^w}({x,y} )- {\Phi ^{wo}}({x,y} ),$$

(29)$$Dfc = \frac{{{M^w}({x,y} )\cdot {I^{wo}}({x,y} )}}{{{M^{wo}}({x,y} )\cdot {I^w}({x,y} )}}.$$

Here, parameters with the superscript w represents the corresponding parameters calculated from raw images with objects. On the contrary, those parameters with the superscript wo indicates the corresponding parameters calculated from raw images without the objects. Regardless of the dose fluctuations and phase-stepping errors, this algorithm is basically the same as the conventional phase-stepping retrieval method. In that case, we believe that the three-step least square iterative algorithm is more universal for the image reconstructions.

3. Simulation conditions and results

The simulation was conducted in order to verify the feasibility and accuracy of the proposed algorithm. It is based on the angular spectrum algorithm, and detailed principles and processes are mentioned in our previous study [24]. The whole forward projection simulation processes are concluded as: free space propagation from the G0 to the G1, interactions with the G1, free space propagation from the G1 to the sample, interactions with the sample. free space propagation from the sample to the G2, interactions with the G2, pixel binning to generate raw image captured by the detector. The X-ray source is simulated as a coherent spherical wave because the line source generated by G0 can be considered coherent. As an X-ray tube with Mo target is used in our experiment, a polychromatic simulation was conducted. The X-ray spectrum data in our experiment was tested by a Si-PIN X-ray detector (Ametek, USA). The Si-PIN X-ray detector was placed after three gratings to record the energy and count of photons. The dose fluctuations in our experiment were also tested, which are illustrated in detail in the experiment results and discussions part. It was proved that a quadratic model was suitable to represent the dose attenuation in our experiment. Therefore, a quadratic dose attenuation model was also chosen in the simulation. Three dose fluctuation coefficients named ${a_0}$, ${a_1}$ and ${a_2}$ were utilized and set as 0.006, -0.1 and 1 separately. It meant that the normalized average light intensity would decrease from 100% to 69.6% in the simulation, as the phase-stepping was performed. The dose fluctuation in the simulation was set larger than in the experiment to prove the robustness of the proposed algorithm. Gratings named G0, G1 and G2 were all chosen as 6µm with a duty cycle of 0.5, which was also similar to the experiment conditions. The G0 and the G2 were both ideal absorption gratings and the G1 was an ideal π-phase grating at the photon energy of 17.48keV, which is the kα of the Molybdenum target. In the simulation, the G0 is completely parallel to the G1. And the G2 was a little tilted compared to the G1 in order to generate fringes with large fringe period. A chest X-ray image from the National Institutes of Health Clinical Center was considered as the sample in the simulation [25]. The pixel values of this image were redefined in order to create a suitable simulation conditions for evaluating the performance of the proposed algorithm. The distance between G1 and G2 was set as 38.1cm, which was the 2^nd Talbot distance for G1 illuminated by X-ray with the photon energy of 17.48keV. In order to guarantee the Lau effect was introduced into the system, the G1-to-G2 distance was also set as 38.1cm, which could improve the visibility of moiré fringes. The tested object was placed 0.1cm behind G1 to achieve high contrast in the DPC image. The field of view (FOV) of the ideal detector was simulated as $40.96cm \times 40.96cm$ with the pixel size of $0.8mm \times 0.8mm$. The number of steps for phase-stepping was chosen as 5 in the simulation. The phase-stepping error was randomly set in each step, but would be no more than 100% of the step size. The given accuracies for iteration were set as ${\varepsilon _1} = 1 \times {10^{ - 6}}$ and ${\varepsilon _2} = 1 \times {10^{ - 6}}$. The Gaussian noise with standard deviation $\sigma = 3$ and mean value $\mu = 0$ and Poisson noise were both added into the raw images in order to simulate the noise conditions in the experiment.

According to the X-ray spectrum data from the Si-PIN X-ray detector, photons with photon energy ranging from 4keV to 40keV was considered. The raw images captured by the detector in the simulation can be calculated by

(30)$$I_{all}^{sim}({x,y} )= \sum\limits_{i = 1}^N {{w_i}I_i^{sim}({x,y} )} ,$$

here, $I_i^{sim}({x,y} )$ is the monochromatic simulated raw images obtained from the detector for a certain photon energy. ${w_i}$ is the weight of the photon energy, which is determined by the X-ray spectrum data. The simulation was conducted and three-step least square iterative algorithm were performed to retrieved the images. The initial values of transmission, differential phase contrast and dark-field were obtained via the traditional phase-stepping retrieval algorithm. Also, the initial phase-stepping errors for iterations were set as 0. As for the dose fluctuation coefficients, the average light intensities of images for each step were calculated and the initial values of three coefficients were obtained by fitting the average light intensity value with the time coordinate. In order to elucidate that both phase-stepping errors and dose fluctuations will cause moiré artifacts in the reconstructed three images, traditional retrieval algorithm, retrieval algorithm only considers phase-stepping error, retrieval algorithm only considers the dose fluctuations and the proposed three-step least square iterative algorithm were all performed for comparison. The retrieved images under these four algorithms and the ground truth are shown in Fig. 2. It can be inferred from Fig. 2 that the presence of the phase-stepping errors and the dose fluctuations will both cause artifacts in the reconstructed transmission, DPC and DF images. The fringe frequency of the moiré artifacts caused by the phase-stepping errors or the dose fluctuations are the same as the captured moiré fringes in transmission image. As for DPC and DF images with the influence of the phase-stepping errors, the fringe frequency of the moiré artifacts will double. However, for these with the influence of the dose fluctuations, the fringe frequency remains the same. This is mainly due to the coefficients of the two different frequency terms. As the preset ${I_0}({x,y} )$ is larger than $M({x,y} )$ in the simulation, the above phenomenon can be explained, which is consistent with the theory. The artifacts caused by dose fluctuations are not obvious in the transmission image, but in DPC and DF images, the artifacts still cannot be ignored. It can be concluded that the artifacts cannot be eliminated if phase-stepping errors or dose fluctuations are considered only. The above results demonstrate that the proposed algorithm is feasible for moiré artifacts reduction in phase-stepping method.

The accuracy of the proposed algorithm was validated by comparing the calculated results of unknown parameters with the preset values. The G2 was set at five different positions for five steps, which means it will be moved four times with different step size. Therefore, parameters including four phase-stepping errors and three dose fluctuation coefficients in the simulations are compared. The comparisons are listed in Table 1, which include both the tested sample and the background. It can be inferred from the relative errors (RE) that the calculated results are pretty close to the preset values, which indicates that the proposed algorithm is of high precision in the simulation. The transmission, DPC and DF images retrieved via the three-step least square iterative algorithm are compared with the ground truth. There are little difference between 2D images of the retrieved images and ground truth just as Figs. 2(j)-(o) shows, the 1D profiles of the row 10 are chosen for comparison. Figures 3(a)-(c) shows the differences between the transmission, DPC and DF images retrieved from the proposed algorithm and the ground truth separately. It is clear that periodic fringes can hardly be seen in the profile of three reconstructed images. The profile of the reconstructed images is highly coincident with the ground truth, which further demonstrates that the moiré artifacts caused by phase-stepping errors and dose fluctuations are mostly eliminated.

Fig. 2. The reconstructed transmission, DPC and DF images via four different image reconstruction algorithms and the ground truth. (a)-(c) From the traditional phase-stepping image reconstruction algorithm. (d)-(f) From the algorithm which only considers the dose fluctuations. (g)-(i) From the algorithm which only considers the phase-stepping errors. (j)-(l) From the proposed three-step least square iterative algorithm. (m)-(o) The ground truth. The green and red dotted lines indicate row 10.

Download Full Size | PDF

Fig. 3. Some important figures in the simulation which shows the performance of the proposed algorithm. (a) The profile of row 10 in the reconstructed transmission image and the ground truth. (b) The profile of row 10 in the reconstructed DPC image and the ground truth. (c) The profile of row 10 in the reconstructed DF image and the ground truth. (d) The stepping curves for pixel (300,500), including only the background and the sample placed in the beam path. (e) The average relative error of four steps which caused by different types of noise. (f) The average relative error of four steps which caused by different standard deviation of Gaussian noise.

Download Full Size | PDF

Table 1. Comparison towards calculated values and preset parameters in the simulation

View Table

The stepping curves for the dose fluctuation corrected pixel (300,500) are shown in Fig. 3(d). It can be inferred that the five points are able to be fitted into a sinusoid curve, which further demonstrates that the proposed algorithm is of high accuracy. Admittedly, there are some errors between the calculated results and the preset values of the parameters in Table 1. In order to study the cause of these errors, different types and intensities of noise are used in the simulation. The Gaussian noise and the Poisson noise are two common types of noise in the X-ray imaging system. In that case, the simulation was first conducted with no noise, only the Poisson noise, only the Gaussian noise ($\mu = 0, \sigma = 3$) and both two noises. The average relative errors of 4 steps are calculated and shown in Fig. 3(e). It is clear that the Gaussian noise contributes more to these errors when the object is placed in the system. While for the background, the influence of these two noise is almost the same. This is because the object absorbs the X-ray and the photons acquired by the detector decreases, which reduces the intensity of Poisson noise. Also, different intensities ($\sigma = 1, \sigma = 2, \sigma = 4, \sigma = 8$) of Gaussian noise were added into the raw image for comparison and the results are displayed in Fig. 3(f). It is obvious that the relative error increase dramatically with the increase of Gaussian noise intensity, which means the SNR of the acquired images is important to guarantee the precision of this algorithm.

4. Experiment results and discussions

In order to verify the effectiveness of the algorithm in the experimental conditions, experiments were carried out in our in-house Talbot-Lau X-ray phase contrast imaging system. The hardware system in this experiment consists of a rotating anode Molybdenum target diagnostic-level X-ray tube (E7290AX, Toshiba Rot anode, Japan), three gratings (G0: 6µm absorption grating, with duty cycle 0.5. G1: 6µm π-phase grating, with duty cycle 0.5. G2: 6µm absorption grating, with duty cycle 0.5. Microworks precision structures, Germany), a motorized translation stage (DDS220/M, Thorlabs, USA) and a CMOS flat-panel X-ray detector with an active area of $57mm \times 64mm$ and the pixel size 49.5µm (Shad-o-Box 1K HS, Teledyne Digital Imaging Inc, Canada). The conventional X-ray tube was operated at 40kV. The G0 was placed close to the X-ray tube. The distance between the G0 and the G1 was set as 38.1cm. And the distance between the G1 and the G2 was also 38.1cm, which was the 2^nd Talbot distance. The X-ray detector was placed close to the G2. The tested samples were placed between the G1 and the G2. The G0 and the G1 were set paralleled, while the G2 was a little tilted compared with the G1 to generate large period fringes which could be captured by the X-ray detector. Also, five steps were adopted for the retrieval of three images, which was the same as those in the simulation. As the minimum step size of our motorized translation stage is 0.05µm, the phase-stepping error was introduced artificially by setting different step size for each step. And each step was performed at the same time interval to insure that the raw data would fit into the dose fluctuation model. The exposure time of each raw image was 3s. The tested samples in the experiments were a hornet specimen and a Gephyrocharax melanocheir specimen, and both were scanned and calculated separately.

The acquisitions of the initial values were the same as those in the simulation, and the iterative algorithm was applied with ${\varepsilon _1} = 1 \times {10^{ - 6}}$ and ${\varepsilon _2} = 1 \times {10^{ - 6}}$. The performance of the three-step least square iterative algorithm is shown in Fig. 4. Herein, the image reconstruction method which only considered the phase-stepping error was also performed to illustrate the influence caused by the dose fluctuations. The reconstructed three images using this algorithm are displayed in Figs. 4(d)-(f) and the three-step iterative algorithm in Figs. 4(e)-(g). The influence caused by dose fluctuations in Figs. 4(d)-(f) is consistent with the simulation results and the theory, which demonstrates that the dose fluctuations cannot be ignored. The average RMSE of five step raw images for algorithm only considers the phase-stepping error is 0.0327, while 0.0051 for the proposed method. This result further illustrates the introduction of dose fluctuation parameters do improve the reconstruction precision. The reconstructed three images of the Gephyrocharax melanocheir specimen are demonstrated in Figs. 5(a)-(c). Abundant information about this fish can be obtained from the reconstructed transmission, DF and DPC images. Figures 5(d)-(f) demonstrates the profiles of row 700 of the reconstructed images, which is a flat area of the image. It can be inferred from the profiles that the periodic fluctuations are mostly reduced via the proposed algorithm, which demonstrates the effectiveness of the algorithm. Also, the reduced moiré artifacts in both 2D images of the hornet specimen and the Gephyrocharax melanocheir specimen further illustrates that the proposed algorithm is still feasible in the experiment. At the same time, the high precision of the proposed method is also proved. As stated above, the phase-stepping errors were artificially introduced. Therefore, it is easy to get the theoretical phase stepping errors. The relative errors between the theoretical values and the calculated results for the hornet specimen and the Gephyrocharax melanocheir specimen are shown in Fig. 6(a). It can be seen that the relative error in each step is higher than that in the numerical simulations. The main cause for this problem is the accuracy of our motorized stage. As the resolution of this stage is 0.05µm, which means there are still some errors between the preset theoretical values and the actual step values that cannot be easily measured. The stepping curves for the dose fluctuation corrected pixel (500,500) of both the hornet specimen and the Gephyrocharax melanocheir specimen in Figs. 6(b), (c), which shows that the five chosen points can be fitted into a sinusoid curve.

Fig. 4. The performance of the three-step least squares iterative algorithm and the influence of the dose fluctuations. (Top: Transmission images. Middle: DPC images. Bottom: DF images.) (a)-(c) The retrieved images via the traditional image reconstruction algorithm in phase stepping. (d)-(f) The retrieved images with image reconstruction algorithm only considering phase-stepping errors, but the dose fluctuations are not considered. It is obvious that the artifacts cannot be eliminated. (g)-(i) The retrieved images via the three-step least squares iterative algorithm.

Download Full Size | PDF

Fig. 5. The performance of the three-step least squares iterative algorithm for the Gephyrocharax melanocheir specimen. (a) Transmission image. (b) DPC image. (c) DF image. (For (a)-(c), the left part of the reconstructed image is from the traditional algorithm, and the right part from the proposed algorithm). Profiles of row 700 in the reconstructed images from the traditional algorithm and the proposed algorithm: (d) In Transmission. (e) In DPC. (f) In DF.

Download Full Size | PDF

Fig. 6. Some figures which illustrates the performance of the proposed three-step iterative algorithm. (a) The relative error between the theoretical values and the calculated results for the hornet specimen and the Gephyrocharax melanocheir (fish) specimen. The stepping curves for pixel (500,500), including only the background and the sample placed in the beam path for (b) the hornet specimen and (c) the Gephyrocharax melanocheir specimen. The average RMSE of different algorithms for (d) the hornet specimen and (e) the Gephyrocharax melanocheir specimen. (f) The dose fluctuation curve in our experiment, which shows that a quadratic fit is proper in our setup.

Download Full Size | PDF

The proposed algorithm is also compared with AIA [19], GPCA [20], Kaeppler’s [21] and Hauke’s [22] algorithm. As the train datasets is an important factor which influences the reconstructed results of Chen’s algorithm [23] and the minimum step requirement for Hashimoto’s algorithm is seven steps for simultaneous estimation of both the stepping errors and the dose fluctuations [17]. These two algorithms are not compared with the proposed algorithm in this paper. The average RMSE of the five step raw images which is calculated using Eq. (17) are used for comparison. The Hauke’s algorithm was performed with two loss functions to optimize DPC and DF images separately. Figure 6(d) demonstrates the average RMSE for the hornet specimen and its background. And Fig. 6(e) for the Gephyrocharax melanocheir specimen and its background. Here, Hauke-DPI refers to the calculated RMSE using Hauke’s algorithm and the loss function to optimize the DPC image. Similarly, Hauke-DFI refers to Hauke’s algorithm with the loss function for DF image optimization. All five algorithms are proved to be effective for moiré artifact reduction in our experiment. However, there are still some visible artifacts in the reconstructed images using AIA and GPCA. The Kaeppler’s algorithm is pretty effective for smooth samples as this algorithm seeks to minimize the optimization problem under total variation. In that case, this algorithm prefers smooth solutions. Some inner part of the hornet in this experiment scattered most of X-rays and the fringe visibility drops to nearly zero in these areas. This leads to large variations in these areas in DF images, which is the cause of high RMSE for Kaeppler’s algorithm in hornet specimen as it searches for the phase-stepping errors which leads to minimum total variation of the transmission and the visibility. The average RMSE of Hauke’s algorithm remains high with either of the two loss functions. This is due to the fact that the Hauke-DPI can reduce or eliminate moiré artifact in the DPC image, but the moiré artifact in the DF image under this loss function may still be remaining. But with the combination of the DPC image in Hauke-DPI and the DF image in Hauke-DFI, this algorithm is still feasible in moiré artifacts reduction. It can be inferred from Fig. 6 that the accuracy of the proposed algorithm is the highest among these algorithm, which demonstrates its high precision in the experiment.

As for the dose fluctuation, experiments were carried out for validation. The experimental setup consists of only the X-ray tube and the X-ray detector. 300 images with the exposure time of 0.2s were captured continuously. By calculating the average light intensity of each image, the dose fluctuation curve can be acquired and is shown in Fig. 6(f). The dose fluctuation verified experiments were performed three times to avoid random effects. It is clear that the X-ray dose attenuated with the time in our experiment. The cause for the continuous dose attenuation is the increasing temperature of the X-ray tube. The equivalent resistance in the filament circuits also increases with the temperature. As DC constant voltage source was utilized to power the filament in our experiment, the filament current would decrease accompanied by the increase of temperature. Therefore, the free electrons overflowing from the filament will be reduced with the time, which is expressed in the form of dose attenuation. In order to reduce the impact of dose attenuation on image reconstruction, the linear, quadratic and cubic fit were adopted to fit the raw data. The goodness of fit (R²) for these three fits are calculated, which is 0.9224 for the linear fit, 0.9933 for the quadratic fit and 0.9946 for the cubic fit. Relatively, three-step iterative algorithm with 2, 3, 4 dose fluctuation parameters were chosen in the experiment and the average RMSE are used for comparison, which are 0.0116, 0.0051 and 0.0051 for the hornet specimen and 0.0149, 0.0059 and 0.0058 for the Gephyrocharax melanocheir specimen. Even though the reconstruction precision with the cubic fit (4 parameters) is a little higher than the quadratic fit (3 parameters), the calculation time of the whole algorithm is increased and the minimum number of the steps required for this algorithm will also be increased due to the introduction of one more variable. At the same time, the reconstruction precision with three dose fluctuation parameters are much higher than with two. Balancing the total calculation time and the accuracy of the reconstructed images, it is rational to choose a quadratic model as the dose fluctuation model in our experiment. As for other Talbot-Lau phase contrast imaging system, the dose fluctuation model may differ. However, TSE can still be utilized to fit the dose fluctuations model and the proposed three-step iterative algorithm still works.

From Fig. 4 and Fig. 5, it is clear that the moiré artifacts caused by the phase-stepping errors and the dose fluctuations are mostly reduced. However, these artifacts cannot be completely eliminated via this algorithm in the experiment. The cause for this phenomenon is that the Eq. (2) is the approximation of the photon distributions at the detector. If higher-order harmonics for the photon distributions at the detector are considered in the expression, the moiré artifacts will be reduced further and may even be eliminated. It is also proved in the simulation that the noise is another important factor which influence the performance of the proposed algorithm. The improvement of the proposed algorithm based on the above two factors will be investigated in future work.

In this work, the influence caused by both the stepping errors and the dose fluctuations on the reconstructed images are theoretically derived in Supplement 1 and first systematically summarized. The dose fluctuation is expressed with Taylor series expansion and optimized along with the phase-stepping errors during iterations for high precision of the reconstructed images, and it is proved to be effective in moiré artifact reduction. The proposed algorithm is demonstrated to be more accurate when compared to some other algorithms. In that case, large step size motorized translation stages with low accuracy and conventional stability X-ray tube can be adopted to reduce cost for Talbot-Lau X-ray phase contrast imaging system.

5. Conclusions

To conclude, phase-stepping error and dose fluctuations are two main factors which will cause moiré artifacts in the reconstructed transmission, differential phase contrast and dark-field images in Talbot-Lau X-ray phase contrast imaging. In this paper, the errors caused by the above two factors on the reconstructed images are mathematically derived. Three-step least square iterative algorithm is proposed for the reduction of the moiré artifacts. The Taylor series expansion is utilized to fit the dose fluctuations model in the experiment. It is demonstrated that the proposed algorithm is feasible and accurate both in the simulation and experiment. This work is able to decrease the restrictions for the stability of the system. Low accuracy motorized translation stages can be utilized in the experiment. We believe that our work will lower the cost for building a Talbot-Lau X-ray phase contrast imaging system, and may be constructive for real X-ray phase contrast imaging clinical use.

Funding

Zhejiang Provincial Ten Thousand Plan for Young Top Talents (2020R52001); Zhejiang Lab (2020MC0AE01); Key Research and Development Program of Zhejiang Province (2020C01116); Major Program of the Natural Science Foundation of Zhejiang Province (LD21F050002); National Natural Science Foundation of China (61735017, 61827825, 62125504).

Acknowledgment

The authors would like to acknowledge the NIH Clinical Center for providing the Chest X-ray image.

Disclosures

The authors declare no conflicts of interest.

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.

Supplemental document

See Supplement 1 for supporting content.

References

1. E. Arendse, O. A. Fawole, L. S. Magwaza, and U. L. Opara, “Non-destructive prediction of internal and external quality attributes of fruit with thick rind: A review,” J. Food Eng. 217, 11–23 (2018). [CrossRef]

2. M. Holler, M. Guizar-Sicairos, E. H. R. Tsai, R. Dinapoli, E. Muller, O. Bunk, J. Raabe, and G. Aeppli, “High-resolution non-destructive three-dimensional imaging of integrated circuits,” Nature 543(7645), 402–406 (2017). [CrossRef]

3. K. Willer, A. A. Fingerle, W. Noichl, et al., “X-ray dark-field chest imaging for detection and quantification of emphysema in patients with chronic obstructive pulmonary disease: a diagnostic accuracy study,” Lancet Digit. Health 3(11), e733–e744 (2021). [CrossRef]

4. U. Bonse and M. Hart, “An X-Ray Interferometer,” Appl. Phys. Lett. 6(8), 155–156 (1965). [CrossRef]

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

6. C. David, B. Nohammer, H. H. Solak, and E. Ziegler, “Differential x-ray phase contrast imaging using a shearing interferometer,” Appl. Phys. Lett. 81(17), 3287–3289 (2002). [CrossRef]

7. A. Momose, S. Kawamoto, I. Koyama, Y. Hamaishi, K. Takai, and Y. Suzuki, “Demonstration of X-Ray Talbot interferometry,” Jpn. J. Appl. Phys. 42(7B), L866–L868 (2003). [CrossRef]

8. 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]

9. F. Pfeiffer, M. Bech, O. Bunk, P. Kraft, E. F. Eikenberry, C. Bronnimann, C. Grunzweig, and C. David, “Hard-X-ray dark-field imaging using a grating interferometer,” Nat. Mater. 7(2), 134–137 (2008). [CrossRef]

10. H. Yoshioka, Y. Kadono, Y. T. Kim, H. Oda, T. Maruyama, Y. Akiyama, T. Mimura, J. Tanaka, M. Niitsu, Y. Hoshino, J. Kiyohara, S. Nishino, C. Makifuchi, A. Takahashi, Y. Shinden, N. Matsusaka, K. Kido, and A. Momose, “Imaging evaluation of the cartilage in rheumatoid arthritis patients with an x-ray phase imaging apparatus based on Talbot-Lau interferometry,” Sci. Rep. 10(1), 6561 (2020). [CrossRef]

11. A. Momose, “X-ray phase imaging reaching clinical uses,” Phys. Medica 79, 93–102 (2020). [CrossRef]

12. N. Bevins, J. Zambelli, K. Li, Z. H. Qi, and G. H. Chen, “Multicontrast x-ray computed tomography imaging using Talbot-Lau interferometry without phase stepping,” Med. Phys. 39(1), 424–428 (2011). [CrossRef]

13. H. Wen, E. E. Bennett, M. A. Hegedus, and S. C. Carroll, “Spatial harmonic imaging of X-ray scattering initial results,” IEEE Trans. Med. Imaging 27(8), 997–1002 (2008). [CrossRef]

14. M. Seifert, V. Ludwig, M. Gallersdorfer, C. Hauke, K. Hellbach, F. Horn, G. Pelzer, M. Radicke, J. Rieger, S. M. Sutter, T. Michel, and G. Anton, “Single-shot Talbot-Lau x-ray dark-field imaging of a porcine lung applying the moire imaging approach,” Phys. Med. Biol. 63(18), 185010 (2018). [CrossRef]

15. L. Romano and M. Stampanoni, “Microfabrication of X-ray Optics by Metal Assisted Chemical Etching: A Review,” Micromachines 11(6), 589 (2020). [CrossRef]

16. C. Hauke, M. Leghissa, G. Pelzer, M. Radicke, T. Weber, T. Mertelmeier, G. Anton, and L. Ritschl, “Analytical and simulative investigations of moire artefacts in Talbot-Lau X-ray imaging,” Opt. Express 25(26), 32897–32909 (2017). [CrossRef]

17. K. Hashimoto, H. Takano, and A. Momose, “Improved reconstruction method for phase stepping data with stepping errors and dose fluctuations,” Opt. Express 28(11), 16363–16384 (2020). [CrossRef]

18. K. Okada, A. Sato, and J. Tsujiuchi, “Simultaneous calculation of phase distribution and scanning phase shift in phase shifting interferometry,” Opt. Commun. 84(3-4), 118–124 (1991). [CrossRef]

19. Z. Y. Wang and B. T. Han, “Advanced iterative algorithm for phase extraction of randomly phase-shifted interferograms,” Opt. Lett. 29(14), 1671–1673 (2004). [CrossRef]

20. J. Vargas, C. O. S. Sorzano, J. C. Estrada, and J. M. Carazo, “Generalization of the Principal Component Analysis algorithm for interferometry,” Opt. Commun. 286, 130–134 (2013). [CrossRef]

21. S. Kaeppler, J. Rieger, G. Pelzer, F. Horn, T. Michel, A. Maier, G. Anton, and C. Riess, “Improved reconstruction of phase-stepping data for Talbot-Lau x-ray imaging,” J. Med. Imaging 4(03), 1 (2017). [CrossRef]

22. C. Hauke, G. Anton, K. Hellbach, M. Leghissa, F. G. Meinel, T. Mertelmeier, T. Michel, M. Radicke, S. M. Sutter, T. Weber, and L. Ritschl, “Enhanced reconstruction algorithm for moire artifact suppression in Talbot-Lau x-ray imaging,” Phys. Med. Biol. 63(13), 135018 (2018). [CrossRef]

23. J. W. Chen, J. T. Zhu, Z. C. Li, W. Shi, Q. Y. Zhang, Z. L. Hu, H. R. Zheng, D. Liang, and Y. S. Ge, “Automatic image-domain Moire artifact reduction method in grating-based x-ray interferometry imaging,” Phys. Med. Biol. 64(19), 195013 (2019). [CrossRef]

24. S. W. Tao, C. X. He, X. Hao, C. F. Kuang, and X. Liu, “Factors Affecting the Spatial Resolution in 2D Grating-Based X-Ray Phase Contrast Imaging,” Front. Phys. 9, 672207 (2021). [CrossRef]

25. X. S. Wang, Y. F. Peng, L. Lu, Z. Y. Lu, M. Bagheri, and R. M. Summers, “ChestX-ray8: Hospital-scale Chest X-ray Database and Benchmarks on Weakly-Supervised Classification and Localization of Common Thorax Diseases,” in 30th Conference on Computer Vision and Pattern Recognition (CVPR) (IEEE, 2017), pp. 3462–3471.

	With the tested object			The background
Parameters	Preset	Calculated	RE	Preset	Calculated	RE
Step 1 (µm)	0.7633	0.7691	0.7673%	1.5317	1.5326	0.1181%
Step 2 (µm)	1.8996	1.8886	-0.5785%	1.1021	1.1017	-0.0900%
Step 3 (µm)	0.4702	0.4730	0.5789%	1.0685	1.0693	0.0914%
Step 4 (µm)	1.8622	1.8580	-0.2222%	0.6780	0.6785	0.1363%
a₀	1.0000	1.0025	0.2458%	1.0000	1.0016	0.1627%
a₁	-0.1000	-0.1013	-1.2762%	-0.1000	-0.1008	-0.8113%
a₂	0.0060	0.0060	-0.1605%	0.0060	0.0060	-0.0986%

Moiré artifacts reduction in Talbot-Lau X-ray phase contrast imaging using a three-step iterative approach

Abstract

1. Introduction

2. Methods

2.1 Phase distributions calculation

2.2 Phase-stepping errors calculation

2.3 Dose fluctuation coefficients calculation

2.4 Iteration

3. Simulation conditions and results

4. Experiment results and discussions

5. Conclusions

Funding

Acknowledgment

Disclosures

Data availability

Supplemental document

References

Supplementary Material (1)

Data availability

Cited By

Figures (6)

Tables (1)

Equations (30)

Optics Express