Comparison of single distance phase retrieval algorithms by considering different object composition and the effect of statistical and structural noise

R. C. Chen; L. Rigon; R. Longo

doi:10.1364/OE.21.007384

1. Introduction

X-ray computed tomography (CT) is a non-destructive technique widely used for visualizing the morphology of samples, and for assessing quantitative information on their three-dimensional (3D) geometries and properties [1, 2]. With the availability of third generation Synchrotron Radiation sources, X-ray phase-sensitive imaging/CT has evolved as an increasingly accepted and utilized technique for characterizing the 3D internal structure of samples in the fields of material science and life sciences [3]. X-ray phase-sensitive imaging/CT utilizes phase shifts rather than absorption information, which is employed in X-ray absorption imaging/CT, as imaging signal [4]. Among X-ray phase-sensitive CT techniques, the X-ray propagation-based phase-contrast computed tomography (PPCT) has a very simple experimental setup, which requires no additional optics in the imaging geometry and is identical to the conventional absorption CT except for providing the beam is sufficiently spatially coherent and increasing the sample-to-detector distance (SDD) [3]. Qualitative PPCT can be performed by applying the standard ðltered back-projection CT reconstruction algorithm to the PPCT data [5]. The results are proportional to the Laplacian of the sample refractive index distribution providing an edge-enhancement, which allows visualizing the boundaries of regions with different refraction properties. Besides, PPCT radiographies fringes contain phase information which could be extracted by means of phase-retrieval [6]. Several phase retrieval algorithms have been developed [6–13]. In general, phase retrieval requires at least two intensity measurements, taken at two different SDDs, such as the transport of intensity equation (TIE) method [6], or the contrast transfer function method [9]. This is a limitation since typically hundreds or thousands of projections will be taken in PPCT experiment. Taking PPCT data at two different SDDs will increase the experiment time and deliver a higher radiation dose to the samples, which could hinder biomedical applications. Definitely, phase retrieval utilizing only one SDD PPCT data is preferable and more feasible, especially when dose is an essential issue in the experiment.

Several phase-retrieval algorithms using a single SDD PPCT data have been proposed, such as the Modified Bronnikov algorithm (MBA) method [13], which modifies the Bronnikov algorithm [11] by introducing an absorption correction factor (ACF) and eliminating the need of additional contact plan projections; the algorithms based on first Born- and Rytov-type approximations proposed by Gureyev [8]; the TIE based method by Paganin [12], which provides a method to reconstruct the sample projected thickness of the homogeneous sample using a single defocused image by solving TIE, thus simultaneously extracting phase and amplitude information; and the phase-attenuation duality algorithm proposed by Wu [14], which also considers the effects of x-ray source coherence and detector resolution. All these algorithms employ only a single SDD PPCT data set and utilize the same two assumptions on material properties: i) that the absorption is weak and homogenous, which means that the intensity in the contact plane can be approximated to unity, and ii) phase-attenuation duality (PAD) property [12, 14], i.e. the refractive index decrement $δ$ and the absorption index $β$ of complex refractive index $n = 1 - δ + i β$ are proportional to each other [15]. Other approaches can be found in [16–18].

Although comparison studies among phase retrieval algorithms, utilizing single- and multi- SDD(s) PPCT data, have been conducted [15, 19–21], it is interesting to know their performance considering the effect of statistical and structural noise, and different object composition. One trend in PPCT experiments is ultra-fast-CT, that aims to reduce the exposure time thus speeding up the data collection and reducing the radiation dose. When the exposure time is reduced, the statistical noise in the collected images increases accordingly and the quality of images becomes poorer. On the other side, a given sample may be surrounded by other uninteresting materials or structures, which have different shape and thickness and will affect the evaluation of the interesting objects. This is the case of structural noise, i.e. the presence in the image of unwanted and unavoidable structures besides the ones of interest. Furthermore, it is also attractive to investigate the performance of phase retrieval algorithms when their assumptions are not strictly satisfied.

In this paper, the performance of different single SDD phase retrieval algorithms will be studied. They are called phase-attenuation duality Born Algorithm (PAD-BA) [8, 22], phase-attenuation duality Rytov Algorithm (PAD-RA) [8], phase-attenuation duality Modified Bronnikov Algorithm (PAD-MBA) [11, 13, 23], phase-attenuation duality Paganin algorithm (PAD-PA) [12], and phase-attenuation duality Wu Algorithm (PAD-WA) [14, 24], respectively in this paper. They are all based on phase-attenuation duality property and on weak absorption of the sample and they employ only a single SDD PPCT data. These algorithms will be investigated via simulated noise-free PPCT data considering the fulfillment of PAD property and weakly absorbing conditions, and with experimental PPCT data of a mixture sample containing absorbing and weakly-absorbing materials, and of a polymer sample considering different degrees of statistical and structural noise.

2. Theory

As shown in Fig. 1 , an object is illuminated by a monochromatic plane x-ray beam, and the PPCT projection images are collected in the image plane. The object can be described by its 3D complex refractive index distribution, $n (x, y, z) = 1 - δ (x, y, z) + i β (x, y, z)$ , where $δ$ and $β$ are the refractive index decrement and the absorption index respectively, and $(x, y, z)$ are the spatial coordinates. Because of the weak interaction of x-rays with matter, the beam propagation path inside the sample can be assumed to be straight and $L$ denotes the linear path in the sample. The wave-object interaction can then be represented as the object transmission function [25]

T_{θ} (x, y) = \exp [- γ_{θ} (x, y) - i ϕ_{θ} (x, y)]

where

θ

represents the CT rotation angle; the sample phase

ϕ_{θ} (x, y)

function and absorption function

γ_{θ} (x, y)

are respectively:

\begin{array}{l} ϕ_{θ} (x, y) = k \int_{L θ}^{} δ (x, y, z) d z \\ γ_{θ} (x, y) = k \int_{L θ}^{} β (x, y, z) d z \end{array}

where

k = 2 π / λ

is the wavenumber,

λ

is the wavelength, and

\int_{L θ}^{}

denotes the line integral over the object along the beam path

L

at CT rotation angle θ.

Fig. 1 Schematic of PPCT scanning geometry.

Download Full Size | PDF

When a sample is quasi-homogeneous, the real and imaginary parts of $n - 1$ are proportional to each other (PAD property) [12, 24]:

δ (x, y, z) = ε β (x, y, z)

where

ε

is a constant. In this paper, a sample which fulfills PAD property will be referred to as a ratio object that is the ratio of the real and imaginary parts of

n - 1

is constant.

In the following, the principle of single SDD phase retrieval algorithms, PAD-BA, PAD-RA, PAD-MBA, PAD-PA and PAD-WA, will be introduced briefly. For describing convenience, let $I_{θ, z}$ and $I_{θ, 0}$ be the intensity distribution at SDD = z and SDD = 0 at rotation angle $θ$ of PPCT; $\hat{ϕ}$ and $\hat{γ}$ denote the Fourier transform of the phase and the absorption function, $ϕ$ and $γ$ , respectively.

2.1 Phase-attenuation duality Born algorithm

Let us assume an object of weak absorption and slowly varying phase shift. According to the Born-type approximation PPCI theory, $I_{θ, z}$ can be approximated by the following equation [8]

F [(I_{θ, z} / I_{θ, 0} - 1) / 2] (ξ, η) = {\hat{γ}}_{θ} \cos χ + {\hat{ϕ}}_{θ} \sin χ

where

(ξ, η)

are the spatial frequencies in the Fourier space corresponding to coordinates

(x, y)

in real space, and

χ = π λ z (ξ^{2} + η^{2})

.

In case the object fulfills PAD condition and is weakly absorbing, i.e. $I_{θ, 0} \approx 1$ , substituting Eq. (3) into Eq. (4) yields the following equation [8]:

ϕ_{θ} (x, y) = F^{- 1} {\frac{F [(I_{z, θ} - 1) / 2]}{ε^{- 1} \cos χ + \sin χ}}

After retrieving the phase function for the entire set of PPCT projections, the 3D refractive index can be reconstructed by applying the standard FBP algorithm to $ϕ_{θ} (x, y)$ [22], that is:

δ (x, y, z) = k^{- 1} \int_{0}^{π} ϕ_{θ} (x, y) * ν d θ

where

*

denotes 1D convolution and

ν

is the CT reconstruction filter.

2.2 Phase-attenuation duality Rytov algorithm

According to the Rytov-Type approximation of the PPCI theory, $I_{θ, z}$ can be approximated as [8]

F [\ln (I_{θ, z} / I_{θ, 0}) / 2] (ξ, η) = {\hat{γ}}_{θ} \cos χ + {\hat{ϕ}}_{θ} \sin χ

If the object satisfies PAD condition and is weakly absorbing ( $I_{θ, 0} \approx 1$ ), combining Eqs. (3) and (7) yields the following equation [8]:

ϕ_{θ} (x, y) = F^{- 1} {\frac{F [\ln (I_{z, θ}) / 2]}{ε^{- 1} \cos χ + \sin χ}}

After retrieving the phase function for the entire set of PPCT projections, the 3D refractive index can be reconstructed by applying the standard FBP algorithm to $ϕ_{θ} (x, y)$ as in Eq. (6).

2.3 Phase-attenuation duality modified Bronnikov algorithm

According to TIE theory [6] and considering the property of 2D and 3D Radon transform, the 3D refractive index decrement $δ (x, y, z)$ of the object can be reconstructed via the following equation that is the main result of the Bronnikov algorithm [11]:

δ (x, y, z) = \frac{1}{4 π^{2} z} \int_{0}^{π} [q (x, y) * * g_{θ} (x, y)] d θ

where

* *

indicate a 2D convolution,

g_{θ} (x, y) = I_{θ, z} / I_{θ, 0} - 1

, and

q (x, y) = \frac{| y |}{x^{2} + y^{2}}

is a filter. Bronnikov algorithm needs double SDDs PPCT data, i.e.

I_{θ, z}

and

I_{θ, 0}

in

g_{θ} (x, y)

, during the reconstruction. For weakly absorbing sample (

I_{θ, 0} \approx 1

), Groso modified the Bronnikov algorithm, by eliminating the need of

I_{θ, 0}

, i.e.

g_{θ} (x, y) = I_{θ, z} - 1

, and introducing an absorption correction factor (ACF: α) in the filter, thus only single SDD PPCT data is needed [13]. The Fourier space form of modified filter is [13]

Q (ξ, η) = \frac{| ξ |}{ξ^{2} + η^{2} + α}

where

α

is the ACF, which is determined by using a semi empirical approach in MBA [13]. However, an inappropriate α value will affect the reconstruction results, which will be blurred with a too small ACF value, while the filter will be eliminated with a too large value. Different groups have performed relative study on this topic and proposed a precise ACF value, which is based on the phase-attenuation duality property of low-Z sample, for MBA [23, 26, 27], namely:

α = \frac{1}{π ε λ z}

The PAD-MBA is the combination of Eq. (9) with $g_{θ} (x, y) = I_{θ, z} - 1$ , and Eqs. (10) and (11).

2.4 Phase-attenuation duality Paganin algorithm

PAD-PA provides a method to reconstruct the projected thickness $t (x, y)$ of the homogeneous sample using a single defocused image by solving the TIE [12], thus it simultaneously extracts phase and amplitude information. The projected thickness of the object $t (x, y)$ is extracted as follows:

t (x, y) = \frac{1}{μ} \ln [F^{- 1} (μ \frac{F (I_{z, θ} / I_{θ, 0})}{4 π^{2} z δ (ξ^{2} + η^{2}) + μ})]

For weakly absorbing sample ( $I_{θ, 0} \approx 1$ ), considering PAD property and substituting $μ = \frac{4 π}{λ} β$ and Eq. (12) into Eq. (2), the object phase function can be obtained via the following equation:

ϕ_{θ} (x, y) = \frac{1}{2} ε \ln {F^{- 1} [\frac{F (I_{z, θ})}{1 + π ε λ z (ξ^{2} + η^{2})}]}

After retrieving the phase function for the entire set of PPCT projections, the 3D refractive index can be reconstructed by applying the standard FBP algorithm to $ϕ_{θ} (x, y)$ as in Eq. (6).

2.5 Phase-attenuation duality Wu algorithm

PAD-WA is proposed by Wu using a single phase-contrast image [14]. Starting from either the paraxial Fresnel–Kirchhoff diffraction theory or the phase-space evolution of the Wigner distributions for x-ray wave fields, object phase function can be retrieved via following formula:

ϕ_{θ} (x, y) = \frac{λ r_{e}}{σ_{K N}} \ln {F^{- 1} [F (I_{z, θ}) / {I_{i n} {\tilde{μ}}_{i n} (λ z) O T F_{\det} [1 + \frac{2 π λ^{2} z r_{e}}{σ_{K N}} (ξ^{2} + η^{2})]}]}

where

r_{e}

is the classical electron radius,

σ_{K N}

denotes the total cross section for x-ray Compton scattering from a single free electron,

{\tilde{μ}}_{i n} (λ z)

presents spatial coherence of the incident x-ray,

O T F_{\det}

is the detector spatial frequency response, and

I_{i n}

denotes the intensity of the incident x-ray upon the object.

Let us assume a perfect coherent, homogeneous incident x-ray source and a perfect detector, i.e. ${\tilde{μ}}_{i n} \approx 1, I_{i n} \approx 1$ , and $O T F_{\det} \approx 1$ . Considering $\frac{r_{e}}{σ_{K N}} = \frac{ε}{2 λ}$ [14,15], the object phase function can be extracted via the following equation:

ϕ_{θ} (x, y) = \frac{1}{2} ε \ln {F^{- 1} [\frac{F (I_{z, θ})}{1 + π ε λ z (ξ^{2} + η^{2})}]}

As mentioned before, the 3D refractive index can be reconstructed by applying the standard FBP algorithm to $ϕ_{θ} (x, y)$ , as in Eq. (6).

It should be noted that, with PAD property and weakly absorbing approximation, the final single SDD phase retrieval equations of PAD-PA and PAD-WA are the same, as shown in Eqs. (13) and (15); thus, in the following study, we will not distinguish between these two algorithms and we will refer to them as PAD-PA/WA.

3. Materials and methods

3.1 Simulations

The performance of all phase retrieval algorithms was first investigated by simulation. The noise-free PPCT data were generated via tomography projection theory and the Fresnel diffraction theory. The simulated phantom is shown in Fig. 2 , in which Fig. 2(a) is the 3D phantom, and Figs. 2(b)-2(d) are phantom slices according to the line positions in Fig. 2(a). As the images show, the 3D phantom is made up of an ellipsoid and two spheres inside. Three different cases are conducted by assigning different complex refractive index values, i.e. the $(δ, β)$ values, to the three phantom regions. The first one is a weakly absorbing ratio phantom ( $ε = 1000$ ), in which the $(δ, β)$ values are (0.0, 0.0) (black, background), (1.0 × 10⁻⁷, 1.0 × 10⁻¹⁰), (2.0 × 10⁻⁷, 2.0 × 10⁻¹⁰), (3.0 × 10⁻⁷, 3.0 × 10⁻¹⁰) (white) respectively; while in the second one, the $(δ, β)$ values are (0.0, 0.0), (1.0 × 10⁻⁷, 1.0 × 10⁻⁹), (2.0 × 10⁻⁷, 2.0 × 10⁻⁹), (3.0 × 10⁻⁷, 3.0 × 10⁻⁹) respectively for the phantom, which simulates an absorbing ratio phantom ( $ε = 100$ ); the final one is a weakly absorbing no-ratio phantom, in which the $(δ, β)$ values are (0.0, 0.0), (1.0 × 10⁻⁷, 1.2 × 10⁻¹⁰), (2.0 × 10⁻⁷, 2.0 × 10⁻¹⁰), (3.0 × 10⁻⁷, 2.0 × 10⁻¹⁰) respectively. The simulation parameters of all PPCT data sets are: energy of 14 keV, SDD = 0.6 m, and an effective pixel size of 9 μm. Within 180 degree CT scan range, 220 PPCT projections were generated.

Fig. 2 Simulation phantom: (a) 3D phantom, (b)-(d) phantom slices according to the line positions in (a).

Download Full Size | PDF

3.2 Experiments

The experimental PPCT data were collected at the SYRMEP beamline [28] at the ELETTRA synchrotron facility, Italy. The SYRMEP beamline employs a bending magnet source with a Si (111) double-crystal monochromator, which provides photon energy ranging from 8.5 to 35 keV.

3.2.1 Mixture sample containing absorbing and weakly-absorbing materials

The assumptions of the single SDD phase retrieval algorithms are the object must fulfill PAD property, i.e. it is a ratio object, and weakly-absorbing. In this case, it is interesting to investigate their performance when the object does not fulfil these conditions. Here, a mixture sample is investigated. It is a PMMA cylinder where five holes have been drilled; four of the five holes are filled with aluminum, Teflon, water and Polyoxymethylene, while the other hole is left empty. The energy of 25 keV, SDD = 0.5 m and a CCD detector with an effective pixel size of 9 μm were used to acquire all PPCT data sets. At 25 keV, the $ε$ values for aluminum, Teflon, water, Polyoxymethylene, and PMMA are 490, 1663, 2152, 2464 and 2787 respectively. In this case, this sample cannot be treated as ratio object since in particular the $ε$ value of aluminum is much smaller than the others. Moreover, aluminum is an absorbing material ( $β =$ 1.8e-09 @ 25keV) compared to other polymer materials ( $β \approx$ 1.0e-10 @ 25keV).

3.2.2 Statistical and structural noise study with a polymer sample

A polymer sample containing wires of nylon ( $\emptyset = 1.6$ mm), polystyrene ( $\emptyset = 1.6$ mm) and PMMA ( $\emptyset = 2$ mm) was investigated. PPCT data sets were collected with different level of statistical and structural noise, which will be explained in the following. The energy of 16 keV, SDD = 0.8 m and a CCD detector with an effective pixel size of 9 μm were used to acquire all PPCT data sets for statistical noise and structural noise study. Within 180 degree CT scan range, 900 PPCT projections were collected.

For statistical noise PPCT data collection, we adjusted the exposure time and the filter thickness in order to obtain different fluence (photons/pixel) reach the detector. Eight PPCT data sets were collected: their parameters are summarized in Table 1 .

Table 1. Fluence (photons/pixel) in the background for the eight PPCT data sets collected for the statistical noise study.

View Table | View all tables in this article

For the structural noise study, the PPCT data sets were acquired by putting a structural noise phantom, Polyamide powder as presented in Fig. 3 , behind the sample during PPCT data collection. As Fig. 3 shows, five different thicknesses of the structural noise phantom were employed in the experiment. Twelve PPCT data sets were collected. The parameters of collected PPCT data sets are summarized in Table 2 .

Fig. 3 Photo of structural noise phantom showing five different thicknesses of Polyamide powder.

Download Full Size | PDF

Table 2. The parameters of collected PPCT data sets for structural noise study.

View Table | View all tables in this article

3.3 PPCT data processing

For PAD-BA, PAD-RA, and PAD-PA/WA, the phantom/object phase function was retrieved by applying the phase retrieval algorithms to the PPCT projections, and then the 3D refractive index was reconstructed by implementing the standard filter back-projection algorithm with Shepp-Logan filter. On the other hand, for PAD-MBA, the processing procedure was firstly filtering the $g_{θ} (x, y) = I_{θ, z} - 1$ using the filter in Eq. (10) with $α = \frac{1}{π ε λ z}$ , and then applying the back-projection to obtain the 3D refractive index. During the phase retrieval, the $ε$ values of 1000, 1000 and 100 are used for simulated weakly absorbing ratio phantom, weakly absorbing no-ratio phantom and absorbing ratio phantom PPCT data respectively, and $ε$ values of 1663 and 1739 are used respectively for mixture sample and sample polymer PPCT data.

All data processing was performed via PITRE software package [29]. PITRE is a freeware which supports phase retrieval for PPCT projections, extracts apparent absorption, refractive and scattering information of diffraction enhanced imaging, and allows parallel beam tomography reconstruction for conventional absorption CT data and for PPCT phase retrieved and DEI-CT extracted information.

3.4 Data analysis

For the quantitative analysis of the results, i.e. slice in this case, the contrast-to-noise ratio (CNR) was measured for the slices reconstructed with the single SDD phase retrieval algorithms. CNR is defined by

CNR = \frac{| S_{a} - S_{b} |}{\sqrt{σ_{a}^{2} + σ_{b}^{2}}}

where

S_{a}

and

S_{b}

are the mean values of two homogeneous region of interest (ROI) in the slices,

a

and

b

representing two different materials; while

σ_{a}^{}

and

σ_{b}^{}

denote the standard deviations of the values in these two ROIs of identical size.

Moreover, the histograms of the reconstructed slices for PAD-BA, PAD-RA, PAD-MBA and PAD-PA/WA obtained from the same PPCT data are calculated. In the histograms, each material in the reconstructed slice will show up as a peak. Histogram entropies have several properties which enable their use for resolution identification. For instance, the width of the peak directly represents the density resolution: the narrower the peak, the higher the density resolution [30, 31].

4. Results

4.1 Simulations

Figures 4(a) -4(d) present the reconstructed phantom slices of noise-free weakly absorbing ratio phantom, the position are same as Fig. 2(c), after implementing the PAD-BA, PAD-MBA, PAD-RA and PAD-PA/WA respectively. Figure 4(f) is profile of Figs. 4(a)-4(d) and of the phantom itself, i.e. Figure 2(c), at line position shown in Fig. 4(a). As the image shows, for the noise-free weakly absorbing ratio phantom PPCT data, after implementing phase retrieval, the refractive index matches the phantom values well for all algorithms. This indicates that PAD-BA, PAD-MBA, PAD-RA and PAD-PA/WA can quantitatively reconstruct the 3D refractive index of a weakly absorbing ratio object by utilizing single SDD PPCT data.

Fig. 4 Noise-free weakly absorbing ratio phantom simulation results: (a)-(d) are phase retrieval results, same slice as Fig. 2(c), of PAD-BA, PAD-MBA, PAD-RA and PAD-PA/WA (ε = 1000) respectively, (f) profile of (a)-(d) and of the phantom itself, i.e. Figure 2(c), at line position shown in Fig. 4(a).

Download Full Size | PDF

Figures 5(a) and 5(b) are profiles of weakly absorbing no-ratio phantom and absorbing ratio phantom simulation results, same slice as Fig. 2(c), after implementing PAD-BA, PAD-MBA, PAD-RA and PAD-PA/WA phase retrieval algorithm, respectively. From the visualization, the phase retrieval slices of all algorithms for these two phantoms are similar to the ones of weakly absorbing ratio phantom, i.e. Figures 4(a)-4(d), therefore they are not presented here.

Fig. 5 (a) Profile of weakly absorbing no-ratio phantom simulation results, same slice as Fig. 2(c), after implementing PAD-BA, PAD-MBA, PAD-RA and PAD-PA/WA (with ε = 1000) and phantom itself, i.e. Fig. 2(c), at same line position shown in Fig. 4(a), (b) Profile of absorbing ratio phantom simulation results, same slice as Fig. 2(c), after implementing PAD-BA, PAD-MBA, PAD-RA and PAD-PA/WA (with ε = 100) and phantom itself, i.e. Fig. 2(c), at same line position shown in Fig. 4(a).

Download Full Size | PDF

From Fig. 5(a), it is clear that although slightly differences exist among different algorithms, they all have problem in reconstructing the correct refractive index for this weakly absorbing no-ratio phantom. These are expected since this phantom does not satisfy the ratio object condition, i.e. PAD property, of phase retrieval algorithms. In this case, although the quantitative result will not be preserved, quantitatively reasonable results may still be obtained. These results improve the contrast and could be used for distinguishing different materials.

When investigating the absorbing ratio phantom, the results of PAD-RA and PAD-PA/WA match the phantom values well as shown in Fig. 5(b), while PAD-BA and PAD-MBA have clearly a problem. When comprising the PAD-BA and PAD-RA phase retrieval algorithms, the only difference between them is the use of $(I_{z, θ} - 1)$ in PAD-BA instead of $\ln (I_{z, θ})$ in PAD-RA; moreover a similar situation exists in PAD-PA/WA and PAD-MBA. When the contrast of recording PPCT projection ( $I_{z, θ}$ ) is low enough, then we can have $\ln (I_{z, θ}) \approx I_{z, θ} - 1$ [8]. Inspired by these points and considering that the contrast of this absorbing ratio phantom is not negligible, we reformulate $g_{θ} (x, y) = I_{θ, z} - 1$ to $g_{θ} (x, y) = \ln (I_{θ, z})$ in PAD-MBA while keeping the other parameters in order to figure out the problem. The same simulated PPCT data was reconstructed with this trial phase retrieval algorithm, which is named as Log-PAD-MBA. The comparison of PAD-MBA and Log-PAD-MBA is shown in Fig. 6 , in which it is clear that after implementing the Log-PAD-MBA, the result has substantially improved compared to PAD-MBA one. This trial indicates that the reason PAD-RA and PAD-PA/WA obtain good results for absorbing ratio phantom is because they use a logarithm considering the contrast of recording PPCT projection is not negligible.

Fig. 6 Profile of absorbing ratio phantom simulation results, same slice as Fig. 2(c), after implementing PAD-MBA and Log-PAD-MBA (with ε = 100) and phantom itself, i.e. Figure 2(c), at same line position shown in Fig. 4(a).

Download Full Size | PDF

4.2 Experiments

Figures 7(a) -7(d) show the reconstructed slices of the mixture sample after applying PAD-BA, PAD-MBA, PAD-RA and PAD-PA/WA phase retrieval algorithm respectively. In the image, it is clearly seen that PAD-BA and PAD-MBA results are corrupted by artefacts, especially near the aluminium rod (the brightest one). On the other hand, the results of PAD-RA and PAD-PA/WA show homogeneous distribution for all materials. Comparing with the simulations results, this is due to the fact that the phantom does not fulfil the PAD property and weakly absorbing assumptions, and the absorption affects the result more than the violation of the ratio object condition.

Fig. 7 Mixture sample reconstructed slices: (a)-(d) are using PAD-BA, PAD-MBA, PAD-RA and PAD-PA/WA (ε = 1663) phase retrieval algorithm respectively.

Download Full Size | PDF

Figures 8(a) -8(d) shows the histogram of reconstructed slices of the mixture sample using PAD-BA, PAD-RA, PAD-MBA and PAD-PA/WA phase retrieval algorithm respectively. The images show that, PPMA and water cannot be distinguished in all cases since their $δ$ values are too close; Polyoxymethylene has a nice sharp peak in PAD-RA and PAD-PA/WA results, and a small peak in PAD-BA result, but not in PAD-MBA one; Teflon and aluminum have clear peak in PAD-RA and PAD-PA/WA results, but these peaks are diminishing sharply in PAD-BA and PAD-MBA ones, especially for aluminum. Moreover, it should be noted that the PAD-MBA reconstructed a negative value for the background/air, which it is not correct in the practical case. This indicates that the PAD-MBA meets problems when treating this kind of sample. The other three algorithms obtained positive value for all materials. These results indicate that PAD-RA and PAD-PA/WA algorithms have more tolerance for samples that do not fulfil their conditions, especially when containing absorbing materials. This could be very useful since the practical sample may be comprised of many different materials and cannot satisfy the ideal requirements. On the other hand, PAD-BA and PAD-MBA algorithms are more sensitive for conditions violations.

Fig. 8 Histogram of reconstructed slices of PAD-BA, PAD-RA, PAD-MBA and PAD-PA/WA (ε = 1663) in Fig. 7.

Download Full Size | PDF

Figure 9 shows the reconstructed slices of statistical noise study of polymer sample. More in detail, Figs. 9(a)-9(c) are phase contrast slices of statistical noise PPCT data index 1, 5 and 8 respectively (see Table 1), while Figs. 9(d)-9(f) are phase retrieval slices obtained with the PAD-BA algorithm of same PPCT data as in Figs. 9(a)-9(c). From the visualization, it is hard to distinguish among phase retrieval slices of PAD-BA, PAD-MBA, PAD-RA and PAD-PA/WA; here only PAD-BA slices are presented. As Fig. 9 shows, with the decrease of counted photons, the result is progressively corrupted by the noise, more obviously in phase contrast slices than in phase retrieval ones. However, the visibilities of all the slices are well preserved in both cases. The reconstructed refractive index values of PMMA, nylon and polystyrene are 1.05e-6, 9.08e-7, 6.74e-7 respectively from slice of PPCT data index 1 with PAD-BA algorithm applied; these values do not match their ideal values that are 1.04e-06, 1.01e-06 and 9.07e-07 respectively for three materials. This could be expected since this polymer sample does not strictly satisfy the required conditions, especially the PAD property, of phase retrieval algorithms. In this case, although the quantitative result will not be preserved, it is obvious the contrast is improved that could be used for distinguishing different materials.

Fig. 9 Polymer sample reconstructed slices: (a)-(c) are phase contrast slices of statistical noise PPCT data index 1, 5 and 8 respectively (see Table 1); (d)-(f) are phase retrieval slices with PAD-BA algorithm and ε = 1739 of the same PPCT data as in (a)-(c)

Download Full Size | PDF

The CNR values for PMMA, nylon and polystyrene of statistical noise reconstructed slices are presented in Table 3 . For calculating the CNR value, material $a$ is one of PMMA, nylon or polystyrene, while material $b$ is air. Their ROIs are the rectangles marked in Fig. 9(d), which applies for both statistical noise and structural noise study of this sample. As expected, the CNR values of three materials decrease as the fluence (photons/pixel) reduces. When close checking the results, PAD-PA/WA gains highest CNR and the trend is ${CNR}_{P A D - P A / W A} {>CNR}_{P A D - R A} {>CNR}_{P A D - M B A} {>CNR}_{P A D - B A}$ . This fact is explained by their histogram, shown in Fig. 10 , in which their density resolution has obvious differences.

Table 3. CNR values for PMMA, nylon or polystyrene of reconstructed slices of polymer sample in the statistical noise study.

View Table | View all tables in this article

Fig. 10 Histogram of reconstructed slices of PPCT data index 8 in Table 1 with different phase retrieval algorithm: (a)-(d) are of PAD-BA, PAD-RA, PAD-MBA and PAD-PA/WA (ε = 1739) respectively. Each histogram identifies all densities on the image in the form of a graph: in this case, X-axis is related to $δ$ value, while Y-axis displays the number of pixels for each $δ$ value.

Download Full Size | PDF

Figure 10 shows the histogram of reconstructed slices of PAD-BA, PAD-RA, PAD-MBA and PAD-PA/WA of PPCT data index 8 in Table 1 of polymer sample. As the image shows, the result of PAD-RA and PAD-PA/WA have narrower peak compared to PAD-BA and PAD-MBA ones, which means PAD-RA and PAD-PA/WA perform better in this statistical noise study than PAD-BA and PAD-MBA. This is likely due to the use of logarithm in PAD-RA and PAD-PA/WA, since they do not need a further assumption, i.e. that the recording PPCT projection ( $I_{z, θ}$ ) should be low enough to fulfill $\ln (I_{z, θ}) \approx I_{z, θ} - 1$ , compared to PAD-BA and PAD-MBA.

Figures 11(a) -11(c) are phase contrast slices of structural noise PPCT data index 3, 7 and 11 respectively (see Table 2), while Figs. 11(d)-11(f) are phase retrieval slices obtained with the PAD-BA algorithm of the same PPCT data as in Figs. 11(a)-11(c). As Fig. 11 shows, with the increase of structural noise (corresponding to increase in phantom thickness), the result is progressively corrupted by the noise. In this case it appears as ring-like artifacts in the slice. This is because the position and structure of structural phantom does not change during a complete PPCT data collection, thus the distribution of structural noise in all the projections is the same. This will generate vertical stripes in the sinogram, and they will appear as ring artifacts superimposed on the slices after back-projection. In the phase contrast slice, the edge enhancement is obviously decreased as compared to the statistical noise study result with roughly the same fluence. This is because the structural noise phantom not only introduces ring-like artifacts, but also acts as a diffuser, i.e. it changes the photon propagation direction, thus reducing the edge enhancement effect.

Fig. 11 Polymer sample reconstructed slices: (a)-(c) are phase contrast slices of structural noise PPCT data index 3, 7 and 11 respectively (see Table 2); (d)-(f) are phase retrieval slices with PAD-BA (ε = 1739) algorithm of the same data as in (a)-(c).

Download Full Size | PDF

The CNR values for PMMA, nylon or polystyrene of structural noise reconstructed slices are presented in Table 4 . The CNR calculation rules are the same used in the statistical noise study. From Table 4, it is clear to see that there are not significant differences of CNR values for all algorithms and small differences with the changes of structural noise phantom thickness. There are two possible reasons: one is the ring-like artifact in the slice, that will increase the value in the denominator of Eq. (16); the other one is the reduced edge enhancement, from which the contrast will be recovered after phase retrieval, caused by the structural phantom.

Table 4. CNR values for PMMA, nylon or polystyrene of structural noise reconstructed slices of polymer sample.

View Table | View all tables in this article

Figure 12 shows the histogram of reconstructed slices of PAD-BA, PAD-RA, PAD-MBA and PAD-PA/WA of PPCT data index 11 in Table 2. Similarly to the statistical study result, PAD-RA and PAD-PA/WA obtain better density resolution than PAD-BA and PAD-MBA ones. Moreover, all the peaks for four algorithms are wider compare to the statistical noise result. This is due to the ring-like artifact caused by the structural noise phantom.

Fig. 12 Histogram of reconstructed slices of PPCT data index 11 in Table 2 with different phase retrieval algorithm: (a)-(d) are of PAD-BA, PAD-RA, PAD-MBA and PAD-PA/WA (ε = 1739) respectively.

Download Full Size | PDF

5. Discussions and conclusions

We investigated the performance of different single SDD phase retrieval algorithms, PAD-BA, PAD-MBA, PAD-RA and PAD-PA/WA, with simulated and experimental PPCT data. These algorithms require that the sample satisfy PAD property and weakly absorbing conditions. When the object fulfills these two assumptions, all algorithms can quantitatively reconstruct the 3D refractive index of a ratio object by utilizing single distance PPCT data. In case the PAD property is violated, they all have problems in reconstructing the correct refractive index although slight differences exist. Though the quantitative result will not be preserved, it is obvious the contrast is improved, compared to phase contrast slices, and this could be used for distinguishing different materials. When the weakly absorbing condition is not satisfied, which means the contrast of recording PPCT projection ( $I_{z, θ}$ ) is too large to satisfy $\ln (I_{z, θ}) \approx I_{z, θ} - 1$ approximation, the PAD-RA and PAD-PA/WA obtain better result than PAD-BA and PAD-MBA because the use of logarithm in their algorithms. This is confirmed in the experimental results of mixture sample.

When considering the statistical noise, the CNR values decrease as the fluence is reduced, and PAD-PA/WA and PAD-RA obtain higher CNR than PAD-BA and PAD-MBA. This is are likely due to the use of logarithm in PAD-RA and PAD-PA/WA, since they do not need the further approximation $\ln (I_{z, θ}) \approx I_{z, θ} - 1$ . The structural noise study shows that with the increase of structural noise, the result is progressively corrupted by ring-like artifact due to the effect of the structural noise phantom. Moreover, because the structural noise phantom acts as a diffuser, the edge enhancement is obviously decreased in the phase contrast slice and thus decreases the CNR values of phase retrieval results.

In summary, we conducted the comparison of single distance phase retrieval algorithms, on the reconstructed 3D slices, by considering their PAD property and weakly absorbing conditions, and the effect of statistical and structural noise. The results show that PAD-RA and PAD-PA/WA obtain better result than PAD-BA and PAD-MBA in the situations considered in this paper.

Acknowledgments

The authors are grateful to Ralf Menk and Fulvia Arfelli for providing the structural noise phantom. They are indebted to Lucia Mancini for preparing the mixture phantom. We wish to thank Liberato De Caro, Sabina Tangaro and Cinzia Giannini for helpful discussion. RCC was supported by the National Basic Research Program (973 Program) (no 2010CB834301).

References and links

1. A. C. Kak and M. Slaney, Principles of Computerized Tomographic Imaging (Society of Industrial and Applied Mathematics, 2001).

2. R. C. Chen, R. Longo, L. Rigon, F. Zanconati, A. De Pellegrin, F. Arfelli, D. Dreossi, R. H. Menk, E. Vallazza, T. Q. Xiao, and E. Castelli, “Measurement of the linear attenuation coefficients of breast tissues by synchrotron radiation computed tomography,” Phys. Med. Biol. 55(17), 4993–5005 (2010). [CrossRef] [PubMed]

3. J. Baruchel, J. Y. Buffiere, P. Cloetens, M. Di Michiel, E. Ferrie, W. Ludwig, E. Maire, and L. Salvo, “Advances in synchrotron radiation microtomography,” Scr. Mater. 55(1), 41–46 (2006). [CrossRef]

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

5. P. Cloetens, M. Pateyron-Salomé, J. Y. Buffière, G. Peix, J. Baruchel, F. Peyrin, and M. Schlenker, “Observation of microstructure and damage in materials by phase sensitive radiography and tomography,” J. Appl. Phys. 81(9), 5878–5886 (1997). [CrossRef]

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

7. X. Z. Wu and H. Liu, “A general theoretical formalism for X-ray phase contrast imaging,” J. XRay Sci. Technol. 11(1), 33–42 (2003). [PubMed]

8. T. E. Gureyev, T. J. Davis, 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]

9. P. Cloetens, W. Ludwig, J. Baruchel, D. Van Dyck, J. Van Landuyt, J. P. Guigay, and M. Schlenker, “Holotomography: Quantitative phase tomography with micrometer resolution using hard synchrotron radiation x rays,” Appl. Phys. Lett. 75(19), 2912–2914 (1999). [CrossRef]

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

11. A. V. Bronnikov, “Theory of quantitative phase-contrast computed tomography,” J. Opt. Soc. Am. A 19(3), 472–480 (2002). [CrossRef] [PubMed]

12. 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(1), 33–40 (2002). [CrossRef] [PubMed]

13. A. Groso, R. Abela, and M. Stampanoni, “Implementation of a fast method for high resolution phase contrast tomography,” Opt. Express 14(18), 8103–8110 (2006). [CrossRef] [PubMed]

14. X. Z. Wu and H. Liu, “X-Ray cone-beam phase tomography formulas based on phase-attenuation duality,” Opt. Express 13(16), 6000–6014 (2005). [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. L. D. Turner, B. B. Dhal, J. P. Hayes, A. P. Mancuso, K. A. Nugent, D. Paterson, R. E. Scholten, C. Q. Tran, and A. G. Peele, “X-ray phase imaging: Demonstration of extended conditions for homogeneous objects,” Opt. Express 12(13), 2960–2965 (2004). [CrossRef] [PubMed]

17. M. A. Beltran, D. M. Paganin, K. Uesugi, and M. J. 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]

18. J. Moosmann, R. Hofmann, A. V. Bronnikov, and T. Baumbach, “Nonlinear phase retrieval from single-distance radiograph,” Opt. Express 18(25), 25771–25785 (2010). [CrossRef] [PubMed]

19. A. M. Yan, X. Z. Wu, and H. Liu, “Robustness of phase retrieval methods in x-ray phase contrast imaging: A comparison,” Med. Phys. 38(9), 5073–5080 (2011). [CrossRef] [PubMed]

20. M. Langer, P. Cloetens, J. P. Guigay, and F. Peyrin, “Quantitative comparison of direct phase retrieval algorithms in in-line phase tomography,” Med. Phys. 35(10), 4556–4566 (2008). [CrossRef] [PubMed]

21. M. N. Boone, W. Devulder, M. Dierick, L. Brabant, E. Pauwels, and L. Van Hoorebeke, “Comparison of two single-image phase-retrieval algorithms for in-line x-ray phase-contrast imaging,” J. Opt. Soc. Am. A 29(12), 2667–2672 (2012). [CrossRef]

22. R. C. Chen, L. Rigon, and R. Longo, “Quantitative 3D refractive index decrement reconstruction using single-distance phase-contrast tomography data,” J. Phys. D Appl. Phys. 44(49), 495401 (2011). [CrossRef]

23. R. C. Chen, H. L. Xie, L. Rigon, R. Longo, E. Castelli, and T. Q. Xiao, “Phase retrieval in quantitative x-ray microtomography with a single sample-to-detector distance,” Opt. Lett. 36(9), 1719–1721 (2011). [CrossRef] [PubMed]

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

25. M. Born and E. Wolf, Principles of Optics, 7th ed. (Cambridge University, New York, 1999).

26. T. E. Gureyev, D. M. Paganin, G. R. Myers, Y. I. Nesterets, and S. W. Wilkins, “Phase-and-amplitude computer tomography,” Appl. Phys. Lett. 89(3), 034102 (2006). [CrossRef]

27. B. D. Arhatari, F. De Carlo, and A. G. Peele, “Direct quantitative tomographic reconstruction for weakly absorbing homogeneous phase objects,” Rev. Sci. Instrum. 78(5), 053701 (2007). [CrossRef] [PubMed]

28. A. Abrami, F. Arfelli, R. C. Barroso, A. Bergamaschi, F. Bille, P. Bregant, F. Brizzi, K. Casarin, E. Castelli, V. Chenda, L. D. Palma, D. Dreossi, A. Fava, R. Longo, L. Mancini, R. H. Menk, F. Montanari, A. Olivo, S. Pani, A. Pillon, E. Quai, S. R. Kaiser, L. Rigon, T. Rokvic, M. Tonutti, G. Tromba, A. Vaseotto, C. Venanzi, F. Zanconati, A. Zanetti, and F. Zanini, “Medical applications of synchrotron radiation at the SYRMEP beamline of ELETTRA,” Nucl Instrum Meth A 548(1-2), 221–227 (2005). [CrossRef]

29. R. C. Chen, D. Dreossi, L. Mancini, R. Menk, L. Rigon, T. Q. Xiao, and R. Longo, “PITRE: software for phase-sensitive X-ray image processing and tomography reconstruction,” J. Synchrotron Radiat. 19(5), 836–845 (2012). [CrossRef] [PubMed]

30. P. Thurner, F. Beckmann, and B. Muller, “An optimization procedure for spatial and density resolution in hard X-ray micro-computed tomography,” Nucl Instrum Meth B 225(4), 599–603 (2004). [CrossRef]

31. E. Hadjidemetriou, M. D. Grossberg, and S. K. Nayar, “Resolution selection using generalized entropies of multiresolution histograms,” Computer Vison Eccv 2350, 220–235 (2002).

PPCT data index	1	2	3	4	5	6	7	8	9	10	11	12
SNPT^a (mm)	0.60	0.60	0.60	0.94	0.94	1.96	1.96	2.97	2.97	5.03	5.03	5.03
Fluence (photons/pixel)	5270	3201	2398	5364	1955	2975	2390	3890	1956	3304	1956	1247

	PAD-BA			PAD-MBA			PAD-RA			PAD-PA/WA
Index	A^a	B^b	C^c	A	B	C	A	B	C	A	B	C
1	73.2	51.0	53.9	88.6	62.0	67.7	162.2	124.7	98.4	181.8	129.3	100.5
2	72.7	50.8	48.5	87.1	61.6	60.1	139.3	109.0	93.0	156.0	107.3	99.9
3	69.0	49.1	43.2	82.6	58.8	54.3	139.1	97.2	88.1	151.1	103.6	90.0
4	59.7	46.9	43.1	70.3	56.4	52.8	124.7	97.2	76.5	137.3	101.3	80.0
5	55.8	46.6	43.0	64.9	56.2	52.2	122.8	93.3	74.4	127.3	98.4	76.8
6	51.3	45.9	42.4	60.1	56.1	51.4	102.9	84.0	67.0	101.5	86.3	66.9
7	47.1	38.4	39.2	54.2	45.2	47.9	95.3	71.4	56.5	97.0	71.1	57.3
8	46.7	31.8	29.3	52.2	35.2	32.0	80.6	50.4	34.5	80.5	50.1	35.6

	PAD-BA			PAD-MBA			PAD-RA			PAD-PA/WA
Index	A^a	B^b	C^c	A	B	C	A	B	C	A	B	C
1	27.4	17.3	13.3	28.1	17.8	13.7	26.9	17.9	13.6	28.2	18.3	13.9
2	38.0	21.3	17.4	39.2	22.6	18.4	37.8	25.4	19.6	39.1	26.0	20.5
3	29.9	21.8	18.3	31.0	23.1	19.2	32.4	25.8	19.8	32.9	26.3	20.2
4	40.9	27.0	21.8	41.3	28.7	23.1	31.0	30.1	21.5	31.5	31.0	22.8
5	38.3	23.8	18.7	38.3	24.9	19.4	31.1	27.2	19.5	32.1	27.6	20.4
6	26.5	27.9	18.1	27.5	29.0	18.9	29.9	28.7	19.7	30.2	28.8	19.7
7	27.6	22.0	16.8	28.9	23.7	18.0	29.2	26.7	19.3	29.4	27.0	19.5
8	27.0	20.8	16.4	27.9	22.1	17.4	28.7	25.6	18.8	29.0	25.9	18.9
9	28.0	22.7	18.9	29.5	23.6	19.9	34.3	23.1	21.6	34.8	24.2	21.9
10	31.4	18.3	16.7	31.8	18.7	17.1	35.6	20.2	18.7	37.2	19.7	18.8
11	26.4	18.9	18.6	27.3	19.3	19.3	30.7	19.8	20.7	29.8	19.7	20.5
12	28.5	20.7	17.8	29.2	21.6	18.6	31.1	24.7	20.1	30.7	24.3	20.8

PPCT data index	1	2	3	4	5	6	7	8	9	10	11	12
SNPT^a (mm)	0.60	0.60	0.60	0.94	0.94	1.96	1.96	2.97	2.97	5.03	5.03	5.03
Fluence (photons/pixel)	5270	3201	2398	5364	1955	2975	2390	3890	1956	3304	1956	1247

	PAD-BA			PAD-MBA			PAD-RA			PAD-PA/WA
Index	A^a	B^b	C^c	A	B	C	A	B	C	A	B	C
1	73.2	51.0	53.9	88.6	62.0	67.7	162.2	124.7	98.4	181.8	129.3	100.5
2	72.7	50.8	48.5	87.1	61.6	60.1	139.3	109.0	93.0	156.0	107.3	99.9
3	69.0	49.1	43.2	82.6	58.8	54.3	139.1	97.2	88.1	151.1	103.6	90.0
4	59.7	46.9	43.1	70.3	56.4	52.8	124.7	97.2	76.5	137.3	101.3	80.0
5	55.8	46.6	43.0	64.9	56.2	52.2	122.8	93.3	74.4	127.3	98.4	76.8
6	51.3	45.9	42.4	60.1	56.1	51.4	102.9	84.0	67.0	101.5	86.3	66.9
7	47.1	38.4	39.2	54.2	45.2	47.9	95.3	71.4	56.5	97.0	71.1	57.3
8	46.7	31.8	29.3	52.2	35.2	32.0	80.6	50.4	34.5	80.5	50.1	35.6

Comparison of single distance phase retrieval algorithms by considering different object composition and the effect of statistical and structural noise

Abstract

1. Introduction

2. Theory

2.1 Phase-attenuation duality Born algorithm

2.2 Phase-attenuation duality Rytov algorithm

2.3 Phase-attenuation duality modified Bronnikov algorithm

2.4 Phase-attenuation duality Paganin algorithm

2.5 Phase-attenuation duality Wu algorithm

3. Materials and methods

3.1 Simulations

3.2 Experiments

3.2.1 Mixture sample containing absorbing and weakly-absorbing materials

3.2.2 Statistical and structural noise study with a polymer sample

3.3 PPCT data processing

3.4 Data analysis

4. Results

4.1 Simulations

4.2 Experiments

5. Discussions and conclusions

Acknowledgments

References and links

Cited By

Figures (12)

Tables (4)

Equations (16)

Optics Express