Non-linear regularized phase retrieval for unidirectional X-ray differential phase contrast radiography

Thomas Thüring; Peter Modregger; Bernd R. Pinzer; Zhentian Wang; Marco Stampanoni

doi:10.1364/OE.19.025545

1. Introduction

X-ray radiography is a standard method of imaging in a variety of applications including medical diagnostics, biological in-situ studies, and materials science. The interaction of X-rays with matter can be described by the complex index of refraction n = 1 – δ + iβ, where δ characterizes the phase shift and β the attenuation properties of the material. In the diagnostic energy range of X-rays (10 – 100keV), δ is typically three orders of magnitude larger than β. The consequently high interaction cross section makes the phase shift measurement a favorable imaging modality. It has been demonstrated that phase images can provide higher contrast than absorption-based images and contain complementary information about the sample [1, 2].

In the past, a variety of phase contrast techniques have been reported. They can be divided into propagation based methods [3–5], interferometric methods [6, 7] and crystal analyzer based methods [8–12]. Among these methods, phase contrast is provided either by the acquisition of the projected phase profile (ϕ), its first spatial derivative (∇ϕ) or its second derivative (∇²ϕ). Many of the reported methods require highly coherent X-rays, only available at synchrotron sources. A recently developed technique based on grating interferometry has established itself as a suitable technique for phase contrast imaging on synchrotron sources [13, 14] as well as on conventional X-ray tubes [15, 16]. The setup consists of an X-ray source, a grating interferometer of two gratings and an X-ray detector. This technique allows the simultaneous measurement of an absorption, a differential phase and a darkfield signal [17]. Due to the acquisition of the differential phase (∇ϕ), the method is also referred to as differential phase contrast (DPC) imaging. Grating interferometry has mainly been demonstrated in unidirectional mode, where the gratings have a periodic line pattern and the phase gradient can only be measured in the perpendicular direction to these lines. In bidirectional grating interferometry, “checkerboard” and “mesh” like grating patterns are used, and the DPC acquisition can be extended to both directions [18].

Compared to conventional absorption imaging, the retrieval of the phase from the raw data is in general not trivial. Paganin et al. derived a general framework for the phase retrieval in coherent imaging systems [19]. Depending on the measured representation of the phase signal (ϕ, ∇φ, ∇²φ), further post processing (integration filter) is necessary to retrieve the quantitative phase signal ϕ. Phase retrieval from the second derivative has mainly been discussed in association with propagation based methods [20, 21]. A huge variety of methods exist for differential phase techniques, where the phase is available in the first derivative. A bidirectional phase gradient acquisition has been performed on a scanning transmission X-ray microscope (STXM) using an annular quadrant detector [22], and a complex filter [23] was used for the quantitative phase retrieval. Also using an STXM, Hornberger et al. [24, 25] quantitatively retrieved the phase by deconvolving the contrast transfer function. Further, Thibault et al. [26] developed a method called scanning diffraction X-ray microscopy (SDXM), which is a combination of STXM and coherent diffractive imaging (CDI) and allows the reconstruction of the complex transfer function of an object by using an iterative algorithm.

For unidirectional grating interferometry, the phase is typically retrieved by a one-dimensional, direct integration. An integration in the image domain is equivalent to dividing by spatial frequency in the Fourier domain, which acts as low-pass filter. Therefore, along the direction of integration, the image undergoes a smoothing operation. However, in the perpendicular direction, high frequencies are unfiltered and amplify the noise. This leads to severe stripe artifacts in the image along the direction of integration. The resulting poor quantitative phase recovery is characterized by a low contrast-to-noise ratio (CNR) in the phase image, often impeding the subsequent analysis of the radiograph.

Kottler et al. also discussed this issue and proposed a bidirectional scanning approach [27] with unidirectional gratings. Two scans are performed, one with horizontally-, the other with vertically-oriented gratings. This yields a two-dimensional phase gradient of the sample and the phase can be retrieved with a complex filter. The method has shown to be well-suited for reducing stripe artifacts. However, it requires further alignment of the gratings or the sample and likewise increases scan time and dose by a factor of two. Stripe removal by acquiring the bidirectional phase gradient has further been demonstrated with a two-dimensional grating interferometer [18]. Similarly, this approach requires a two-dimensional stepping scheme and thus quadratically increases scan-time and dose compared to the unidirectional case. For diffraction enhanced imaging (DEI) and multiple image radiography (MIR), a linear least squares approach for the reduction of stripe artifacts has been reported [28].

Here, an iterative method for solving the problem of stripe artifacts from direct integration in unidirectional DPC images is proposed. The method is based on non-linear constrained optimization and is purely algorithmic. It works on the raw, noisy DPC measurements and thus does not need longer exposure times. Apart from an initial grating alignment, no further adjustments of the setup hardware or the sample are required. The stripes from the integration can be suppressed and the visibility of image details is enhanced. The performance of the method is evaluated by calculating the root mean squared error (RMSE) of the images compared to a simulated ground-truth image for phantom data and the image CNR using numerical and experimental datasets. The significant increase in CNR and the improved quantitative phase recovery makes this method a valuable tool in differential phase contrast radiography.

2. Methods

2.1. Experimental setup and phase retrieval

Figure 1 shows a typical grating interferometer setup for parallel beam geometry. The first grating of the interferometer is a phase grating, acting as a beam splitter by introducing a periodic phase shift of zero or π to the wavefront. This generates an interference pattern downstream, with maximum intensity differences at fractional Talbot distances, given by $d_{m} = m g_{1}^{2} / 8 λ$ , where m is an odd integer, g₁ the period of the phase grating and λ the wavelength [29].

Fig. 1 (a) DPC setup for a parallel beam configuration. The phase grating generates an interference pattern (beam splitter). A sample in front of the phase grating causes a lateral (x-direction) shift of the interference fringes at the position of the absorption grating. The absorption grating is used to analyze the interference fringes and determine this shift, which corresponds to the DPC signal. (b) Reference (blue) and object (light red) phase stepping curve (PSC). From these curves, the absorption, phase and a scattering signal can be calculated.

Download Full Size | PDF

A sample positioned in front of the phase grating introduces a phase shift in the wavefront, generating a phase projection profile ϕ(x, y), which is given by the line integral

ϕ (x, y) = \frac{2 π}{λ} \int δ (x, y, z) d z .

In this equation, λ is the wavelength and δ(x, y, z) is the refractive index decrement, associated to the real part of the complex index of refraction n = 1 – δ+ iβ. The partial derivative of the phase profile with respect to x is proportional to the beam refraction angle α [29]:

α (x, y) = \frac{λ}{2 π} \frac{\partial ϕ (x, y)}{\partial x} .

Beam refraction causes a shift of the interference pattern in x-direction (see Fig. 1). This fringe shift, measured in radians with respect to the fringe period g₂, corresponds to the DPC signal and is given by

φ (x, y) = 2 π \frac{d}{g_{2}} α (x, y) = \frac{λ d}{g_{2}} \frac{\partial ϕ (x, y)}{\partial x} .

Here, d is the propagation distance downstream of the phase grating and Eq. (2) was used to obtain the right part of Eq. (3). The fringe period g₂ is typically a few microns and therefore, the interference pattern cannot generally be resolved sufficiently by the X-ray detector. For the determination of φ, an absorption grating, having the same period as the fringes (g₂) and acting as an analyzer mask, is positioned at a fractional Talbot order distance. The acquisition protocol consists in a so called “phase stepping scan”, moving the absorption grating in equidistant steps in x direction and acquiring an image at each position [29]. This is done in a reference (without sample) and an object phase stepping scan (Fig. 1(b)), respectively, to account for an inhomogeneous detector illumination. From the phase stepping curves (PSC), available in each pixel, the absorption (a), the differential phase (φ), and a scattering signal (v) can be calculated by a Fourier analysis:

a = \frac{a_{0, obj}}{a_{0, ref}}

φ = φ_{1, obj} - φ_{1, ref}

v = \frac{a_{1, obj}}{a_{0, obj}} \cdot \frac{a_{0, ref}}{a_{1, ref}} = \frac{V_{obj}}{V_{ref}},

where a_i is the magnitude and φ_i is the phase of the i-th Fourier component of the PSC. Signals from reference and object scan are labeled by “obj” and “ref”, respectively. The scattering signal v corresponds to the reduced fringe visibility V_obj/V_ref.

In this paper, we focus on the recovery of the phase signal from the DPC measurement φ. According to Eq. (3), the retrieval of the phase image ϕ(x, y) requires a one-dimensional integration in x-direction, given by

ϕ (x, y) = \frac{g_{2}}{λ d} \int_{0}^{x} φ (x^{'}, y) d x^{'} .

Due to noise in the DPC image, this integration accumulates the errors and therefore the noise variance. The typical result are strong stripe artifacts, reducing image contrast and impeding an exact phase retrieval. Figure 2 shows a noisy DPC image generated from the modified Shepp-Logan phantom [30] and the phase retrieval obtained from an integration in the horizontal direction. Prior to the integration step, the mean value was removed in every line of the DPC image to enforce a cumulative sum of zero. However, this simple correction approach could not remove strong horizontal stripe artifacts which still appear across the image.

Fig. 2 Simulation of the integration of noisy DPC data using the modified Shepp-Logan phantom. a) Original data (some low-pass filtering applied), b) noisy DPC measurement, c) phase retrieval by direct integration

Download Full Size | PDF

2.2. Measurement and noise model

The method described in the next section is an algorithmic approach for the suppression of the aforementioned stripe artifacts arising from direct integration. It is based on constrained optimization and thus requires a measurement and a noise model, respectively. The DPC measurement model is a standard linear regression model, which can be written in operator notation, given by

φ = D_{x} ϕ + w .

D_x is a forward operator which models the relation between ϕ and φ in absence of noise and w is a random image modelling the noise in measured data. The noise standard deviation derived by Engel et al. [31] (Eq. 25) is given by

σ_{DPC} \propto \frac{1}{V \cdot \sqrt{N}},

where V is the visibility and N is the number of photons in the measurement.

Equation (6) is the discretization of Eq. (3), whereas the constant of proportionality ( $\frac{g_{2}}{λ d}$ ) has been omitted for the sake of simplicity. D_x models the first derivative and can be implemented as a finite differences transform operator, given by

D_{x} ϕ (i, j) = {\begin{array}{l} ϕ (i + 1, j) - ϕ (i, j) & if 1 \leq i < N_{i} \\ 0 & if x = N_{j}, \end{array}

where i and j are now discrete coordinates (pixel coordinates) and N_i is the image size in x-direction.

2.3. Regularized image integration method

A direct inversion of Eq. (6) by integrating the noisy measurement φ causes the typical horizontal stripe artifacts. The phase retrieval can be improved by solving a constrained optimization problem. While maintaining consistency with the measured DPC data, the solution is retrieved by minimizing a cost function. The constrained optimization problem for the above model can be written as

\begin{matrix} minimize {‖ D_{y} f ‖}_{ℓ_{p}} \\ subject to : {‖ W (D_{x} f - φ) ‖}_{ℓ_{2}} < ɛ . \end{matrix}

Essentially, this optimization problem seeks for the minimal ℓ_p norm of D_y f for any f being consistent with the measurement data φ in a least squares sense. W is a weighting operator, multiplying each pixel with its inverse noise standard deviation (1/σ) to account for the noise model. This approach is also known under the name penalized weighted least squares (PWLS) [32]. D_x and D_y are finite differences transforms in the x- and y-direction, respectively, f is the solution of the optimization problem and ɛ is a boundary for the noise. The ℓ_p norm of an image is defined as

{‖ a ‖}_{ℓ_{p}} = {(\sum_{i = 1}^{N_{i}} \sum_{j = 1}^{N_{j}} {| a (i, j) |}^{p})}^{1 / p} .

For instance, the most often used ℓ₂ norm corresponds to the Euclidean norm. In this case, problem (9) is linear and an explicit inversion for f exists. Here, we focus on the use of the ℓ₁ norm, which is the sum of the magnitude of all pixel values. The ℓ₁ norm minimization of the finite differences transform is well-known in denoising applications and often referred to as total variation denoising (two-dimensional finite differences) [33]. The ℓ₁ norm minimization typically shows strong edge preserving characteristics and causes less blurring than the ℓ₂ norm. The better performance comes at the expense of a non-linear optimization problem, having no explicit inversion formula.

For the solution of (9) with p = 1, the problem is first reformulated to an unconstrained form:

minimize F (f) = {‖ W (D_{x} f - φ) ‖}_{ℓ_{2}}^{2} + λ {‖ D_{y} f ‖}_{ℓ_{p}} .

The data consistency constraint and the optimization function are now expressed in a single objective function F(f). A Lagrange multiplier λ (regularization parameter) is used for weighting the two terms. The parameter ɛ from Eq. (9) is now implicitly chosen by adjusting λ. This parameter is strongly data-dependent and in general difficult to determine. Here, it will be determined by evaluating well-known image quality metrics.

Equation (11) can further be generalized to hold multiple regularization terms. The problem formulation is then given by

\begin{array}{l} minimize F (f) = {‖ W (D_{x} f - φ) ‖}_{ℓ_{2}}^{2} + λ_{1} {‖ T_{1} f ‖}_{ℓ_{p_{1}}}^{p_{1}} + \\ + λ_{2} {‖ T_{2} f ‖}_{ℓ_{p_{2}}}^{p_{2}} + \dots \end{array}

In principle, this allows for any number of regularization terms to be added, which is particularly useful if more a-priori knowledge about the sample is available. Minimum norm constraints with any norm number p and for any linear transform Tf can be added. An example for regularized integration using two regularization terms will be given in the next section.

For the numerical solution of the unconstrained optimization problem, a non-linear Conjugate Gradient (NLCG) algorithm [34] has been used. NLCG is characterized by a fast convergence for the inversion of large-scale linear systems. Moreover, this algorithm can handle both, the linear (p = 2) and the non-linear (p = 1) case and thus provides a convenient way for a comparison. The implementation was done in Matlab (MathWorks, Version R2010a) and run on a 64-bit Linux machine (Intel Xeon 2.83GHz processor).

3. Results

3.1. Numerical data

Numerical phantom data have been used for the theoretical assessment of the method. The data were generated according to Eq. (6), where ϕ was the known (ground-truth) image, i.e. the Shepp-Logan phantom, representing the line integral of δ (real part of the refractive index) through an object. Noise has been generated by using the previously discussed model for DPC measurements (Eq. (7)). The calculation of a noise standard deviation map required the simulation of a transmission and a visibility image. For the transmission image, we assumed a constant ratio δ/β = 10 (β is the imaginary part of the refractive index), yielding the transmission image directly from the known phase image. For the visibility image, a high scattering signal at strong edges of the object was assumed and small angle scattering was neglected. Therefore, an edge enhanced image was generated by applying a discrete two-dimensional Laplace filter to the phantom image, from which the visibility map was determined. Figure 3 shows the simulated transmission and visibility image and the combined noise standard deviation map. Using the noise standard deviation map, noisy DPC data, as shown in Fig. 2(b), was generated as input for the regularized image integration.

Fig. 3 Simulation of (a) the transmission image and (b) the visibility map of the Shepp-Logan phantom. (c) noise standard deviation map (normalized), obtained by combining (a) and (b) according to Eq. (7).

Download Full Size | PDF

In the following, regularized and direct integration are evaluated in terms of well known image metrics, including root mean squared error (RMSE) compared to the ground-truth image, contrast-to-noise ratio (CNR) and the noise pattern standard deviation. The ground-truth image is defined as the radiographic phase image of the sample in absence of errors (e.g., noise, material imperfections or measurement errors). RMSE and CNR are computed as a function of the regularization parameter λ. First, this allows a direct comparison of regularized and conventional direct integration and second, it provides support for the choice of an optimal λ, although this is of course dependent on the chosen metric. In order to compare the the non-linear ℓ₁ and the linear ℓ₂ norm minimization, the simulations are performed for p ∈ {1, 2}. For all simulations, the same convergence criterion in the NLCG algorithm is applied.

Figure 4(a) shows the phase image obtained by direct integration and Fig. 4(b) displays the optimal result of the regularized integration in terms of the RMSE. The RSME of two images f₁ and f₂ is defined as

RMSE (f_{1}, f_{2}) = \frac{{‖ f_{1} - f_{2} ‖}_{ℓ_{2}}}{\sqrt{N_{i} N_{j}}},

where N_i and N_j are the horizontal and vertical image sizes, respectively. Figure 4(e) displays the RMSE as a function of the regularization parameter λ and with p ∈ {1, 2}, for regularized integration (RI), and the RMSE for direct integration (DI). For all values of λ within the investigated range, a smaller RMSE could be measured for regularized integration. The minimum was achieved with the ℓ₁ norm minimization at λ_RMSE = 5 · 10⁻³, yielding the optimal solution (displayed in Fig. 4(b)).

Fig. 4 On the left side, simulation results on the Shepp-Logan phantom are shown. (a) Direct integration of the numerical phantom data of Fig. 2(b), causing the typical stripe artifacts. (b) Regularized image integration with p = 1 and λ = 5·10⁻³, containing no stripe artifacts. (c) noise pattern for the direct integration and (d) for the regularized integration, respectively, calculated by subtracting the output image from the ground-truth image. The noise standard deviations are displayed on the images. The right side of the figure shows the evaluation of (e) the RMSE and (f) the CNR as a function of the regularization parameter and for p ∈ {1, 2}. The red boxes in (a) show the regions used for the calculation of the CNR. The optimal regularization parameter for both metrics is indicated with a vertical line (λ_RMSE and λ_CNR). Figure (g) and (h) show vertical and horizontal line profiles, respectively, through the images in (a), (b) and the ground-truth image. The locations of the extracted line profiles are indicated with the dashed line in (b). In these profiles, the reduction of high vertical variations and the improved quantitative phase recovery in the x-direction compared to direct integration is clearly visible.

Download Full Size | PDF

Regularized image integration changes the noise pattern of the images. Apparently, direct integration introduces a horizontally oriented stripe pattern, as shown in Fig. 4(c). After regularization, the pattern is stripe-free and more homogeneous (Fig. 4(d)), although a few new artifacts, mainly at strong horizontal edges of the sample (top and bottom of phantom skull), have been introduced. Artifacts at horizontal edges appear due to the low DPC signal at these locations, which is a result of the unidirectional sensitivity (horizontal direction) of the grating interferometer. The standard deviation of the noise patterns, indicated by σ in Fig. 4(c) and (d), is smaller for the regularized integration, confirming the improved image quality.

The RMSE metric requires the knowledge of a ground-truth image for the computation, which is in practice unknown. Therefore, the method has further been evaluated using the image CNR, which represents a well-known metric for the quantification of the image contrast and is independent of ground-truth data. CNR is defined as

CNR = 2 \frac{| S_{obj} - S_{bg} |}{σ_{obj} + σ_{bg}},

where S and σ are the mean value and the standard deviation, respectively, within a region-of-interest (ROI) of the object (obj) and the background (bg). Figure 4(f) shows the CNR as a function of λ, calculated by using the object and background regions marked with the red boxes in Fig. 4(a). The ℓ₁ norm minimization again performed better than the ℓ₂ norm and the optimal regularization parameter is λ_CNR = 4 · 10⁻², which is almost an order of magnitude larger than the optimal value in terms of RMSE. This is mainly because high regularization parameters significantly suppress noise, resulting in a higher CNR. On the other hand, a high λ can smooth the image and reduce edged sharpness. For this reason, optimizing the CNR is not always the best approach and may lead to an overestimation of λ.

The previous results showed that regularized image integration outperformed direct integration of differential phase contrast data in terms of RMSE, CNR and noise variance. The vertical and horizontal line profiles through the image of Fig. 4(b), as shown in Fig. 4(g) and (h), further illustrate how well the vertical variations could be suppressed and that the regularized image is in significantly higher accordance with the ground-truth image. Moreover, the non-linear ℓ₁ norm minimization performed clearly better than the ℓ₂ norm for both metrics. For this reason, in the examples of the following sections, only the ℓ₁ norm minimization will be used for the regularized phase retrieval.

3.2. Experimental data

The method has further been tested with experimental DPC data. In order to investigate the quantitative robustness of the stripe removal method also for experimental data, a phantom made of materials with known indices of refraction was fabricated and imaged. The phantom has been acquired with the DPC setup of the TOMCAT beamline at the Swiss Light Source [35]. The X-ray source is a 2.9T superbend bending magnet. A double crystal multilayer monochromator delivers monochromatic radiation. The photon energy was set to 25keV and the grating interferometer was designed for an inter-grating distance of the 3^rd fractional Talbot order. The source-to-phase grating distance was 25m, the inter-grating distance was d = 120mm and the grating periods were g₁ = 3.98μm and g₂ = 2.0μm. The materials of the absorption and phase grating were gold and silicon, respectively [36]. The total exposure time was 1.6 seconds, using 8 phase steps.

The phantom consists of three concentric rods, which are surrounded by a circular wall. The wall is made of polymethylmetharcylate (PMMA) and the materials for the rods are polypropylen (PP), polystyrene (PS) and polyoxymethylene (POM), respectively. The space between the rods and the wall was filled with de-ionized water. The complex refractive index n = 1 – δ + iβ at 25keV of each material was known (available at http://sergey.gmca.aps.anl.gov), allowing for a calculation of the ground-truth image (calculated analytically by assuming a known sample profile).

Figure 5(a)–(d) shows the DPC, integrated, regularized and ground-truth image, respectively, of a projection of the phantom and Fig. 5(e) shows the normalized horizontal profiles after summing up the image rows in vertical direction. In this example, an additional image constraint was applied to demonstrate the possibilities of the generalized problem formulation introduced in the previous section. The additional constraint is based on a-priori knowledge about the image boundaries, assuming a background signal of zero. This assumption holds because the DPC image is flat-field corrected, essentially meaning that in a background pixel, the phase signals of the object and the reference scan should be the same. The optimization problem, consisting now of two regularization terms, is defined as:

minimize F (f) = {‖ W (D_{x} f - φ) ‖}_{ℓ_{2}}^{2} + λ_{1} {‖ D_{y} f ‖}_{ℓ_{1}} + λ_{2} {‖ M f ‖}_{ℓ_{2}}^{2} .

In the second regularization term, the matrix M operates as a binary mask for the pixels in the image f at the left and right edges. For the norm parameter, p = 2 was chosen, because this constrains the image to have low energy at the edge pixels. While the first regularization parameter λ₁ was varied, the second regularization parameter was kept constant at λ₂ = 10⁻². Since the second regularization term only constrains parts of the image, the phase retrieval is not very sensitive to λ₂, making an exact determination of an optimum for this parameter unnecessary. There is a high certainty that a background signal of zero is a true assumption, therefore the value should be kept high. On the other hand, in order to avoid a too strong influence of this term to the result, λ₂ should at the same time be kept small. The choice of λ₂ = 10⁻² is purely empirical and does not represent a unique optimum.

Fig. 5 Projection images of a phantom, acquired at 25keV at the DPC setup of the TOM-CAT beamline at the Swiss Light Source. (a) DPC image, (b) direct integration, (c) regularized integration using p₁ = 1, λ₁ = 6 · 10⁻⁴, p₂ = 2 and λ₂ = 10⁻² (d) ground-truth image. The colormaps of the images (b)–(d) all have the same minimum and maximum value. (e) shows the evaluation of RMSE and CNR for DI and RI and (f) a normalized horizontal profile after summing up the images in the vertical direction. CNR has been calculated using the ROIs marked with the red boxes in (b).

Download Full Size | PDF

Figure 5 shows the evaluation of the RMSE and the CNR for the regularized and the directly integrated phase recovery of the phantom. In this example, the optimal values for the regularization parameter match exactly (λ_1,RMSE = λ_1,CNR = 5 · 10⁻⁴). Again, the suppression of stripe artifacts in the regularized integration is clearly visible, if Fig. 5(b) and (c) are compared. According to Fig. 5(e), the quantitative phase recovery could be improved over the whole investigated range of the regularization parameter, although a complete agreement to the calculated ground-truth was not achieved (Fig. 5(f)). However, it is important to notice that the result of Fig. 5(f) must be taken with care due to the following important facts: In this example, RI was compared with DI after a vertical summation of the images, not by looking at line profiles (as it was done in the simulations). This summation also smoothes out the stripe artifacts in DI and therefore, the DI line in Fig. 5(f) would be expected to highly match the ground-truth line, which is not the case. This inconsistency might result from several random effects, as for example inhomogeneities in the material, inexact numbers for the material densities and the refractive indeces or inconsistencies of the DPC measurements resulting from grating drift or grating imperfections from the fabrication process. Due to these facts, it is impossible to state that RI was indeed quantitatively more accurate than DI, simply because the ground-truth is actually unknown. On the other hand, due to the removal of the stripe artifacts, it is very likely that a line profile from RI is quantitatively more accurate than a line profile (no vertical summation) of DI.

Figure 6(a) shows a DPC image of a mouse (scanned within a sample holder). The scan was taken on a grating interferometer setup using an X-ray tube source, operating at 40kV, 25mA. The pixel size of the detector was 50 × 50μm². The anode material of the tube is tungsten and the spotsize is approximately 0.5mm. Due to the low spatial coherence of the beam, a source grating was used [15]. 16 phase steps were acquired and the total exposure time was 120 seconds. The full body of the mouse was scanned by stitching multiple scans together. The total size of this image is 2360 × 1200 pixels. The standard problem formulation of Eq. (11) with p = 1 was used to retrieve the regularized phase image. According to a CNR analysis, the optimal regularization parameter was λ_CNR = 4 · 10⁻³. The integrated and regularized images are shown in Fig. 6(b) and 6(c), respectively. Significant improvements compared to direct integration can be observed in the images. Features that are hidden in the directly-integrated image become visible in the regularized integration (see ribs in zoomed area). Due to the large image size, the algorithm took 450 iterations with a total run time of approximately 11 minutes.

Fig. 6 (a) DPC image of a mouse taken at an X-ray tube setup (40kV, 25mA), (b) integrated image, (c) regularized integration for λ = 4 · 10⁻³. In the zoomed area below, the ribs of the mouse can clearly be identified in the regularized integration, while these details vanished in the direct integration. The colormaps of the images in (b) and (c) with the same zoom level have the same minimum and maximum value.

Download Full Size | PDF

4. Discussion

Non-linear, ℓ₁ regularized image integration can significantly improve the quality of phase images obtained from noisy, unidirectional DPC measurements without increasing the radiation dose. Longer exposure times, additional phase steps, or bidirectional measurements can be avoided. The method outperformed direct integration in the given examples in terms of RMSE compared to ground-truth data, noise pattern and CNR. The improved RMSE promises more accurate results for a quantitative data analysis, while a higher CNR will facilitate data post processing as for instance image segmentation.

The application of a regularization technique to a problem, which is in principle well-posed, might seem a bit unusual at a first glance. However, since the integration of differential data is highly sensitive to noise and generates strong image artifacts, regularization provides a convenient way to specify data fidelity. This is done using the regularization parameter, λ, being the key parameter of the method. It determines the weighting between the data consistency and the optimization function and can be interpreted as a tuning parameter for the resulting image quality. If λ is chosen too high, overregularization occurs, generating new image artifacts. A typical overregularization effect when using the ℓ₁ norm minimization of the finite differences transform is the generation of piecewise constant areas in the image, looking similar to median filtering. If λ is below its optimum, the regularization term is too weak and the stripe artifacts will not vanish. Naturally, λ is data dependent and the optimum is not always trivial to determine. The evaluation of image metrics (CNR, RMSE) for the optimization of λ is one possible approach amongst others. Another well known method is the L-curve analysis, where the data constraint term and the minimization terms are evaluated after each iteration [37]. Different methods usually yield different optima, and therefore, an objective optimum does usually not exist. CNR and RMSE turned out to be convenient metrics for the image evaluation and the results showed, that the regularized images have been superior to the directly integrated images over a wide range of λ. Fortunately, this relaxes the requirement to exactly find an optimum value for λ. Eventually, finding the optimum solution will remain a task dependent issue.

Another important parameter is the norm parameter p. It has been demonstrated, that for the minimization of the finite differences transform, the ℓ₁ norm performed better than the ℓ₂ norm in terms of the evaluated metrics. The ℓ₁ norm has also been subject to many investigations in the field of tomographic image reconstruction, where an image is assumed to have a sparse representation in a linear transform domain. The minimization of the ℓ₁ norm in such a domain enforces sparsity and may improve image reconstruction. This concept is also consistent to our problem, since the finite difference transform in the y-direction is expected to yield a sparse image representation. The benefit of minimizing the ℓ₁ norm instead of the ℓ₂ norm comes at the expense of a non-linear problem formulation. For the finite difference transform, the ℓ₁ norm has been shown to perform clearly better, however for other transform domains, the ℓ₂ norm may be a better choice.

The possibility of introducing theoretically any number of regularization terms improves the flexibility of the method to handle many kinds of samples. In particular sample specific constraints, which should be added if a-priori knowledge about the sample is available, are expected to further improve the phase retrieval. For instance, if the sample is known to consist mainly of piecewise constant parts, a total variation minimization would be appropriate. Another example, which is less sample specific but also often applicable, would be an “isolated specimen” constraint, meaning that there is a zero background signal all around the sample. Essentially, this would be an generalization of the proposed zero-background constraint, which assumed the signal to be zero only at the left and right edges. In any case, the potential improvement always comes at the expense of an additional free regularization parameter and an increased computational complexity. Additional regularization terms are certainly an important subject for future investigations.

The iterative NLCG algorithm has proven to be a suitable method for solving the non-linear optimization problem, even though the processing time can rapidly increase for very large image sizes. Code optimization, which has not yet been performed, should lead to a significant decrease in computational time.

It is important to point out that non-linear, regularized image integration provides a general signal recovery method for unidirectional differential data. In this work, it has only been applied to data obtained from grating interferometry, although other techniques as for example crystal analyzer based methods would likewise benefit from our approach.

We strongly believe that DPC radiography is a valuable imaging method for various applications in radiology. In order to effectively exploit the additional information obtained from this technique, the proposed method provides an important tool for obtaining the quantitative phase projection. This is crucial for improving the analysis of structures in the examined materials (e.g., in medical images).

Acknowledgments

We would like to acknowledge S. Rutishauser and C. David (PSI Villigen) for providing the grating interferometer setup and M. Cacquevel (Brain Mind Institute, École Polytechnique Fédérale de Lausanne, Switzerland) for providing the mouse sample. All procedures involving the mouse sample were approved by the Committee on Animal Experimentation for the canton of Vaud, Switzerland, in accordance with Swiss Federal Laws on Animal Welfare and the European Community Council directive (86/609/EEC) for the care and use of laboratory animals. This study was supported by Centre d’Imagerie BioMédicale (CIBM) of the UNIL, UNIGE, HUG, CHUV, EPFL and the Leenaards and Jeantet Foundations.

References and links

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

2. A. Momose and J. Fukuda, “Phase-contrast radiographs of nonstained rat cerebellar specimen,” Med. Phys. 22, 375–379 (1995). [CrossRef] [PubMed]

3. A. Snigirev, I. Snigireva, V. Kohn, S. Kuznetsov, and I. Schelokov, “On the possibilities of X-ray phase contrast microimaging by coherent high-energy synchrotron radiation,” Rev. Sci. Instrum. 66, 5486–5492 (1995). [CrossRef]

4. S. Wilkins, T. Gureyev, D. Gao, A. Pogany, and A. Stevenson, “Phase-contrast imaging using polychromatic hard X-rays,” Nature 384, 335–338 (1996). [CrossRef]

5. P. Cloetens, R. Barrett, J. Baruchel, J. Guigay, and M. Schlenker, “Phase objects in synchrotron radiation hard X-ray imaging,” J. Phys. D. 29, 133–146 (1996). [CrossRef]

6. U. Bonse and M. Hart, “An X-ray interferometer,” Appl. Phys. Lett. 6, 155–156 (1965). [CrossRef]

7. A. Momose, T. Takeda, Y. Itai, and K. Hirano, “Phase-contrast X-ray computed tomography for observing biological soft tissues,” Nat. Med. 2, 473–475 (1996). [CrossRef] [PubMed]

8. T. Davis, D. Gao, T. Gureyev, A. Stevenson, and S. Wilkins, “Phase-contrast imaging of weakly absorbing materials using hard X-rays,” Nature 373, 595–598 (1995). [CrossRef]

9. D. Chapman, W. Thomlinson, R. Johnston, D. Washburn, E. Pisano, N. Gmür, Z. Zhong, R. Menk, F. Arfelli, and D. Sayers, “Diffraction enhanced X-ray imaging,” Phys. Med. Biol. 42, 2015–2025 (1997). [CrossRef] [PubMed]

10. W. Boettinger, H. Burdette, and M. Kuriyama, “X-ray magnifier,” Rev. Sci. Instrum. 50, 26–30 (1979). [CrossRef] [PubMed]

11. M. Stampanoni, G. Borchert, and R. Abela, “Progress in microtomography with the Bragg Magnifier at SLS,” Rad. Phys. Chem. 75, 1956–1961 (2006). [CrossRef]

12. P. Modregger, D. Lübbert, P. Schäfer, and R. Köhler, “Two dimensional diffraction enhanced imaging algorithm,” Appl. Phys. Lett. 90, 193501 (2007).

13. C. David, B. Nöhammer, H. Solak, and E. Ziegler, “Differential X-ray phase contrast imaging using a shearing interferometer,” Appl. Phys. Lett. 81, 3287–3289 (2002). [CrossRef]

14. A. Momose, S. Kawamoto, I. Koyama, Y. Hamaishi, K. Takai, and Y. Suzuki, “Demonstration of X-Ray Talbot Interferometry,” Jpn. J. Appl. Phys. 42, L866–L868 (2003). [CrossRef]

15. F. Pfeiffer, T. Weitkamp, O. Bunk, and C. David, “Phase retrieval and differential phase-contrast imaging with low-brilliance X-ray sources,” Nature Phys. 2, 258–261 (2006). [CrossRef]

16. M. Engelhardt, J. Baumann, M. Schuster, C. Kottler, F. Pfeiffer, O. Bunk, and C. David, “High-resolution differential phase contrast imaging using a magnifying projection geometry with a microfocus X-ray source,” Appl. Phys. Lett. 90, 224101 (2007). [CrossRef]

17. F. Pfeiffer, M. Bech, O. Bunk, P. Kraft, E. Eikenberry, C. Brönnimann, C. Grünzweig, and C. David, “Hard-X-ray dark-field imaging using a grating interferometer” Nature Mater. 7, 134–137 (2008). [CrossRef]

18. I. Zanette, T. Weitkamp, T. Donath, S. Rutishauser, and C. David, “Two-Dimensional X-Ray Grating Interferometer,” Phys. Rev. Lett. 105 (2010). [CrossRef]

19. D. Paganin, “Phase retrieval using coherent imaging systems with linear transfer functions,” Opt. Commun. 234, 87–105 (2004). [CrossRef]

20. D. Paganin, S. Mayo, T. Gureyev, P. Miller, and S. Wilkins, “Simultaneous phase and amplitude extraction from a single defocused image of a homogeneous object,” J. Microsc. 206, 33–40 (2002). [CrossRef] [PubMed]

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

22. M. de Jonge, B. Hornberger, C. Holzner, D. Legnini, D. Paterson, I. McNulty, C. Jacobsen, and S. Vogt, “Quantitative Phase Imaging with a Scanning Transmission X-Ray Microscope,” Phys. Rev. Lett. 100, 163902 (2008). [CrossRef] [PubMed]

23. M. Arnison, K. Larkin, C. Sheppard, N. Smith, and C. Cogswell, “Linear phase imaging using differential interference contrast microscopy” J. microsc. 214, 7–12 (2004). [CrossRef] [PubMed]

24. B. Hornberger, M. Feser, and C. Jacobsen, “Quantitative amplitude and phase contrast imaging in a scanning transmission X-ray microscope” Ultramicroscopy 107, 644–655 (2007). [CrossRef] [PubMed]

25. B. Hornberger, M. de Jonge, M. Feser, P. Holl, C. Holzner, C. Jacobsen, D. Legnini, D. Paterson, P. Rehak, L. Strüder, and S. Vogt, “Differential phase contrast with a segmented detector in a scanning X-ray microprobe” J. Synchrotron Radiat. 15, 355–362 (2008). [CrossRef] [PubMed]

26. P. Thibault, M. Dierolf, A. Menzel, O. Bunk, C. David, and F. Pfeiffer, “High-resolution scanning X-ray diffraction microscopy” Science 321, 379–382 (2008). [CrossRef] [PubMed]

27. C. Kottler, C. David, F. Pfeiffer, and O. Bunk, “A two-directional approach for grating based differential phase contrast imaging using hard X-rays,” Opt. Express 15, 1175–1181 (2007). [CrossRef] [PubMed]

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

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

30. L. Shepp and B. Logan, “The Fourier reconstruction of a head section,” IEEE Trans. Nucl. Sci. 21, 21–43 (1974).

31. K. Engel, D. Geller, T. Köhler, G. Martens, S. Schusser, G. Vogtmeier, and E. Rössl, “Contrast-to-noise in X-ray differential phase contrast imaging,” Nucl. Instrum. Methods A 648, 202–207 (2010). [CrossRef]

32. K. Sauer and C. Bouman, “A local update strategy for iterative reconstruction from projections,” IEEE Trans. Signal Process. 41, 534–548 (1993). [CrossRef]

33. L. Rudin, S. Osher, and E. Fatemi, “Nonlinear total variation based noise removal algorithms,” Phys. D Nonlinear Phenom. 60, 259–268 (1992). [CrossRef]

34. J. Nocedal and S. Wright, Numerical optimization (Springer verlag, 1999). [CrossRef]

35. M. Stampanoni, A. Groso, A. Isenegger, G. Mikuljan, Q. Chen, A. Bertrand, S. Henein, R. Betemps, U. Frommherz, P. Böhler, D. Meister, M. Lange, and R. Abela, “Trends in synchrotron-based tomographic imaging: The SLS experience,” Proc. SPIE 6318, 63180M (2006). [CrossRef]

36. C. David, J. Bruder, T. Rohbeck, C. Grünzweig, C. Kottler, A. Diaz, O. Bunk, and F. Pfeiffer, “Fabrication of diffraction gratings for hard X-ray phase contrast imaging,” Microelectron. Eng. 84, 1172–1177 (2007). [CrossRef]

37. P. Hansen, “Analysis of Discrete Ill-Posed Problems by Means of the L-Curve,” SIAM Rev. 34, 561–580 (1992). [CrossRef]

Non-linear regularized phase retrieval for unidirectional X-ray differential phase contrast radiography

Abstract

1. Introduction

2. Methods

2.1. Experimental setup and phase retrieval

2.2. Measurement and noise model

2.3. Regularized image integration method

3. Results

3.1. Numerical data

3.2. Experimental data

4. Discussion

Acknowledgments

References and links

Cited By

Figures (6)

Equations (17)

Optics Express