1. Introduction

Conventional x-ray imaging for soft tissues (such as breast, brain, liver, etc.) is limited in its sensitivity for detecting subtle tissues pathological changes, since the imaging relies on the small difference in x-ray attenuation between the lesions and soft tissues of variable structure. However, as x-ray wave passes through tissues, x-ray undergoes phase-shifts as well. The amount of x-ray phase shift ϕ by tissues is

where re is the classical electron radius, h is the Plank constant and c the speed of light, E is x-ray photon energy. Here ρe denotes tissue electron density, and the integration of electron density ρ_e is over the x-ray path and the integral is the projected electron density ρ_e,p [1, 2, 3]. It turns out that x-ray phase-shift differences between tissue and lesions are about one thousand times larger than their attenuation differences [1, 2, 3]. Hence the phase-sensitive imaging has the potential to greatly enhance the lesion detection sensitivity and specificity, and reduce the radiation dose associated with x-ray imaging. In the inline phase contrast imaging as studied in previous works, 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, and Laplacian and gradients of x-ray phase-shifts caused by the object, and gradients of the objects attenuation [4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14]. The procedure of disentangling tissue phase-shifts from the mixed contrast mechanism and retrieving the phase maps from acquired phase contrast images is called the phase retrieval. Phase retrieval technique plays a central role in phase-sensitive x-ray imaging. By means of phase retrieval, one can reconstruct a quantitative map of phase-shifts, a phase image of the imaged object [3, 6, 13, 15, 16]. According to Eq. (1), a retrieved phase image is equivalently a map of imaged objects quantitative projected electron densities. Therefore, tissue phase images can provide quantitative tissue characterization, which potentially can increase the sensitivity and specificity of a diagnosis. Moreover, phase retrieval is a necessary procedure for phase-sensitive volumetric imaging such as the phase-sensitive tomography and tomosynthesis [7, 14, 15]. In the phase-sensitive imaging, x-ray does NOT propagate along straight lines because of x-ray refraction and diffraction at tissue interface and boundaries. The diffraction fringes generated by phase contrast at tissue interfaces will obviously render conventional volumetric image reconstruction methods erroneous, since conventional volumetric reconstruction methods assume that x-rays propagate along straight lines. The phase-sensitive CT experiments showed that conventional tomography reconstruction algorithms applied to phase-contrast projections can lead to erroneous tomography images with artifacts such as the streaks and negative densities [7]. The exact mechanism of the phenomenon was revealed in one of our recent work. In fact, if we denote the imaged objects linear attenuation coefficient by µ(r⃗) and its electron density by ρ_e(r⃗), with the conventional filtered backprojection method the reconstructed “apparent attenuation coefficients” µ _recon(r⃗) was found as [17]:

where R ₂ denotes the object-detector distance, λ the x-ray wavelength. Eq. (2) shows that the reconstructed tomogram consists of three sub-tomograms. The first sub-tomogram is the tomogram of the tissue linear attenuation coefficients µ(r⃗). The second sub-tomogram is a scaled map of the 3-D Laplacian of tissue electron densities. The third sub-tomogram µ _mixed(r⃗) represents the artifacts reflecting the global variation of both µ(r⃗) and ρ_e(r⃗), but is not representing the local variations of µ(r⃗) and ρ_e(r⃗). It is this term that causes artifacts such as the streaks in reconstruction images. In order to correctly deal with volumetric reconstruction issue, as a general strategy one should first disentangle tissue phase-shifts from the mixed contrast mechanism from acquired phase-sensitive projections and retrieve the phase maps for all projections, then reconstruct the 3D images of object attenuation coefficients and electron densities, respectively.

For phase retrieval one in general needs to acquire at least two and mostly four to five images with different object-detector distances [3, 15]. This requirement is imposed by the x-ray wave nature and its wave equation [3, 8, 9, 10]. Multiple image-acquisitions for a phase retrieval increase radiation dose associated with the technique. For phase sensitive tomography or tomosynthesis, multiple image-acquisitions for each projection-angle make phase-sensitive tomography and tomosynthesis implementation especially cumbersome: multiple phase-contrast images have to be acquired with different object-detector distances for all the projections during a scan, and the number of required projections can be as high as one thousand. In addition, the image data alignment procedure is quite time-consuming. Moreover, radiation dose is multiplied for multiple image acquisition per projection view as well. Therefore, an important challenge is to find an effective and low-dose approach for phase retrievals for phase-sensitive imaging. Many studies were conducted addressing this challenge [6, 13, 14]. We proposed a unique solution for phase-sensitive imaging with high energy x-rays. We noted that with increasing x-ray energy tissues photoelectric absorption decreases as 1/Eⁿ with n≈3-4. but tissues phase shift decreases much slower, and only as 1/E as is shown in Eq. (1). Therefore, phase-sensitive imaging with x-rays can fully take advantage of phase contrast and greatly reduce radiation doses associated with x-ray imaging, by taking advantage of the dual relation between the high energy x-rays phase-shifts and attenuation for soft tissues. Furthermore, analyzing the tissue attenuation data, we observed that the soft tissue attenuation cross sections are very well approximate by that of x-ray incoherent scattering for x-rays of about 60 keV to 500 keV. Under this circumstance, we noted that soft-tissue attenuation and soft-tissue phase are all related to the projected electron density, we called this new notion as the phase-attenuation duality. Based on this concept of phase-attenuation duality, we found a way to retrieve a phase-map from only one single recorded image, that is, the retrieved phase map is given by [13]:

where R ₁ is the distance from the point x-ray source to object, R ₂ is the distance from object to detector, M=(R ₁+R ₂)/R ₁ is the geometric magnification, r⃗ is the position vector (see Fig. 1.) In Eq. 3, I(Mr⃗;R ₁+R ₂) is the measured image intensity at detector plane, and I _in is the intensity at entrance, the operator ∇² denotes the Laplacian operator. In Eq. (3) σ_KN is the Klein-Nishina total cross section, that is, the total cross section for x-ray Compton scattering with a free electron [17]. Note that σ_KN is x-ray energy dependent, as is shown below. We point out that the high sensitivity of x-ray phase change is realized as the high ratio λr_e/σ_KN≫1 in this circumstance. In fact, for clinical imaging photon energy is no higher than 150 keV, and in this case λr_e/σ_KN~10³. In other words, with this duality the high sensitivity of phase imaging manifests itself as the high ratio of x-ray wavelength over the classical electron radius. Compared to the common multiple-acquisition based phase retrieval approaches in the literature [3, 6, 15], our duality-based phase retrieval method is a conveniently implemented method with tremendous dose-saving for many potential applications such as mammography. These advantages will be critical for phase-sensitive volumetric imaging such as breast tomosynthesis and computed tomography.

 

Fig. 1. A shematic diagram for the measuring system

Download Full Size | PDF

However, there is a limitation on applicability of the phase retrieval formula Eq. (3). In fact. Eq. (3) is not valid for the cases with very high spatial resolution and with large object-detector distances. The physics behind this can be intuitively explained as follows. In fact, in inline phase-sensitive imaging the phase contrast is formed by the diffraction from the object to the detector of x-rays with shifted phases. The Fresnel diffraction can be thought as a x-ray wave-front filtering with the Fresnel propagator exp(iπλ zu⃗²), where λ is the wavelength and u⃗ is the wave-front transverse spatial frequency vector, and z=R ₂/M in terms of the parameters in Eq. (3). So the phase of the propagator, i.e. πλzu⃗², determines how strongly the wave-front exit from the object will be diffracted in x-ray propagation. The phase retrieval formula Eq. (3) was derived by linearizing Fresnel propagator phase under the assumption of πλ zu⃗²≪1 [8, 13]. This linearizing is equivalent to simplifying x-ray wave equation of motion to the transport of intensity of equation (TIE) [8]. However, for cases with very high spatial resolution and with large object-detector distances, such as cases with phase-sensitive imaging experiments with synchrotron based x-ray sources, or the high resolution imaging experiments with biological specimen or small animals, the high spatial resolution or large propagation distances make Fresnel propagator phases πλzu⃗² larger than 1 for many transverse spatial frequencies of x-ray wave-front. For example, for 60 keV x-rays, if M=1 and R ₂=2m, detector pixel is 5 µm, then πλzu⃗²=1.299. For these cases one cannot linearize Fresnel propagators and should consider the full exponential propagators instead. In section 2, we will first briefly review the concept of the phase-attenuation duality and how it greatly simplify the phase retrieval for soft-tissue objects. We will then derive a new and general PA-Duality based phase retrieval formula, as the main result of this work. With this general formula, one can retrieve the phase from one single phase contrast image and is valid for cases involved with large Fresnel propagator phases. Hence this new phase retrieval formula extends its applicability to cases with very high spatial resolution and with large object-detector distances. In section 3, we show the simulation results with the general formula for phase retrieval.

2. Phase-attenuation duality and phase retrieval from one single phase contrast image for soft tissue objects

As is explained in previous works [3, 4, 8, 9, 10], in phase-sensitive imaging the imaged object can be modeled by a two-dimensional transmission function $T\left(\overrightarrow{r}\right)=A_0\left(\overrightarrow{r}\right)e^{i{\varphi }_0\left(\overrightarrow{r}\right)}$, where ϕ ₀(r⃗) denotes the x-ray phase-shift caused by the object, and A ₀(r⃗) denotes the x-ray transmission, or the attenuation-map of object. For sake of clarity in exploring the effects of non-linearity of the Fresnel propagators, we assume the x-ray source is a quasi-monochromatic point source in following derivations. We will discuss more general x-ray sources in section 3.

Let R ₁ denote the source-object distance and R ₂ the object-detector distance. The geometric magnification factor M is equal to (R ₁+R ₂)/R ₁. We also assume that one encounters only moderate variation of the object phase function ϕ(r⃗) such that |ϕ(r⃗+λR ₂ u⃗/2M)−ϕ(r⃗−λR ₂ u⃗/2M)|≪1 [3]. Transmitting from the object, x-rays undergo attenuation and phase-shift, and diffract over a distance R2 to the detector. An important task is to find out the formula of the image intensity as function of the objects transmission function and imaging geometry. A general solution had been found in one of our previous works [3]. Starting from the paraxial Fresnel-Kirchhoff diffraction theory, we found out in this work that the Fourier transform of x-ray irradiance at detector plane is given by:

where I(r⃗_D) denotes the x-ray irradiance at the detector entrance, and the symbol 𝓕̂ (·) denotes the 2-D Fourier transform and I _in is the intensity of the incident x-ray upon the object, and u⃗ is the spatial frequency vector in the object plane. Note that the Fresnel propagator phases πλR ₂ u⃗²/M enter Eq. (4) as the argument of sinusoidal transfer functions. For cases with moderate spatial resolutions u⃗ and R ₂, Eq. (4) is much simplified, and the inverse Fourier transformation of the simplified Eq. (4) can be written as [8]:

Equation (5) is the TIE-based formula for image irradiance, which was first derived by Nugent and colleagues[3]. Recently Guigay et al showed the limitations of Eq. (5) with theoretical analysis and experiments [17]. They found that for in-line phase-sensitive imaging with long R ₂ and high resolutions such that πλR ₂ u⃗²/M>1, the phase retrieval based on Eq. (5) would be inaccurate, one has to use the general formula Eq. (4) for accurate quantitative phase imaging.

Furthermore, we note that soft tissues encountered in clinical imaging are composed almost exclusively by atoms of the light elements with the atomic numbers Z<10. For example, the composition of breast tissue consists mainly of light elements such as hydrogen, carbon, nitrogen and oxygen. In contrast, the sum of weight fractions of other heavier elements in breast tissue is only less than 1% [18]. The same is true for other soft tissues such as the brain grey/white matter [18]. Analyzing the tissue attenuation data, we observed that the soft tissue attenuation cross sections are very well approximated by that of x-ray incoherent scattering for x-rays of about 60keV to 500keV[19]. In addition, the incoherent scattering functions become linearized as well [19]. Under this circumstance, while the phase-shift is caused by the coherent x-ray scattering, but soft tissue attenuation is determined by the incoherent x-ray scattering for soft tissues, so both the x-ray phase-shift and attenuation by soft tissues are all determined by the projected electron density for these high energy x-rays. We called this complementary relationship between phase-shift and attenuation for soft tissues as the phase-attenuation duality. When the P-A duality holds, x-ray attenuation and phase shift by the object is related in following way:

where σ_KN is the total cross section for x-ray Compton scattering with a free electron:

where η=E _photon/m_ec ², here we denote the photon energy of the primary x-ray beam by E _photon, and m_ec ² is the rest electron energy and equal to 511keV. The duality condition can be used to simplify Eq. (4). In fact, we can rewrite Eq. (4) as

Especially with the duality condition we found

In Eq. (9) we made an approximation of dropping the negligible contribution of the higher order terms of σ_KN/λr_e. The approximation is based on σ_KN/λr_e≈10⁻³, and we assumed |ϕ(r⃗+λR ₂ u⃗/2M)−ϕ (r⃗−λR ₂ u⃗/2M)|≪1, moderate variations of object phase shifts. Substituting Eqs. (9) into Eq. (8), we found:

Equation (10) is the central result of this work, it is the formula for x-ray irradiance for high resolution soft tissue imaging with high energy x-rays. In the derivation the parameter λr_e/σ_KN plays a significant role. It is interesting to note that Gureyev et al obtained a similar result in ref.[20] (Eq.(24)) with the exception that the term πλR ₂ u⃗²/M in Eq. (10) was absent. The differences may lie in the different assumptions made regarding the object transmission functions and x-ray energies. Gureyev and coauthors assumed the objects as single material homogeneous objects [20], but our results are for heterogeneous objects of the light elements imaged with high energy x-rays (60 keV or higher), i.e., for the cases the phase-attenuation duality holds. In addition, in their derivation they separated the slowly and rapidly varying components of the amplitude and phase, and they assumed that the magnitude of the rapidly varying components are much smaller than 1, but we did not make such limiting assumption in derivation of Eq. (10). We stress that the Klein-Nishina cross-section is energy dependent shown in Eq. (7) and λr_e/σ_KN is about 1000 for inline x-ray energy range used in diagnostic imaging. This large λr_e/σ_KN represents the high sensitivity of phase contrast in tissue differentiation, but it also makes the x-ray irradiance formula Eq. (10) getting very much simplified. It is easy to verify that for the low resolution cases such as those with πλR ₂ u⃗²/M≪1, Eq. (10) is reduced to

we recover our previous result as shown in [13]. On the other hand, for general cases πλR ₂ u⃗²/M~1 or πλR ₂ u⃗²/M>1 Eq. (11) cannot be applied, and one should use Eq. (10) for phase retrieval as is shown below.

3. Simulation results

Comparing Eq. (10) and Eq. (11) for use in performing phase retrieval, the only difference lies in the inverse filters derived for retrieving 𝓕̂ (A ² ₀). For phase retrieval based on Eq. (10) the inverse filter is found as

where α=πλR ₂ u⃗²/M, the FT-space Fresnel Propagator phase as defined in section 1, and γ=λr_e/σ_KN, which represents the phase contrast enhancement as explained earlier. $\sqrt{1+{\left(2\gamma +\alpha \right)}^2}$ and α ₁=α−sin⁻¹[(2γ+α)/δ]. On the other hand, for phase retrieval based on Eq. (11) the inverse filter will be

Since γ~10³, the filters in Eq. (12) and Eq. (13) are close when α is small. When α is getting bigger, close to or greater than π (α ₁ close to π/2), the difference between the two inverse filters is getting bigger. Figure 2 shows a plot of the two filters, with respect to the change of α. The inverse filter in Eq. (12), depicted in Fig. 2 as the solid line, changes the sign when α passes through kπ (and thus α ₁ passes through (k−1/2)π) and diverges to ∞ when α ₁=(k−1/2)π, k=1,2,⋯. On the other hand, the filter in Eq. (13), depicted as the dashed line in Fig. 2, keeps positive and approaches zero as α gets bigger. So the difference of the retrieved phase maps using Eq. (10) and Eq. (11) will be highly dependent on the value of the Fresnel propagator phase α=πλR ₂ u⃗²/M.

To demonstrate this phenomenon, we performed two simulations. First, we set the phase shift ϕ to be a sine-function in x direction with period of 32 pixel size and constant in y direction. Assuming the pixel pitch is p µm, the spatial frequency corresponding to ϕ, in x direction, is µ _x0=1/32p 1/µm. In our simulations, for the purpose of simplification, the magnification M as well as the intensity at entrance, I _in, are assumed to be 1. The x-ray photon energy is set to 60 keV. This gives the x-ray wavelength λ=2.0667×10⁻¹¹ m and the total cross section for x-ray Compton scattering with a free electron σKN=5.4526×10⁻²⁹ m² as determined from Eq. (7). The classical electron radius r_e is set to 2.82×10⁻¹⁵ m. By assuming the detectors resolution p=1 µm and the phase shift ϕ(x,y)=0.5sin(32px)-1, according to the phase-attenuation duality discussed earlier, the attenuation map is related to the phase map as A ² ₀ (x,y)=exp[σ_KN/λr_e·ϕ]=exp[(sin(32px)-2)/2γ], with γ=λr_e/σ_KN=1.0726×10³. We then performed Fresnel diffraction to determine the phase-contrast image irradiances at distances downstream[21], and the retrieved maps of A ² ₀ were found by using the inverse filters of Eq. (12) and Eq. (13), which were derived based on Eq. (10) and Eq. (11), respectively. Because of the large value of γ, phase contrast image irradiances {I} were found be maps with each I being dominated by two single frequency component at u_x=±u _x0 (u _x0=1/32p=1/32×10⁶ m) plus a zero frequency component. This being so, the shape and variation of the retrieved A20 - maps were determined by their dominant frequency component at frequencies u_x=±u _x0. In order to clearly show the differences between phase retrievals based on Eq. (10) and Eq. (11), we examined the ratio of the dominant frequency components of the two retrieved A ² ₀ -maps based on Eq. (10) and Eq. (11), respectively. We found the ratio being r=(1+2γα)/(δ · cos(α ₁)) for a given object-detector distance R ₂. This ratio reflects the capability of achieving accurate phase retrieval with Eq. (10) against Eq. (11). For small α≪1, that is, for small Fresnel propagator phases, this ratio is 1, so Eq. (10) against Eq. (11) both are equally good for the accurate phase retrieval. However, when α is chosen such that cos(α ₁) is close to 0, this ratio can be very large. For example by setting R ₂=49.55 m, we have α mentioned above equals to π, the ratio r=(1+2γ)/(δ · cos(α ₁))≈6700. This large ratio means that the retrieved A ² ₀ -map based on Eq. (11) has diminished image-contrast and erronous values, as is shown in Fig.3(c). The phase retrieval result based on Eq. (10), shown in Fig. 3(d), however is much more accurate. When R ₂=49.55 m, the FT-space Fresnel Propagator phase α=πλR ₂ u⃗²/M≈π. So we have, from Eq. (10), 𝓕̂((M ²/I _in)I)≈−𝓕̂(A ² ₀), which is different from the result 𝓕̂((M ²/I _in)I)≈[1+2π(λr_e/σ_KN)]𝓕̂(A ² ₀) derived from the conventional TIE prediction (Eq. (11)), which does not hold for such large Fresnel Propagator phases. This phenomenon is shown in Fig. (3)(b), where the intensity I is acquired with Fresnel diffraction, in which the dark-white area in A ² ₀ is reversed. Note that the inverse filter in Eq. (12) is singular at α ₁=(k−1/2)π, k=1,2,⋯, because of the zero-crossing of cos(α ₁). To deal with the singularity associated with the zero-crossing of cos(α ₁), one should regularize the inverse filter. A common regularization technique is the Tikhonov regularization, which modifies 1/cos(α ₁) in the filter of Eq. (12) to cos(α ₁)/(cos(α ₁)²+κ ²), where κ is the Tikhonov regularization parameter. In this specific example of the sinusoidal object, however, since we know 𝓕̂ (I) is dominated by frequency-components only at u_x=±u _x0 and (u,v)=(0,0), by setting 𝓕̂ (I)=0 at all other frequency points purposely, we can avoid the problem caused by the zero-crossing of the inverse filter.

 

Fig. 2. Plots of filters (13) and (12)

Download Full Size | PDF

 

Fig. 3. Comparation of the quality of the retrieved results A ² ₀ using Eq. (11) and Eq. (10). A ² ₀=exp[ϕ/γ], where phase shift ϕ=0.5sin(32px)-1, pixel pitch p=1µm, and object to detector distance R ₂=49.55m. (a) is the true A ² ₀ ; (b) is the phase-contrast intensity image I acquired from Fresnel diffraction; (c) is retrieved A ² ₀ using Eq. (11); (d) is retrieved A ² ₀ using Eq. (10).

Download Full Size | PDF

In another simulation, we set all the parameters the same as the previous simulation except replacing the test-object with the well-known Lena image as is shown in Fig. 4(a). The simulation results for the case of R ₂=5 m are shown in Figures 4(c)–4(d). The phase-contrast image I, shown in Fig. 4(b), is acquired with Fresnel diffraction from the object downstream to the detector [21]. Figure 4(c) shows the retrieved image with the inverse filter of Eq. (13), which is based on Eq. (11). Since the filter of Eq. (13) is incorrect when the Fresnel propagator phase α is big, say close or greater than π, the high frequency components of the object was unduly suppressed in the phase-retrieval. This deficiency is reflected in Fig. 4(c) such that the retrieved image gets blurred in locations with high contrast changes, such as the hair, the edge of the hat and body. In contrast, the phase retrieval using the inverse filter of Eq. (12), which is based on Eq. (10), preserves objects high frequency components, as is shown in Fig. 4(d).

As stated above, a regularization must be used to overcome the difficulty associated with the zero-crossing in the denominator of inverse filter in Eq. (13). The zero-crossing occurs obviously at α ₁=(k−1/2)π, k=1,2,⋯. The Tikhonov regularization would replace 1/cos(α ₁) in the inverse filter of Eq. (13) by the following

where κ is a small constant. However, this form of regularization is not optimal since it will deviate from the original inverse filter when α ₁ is small. In other words, the test-object’s low frequency components would not be fully reconstructed in the retrieval. In order to get around this drawback, we note that the zero-crossing occurs only when α1 is close or greater than π/2, while the most low frequency components have α<1. Hence we set the regularization parameter κ=0 when α is less than some threshold value α ₀, and κ>0 when α is greater than α ₀. In our simulation, we set this threshold α ₀=0.8, so κ=0 when α<0.8, and κ=0.12 when α>0.8. We can see from the image in Fig. 4(c) that most of the low frequency information is preserved. On the other hand, the quality of the results can also be numerically determined by computing the error comparing to the true A ² ₀ image. By taking the L ₂-measure, which is defined as the integrated squared errors, we have the absolute errors for methods (10) and (11) being 0.014 and 0.021 respectively.

 

Fig. 4. Fig. 4. Comparation of the retrieved results using Eq. (11) and Eq. (10) with R2=5m. (a) is the true A20; (b) is the phase-contrast image I acquired from Fresnel diffraction; (c) is retrieved using Eq. (11); (d) is retrieved using Eq. (10).

Download Full Size | PDF

4. Conclusions

Phase retrieval is required for achieving artifact-free x-ray phase-sensitive 3D imaging. Multiple image-acquisitions per projection are required in general for phase retrieval [3, 4, 15, 17, 22, 23, 24]. This requirement impedes the implement of x-ray phase-sensitive 3D imaging, especially for phase-sensitive 3D imaging in clinical applications. A phase-retrieval approach based on the phase-attenuation duality with high energy x-rays can greatly facilitate for phase sensitive imaging by allowing phase retrieval from only one single projection image. However, the previously derived phase retrieval formula is valid only for small Fresnel propagator phases corresponding to common clinical imaging tasks. In this work we generalized the previous phase retrieval formula to Eq. (10), which can be applied for cases with large Fresnel propagator phases corresponding to high spatial resolution imaging. The computer simulation demonstrated superiority of this new formula over the previous duality-based phase-retrieval formula in reconstructing the high frequency components of imaged objects. In addition, a modified Tikhonov regularization technique has been devised for phase retrieval in cases of very high resolution and large object-detector distance such that some Fresnel propagator phases may be close or greater than π. This new phase retrieval formula of Eq. (10) lays the foundation for implementing high-resolution phase-sensitive 3D imaging of soft tissue objects.