Phase retrieval in in-line x-ray phase contrast imaging based on total variation minimization

Alexander Kostenko; K. Joost Batenburg; Heikki Suhonen; S. Erik Offerman; Lucas J. van Vliet

doi:10.1364/OE.21.000710

1. Introduction

Recently, the field of X-ray phase-contrast imaging (PCI) has been growing rapidly. X-ray PCI found applications in materials science, ranging from investigating the microstructure of carbon-based materials [1, 2] to in-situ measurements of dynamic processes taking place in metal alloys and semiconductors [3–5]. X-ray PCI is also entering the field of pre-clinical biomedical research, namely, small animal imaging and various ex-vivo/in-vitro studies [6–11]. The increasing availability of the X-ray PCI techniques over the last years was stimulated by advances in instrumentation and phase retrieval algorithms.

The scope of our current investigation lays within phase retrieval for propagation-based X-ray PCI. For more than a decade various algorithms were developed to permit accurate reconstruction of the specimen’s phase and attenuation images from phase-contrast data acquired using the propagation-based approach [12]. A major effort was aimed at the development of linear approximations to the image formation of PCI that would permit a stable solution of the resulting inverse problem [12–16]. Using these approximations, the phase and attenuation images of the specimen can be calculated from a series of phase-contrast images acquired at different propagation distances. To allow phase retrieval from a single phase-contrast image, methods based on prior information about the specimen’s composition were developed [17–19]. They are referred to as the so called phase-attenuation duality models.

In all linear phase retrieval models the accuracy of the reconstruction is a function of spatial frequency. Depending on the acquisition conditions, the signal-to-noise ratio (SNR), and the fitness of the linear approximation, the phase image of the specimen will be irrecoverable within a particular set of spatial frequencies [20, 21]. In the case of multi-distance phase retrieval [16] this can lead to large errors at low spatial frequencies while in single-distance approaches [17–19] artifacts are produced at middle and high spatial frequencies. In order to avoid large errors in the reconstructed images, most of the phase retrieval approaches rely on a so called L2-norm based regularization also known as Tikhonov regularization. When L2 regularization is used, the solution that has the minimum L2-norm (i.e. Euclidean norm) is promoted. This leads to a suppression of spatial frequencies that are ill-determined by the phase retrieval model or heavily corrupted by noise. Such a solution may not be optimal, especially when it results in a strong suppression of a large band of low frequencies in multi-distance retrieval methods.

Another regularization approach that is currently used in an increasing number of image reconstruction applications is called Total Variation (TV) minimization. It was initially developed for image denoising [22] and recently been introduced in such fields as deblurring, super-resolution, inpainting and tomography [23–26]. The idea underlying TV minimization is to promote a solution that has the sparsest gradient. It was theoretically proven [27] that under certain conditions TV minimization allows exact reconstruction of signals with a sparse gradient from highly incomplete sets of observations. In cases where the gradient of the reconstructed image is not exactly sparse, TV minimization is nevertheless preferred to L2 regularization in many applications [28].

In this paper we introduce a TV minimization approach for solving the inverse problem of phase retrieval in propagation-based X-ray PCI based on various linear models. Any implementation of TV minimization can be chosen from a wide range of algorithms [25]. Here we will present only results acquired with an algorithm based on the so called Fast Iterative Shrinkage-Thresholding Algorithm (FISTA) [29]. The original implementation of the algorithm was modified to include frequency weighting into the minimization scheme. Frequency weighting permits to account for the frequency-dependent nature of the signal-to-noise ratio and is shown to have a significant influence on the accuracy of the phase retrieval.

2. Materials and methods

In the following subsections we will describe how the phase retrieval problem of phase retrieval can be solved using iterative least-squares minimization and TV-minimization. In order to do so, we will introduce a matrix formalism which is uncommon in the field of phase retrieval. To keep our description compact we will refer the reader to the paper of M. Langer [16] for details regarding the theoretical background of X-ray phase retrieval algorithms.

2.1. Matrix formalism for phase propagation model

It is essential for the purpose of this paper that the convolution integral can be expressed using matrix formalism. Let us, for instance, define the propagated X-ray wavefront as a matrix product. Let 𝒜 denote the set of square integrable functions ℝ² → ℂ with bounded support. Within the paraxial approximation, the X-ray field H_D ∈ 𝒜 propagated to a distance D from the object can be expressed as the convolution (★) of the unpropagated field T ∈ 𝒜 with the Fresnel propagator P_D ∈ 𝒜 :

H_{D} = P_{D} ★ T .

If we denote the Fourier domain representations of T, P_D, and H_D by T̃, P̃_D, and H̃_D, respectively, we find that

{\tilde{H}}_{D} (f) = {\tilde{P}}_{D} (f) \cdot \tilde{T} (f),

where f denotes the spatial frequency. We now discretize the Fourier domain, representing a spectrum by its values in a discrete set {f₁,..., f_k} of basis elements

\tilde{t} = (\begin{matrix} \tilde{T} (f_{1}) \\ ⋮ \\ \tilde{T} (f_{k}) \end{matrix}) and {\tilde{P}}_{D} = (\begin{matrix} {\tilde{P}}_{D} (f_{1}) & \emptyset \\ {\tilde{P}}_{D} (f_{2}) \\ ⋱ \\ \emptyset & {\tilde{P}}_{D} (f_{k}) \end{matrix}) .

where P̃_D is a diagonal matrix containing a discrete representation of the propagator P̃_D(f_k) and t̃ is a vector that contains a discrete representation of T̃(f_k). Consequently, a discretized representation h̃_D of the propagated field H̃_D can be defined in frequency domain:

{\tilde{h}}_{D} = {\tilde{P}}_{D} \cdot \tilde{t},

where the symbol (·) denotes the matrix-vector product. If we introduce a vector t corresponding to the discrete representation of the unpropagated field T, it is possible to construct a Toeplitz matrix P_D such that the discretized propagated field h_D can be expressed in the spatial domain:

h_{D} = P_{D} \cdot t

As matrix P_D is dense, it is hard to compute Eq. (5) directly or solve the inverse problem for a large vector t in the spatial domain. However, the fact that matrix P_D has a complimentary frequency space representation P̃_D, both play an important role in computing the explicit deconvolution and iterative least-squares inversion. Since P̃_D is a diagonal matrix, its eigenvalues correspond to the discrete representation of the propagator P_D(f_k). That property allows us to construct the well-known explicit deconvolution Eq. (8) in Section 2.2 and to analyse the properties of iterative least-squares inversion in Section 2.3.

2.2. Linear phase retrieval algorithms: L2-norm regularization

The standard approach to phase retrieval in in-line x-ray PCI relies on the fact that the observed phase-contrast image can be approximated by a linear transformation acting on some unknown image (or a combination of a phase and an attenuation image). This problem can be viewed as a system of linear equations:

b = A \cdot x

where A is a matrix that represents one of the linear phase models, b a vector containing the observed phase-contrast images, and x a vector containing the unknown images. Since the observed images are corrupted by noise and inversion of A is usually ill-posed, Eq. (6) is often replaced by an L2-norm regularized least-squares problem (a.k.a. a Tikhonov regularization problem):

\underset{x}{argmin} : {‖ A \cdot x - b ‖}_{2}^{2} + ε_{L 2} {‖ x ‖}_{2}^{2},

where

{‖ \dots ‖}_{2}^{2}

denotes the L2-norm and ε the regularization weight. The minimization of the first term guarantees the best fit of the linear phase-contrast model to the observed images, while the second term promotes solutions with the smallest L2-norm thereby suppressing noise and outliers. Given that matrix A has a diagonal representation in frequency domain Ã, Eq. (7) has an analytical solution according to which each j^th frequency component x̃_j can be calculated using the following expression:

{\tilde{x}}_{j} = \frac{{\tilde{A}}_{j, j}^{*} {\tilde{b}}_{j}}{{| {\tilde{A}}_{j, j} |}^{2} + ε_{L 2}},

where

{\tilde{A}}_{j, j}^{*}

and |Ã_j,j|² are respectively the conjugate and the squared magnitude of the diagonal matrix Ã. When in underdetermined or ill-conditioned system |Ã_j,j| → 0, the corresponding frequency component of x̃ will be suppressed by L2 regularization:

{\tilde{x}}_{j} \approx \frac{{\tilde{A}}_{j, j}^{*} {\tilde{b}}_{j}}{ε_{L 2}} .

2.3. Linear phase retrieval algorithms: TV minimization

Similar to the L2-regularized approach Eq. (6), Eq. (7) can be replaced by a TV-regularized version:

\underset{x}{argmin} : {‖ A \cdot x - b ‖}_{2}^{2} + ε_{T V} {‖ x ‖}_{T V},

where the second term ||x||_TV denotes the so called Total Variation norm of x which is defined by the gradient magnitude of the image x

{‖ x ‖}_{T V} = \sum \sqrt{{(\nabla_{h} x)}^{2} + {(\nabla_{v} x)}^{2}},

where ∇_h and ∇_v are the horizontal and vertical finite difference operators. It has to be noted that the TV norm defined by Eq. (11) has the dimension of 1/length (assuming x dimensionless), so ε_TV has the dimension of length unlike the dimensionless ε_L₂. Therefore the ratio between the optimal ε_TV and ε_L₂ will in general depend on the choice of the pixel size. All values of ε_TV given in this article implicitly have the dimension of the pixel size that is used in the calculation.

In contrast with the L2-regularized Eq. (7) the second term of Eq. (10) promotes solutions with sparse gradient magnitude. Expression (10) represents a non-smooth convex minimization problem and can be calculated numerically using one of the iterative TV minimization algorithms [25]. During this study we used an implementation of TV minimization based on FISTA [29] which was introduced by Beck and Teboulle in [30].

The Beck and Teboulle algorithm can be viewed as an extension of dual gradient minimization [28]. Minimizing Eq. (10) is achieved by splitting each iteration into two sub-problems (or steps), the so called gradient step and denoising step:

Gradient step: finding an image x₀ by reducing the unregularized L2 residual term ${‖ A \cdot x - b ‖}_{2}^{2}$ in the beginning of each iteration,
Denoising step: followed by regularizing this intermediate image x₀ by minimizing: $\underset{x}{argmin} : {‖ x_{0} - x ‖}_{2}^{2} + ε_{T V} {‖ x ‖}_{T V} .$

Alternating gradient and denoising steps has been shown previously to speed up the convergence without sacrificing accuracy [31].

The gradient step is based on FISTA where, in order to reach a high rate of convergence, the current guess is updated using information from two previous iterations:

\begin{array}{l} 1 : & y^{n - 1} = x^{n - 1} + (\frac{t_{n - 1} - 1}{t_{n}}) (x^{n - 1} - x^{n - 2}); \\ 2 : & y^{n} = y^{n - 1} - \frac{2}{L} A^{T} \cdot (A \cdot y^{n - 1} - b); \\ 3 : & x^{n} = D_{L, ε} (y^{n}) . \end{array}

Here an intermediate vector y is introduced to take into account solution updates from two previous iterations, t_n is a scalar that is determined at each iteration as

t_{n} = \frac{1 + \sqrt{1 + 4 t_{n - 1}^{2}}}{2}

, t₀ = 1, and L denotes a so called Lipschitz constant that can be calculated as the maximum eigenvalue of the product A^* · A (see [29] for a detailed description). Operator D_L,ε signifies the denoising step which can be implemented using Fast Gradient Projection (FGP) [23]. During the denoising step, the TV norm of the current guess is minimized depending on the regularization weight ε_TV and the Lipschitz constant L.

If matrix Ã is diagonal, all frequency components of the solution ${\tilde{y}}_{j}^{n}$ can be updated during the gradient step (13) independently from each other:

{\tilde{y}}_{j}^{n} = {\tilde{y}}_{j}^{n - 1} - \frac{2 ω_{j}}{L} {\tilde{A}}_{j, j}^{*} ({\tilde{A}}_{j, j} {\tilde{y}}_{j}^{n - 1} - {\tilde{b}}_{j}) .

A frequency weighting vector ω_j ≤ 1 is introduced to control the convergence of the algorithm.

In general, convergence properties might vary among different minimization algorithms. However, the following observation is likely to be correct for algorithms similar to the one described above: frequency components of the solution x̃_j that correspond to small matrix elements |Ã_j,j| → 0 will only be modified in the denoising step. Hence, its final value is determined solely by the TV-term of Eq. (10). The latter ensures a significant difference between how TV minimization and L2 regularization (see Section 2.4) are computing the frequency components of the solution that can not be retrieved from the observations.

2.4. Linear phase retrieval algorithms: models

As was explained in the previous section, TV minimization can be used for phase retrieval as long as the phase-contrast model can be expressed as a linear system Eq. (6).

Let us introduce matrices C_D and S_D that represent the convolution with the real and imaginary parts of the Fresnel propagator P_D respectively. Their frequency domain representations have the following form:

{\tilde{C}}_{D} = (\begin{matrix} cos (α_{1}) & \emptyset \\ ⋱ \\ \emptyset & cos (α_{k}) \end{matrix}), {\tilde{S}}_{D} = (\begin{matrix} sin (α_{1}) & \emptyset \\ ⋱ \\ \emptyset & sin (α_{k}) \end{matrix})

where α_j = πλD|f_j|², λ stands for the X-ray wavelength and D is the effective propagation distance. Using matrices C_D and S_D we can construct matrix A and the corresponding update rule for the gradient step for the following linear phase retrieval models:

CTF model

The CTF model is widely used for phase retrieval in cases when the specimen yields negligible attenuation and slowly varying phase [12]. The Fourier transform of the phase-contrast image Ĩ_D(f) is approximated by the Fourier transforms μ̃ (f) and φ̃(f) of the projected attenuation and projected phase images of the specimen:

{\tilde{I}}_{D} (f) = δ (f) - 2 cos (α) \tilde{μ} (f) + 2 sin (α) \tilde{φ} (f) .

Linear systems that are formed by combining Eqs. (16) for a set of m phase-contrast images {I_D₍₁₎,..., I_D₍_m₎} can easily be represented in matrix form Eq. (6) (in frequency domain for |f_j| > 0) by construction of Ã, x̃ and b̃:

\tilde{A} = (\begin{matrix} - 2 {\tilde{C}}_{D (1)} & 2 {\tilde{S}}_{D (1)} \\ \dots \\ - 2 {\tilde{C}}_{D (m)} & 2 {\tilde{S}}_{D (m)} \end{matrix}), \tilde{x} = (\begin{matrix} \tilde{μ} \\ \tilde{φ} \end{matrix}), \tilde{b} = (\begin{matrix} {\tilde{I}}_{D (1)} \\ \dots \\ {\tilde{I}}_{D (m)} \end{matrix}) .

Here the unknown vectors μ̃ and φ̃ are concatenated into a single vector x̃, vector b̃ contains all observed images Ĩ_D_(i) and matrix Ã is obtained by concatenating pairs of matrices C̃_D_(i) and S̃_D_(i), where each pair corresponds to a particular propagation distance D(i). Using this representation we can find an update rule similar to Eq. (14) for each j^th frequency component of the unknown vector x̃. We will separate it into the update rules for μ̃ and φ̃ as follows:

\begin{array}{l} {\tilde{μ}}_{j}^{n} = {\tilde{μ}}_{j}^{n - 1} + \frac{4 ω_{j}}{L m} \sum_{i = 1}^{m} {\tilde{C}}_{D (i)}^{(j, j)} (- 2 {\tilde{C}}_{D (i)}^{(j, j)} {\tilde{μ}}_{j}^{n - 1} + 2 {\tilde{S}}_{D (i)}^{(j, j)} {\tilde{φ}}_{j}^{n - 1} - {\tilde{I}}_{D (i)}^{j}), \\ {\tilde{φ}}_{j}^{n} = {\tilde{φ}}_{j}^{n - 1} - \frac{4 ω_{j}}{L m} \sum_{i = 1}^{m} {\tilde{S}}_{D (i)}^{(j, j)} (- 2 {\tilde{C}}_{D (i)}^{(j, j)} {\tilde{μ}}_{j}^{n - 1} + 2 {\tilde{S}}_{D (i)}^{(j, j)} {\tilde{φ}}_{j}^{n - 1} - {\tilde{I}}_{D (i)}^{j}) . \end{array}

Calculation of Eq. (18) has to be carried out at each iteration of the gradient and denoising steps. Since the denoising step must be computed in the spatial domain, two inverse Fourier transforms (for μ̃ and φ̃) must be calculated at each iteration before the denoising step and two Fourier transforms after the denoising step.

Mixed model

The Mixed model [14] is used for phase retrieval in cases with (significant) attenuation. In an approximated version of this model (assuming only the first two terms), the phase-contrast image Ĩ_D(f) becomes:

{\tilde{I}}_{D} (f) = cos (α) {\tilde{I}}_{0} (f) + 2 sin (α) (\tilde{I_{0} φ}) (f) .

Here Ĩ₀(f) denotes the Fourier transform of the intensity image at zero distance. Ĩ₀(f) is fully determined by the attenuation image of the specimen and can be expressed in the spatial domain as I₀ = e^−2μ. The linear system that describes a set of phase-contrast images {Ĩ_D₍₁₎,..., Ĩ_D_(m)} based on Eq. (19) can be expressed through the following Ã, x̃ and b̃:

A = (\begin{matrix} {\tilde{C}}_{D (1)} & 2 {\tilde{S}}_{D (1)} \\ \dots \\ {\tilde{C}}_{D (m)} & 2 {\tilde{S}}_{D (m)} \end{matrix}), \tilde{x} = (\begin{matrix} {\tilde{I}}_{0} \\ \tilde{I_{0} φ} \end{matrix}), \tilde{b} = (\begin{matrix} {\tilde{I}}_{D (1)} \\ \dots \\ {\tilde{I}}_{D (m)} \end{matrix}) .

Here we treat the element-wise product of the intensity image and the phase image (I₀φ) as an unknown image independent from the unknown intensity image I₀. Using this representation we can write down the update rule Eq. (14) for the gradient step as:

\begin{array}{l} {\tilde{I}}_{0, j}^{n} = {\tilde{I}}_{0, j}^{n - 1} - \frac{2 ω_{j}}{L m} \sum_{i = 1}^{m} {\tilde{C}}_{D (i)}^{(j, j)} ({\tilde{C}}_{D (i)}^{(j, j)} {\tilde{I}}_{0, j}^{n - 1} + 2 {\tilde{S}}_{D (i)}^{(j, j)} {(\tilde{I_{0} φ})}_{j}^{n - 1} - {\tilde{I}}_{D (i)}^{j}) \\ {(\tilde{I_{0} φ})}_{j}^{n} = {(\tilde{I_{0} φ})}_{j}^{n - 1} - \frac{4 ω_{j}}{L m} \sum_{i = 1}^{m} {\tilde{S}}_{D (i)}^{(j, j)} ({\tilde{C}}_{D (i)}^{(j, j)} I_{0, j}^{n - 1} + 2 {\tilde{S}}_{D (i)}^{(j, j)} {(\tilde{I_{0} φ})}_{j}^{n - 1} - {\tilde{I}}_{D (i)}^{j}) . \end{array}

Phase-attenuation duality models

These models can be used when the specimen has a homogeneous composition or, in the limited range of X-ray energies, when the specimen is composed of light elements [17, 19]. In duality models the phase and attenuation images of a specimen are assumed to be proportional to each other, i.e. $σ = \frac{φ}{μ}$ , permitting to reduce the number of unknown variables. In this study we consider two duality models:

{\tilde{I}}_{D} (f) = 2 (σ sin (α) - cos (α)) \tilde{μ} (f),

which was derived from the CTF model [18] and :

{\tilde{I}}_{D} (f) = (cos (α) + (α - σ) sin (α)) {\tilde{I}}_{0} (f),

which was derived from the Mixed model [19]. Both models can be used to retrieve the projected attenuation image of the specimen based on a single phase contrast image Ĩ_D(f) acquired at a suitable distance D > 0 as well as using a set of m phase-contrast images {I_D₍₁₎,..., Ĩ_D_(m)} recorded at different distances. The update rule for the gradient step based on Eq. (22) is then simplified into:

{\tilde{μ}}_{j}^{n} = {\tilde{μ}}_{j}^{n - 1} - \frac{4 ω_{j}}{L m} \sum_{i = 1}^{m} {\tilde{B}}_{i} (2 {\tilde{B}}_{i} {\tilde{μ}}_{j}^{n - 1} - {\tilde{I}}_{D (i)}^{j}),

where

{\tilde{B}}_{i} = (σ {\tilde{S}}_{D (i)}^{(j, j)} - {\tilde{C}}_{D (i)}^{(j, j)})

. The update rule for the gradient step based on Eq. (23) has the following form:

{\tilde{I}}_{0, j}^{n} = {\tilde{I}}_{0, j}^{n - 1} - \frac{2 ω_{j}}{L m} \sum_{i - 1}^{m} {\tilde{B}}_{i} ({\tilde{B}}_{i} {\tilde{I}}_{0, j}^{n - 1} - {\tilde{I}}_{D (i)}^{j}),

where

{\tilde{B}}_{i} = ({\tilde{C}}_{D (i)}^{(j, j)} + (α - σ) {\tilde{S}}_{D (i)}^{(j, j)})

.

3. Simulations

As indicated before, TV minimization permit an accurate solution for a class of inverse problems based on severely incomplete sets of observations. The underlying assumption of all TV minimization approaches is that the unknown signal must have a sparse gradient magnitude. In the field of image reconstruction such assumption is fulfilled when the reconstructed image is piecewise constant, i.e. the intensity only changes in a step-like manner.

3.1. Phantom image with sparse gradient magnitude

Our first demonstration of phase retrieval based on TV minimization uses a ’flat’ piece-wise constant phantom. Here TV minimization is expected to yield high accuracy. However, in the following subsections we will abandon this restriction and use a ’spheres’ phantom to demonstrate the performance of TV minimization in more realistic cases.

Figure 1 shows the comparison between L2-regularized and TV-regularized phase retrieval for the CTF model. The ground truth projected attenuation and phase images (256 × 256 pixels) were computed for a randomly generated composition of overlapping polyethylene disks immersed in a layer of water with a total thickness d = 0.1 mm (Fig. 1(a)). Subsequently phase-contrast images I_D₁ and I_D₂ are generated using Fresnel propagation for a monochromatic X-ray energy of 20KeV and a pixel size of 1 μm (Fig. 1(b, c)). In the current simulation we used propagation distances {0 m, 1 m}. Gaussian noise with standard deviation of 0.02 was added to both images mimicking acquisition with poor SNR. The optimal regularization parameters ε_L₂ and ε_TV were chosen such that the overall Root Mean Square Error (RMSE) is minimal for the listed conditions (see Section 3.3). In practice, they must be derived from an estimate of the SNR of the measured images. No stopping criteria was used in the gradient step of the TV minimization, instead all reconstructions were computed using 1000 iterations in order to guarantee convergence.

Fig. 1 Phase reconstructions based on the simulated data of the ’flat’ phantom for TV and L2 regularization with and without frequency weighting (see Sec. 3.1 for simulation parameters). (a) Ground truth. (b) Intensity image at zero distance with Gaussian noise. (c) Propagated phase-contrast image at 1m with Gaussian noise. (d) L2-regularized solution, no frequency weighting. (e) TV-regularized solution, no frequency weighting. (f) L2-regularized solution, with frequency weighting. (g) TV-regularized solution with frequency weighting.

Download Full Size | PDF

The resulting solutions are shown in Fig. 1(d, e). The frequency spectrum (we will use this term for the angular average of the magnitude of 2D Fourier transform of the image) of the error magnitude of the solutions can be seen on Fig. 2 (right). The frequency spectrum of the error was normalized against the frequency spectrum of the ground truth image. It is evident from this graph that when the optimal regularization weights are used, TV regularization is significantly more accurate then L2 regularization at all frequencies where the direct phase retrieval problem is undetermined (for frequencies components with |Ã_j,j| → 0). Note that TV minimization can yield higher error within particular bands of frequencies to promote a sparse gradient. We have already pointed out in Section 2.2 and Section 2.3 that there are major differences in how TV minimization treats frequency components of the solution that correspond to |Ã_j,j| → 0. In order to confirm this, we compared the results of the TV method with ε_TV = 0.02 and ε_TV = 0 (no TV minimization) when applied to noise-free data. Figure 2(left) shows that the frequency spectrum of the normalized difference between the regularized and non-regularized solutions. It is evident that the regularization strength depends heavily on |Ã_j,j|.

Fig. 2 The radial frequency spectrum (angular averaged) of the reconstruction. Dotted line shows |Ã_j,j| on both graphs. Left: normalized difference between the solutions with and without TV-regularization, i.e. (ε_TV = 0.02) vs. (ε_TV = 0). Right: normalized error-magnitude of the TV-regularized solution (solid line) and the L2-regularized solution (dashed line). The graphs are normalized against the radial frequency spectrum of the ground truth image.

Download Full Size | PDF

So far we used scalar regularization weights ε_L₂ and ε_TV in both regularization approaches. One can achieve a significantly improved accuracy by using a frequency dependent regularization weighting instead. In L2 regularization, the scalar ε_L₂ can be replaced by the frequency dependent factor ε_L2,j ∼ SNR⁻¹(f_j) (i.e. Wiener deconvolution). In our investigation we assumed the SNR of the reconstructed image to be proportional to its frequency spectrum. The latter can be estimated from a preliminary reconstruction based on a scalar regularization weight using either L2 or TV regularization. Such estimation is demonstrated in Fig. 3(left). The estimate of the image SNR was introduced into TV minimization approach using the frequency weighting vector ω_j ∼ SNR(f_j) in Eq. (14). Resulting L2- and TV-regularized solutions are depicted in Fig. 1(f) and Fig. 1(g) respectively. It can be seen in Fig. 3(right) that the error of TV-regularized solution is significantly lower than the error of L2-regularized solution in a wide frequency range.

Fig. 3 Radial frequency spectrum (angular averaged) of the retrieval results using the frequency dependent regularization weights (i.e. frequency weighting). Dotted line shows |Ã_j,j| on both graphs. Left: estimated frequency spectrum of the specimen (solid line). Right: normalized error-magnitude of the TV-regularized solution (solid line) and L2-regularized solution (dashed line) with frequency weighting.

Download Full Size | PDF

3.2. Realistic phantom

In order to demonstrate the performance of the TV-minimization in realistic cases we tested L2-regularization and TV-minimization phase retrieval approaches for phantoms with a non-sparse gradient. Figure 4 shows reconstructions of the projected phase image of a composition of randomly sized and positioned polyethylene spheres immersed in water (the rest of simulation parameters match those from Section 3.1). Reconstructed images are obtained using both L2- and TV-regularization based on the CTF model, the Mixed model as well as their phase-attenuation duality modifications. In phase retrieval based on duality models we have used a single simulated phase-contrast image with propagation distance 1 m. The corresponding normalized frequency spectra of the error-magnitude are depicted in Fig. 5. It is evident that the TV-regularized solutions yield lower error in comparison with the L2-regularized ones in a broad range of frequencies. It is also apparent that given the parameters used in the current simulation (low attenuation and moderate phase changes), reconstructions obtained with the CTF model and the Mixed model are virtually indistinguishable, while there is some discrepancy with the models that are based on phase-attenuation duality assumption.

Fig. 4 Reconstructions based on the simulated data of the ’spheres’ phantom for different propagation models. All results are obtained using frequency weighting. Top row of images (a, c, e, g) show L2-regularized solutions. The bottom row of images (b, d, f, h) show TV-regularized solutions. Solutions (a, b) are based on the CTF model, (c, d) on the Mixed model, (e, f) on the dual-CTF model and (g, h) on the dual-Mixed model.

Download Full Size | PDF

Fig. 5 The radial frequency spectrum (angular averaged) of the reconstructions for different propagation models. Dotted line shows |Ã_j,j| on all graphs. Normalized error magnitude of the TV-regularized solution (solid line) and L2-regularized solution (dashed line).

Download Full Size | PDF

3.3. Optimal regularization weights

The accuracy of both regularization methods considered in this paper strongly depend on the choice of the regularization weights ε_L₂ and ε_TV. In practice, these parameters are estimated from the measured data or chosen in some heuristic manner. In order to investigate how the choice of the regularization weight affects the accuracy of phase retrieval, we measured the total RMSE of the reconstructed phase images varying two parameters: standard deviation (STD) of the Gaussian noise and the regularization weight of the phase retrieval algorithm. The phase-contrast images were simulated using the polyethylene spheres phantom using the simulation parameters described in Section 3.1 and 3.2. Fig. 6(left) shows the values of RMSE for phase retrieval based on the CTF model. Curves L2-1, L2-2, L2-3 and L2-4 represent the RMSE of the L2-regularized solutions with noise STD = 0.01, 0.02, 0.05 and 0.1 respectively. Curves TV-1, TV-2, TV-3 and TV-4 represent the RMSE of the TV-regularized solutions with corresponding noise STD. It can be seen that TV-regularization yields solutions with a 4–7 times lower total RMSE in comparison to those obtained with L2 regularization. It is also evident that estimation of the optimal regularization parameter is more important for TV-regularization then for L2 regularization as it affects the RMSE to a larger extent.

Fig. 6 Influence of the regularization weight on samples with increasing noise levels and on samples of increasing thickness. RMSE of the solution is shown against the magnitude of regularization parameter for L2 regularization (crosses) against TV regularization (plus signs). Left: each curve shows errors for different noise levels, STD = 0.01, 0.02, 0.05, 0.1. Right: each curve shows errors for different specimen thickness, 0.1 mm, 0.2 mm, 0.5 mm, 1 mm.

Download Full Size | PDF

Along with various types of additive noise, non-linear effects that are not taken into account by the phase retrieval models can become an important source of errors. The fraction of signal to error due to non-linearity of the CTF model is, in a certain range of conditions, proportional to |sin(α)| [20]. That property allows to treat the non-linearity error as another form of noise. Figure 6(right) demonstrates the values of RMSE of the phase retrieval based on the CTF model applied to noise-free phase-contrast images simulated for phantoms of different thicknesses. Curves L2-1, L2-2, L2-3 and L2-4 show the RMSE of the L2-regularized phase retrieval for phantoms with a total thickness = 0.1 mm, 0.25 mm, 0.5 mm and 1 mm respectively. Curves TV-1, TV-2, TV-3 and TV-4 represent the RMSE of the corresponding TV-regularized solutions. The thickness of the phantom has a dramatic effect on the accuracy of the phase retrieval based on the CTF model since it introduces larger variations in projected attenuation and phase images which leads to greater non-linearity of the observed phase-contrast image. It is evident that TV-regularized phase retrieval yields similar accuracy with the L2-regularized solution in the case of thin phantom. However the advantage of TV regularization becomes significant for thicker phantoms.

4. Experiment

To test the TV-regularized phase retrieval on experimental data we have used X-ray phasecontrast images of a test pattern designed to assess the resolution of the X-ray imaging system. The pattern consisted of 700 nm high lithographic gold structures on top of a Si substrate. Phase-contrast data was acquired at the beamline ID22NI of the European Synchrotron Radiation Facility (Grenoble, France). The incoming x-rays were focused using Kirkpatrick-Baez mirror system, which gave a point focus with sub 100 nm with and height (FWHM) [32]. The mean energy of the x-rays was 17.5 KeV (∼ 1.5% bandwidth) with 10¹² photons/s in the focus. The focus was used as a point source for projection microscopy [33] giving effective pixel size of 53 nm. Images were acquired at four propagation distances {27.4 mm, 28.3 mm, 31.8 mm, 40.3 mm}.

Figure 7 shows reconstructed phase images of the lines and dots pattern and so called Siemens star pattern. The L2-regularized solutions are shown in the top row, while the bottom row shows the TV-regularized solutions. Both phase-retrieval approaches were based on phase-attenuation duality CTF model. The reconstructions depicted in sub-figures (a, b, e and f) are based on phase-contrast images acquired at four different distances. The reconstructions from sub-figures (c, d, g and h) are based on a single phase-contrast image acquired at a propagation distance 27.4 mm.

Fig. 7 Phase retrieval from experimental data using the dual-CTF model. (a–d) Phase retrieval of a ’dots and lines’ pattern. (e–h) Phase retrieval of a ’star’ pattern. (a, e) L2-regularized solution based on 4 images recorded at different propagation distances. (b, f) TV-regularized solution based on 4 images recorded at different propagation distances. (c, g) L2-regularized solution based on a single recorded image. (d, h) TV-regularized solution based on a single recorded image.

Download Full Size | PDF

It is evident that TV regularization permits a very high quality reconstruction of lines and dots pattern based only on a single phase-contrast image. The high frequency ripple-like artifacts are efficiently suppressed while the pattern is accurately reconstructed. That can be explained by the fact that the reconstructed image has a sparse gradient magnitude and that spatially localized structures such as dots and lines have a broad footprint in frequency space. The problem of phase retrieval based on a single phase-contrast image is ill-posed within a set of spatial frequencies that depends on the acquisition parameters. TV regularization allows to fill-in these gaps in frequency space by applying the constraint of sparse gradient magnitude. Phase retrieval of the Siemens star pattern shows that the structures which are periodic in space are harder to reconstruct. Their footprint in frequency space is localized and might be completely irrecoverable from a single phase-contrast image with the given acquisition parameters.

5. Conclusion

Phase retrieval in propagation based X-ray PCI can be improved using iterative TV minimization algorithms. Reconstructions based on simulated and experimental data show that phase retrieval based on TV minimization can significantly outperform the current practice, a deconvolution approach with L2 regularization. TV minimization can be used with different linear phase retrieval models including the CTF model, the Mixed model and the phase-attenuation duality models. Although the method works best for specimen that adhere to the constraints, the method also improves the phase reconstruction for specimen that are not exactly sparse in their gradient magnitude. TV minimization provides an effective regularization instrument for solving an underdetermined linear systems of equations. Analysis of the Fourier spectrum of the error of the reconstructed phase images clearly demonstrates that TV minimization allows partial recovery of the solution within the frequency bands that are undefined by the particular phase-contrast model or corrupted by noise. That feature of TV minimization allows effective suppression of the high frequency artifacts in single-distance phase retrieval based on phase-attenuation duality models and permits more accurate reconstruction of low spatial frequencies in multi-distance approaches.

TV minimization finds the solution in the form of a TV-regularized least-squares fit and it does not require the knowledge of the attenuation part of the specimen. A frequency dependent estimation of the signal-to-noise ratio can be used for each observed image, dramatically improving the accuracy of reconstruction. Simulations show that TV regularization can suppress errors that occur both due to additive noise and due to non-linearity of the phase propagation. Experimental data has shown that TV minimization can also significantly improve the accuracy of phase reconstruction of real specimen that comply with gradient sparsity condition and imaged under realistic circumstances. These results could allow to decrease the number of phase-contrast images that are needed for the accurate image reconstruction in some applications of X-ray PCI. This is particularly important for in situ experiments or to reduce radiation damage to the specimen.

Acknowledgments

This work was partially supported by the Care4U project with financial support of Point One, an innovation program of the Ministry of Economic Affairs in The Netherlands.

References and links

1. O. Coindreau, C. Mulat, C. Germain, J. Lachaud, and G. L. Vignoles, “Benefits of x-ray CMT for the modeling of C/C composites,” Adv. Eng. Mat. 13, 178–185 (2011). [CrossRef]

2. F. Cosmi, A. Bernasconi, and N. Sodini, “Phase contrast micro-tomography and morphological analysis of a short carbon fibre reinforced polyamide,” Compos. Sci. Technol. 71, 23–30 (2011). [CrossRef]

3. M. Herbig, A. King, P. Reischig, H. Proudhon, E. M. Lauridsen, J. Marrow, J. -Y. Buffire, and W. Ludwig, “3-D growth of a short fatigue crack within a polycrystalline microstructure studied using combined diffraction and phase-contrast x-ray tomography,” Acta Mater. 59, 590–601 (2011). [CrossRef]

4. A. Kostenko, H. Sharma, E. G. Dere, A. King, W. Ludwig, W. van Oel, S. Stallinga, L. J. van Vliet, and S. E. Offerman, “Three-dimensional morphology of cementite in steel studied by x-ray phase-contrast tomography,” Scr. Mater. 67, 261–264 (2012). [CrossRef]

5. T. Argunova, M. Gutkin, J. H. Je, E. Mokhov, S. Nagalyuk, S. Sergey, and Y. Hwu, “SR phase-contrast imaging to address the evolution of defects during SiC growth,” Phys. Status Solid. A Appl. Mat. 208, 819–824 (2011). [CrossRef]

6. J. Xu, A. W. Stevenson, D. Gao, M. Tykocinski, D. Lawrence, S. W. Wilkins, G. M. Clark, E. Saunders, and R. S. Cowan, “The role of radiographic phase-contrast imaging in the development of intracochlear electrode arrays,” Otology and Neur. 22, 862–868 (2001). [CrossRef]

7. P. Weiss, L. Obadia, D. Magne, X. Bourges, C. Rau, T. Weitkamp, I. Khairoun, J.M. Bouler, D. Chappard, O. Gauthier, and G. Daculsi, “Synchrotron x-ray microtomography (on a micron scale) provides three dimensional imaging representation of bone ingrowth in calcium phosphate biomaterials,” Biomat. 24, 4591 (2003). [CrossRef]

8. P. Coan, F. Gruener, C. Glaser, T. Schneider, A. Bravin, M. Reiser, and D. Habs, “Phase contrast medical imaging with compact x-ray sources at the Munich-centre for advance photonics,” Nucl. Instr. and Meth. Phys. Res. A 608, S44–S46 (2009). [CrossRef]

9. E. C. Ismaila, W. Kaabara, D. Garritya, O. Gundogdua, O. Bunkb, F. Pfeifferb, M. J. Farquharsond, and D. A. Bradleya, “X-ray phase contrast imaging of the bone-cartilage interface,” Appl. Rad. and Isot. 68, 767–771 (2010). [CrossRef]

10. Q. Tao, D. li, L. Zhang, and S. Luo, “Using x-ray in-line phase-contrast imaging for the investigation of nude mouse hepatic tumors,” PLoS ONE 7, e39936, (2012). [CrossRef] [PubMed]

11. A. A. Appel, J. C. Larson, S. Somo, Z. Zhong, P. P. Spicer, F. K. Kasper, A. B. Garson III, A. M. Zysk, A. G. Mikos, M. A. Anastasio, and E. M. Brey, “Imaging of poly(α-hydroxy-ester) scaffolds with x-ray phase-contrast microcomputed tomography,” Tissue Eng. Part C 18, (2012).

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

13. T.E. Gureyev, T.J. Davis, A. Pogany, S.C. Mayo, and S.W. Wilkins, “Optical phase retrieval by use of first Born- and Rytov-type approximations,” Appl. Opt. 43, 2418–2430 (2004). [CrossRef] [PubMed]

14. X. Wu and H. Liu, “A general theoretical formalism for x-ray phase contrast imaging,” J. of X-ray Sci. and Techn. 11, 33–42 (2003).

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

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

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

18. L. Turner, B. Dhal, J. Hayes, A. Mancuso, K. Nugent, D. Paterson, R. Scholten, C. Tran, and A. Peele, “X-ray phase imaging: demonstration of extended conditions with homogeneous objects,” Opt. Express 12, 2960–2965 (2004). [CrossRef] [PubMed]

19. X. Wu and A. Yan, “Phase retrieval from one single phase contrast x-ray image,” Opt. Express 17, 11187 (2009). [CrossRef] [PubMed]

20. R. Hofmann, J. Moosmann, and T. Baumbach, “Criticality in single-distance phase retrieval,” Opt. Express 19, 25881–25890 (2011). [CrossRef]

21. M. Langer, P. Cloetens, and F. Peyrin, “Fourier-wavelet regularization of phase retrieval in x-ray in-line phase tomography,” J. Opt. Soc. Am. A 26, 1876–1881 (2009). [CrossRef]

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

23. A. Chambolle, “An algorithm for total variation minimization and applications,” J. Math. Imaging Vision 20, 89–97 (2004). [CrossRef]

24. E. Y. Sidky, M. A. Anastasio, and X. Pan, “Image reconstruction exploiting object sparsity in boundary-enhanced x-ray phase-contrast tomography,” Opt. Express 18, 10404 (2010). [CrossRef] [PubMed]

25. J. Dahl, P. C. Hansen, S. H. Jensen, and T. L. Jensen, “Algorithms and software for total variation image reconstruction via first-order methods,” Num. Alg. 53, 67–92 (2010). [CrossRef]

26. A. W. M. van Eekeren, K. Schutte, and L. J. van Vliet, “Multiframe super-resolution reconstruction of small moving objects,” IEEE Trans. Im. Proc. 19, 2901–2912 (2010). [CrossRef]

27. E. Candes, J. Romberg, and T. Tao, “Robust uncertainty principles: exact signal reconstruction from highly incomplete frequency information,” IEEE Trans. Inform. Theory 52, 489–509 (2006). [CrossRef]

28. X. Bresson and T. F. Chan, “Fast dual minimization of the vectorial total variation norm and applications to color image processing,” Inverse Probl. and Imag. 2, 455–484 (2008). [CrossRef]

29. A. Beck and M. Teboulle, “A fast iterative shrinkage-thresholding algorithm for linear inverse problems,” SIAM J. Im. Sci. 2, 183–202 (2009). [CrossRef]

30. A. Beck and M. Teboulle, “Fast gradient-based algorithms for constrained total variation image denoising and deblurring problems,” IEEE Trans. Im. Proc. 18, 2419–2434 (2009). [CrossRef]

31. T. Q. Pham, L. J. van Vliet, and K. Schutte, “Robust super-resolution by minimizing a gaussian-weighted l2 error norm,” J. Phys.: Conf. 124, 012037 (2008). [CrossRef]

32. R. Barrett, R. Baker, P. Cloetens, Y. Dabin, C. Morawe, H. Suhonen, R. Tucoulou, A. Vivo, and L. Zhang. “Dynamically-figured mirror system for high-energy nanofocusing at the ESRF, Proc. SPIE 12, 813904 (2011). [CrossRef]

33. R. Mokso, P. Cloetens, E. Maire, W. Ludwig, and J.-Y. Buffire, “Nanoscale zoom tomography with hard x-rays using Kirkpatrick-Baez optics, Appl. Phys. Lett. 90, 144104 (2007). [CrossRef]

Phase retrieval in in-line x-ray phase contrast imaging based on total variation minimization

Abstract

1. Introduction

2. Materials and methods

2.1. Matrix formalism for phase propagation model

2.2. Linear phase retrieval algorithms: L2-norm regularization

2.3. Linear phase retrieval algorithms: TV minimization

2.4. Linear phase retrieval algorithms: models

CTF model

Mixed model

Phase-attenuation duality models

3. Simulations

3.1. Phantom image with sparse gradient magnitude

3.2. Realistic phantom

3.3. Optimal regularization weights

4. Experiment

5. Conclusion

Acknowledgments

References and links

Cited By

Figures (7)

Equations (25)

Optics Express