Dual-energy phase retrieval algorithm for inline phase sensitive x-ray imaging system

Muhammad U. Ghani; Aimin Yan; Laurie. L. Fajardo; Xizeng Wu; Hong Liu

doi:10.1364/OE.431623

1. Introduction

Conventional X-ray imaging for soft tissues (such as breast, brain, liver, etc.) is limited in its sensitivity for detecting subtle tissue pathological changes, since the imaging relies on small differences in X-ray attenuation between the lesions and soft tissues of variable structure. In contrast to conventional X-ray imaging, X-ray phase contrast imaging, which has been of interest in the last couple of decades for potential applications in medical imaging and material science, relies on the X-ray phase shifts that tissue can generate. The amount of X-ray phase shift is determined by the summation of tissue electron density over the ray path, i.e.

\phi(\vec{\mathbf{r}})={-}\lambda {r_{\mathrm{e}}}\int \rho_{\mathrm{e}}(\vec{\mathbf{r}}) ds={-}\lambda {r_{\mathrm{e}}} \rho_{\mathrm{ep}},

where $\lambda$ is the X-ray wavelength, ${r_{\mathrm {e}}}$ the classical electron radius and $\rho _{\mathrm {ep}}$ [1–3], the integral of tissue electron densities, $\rho _{\mathrm {e}}$, over the ray path, is frequently called the projected electron density. In the X-ray energy range of 5–200 keV the differences in X-ray phase shifts between soft tissues are about 2-3 orders of magnitude greater than the differences in the projected linear attenuation coefficients [1–3]. Therefore, X-ray phase-contrast imaging has great potential to increase detection sensitivity. Different from the interferometry-based setups that use elaborate X-ray optics (such as gratings or an analyzer crystal) for phase measurements [4,5], inline propagation method utilizes interference patterns at longer sample-to-detector distances with an X-ray source of sufficiently large lateral coherence [2,6]. In the inline phase contrast imaging, the X-rays with shifted phases diffract from the object to the detector, and the diffracted X-rays form bright and dark fringes at tissue boundaries. The edge enhancement thus generated relies on the spatial coherence of the X-ray source, the Laplacian and gradients of X-ray phase shifts caused by the object’s, and the gradients of the objects attenuation [7–17].

On a recorded phase contrast image, tissue attenuations and phase shifts are entangled, which are both unknown in practice. The procedure of disentangling tissue phase shifts from the mixed contrast mechanism and retrieving the phase maps from acquired phase contrast images is called phase retrieval. Phase retrieval technique plays a central role in phase-sensitive X-ray imaging. In general, at least two images with varying object-detector distances or at different X-ray energies per projection are required to correctly reconstruct the phase shifts and attenuations [18–21]. This adds complexity to the experimental setup and increases the radiation dose delivered to the tissue. In order to perform accurate reconstruction with one image, additional validity assumptions must be added. Different from the other methods [9,22–24], in which a priori knowledge about the imaged sample is assumed, Wu et al. [16] developed a Phase-Attenuation Duality (PAD) algorithm for soft tissues of low atomic numbers (Z<10) with high X-ray energies ($E\geq 60$ keV). This algorithm is valid because for soft tissues which are composed mainly of low atomic number elements (Z<10), a condition, called the phase-attenuation duality (PAD), is satisfied when X-ray energy is $60$ keV or higher. In other words, at this higher X-ray energy range, the effects related to the photoelectric effect (PE) are negligible. This PAD algorithm is robust add has the capability of (high frequency) noise reduction [16]. The PAD algorithm also has its deficiencies. For the high level of X-ray energy, the absorption contrast is lost, which in principle yields variations in atomic number composition (the current gold standard in radiology) [25].

Current scientific advancements have enabled the application of Photon-Counting Detectors (PCDs) in medical imaging research [26,27]. Unlike the energy integrating detectors, PCDs have advanced application specific integrated circuits (ASICs) designed for readout which allows photons energy classification while suppressing electronic noise for improved image quality [28]. By selecting appropriate energy thresholds, the incident photons are registered and counted into different bins based on their energies. Hence, with a PCD, one can get multiple low dose energy bin images from a single X-ray exposure. This unique ability of the PCDs can be potentially utilized for the phase retrieval and accurate quantification of tissues with multiple images without increasing the x-ray dose.

Furthermore, this photon counting technique also eliminates the chance of misalignment frequently encountered with multiple exposures. Gureyev et al. [29,30] presented algorithms that use a material-independent scaling of attenuation and phase effects with respect to the X-ray energy. This model works only for the low X-ray energies (1–20 keV) where the PE effects dominate [29,30]. Recently, Gürsoy & Das [25] and Vazquez et al. [31] employed the Alvarez-Macovski (AM) model [32] to separate attenuation as a combination of PE absorption and Compton scattering and derived a linear equation in the frequency domain. The problem is then to find the solution with some proper linear equation solver, such as the least square approximation. This one step solution in Fourier space is efficient and accurate if the input data has limited quantum noise. Schaff et al. [33] have shown that spectral phase retrieval using the Alvarez-Macovski basis is equivalent to using well-behaved filters in Fourier-space. The robustness of a linear equation depends on its condition number, i.e. the ratio of the largest eigenvalue to the smallest eigenvalue. In medical imaging system, the data size is considerably large. For example, with a setup of source-to-object and object-to-detector distance $R_1=R_2=1$m, and X-ray energy $E=20$ keV, when the detector size is $2048\times 2048$ with pixel pitch = $50\mu$m, the condition number for matrix $\sigma _{\mathrm {kn}}+2\pi (\lambda ^{2} {r_{\mathrm {e}}}R_2)/M\cdot \vec {\mathbf {u}}^{2}$ can be as large as $10^{36}$. Special care must be taken if high frequency noise are present. Li et al. [34,35] employed the same PE absorption and Compton scattering partition with the Alvarez-Macovski (AM) model [32] but with an iterative procedure using finite difference method to estimate the residue to improve the solution. The complexity of the algorithm makes it hard to be widely accepted.

In this paper, we adopt similar PE absorption and Compton scattering partition of the attenuation. With this partition, we separate the phase contrast image into two contributions $I={A_{\mathrm {pe}}}^{2}\times {I_{\mathrm {kn}}}$, where ${A_{\mathrm {pe}}}^{2}$ is the attenuation generated by PE absorption, and ${I_{\mathrm {kn}}}$ is a “pseudo” phase contrast image caused by electron density. Since ${I_{\mathrm {kn}}}$ is determined only by electron density, the phase-attenuation duality condition is satisfied, so the PAD algorithm is valid for phase retrieval from ${I_{\mathrm {kn}}}$ [16]. Furthermore, the known relationship of linear attenuation coefficients with respect to X-ray energy helps the transformation of attenuation caused by PE absorption between different X-ray energies. This way, one can develop a simple iterative algorithm. A flow chart of the algorithm is shown in Fig. 1. Briefly, with given ${A_{\mathrm {pe}}}^{2}$ at energy $E_1$, one can solve for the energy independent projected electron density $\rho _{\mathrm {ep}}$ through PAD algorithm. Then the wavefront ${A_{\mathrm {kn}}}\cdot \exp \left [i\phi \right ]$ at energy $E_2$ is acquired with the Compton scattering effect. By applying Fresnel propagation, one gets the pseudo phase contrast image ${I_{\mathrm {kn}}}$, and thus the attenuation arising from PE absorption ${A_{\mathrm {pe}}}^{2}=I_2/{I_{\mathrm {kn}}}$ at energy $E_2$. Finally, using the relationship of X-ray energy with respect to linear attenuation coefficients caused by PE absorption, one gets ${A_{\mathrm {pe}}}$ at energy $E_1$. The well-conditioned PAD algorithm and Fresnel propagation guarantees the convergence and stability of the iterations.

Fig. 1. Flow chart of the dual-energy iterative algorithm. $I_1$, and $I_2$ are input phase contrast images with corresponding X-ray energies $E_1$, and $E_2$ ($E_1>E_2$) respectively. Here $\mathfrak {D}$ and $\mathfrak {F_r}$ represents the PAD and Fresnel propagation methods respectively.

Download Full Size | PDF

The paper is organized as follows. In section 2., we first separate the phase contrast image into two contributions, the attenuation caused by PE absorption and the pseudo phase contrast image caused by electron density. We then develop the iterative algorithm by adopting the PAD algorithm and the forward Fresnel propagation method in section 3. In sections 4. and 5., we perform experimental validation of the proposed algorithm with the inline phase contrast images acquired by a photon counting detector. We present our conclusion in section 6.

2. Method

Phase retrieval for the inline phase contrast imaging is based on the X-ray propagation equations derived from the Fresnel diffraction theory or the Wigner distribution based on the parabolic wave equation for optical phase-space formalism [6,7,11,36,37]. To be specific, let $\phi (\vec {\mathbf {r}})$ be the X-ray phase shift caused by the imaging object, and $A(\vec {\mathbf {r}})$ be the X-ray transmission or the attenuation-map of the object. Then the Fourier transform of the X-ray intensity,

\hat{\mathcal{F}}\hspace{-2.pt}\left\lgroup {I}\right\rgroup\hspace{-2.pt} (\vec{\mathbf{u}})=\int_{\mathbb{R}^{2}} I(\vec{\mathbf{r}})\exp\left[2\pi i \vec{\mathbf{r}}\cdot\vec{\mathbf{u}}\right]\,\mathrm{d}\vec{\mathbf{r}},

at the detector entrance away from the object with distance $R_2$, of the monochromatic point X-ray source at a place away from the object $R_1$ distance, can be modeled by

(1)$$\begin{aligned}&\hat{\mathcal{F}}\hspace{-2.pt}\left\lgroup {I}\right\rgroup\hspace{-2.pt} \left(\frac{\vec{\mathbf{u}}}{M};R_1+R_2\right) = I_{\mathrm{in}}\left\{\cos\left(\frac{\pi\lambda R_2}{M}\vec{\mathbf{u}}^{2}\right)\cdot\hat{\mathcal{F}}\hspace{-2.pt}\left\lgroup {A^{2}}\right\rgroup\hspace{-2.pt} {\vec{\mathbf{u}}} + \right.\\ &\qquad+\left[2\sin\left(\frac{\pi\lambda R_2}{M}\vec{\mathbf{u}}^{2}\right)-2\left(\frac{\pi\lambda R_2}{M}\vec{\mathbf{u}}^{2}\right)\cdot\cos\left(\frac{\pi\lambda R_2}{M}\vec{\mathbf{u}}^{2}\right)\right]\cdot\hat{\mathcal{F}}\hspace{-2.pt}\left\lgroup {A^{2}\phi}\right\rgroup\hspace{-2.pt} - \\ &\qquad\qquad- \cos\left(\frac{\pi\lambda R_2}{M}\vec{\mathbf{u}}^{2}\right)\cdot\frac{\lambda R_2}{2\pi M}\cdot\hat{\mathcal{F}}\hspace{-2.pt}\left\lgroup {\nabla(A^{2}\nabla\phi)}\right\rgroup\hspace{-2.pt} - \\ &\qquad\qquad\qquad\left.-\frac{\lambda R_2}{4\pi M} \sin\left(\frac{\pi\lambda R_2}{M}\vec{\mathbf{u}}^{2}\right)\cdot\hat{\mathcal{F}}\hspace{-2.pt}\left\lgroup {\nabla^{2}A^{2}}\right\rgroup\hspace{-2.pt} \right\},\end{aligned}$$

where ${I_{\mathrm {in}}}$ is the incident X-ray intensity at $R_1$, $\lambda$ is the wavelength of the monochromatic point X-ray source and $M = (R_1 + R_2 )/R_1$ is the geometric magnification. When the FT-Space Fresnel propagator $\pi \lambda R_2 \vec {\mathbf {u}}^{2}/M\ll 1$, Eq. (1) can be simplified to the Transport of Intensity Equation (TIE) [11,38,39]

(2)$$I(\vec{\mathbf{r}}; R_1+R_2; E) = \frac{I_{\mathrm{in}}}{M^{2}}\cdot\left\{A^{2}\left(\frac{\vec{\mathbf{r}}}{M}; E\right) - \frac{\lambda(E) R_2}{2\pi M}\nabla \left[A^{2}\left(\frac{\vec{\mathbf{r}}}{M}; E\right)\cdot\nabla \phi\left(\frac{\vec{\mathbf{r}}}{M}; E\right)\right]\right\}.$$

The requirement of $\pi \lambda R_2 \vec {\mathbf {u}}^{2}/M\ll 1$ is necessary since for high-resolution images the FT-Space Fresnel propagator $\pi \lambda R_2 \vec {\mathbf {u}}^{2}/M$ is close or greater than $\pi$, any phase retrieval algorithm based on Eq. (2) needs to be substituted to Eq. (1) [20,40,41]. In this paper, the algorithm is also based on Eq. (2), i.e. the moderate resolution requirements are satisfied.

It is well known that when the X-rays pass through a tissue, the X-ray phase shift arises from the coherent X-ray scattering by tissues. That is $\phi (\vec {\mathbf {r}})=-\lambda {r_{\mathrm {e}}} \rho _{\mathrm {ep}}(\vec {\mathbf {r}})$, where $\rho _{\mathrm {ep}}(\vec {\mathbf {r}})=\int \rho _{\mathrm {e}}(\vec {\mathbf {r}},z)\,\mathrm {d} z$ is the integral of tissue electron density along the ray path, and is called the projected electron density. For the diagnostic X-rays, the tissue attenuation arises from three X-ray tissue interactions: the PE absorption, coherent scattering and incoherent scattering (also known as Compton scattering) [11,16,42,43]. In other words, the X-ray attenuation is given as $A^{2}(\vec {\mathbf {r}}) = \exp \left [-\int \mu (\vec {\mathbf {r}},z)\,\mathrm {d} z\right ]$, where the linear attenuation coefficient is written as

\mu(\vec{\mathbf{r}},z) = \mu_{\mathrm{pe}}(\vec{\mathbf{r}},z) + \mu_{\mathrm{coh}}(\vec{\mathbf{r}},z) + \mu_{\mathrm{incoh}}(\vec{\mathbf{r}}, z).

Here $\mu _{\mathrm {pe}}(\vec {\mathbf {r}},z)$, $\mu _{\mathrm {coh}}(\vec {\mathbf {r}},z)$ and $\mu _{\mathrm {incoh}}(\vec {\mathbf {r}}, z)$ are the linear attenuation coefficients generated by PE absorption, coherent scattering and incoherent scattering respectively. Of the three contributions, $\mu _{\mathrm {pe}}$ has the most contribution when low energy X-rays (1–20 keV) are employed for imaging [29,30]. On the other hand, when high-energy X-rays (60-500 keV) are used, the effects related to PE and coherent scattering are minimum for light elements with atomic numbers $Z<10$, thus leaving $\mu _{\mathrm {incoh}}$ as the primary contributor. Alvarez-Macovski (AM) attenuation model describes the attenuation as a combination of PE absorption $\mu _{\mathrm {pe}}$, and Compton scattering $\mu _{\mathrm {incoh}}$ for the intermediate energy range [32]. Our phase retrieval algorithm adopts the AM model to describe the tissue attenuation as

(3)$$A^{2}\left(\vec{\mathbf{r}}\right) = \exp\left[-\int \mu(\vec{\mathbf{r}},z)\,\mathrm{d} z\right] = \exp\left[-\int \mu_{\mathrm{pe}}(\vec{\mathbf{r}},z)\,\mathrm{d} z\right]\cdot \exp\left[-\int \mu_{\mathrm{incoh}}(\vec{\mathbf{r}},z)\,\mathrm{d} z\right].$$

Note $\mu _{\mathrm {pe}}$ and $\mu _{\mathrm {incoh}}$ relate X-ray energy in the form [32]

(4)$$\left\{\begin{array}{l} \mu_{\mathrm{pe}}(\vec{\mathbf{r}},z) \propto \frac{1}{E^{3}}, \\ \mu_{\mathrm{incoh}}(\vec{\mathbf{r}},z) = \sigma_{\mathrm{kn}}(E) \rho_{\mathrm{e}}(\vec{\mathbf{r}},z), \end{array}\right.$$

where $\sigma _{\mathrm {kn}}$ is the total cross-section for X-ray photon Compton scattering from a single free electron, which can be determined by the time-honored Klein-Nishina formula:

(5)$$\sigma_{\mathrm{kn}}(E) = 2\pi r_e^{2}\left\{\frac{1+\eta}{\eta^{2}}\left[\frac{2(1+\eta)}{1+2\eta} - \frac{1}{\eta} \log(1+2\eta)\right] + \frac{1}{2\eta} \log(1+2\eta) - \frac{1+3\eta}{(1+2\eta)^{2}}\right\},$$

where $\eta = E/m_ec^{2}$, with $E$ denotes the photon energy of the primary X-ray beam and $m_ec^{2}$ the resting electron energy. Let $\mu _0(\vec {\mathbf {r}})$ be the ray-integral of the linear attenuation coefficient generated by PE absorption at energy $E_0=1$ keV. One can split the attenuation into two contributions

(6)$$A^{2}(\vec{\mathbf{r}};E) = {A_{\mathrm{pe}}}^{2}(\vec{\mathbf{r}}; E)\cdot {A_{\mathrm{kn}}}^{2}(\vec{\mathbf{r}}; E),$$

where

(7)$$\left\{\begin{array}{l} {A_{\mathrm{pe}}}^{2}(\vec{\mathbf{r}}; E) = \exp\left[-\mu_0\cdot \left(E_0/E\right)^{3} \right],\\ {A_{\mathrm{kn}}}^{2}(\vec{\mathbf{r}}; E) = \exp\left[-\sigma_{\mathrm{kn}}(E)\rho_{\mathrm{ep}}(\vec{\mathbf{r}})\right]. \end{array}\right.$$

Assuming the fluctuation of ${A_{\mathrm {pe}}}$ is negligible compared to the phase change and by substituting Eq. (6) into Eq. (2), one arrives

(8)$$I\left(\vec{\mathbf{r}}; R_1+R_2; E\right) = {A_{\mathrm{pe}}}^{2}\left(\frac{\vec{\mathbf{r}}}{M};E\right)\cdot {I_{\mathrm{kn}}}\left(\vec{\mathbf{r}};R_1+R_2;E\right),$$

where

(9)$$\begin{aligned} {I_{\mathrm{kn}}}\left(\vec{\mathbf{r}};R_1+R_2;E\right) =& \frac{{I_{\mathrm{in}}}}{M^{2}}\cdot \left\{{A_{\mathrm{kn}}}^{2}\left(\frac{\vec{\mathbf{r}}}{M};E\right) - \frac{\lambda(E) R_2}{2\pi M}\nabla \left[{A_{\mathrm{kn}}}^{2}\left(\frac{\vec{\mathbf{r}}}{M};E\right)\cdot\nabla \phi\left(\frac{\vec{\mathbf{r}}}{M};E\right)\right]\right\}\\ =& \frac{{I_{\mathrm{in}}}}{M^{2}}\cdot \left\{{A_{\mathrm{kn}}}^{2}\left(\frac{\vec{\mathbf{r}}}{M};E\right) + \frac{\lambda(E) R_2}{2\pi M} \cdot \frac{\lambda(E) {r_{\mathrm{e}}}}{\sigma_{\mathrm{kn}}(E)} \cdot \nabla^{2} {A_{\mathrm{kn}}}^{2}\left(\frac{\vec{\mathbf{r}}}{M};E\right)\right\}.\end{aligned}$$

Many researchers chose to simplify the TIE equation by factoring $A^{2}(\vec {\mathbf {r}})$ out of the differential operator as an approximation for the TIE [25,31,44]. Here we keep ${A_{\mathrm {kn}}}(\vec {\mathbf {r}})$ untouched as this approach entails a better approximation and the above PDE (Eq. (9)) is well-conditioned. Note that Eq. (9) satisfies the phase-attenuation duality (PAD) condition [16]. Thus, the PAD algorithm can be used to retrieve ${A_{\mathrm {kn}}}$ and thus the projected electron density $\rho _{\mathrm {ep}}$ using Eq. (7) provided that ${I_{\mathrm {kn}}}$ is known. One can see from the Eq. (9) that ${I_{\mathrm {kn}}}$ is the phase contrast image arising from Compton scattering. Thus, ${I_{\mathrm {kn}}}$ can be computed once ${A_{\mathrm {kn}}}$ is known. Li et al. simulated the forward propagation using the finite difference method [34]. This method is complicated and most importantly, it usually has difficulty when dealing with noisy data. We choose a different approach as follow to simulate the forward propagation. Recall the Fresnel wave propagation equation from $R_1$ to $R_2$, for a point source, can be expressed as

(10)$$W(\xi; R_1+R_2) = \frac{\sqrt{I_{\mathrm{in}}}}{\lambda R_{2}}\exp\left[{-}2\pi i\frac{R_{1}+R_{2}}{\lambda}\right]\cdot\exp\left[i\pi\frac{\xi^{\,2}}{\lambda(R_{1}+R_{2})}\right]\cdot \left(G\circledast T\right)(\xi),$$

where $I_{\mathrm {in}}$ denotes the incident X-ray intensity at $R_1$, $\lambda$ is the X-ray wave length,

G(x)=\exp\left[i \pi x^{2}\cdot M/(\lambda R_2)\right],

$M=(R_1+R_2)/R_1$, and $T(x)$ is the wavefront at $R_1$ plane. The symbol $\circledast$ in the above equation represents the convolution operator. The wavefront $W$ can be numerically simulated through Fast Fourier Transform (FFT). Since $G(x)$ is the Fresnel propagator, which is a Gaussian-like function, the convolution operation indicates it plays a role of spreading and smoothing the wavefront in the evolution process. With this method, the noise can be effectively restrained.

Keeping these considerations in mind, we developed an iterative algorithm for dual-energy phase retrieval that is explained in the next section.

3. Dual-energy phase reconstruction algorithm

Assume that we have two phase contrast images $I_{1}$, and $I_{2}$ which are acquired with the same geometric setup and a perfect detector. Assume that the X-ray energies for $I_{1}$ and $I_{2}$ images are $E_1$, and $E_2$ with $E_1>E_2$. The algorithm is as follow.

Step 1: First assume that the attenuation arising from PE absorption at energy $E_1$ is ${A_{\mathrm {pe}}}^{2}(E_1)=1$.
Step 2: From Eq. (8) we know that the phase contrast image arises from electron density is given by ${I_{\mathrm{kn}}}(E_1)=I_1(E_1)/{A_{\mathrm{pe}}}^{2}(E_1).$ Using ${I_{\mathrm {kn}}}$, compute the projected electron density $\rho _{\mathrm {ep}}$ by applying the already developed PAD based method [16].
Step 3: Compute the X-ray wave amplitude attenuated by Compton scattering ${A_{\mathrm {kn}}}$, and phase shift $\phi$ at energy $E_2$, based on the relationship of ${A_{\mathrm {kn}}}$ and $\phi$ with respect to $\rho _{\mathrm {ep}}$. Then Fresnel wave propagate the sample’s transmittance $T(\vec{\mathbf{r}}; E_2)={A_{\mathrm{kn}}}(\vec{\mathbf{r}}; E_2)\cdot\exp\left[i\phi(\vec{\mathbf{r}}; E_2)\right]$ from object plane $R_1$ to detector entrance at $R_2$, one gets the phase contrast image ${I_{\mathrm {kn}}}(E_2)$.
Step 4: Using Eq. (8), the X-ray wave amplitude attenuated by PE absorption is computed as ${A_{\mathrm{pe}}}^{2}(E_2)=I_2(E_2)/{I_{\mathrm{kn}}}(E_2).$ Then the ray-integral of the linear attenuation coefficient generated by PE absorption at energy $E_0=1$ keV is computed from $\mu_0={-}\left(E_2/E_0\right)^{3}\cdot\log({A_{\mathrm{pe}}}^{2}(E_2)).$
Step 5: If the ($L_2$-norm) error of $\mu _0$ between two adjacent iteration does not decrease, the iteration stops. Otherwise, compute ${A_{\mathrm {pe}}}^{2}$ at energy $E_1$ with ${A_{\mathrm{pe}}}^{2}(E_1)=\exp\left[-\mu_0\cdot\left(E_0/E_1\right)^{3}\right],$ and go back to Step 2.

Remarks:

1. The algorithm can be readily extended to multiple phase contrast images. To replace the two input phase contrast images $I_1$ and $I_2$ with two-phase contrast image groups $I_{g1}$ and $I_{g2}$, the energies in group $I_{g1}$ must be greater than the energies in group $I_{g2}$. The corresponding $\mu _0$ and $\rho _{\mathrm {ep}}$ are taken as the average of each individual $\mu _0$ and $\rho _{\mathrm {ep}}$. This way, the noise can be further constrained, and a better quality of the retrieved images is expected.
2. The requirement of $E_1>E_2$ is necessary. This is because the conversion from ${A_{\mathrm {pe}}}(E_2)$ to ${A_{\mathrm {pe}}}(E_1)$: ${A_{\mathrm{pe}}}^{2}(E_1)=\exp\left[(E_2/E_1)^{3}\log({A_{\mathrm{pe}}}^{2}(E_2)\right],$ can easily diverge if $E_1<E_2$.
3. Once $\mu _0$ and $\rho _{\mathrm {ep}}$ are retrieved, one can get the tissue’s attenuation change $A$ and phase shift $\phi$ at any energy: $(11)$$\left\{\begin{array}{lll} A^{2}(E) & = &\exp\left[-\mu_0/E^{3} - \sigma_{\mathrm{kn}}(E)\cdot \rho_{\mathrm{ep}}\right],\\ \phi(E) & = &{-}\lambda(E) r_{\mathrm{e}} \cdot \rho_{\mathrm{ep}}. \end{array}\right.$$$

4. Experimental validation

To validate the proposed phase retrieval algorithm, a multi-material phantom was imaged with an inline phase sensitive x-ray imaging prototype employing a photon counting detector (PCD). The prototype consists of a microfocus x-ray tube (Hamamatsu Photonics, Japan) that has a tungsten (W) anode and 0.5 mm thick Beryllium (Be) output window with tube voltage and tube current ranging from 40–130 kV and 10–300 $\mu$A, respectively. The source focal spot size varies from 16–50 $\mu$m as the output power varies from 8–39 W. The x-ray emission beam angle is approximately $100^{\circ }$, and a collimator is employed to limit the x-ray beam size and scattering. A 3.5 mm thick aluminum (Al) filter was placed in front of the output window for beam hardening. The prototype incorporates a high-resolution PCD (Advacam s.r.o, Czech Republic) that has a 1-mm thick cadmium telluride (CdTe) semiconducting layer for the direct detection under a bias voltage of 300 V. The PCD has five Medipix3RX chips with a pixel array of $256 \times 1280$ and a pixel pitch of 55 $\mu$m. Each pixel has two energy discrimination thresholds and two integrated 12-bit digital counters. Both counters can be joined to a single 24-bit counter, providing an improved dynamic range. The PCD can operate in imaging mode and spectroscopic mode. In imaging mode, the threshold values in energy (keV) are fixed values when detecting X-ray photons. In spectroscopic mode, the detector sweeps the whole energy range from low energy threshold to high energy threshold through sequentially increasing the threshold by a preset threshold step. The source to object distance ($R_1$) was 68.58 cm, while the source to detector distance ($R_1+R_2$) was 137.16 cm. The imaging setup used in this study fulfills the partial coherence requirements dictated by Wu et al. [45], which are necessary for the visibility of the phase contrast effects. The schematic of the prototype is given in Fig. 2(a). The imaging phantom consists of five cylindrical rods simulating the X-ray interaction of with the adipose, 25 mg/cc hydroxyapatite, muscle, lung and breast tissues. Each rod had a diameter of 5 mm and a thickness of approximately 2.25 mm. All the cylindrical rods were fixated on a 7.75 mm thick acrylic substrate as shown in the schematics of Fig. 2(b). All of the materials $Z_{\mathrm {eff}}$ values are well below 10. Spectroscopic scans of the phantom for each energy threshold were collected at 100 kV, 300 $\mu$A and 25 sec.

Spectral data was generated by threshold scanning followed by energy binning. First, a series of phase contrast scans were acquired, each with a 12 keV energy threshold (TH) window. Each TH scan comprises the photon counts with energies above the TH value. The next step involves the difference of the two adjacent scans and then divided by the difference of the two corresponding adjacent flat field images to get the resultant pseudo monochromatic phase contrast images bounded between 16-28, 29-40, 41-52, 53-64, 65-76, 77-88, 89-100 keV bins. The corresponding (weighted averaging) effective energies of these bins were found to be 24.67, 34.22, 45.74, 57.92, 69.26, 81.38, 92.34 respectively, based on the measured x-ray spectrum at 100 kV and 3.5 mm Al filter as shown in Fig. 2(c).

Fig. 2. (a) Schematics of the inline phase sensitive x-ray imaging prototype with core components and geometric conditions; (b) Schematics of the phantom consisting of 5 cylindrical rods simulating the adipose, 25 mg/cc hydroxyapatite, muscle, lung and breast tissues; (c) X-ray spectrum for a tungsten (W) target source acquired at 100 kV, 3.5 mm Al filter showing various energy bins used in this study.

Download Full Size | PDF

5. Results

Two tests were performed to evaluate the performance of the iterative dual-energy phase retrieval algorithm. In the first test, five input phase contrast images shown in Fig. 3 with energy bounded between 16-28, 29-40, 41-52, 53-64, and 65-100 keV were chosen. The five images were set into two groups: $I_{g1}=\{I_{41-52}, I_{53-64}, I_{65-100}\}$, and $I_{g2}=\{I_{16-28}, I_{29-40}\}$. The effective energies corresponding to these five images are 45.74, 57.92, 75.68, 24.67, and 34.22 keV respectively. Figure 4 shows the reconstructed projected electron density $\rho _{\mathrm {ep}}$ and the ray-integral of the linear attenuation coefficient generated by PE absorption, $\mu _0$, at energy $E_0=1$ keV after 72 iteration steps. Comparing the relative errors with the true expected value shown in Tab. 1, one concludes the results are reasonably accurate. For example, for the background acrylic, the error is as low as $0.973\%$. Even for the largest value (the muscle rod), the error is below $7.3\%$. This indicates the algorithm is accurate and robust.

Fig. 3. Input phase contrast images for the test. From top to bottom, the effective energies are 24.67, 34.22, 45.74, 57.92, 75.68 keV respectively. The input phase contrast images are grouped to $I_{g1}=\{I_{\mathrm {41-52}}, I_{\mathrm {53-64}}, I_{\mathrm {65-100}}\}$, $I_{g2}=\{I_{\mathrm {16-28}}, I_{\mathrm {29-40}}\}$. The profiles (11 pixel averaged vertically) along the central lines are shown to the right plots.

Download Full Size | PDF

Fig. 4. The retrieved projected electron density $\rho _{\mathrm {ep}}$ (top) and the ray-integral of the linear attenuation coefficient $\mu _0$ generated by PE absorption at energy $E_0=1$ keV. The input phase contrast images are $I_{g1}=\{I_{\mathrm {41-52}}, I_{\mathrm {53-64}}, I_{\mathrm {65-100}}\}$, $I_{g2}=\{I_{\mathrm {16-28}}, I_{\mathrm {29-40}}\}$. The profiles (11 pixel averaged vertically) along the central lines are shown to the right plots.

Download Full Size | PDF

Table 1. $\rho _{\mathrm {ep}}$ accuracy: $I_{g1}=\{I_{\mathrm {40-52}}, I_{\mathrm {52-64}}, I_{\mathrm {64-100}}\}$, $I_{g2}=\{I_{\mathrm {16-28}}, I_{\mathrm {28-40}}\}$.

View Table | View all tables in this article

To demonstrate the noise reduction capability of the algorithm, we compute the attenuation maps with Eq. (11) at the five effective energies 24.67, 34.22, 45.74, 57.92, and 75.68 keV. The computed images are shown in Fig. 5. The quality improvement is intuitive by comparing Fig. 5 and Fig. 3. The quantitative improvements of the contrast noise ratio (CNR) for the five synthetic tissue rods are shown in Tab. 2. The improvements are more for higher x-ray energies compared to lower ones. This phenomenon is sample-dependent. Since the sample chosen here is composed of light elements, the attenuation generated by PE absorption is relatively low compared to the attenuation generated by Compton scattering. At higher X-ray energies, the contributions of Compton scattering are greater in addition to the high photon counts.

Fig. 5. Display of the computed attenuation maps. The effective energies (from top to bottom) are 24.67, 34.22, 45.74, 57.92, and 75.68 keV respectively. The five input phase contrast images for the test are $I_{g1}=\{I_{\mathrm {40-52}}, I_{\mathrm {52-64}}, I_{\mathrm {64-100}}\}$, $I_{g2}=\{I_{\mathrm {16-28}}, I_{\mathrm {28-40}}\}$. The profiles (11 pixel averaged vertically) along the central lines are shown to the right plots.

Download Full Size | PDF

Table 2. Contrast Noise Ratio (CNR) of the input phase contrast images and the reconstructed attenuation images. $I_{g1}=\{I_{\mathrm {41-52}}, I_{\mathrm {53-64}}, I_{\mathrm {65-100}}\}$, $I_{g1}=\{I_{\mathrm {16-28}}, I_{\mathrm {29-40}}\}$.

View Table | View all tables in this article

In the second test, we combine images $I_{\mathrm {41-52}}$, $I_{\mathrm {53-64}}$, and $I_{\mathrm {65-100}}$ together to form one image $I_{\mathrm {41-100}}$, and $I_{\mathrm {16-28}}$, $I_{\mathrm {29-40}}$ together to form image $I_{\mathrm {16-40}}$. With these two images as inputs, the corresponding effective energies are $59.31$ and $31.84$ keV, respectively. This time the iteration stops after 132 steps. Tables 3 and 4 lists the error in accuracy and improvement with CNR. The results are accurate, although relatively less accurate compared to the five-image input results. In the tests above, we let the algorithm to iterate as many steps as it can, and just to demonstrate the algorithm is robust. In practical situations, one doesn’t need such extensive iterations, as can be seen from Fig. 6, with 5 iteration steps, the relative error for the projected electron density drops from $3.9\%$ to $1.3\%$ for the five-image input test (the blue curve); while for the two-image input test, the relative error drops from $5.5\%$ to $3.1\%$ after 5 iteration steps. After that, the improvements are almost negligible. Figure 6 shows again that better retrieving results are expected with more input images.

Fig. 6. The log-log plot of the relative errors of the retrieved projected electron density $\rho _{\mathrm {ep}}$ in each iteration step. The solid blue curve corresponding to the five-image input test, while the dashed red curve corresponding to the two-image input test.

Download Full Size | PDF

Table 3. Accuracy of reconstructed projected electron density $\rho _{\mathrm {ep}}$. The input phase contrast images are $I_{g1}=\{I_{\mathrm {16-40}}\}$, $I_{g2}=\{I_{\mathrm {41-100}}\}$.

View Table | View all tables in this article

Table 4. Contrast Noise Ratio (CNR) of the input phase contrast images and the reconstructed attenuation images. $I_{g1}=\{I_{\mathrm {16-40}}\}$, $I_{g2}=\{I_{\mathrm {41-64}}\}$.

View Table | View all tables in this article

6. Discussion and conclusion

Phase retrieval is a crucial step for quantitative imaging such as reconstructing the 3-D distribution of tissue linear attenuation coefficients and refraction indices. However, images acquired in medical applications are relatively noisy due to the radiation dose constraints. Thus, an accurate and robust algorithm to reduce quantum noise that appears in acquired images is essential. In [16], the authors developed a robust algorithm (PAD) that can accurately retrieve the tissue’s projected electron density. To apply this algorithm, a fundamental condition, the phase-attenuation duality condition must be satisfied. Fortunately for soft tissues, which are mainly composed of light elements($Z<10$), when X-ray energies are 60 keV or higher, the PAD condition is satisfied. For dense materials, such as bone or lower X-ray energies, this PAD condition fails to be fulfilled. In general, one needs at least two phase contrast images at different X-ray energies or different geometric setups to completely retrieve the phase values. Images with different geometric setups are usually not permissible in the medical imaging since they tend to increase the X-ray radiation. However, with the advancements made in the photon counting detection technology, multiple X-ray images at different energy levels from a single X-ray exposure can be acquired. Thus, without increasing radiation dose levels, accurate reconstruction of the 3-D distribution of tissue linear attenuation coefficients and refraction indices becomes possible.

In this paper, the authors proposed an iterative method of simultaneously retrieving tissue projected electron density $\rho _{\mathrm {ep}}$ and the ray-integral of the linear attenuation coefficient $\mu _0$, at energy $E_0=1$ keV, generated by PE absorption using two or more phase contrast images at different energy levels. By factoring out the attenuation caused by Compton scattering from the attenuation map, and with a restraint of small fluctuations of the attenuation caused by PE absorption, we separate the phase contrast image into two factors: One is a pseudo phase contrast image caused by electron density, the other is the attenuation map caused by PE absorption, Eq. (8). After this decomposition, the noise reduction capability of the PAD algorithm allows the implementation of a stable iterative algorithm. With a multimaterial phantom imaged with a PCD based X-ray inline phase sensitive imaging prototype, we demonstrate the algorithm’s accuracy and its robustness against noise. It was shown that the accuracy of the retrieval and improvement in CNR with five-image input was higher than the two-input images. It is worth mentioning that all the images were acquired using the spectroscopic mode of the detector. This was mainly due to the presence of only two energy discrimination thresholds in the imaging mode, limiting the ability to acquire more than two images. Unlike the five-input images, the two-input images can also be acquired using the imaging mode of the PCD by setting the first threshold (TH0) to 16 keV, and second threshold (TH1) to 40 keV. Thus, the practicality of the two-input images for phase retrieval suits well the current PCD in the imaging mode. Future investigations will be performed using the imaging mode of the PCD to further test the proposed method for accurate retrieval purposes and improved image quality.

In the simulation, we assumed the linear attenuation arising from PE is inversely proportional to the cubic of X-ray energy, Eq. (4). This is an approximation. The true relationship is material dependent. An incorrect choice of the index number will certainly affect the quality of the results. Also, for images with multiple materials, a universal index number may also affect the results. In our experimental test, one can see that the edge enhancements in the retrieved projected electron densities $\rho _{\mathrm {ep}}$ (Fig. 4(top)) are completely removed. But for the $\mu _0$ image, the edge enhancements are still present. We believe this is due to Step 4 in the algorithm where although the edge enhancements can be removed partially by the division $I_2(E_2)/{I_{\mathrm {kn}}}(E_2)$, the difference in experimental data and software simulation cannot remove all of them. The noise in the experimental data is also piled up in the iterations. That is why the retrieved $\mu _0$ image is noisier than the $\rho _{\mathrm {ep}}$ image. When converting to attenuation images the noise and edge enhancements present in $\mu _0$ will add up to the attenuation images, especially for low energy images. This is because the composition of the attenuation arises from PE is weighted more compared to images of higher energy attenuation images. Further investigation is needed about this issue.

This algorithm does not require a priori information about the composition of materials within the object. As a result, the algorithm can be applied universally to any sample X-ray image, including dense materials, for a wide X-ray ranges. Once the projected electron density $\rho _{\mathrm {ep}}$ and the ray-integral of the linear attenuation coefficient $\mu _0$ at energy $E_0=1$ keV, caused by PE absorption, are reconstructed, further analysis based on attenuation or electron density, such as material decomposition, or 3-D quantitative analysis (combined with CT reconstruction) can be performed.

Funding

National Institutes of Health (F32CA250300, P30CA225520).

Disclosures

The authors declare no conflicts of interest.

Data availability

The data that support the findings of this study are available from the corresponding author upon reasonable request.

References

1. A. Snigirev, I. Snigireva, V. Kohn, S. Kuznetsov, and I. Schelokov, “On the possibilities of x-ray phase contrast micro-imaging by coherent high-energy synchrotron radiation,” Rev. Sci. Instrum. 66(12), 5486–5492 (1995). [CrossRef]

2. S. Wilkins, T. Gureyev, D. Gao, A. Pogany, and A. Stevenson, “Phase contrast imaging using polychromatic hard x-ray,” Nature 384(6607), 335–338 (1996). [CrossRef]

3. K. A. Nugent, “Coherent methods in the x-ray sciences,” Adv. Phys. 59(1), 1–99 (2010). [CrossRef]

4. T. J. Davis, D. Gao, T. E. Gureyev, A. W. Stewenson, and S. W. Wilkins, “Phase-contrast imaging of weakly absoprbing materials using hard x-rays,” Nature 373(6515), 595–598 (1995). [CrossRef]

5. P. Cloetens, J. P. Guigay, C. D. Martino, J. Baruchel, and M. Schlenker, “Fractional talbot imaging of phase gratings with hard x rays,” Opt. Lett. 22(14), 1059–1061 (1997). [CrossRef]

6. D. Paganin and K. Nugent, “Noninterferometric phase imaging with partially coherent light,” Phys. Rev. Lett. 80(12), 2586–2589 (1998). [CrossRef]

7. A. Pogany, D. Gao, and S. Wilkins, “Contrast and resolution in imaging with a microfocus x-ray source,” Rev. Sci. Instrum. 68(7), 2774–2782 (1997). [CrossRef]

8. F. Arfelli, V. Bonvicini, A. Bravin, G. Cantatore, E. Castelli, L. D. Palma, M. D. Michiel, M. Fabrizioli, R. Longo, R. H. Menk, A. Olivo, S. Pani, D. Pontoni, P. Poropat, M. Prest, A. Rashevsky, M. Ratti, L. Rigon, G. Tromba, A. Vacchi, E. Vallazza, and F. Zanconati, “Mammography with synchrotron radiation: phase-detection techniques,” Radiology 215(1), 286–293 (2000). [CrossRef]

9. D. Paganin, S. Mayo, T. Gureyev, P. Miller, and S. Wilkins, “Simultaneous phase and amplitude extraction from a single defocused image of a homogeneous object,” J. Microsc. 206(1), 33–40 (2002). [CrossRef]

10. S. Mayo, T. Davis, T. Gureyev, P. Miller, D. Poganin, A. Pogany, A. Stevenson, and S. Wilkins, “X-ray phase-contrast microscopy and microtomography,” Opt. Express 11(19), 2289–2302 (2003). [CrossRef]

11. X. Wu and H. Liu, “A general theoretical formalism for x-ray phase contrast imaging,” J. X-Ray Sci. Technol. 11(1), 33–42 (2003).

12. X. Wu and H. Liu, “Clinical implementation of phase-contrast x-ray imaging: Theoretical foundations and design considerations,” Med. Phys. 30(8), 2169–2179 (2003). [CrossRef]

13. X. Wu and H. Liu, “A new theory of phase-contrast x-ray imaging based on wigner distributions,” Med. Phys. 31(9), 2378–2384 (2004). [CrossRef]

14. E. Donnelly, R. Price, and D. Pickens, “Experimental validation of the wigner distributions theory of phase-contrast imaging,” Med. Phys. 32(4), 928–931 (2005). [CrossRef]

15. D. Zhang, M. Donvan, L. Fajardo, A. Archer, X. Wu, and H. Liu, “Preliminary feasibility study of an in-line phase contrast x-ray imaging prototype,” IEEE Trans. Biomed. Eng. 55(9), 2249–2257 (2008). [CrossRef]

16. X. Wu, H. Liu, and A. Yan, “X-ray phase-attenuation duality and phase retrieval,” Opt. Lett. 30(4), 379–381 (2005). [CrossRef]

17. X. Wu and H. Liu, “X-ray cone-beam phase tomography formulas based on phase-attenuation duality,” Opt. Express 13(16), 6000 (2005). [CrossRef]

18. A. Bravin, P. Coan, and P. Suortti, “X-ray phase-contrast imaging: from pre-clinical applications towards clinics,” Phys. Med. Biol. 58(1), R1–R35 (2013). [CrossRef]

19. P. Tafforeau, R. Boistel, E. Boller, A. Bravin, M. Brunet, Y. Chaimanee, P. Cloetens, M. Feist, J. Hoszowska, J. J. Jaeger, R. Kay, V. Lazzari, L. Marivaux, A. Nel, C. Nemoz, X. Thibault, P. Vignaud, and S. Zabler, “Applications of x-ray synchrotron microtomography for non-destructive 3d studies of paleontological specimens,” Appl. Phys. A 83(2), 195–202 (2006). [CrossRef]

20. X. Wu and A. Yan, “Phase retrieval from one single phase contrast x-ray image,” Opt. Express 17(13), 11187–11196 (2009). [CrossRef]

21. L. J. Allen and M. P. Oxley, “Phase retrieval from series of images obtained by defocus variation,” Opt. Commun. 199(1-4), 65–75 (2001). [CrossRef]

22. A. Burvall, U. Lundström, P. A. C. Takman, D. H. Larsson, and H. M. Hertz, “Phase retrieval in x-ray phase-contrast imaging suitable for tomography,” Opt. Express 19(11), 10359 (2011). [CrossRef]

23. M. Beltran, D. Paganin, K. Uesugi, and M. Kitchen, “2d and 3d x-ray phase retrieval of multi-material objects using a single defocus distance,” Opt. Express 18(7), 6423 (2010). [CrossRef]

24. L. C. P. Croton, K. S. Morgan, D. M. Paganin, L. T. Kerr, M. J. Wallace, K. J. Crossley, S. L. Miller, N. Yagi, K. Uesugi, S. B. Hooper, and M. J. Kitchen, “In situ phase contrast x-ray brain ct,” Sci. Rep. 8(1), 11412 (2018). [CrossRef]

25. D. Gürsoy and M. Das, “Single-step absorption and phase retrieval with polychromatic x rays using a spectral detector,” Opt. Lett. 38(9), 1461 (2013). [CrossRef]

26. K. Taguchi and J. S. Iwanczyk, “Vision 20/20: Single photon counting x-ray detectors in medical imaging,” Med. Phys. 40, 100901 (2013). [CrossRef]

27. M. J. Willemink, M. Persson, A. Pourmorteza, N. J. Pelc, and D. Fleischmann, “Photon-counting ct: Technical principles and clinical prospects,” Radiology 289(2), 293–312 (2018). [CrossRef]

28. L. Ren, B. Zheng, and H. Liu, “Tutorial on x-ray photon counting detector characterization,” J. X-Ray Sci. Technol. 26(1), 1–28 (2018). [CrossRef]

29. T. E. Gureyev, S. Mayo, S. W. Wilkins, D. Paganin, and A. W. Stevenson, “Quantitative in-line phase-contrast imaging with multienergy x-rays,” Phys. Rev. Lett. 86(25), 5827–5830 (2001). [CrossRef]

30. S. C. Mayo, P. R. Miller, S. W. Wilkins, T. J. Davis, D. Gao, T. E. Gureyev, D. Paganin, D. J. Parry, A. Pogany, and A. W. Stevenson, “Quantitative x-ray projection microscopy: phase-contrast and multi-spectral imaging,” J. Microsc. 207(2), 79–96 (2002). [CrossRef]

31. I. Vazquez, I. E. Harmon, J. C. R. Luna, and M. Das, “Quantitative phase retrieval with low photon counts using an energy resolving quantum detector,” J. Opt. Soc. Am. A 38(1), 71–79 (2021). [CrossRef]

32. R. E. Alvarez and A. Macovski, “Energy-selective reconstructions in x-ray computerised tomography,” Phys. Med. Biol. 21(5), 733–744 (1976). [CrossRef]

33. F. Schaff, K. S. Morgan, D. M. Paganin, and M. J. Kitchen, “Spectral x-ray imaging: Conditions under which propagation-based phase-contrast is beneficial,” Phys. Med. Biol. 65(20), 205006 (2020). [CrossRef]

34. H. T. Li, A. M. Kingston, G. R. Myers, L. Beeching, and A. P. Sheppard, “Linear iterative near-field phase retrieval (lipr) for dual-energy x-ray imaging and material discrimination,” J. Opt. Soc. Am. A 35(1), A30 (2018). [CrossRef]

35. H. T. Li, F. Schaff, L. C. P. Croton, K. S. Morgan, and M. J. Kitchen, “Quantitative material decomposition using linear iterative near-field phase retrieval dual-energy x-ray imaging,” Phys. Med. Biol. 65(18), 185014 (2020). [CrossRef]

36. X. Wu and H. Liu, “A dual detector approach for x-ray attenuation and phase imaging,” J. X-Ray Sci. Technol. 12(1), 35–42 (2004).

37. X. Wu and H. Liu, “Phase-space evolution of x-ray coherence in phase-sensitive imaging,” Appl. Opt. 47(22), E44–E52 (2008). [CrossRef]

38. M. R. Teague, “Deterministic phase retrieval: a green’s function solution,” J. Opt. Soc. Am. 73(11), 1434–1441 (1983). [CrossRef]

39. K. A. Nugent, T. E. Gureyev, D. F. Cookson, D. Paganin, and Z. N. M. N. Barnea, “Quantitative phase imaging using hard x rays,” Phys. Rev. Lett. 77(14), 2961–2964 (1996). [CrossRef]

40. J. Guigay, M. Langer, R. Boistel, and P. Cloetens, “Mixed transfer function and transport of intensity approach for phase retrieval in the fresnel region,” Opt. Lett. 32(12), 1617–1619 (2007). [CrossRef]

41. T. Gureyev, Y. Nesterets, D. Paganin, A. Pogany, and S. Wilkins, “Linear algorithms for phase retrieval in the fresnel region. 2. partially coherent illumination,” Opt. Commun. 259(2), 569–580 (2006). [CrossRef]

42. N. A. Dyson, X-rays in atomic and nuclear physics (Cambridge University Press, 1990).

43. X. Wu, A. Dean, and H. Liu, Biomedical Photonics Handbook, chap. 26 (CRC Press, Tampa, Fla., 2003).

44. F. Schaff, K. S. Morgan, J. A. Pollock, S. B. H. Hooper, L. C. P. Croton, and M. J. Kitchen, “Material decomposition using spectral propagation-based phase contrast x-ray imaging,” IEEE Trans. Med. Imaging 39(12), 3891–3899 (2020). [CrossRef]

45. X. Wu and H. Liu, “Clarification of aspects in in-line phase-sensitive x-ray imaging,” Med. Phys. 34(2), 737–743 (2007). [CrossRef]

	Retrieved $ρ_{e p}$	Expected $ρ_{e p}$	Difference (%)	Rod only $ρ_{e p}$	Expected $ρ_{e p}$	Difference (%)
	(electrons/m $^{2}$ )	(electrons/m $^{2}$ )		(electrons/m $^{2}$ )	(electrons/m $^{2}$ )
	$\times 10^{27}$	$\times 10^{27}$		$\times 10^{26}$	$\times 10^{26}$
Background (Acrylic)	3.028	2.999	0.973	–	–	–
Rod 1 (Adipose)	3.695	3.710	−0.404	6.735	7.103	−5.182
Rod 2 (25mg/cc)	3.839	3.811	0.729	8.007	8.118	−1.371
Rod 3 (Muscle)	3.900	3.793	2.807	8.516	7.941	7.238
Rod 4 (Lung)	3.191	3.150	1.307	1.586	1.510	5.076
Rod 5 (Breast)	3.761	3.733	0.740	7.327	7.337	−0.143

energy		Rod 1:	Rod 2:	Rod 3:	Rod 4:	Rod 5:
(keV)		(Adipose)	(25mg/cc)	(Muscle)	(Lung)	(Breast)
24.66	$I$	0.84	1.00	1.12	0.20	0.93
	$A^{2}$	0.97	1.20	1.39	0.22	1.11
	ratio	1.16	1.21	1.24	1.14	1.19
34.22	$I$	1.04	1.29	1.45	0.25	1.17
	$A^{2}$	2.27	2.77	3.09	0.55	2.53
	ratio	2.19	2.15	2.12	2.21	2.16
45.74	$I$	0.96	1.19	1.31	0.23	1.08
	$A^{2}$	5.17	6.29	6.87	1.27	5.72
	ratio	5.39	5.29	5.26	5.55	5.28
57.92	$I$	0.98	1.17	1.26	0.26	1.04
	$A^{2}$	9.06	11.00	11.92	2.23	9.98
	ratio	9.24	9.37	9.47	8.63	9.58

	Retrieved $ρ_{e p}$	Expected $ρ_{e p}$	Difference (%)	Rod only $ρ_{e p}$	Expected $ρ_{e p}$	Difference (%)
	(electrons/m $^{2}$ )	(electrons/m $^{2}$ )		(electrons/m $^{2}$ )	(electrons/m $^{2}$ )
	$\times 10^{27}$	$\times 10^{27}$		$\times 10^{26}$	$\times 10^{26}$
Background (Acrylic)	3.085	2.999	2.868	–	–	–
Rod 1 (Adipose)	3.767	3.710	1.559	6.888	7.103	−3.034
Rod 2 (25mg/cc)	3.919	3.811	2.821	8.238	8.118	1.472
Rod 3 (Muscle)	3.984	3.793	5.024	8.800	7.941	10.816
Rod 4 (Lung)	3.252	3.150	3.238	1.661	1.510	10.021
Rod 5 (Breast)	3.840	3.733	2.856	7.525	7.337	2.564

energy		Rod 1:	Rod 2:	Rod 3:	Rod 4:	Rod 5:
(keV)		(Adipose)	(25mg/cc)	(Muscle)	(Lung)	(Breast)
24.66	$I$	0.84	1.00	1.12	0.20	0.93
	$A^{2}$	0.92	1.12	1.31	0.22	1.04
	ratio	1.10	1.13	1.17	1.10	1.12
34.22	$I$	1.04	1.29	1.45	0.25	1.17
	$A^{2}$	2.17	2.63	2.95	0.53	2.42
	ratio	2.09	2.05	2.03	2.13	2.07
45.74	$I$	0.96	1.19	1.31	0.23	1.08
	$A^{2}$	5.04	6.13	6.72	1.24	5.59
	ratio	5.25	5.16	5.15	5.44	5.16
57.92	$I$	0.98	1.17	1.26	0.26	1.04
	$A^{2}$	8.84	10.76	11.69	2.19	9.79
	ratio	9.02	9.17	9.30	8.45	9.40

	Retrieved $ρ_{e p}$	Expected $ρ_{e p}$	Difference (%)	Rod only $ρ_{e p}$	Expected $ρ_{e p}$	Difference (%)
	(electrons/m $^{2}$ )	(electrons/m $^{2}$ )		(electrons/m $^{2}$ )	(electrons/m $^{2}$ )
	$\times 10^{27}$	$\times 10^{27}$		$\times 10^{26}$	$\times 10^{26}$
Background (Acrylic)	3.028	2.999	0.973	–	–	–
Rod 1 (Adipose)	3.695	3.710	−0.404	6.735	7.103	−5.182
Rod 2 (25mg/cc)	3.839	3.811	0.729	8.007	8.118	−1.371
Rod 3 (Muscle)	3.900	3.793	2.807	8.516	7.941	7.238
Rod 4 (Lung)	3.191	3.150	1.307	1.586	1.510	5.076
Rod 5 (Breast)	3.761	3.733	0.740	7.327	7.337	−0.143

Dual-energy phase retrieval algorithm for inline phase sensitive x-ray imaging system

Abstract

1. Introduction

2. Method

3. Dual-energy phase reconstruction algorithm

4. Experimental validation

5. Results

6. Discussion and conclusion

Funding

Disclosures

Data availability

References

Data availability

Cited By

Figures (6)

Tables (4)

Equations (21)

Optics Express