1. Introduction

X-ray phase contrast imaging (PCI) consists in a number of methods that allow to detect and quantify changes in the X-ray wave phase after its interaction with an object (for a brief overview and a comprehensive review see [1, 2] correspondingly). Most of the PCI experiments are done in such a way that the intensity of phase-contrast images depends either on the Laplacian of the wave phase [3] or on its first partial derivative [4–6]. The X-ray wave phase in the object exit plane can be computed from phase contrast images. This procedure is often called phase retrieval. In this work we consider the problem of the phase retrieval from X-ray differential phase contrast (DPC) images, which contain information about the first derivative of the phase over one spatial coordinate (so called unidirectional DPC images [7]). These data are commonly acquired in experiments with analyzer crystals [8], grating interferometry setups [5, 9], coded aperture technique [6], and speckle-based approaches [10].

Phase retrieval by means of a straightforward one-dimensional integration produces strong streak artifacts due to inevitable presence of noise in experimental data. This problem was addressed in several works [7, 8, 11], where the integration problem was converted to a minimization problem and regularization terms were introduced in order to suppress the streak artifacts. The general idea behind the regularization is to enforce a coupling between spatial coordinates, which physically reflects an integrity of the object under investigation. The smoothness and continuity of the desired solution is controlled by a stabilizing (regularizing) functional. Wernick et al. used the L2 norm of the second derivative as a stabilizing functional, which essentially corresponds to a priory assumption that a linear function is a good approximation of the phase continuity in the direction perpendicular to the direction of integration [8]. The solution can be found analytically and a very fast computation method, based on fast Fourier transform, can be applied. Thuering et al. and, more recently, Sperl et al. suggested to use L1 and L2 norms of the first derivative as a stabilizing functional. Such regularization approaches are well-known in image processing. It is considered to have the advantage to preserve the spatial resolution, as physically this regularization term implies that the sought function can be represented as a piecewise constant function [12]. The solution of the resulting minimization problem can not be obtained directly and requires the application of iterative minimization methods. The computation time for a large data set might be orders of magnitude higher than in case of direct inversion [8]. In the presence of strong noise and for large data sets, advanced computation methods should be used to find a reliable solution in a reasonable time [11].

It was also shown that the basic phase retrieval problem can be converted to an elliptic partial differential equation [13, 14]. In this work we explore this approach in details. First, the suppression of streak artifacts is verified on a numerical phantom similar to those used in other works [7, 11]. Second, the accuracy of the phase retrieval is tested using a simple experimental phantom of known composition. A possible practical application is also demonstrated using a biomedical sample with a complex internal structure. For this purpose, a tomographic set of X-ray DPC images was processed by the suggested phase retrieval algorithm. The obtained phase projections were used to perform CT reconstruction of the index of refraction inside the object. This reconstruction yields absolute magnitudes of the refractive index decrement, while typically a reference material has been used in order to normalize the reconstructed data and to obtain such quantitative information [15, 16]. Finally, the differential equation for the phase retrieval is equivalent to the heat conduction one. There exists a number of software packages optimized to solve heat transfer problems, which means that an existing software can be adapted to perform the phase retrieval. In the end of this work we present a brief description of how the phase retrieval can be performed with the aid of the ANSYS^® software [17].

2. Derivation of the equation

Figure 1 shows a schematic of a PCI experiment and of the coordinate notations in the general case of CT imaging. Here we assumed that phase-contrast imaging can be done using various approaches. The image intensity depends on both the attenuation and the refraction (i.e. angular deflection) of the X-ray beam induced by the interaction with an object. Depending on the used experimental technique, the phase-stepping [5], diffraction enhanced imaging [18], or other algorithms [6, 10] are used to calculate deflection angles of X-rays in the object exit plane. Two cases are considered: (1) the object rotation axis is perpendicular to the phase contrast sensitivity axis, so that the differential phase signal depend on the amount of X-ray deflection in the CT reconstruction plane (in-plane geometry); (2) the object rotation axis is parallel to the phase contrast sensitivity axis, so that the differential phase signal is measured in the plane orthogonal to the CT reconstruction plane (out-of-plane geometry). Following derivations are done for the former case. The latter case is related to the results presented in Section 4.1.

 

Fig. 1 Sketch of the PCI experiment with the definition of the coordinate system and parameters used in derivations. Two possible PCI experiments are distinguished depending on the orientation of the object rotation axis with respect to the phase-contrast sensitivity axis: panel (a) shows the in-plane geometry, panel (b) is the out-of-plane acquisition geometry.

Download Full Size | PDF

Relations between deflection angles α, wave phase ϕ and distribution of the refractive index decrement δ inside the object can be derived using the geometrical optics approximation and the paraxial ray equation:

Equation (2) relates deflection angles to the Radon transform of the refractive index along the ray path. The phase delays induced by the object can be expressed as:

It is seen, that the X-ray angular deflection at a point in the object exit plane is linearly proportional to the local spatial derivative of the wave phase.

Phase retrieval can be considered as the 1D Cauchy problem [19], providing that phase delays are zero at the reference edge of the image, so that the initial condition for Eq. (4) can be defined. However, the right hand side of the Eq. (4) is the quantity measured in PCI experiments and it always contains noise. If the direct one-dimensional integration is performed the extracted phase images will be severely distorted by streak artifacts [see Fig. 2(c)].

 

Fig. 2 Definition of the phase and its derivative at object boundaries. Case (a) differs from case (b) in direction of the phase contrast sensitivity axis. It is assumed that the object (the gray area) exceeds the field of view (which is shown as a black rectangle) in the vertical direction. For the sake of compactness (x,y) coordinates of lines that define object boundaries are indicated with subscripts left, right, top, bottom. The Dirichlet boundary condition can be then set as following: ϕ_left = ϕ_right = 0, since in the horizontal direction a part of the unperturbed wavefront is taken in the image; the phase can be constrained at the top and bottom boundaries assuming that object is homogeneous (ϕ : hom. obj. approx.) and its thickness is approximately known. Measured DPC values can be used in both cases as the Neumann condition for the phase derivative along the axis parallel to the DPC axis (∂_x,yϕ : exp.). The phase derivative in the orthogonal direction is unknown, as, for instance, ∂_yϕ at the top and bottom edges in case (a). The periodic boundary condition can be used in this case. In case (b) ∂_xϕ is unknown, but a reasonable approximation (∂_xϕ : approx.) can be made in order to reduce the artifacts.

Download Full Size | PDF

Only a small fraction of a priori information can be utilized if the phase is reconstructed using the 1D Cauchy problem. In practice more information about an object is typically available and can be used to reduce reconstruction artifacts. For instance, the sample is often partially isolated so that the phase should be zero at all sample boundaries with air, not only at a single reference line. Next, it is common to embed samples (in particular biomedical ones) in plastic cylindrical containers. An assumption can be made about phase gradients at the edges of the container and about the phase in image regions where the sample composition and geometry are known. In principle solution can be constrained at all points, where the phase and(or) its gradient are known or can be estimated analytically. This a priori information can be exploited if the phase retrieval is transformed to a boundary value problem, which requires a larger set of additional constraints than an initial value problem.

Equation (4) can be converted to a boundary value problem in the following way. Taking the first partial derivative of both sides along the y-coordinate gives the 1D second order partial differential equation:

The weighted regularization term can be then added to the left-hand side in order to enforce the smoothness of the retrieved phase in the direction perpendicular to the DPC axis:

The obtained expression is a 2D elliptic partial differential equation similar to those describing 2D heat transfer problems. In order to solve Eq. (6) the boundary condition should be defined.

2.1. Phase-retrieval boundary value problem and a priori information

Boundary condition of the first-type (Dirichlet) is required in order to find the unique solution of Eq. (6) and to retrieve the absolute value of the phase delay. However, a simple homogeneous Dirichlet condition of zero phase delay at image boundaries (suggested, for instance, in [13]) is a rarely encountered condition in practice. It is often impossible to isolate an object in the center of the camera field of view so that it is surrounded by air, where X-ray phase remains unperturbed. In an experiment, for instance in CT imaging, object intersects at least one of image boundaries. At this edge the phase delays can be as large as hundreds of radians and the reconstruction under assumption of zero phase can produce erroneous results. Although exact values of the phase delay at the reconstruction region boundaries are unknown, a reasonable approximation to the expected phase delays will greatly reduce the reconstruction artifacts.

Large difference in the index of refraction between a material and the air is also the typical source of artifacts. These errors can be reduced if a priory information about the phase derivative (the second-type or the Neumann boundary condition) is added in the reconstruction problem. The phase derivative in the phase-contrast sensitivity direction is readily available as it is measured in the experiment. In the orthogonal direction a reasonable assumption about the expected phase derivative can be made to reduce the reconstruction artifacts. Generally, the magnitude of the phase gradient at the object edge will determine the rate, at which phase delay changes near the object edge. The value of the phase derivative can be adjusted according to the object shape near the edge. For instance, a very large value of the phase gradient can approximate a sharp discontinuity at the object edge.

Generally the amount of a priori information about the sample allows to define the inhomogeneous mixed boundary conditions and in this way to constrain both the phase and its derivative. Figure 2 illustrates how boundary conditions were handled in this work. A case frequently encountered in PCI CT experiments is considered when sample dimensions exceed the field of view in one direction while in the other direction both object edges fit inside the DPC image. Without loss of generality it can be assumed that the object intersects the top and bottom image boundaries. Dirichlet boundary condition was defined as follows: (1) zero phase delay along the interfaces with air was set; (2) at the top and bottom boundaries the phase delay was approximated assuming a homogeneous object consisting of a single material. Our work concerned with biomedical samples and we estimated the phase at boundaries assuming a material equivalent to water. The projected object thickness can be estimated using, for instance, complementary absorption data typically obtained in PCI CT imaging. The Neumann boundary condition can be readily defined along the axis, which coincides with the phase contrast sensitivity axis. In the orthogonal direction, the phase derivative can be approximated either by zero, the periodic boundary condition, or a finite value. Several trials might be necessary in order to find Neumann boundary condition which minimizes artifacts. Note that object edges, at which boundary conditions are defined, should not necessarily coincide with the coordinate lines.

Finite element method was used to retrieve the phase using Eq. (6) and inhomogeneous mixed boundary conditions. Finite element discretization of the boundary value problem produces a linear system of equations, which has a unique solution [20]. Although the rank of this system of linear equations is large (the number of equations is equal to the number of pixels in the reconstruction area), the left-hand side coefficient matrix is very sparse and a solution can be obtained extremely quickly. The addition of the a priori information and the corresponding reduction of the number of unknowns and rearranging of the linear system of equations (which is sometimes called reduction) is the process that requires optimization. In our own implementation of the finite element method in C++ it is the most time consuming part of the computation. We have seen that a commercial software can set the boundary conditions and yield the solution of Eq. (6) in a matter of seconds at regular desktop computers even for a very large image. It is worth noting that we have also tried to solve the phase retrieval boundary value problem using the finite difference method, however the resulting coefficient matrix appeared to be ill-conditioned and it was often impossible to retrieve the phase.

3. Tests with numerical and experimental phantoms

3.1. Numerical data

A numerical phantom was used both for the general assessment of the method and for studying the dependence of the solution on the magnitude of the regularization parameter γ. For the sake of comparison with results of [7, 11], the Shepp-Logan phantom was generated on a 256 × 256 grid. It was assumed that the generated projection shown in Fig. 3(a) represents the line integral of the index of refraction decrement δ. The DPC image was calculated according to Eq. (2) by a simple forward difference approximation to the 1D derivative operator. It is difficult to derive the exact noise model, since several mathematical operations are always performed to calculate the DPC image from a sequence of raw phase-contrast images. Therefore, we added a Gaussian noise with the root mean square value, which is equal to 5% of the maximum intensity value in the generated DPC image. The appearance of the resulting noisy DPC image (Fig. 3(b)) is similar to the typical data obtained in our experiments, where about 1000 photons contributed to the signal in each image pixel within 100 msec exposure time. Note that neither low-pass filtering (smoothing) was applied, nor values of bony tissues were decreased when modeling the DPC signal.

 

Fig. 3 Tests with the Shepp-Logan phantom: (a) ground truth image; (b) simulated noisy DPC image; (c) direct 1D integration of DPC image displays characteristic stripe artifacts. Phase retrieved using Eq. (6) with (d) γ = 0.001, own software; (e) γ = 0.05, own software; (f) γ = 0.05, ANSYS software. Insets (g) and (h) show 1D phase profiles over row 118 and column 118 respectively. Profiles were extracted from images (a,c,e) and there location is indicated by dashed red lines in Fig. 3(e).

Download Full Size | PDF

Typical streak artifacts arising after 1D integration are seen in Fig. 3(c). Figures 3(d) and 3(e) present the phase projection obtained using the suggested Eq. (6) with γ = 0.001 and γ = 0.05 correspondingly. For γ = 0.05 the result is also presented that was obtained with the ANSYS software, Fig. 3(f). It can be seen that when γ is increased the residual stripe pattern is vanishing. At large γ the appearance of another artifact can be noticed in the vicinity of the upper and lower vertices of the outer ellipse. These areas are the most complicated as the phase-contrast sensitivity direction (along the image rows in this case) coincides there with the surface tangent. This error is somewhat amplified in the described case due to the small dimensions of the computation grid and, as a consequence, lesser accuracy of the finite difference approximation to the derivative operator. Nevertheless, this effect will always appear in practice, in particular at interfaces between materials, where the refractive index experiences the largest fluctuations (f.i. bone-soft tissue interface).

We quantified the root mean square error (RMSE) and contrast-to-noise ratio (CNR) metrics of the reconstructed phase images. The RMSE of two images was calculated as the ratio of the L2 norm of the residual image to the square root of the total number of pixels in the image. The CNR factor was found by selecting two regions (object and background, same as in [7]) and calculating the double ratio of the difference between mean values in these two regions to the sum of standard deviations in them. Results are gathered in the Table 1.

Table 1. Root mean square error and contrast to noise ration for several γ values.

View Table | View all tables in this article

It can be seen, that the CNR continues to increase as γ grows, however the RMSE metric increases much faster due to excessive smoothing introduced by the stabilizing term in Eq. (6). We were not able to reach high CNR values of ∼40, that were achieved in the work of Thuering et al., however reasonably good image metrics were obtained with our method at γ ≈ 0.05. The good agreement between extracted and ground truth values is also illustrated in Figs. 3(g) and 3(h) where 1D phase profiles are plotted. Streak artifacts are largely suppressed while excessive blurring is avoided at sharp edges. For this reason we always selected the regularization multiplier in the range 0.01–0.05 when reconstructing the experimental data.

3.2. Experimental data

Data was acquired with the analyzer based imaging setup located at the biomedical beamline ID17 of the European Synchrotron Radiation Facility (France). The Bragg reflection from [333] planes of a perfect Si crystal was used to prepare the monochromatic beam with a very small divergence on the sample’s entrance plane. After the wave passes through the sample, the phase-contrast is generated by reflecting the wave from another Si crystal, identical to the first one. Detailed specifications of the imaging setup and acquisition procedure can be found in [15]. Each differential phase contrast projection was calculated from a pair of raw phase-contrast images by using a non-linear extension for the diffraction enhanced imaging algorithm [21]. All experiments were done at 51 keV photon energy.

In order to investigate the accuracy of the phase retrieval algorithm, the last was applied to the DPC image of a simple plastic phantom. It is composed of two coaxial cylinders. The smaller polyethylene (PE) cylinder with diameter of 6 mm was tightly inserted in a hole drilled in a 15 mm-diameter Plexiglas (PMMA) cylinder. Inside the sample there is also a small air cavity remained after the drill lip (see Fig. 4(a), triangular shaped object at the bottom of the PE cylinder). Images of the phantom were obtained with a detector pixel size of 96 × 96 μm².

 

Fig. 4 Experimental images of a phantom: (a) DPC projection, (b) phase projection obtained with Eq. 6, (c) 1D phase profiles taken over line 280 (indicated with the dashed red line in panel (b)) together the expected (theoretical) phase profile.

Download Full Size | PDF

Figure 4(b) shows the phase image retrieved from the DPC image using Eq. (6). The inner PE cylinder, which is clearly visible owing to the edge enhancement in the DPC image, has a very small contrast in the phase image, since the refractive index decrements of both materials have a similar magnitude. Still, the inner cylinder can be recognized in the 1D phase profile (Fig. 4(c), solid red line). Using Eq. (3) the experimental phase profile can be fitted very well (Fig. 4(c), dotted line) with the values of the refractive index decrement for PMMA and PE materials taken from the Evaluated Photon Data Library (EPDL97) tables: δ_PMMA = 10.21 × 10⁻⁸, δ_PE = 8.44 × 10⁻⁸. The maximum experimental phase delay is in a good agreement with the tabulated values. However the deviation from the expected phase profile becomes noticeable as the distance from the center of the cylinder increases. The discrepancy between experimental and expected data is most severe near the PMMA-air interface. This can be explained by a very large difference (a discontinuity) between the index of refraction of the object and air. The magnitude of the refraction angles can significantly exceed the width of the rocking curve in these regions. This gives rise to the spatially localized systematic error in the intensity of analyzer based images. This error can be suppressed by immersing the sample in a water tank. Also, the profile of the real experimental rocking curve slightly deviates from a Gaussian profile, but this fact is not taken into account in the algorithm that we used to calculate DPC images. Finally, it was assumed that the incident X-ray beam is a perfect plane wave, which was not the case in reality. This was not taken into account either in our calculations. In order to eliminate these errors, further optimization of the imaging setup and reconstruction algorithms is required.

4. Refraction-based CT using phase projections obtained from DPC images

In this section we demonstrate a possible practical application of our phase retrieval method combined with refraction-based X-ray CT imaging. The experiment was made using a rabbit knee joint sample. The expected benefit of imaging the refractive index inside such sample is the simultaneous visualization of the bone and the soft tissue surrounding the bone. This is very important for the study of the osteoarthritis disease: a good contrast of both bony and soft tissues can not be achieved neither with MRI diagnostic nor with absorption CT imaging [22]. To this end and to further validate the potential of the suggested phase retrieval algorithm, we applied it to a tomographic data set of DPC images of the rabbit knee sample. We then used the retrieved phase projections to reconstruct the index of refraction inside the sample.

The rabbit’s knee joint was extrarticulated and placed in a cylindrical plastic container with a inner diameter or 50 mm. The residual volume in the container was filled in with the solidified agarose gel. The thickness of the container wall was about 1 mm. A tomographic set of 500 DPC projections was obtained under the conditions described in Sec. 3.2. The only difference was the detector pixel size, which was changed to 46 × 46 μm² in order to achieve higher spatial resolution. The size of the projections (computation grid) processed by the phase retrieval algorithm was 1189 × 499 pixels.

Figure 5(a) shows a DPC projection of the rabbit knee joint. Despite the large maximum projected thickness of the sample, the DPC signal can be clearly distinguished at the interfaces between different tissues and details. On the other hand, very few details can be seen in the retrieved phase projection presented in Fig. 5(e). The maximum phase delay in the central part of the sample (where the total object thickness equals to 50 mm) exceeds 1000 radians. That is why fine inner features can be barely seen in the phase projection image. More details appear if the contrast enhancement is applied to a fragment of the projection, as shown in Fig. 5(f).

 

Fig. 5 Reconstruction of the refractive index distribution inside the rabbit knee leg. Results obtainable with the presented technique (shown in the bottom row) are compared with data reconstructed using FBP algorithm for gradient projections (FBP with Hilbert filter [9]). Images present the following results: one of the DPC projections obtained in the experiment (a), the colorbar on its left-hand side is proportional to the X-ray deflection angles in radians; axial and sagittal images of refractive index distribution obtained using the FBP algorithm for gradient projections (c) and (d); the phase projection (e) retrieved from the DCP image (a), colorbar on its left-hand side is proportional to the X-ray phase delay in the sample exit plane in radians; axial and sagittal views obtained using FBP algorithm for ordinary projections (g) and (h), the colorbar on the right-hand side shows the magnitude of the refractive index decrement. Insets (b) and (f) show magnified fragments of images (a) and (e) with automatic contrast adjustment made in ImageJ software.

Download Full Size | PDF

The internal composition of the sample is clearly depicted when 500 phase projections are used for CT reconstruction. Simplest filtered backprojection algorithm with Ram-Lak filter (FBP, [23]) was applied. Images of the refractive index in both axial [Fig. 5(g)] and sagittal [Fig. 5(h)] planes allow to distinguish various features inside the rabbit knee leg. Arrows and numbers in Fig. 5(g) indicate: (1) the cartilage next to the bone, (2) a feature inside the muscle tissue, (3) stitches, remained after surgical intervention, and (4) a fatty tissue. These structures can be also seen in images shown in Figs. 5(c) and 5(d), which were obtained by means of the FBP algorithm for gradient projections. Intensity of these images is proportional to the local deviation of the refractive index from its mean value [24]. In order to retrieve the absolute quantities further normalization or integration should be done.

On the other hand images obtained using integrated DPC projections are proportional to the absolute value of the index of refraction decrement δ. Table 2 contains values of δ for several materials and tissues inside the rabbit leg. Although theoretical data are not available for the studied case, reference values of the corresponding human tissues derived from tables [25] were indicated if possible for the sake of comparison. The refractive index decrement δ is linearly proportional to the physical density of materials at X-ray energies well above the K-edges of the materials composing the sample. Therefore presented δ values can give an idea about the relative density of the different tissue and materials [15, 16]. For instance, we found that the muscle tissue in the rabbit leg was on average slightly denser than the human one: δ_rabbit = (9.2±0.1)×10⁻⁸ against δ_human = 8.9×10⁻⁸. On the other hand the density of the rabbit bone is less than the typical density of the human cortical bone: δ_rabbit = (13.5±0.1)×10⁻⁸ against δ_human = 15.0 × 10⁻⁸. The measured average density of the material filling the spongy inner part of the rabbit bone is also smaller: δ_rabbit = (7.8±0.1)×10⁻⁸ against δ_human = 8.1×10⁻⁸. The solidified agarose appeared to be rather inhomogeneous. This can be easily explained by the fact that during the final stage of the solidification process it was not possible to agitate the gel. The agarose stratification is also visible at the bottom of Fig. 5(a).

Table 2. Measured and expected index of refraction decrement of materials and tissues inside the biological sample.

View Table | View all tables in this article

4.1. Refractive index CT in imaging systems with rotation axis parallel to the direction of phase sensitivity

Suggested technique was further tested in the refractive index CT reconstruction using the DPC images acquired in the setup, in which the rotation axis is parallel to the plane in which DPC is measured. In this case DPC signal depends on the amount of the X-ray deflection in the direction orthogonal to the CT reconstruction plane [the out-of-plane geometry shown in Fig. 1(b)]. The common reconstruction approach based on the FBP with a specialized gradient filter function (Hilbert filter [9]) can not be used. We show that the reconstruction can be performed in two steps: first the phase delay projections are obtained from DPC images and then an ordinary FBP algorithm is used for CT.

The experiment was performed using a human breast sample embedded in a PMMA container. The effective pixel size was 96 × 96 μm², the photon energy was set to 51 keV. First, DPC projections were acquired in the in-plane geometry. Then the object rotation axis was aligned along the phase contrast sensitivity plane and the data acquisition was repeated. Phase images were retrieved using mixed boundary conditions shown in Figs. 2(a) and 2(b) for the in-plane and out-of-plane setups accordingly. The regularization term used to suppress the stripe artifacts was set to γ = 0.05 for the in-plane and γ = 0.02 for the out-of-plane geometry. FBP algorithm with the Ram-Lak filter was used for tomographic reconstruction.

Examples of the index of refraction distribution reconstructed in the axial and sagittal planes are shown in Fig. 6. Reconstruction artifacts are clearly visible in the left upper quadrant of the axial distribution reconstructed in the out-of-plane geometry [Fig. 6(b)]. They are localized and do not spoil the results completely. We attribute these artifacts to the presence of air bubbles at the top of the sample in the region where the boundary condition were defined. Sagittal view taken in a plane, which does not coincides with these artifacts [Fig. 6(d)], depicts the inner structures even better than the corresponding image obtained in the in-plane geometry [Fig. 6(c)]. This can be expected considering that the direction of sensitivity in the out-of-plane geometry is parallel to the sagittal plane.

 

Fig. 6 Reconstructed distributions of the index of refraction in the breast sample. Axial and sagittal views obtained in the in-plane geometry are shown in panels (a,c) correspondingly. Panels (b,d) show axial and sagittal views reconstructed in the out-of-plane geometry. Dashed red lines in Figs. 6(a) and (b) show the location of sagittal sections. The right hand side and left hand side gray bars show tissue density in g/mm³ (derived from δ values) for images (a,c) and (b,d) correspondingly.

Download Full Size | PDF

A simple quantitative evaluation of the reconstruction results was made. Table 3 shows the density of tissues derived from δ and the signal to noise ratio (SNR) at several image regions with dimensions 20 × 20 pixels. The SNR was calculated as $\mathrm{SNR}=\overline{\rho }/\sigma$, where $\overline{\rho }$ is the calculated mean value in the considered region and σ is the standard deviation (STD). Although values measured in the sagittal plane [Fig. 6(d)] are in a decent agreement with the expected ones, in the image areas spoiled by artifacts [as the left upper part in the Fig. 6(b)] the relative deviation from the correct value can be as large as 50%. On the other hand reconstruction made in the in-plane geometry provides a good estimation of the tissue density in the entire volume.

Table 3. Quantitative results obtained from sagittal slices shown in Figs. 6(c) and 6(d). The values are measured in the regions marked by a rectangle (adipose tissue) and triangle (skin layer).

View Table | View all tables in this article

Concluding we would like to note that in PCI CT experiments at synchrotron radiation sources imaging in the out-of-plane setup might be preferred over the in-plane configuration. Typically the vertical coherence length of the synchrotron X-ray beam is at least order of magnitude larger than the horizontal one. In order to exploit larger spatial coherence length in the unidirectional X-ray DPC imaging, one would need to align the phase-contrast sensitivity axis in the vertical plane. However, in CT experiments it is also much more convenient to set the axis of rotation of the sample vertically, rather than to rotate it around the horizontal axis. In this case the DPC CT projections are acquired in the out-of-plane geometry and the index of refraction CT can be performed through the phase retrieval.

5. Discussion

We have described an approach to the problem of the phase retrieval from unidirectional X-ray DPC measurements. The method largely suppresses streak artifacts, which arise in the phase image if a straightforward 1D integration of the noisy DPC data is used. Contrary to existing methods, the considered approach does not make use of a minimization problem, but it converts the phase retrieval to the second order partial differential equation. Boundary conditions should be imposed in order to find the unique solution of a differential equation. In this work we propose to use the mixed inhomogeneous boundary conditions. The Dirichlet boundary condition allowed us to retrieve absolute values of the phase, without necessity to find the reference material in the image and to re-normalize the reconstructed values. The Neumann boundary condition was imposed in order to suppress reconstruction artifacts at the object boundaries. It is very important for the phase contrast CT, when many phase projections of the sample should be processed together. Our example illustrates that CT reconstruction of the refractive index inside an object can be made if the phase projections are first retrieved from a set of DPC projections. It is worth noting that the phase projections are ordinary tomographic projections [24] and any CT algorithm can be utilized then for the reconstruction. The method could be particularly useful if the differential phase is measured in the plane parallel to the CT rotation axis. In this case the refraction index tomography cannot be performed by means of FBP with a specialized gradient filter function (Hilbert filter).

The absolute accuracy of the suggested phase retrieval algorithm and index of refraction tomography can be tested through a direct comparison of the retrieved δ with a directly measured value. In order to do it several tissue samples should be extracted from the complex object and measurements of their refractive index or density should be performed. Another possibility is to embed a simple reference object composed of known material into a complex sample as it was done in [16]. However the estimated accuracy of the theoretically calculated δ is about 2% [26]. Also very little experimental data about the index of refraction decrement in the energy range above 30 keV are available. An exact comparison would require preliminary measurements of δ for several reference materials as it was done in [27]. The same method can be adapted to measure the refractive index of biological tissues. This work is the part of our future research.

If one applies the finite element method to solve Eq. (6) together with, for instance, the Dirichlet boundary condition, a non-degenerate, consistent, sparse linear system of equations is obtained. The problem can be scaled to big computation grids (large image sizes), and the solution should be unique regardless of the noise level. In principle, other methods can be also applied to solve Eq. (6). However the finite element discretization of the phase retrieval boundary value problem can be particularly handy in cases when additional a priori information about the sample is available and can be taken into account in the reconstruction.

It is necessary to note that in the current implementation of the method, the solution error is higher around image points where the phase projection has large gradients. Typically these are interfaces between the sample edges and the air, and between soft tissues and the bone. Improving the accuracy of the phase retrieval is the aim of our further research. One possible way to achieve this is to introduce a spatially dependent regularization parameter. It is relatively easy to automatically identify edges in the image. Then the regularization parameter can be decreased in the points close to the sharpest edges, in order to avoid excessive blurring introduced by the regularization with the second spatial derivative. In principle, several regularization terms can be added in the Eq. (6) and their contribution in each image point can be controlled by spatially dependent regularization multipliers [Eq. (7)]. The finite element method readily allows to do this without any loss of computation efficiency.

We also demonstrate that an existing software can be adapted and utilized to perform the phase retrieval using Eq. (6). This might not provide the greatest flexibility in the definition of custom boundary conditions, but it gives the opportunity to try the phase retrieval without investing time in the computer implementation of numerical techniques.

6. Conclusions

It is shown that the phase can be retrieved from unidirectional DPC images using a well-posed boundary value problem, which is essentially the 2D stationary heat conduction equation. We demonstrate that this method can suppress the typical streak artifacts associated with the phase retrieval problem. It is also shown that the suggested approach to the phase retrieval can be used for the quantitative evaluation of the results, for instance, in order to estimate magnitudes of the refractive index decrement of materials composing the sample. We also discuss a possible practical application of the phase retrieval method to refraction-based CT. As an example, we present the reconstruction of the index of refraction inside a complex biological sample. The quality of the obtained images allows to use them for non destructive examinations in preclinical studies and biomedical research as well as in the material science.

Appendix: retrieving the phase with the ANSYS^® software

There is a variety of both open source and commercial finite element software packages that can be potentially used to perform the phase retrieval using Eg. (6). Here we show how the problem can be set up and solved by means of the Workbench Mechanical component from ANSYS^® Academic Research, Release 15.0. The procedure includes several steps.

1. Calculate the right-hand side of the Eq. (6) and export the result in comma separated ASCII file as following: 1st, 2nd and, 3d columns are x,y, and z coordinates of the pixel and 4th column contain the calculated quantity. Assuming that reconstruction takes place in xy-plane, z coordinate can be set to zero everywhere.
2. Launch the workbench and select a Steady-state thermal project. Right click on Geometry and in Properties change the Analysis Type to 2D. From the component system list add the External data into the project.
3. Edit properties of the External data. Add the file with right-hand side values in the Data source panel, and in the Table of file panel assign data types according to the columns, i.e. X,Y, and Z coordinates for the first three columns and Heat generation type for the experimental data (4th column). Illustrated description can be found at http://www.edr.se/blogg/blogg/ansys_tutorial_external_data.
4. In the section Engineering data it is necessary to create a new material. It is sufficient to define the Orthotropic thermal conductivity in the Thermal properties of the created material. Unit conductivity should be selected for the DPC signal direction (y-axis in this case). The conductivity in the orthogonal direction is equal to the desired value of the regularization parameter γ.
5. In the Model section it is necessary to sketch in the xy-plane a rectangle with dimensions equal to the physical size of the DPC image. Coordinates of vertices should match the values exported to the comma separated file in the Step 1. Creating of the Surface body around the sketched rectangle concludes this step. Tutorial called 2D Steady Conduction located at https://confluence.cornell.edu/pages/viewpage.action?pageId=146918509 can provide a very good guidance for the beginner.
6. In the ANSYS Mechanical application (double click Model to launch it) a regular mesh should be mapped to the rectangular domain created in the previous step. Mesh sampling should roughly match the pixel size. Material created at Step 4 should be assigned to the surface body.
7. In the section B5 the boundary condition and the right-hand side of the equation should be defined. Temperature corresponds to the phase, so that the Dirichlet boundary condition can be set by assigning the zero temperature to all image edges. Correspondingly, the heat flow reflects information about the first phase derivative. Our right-hands side values $k\frac{\partial \alpha (x,y)}{\partial y}$ are equivalent to heat sources and sinks. To add them, right click on the Imported load folder and select Heat generation. The phase derivative imported at step 3 will be interpolated on the computation grid (see Fig. 7).
8. Solve the boundary value problem to obtain the Temperature. It is the 2D distribution of phase delays in the object image. One dimensional phase profiles can be obtained by mapping the solution to a path.

 

Fig. 7 Workspace in the ANSYS mechanical. Note that: two Edge Sizing properties control the mesh sampling; Heat flow and Temperature are defined as boundary conditions. Imported phase derivative plays a role of the Heat Generation (note contours of the Shepp-Logan phantom in the right panel). Temperature is added in the Solution list.

Download Full Size | PDF

In a similar way the phase can be retrieved with ANSYS using other boundary conditions.

Acknowledgments

Authors would like to acknowledge the financial support from Deutsche Forschungsgemeinschaft cluster of excellence - Munich Center for Advanced Photonics ( EXE158), and Open Access Publishing Fund of Karlsruhe Institute of Technology; E. Brun, H. Requardt, T. Brochard, C. Nemoz and for the technical support at the beamline, as well as V. Gasilov for helpful conversations about mathematical physics.

References and links

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

2. A. Bravin, P. Coan, and P. Suortti, “X-ray phase-contrast imaging: from pre-clinical applications towards clinics,” Phys. Med. Biol 58, R1–R35 (2013). [CrossRef]  

3. 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–337 (1996). [CrossRef]  

4. V. N. Ingal and E. A. Beliaevskaya, “X-ray plane-wave topography observation of the phase contrast from a non-crystalline object,” J. Phys. D 28, 2314–2318 (1995). [CrossRef]  

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

6. P. R. T. Munro, L. Rigon, K. Ignatyev, F. C. M. Lopez, D. Dreissi, R. D. Speller, and A. Olivo, “A quantitative, non-interferometric X-ray phase contrast imaging techniques,” Opt. Express 21, 647 (2012). [CrossRef]  

7. T. Thuering, P. Modregger, B. Pinzer, Z. Wang, and M. Stampanoni, “Non-linear regularized phase retrieval for unidirectional X-ray differential phase contrast radiography,” Opt. Express 19, 25545–25558 (2011). [CrossRef]  

8. M. Wernick, Y. Yang, I. Mondal, D. Chapman, M. Hasnah, C. Parham, E. Pisano, and Z. Zhong, “Computation of mass-density images from X-ray refraction-angle images,” Phys. Med. Biol. 51, 1769–1778 (2006). [CrossRef]   [PubMed]  

9. F. Pfeiffer, C. Kottler, O. Bunk, and C. Davis, “Hard X-ray phase tomography with low-brilliance sources,” Phys. Rev. Lett. 98, 108105 (2007). [CrossRef]   [PubMed]  

10. S. Berujon, E. Ziegler, R. Cerbino, and L. Peverini, “Two-dimensional X-Ray beam phase sensing,” Phys. Rev. Lett. 108, 158102 (2012). [CrossRef]   [PubMed]  

11. J. Sperl, D. Bequ, G. Kudielka, K. Mahdi, P. Edic, and C. Cozzini, “A Fourier-domain algorithm for total-variation regularized phase retrieval in differential X-ray phase contrast imaging,” Opt. Express 22, 450–462 (2014). [CrossRef]   [PubMed]  

12. W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, Numerical Recipes in C, 2. (Cambridge University, 1992, pp. 810).

13. J. M. Di Martino, J. L. Flores, F. Pfeiffer, K. Scherer, G. A. Ayubi, and J. A. Ferrari, “Phase retrieval from one partial derivative,” Opt. Lett. 38, 4813–4816 (2013). [CrossRef]  

14. S. Gasilov, A. Mittone, E. Brun, A. Bravin, S. Grandl, A. Mirone, and P. Coan, “Tomographic reconstruction of the refractive index with hard X-rays: an efficient method based on the gradient vector-field approach,” Opt. Express 22, 5216–5227 (2014). [CrossRef]   [PubMed]  

15. S. Gasilov, A. Mittone, E. Brun, A. Bravin, S. Grandl, and P. Coan, “On the possibility of quantitative refractive-index tomography of large biomedical samples with hard X-rays,” Biomed. Opt. Express 4, 1512–1518 (2013). [CrossRef]   [PubMed]  

16. A. Sarapata, M. Willner, M. Walter, T. Duttenhofer, K. Kaiser, P. Meyer, C. Braun, A. Fingerle, P. B. Nooel, F. Pfeiffer, and J. Herzen, “Quantitative imaging using high-energy X-ray phase-contrast CT with a 70 kVp polychromatic X-ray spectrum,” Opt. Express 23, 523–535 (2015). [CrossRef]   [PubMed]  

17. http://www.ansys.com/Products/ANSYS+15.0+Release+Highlights/ANSYS+Workbench

18. F. A. Dilmanian, Z. Zhong, B. Ren, X. Y. Wu, L. D. Chapman, I. Orion, and W. C. Thomlinson, “Computed tomography of x-ray index of refraction using the diffraction enhanced imaging method,” Phys. Med. Biol. 45, 933–946 (2000). [CrossRef]   [PubMed]  

19. M. O. Hasnah, C. Parham, E. D. Pisano, Z. Zhong, O. Oltulu, and D. Chapman, “Mass density images from the diffraction enhanced imaging technique,” Med. Phys. 32, 549–552 (2005). [CrossRef]   [PubMed]  

20. G. Birkhoff and R. E. Lynch, Numerical Solution of Elliptic Problems (SIAM Studies in Applied Mathematics1984, pp. 187).

21. A. Maksimenko, “Nonlinear extension of the X-ray diffraction enhanced imaging,” Appl. Phys. Lett. 90, 154106 (2007). [CrossRef]  

22. A. Horng, E. Brun, A. Mittone, S. Gasilov, L. Weber, T. Geith, S. Adam-Neumair, S. D. Auweter, A. Bravin, M. F. Reiser, and P. Coan, “Cartilage and soft tissue imaging using X-rays propagation-based phase-contrast computed tomography of the human knee in comparison with clinical imaging techniques and histology,” Investigative Radiology 49, 627–634 (2014). [CrossRef]  

23. A. C. Kak and M. Slaney, Principles of Computerized Tomographic Imaging (IEEE, 1988, Chap. 5).

24. G. W. Faris and R. L. Byer, “Three-dimensional beam deflection optical tomography of a supersonic jet,” Appl. Opt. 27, 5202–5212 (1988). [CrossRef]   [PubMed]  

25. Tissue Substitutes in Radiation Dosimetry and Measurement (ICRU Report 44, 1989).

26. C. T. Chantler, “Detailed tabulation of atomic form factors, photoelectric absorption and scattering cross section, and mass attenuation coefficients in the vicinity of absorption edges in the soft X-ray (Z=3036, Z=6089, E=0.1 keV10 keV),” J. Phys. Chem. Ref. Data 29, 597–1048 (2000). [CrossRef]  

27. C. M. Slack, “The refraction of X-Rays in prisms of various materials,” Phys. Rev. 27, 691–695 (1926). [CrossRef]