1. Introduction

X-ray phase-contrast imaging has emerged as a powerful imaging technique during the last two decades. Its strength lies mainly in the improved contrast for low-absorbing materials, such as soft biological tissue, compared to conventional absorption-based X-ray imaging [1]. The improvement stems from the fact that phase shift, which is the origin of phase contrast, can give a large contrast despite small differences in absorption. As the phase shift cannot be directly detected, the challenge of phase-contrast imaging is to obtain the phase through various indirect measures. Several methods have been developed [2], most notably grating-based imaging [3,4], propagation-based imaging (PBI) [5–7], and speckle-based imaging [8,9]. All three have been implemented both at synchrotron facilities and in compact laboratory arrangements.

An advantage of propagation-based phase-contrast imaging, compared to other phase-contrast techniques, is the simplicity of the laboratory arrangement. Furthermore, PBI poses a lower dose requirement compared to grating-based imaging when imaging smaller features using polychromatic, i.e., lab sources [10]. Provided partial spatial coherence and a longer propagation distance, the phase shift in the object appears as edge enhancement. As a first approximation, this edge-enhanced image is proportional to the transverse laplacian of the phase shift, not the phase shift itself [11]. Unprocessed images, with their strong edge enhancement, are ideal for feature detection, but further quantitative processing and analysis requires retrieval of the phase [12–14]. The phase retrieval aims to obtain the complex refractive index of the sample, n = 1 − δ + iβ. The two components of the complex refractive index, δ and β, represent phase shift and absorption respectively. Without any a priori information, one image is not sufficient to retrieve the complex refractive index. Multiple images, acquired at different propagation distances, provide sufficient information but pose a practical challenge, in particular for laboratory systems. For propagation-based imaging with compact lab sources, single-distance imaging is still the prevalent method. To circumvent the lack of information, a priori assumptions are used in the phase-retrieval process. For samples consisting of one or two materials, several such schemes exist [15]. However, most real samples consist of multiple materials, and are not properly retrieved using those techniques.

In this work, we compare the two existing analytical single-distance phase-retrieval methods designed for samples consisting of three or more materials. The two methods, by Beltran et al. [16, 17] and Ullherr & Zabler [18], are here referred to as the parallel method and the linear method, respectively. Through derivation of the phase-retrieval formula and quantitative simulations we show that the linear method returns quantitative values of the attenuation coefficient, μ, as already demonstrated for the parallel method [16]. In section 2 we discuss and clarify the derivation of the two methods, relying on the full derivation of the linear method as outlined in Appendix A. For completeness the derivation of the parallel method is also included, as previously outlined by Beltran et al. [16]. We show that the correction term for sample thickness used by Beltran et al. is not required in the linear method. Section 3 contains simulations to assess the methods qualitatively and quantitatively. Section 4 shows both methods used on experimental data from a laboratory arrangement. Section 5 highlights a previously neglected effect, inverted contrast, which appears for certain multi-material combinations and has important implications for phase retrieval. The discussion and conclusion in section 6 summarizes differences, characteristics, and challenges of the two methods.

2. Phase-retrieval formulas and algorithms

2.1. Phase-retrieval formulas

A widely used single-material phase-retrieval method is Paganin’s method [14]. It is derived from the transport of intensity equation (TIE) [19]. In the derivation of the TIE two assumptions on the free-space propagation are made. First, the Fresnel number F = a²/(dλ), where a is the smallest visible feature, d is the effective propagation distance [15] and λ is the wavelength, should fulfill the condition F ≫ 1. Equivalently, the smallest observable feature should be larger than $\sqrt{d\lambda }$. Second, the contrast in the intensity distribution must be low. The intensity is proportional to the transverse laplacian of the phase so the second assumption is equivalent to assuming a slowly varying phase. Paganin derives the formula [14]

The parallel method’s phase-retrieval formula for two materials [16] is an extension of Paganin’s formula:

For the linear method, the full derivation is outlined in Appendix A. The derivation is based on the assumption of slowly varying total projected thickness, resulting in the phase-retrieval formula [18]

In practice, Paganin’s methods is frequently used on multi-material samples, i.e., on samples not consisting of a single material in air or vacuum. To obtain a good phase retrieval despite this difference, the parameter δ/μ is tuned until edge enhancement is removed. The validity of this approach, though widely used, has not been mathematically confirmed. The optimal value of δ/μ at an interface between two material turns out to be Δδ/Δμ, i.e., the same as Eq. (3). We have derived this formula in Appendix A, and thereby shown that the commonly used approach of adjusting Eq. (1) actually has a solid mathematical foundation. Simply using equation Eq. (3), without the subsequent steps in the linear algorithm described in section 2.2, will be referred to as extended Paganin.

2.2. Algorithms

Figure 1 shows a flowchart of the parallel method’s algorithm, where numbers refer to the processing steps in the figure. It is assumed that A(r_⊥) is known. If not known a priori, it can be approximated from the data set using a conventional single-material phase retrieval (Paganin), tomographic reconstruction, and forward projection using the Radon transform [20]. When this pre-step is finished two or more phase retrievals, using different phase-retrieval parameters, are performed on the raw projection data. In general, the sample is assumed to consist of one or more materials inside an encasing material. The sample in Fig. 1 consists of two materials and air, requiring two separate phase retrievals.

Following the notation in Fig. 1, the parallel formula in Eq. (2), with the A(r_⊥)-factor, is used for all materials in the encasing (step 1). Δδ/Δμ is optimized to retrieve the interface between the encasing and the material in question. A normal Paganin phase retrieval is also applied to one copy of the raw images to retrieve the encasing (step 2). This can in theory be extended to an arbitrary number of materials. Next, tomographic reconstruction is performed on the phase-retrieved images. One reconstruction per sample material must be performed (steps 3 & 4). Finally, an intensity-based threshold segmentation (steps 5 & 6) is used to merge the data sets together (step 7). All materials are added to the encasing image. As these images contain only the difference between encased and encasing materials, μ₂ − μ₁, the constant value of μ₁ is added to the cut out G (see Fig. 1).

 

Fig. 1 The parallel algorithm: 1: The projection image (B) together with the projected thickness A(r_⊥) (A) are used in Beltran’s formula to retrieve (μ₂ − μ₁)T₂ (C). 2: D is also retrieved from B using Eq. (1). 3 & 4: C and D are tomographically reconstructed. Slices are shown in E and F, respectively. 5 & 6: The segmentation is performed with an intensity threshold in F. This is used to obtain only material 2 from E, and to remove material 2 from F. 7: The cut out G is combined with H. The value of the encasing (μ₁) must be added to G before merging to get μ in the final image I.

Download Full Size | PDF

The parallel method’s final steps can be expressed as

 

Fig. 2 The linear method’s algorithm: 1: phase retrieval on projection images (A), 2: tomographic reconstruction, 3: segmentation (highly absorbing part D removed), 4: additional phase retrieval on volume image (E), 5: merge of D and F to form final result G.

Download Full Size | PDF

The linear method outlined in Fig. 2 is constructed for one highly absorbing material, such as bone, in the same sample as several less-absorbing materials, such as different kinds of tissue or air. Here the so-called forward correction is used. The linear method also exists as a backward correction, though not treated here due to its iterative nature and less successful result [18]. The raw projection images are immediately phase retrieved (step 1). This phase retrieval is optimized for the interface between highly and less absorbing materials. The second step is the tomographic reconstruction (step 2). In the volume image an intensity based threshold segmentation is employed to cut out the highly absorbing part (step 3). The less absorbing part will have remaining edge enhancement as the first phase retrieval was incomplete for this interface. A new phase retrieval is therefore applied to this part (4). This phase retrieval is carried out in the volume image, not on a projection image. To keep the values correct after the second phase retrieval the image’s segmentation mask is filtered like the image and then used to normalize the image. After the volume phase retrieval the images are merged (step 5). Mathematically, the volume image manipulations (step 3–5 in Fig. 2) can be expressed as

The term Δδ₁/Δμ₁ in the nominator is the same as the one used in the first phase retrieval, e.g., the interface between tissue and bone. The ratio in the denominator Δδ₂/Δμ₂ is calculated from δ and μ for the less absorbing materials, e.g., air and soft tissue. The K_L filter will in most cases be between zero and one (for exceptions see section 5), since the remaining edge enhancement should be removed.

2.3. Validity of phase retrieval on volume images

In the linear algorithm phase retrieval is performed both before and after the tomographic reconstruction. It is not immediately obvious that phase retrieval can be performed in the volume image, at least with the currently used formula. It is, however, possible because both phase retrieval and normal tomographic reconstruction (filtered back projection) are filters in Fourier space. The commutativity of two consecutive Fourier filters makes it possible to change the order. Ruhlandt et al. have shown this formally [21]. The phase retrieval formulas contain, however, a logarithm. In theory this makes it impossible to change the order, but in practice the input in phase retrieval, I/I₀, is around 1 due to the low contrast. Hence the logarithm is approximately linear. This corresponds to the condition of an optically weak object, very often the case in phase-contrast imaging.

3. Quantitative simulations

To quantitatively test the methods, a simple phantom containing three materials was simulated for monochromatic radiation (20 keV) using the simulation scheme described by Lundström [22]. The simulated system source spot and detector PSF had a FWHM of 10 μm and 25 μm respectively. The simulated images were sampled to model a detector with 9 μm pixels and a realistic noise power spectrum. The simulation sampling was 0.9 μm to avoid aliasing between pixels. The source-to-object distance was 1 m and the source-to-detector distance was 2.5 m giving an effective propagation distance of 0.6 m. The phantom was a 1 mm thick cylinder, in air, with a 0.4 mm core. The outer layer was adipose tissue and the core soft tissue. Though it might appear trivial, this kind of sample cannot be well retrieved using single-material phase retrieval. Phase retrieval using both the parallel and the linear method was applied to the simulated data. The theoretical values of Δδ/Δμ for the two interfaces were used in the reconstruction. Repeated trials to observe which Δδ/Δμ removes the edge enhancement can be used if the materials are unknown, albeit unpractical for large data sets. The theoretical values used in the simulation and reconstruction were, in this case, ratios Δδ₁/Δμ₁ = 1.7 · 10⁻⁹ m and Δδ₂/Δμ₂ = 1 · 10⁻⁸ m.

Phase retrieval using the linear method is shown in Fig. 3, on noise-free data. Intensity threshold segmentation cannot be directly applied to image 3(b) as the remaining edge enhancement is strong, as seen in the corresponding profile in Fig. 3. The problem is solved by creating another mask. By performing a strong phase retrieval on image 3(b) (ratio 1 · 10⁻⁸ m) followed by a threshold segmentation (threshold = 56 m⁻¹) a rough mask is obtained. This mask does not accurately follow the highly-less absorbing interface, but it singles out the right region. If a normal mask (ratio 1.7 · 10⁻⁹ m, threshold 60 m⁻¹) is multiplied with this rough mask all regions with remaining edge enhancement are removed since they are not part of the rough mask. The simulations, as illustrated in Fig. 4, are done both with and without noise. The noise level is estimated in the raw projection images by the differential signal-to-noise ratio [23, p. 169]

 

Fig. 3 Phase retrieval using the linear method on simulated cylindrical object. (a) Tomographic slice without phase retrieval. (b) Slice after partial phase retrieval. (c) Slice after second phase retrieval and merge. (d–f) Profiles of the images above.

Download Full Size | PDF

The simulation shows that the linear method successfully retrieves μ quantitatively without leaving edge enhancement or causing any substantial blurring. Values for the different materials and noise levels are shown in Table 1. The error is less than 4% for the linear method.

Table 1. Agreement between simulated and phase retrieved μ-values [m⁻¹]. All values are given as mean ± standard deviation. Differential signal-to-noise ratio (SNR_diff) for simulated projections is given.

View Table | View all tables in this article

The parallel method is quantitative as well with an error of less than 4%. One clear difference is, however, visible in the result (see Fig. 5). The material-material edge in the linear method is sharper. This is because the more absorbing part (here: soft tissue) is removed before the phase retrieval in the volume image is applied. In the parallel method, all phase retrieval is performed on the projection images and the highly-absorbing material will therefore affect the surrounding region. We can quantify this difference in sharpness by defining a measure similar to the rise time used in signal processing. The rise distance for an edge is the number of pixels having values between 10% and 90% of the final value. The image retrieved with the linear method has a rise distance that is comparable to the simulated system resolution, whereas the rise distance for the parallel method is 3–4 times longer. The parallel method can be improved by adjusting the segmentation. Instead of setting the mask boundary at the sharp transition between soft tissue and adipose tissue it can be moved further away. This will reduce the negative effect but the solution has limitations. If there is an edge in close proximity to the sharp transition, as can be seen in Fig. 5, there is a strict limitation on how far the segmentation can be adjusted.

The quality of the reconstruction can also be quantified with a χ² test that measures the difference between the object and the reconstruction: ${\chi }^2=\sum {\left|\text{reconstructed}-\text{object}\right|}^2/\sum \text{object}$. A perfect reconstruction would result in a χ² of zero [24]. We note that the absolute results of χ² are not relevant in this case, but rather the relative differences between reconstructions. Values for the multi-materials methods are shown in Table 2. Without additional noise Paganin’s method has a χ² of 126. The result follows previous observations in that the multi-material methods show a clear improvement over single-material methods, and the linear method produces slightly more accurate results.

Table 2. χ² values for the multi-material reconstructions with different levels of noise.

View Table | View all tables in this article

In the presence of noise there is a clear reduction in accuracy (see Fig. 4 and Table 1 & 2). The mean values decrease and the standard deviation increases. The linear method displays a smaller standard deviation at high noise. The large drop in μ that the parallel method displays for soft tissue is partly due to the blurry edge.

Repeating the simulation and phase retrieval for polychromatic radiation does not result in a quantitative reconstruction, nor correct proportions between reconstructed values. For single-material samples there have been extensions to conventional phase-retrieval methods that allow quantitative reconstruction [25,26]. However, these extensions do not solve the problem in multi-material phase retrieval. Using the full source spectrum could potentially provide quantitative reconstruction.

 

Fig. 4 Simulated slices after the linear method’s partial phase retrieval showing three noise levels: (a) no added noise, (b) SNR_diff = 1.4, and (c) SNR_diff = 0.9.

Download Full Size | PDF

 

Fig. 5 Phase retrieval using the parallel method on simulated cylindrical object. (a) Tomographic slice after phase retrieval with the parallel formula. (b) Slice after phase retrieval with Paganin’s formula. (c) Merge of (a) and (b) according to the parallel algorithm. (d–f) Profiles of the images above.

Download Full Size | PDF

4. Phase retrieval on experimental data

We demonstrate the phase-retrieval methods on a tomography of an ex vivo human coronary artery. The experimental data was acquired in a previous work [27]. This coronary artery suffers from atherosclerosis, the cause of ischemic heart disease which is the leading cause of death in the world [28]. In atherosclerosis, depositions of calcium and/or lipids are accumulated within the artery wall. Discriminating the less absorbing lipid depositions and the artery wall smooth muscle tissue in the presence of highly absorbing calcifications may provide further understanding of the biological mechanisms. Fig. 6 shows a tomographic slice where phase retrieval has been applied using different methods. The regions visible in the images are calcified plaque (A), the air-filled lumen (B), adipose tissue (C), and muscle tissue (D). In (a) the plaque-muscle tissue interface AD was retrieved with δ/μ = 2 · 10⁻¹⁰ m using the extended Paganin method. It is clear that, while showing all features well, the image has an incomplete phase retrieval with remaining edge enhancement at the tissue-air interface. In (b) the muscle tissue-air interface BD is retrieved with δ/μ = 1 · 10⁻⁸ m using the extended Paganin method. The result is no remaining edge enhancement, but the image suffers from blurring around the highly absorbing plaque. Applying the linear and the parallel methods in (c) and (d), respectively, we get a different result. The blurring is reduced while edge enhancement is removed. We note that the interface AD in (d) displays a slight blurring. This is the material-material interface blurring in the parallel method described in section 3.

 

Fig. 6 Phase retrieved images of coronary artery with (a) partial and (b) full extended Paganin phase retrieval as well as (c) the linear method, and (d) the parallel method.

Download Full Size | PDF

To quantify the difference between the phase retrievals we have defined two simple metrics: overshoot and rise distance. The overshoot estimates the remaining edge enhancement, and the rise distance (as defined in Section 3) estimates the interface blurring. A low value is desirable for both metrics. Measurements were done for the BD and AD interfaces, as marked in Fig. 6(a) by a white line and a black line, respectively. The results are shown in Table 3. As expected the partial retrieval has a clear overshoot due to remaining edge enhancement, while the full phase retrieval shows no overshoot. The multi-material methods have both some overshoot for interface AD. Further phase retrieval could remove the overshoot, but it would also increase the rise distance (blurring). The rise distance for the AD interface is clearly reduced using the multi-material methods. The interface BD is retrieved last and is therefore not affected much by the multi-material methods.

Table 3. Overshoot and rise distance at interfaces BD and AD for all images in Fig. 6. An advantage of multi-material methods as well as clear trade-off between overshoot and rise distance is evident from the numbers.

View Table | View all tables in this article

5. Low-frequency contrast reversal due to negative Δδ/Δμ - ratio

The multi-material methods, as well as the widely used extension of the Paganin method, share a problem, that in certain cases complicates the phase retrieval. In Beltran et al.’s derivation the following expression appears (see Eq. (24), Appendix A),

The right-hand side is known from experimental data and the left-hand side is unknown. Dividing by $\left(d({\delta }_2-{\delta }_1)/({\mu }_2-{\mu }_1)4{\pi }^2{\mathbf{w}}_{\perp }^2+1\right)$ ultimately yields the sought parameter. The same step exists in the linear method (see Eq. (36), Appendix A). In most cases Δδ/Δμ is positive and this operation is trivial. The term Δδ/Δμ can, however, be negative. If so, at the spatial frequency

the Paganin method implies division by zero. This occurs at relatively low spatial frequencies. The zero implies a loss of information not during the phase retrieval, but during the imaging process since the transfer function is zero. This issue can not be resolved simply by taking the absolute value of Δδ/Δμ. Instead, the negative value should be used, and a Tikhonov regularization implemented to limit the noise amplification near the zero division [29,30]. This type of modification, shown in Fig. 7, is also referred to as a Wiener filter [31].

 

Fig. 7 Effect of a Wiener filter on a phase retrieval filter H as a function of frequency ω. The dashed dotted line in blue is a conventional filter obtained for positive Δδ/Δμ. The solid red line shows a filter with zeros in the denominator caused by negative Δδ/Δμ and the dashed black line the result of applying a Wiener filter to the red filter.

Download Full Size | PDF

If a sinusoidal pattern, with period corresponding to the lost frequency, is imaged the contrast becomes zero (see Fig. 8). For frequencies above this zero frequency the contrast is reversed. The zero contrast occurs when the intensity contributions from phase and attenuation are equally strong and shifted by π. A similar interference effect has been shown in [32]. The image in Fig. 8(c) has been retrieved with a Wiener filter and the reversion of contrast has been removed. The region with zero contrast is however largely unchanged.

The implication of a negative Δδ/Δμ ratio is clearly visible in a tomographic image. Figure 9 shows the slice of a simulated cylindrical object. The cylinder is Polyvinyl chloride (PVC) and the surrounding material is Polytetrafluoroethylene (PTFE). For most material combinations the phase contrast adds constructively to the absorption contrast to produce a positive edge enhancement. For the PVC-PTFE interface this is not the case. Here the edge is reversed and the phase contrast produces a negative edge enhancement (see Fig. 9(a)). The result of a phase retrieval using a positive Δδ/Δμ is shown in Fig. 9(b). Despite the improved contrast, a new kind of edge enhancement has appeared. If the sign is correct, and a Wiener filter is applied Fig. 9(c) is obtained. The edge is retrieved quite well, but low frequencies close to w₀ are amplified, thus degrading the result.

 

Fig. 8 Sinusoidal Siemens star with negative Δδ/Δμ material combination. The diameter in the detector plane is 6.3 mm. (a) Object, (b) simulated phase contrast image, and (c) phase retrieved image.

Download Full Size | PDF

As outlined so far, this phenomenon affects phase retrieval of single-interface samples. For phase retrieval of several different interfaces, such as the methods tested in this paper, there is an additional problem; samples can now contain two different interfaces, one that should be retrieved with a positive Δδ/Δμ and one with a negative. Performing either will negatively affect the other. This will not be solved with a simple Wiener filter.

 

Fig. 9 (a) Tomographic slice of a sample containing two materials. The interface has negative edge enhancement (compare with Fig. 3). (b) Result of erroneous phase retrieval with positive Δδ/Δμ if applied to image (a). (c) Result of phase retrieval with negative Δδ/Δμ and a Wiener filter. Although the retrieval of the edge is improved, amplification of low frequencies close to w₀ are visible.

Download Full Size | PDF

6. Discussion and conclusion

From both simulations and experimental data it is clear that the linear and parallel methods are comparable regarding image quality. Both accurately produce quantitative results for monochromatic radiation, but not for polychromatic radiation, at least not without additional information and processing. An advantage of the linear method is that the parameters may be set one at a time whereas in the parallel method both the Δδ/Δμ ratios and μ of the encasing must be known. The linear method also removes the highly absorbing material before the phase retrieval is complete which improves the final image by reducing material-material interface blurring.

Considering the potential data size, computational performance is of importance. Both methods perform multiple phase retrievals increasing the computational time by a factor of three or four compared to single-material algorithms such as Paganin’s method. However, the computational complexity, i.e. how the computational time scales with data size, is not changed. The major difference is in the tomographic reconstruction. The linear method, like any other analytical method, requires a single tomographic reconstruction. The parallel method, however, requires several such reconstructions. Again there is no change in terms of complexity, but the increase in time can be substantial. Unless the total projected thickness is known through laser profilometry or some other method, it must be calculated, adding not only another tomographic reconstruction, but also a forward projection. We note that if reconstruction is done in a dedicated software, forward projection may not be available. Performing the forward projection, care must be taken in alignment and resizing so that the calculated projected thickness images correspond well to an original projection image. All this makes the parallel method more time-consuming and also impractical to work with.

Another critical issue is the segmentation. High noise results in faulty segmentation. Segmentation is also made more complicated as remaining edge enhancement can have as high pixel values as the highly absorbing materials. This second issue is particularly relevant for the linear method. In the parallel method segmentation usually also has to be applied in the pre-step, where the total projected thickness is calculated. This segmentation is done to differentiate between air and sample, which can be difficult, partially due to artifacts.

We have shown that for certain material combinations and energies the edge enhancement is reversed. For single-interface samples a Wiener filter partly solves this issue, but for multi-material samples this approach will not work. However, as δ often increases with μ it is rare to come across a material combination that has Δδ/Δμ < 0 in the normal energy range. This phase retrieval problem should therefore be a minor issue in biomedical applications, while for imaging of, e.g., composite plastics it might pose a problem.

In conclusion, we have shown theoretically and by simulations that the linear method is derived from the TIE and can return quantitative images. We have also shown that multi-material samples can exhibit contrast inversion. Both compared methods produce quantitatively correct images of the samples given monochromatic radiation and low noise. The linear method has a small advantage in quality, and a larger advantage in computational time and implementation.

Appendix A: Derivation phase-retrieval formulas

Below follows the mathematical derivations of the phase retrieval methods in this paper. Both methods follow essentially Paganin’s derivation from 2002 [14], but they handle the issue of a multi-material samples differently. The derivation of the parallel method is outlined in [16], while the full derivation of the linear method has not been previously published.

Derivation of the parallel phase-retrieval formula

The phase is related to the measured intensity by the Transport of intensity equation (TIE) [19]

(10)$$\nabla \cdot \left[I(\mathbf{r},z)\nabla \varphi (\mathbf{r},z)\right]=-k{\partial }_{z}I(\mathbf{r},z),$$

where I is the intensity, φ is the phase shift, k is the wave number, z is the distance in the propagation direction and r is the spatial coordinate in the plane perpendicular to z. The transmitted intensity is given by the Beer–Lambert law for a sample with two materials,

(11)$$I(\mathbf{r},z=0)=I_0\text{exp}\left(-{\mu }_1{T}_1(\mathbf{r})-{\mu }_2{T}_2(\mathbf{r})\right),$$

where I₀ is the incident intensity, μ is the attenuation coefficient, and T(r) is the projected thickness of the respective materials. The phase shift depends on the thickness and real decrement δ of the refractive index n = 1− δ + iβ as

(12)$$\phi (\mathbf{r},z=0)=-k\left[{\delta }_1{T}_1(\mathbf{r})+{\delta }_2{T}_2(\mathbf{r})\right].$$

The total projected thickness is defined as

(13)$$A(\mathbf{r})=T_1(\mathbf{r})+T_2(\mathbf{r}).$$

Substituting Eqs. (11) and (12) into the TIE yields

(14)$$\nabla \cdot \left[I_0\text{exp}\left(-{\mu }_1{T}_1(\mathbf{r})-{\mu }_2{T}_2(\mathbf{r})\right)\nabla \left[-k\left({\delta }_1{T}_1(\mathbf{r})+{\delta }_2{T}_2(\mathbf{r})\right)\right]\right]=-k{\partial }_{z}I(\mathbf{r},z),$$

and removing the factor −k gives

(15)$$\nabla \cdot \left[I_0\text{exp}\left(-{\mu }_1{T}_1(\mathbf{r})-{\mu }_2{T}_2(\mathbf{r})\right)\nabla \left[\left({\delta }_1{T}_1(\mathbf{r})+{\delta }_2{T}_2(\mathbf{r})\right)\right]\right]={\partial }_{z}I(\mathbf{r},z).$$

The two parts containing μ, δ and the thicknesses are rewritten as

(16)$${\delta }_1{T}_1+{\delta }_2{T}_2={\delta }_1(T_1+T_2)+({\delta }_2-{\delta }_1)T_2={\delta }_{1}A+({\delta }_2-{\delta }_1)T_2$$

(17)$$-{\mu }_1{T}_1-{\mu }_2{T}_2=-{\mu }_1(T_1+T_2)-({\mu }_2-{\mu }_1)T_2=-{\mu }_{1}A-({\mu }_2-{\mu }_1)T_2.$$

Substituting Eqs. (16) and (17) into Eq. (15) yields

(18)$$\nabla \cdot \left[I_0\text{exp}\left(-{\mu }_{1}A(\mathbf{r})\right)\text{exp}\left(-\left({\mu }_2-{\mu }_1\right)T_2(\mathbf{r})\right)\nabla \left[{\delta }_{1}A(\mathbf{r})+({\delta }_2-{\delta }_1)T_2(\mathbf{r})\right]\right]={\partial }_{z}I(\mathbf{r},z).$$

The projected thickness of the entire object is assumed to vary slowly, i.e. ∇A(r) ≈ 0, yielding

(19)$$I_0\text{exp}\left(-{\mu }_{1}A(\mathbf{r})\right)\nabla \cdot \left[\text{exp}\left(-\left({\mu }_2-{\mu }_1\right)T_2(\mathbf{r})\right)\nabla \left[({\delta }_2-{\delta }_1)T_2(\mathbf{r})\right]\right]={\partial }_{z}I(\mathbf{r},z).$$

Using the chain rule, the left-hand side of Eq. (19) may be written

(20)$$({\delta }_2-{\delta }_1)I_0\text{exp}\left(-{\mu }_{1}A(\mathbf{r})\right)\left(-\frac{1}{{\mu }_2-{\mu }_1}{\nabla }^2\text{exp}\left(-(-{\mu }_2-{\mu }_1)T_2(\mathbf{r})\right)\right).$$

Consider only the right-hand side of Eq. (20). The approximate derivative may be written as

(21)$$\begin{array}{ll}\frac{\partial }{\partial z}I(\mathbf{r},z)\hfill & \approx \frac{I(\mathbf{r},z=d)-I(\mathbf{r},z=0)}{d}\hfill \\ \hfill & =\frac{I(\mathbf{r},z=d)-I_0\text{exp}\left(-({\mu }_2-{\mu }_1)T_2(\mathbf{r})\right)\text{exp}\left(-{\mu }_{1}A(\mathbf{r})\right)}{d}.\hfill \end{array}$$

The full equation is then

(22)$$\begin{array}{c}({\delta }_2-{\delta }_1)I_0\text{exp}(-{\mu }_{1}A(\mathbf{r}))\left(-\frac{1}{{\mu }_2-{\mu }_1}{\nabla }^2\text{exp}\left(-({\mu }_2-{\mu }_1)T_2(\mathbf{r})\right)\right)\\ =\frac{I(\mathbf{r},z=d)-I_0\text{exp}\left(-({\mu }_2-{\mu }_1)T_2(\mathbf{r})\right)\text{exp}(-{\mu }_{1}A(\mathbf{r}))}{d}.\end{array}$$

Rearranging yields

(23)$$\left[-\frac{d({\delta }_2-{\delta }_1)}{({\mu }_2-{\mu }_1)}{\nabla }^2+1\right]\text{exp}\left(-({\mu }_2-{\mu }_1)T_2(\mathbf{r})\right)=\frac{I(\mathbf{r},z=d)}{I_0\text{exp}\left(-{\mu }_{1}A(\mathbf{r})\right)}.$$

Taking the Fourier transform of Eq. (23), the Laplacian operator, acting on the space coordinate r_⊥, produces a factor $-4{\pi }^2{\mathbf{w}}_{\perp }^2$. The equation can be writen as

(24)$$𝔉\left\{\text{exp}\left(-({\mu }_2-{\mu }_1)T_2(\mathbf{r})\right)\right\}\left(\frac{d({\delta }_2-{\delta }_1)}{({\mu }_2-{\mu }_1)}4{\pi }^2{\mathbf{w}}_{\perp }^2+1\right)=𝔉\left\{\frac{I(\mathbf{r},z=d)}{I_0\text{exp}\left(-{\mu }_{1}A(\mathbf{r})\right)}\right\}.$$

Rearranging and applying the inverse Fourier transform and logarithm gives the final result as

(25)$$({\mu }_2-{\mu }_1)T_2(\mathbf{r})=-\text{ln}\left({𝔉}^{-1}\left\{\frac{1}{\frac{d({\delta }_2-{\delta }_1)}{({\mu }_2-{\mu }_1)}4{\pi }^2{\mathbf{w}}_{\perp }^2+1}𝔉\left\{\frac{I(\mathbf{r},z=d)}{I_0\text{exp}\left(-{\mu }_{1}A(\mathbf{r})\right)}\right\}\right\}\right).$$

Derivation of the linear phase-retrieval formula

Again we start from Beer-Lambert’s law, though slightly rewritten:

(26)$$I(\mathbf{r},z=0)=I_0\text{exp}\left(-{\mu }_2{T}_2(\mathbf{r})-{\mu }_1{T}_1(\mathbf{r})\right)=I_0\text{exp}\left(-\int \mu (\mathbf{r},z)\text{d}z\right).$$

The phase depends upon thickness and real increment δ as

(27)$$\phi (\mathbf{r},z=0)=-k\left[{\delta }_2{T}_2(\mathbf{r})+{\delta }_1{T}_1(\mathbf{r})\right].$$

The entity

(28)$$M(\mathbf{r})={\mu }_2{T}_2(\mathbf{r})+{\mu }_1{T}_1(\mathbf{r}),$$

representing the total attenuation of the sample is defined. Looking at two pairs of δ and μ in a 2-dimensional δ − μ space, as illustrated in Fig. 10, one can draw a line between the pair and write an expression for the line as

(29)$${\delta }_1={\delta }_0+\frac{\mathrm{\Delta }\delta }{\mathrm{\Delta }\mu }{\mu }_1$$

where Δδ/Δμ is the slope and δ₀ is the offset. Inserting the expressions in Eq. (29) into Eq. (27) yields

(30)$$\frac{\phi (\mathbf{r},z=0)}{-k}={\delta }_2{T}_2(\mathbf{r})+{\delta }_1{T}_1(\mathbf{r})=\left({\delta }_0-\frac{\mathrm{\Delta }\delta }{\mathrm{\Delta }\mu }{\mu }_2\right)T_2(\mathbf{r})+\left({\delta }_0-\frac{\mathrm{\Delta }\delta }{\mathrm{\Delta }\mu }{\mu }_1\right)T_1(\mathbf{r})$$

(31)$$={\delta }_0\left(T_2(\mathbf{r})+T_1(\mathbf{r})\right)+\frac{\mathrm{\Delta }\delta }{\mathrm{\Delta }\mu }M(\mathbf{r})={\delta }_{0}A(\mathbf{r})+\frac{\mathrm{\Delta }\delta }{\mathrm{\Delta }\mu }M(\mathbf{r})$$

Insertion into the TIE yields

(32)$$\nabla \cdot \left[I_0\text{exp}(-M(\mathbf{r}))\nabla \left[-k\left({\delta }_{0}A(\mathbf{r})+\frac{\mathrm{\Delta }\delta }{\mathrm{\Delta }\mu }M(\mathbf{r})\right)\right]\right]=-k{\partial }_{z}I(\mathbf{r},z).$$

Removing the factor −k and using ∇A(r) ≈ 0 due to slow variations of the total thickness gives

(33)$$\frac{\mathrm{\Delta }\delta }{\mathrm{\Delta }\mu }\nabla \cdot \left[I_0\text{exp}\left(-M(\mathbf{r})\right)\nabla M(\nabla )\right]={\partial }_{z}I(\mathbf{r},z).$$

Use the same steps as Paganin [14] the equation may be rewritten as

(34)$$-\frac{\mathrm{\Delta }\delta }{\mathrm{\Delta }\mu }{I}_0{\nabla }^2{\text{e}}^{-M(\mathbf{r})}=\frac{I(\mathbf{r},z=d)-I_0\text{exp}\left(-M(\mathbf{r})\right)}{d}$$

or

(35)$$\left(-\frac{\mathrm{\Delta }\delta }{\mathrm{\Delta }\mu }d{\nabla }^2+1\right)I_0\text{exp}\left(-M(\mathbf{r})\right)=I(\mathbf{r},z=d)$$

In a manner analogous to the step between Eqs. (23) and (24), performing a Fourier transform yields

(36)$$𝔉\left\{\text{exp}\left(-M(\mathbf{r})\right)\right\}\left(\frac{\mathrm{\Delta }\delta }{\mathrm{\Delta }\mu }4{\pi }^{2}d{\mathbf{w}}_{\perp }^2+1\right)=𝔉\left\{\frac{I(\mathbf{r},z=d)}{I_0}\right\}$$

Finally, applying an inverse Fourier transform and using the definition of M(r) in Eq. (28) yields the phase-retrieval formula used in the linear phase retrieval algorithm:

(37)$$\int \mu (\mathbf{r},z)\text{d}z=M(\mathbf{r})=-\text{ln}\left({𝔉}^{-1}\left\{\frac{1}{\frac{\mathrm{\Delta }\delta }{\mathrm{\Delta }\mu }4{\pi }^{2}d{\mathbf{w}}_{\perp }^2+1}𝔉\left\{\frac{I(\mathbf{r},z=d)}{I_0}\right\}\right\}\right)$$

 

Fig. 10 For two materials there is a linear relationship between δ and μ at a given energy.

Download Full Size | PDF

Funding

Knut and Alice Wallenberg Foundation; Swedish Research Council.

Acknowledgments

We thank Maximilian Ullherr for sharing the linear method code which we used as a starting point to develop our own.

Disclosures

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

References and links

1. R. Fitzgerald, “Phase-sensitive X-ray imaging,” Phys. Today 53(7), 23–26 (2000). [CrossRef]  

2. A. Momose, “Recent Advances in X-ray Phase Imaging,” Jpn. J. Appl. Phys. 44, 6355–6367 (2005). [CrossRef]  

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

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

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

6. P. Cloetens, R. Barrett, J. Baruchel, J-P. Guigay, and M. Schlenker, “Phase objects in synchrotron radiation hard x-ray imaging,” J. Phys. D 29, 133–146 (1996). [CrossRef]  

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

8. K. S. Morgan, D. M. Paganin, and K. K. W. Siu, “X-ray phase imaging with a paper analyzer,” Appl. Phys. Lett. 100(12), 124102 (2012). [CrossRef]  

9. I. Zanette, T. Zhou, A. Burvall, U. Lundström, D. H. Larsson, M.-C. Zdora, P. Thibault, F. Pfeiffer, and H. M. Hertz, “Speckle-based x-ray phase-contrast and dark-field imaging with a laboratory source,” Phys. Rev. Lett. 112, 253903 (2014). [CrossRef]   [PubMed]  

10. T. Zhou, U. Lundström, T. Thüring, S. Rutishauser, D. H. Larsson, M. Stampanoni, C. David, H. M. Hertz, and A. Burvall, “Comparison of two x-ray phase-contrast imaging methods with a microfocus source,” Opt. Express 21(25), 30183–30195 (2013). [CrossRef]  

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

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

13. A. V. Bronnikov, “Reconstruction formulas for phase-contrast imaging,” Opt. Commun. 171, 239–244 (1999). [CrossRef]  

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

15. 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–10376 (2011). [CrossRef]   [PubMed]  

16. 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–6436 (2010). [CrossRef]   [PubMed]  

17. M. Beltran, D. Paganin, K. Siu, A. Fouras, S. Hooper, D. Reser, and M. Kitchen, “Interface-specific x-ray phase retrieval tomography of complex biological organs,” Phys. Med. Biol. 56, 7353–7369 (2011). [CrossRef]   [PubMed]  

18. M. Ullherr and S. Zabler, “Correcting multi material artifacts from single material phase retrieved holo-tomograms with a simple 3D Fourier method,” Opt. Express 23(25), 32718–32727 (2015). [CrossRef]   [PubMed]  

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

20. J. Radon, “Über die Bestimmung von Funktionen durch ihre Integralwerte längs gewisser Mannigfaltigkeiten,” Berichte über die Verhandlungen der Königlich-Sächsischen Akademie der Wissenschaften zu Leipzig 69, 262 (1917).

21. A. Ruhlandt and T. Salditt, “Three-dimensional propagation in near-field tomographic X-ray phase retrieval,” Acta Crystallographica Section A 72(2), 215–221 (2016). [CrossRef]  

22. U. Lundström, Phase-Contrast X-ray Carbon Dioxide Angiography, Ph.D. thesis (KTH Royal Institute of Technology, 2014).

23. J. Beutel, H. L. Kundel, and R. L. Van Metter, Handbook of Medical Imaging: Volume 1. Physics and Psychophysics (SPIE Press, 2000), Chap. 3.

24. J. L. Devore and K. N. Berk, Modern Mathematical Statistics with Applications, 2nd ed. (Springer, 2012), Chap. 13. [CrossRef]  

25. G. R. Myers, S. C. Mayo, T. E. Gureyev, D. M. Paganin, and S. W. Wilkins, “Polychromatic cone-beam phase-contrast tomography,” Phys. Rev. A 76, 045804 (2007). [CrossRef]  

26. B. D. Arhatari, K. Hannah, E. Balaur, and A. G. Peele, “Phase Imaging Using A Polychromatic X-ray Laboratory Source,” Opt. Express 16(24), 19950–19956 (2008). [CrossRef]   [PubMed]  

27. W. Vågberg, J. Persson, L. Szekely, and H. M. Hertz, “Cellular-resolution 3D virtual histology of human coronary arteries using x-ray phase tomography,” Submitted.

28. World Health Organization, Health in 2015: from MDGs, Millennium Development Goals to SDGs, Sustainable Development Goals (World Health OrganizationFrance, 2015).

29. T. E. Gureyev, T. J. Davids, A. Pogany, S. C. Mayo, and S. W. Wilkins, “Optical phase retrieval by use of first Born-and Rytov-type approximations,” Appl. Opt. 43(12), 2418–2430 (2004). [CrossRef]   [PubMed]  

30. A. N. Tikhonov and V. Y. Arsenin, Solutions of Ill-Posed Problems (Wiley, 1977).

31. J. W. Goodman, Introduction to Fourier optics, 3rd ed. (Roberts & Company Publishers, 2005), Chap. 8.

32. S. Zabler, H. Riesemeier, P. Fratzl, and P. Zaslansky, “Fresnel-propagated imaging for the study of human tooth dentin by partially coherent x-ray tomography,” Opt. Express 14(19), 8584–8597 (2006). [CrossRef]   [PubMed]