Fast iterative reconstruction of differential phase contrast X-ray tomograms

: Differential phase-contrast is a recent technique in the context of X-ray imaging. In order to reduce the specimen’s exposure time, we propose a new iterative algorithm that can achieve the same quality as FBP-type methods, while using substantially fewer angular views. Our approach is based on 1) a novel spline-based discretization of the forward model and 2) an iterative reconstruction algorithm using the alternating direction method of multipliers. Our experimental results on real data suggest that the method allows to reduce the number of required views by at least a factor of four.


Introduction
X-ray phase-contrast imaging (PCI) is a promising alternative to absorption-based computed tomography (CT) for visualizing many structures in biological samples, in particular soft-tissues. PCI is based on the phase shift induced by the propagation of a coherent wave through the investigated object. Various PCI methods have been developed including analyzer based [1][2][3], interferometric [4][5][6] and free space propagation methods [7][8][9]. These methods differ substantially in terms of the physical signal that is measured and the required experimental setup.
Weitkamp et al [6] and Momose et al [10] recently developed a new X-ray phase-contrast imaging method based on grating interferometry (GI). It has attracted increasing interest in a variety of fields owing to its unique combination of imaging characteristics. First, GI provides a high sensitivity to electron-density variations, down to 0.18e/nm 3 [11]. This makes the technique particularly suitable for soft-tissue specimens. It has been applied successfully to biological samples such as insects [6,11], rat-brain tissue and even human as breast tissue [12]. Second, GI produces three complementary types of information; namely, attenuation, phase-shift and dark-field measurements. In differential phase-contrast imaging (DPCI), one focuses exclusively on the phase information, which in principle allows for reconstructing the refractive index distribution of the object. Third, GI does not require a highly monochromatic source, which means that conventional laboratory X-ray tubes can be used. The combination of the aforementioned characteristics makes GI suitable for a broad range of applications, such as material sciences (e.g., material testing), biomedical research (e.g., monitoring drug effects), or even clinical diagnostics (e.g., mammography).
DPCI essentially yields the derivative of the Radon transform of the object's refractive-index map. This map can therefore be retrieved using a suitable variant of the filtered back-projection (FBP) algorithm known from conventional tomography. However the FBP method requires a large number of viewing angles, and thus the acquisition time is very long. The long radiation time can also damage biological samples. This provides a strong motivation for developing algorithms that can handle fewer views, so as to reduce the exposure time.

Contributions
To this end, it is natural to formulate the reconstruction as an inverse problem that is regularized using prior information on the sample. The first step in such a formulation is the discretization of the forward model. The contributions of this paper are summarized as follows: • The proposal of a new discretization scheme for derivatives of the Radon transform based on B-spline calculus. We show that the model lends itself to an analytical treatment, yielding a fast and accurate implementation.
• The design of an iterative reconstruction algorithm taking advantage of the proposed discretization and using a combination of TV and Tikhonov regularization. Our method follows an augmented-Lagrangian optimization principle and makes use of the conjugate gradient method to solve the linear step in the alternating direction method of multipliers (ADMM).
• The application of a problem-specific preconditioner that speeds up the convergence of the linear optimization step considerably.
• The experimental evaluation of the method using real data.
An early version of the proposed algorithm was presented at ISBI 2012 [13].

Related work
While ADMM has been adapted for several imaging modalities, e. g., magnetic resonance imaging [14] and CT [15], it has not yet been used in the context of DPCI. In fact, only a small number of contributions have proposed to apply iterative reconstruction techniques to DPCI [16][17][18]. While these techniques also rely on TV regularization, they differ from the present work in the way the forward model is discretized and in terms of the optimization algorithm. We now briefly review the characteristics of the aforementioned works. Qi et al. [16] express the derivative in the projection domain using derivatives along the two axes in the image domain. Then they implement these derivatives with finite differences. A basic steepest descent algorithm is used to solve the corresponding TV problem.
Kohler et al. [17] and Xu et al. [18] discretize the forward model using Kaiser-Bessel window functions. The advantage of these functions is that the derivative of their Radon transform is known analytically. However, they fail to satisfy the partition of unity, which is a necessary condition for controlling the error of approximation [19]. Note that [18] also presents results obtained with rectangular basis functions (splines of degree zero). However, the Radon transform of these functions is not differentiable analytically. Thus the authors have to use a numerical approximation of the derivative. In contrast, our work is based on higher-degree splines, which allows for fully analytic derivations. Concerning the optimization algorithm, [17] uses the Newton-Raphson method with a separable-paraboloidal-surrogate, while [18] relies on an adaptive-steepest-descent algorithm combined with projections onto convex sets [20].

Outline
In Section 2, we briefly review the physics of differential phase-contrast imaging. The endresult is that DPCI can be described mathematically in terms of derivative variants of the Radon transform. We discuss the properties of these operators in Section 3. We then propose a discretization of the forward model based on B-spline functions in Section 4 as well as a fast and accurate implementation. In Section 5, we describe our reconstruction scheme based on ADMM to solve regularized reconstruction problem. The framework and the proposed reconstruction method are validated in Section 6 through experiments with real data. We summarize and conclude our work in Section 7.

Physical model
An X-ray plane wave can be characterized by its intensity and phase. However the intensity is the only directly measurable part. Therefore, to extract the phase it is necessary to map this information into intensity patterns. In DPCI, this is achieved by using grating interfereometers. The principle has been described in detail in [6,21,22], and we briefly review it here. Figure 1 shows the physical setup of DPCI. It consists of an object on a rotation stage, a phase grating (G1) and an analyzer absorption grating (G2) which are behind the sample. Typically the phase grating produces a phase shift of π on the incident wave.
In order to specify the physical model of DPCI, we first need to set the geometry: θ θ θ = (cos θ , sin θ ) is a unit vector that is orthogonal to the projection direction and y is the projection coordinate depicted in Figure 1. The spatial coordinate of the 2-D object is denoted by The phase shift induced by the object on the transmitted beam is given by where is the real part of the object's refractive index, θ is the angle of view of the monochromatic wave illuminating the object and R{ f (x x x)}(y, θ ) is the Radon transform of f (x x x).
) maps a function on R 2 into the set of its integrals over the lines of R 2 . Specifically, if θ θ θ = (cos θ , sin θ ) and y ∈ R, then where θ θ θ ⊥ is perpendicular to θ θ θ .  The phase grating introduces a phase shift of π in the transmitted wave. The absorption grating is necessary for measuring the received wave given the limited resolution of the detector. The diagram on the right shows the received intensity for one pixel of the CCD camera with respect to the position of the grating. The red curve corresponds to the situation where the object is in the illumination path. The black curve shows the received intensity when there is no object.
The object introduces changes in the propagation direction of the illumination wave. The refraction angle is proportional to the derivative of the phase of the output wave with respect to y: where α(y, θ ) is the refraction angle in radians, and λ is the wavelength After passing through the object, the X-ray reaches the phase grating which introduces a periodic phase modulations. It essentially splits the beam into its first two diffraction orders. A periodic interference pattern perpendicular to the optical axis is formed. To measure this pattern with high resolution, one then uses a phase stepping technique (PST) [6]. To that end, an absorption grating is placed in front of the detector with the same periodicity and orientation as the pattern. In this technique, the absorption grating is moved perpendicular to the optical axis and the intensity signal in each pixel in the detector plane is recorded as a function of the grating position, x g .
When an object is placed on the optical path, the illumination wave is refracted which causes a local shift of the interference pattern. This displacement is given by where d is the distance between the two gratings. Combining this relation with (1) and (3) yields For our purpose, this equation characterizes the forward model of DPCI: f is the quantity of interest and g is the physically measurable signal.

Review of the derivative variants of the Radon transform and generalized filtered back-projection
Since the measurements are related to the Radon transform and its derivatives, we start by establishing some higher-level mathematical properties of these operators. The nth derivative of the Radon transform of a function f (x x x) is denoted by The derivatives of the Radon Transform are linear operators with the following properties: • Scale invariance: • Pseudo-distributivity with respect to convolution: • Projected translation invariance: To derive the necessary relations for the rest of the paper, we define a new operator, the Hilbert transform along the second coordinate.

Proposition 1.
The sequential application of the Radon transform, the nth derivative operator and the adjoint of the Radon transform on function where (− ) 1 2 is the fractional Laplace operator with transfer function ||ω ω ω|| and (− ) n 2 is n times application of this operator.
(11) is a key equation. It implies that if one is interested in solving the inverse problem the direct solution is to first apply the Radon adjoint on the measured data and then apply the inverse of the operator 2π (−1) n H n 2 (− ) (n−1)/2 . The transfer function of this inverse is 1 2π i n sgn(ω 2 ) n ω ω ω −(n−1) . An equivalent form of (11) using the fact that where R (n) is the adjoint of the n-th derivative of the Radon transform and the transfer function of q(y) is: (17) is the basis for the generalized filtered back projection (GFBP). The full procedure is described in Algorithm 1.
Filter the input data with the transfer function: 4: end for 5: Apply the adjoint of the n-th derivative of the Radon transform on the output of the previous stage.

Discretization scheme and implementation of the forward model
In this section we define discrete versions of the Radon transform and its derivatives which are consistent with the continuous-domain operators. This is the basis for a fast and accurate implementation of the DPCI forward model.

Discrete model of the Radon transform and its derivatives
We propose a formulation that relies on the practical generalized sampling scheme [19]. Let ϕ(x x x) ∈ L 2 (R 2 ) denote a generating Riesz basis function for the approximation space V as: The orthogonal projection of the object Using the linearity and translation invariance properties of the Radon transform, the Radon where and R{ϕ}(y, θ ) = F −1 { ϕ(ωθ θ θ )}(y). Likewise, we obtain the derivatives of the Radon transform: where Based on (22) and (24) we need to choose an appropriate generating function and to determine its Radon transform analytically. In this work we use a tensor-product B-spline, which provides a good compromise between support size and approximation order [24]: Here β m is a centered univariate B-spline of degree m: where

Proposition 2.
The n-th derivative of the Radon transform of the tensor-product B-spline is given by Proof. We define g(y, θ ) = R (n) {β }(y, θ ). By an application of the Fourier-slice theorem let us now focus on the Fourier-domain expression Taking the inverse Fourier transform of this expression yields the desired result.
To discretize our forward model, we need to calculate the derivatives of the Radon transform of f at the location (y j , θ i ) where y j = jΔy and θ i = iΔθ The values R (n) ϕ k k k (y j , θ i ) are then determined using Proposition 2.

Matrix formulation of the forward model
To describe the forward model in a more concise way, we now introduce an equivalent matrix formulation of (30) g = Hc , where c is a vector whose entries are the B-spline coefficients c k k k , g is the output vector and H is the system matrix defined by: The corresponding adjoint operator is represented by the conjugate transpose of the system matrix.
In our implementation we use cubic B-spline basis functions and take advantage of the closed-form expression (28). In the next section, we present an effective algorithm for calculating the forward model and its adjoint.

Fast implementation
The calculation of the Radon transform (or its variants) involves the application of the system matrix on the B-spline coefficients of the object. However the system matrix H cannot be stored explicitly due to its size. To circumvent this problem we exploit the translation-invariance of the Radon transform. This property implies that all the matrix entries in (32) can be derived from a single derivative of the Radon transform, namely that of the generating function ϕ: To improve the speed, we oversample R (n) ϕ(y, θ ) on a fine grid Y × Θ with 100 samples along each angular direction and store the values in a lookup table L. To compute the matrix entries we define a mapping with K is the number of samples along each direction, P is the number of projections and (Y ( j), Θ(i)) is the sample in Y × Θ that is nearest to (y, θ ). Therefore, we have Note that the algorithm can easily be parallelized, since projections corresponding to different angles are completely independent of each other. We designed a multithreaded implementation for an 8-core workstation which allows for 8 simultaneous projection computations. Similarly, for the adjoint of the forward model, the computation can be parallelized with respect to each object point.
In summary, our implementation is based on an accurate continuous-to-discrete model. Moreover it is fast thanks to the use of look-up tables and multi-threading. Note that our method could also be adapted to a fan-beam geometry, which can always be mapped back to the parallel one. This would lead to a non-uniform sampling pattern but our method can account for this at no additional cost (thanks to our look-up-table-based implementation).

B-spline-based discretization scheme vs Kaiser Bessel window functions
To compare the performance of the two families of basis functions, we test their consistency with a continuously-defined sample, both with respect to the forward model and with respect to the reconstruction problem. It is difficult to do this with real data for two reasons. First, we do not have access to the ground truth. Second, owing to the limited number of viewing angles and the presence of noise, one needs to apply some regularization which has a significant effect on the quality of the reconstruction when the problem is ill-posed (see Section 5). To factor out this effect, we introduce a synthetic sample for which the derivative of the Radon transform can be computed analytically.

Proposition 3. The derivative of the Radon transform of the isotropic function
where Proof. Since the function is isotropic, its Radon transform is independent of the direction. For simplicity we consider θ = 0 whose projection coordinate is parallel to the x 1 axis. Thus we have which yields the desired result.  The linearity and pseudo shift invariance of the derivative of the Radon transform then allows for an analytic computation of the projection of any function of the form where α ∈ R and x k ∈ R 2 . Figure 2(a) shows the analytic phantom that we consider as the object of interest. It was generated using 10 random radiuses a k , 10 random locations x k and 10 random intensities α k . Its differentiated sinogram is displayed in Figure 2(b) with 1800 viewing angles that are chosen uniformly between 0 and π. In our first experiment we computed the derivative of the Radon transform of the object using our discretization scheme and the one proposed in [18], which is based on Kaiser-Bessel functions. The signal-to-noise ratios (SNR) of both methods were computed with respect to the analytical projection of the phantom. The results are shown in Table 1 and demonstrate that the projections obtained with B-spline functions are more accurate than those obtained with the Kaiser-Bessel functions used in [18], and this despite the fact that the initial object is perfectly isotropic and the B-splines are not.
In our second experiment, we used the analytical projections along 1800 viewing angles as the measurement and we reconstructed the object by minimizing the quadratic cost function where H is the discretized forward operator, c contains the expansion coefficient in the chosen basis and g is the measurement vector. We use the sampled version of the object on a 1024 × 1024 grid as the ground truth. Since the data is correctly sampled and noise-free, there is not need for regularization and the reconstruction was performed using the standard CG algorithm. We used the same basis functions as in our first experiment. The results are shown in the last row of Table 1 and again demonstrate the advantage of our spline-based formulation. They confirm that the partition-of-unity condition provides a significant advantage in terms of reconstruction quality. This condition is necessary and sufficient for the approximation error to vanish as the grid size tends to zero [19]. Moreover, if we constrain the support of the basis function, polynomial B-splines are known to provide the largest-possible order of approximation, and hence the minimal approximation error [26].

Image reconstruction
We have already presented a direct method for inverting (31), namely the generalized filtered back-projection (GFBP) described in Algorithm 1. However, for large images with a limited number of measurements, direct methods such as FBP are not accurate enough. Thus, in this section, we formulate the reconstruction as a penalized least-squares optimization problem. The cost function we want to minimize is: where g is the measurement vector and the regularization term Ψ(c) is chosen based on the following considerations. First we observe that the null space of the derivative of the Radon transform contain the zero frequency (constant). To compensate for this we incorporate the energy of the coefficients into the regularization term. Second, to enhance the edges in the reconstructed image we impose a sparsity constraint in the gradient domain. This leads to where λ 1 and λ 2 are regularization parameters and {Lc} k k k ∈ R 2 denotes the gradient of the image at position k k k. More precisely, the gradient is computed based on our continuous-domain image model using the following property.
The gradient of f on the Cartesian grid is where High-dimensional sparsity-regularized problems are typically solved using the family of iterative-shrinkage/thresholding algorithms (ISTA). The convergence speed of these methods depend on the conditioning of H T H. In our case, the derivative of the Radon transform is a very ill-conditioned operator leading to very slow convergence. To overcome this difficulty, we use a variable-splitting scheme so as to map our general optimization problem into simpler ones. Specifically, we introduce the auxiliary variable u [14,[27][28][29][30] and reformulate the TV problem (42) as a linear equality-constrained problem c = argmin c,u We can map this into an unconstrained problem by considering the augmented-Lagrangian functional where α α α is a vector of Lagrange multipliers. To minimize (45), we adopt a cyclic update scheme also known as Alternating-Direction Method of Multipliers (ADMM), which consists of the iterations L μ (c, u k , α α α k ) is a quadratic function with respect to c whose gradient is We minimize L μ (c, u k , α α α k ) iteratively using the conjugate gradient (CG) algorithm to solve the equation Ac = b. Since A has a large condition number, it is helpful to introduce a preconditioning matrix M. This matrix is chosen such that M −1 H T H + μL T L + λ 1 I has a condition number close to 1. To design this problem-specific preconditioner, we use the following property of the forward model.

Proposition 5.
The successive application of the derivatives of the Radon transform and their adjoint is a highpass filter with frequency response ω ω ω 2n−1 : Proof. This follows from ( ∂ ∂ y ) * = − ∂ ∂ y and Proposition 1.
This shows that R (1) * R (1) has a Fourier transform that is proportional to ω ω ω while L T L is a discretized Laplace operator whose continuous-domain frequency response is ω ω ω 2 . Thus, the preconditioner M −1 that we use in the discrete domain is the discrete filter that approximates the frequency response 1 ω ω ω +μ ω ω ω 2 +λ 1 . A standard result for FISTA-type algorithms is that the solution of the minimization of L μ (c k , u, α α α k ) with respect to u is given by the shrinkage function [31,32] The complete reconstruction method is summarized in Algorithm 2 below. Note that in standard implementations of the CG algorithm the initial estimate is often chosen to be zero. Here we use a warm initialization in the sense that the starting point of each inner PCG iteration is the outcome of the previous PCG iteration (see step 4).

Performance metrics
We use the structural similarity measure (SSIM) [33,34] and signal-to-noise ratio (SNR) for measuring the quality of the reconstructed image. SSIM is a similarity measure proposed by Z. Wang, et al which compares the luminance, contrast and structure of images. SSIM is computed for a window of size R × R around each image pixel. The SSIM measure for two images x and x for the specified window is where C 1 and C 2 are small constants values to avoid instability. μ x and μx denote the empirical mean of image x and x in the specified window, respectively. σ x and σx are the empirical variance of the corresponding images and the covariance of two images is denoted by σ xx for the corresponding window. In our experiments, we choose C 1 = C 2 = (.001 * L) 2 where L is the dynamic range of the image pixel values. SSIM for the total image is obtained as the average of SSIM over all pixels. It takes values between 0 and 1 with 1 corresponding to the highest similarity.
Our other quality measure is the SNR. If x is the oracle and x is the reconstructed image we have Higher values of the SNR correspond to a better match between the oracle and the reconstructed image.

Parameter selection
The proposed method has several parameters which are tuned as follow. From a statistical point of view, the Baysian estimator for the additive gaussian white noise inverse problem, g = Hc + n, can be written aŝ where σ 2 n is the noise variance and p(c) is the prior density of the signal. Assuming that the gradient sample values u k k k = {Lc} k k k are i.i.d. Laplace distributed, we have where σ 2 u is the variance of the gradient values. On the other hand, by definition SNR = 10 log( Hc 2 2 /Nσ 2 n ) where N is the number of pixels in the image. Therefore σ 2 n = 10 −.1SNR Hc 2 2 /N. In practice Hc 2 2 can be approximated by g 2 2 and σ 2 u is proportional to ∑ k k k u k k k 2 2 which is empirically estimated using g 2 2 . This leads to the following rule of thumb for setting the TV regularization parameter: λ 2 ∝ 10 −.1SNR g 2 . Therefore, the TV parameter is proportional to the norm of the measurement. The proportionality constant is related to its signal-to-noise ratio.
The Thikhonov parameter is chosen as small as possible and is set to λ 1 = 10 −5 . Here we use λ 2 = ||g|| 2 × 10 −3 . Based on our experience, the parameter μ can be chosen ten times larger than λ 2 . The number of inner iterations for solving the linear step plays no role in the convergence of the proposed technique, but affects speed; we suggest to choose it as small as possible, e.g., 2 or 3. We use cubic B-splines with m = 3 as the basis functions.

Experimental result
To validate our reconstruction method, we conducted experiments with real data acquired using the TOMCAT beam line of the Swiss Light Source at the Paul Scherrer Institute in Villigen, Switzerland. The synchrotron light is delivered by a 2.9 T super-bending magnet. The energy of the X-ray beam is 25keV [35]. We used nine phase steps over two periods to measure the displacement of the diffraction pattern described in Section 2. For each step a complete tomogram was acquired around 180 degrees; we used 721 uniformly distributed projection angles. Image acquisition was performed with a CCD camera whose pixel size was 7.4μm.
For our experiments we used a rat-brain sample. The sample is embedded in liquid paraffin at room temperature. This is necessary to match the refractive index of the sample with its environment, so that the small-refraction-angle approximation holds. Finally the projections were post-processed, including flat-field and dark-field corrections, for the extraction of the phase gradient. Figure 3 gives insights into the convergence behavior of our algorithm. In comparison with the established FISTA algorithm one can observe a drastic acceleration. This results from the combination of 1) the ADMM solver, 2) our problem-specific preconditioner and 3) our warminitialization method for the inner PCG iterations. One can also observe that the convergence speed is not very sensitive to the number of inner iterations when we use warm initialization.
In practice our reconstruction scheme appears to converge after 5 outer iterations, each involving 2 inner PCG iterations. A PCG iteration essentially requires one application of the forward operator and one application of its adjoint. Since the computational complexities of the two operators are comparable, the overall cost is equivalent to 2 × 2 × 5 = 20 evaluations of the forward operator. Based on these considerations, we expect our algorithm to be about 6 times faster than the method described in [18]. This figure is based on a personal communication of the main authors of [18], who state that their algorithm converges in about 65 iterations, each being equivalent in cost to one of our inner iterations.
We further investigated the performance of the direct filtered back-projection and of the proposed iterative reconstruction techniques on a coronal rat-brain section. The reconstructed images for 721 and 181 viewing angles using GFBP and our method are shown in Figure 4. The figure of merits (SNR and SSIM) are indicated below each image. As the number of views is being reduced, the influence of regularization becomes predominant compared to the choice of a specific basis function. This was also confirmed to us by the fact that the reconstructions obtained by running the same algorithm with a different set of basis functions (Kaiser-Bessel with α = 10.4) were essentially indistinguishable for those in Figure 4(b, e). On the other hand, the effect of choosing an optimized set of basis functions is more significant in the well-posed scenario without regularization as demonstrated in Section 4.4.
The GFBP reconstruction contains artifacts at the boundary of the image and with respect to specific anatomical features. For example, the bottom-right sub images in the reconstructions of Figure 4 show the mammal-thalamic tract in this coronal section. One can clearly see oscillatory artifacts in the GFBP reconstruction. This could confuse the biologist or automated diagnostic systems for determining the nucleus and immoneurins in that region. The middle-right and top-right images show a part of the thalamus and the region between the thalamus and the hippocampus, respectively.
To reduce the aforementioned artifacts, we also implemented a smoothed version of the GFBP algorithm. Specifically, we modified the filter in the third step of Algorithm 1 as follows: q(ω y ) = 1 2π where h(ω y ) is a lowpass filter and k is an exponent that acts as a smoothing parameter. We chose h(ω) to be the standard Hamming window. The reconstructions are shown in the bottom row of Figure 4. They suggest that with this smoothed version of GFBP there is a trade-off between artifacts and image contrast. Note that for these experiments the parameter k was optimized so as to achieve the best SNR. Visually the reconstructed image with 721 angles using our proposed technique is more faithful in comparison with the GFBP approach. Therefore, we consider it as our gold standard for investigating the dependence of our algorithm on the number of views as shown in Figure 5. The SNR and SSIM values are computed for the main region of the sample which includes the brain. They suggest that we can reduce the number of views with our method by at least a factor of four while essentially maintaining the quality of the standard reconstruction method (FBP with a complete set of views).

Conclusion
Up to now, in-vivo tomography with grating interferometry faces the challenge of large dose deposition, which potentially harms the specimens e.g., in small rodent scanners. To reduce the total scanning time, we presented an exact B-spline formulation for the discretization of the derivative variants of the Radon transform that avoids any numerical differentiation. We formulated the reconstruction problem as an inverse problem with a combination of regularization terms. We introduced a new iterative algorithm that uses the alternating direction method of multipliers to solve our TV-regularized reconstruction problem. An important practical twist is the introduction of a problem-specific preconditioner which significantly speeds up the quadratic optimization step of the algorithm. Finally, we presented experimental results to validate the proposed discretization method and corresponding iterative technique. Our finding confirms that ADMM is quite competitive for solving TV-regularized problems. Moreover, our method allows for a substantial dose reduction while preserving the image quality of FBP-type methods. This is a crucial step towards the diffusion of DPCI in medicine and biology.