1. Introduction

X-ray computer tomography (CT) is a well-established technique for non-destructive three-dimensional (3D) imaging of internal structure of optically opaque samples [1–3]. Conventional CT is based on the differential attenuation of transmitted X-rays by different components of the sample. This contrast mechanism is effective for distinguishing between components with significant differences in atomic number or density, e.g. between flesh and bones in the case of medical CT. However, the difference in X-ray attenuation by different types of soft tissues (e.g. healthy and malignant ones) is rather weak, which typically results in poor image contrast. It has been suggested that X-ray phase contrast can be utilized for improvement of the contrast in transmission images of non-crystalline samples consisting predominantly of light chemical elements [4–6]. Subsequently, phase-contrast X-ray CT (PCT) has been implemented in several forms, including crystal-based and grating-based X-ray interferometry [5,7,8], analyser-based [9,10] and in-line (also called propagation-based) phase contrast [11–15] and others. The present paper investigates a particular method of object reconstruction in in-line PCT.

In-line X-ray images of most samples usually display conventional absorption contrast as well as phase contrast [4,6]. The in-line phase contrast appears as the result of free-space propagation of the transmitted beam, which transforms phase variations in the plane located immediately after the sample into measurable intensity variations in the detector plane, provided the propagation distance is large enough for the given X-ray energy, imaging configuration and sample composition and structure. For hard X-rays (with the wavelength of around 1 Å or shorter), the complex refractive index, $n(r)\equiv 1-\Delta (r)+i\beta (r)$, of most materials is very close to unity (e.g., the decrement, $\Delta (r)$, and the imaginary part, $\beta (r)$, are typically of the order of 10⁻⁶ and 10⁻⁸, respectively), hence the refraction angles are small too and it is usually possible to approximate the X-ray trajectories through a sample by straight lines. Here we also do not consider such diffraction effects as small-angle and wide-angle scattering as it is done for instance in diffraction tomography [3,16–18]. This simplifies the mathematical formalism of PCT and leads to relatively simple linear reconstruction algorithms [12,19–21]. The algorithms are particularly simple in the case of objects exhibiting negligible absorption [22] and, more generally, in the case of so-called “homogeneous” objects [21].

Recently, considerable effort has been applied to developing successful practical strategies for PCT of “general” objects, i.e. the objects that are not “homogeneous”. X-ray absorption and refractive properties of such objects are both non-trivial and “uncorrelated”. It is intuitively clear that in order to reconstruct the 3D distribution of the complex refractive index in such objects, one needs to measure 2D intensity distributions in at least two X-ray projections at each rotational position of the sample, with different values of some essential parameter, such as the sample-to-detector distance or the X-ray energy [23,24]. The main difficulty in the reconstruction of the complex refractive index distribution from such data has so far been presented by the spatial and temporal instability of X-ray beams produced by laboratory sources and synchrotron beamlines. This instability often prevents one from accurate co-registration of different images collected (at different points of time) at a fixed view angle. Recently some techniques have been suggested to overcome this difficulty [25–29], but a general practical solution is yet to be found. Other reconstruction approaches successfully applied in propagation-based PCT include the ones based on the first Born approximation, rather than on the projection approximation [15,30,31]. These methods can account for higher-order Fresnel diffraction, but require the object to be “weak” in the sense that diffracted wave is much lower in amplitude compared to the primary transmitted one.

In the present paper we investigate the options for reconstruction of a newly introduced quantity, which represents a particular combination of the real and imaginary parts of the refractive index, from X-ray PCT data. This quantity is theoretically derived from the phase-retrieval approach which was proposed by Teague [23] and further developed by Paganin and Nugent [32]. We show that in the case of objects with negligible absorption the new quantity coincides with the decrement of the refractive index of the object. In the general case, the quantity is affected by the imaginary part of the refractive index as well, with the dependence displaying a relatively simple and predictable character, as demonstrated by our numerical experiments. The advantage of the proposed approach is in the significantly simplified form of the PCT algorithm that can be applied for the reconstruction of the introduced quantity. This method can be viewed as an extension of the single-step PCT algorithm originally proposed in [22].

We present the theoretical overview of the proposed method in Section 2 of the paper, the results of the numerical simulations are found in Section 3, while the conclusions are discussed in Section 4.

2. Theoretical background

Consider an experiment with an object illuminated by a plane monochromatic X-ray wave with wavelength $\lambda$ and intensity $I_{in}$ and with the transmitted wave registered by a position-sensitive detector. The interaction of X-rays with the object is described by the spatial distribution of the complex refractive index in the object, $n(r)$, where $r=(x,y,z)$ are the usual Cartesian spatial coordinates (we omit the dependence of all quantities on the wavelength for the sake of brevity). We only consider here the case of fully coherent X-rays, however, the theory can be extended to the partially-coherent case as in [33,34].

The propagation direction of the incident X-ray wave, $I_{\text{in}}^{\text{1/2}}\exp (ik{z}^{\prime }$, where $k=2\pi /\lambda$, makes an angle ${\theta }^{\prime }$ with the z axis in the (x, y) plane, where $-\pi /2\le {\theta }^{\prime }<\pi /2$. It will be also convenient to use the angle $\theta ={\theta }^{\prime }+\pi /2$ (Fig. 1 ).

 

Fig. 1 Generic layout of the in-line CT experimental setup.

Download Full Size | PDF

We assume as usual that the projection approximation [3] can be applied to calculate the phase and intensity, respectively, of the wave after transmission through the object,

As conventional X-ray detectors can only register the intensity distribution in the transmitted beam, the transmission images collected in the plane immediately after the object can provide information only about the imaginary part of the refractive index, as follows from Eq. (2). In-line phase-contrast imaging allows one to register the phase-contrast effects caused by the real part of the refractive index, related to the phase of the transmitted beam via Eq. (1), and consequently, to reconstruct the distribution of the real part of refractive index in the sample. In order to accomplish such a reconstruction one needs to have an accurate quantitative model relating the transmitted phase distribution, ${\phi }_{\theta }(x\text{'},y)$, to the in-line intensity in the image, $I_{\theta }^R(x\text{'},y)$, acquired at a free-propagation distance R from the object. One such model is given by the finite-difference form of the Transport of Intensity equation (TIE) [23], which can be written as follows:

Teague [23] suggested to solve Eq. (4) via the introduction of an auxiliary function ${\psi }_{\theta }$, which satisfies

After the phase distribution ${\phi }_{\theta }$ is reconstructed at each projection angle $\theta$, it is possible to reconstruct the 3D distribution of the real part of refractive index by CT techniques. Taking the 2D Fourier transform of Eq. (3) and using the Fourier slice theorem [2] one can easily obtain that

Note that in the case of Eq. (2) the same approach yields the conventional (“absorption”) CT reconstruction formula:

These problems can be partially eliminated in the case of objects with negligible absorption, by means of a single-step reconstruction formula due to Bronnikov [22]:

A similar (but much stronger) stabilization of the PCT reconstruction with respect to noise in the input data is achieved in the case of so-called “homogeneous” objects [21,36]. The latter class of objects includes objects consisting predominantly of a single material, and, in the case of X-rays with energies between approximately 60 keV and 500 keV, any objects containing only light chemical elements with Z<10 [37].

Using Eq. (13) and Eq. (10) one can easily verify the following identity:

Using Eq. (14), the single-step reconstruction algorithm, Eq. (13), can be re-written in the following form:

Unfortunately, in the cases with non-trivial absorption the approach which led to Eq. (15) is usually no longer valid, as, unlike the Laplace operator, the general TIE operator, $D_{\text{TIE}}[{\phi }_{\theta }](x^{\prime },y)=-{\nabla }_{\perp }•[I_{\theta }(x^{\prime },y){\nabla }_{\perp }{\phi }_{\theta }(x^{\prime },y)]$, does not commute with the Radon transform. This is why we propose the following approach based on the Teague's method for solution of the TIE described above.

Let us assume that there exists a function $\tilde{\Delta }(x,y,z)$ such that

By an exact analogy with the derivation of Eq. (15) from Eq. (13), one can now also derive the following reconstruction formula from Eq. (17):

An important question that needs to be answered with respect to quantity $\tilde{\Delta }(x,y,z)$ is that about its physical nature. Firstly, we note that in the case of objects with negligible absorption, i.e. when it is possible to use the approximation $I_{\theta }(x\text{'},y)\cong I_{\text{in}}$, Eq. (17) reduces to Eq. (13), and hence $\tilde{\Delta }(x,y,z)\equiv \Delta (x,y,z)$ in this case. In the general case, when absorption inside the object cannot be ignored, the quantity $\tilde{\Delta }(x,y,z)$ depends both on the real and imaginary parts of the complex refractive index. The following formal expression can be derived from Eqs. (17), (6), (4), (1) and (2):

We carried out computer experiments with numerically simulated objects (phantoms) that generate different amounts of X-ray absorption and refraction at a given X-ray energy in order to investigate the behavior of the quantity $\tilde{\Delta }(x,y,z)$ as a function of the real and imaginary parts of the refractive index of the phantoms. The results of the simulations presented in the next section of the paper shed some light on the nature of physical information about the object that can be obtained by means of in-line phase-contrast CT experiments combined with the simplified reconstruction algorithm described by Eq. (17) or Eq. (18). The simple formula (Eq. (22)) connecting $\tilde{\Delta }$, $\Delta$ and $\beta$ for certain types of objects was also verified by our simulations.

3. Methods and results

The simulated object, as depicted in Fig. 2 , was a collection of 12 non-overlapping homogeneous spheres of radius 100 µm, one in each of the eight corners of the bounding cube and one on each face excluding the top and bottom. It was represented as two 1024 × 1024 × 1024 arrays of 32-bit floating-point numbers, corresponding to the $\beta$ and $\Delta$ distributions of the complex refractive index $n=1-\Delta +i\beta$, which was different for each of the spheres. Both arrays were stored and processed as 1024 square planar cross-sections parallel to the x-z plane. The voxel size was 1 × 1 × 1 $\mu m^3$.

 

Fig. 2 The structure of the object—and its bounding cube—that was simulated. It consisted of twelve spheres arranged in three layers of four, each parallel to the base. The bottom and top layers had the spheres situated in the corners, and the middle layer had them in the centre of each face. The lightness of the spheres decreases with increasing $\beta$ for the constant $\Delta$ objects, $\Delta$for the constant $\beta$ objects, and $\beta$ and $\Delta$ for the objects with both varying.

Download Full Size | PDF

The simulations were performed using X-TRACT [39], a program that amongst other things allows one to calculate parallel-beam X-ray projections of objects given the object's $\beta$ and $\Delta$ slices. X-TRACT allows one to perform a variety of useful image operations such as pixel-wise addition and inverse Laplacian, to process image files into sinograms, and to perform CT reconstructions by calculating the inverse Radon transform of sinogram images. All of these operations were required for the simulations.

First, X-TRACT was used to generate parallel-beam projections of the object at two different sample-to-detector propagation distances, namely 0 cm and 10 cm. These projections were simulated for the X-ray wavelength of 1 Å. A Gaussian filter with the width 2$\sigma$ = 4 pixels was applied to each stack of slices before propagating. This models roughness at the interfaces of the object and the implicit convolutions with point-spread functions of the source and the detector. In total, 180 of these projections were generated, with each projection corresponding to a one degree rotation of the object from the previous position.

The derivative of the intensity in the direction of propagation, ${\partial }_z{I}_{\theta }(x\text{'},y)$, was then approximated by performing a pixelwise subtraction of the 10 cm projection from the 0 cm projection and then dividing by the propagation distance, i.e 10 cm. Next an inverse 2D Laplacian operation was performed on each derivative image, and because the boundary conditions make the inverse Laplacian unique only up to an additive constant, this constant was approximated—by averaging the background values—and subtracted in order to force the background values to be as close to zero as possible. The images were also multiplied by the quantity $-k/R$ in order to generate the 2D distributions of the function ${\psi }_{\theta }(x^{\prime },y)$ in accordance with Eq. (6). The obtained images were then used to create a series of sinograms from “projections” ${\psi }_{\theta }(x^{\prime },y)$, and the CT reconstruction was performed to these sinograms in the same way as normally applied for the solution of Eq. (1). This process yielded the 3D distribution of the quantity $\tilde{\Delta }(x,y,z)$ in accordance with Eq. (16).

In order to explore the relationship between the complex refractive index $n=1-\Delta +i\beta$and $\tilde{\Delta }$ with regards to CT reconstruction, six different objects were used: three with constant $\Delta$ and three with constant $\beta$. Then, to more closely model experimental results, the constant $\beta$ simulations were performed again with Poisson noise being added to the projections, with a relative standard deviation, ${\sigma }_{rel}^{}$, of 5% calculated for the average intensity.

The three constant $\beta$ objects had $\beta$ values of 1.5 × 10⁻⁹, 1 × 10⁻⁸, and 2 × 10⁻⁸, each with a Δ range of 5 × 10⁻⁸ to 1 × 10⁻⁶. The three constant $\Delta$ objects had $\Delta$ values of 5 × 10⁻⁷, 1 × 10⁻⁶, and 2 × 10⁻⁶, each with a β range of 1 × 10⁻⁹ to 2 × 10⁻⁸. These values were chosen close to refractive indices of certain real materials; at a wavelength of 1 Å [40] are shown in Table 1 . The results are shown in Figs. 3 , 4 and 5 , and all confirmed the theoretical prediction that $\tilde{\Delta }$ would approach Δ as $\beta$ approached 0.

Table 1. The Components of Refractive Indices of Certain Real Materials at a Wavelength of 1 Å [35].

View Table | View all tables in this article

 

Fig. 3 The reconstructed values using $\tilde{\Delta }$ for the three objects with constant $\beta$ (closed circles, squares and diamond symbols). The dashed line is $\tilde{\Delta }=\Delta$. The uncertainty is taken to be 1 standard deviation of the values in the centre of the spheres in the reconstructed slices. Also shown are the curves of best fit $\tilde{\Delta }=\Delta \exp \{-2k\frac{<P\beta >}{\beta }\beta \}$ for closed circles, squares and diamond symbols, which parameters are shown in Table 2. The triangular symbols show the result of a “homogenous” TIE reconstruction.

Download Full Size | PDF

 

Fig. 4 The reconstructed values using $\tilde{\Delta }$ for the three objects with constant β and 5% Poisson noise. The dashed line is $\tilde{\Delta }=\Delta$. The uncertainty is taken to be 1 standard deviation of the values in the centre of the spheres in the reconstructed slices. Also shown are the same curves of best fit as in Fig. 3.

Download Full Size | PDF

 

Fig. 5 $\tilde{\Delta }$ (triangular, open circles and inverted triangular symbols) and the reconstructed values of $\Delta$ (see Eq. (22)) (closed circles, diamond and square symbols) using $\tilde{\Delta }$ for the three objects with constant Δ. The uncertainty is taken to be 1 standard deviation of the values in the centre of the spheres in the reconstructed slices. Also shown are the curves of best fit $\tilde{\Delta }=\Delta \exp \{-2k\frac{<P\beta >}{\beta }\beta \}$, which parameters are given in Table 3.

Download Full Size | PDF

Figure 3 suggests that $\tilde{\Delta }$ can be approximated by a linear function $\tilde{\Delta }=m\Delta$, with $0<m<1$, in agreement with Eq. (22) for constant $\beta$. This is also supported by Fig. 5, and both figures show the curves of best fit for each simulation's data series. Tables 2 and 3 show the values of $\frac{<P\beta >}{\beta }$ corresponding to each curve of best fit. Figure 5 also demonstrates that Eq. (22) may allow one to accurately obtain the function $\Delta (x,y,z)$ from the reconstructed functions $\tilde{\Delta }(x,y,z)$ and $\beta (x,y,z)$ in the cases with significant spatial variation of both $\beta$ and $\Delta /\beta$. It should be noted that while the slow variation of $(P_{\theta }\beta )(x\text{'},y)$ was demonstrated to be sufficient for the validity of Eq. (22), it is not a necessary condition. As our numerical results show, Eq. (22) remains a good approximation for a broader range of objects than that satisfying the condition of slow variation, including certain types of objects with localized areas of rapid variation of $(P_{\theta }\beta )(x\text{'},y)$.

Table 2. The Best-Fit Parameters Obtained from Figs. 3 and 4.

View Table | View all tables in this article

Table 3. The Best-Fit Parameters Obtained from Fig. 5.

View Table | View all tables in this article

For comparison, we have also performed the reconstruction of the real decrement of the refractive index $\Delta (x,y,z)$ in the sample with constant $\beta =1.5\times {10}^{-9}$ using the PCT reconstruction formulae for homogeneous objects [21] with $\Delta =5\times {10}^{-7}$ (which approximately corresponded to an average Δ value for this sample) . The appropriate results are shown in Fig. 3 by green triangular symbols. This data demonstrates that the “homogeneous” TIE method, although producing reasonable looking CT reconstructions, cannot quantitatively reconstruct with satisfying accuracy the distribution of the real decrement of the refractive index in objects with significant spatial variability of the $\Delta /\beta$ ratio.

Figure 4 shows the constant $\beta$ reconstructions with 5% noise. Whilst being visibly more unstable than their noiseless counterparts, the quantity $\tilde{\Delta }(x,y,z)$ is still reconstructed with remarkable stability and closely followed the same trend lines.

The graphs of Fig. 5, whilst being largely linear, are all marred by a small but systematic 'jump' between the fourth and fifth data points. This jump is likely caused by the geometry of the simulated object, as the four spheres in the middle layer are closer together than those in the top and bottom layers.

Figure 6 shows the projection at $\theta ={26}^{\circ }$ of the object with $\Delta =5\times {10}^{-7}$ at an object-to-detector propagation distance of 10 cm, along with a reconstructed slice of $\Delta$ through the bottom layer of spheres. Visible are the bright diffraction fringes encircling each sphere (Fig. 6(a)), as typical for the in-line phase contrast images.

 

Fig. 6 (a) Typical projection with an object-to-detector projection distance of 10 cm; (b) reconstructed slice of $\Delta$ through the bottom layer of spheres and (c) a cross-section indicated by the solid line through the bottom row in (a). This particular object had spheres with constant $\Delta =5\times {10}^{-7}$ and $\beta$ varying from 1 × 10⁻⁹ to 2 × 10⁻⁸.

Download Full Size | PDF

4. Conclusion

We presented a variant of in-line phase-contrast tomography based on the Teague’s method. Our newly introduced quantity $\tilde{\Delta }(x,y,z)$ is uniquely related to the complex refractive index in the object and can be reconstructed using two projections at different propagation distances collected at each view angle. The proposed algorithm for the reconstruction of $\tilde{\Delta }$ from the experimental X-ray projection data is significantly simpler compared to that of the real decrement of the refractive index, $\Delta$. We also showed that, similarly to the case of non-absorbing objects, but unlike the case of $\Delta$ in objects with significant and non-uniform X-ray absorption, the quantity $\tilde{\Delta }$ can always be reconstructed by means of inverse 3D Laplacian applied to the difference between the “conventional” CT reconstructions obtained separately from the two projection data sets collected at different object-to-detector distances. We argue that this approach may present certain advantages compared to the more conventional PCT approach that requires phase retrieval to be applied at each projection angle separately. Furthermore, knowing the distributions of $\tilde{\Delta }$ and $\beta$ in the object, may allow one to find the real decrement of the refractive index, $\Delta$of materials in the object. Namely, we showed that for a certain class of objects the connection between$\tilde{\Delta }$, $\Delta$ and $\beta$ can be approximated by a simple expression, $\Delta (x,y,z)\cong \tilde{\Delta }(x,y,z)\exp \{2k<P\beta >(y)\}$, where $<P\beta >(y)$ is the projection $(P_{\theta }\beta )(x^{\prime },y)$ averaged over the coordinates $x^{\prime }$ and $\theta$. Our computer simulations have shown that this significantly simplified reconstruction procedure is stable with respect to noise. We intend to conduct suitable experimental tests of the proposed method in the near future. The outcomes of such experimental tests will be presented in a subsequent publication.