1. Introduction

Absorption based x-ray computed tomography (CT) is nowadays a standard technique with spatial resolution around one micron [1–4]. On the other hand, a wide range of samples studied in materials science, biology and medicine show very weak absorption contrast, yet produce significant phase shifts of the x-ray beam. The use of phase information can therefore provide new and otherwise not accessible information. Phase contrast imaging techniques have two key advantages: first, light elements - showing poor contrast in absorption radiography - can be easily detected and second, phase-contrast radiography helps to reduce the radiation dose deposited on the object under investigation.

As far as phase tomography (quantitative 3D reconstruction of the phase or the refractive index from 2D phase images) is concerned, several attempts were made by using interferometric [5–8] and non-interferometric phase retrieval methods [9–17]. All these methods are based on a two step approach: first, the projections of the phase are determined in the form of Radon projections and then the object function, i.e. the refractive index δ is reconstructed applying a conventional filtered backprojection algorithm. Since far from absorption edges δ is linearly related to the electron density ρ, which is - except for hydrogen rich materials - proportional to the mass density, the resulting reconstruction represents approximately the distribution of the mass density in the sample [18].

The reconstruction algorithm suggested by Bronnikov in his theoretical works [19, 20] presents an alternative approach which requires no intermediate step of 2D phase retrieval and provides a direct 3D reconstruction of the refractive index from the intensity distributions (radiographic projections at different angles) measured in a single plane (pure phase object) and from the intensity distributions measured in two planes in the case of mixed phase-amplitude object.

In a previous work [21] we have simulated and validated the quantitative phase tomography method described above. We have adapted the original algorithm in order to get 3D reconstruction of the refractive index from the intensity distributions measured in a single plane also in the case of mixed weakly absorbing-phase object. When using this one step approach in phase reconstruction in 3D the experimental setup can be kept extremely simple and the radiation damage is kept to the minimum, which is of special importance for biological specimens. Besides, phase retrieval as well as data alignment (necessary in case of measurements in multiple planes) are superfluous.

In this work we describe in detail the filter adaptation and the means of choosing the filter parameters. The method is applied to real experimental data and also validated by the Differential Phase Contrast Imaging Technique [8].

2. Theoretical background

When a monochromatic plane wave propagates along the positive z-axis and impinges upon a thin mixed phase-amplitude object (characterized by the linear absorption coefficient μ(x,y,z) and the real part of the index of refraction δ(x, y, z)), the intensity distribution at distance z = d and angle of rotation θ is approximated by the following expression, also known as Transport of Intensity equation [9–16]:

Where ϕ_θ (x, y) is the phase function of the object. This is valid in the near field Fresnel region, where d << a ²/λ (a is the transversal size of the smallest structure in the object and λ is the x-ray wavelength). For a mixed phase amplitude object with weak and almost homogeneous absorption, i.e. $\frac{\partial {\mu }_{\theta }}{\partial x},\frac{\partial {\mu }_{\theta }}{\partial y}\approx 0,$, Eq. (1) simplifies further to [19, 20]:

with ${\nabla }^2\equiv \frac{{\partial }^2}{\partial x^2}+\frac{{\partial }^2}{\partial y^2}.$ The goal of quantitative phase tomography is to reformulate Eq. (2) to obtain δ(x, y, z) from the knowledge of I _θ,z=d (x,y) for θ∈[0,π]. Expressing Eq. (2) as ${\nabla }^2{\varphi }_{\theta }(x,y)=-\frac{2\pi }{\mathit{\lambda d}}{g}_{\theta }(x,y)$ with $g_{\theta }(x,y)=\frac{I_{\theta ,z=d}(x,y)}{I_{\theta ,z=0}(x,y)}-1$, applying then the 3D Radon transform (denoted by the symbol ∧) and calculating finally the second derivative with respect to the variable s = xsinω + ycosω one gets:

Expression (3) is a theorem which states that from the 2D Radon transform of the measured value g, one can directly find the 3D Radon transform of δ. An explicit inversion formula for the 3D Radon transform was given by Radon and Lorentz already at the beginning of the previous century, and one therefore has a direct solution to the reconstruction problem for quantitative phase contrast CT [22]:

where the stars indicate a 2D convolution and

This convolution integral can be computed in the Fourier domain by taking the two dimensional Fourier transform. In the Fourier domain, Eq. (5) has the low-pass filter form given by [20]:

Expressed in this (filtered backprojection) form [Eq. (4)], this reconstruction approach becomes very interesting since in case of a pure phase object (I _θ,z=0 (x,y) = 1) it allows recovering the refractive index in 3D from only one single tomographic data set:

This is a significant improvement since (a) the experimental setup can be kept extremely simple (it is actually the same as for standard, absorption based tomography) and (b) the radiation damage is kept to the minimum, which is of special importance for biological specimens and for potential future in vivo studies. In addition, time consuming data alignment routines as well as phase retrieval processes become superfluous.

3. Results and analysis

In a previous work [21], we have validated the theory by implementing the aforesaid reconstruction algorithm and successfully recovered the refractive index of the simulated object. In the second step, the original algorithm [19, 20] has been used to reconstruct phase tomograms from experimental data sets, obtained at the Tomography Station of the Materials Science beamline at the Swiss Light Source [3]. The reconstructed images using one distance data set have shown artifacts due to residual absorption (see section 3.1). We have found a solution to this problem by modifying the original algorithm (section 3.2).

3.1. Experimental results

As already pointed out, the experiment setup is equivalent to the one for standard absorption tomography. Nevertheless, two additional conditions have to be satisfied:

After traversing the sample, the x- rays are converted into visible light by a 20 microns thick YAG scintillator and then magnified (10× or 20×) by an optical microscope and detected with a 14 bit, 2048 × 2048 pixel CCD camera, giving an effective pixel size of 1.4 and 0.7 microns respectively.

As a test, we investigated the internal microstructure of 350 microns thin wood samples. Even when low photon energy (10 keV) is used there is a weak absorption contrast and only the edges of internal structures are visible [see Fig. 1(a)].

In order to apply the phase tomography method in a pure phase object approximation, a single data set (I _θ,z=d (x,y)) of 501 angular projections at sample detector distance of 32 mm and 25 keV photon energy was acquired. The nonuniform illumination intensity was normalized by dividing through reference (sample-out of the beam) images. For the practical implementation of the filter in Eq. (6) a sufficiently small sampling interval has to be used in the frequency domain. This can be achieved by zero padding the discrete data of the experimental function g_θ(x,y). This is a standard technique to avoid the wraparound error in digital signal processing. Fig. 1(b) presents a slice through the sample reconstructed using Eq. (6) and Eq. (4).

 

Fig. 1. Tomographic reconstruction of a 350 microns thin wood sample: (a) absorption contrast obtained with 10 keV photon energy, (b) phase contrast, obtained with 25 keV. Effective pixel size is 0.7 microns. The length of the scale bar is 50 μm. Sample by courtesy of P.Trtik, ETHZ

Download Full Size | PDF

3.2. Analysis

The reconstructed slice shown in Fig. 1(b) is severely corrupted by the residual absorption artifact (see cap in the center). This means that, even though the absorption level for the (wood) sample is calculated to be in the range of only 2 %, pure phase object condition is not satisfied and an intensity measurement at ‘zero’ distance (Iθ,z=o (x,y)) is also required. We have adopted another approach - the reconstruction artifacts are corrected by amending the original method [20] leading to what we named Modified version of Bronnikov’s Algorithm (MBA).

The correction consists of adding, in the denominator of the low-pass filter given in Eq. (6), an absorption dependent correction factor α_exp. Consequently, the new filter has the following form:

The values of α_exp to be used are found by using a semi empirical (simulations-experiment) approach. First, we performed simulations using a mathematical phantom made up of three spheres of different size inside a cube. Two cases, a pure phase object (β=0, where β is the imaginary part of the refractive index) and mixed, weakly absorbing phase object (β of the order of 10^-9) were studied. The phase function ϕ was calculated analytically from the knowledge of δ and sets of phase contrast projection data (at different rotation angles) were generated using the Fresnel propagator [23, 24] with λ = 1 Å, sample detector distance z = 25 mm, 1852 pixels detector array and an isotropic pixel size of 1 μm. For the pure phase objects, applying Eq. (4) on one data set (g_θ(x, y) = I_θ,z=d (x, y) -1) the phase information (refractive index) has been successfully reconstructed with errors around 1%. For mixed phase-amplitude object, reconstruction using two distances data sets (g_θ(x,y) = [I_θ,z=d(x,y)/ I_θ,z=0(x,y)]-1) leads to similar result. If only a data set of angular projections at one distance is used for the reconstruction of a simulated mixed object, one gets a cap similar to the one shown in Fig. 1(b). The right correction factor in case of simulations (α_simul) has been obtained by the iterative minimization of the sum square error (SSE) defined as:

The achieved values for SSE were of the order of 10^-4. The same iterative procedure was applied for different values (in the range of 1–10 %) of absorption and the corresponding correction factors are shown in Fig. 2.

 

Fig. 2. Values for the correction factor α_simul as function of transmission found by simulations.

Download Full Size | PDF

For a specific experimental data set, the correction factor α_exp to be used in the filter of Eq. (8) is then calculated in the following way:

Average absorption/transmission level of the object is measured/calculated based on a ‘zero’ distance intensity projection. The corresponding α_simul is found from the graph depicted in Fig. 2.

Based on α_simul, the correction factor α_exp is found according to the following expression:

Where pixel_simul and array_simul are the pixel size and the detector array size (in pixels) used for simulations and pixel_exp and array_exp are the corresponding values for a specific experimental data set.

Applying the MBA method, we reconstructed the same data set, for the wood sample. The result is presented in Fig. 3. It appears clearly that the absorption artifact has completely disappeared.

 

Fig. 3. Phase tomography reconstruction (a) and the 3D rendering (b) of a 350 microns thin wood sample using modified filter given in the Eq. (8). The length of the scale bar is 50 μm.

Download Full Size | PDF

3.3. Validation

In order to validate our approach, we have compared our results with the ones obtained using an established phase contrast method, the Differential Phase Contrast (DPC) Technique [8]. To make a direct evaluation, experiments using the two approaches have been performed on the same samples and with the same detector effective pixel size (1.4 microns). The used photon energies were different, optimized for the individual techniques (17.5 keV for the DPC and 25 keV for the MBA).

 

Fig. 4. Validation of the MBA method: (a) Phase tomographic reconstruction of sample consisting of polyacrylate, starch and cross-linked rubber matrix obtained using DPC and (b) using MBA. The length of the scale bar is 100 μm.

Download Full Size | PDF

4. Discussion

We have demonstrated that the 3D distribution of the phase (refractive index) of a weakly absorbing object can be directly reconstructed from a single tomographic data set (one step). In order to reduce the residual absorption artifact the method presented in Ref. [20] has been amended leading to what we named Modified Bronnikov Algorithm (MBA). The results show that the contrast is increased, while keeping dose minimal and spatial resolution comparable to the conventional absorption based technique.

Depending on the object structure and composition, the approximation of a constant absorption over the entire object (that we used) can induce some artifacts. Additionally, the simulations, which are at the basis for the calculation of the filter correction factor, are performed without taking detector and beam characteristics into account. However, as well shown in Fig. 3(a) and Fig. 4(b), the quality of the reconstruction is largely sufficient to perform optimal segmentation and further post processing.

The MBA was validated using the Differential Phase Contrast Imaging technique [8]. For the particular example shown in the Fig. 4 the contrast obtained using MBA is even better. In general, the two methods are complementary. The full 3D phase imaging approach (MBA) is experimentally simple, very fast and it is ideally suited for small objects when resolution around 1 micron is needed. Deposited dose is equivalent (or lower, since higher photon energy is used) to the one deposited in conventional tomography. The MBA spatial resolution is limited, equivalently to standard absorption tomography, by the detector resolution to around 1 micron [3]. The DPC is more demanding in terms of instrumentation and acquisition time (therefore also deposited dose) but it is more sensitive and can be scaled up to large fields of view. Spatial resolution is limited to two period of analyzer grating [8] e.g. to 2 microns with present fabrication technology.

Future developments will focus on the fine tuning of the filtering step and on the streamlining of the whole process in order to allow high-resolution and high throughput imaging of soft biological and materials science specimens.