## Abstract

Propagation-based phase-contrast X-ray imaging provides high-resolution, dose-efficient images of biological materials. A crucial challenge is quantitative reconstruction, referred to as phase retrieval, of multi-material samples from single-distance, and hence incomplete, data. In this work, the two most promising methods for multi-material samples, the parallel method, and the linear method, are analytically, numerically, and experimentally compared. Both methods are designed for computed tomography, as they rely on segmentation in the tomographic reconstruction. The methods are found to result in comparable image quality, but the linear method provides faster reconstruction. In addition, as already done for the parallel method, we show that the linear method provides quantitative reconstruction for monochromatic radiation.

© 2017 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

## 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*^{2}/(*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]

*T*is the projected thickness of the object,

*μ*is the attenuation coefficient of the object material,

*d*is the effective propagation distance,

*I*is the intensity in the phase-contrast image,

*I*

_{0}is a flat-field image, and

**w**are the spatial frequencies corresponding to

**r**

_{⊥}, the set of spatial coordinates. 𝔉

_{2D}and ${\U0001d509}_{2D}^{-1}$ are the two-dimensional Fourier and inverse Fourier transforms, respectively. Here it is assumed that the sample consists of one material surrounded by air/vacuum, and that

*δ*and

*β*are known. The formulas used for the parallel method by Beltran

*et al.*[16] and for the linear method by Ullherr & Zabler [18] both reduce to Eq. (1) for samples containing one material and air/vacuum.

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

*μ*

_{1}and

*δ*

_{1}, and

*μ*

_{2}and

*δ*

_{2}respectively, and their projected thicknesses are

*T*

_{1}(

**r**

_{⊥}) and

*T*

_{2}(

**r**

_{⊥}). To arrive at this result Beltran

*et al.*[16] assumes that the total sample thickness varies slowly. The main difference from Eq. (1) is that the

*δ/μ*ratio is replaced by Δ

*δ*/Δ

*μ*= (

*δ*

_{2}−

*δ*

_{1})/(

*μ*

_{2}−

*μ*

_{1}), i.e. the ratio between the difference in

*δ*and

*μ*for two materials. A factor exp (−

*μ*

_{1}

*A*(

**r**

_{⊥})) is also included in the initial normalization step.

*A*(

**r**

_{⊥}) =

*T*

_{1}(

**r**

_{⊥}) +

*T*

_{2}(

**r**

_{⊥}) is the total projected thickness of the sample, that is, the thickness of both materials minus any internal voids.

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]

*μ*(

*z*) now varies with position. However, it differs from the quantity retrieved in Eq. (2). In the parallel method, the retrieved quantity is the difference (

*μ*

_{2}−

*μ*

_{1})

*T*

_{2}(

**r**) in projected attenuation between material 2 and material 1. The linear method, on the other hand, retrieves the total attenuation caused by all materials. This is also the reason for the correction coefficient

*A*(

**r**

_{⊥}) in the parallel method. This correction is missing in the linear method, not because of any additional approximations, but because it should not be present if the full attenuation is reconstructed (see Appendix A). Even though Eqs. (2) and (3) retrieve different quantities both the linear and the parallel algorithms will, in the end, return

*μ*.

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, *μ*_{2} − *μ*_{1}, the constant value of *μ*_{1} is added to the cut out G (see Fig. 1).

The parallel method’s final steps can be expressed as

where*M*are binary masks and

*M*

_{2}= 1 −

*M*

_{1}.

*V*are the reconstructed volumes. Indices 1 and 2 correspond to encasing material and encased material, respectively.

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

*V*is the volume image from the reconstruction,

*M*

_{H}and

*M*

_{L}are binary masks for segmentation of highly and less absorbing parts respectively, and

*K*

_{L}, the filter used in step 4, i.e. on the less absorbing part, is

The term Δ*δ*_{1}/Δ*μ*_{1} 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 Δ*δ*_{2}/Δ*μ*_{2} 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*_{0}, 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}/Δ*μ*_{1} = 1.7 · 10^{−9} m and Δ*δ*_{2}/Δ*μ*_{2} = 1 · 10^{−8} 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^{−8} m) followed by a threshold segmentation (threshold = 56 m^{−1}) 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^{−9} m, threshold 60 m^{−1}) 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]

*S̄*

_{A}is the mean pixel value over the whole object,

*S̄*

_{B}is the mean pixel value over the entire background, and

*σ*

_{B}is the standard deviation of the background. The same regions are hence compared for all images.

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.

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 *χ*^{2} 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 *χ*^{2} of zero [24]. We note that the absolute results of *χ*^{2} 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 *χ*^{2} 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.

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.

## 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^{−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^{−8} 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.

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.

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

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*_{0} are amplified, thus degrading the result.

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.

## 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]

*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,

*I*

_{0}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

*k*gives

*μ*,

*δ*and the thicknesses are rewritten as

*A*(

**r**) ≈ 0, yielding

**r**

_{⊥}, produces a factor $-4{\pi}^{2}{\mathbf{w}}_{\perp}^{2}$. The equation can be writen as

## Derivation of the linear phase-retrieval formula

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

*δ*as

*δ*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 where Δ

*δ*/Δ

*μ*is the slope and

*δ*

_{0}is the offset. Inserting the expressions in Eq. (29) into Eq. (27) yields

*k*and using ∇

*A*(

**r**) ≈ 0 due to slow variations of the total thickness gives

*M*(

**r**) in Eq. (28) yields the phase-retrieval formula used in the linear phase retrieval algorithm:

## 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]