## Abstract

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.

© 2013 Optical Society of America

## 1. 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–3], interferometric [4–6] and free space propagation methods [7–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.18*e*/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.

#### 1.1. 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].

#### 1.2. 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–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].

#### 1.3. Outline

In Section 2, we briefly review the physics of differential phase-contrast imaging. The end-result 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.

## 2. 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

**= (**

*x**x*

_{1},

*x*

_{2}).

The phase shift induced by the object on the transmitted beam is given by

*f*(

**) =**

*x**n*(

**) − 1 with**

*x**n*(

**) is the real part of the object’s refractive index,**

*x**θ*is the angle of view of the monochromatic wave illuminating the object and

*ℛ*{

*f*(

**)}(**

*x**y*,

*θ*) is the Radon transform of

*f*(

**).**

*x***Definition 1.**[23] The (2-dimensional) Radon transform *ℛ : L*_{2}(ℝ^{2}) → *L*_{2}(ℝ× [0, *π*]) maps a function on ℝ^{2} into the set of its integrals over the lines of ℝ^{2}. Specifically, if ** θ** = (cos

*θ*, sin

*θ*) and

*y*∈ ℝ, then

*θ*^{⊥}is perpendicular to

**.**

*θ*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*:

*α*(

*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

*f*is the quantity of interest and

*g*is the physically measurable signal.

## 3. 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 *n*th derivative of the Radon transform of a function *f*(** 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:$${\mathcal{R}}^{\left(n\right)}\left\{\left(f*g\right)\left(\mathit{x}\right)\right\}\left(y,\theta \right)=\left({\mathcal{R}}^{\left(n\right)}f\left(\cdot ,\theta \right)*\mathcal{R}g\left(\cdot ,\theta \right)\right)\left(y,\theta \right)=\left(\mathcal{R}f\left(\cdot ,\theta \right)*{\mathcal{R}}^{\left(n\right)}g\left(\cdot ,\theta \right)\right)\left(y,\theta \right).$$

To derive the necessary relations for the rest of the paper, we define a new operator, the Hilbert transform along the second coordinate.

**Definition 2.** The Hilbert transform along the *x*_{2} axis, *ℋ*_{2} : *L*_{2}(ℝ^{2}) → *L*_{2}(ℝ^{2}), is defined in the Fourier domain as

*ω*

_{1},

*ω*

_{2}) are the spatial frequency coordinates corresponding to (

*x*

_{1},

*x*

_{2}).

**Proposition 1.** The sequential application of the Radon transform, the *n*th derivative operator and the adjoint of the Radon transform on function *f*(** x**) ∈

*L*

_{2}(ℝ

^{2}) is

**‖ and ${\left(-\mathrm{\Delta}\right)}^{\frac{n}{2}}$ is**

*ω**n*times application of this operator.

*Proof.* Let *g*(*y*, *θ*) = *ℛf*(*y*, *θ*). The Fourier Slice Theorem states that:

The Fourier transform of the *n*th derivative of *g*(*y*, *θ*) with respect to *y* is (*iω*)* ^{n}ĝ*(

*ω*,

*θ*). Thus,

*ω*| = ‖

**‖ with**

*ω***= (**

*ω**ω*

_{1},

*ω*

_{2}) =

*ω*(cos

*θ*, sin

*θ*). For

*θ*∈ (0,

*π*), sgn(

*ω*) = sgn(

*ω*sin

*θ*) = sgn(

*ω*

_{2}). The space-domain equivalent is

*ℛ*

^{*}

*ℛ*= (−Δ)

^{1/2}, yields the desired results.

(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\pi {\left(-1\right)}^{n}{\mathscr{H}}_{2}^{n}{\left(-\mathrm{\Delta}\right)}^{\left(n-1\right)/2}$. The transfer function of this inverse is $\frac{1}{2\pi {i}^{n}}\text{sgn}{\left({\omega}_{2}\right)}^{n}{\Vert \mathit{\omega}\Vert}^{-\left(n-1\right)}$.An equivalent form of (11) using the fact that ${\left(\frac{\partial}{\partial y}\right)}^{*}=-\frac{\partial}{\partial y}$ is

*ℛ*

^{(n)}is the adjoint of the

*n*-th derivative of the Radon transform and the transfer function of

*q*(

*y*) is:

## 4. 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.

#### 4.1. Discrete model of the Radon transform and its derivatives

We propose a formulation that relies on the practical generalized sampling scheme [19]. Let *φ*(** x**) ∈ L

_{2}(ℝ

^{2}) denote a generating Riesz basis function for the approximation space

**V**as:

The orthogonal projection of the object *f*(** x**) on the reconstruction space

**V**is

*φ*

_{k}(

**) =**

*x**φ*(

**−**

*x***) and**

*k***c**

_{k}= 〈

*f*,

*φ̃*(· −

**)〉 in which**

*k**φ̃*(

**) is the dual of**

*x**φ*(

**) with Fourier transform $\widehat{\tilde{\phi}}\left(\mathit{\omega}\right)\widehat{\phi}\left(\mathit{\omega}\right)/{\sum}_{\mathit{k}\in {\mathbb{Z}}^{2}}{\left|\widehat{\phi}\left(\mathit{\omega}-2\pi \mathit{k}\right)\right|}^{2}$.**

*x*Using the linearity and translation invariance properties of the Radon transform, the Radon
transform of *f*(** x**) is

*ℛ*{

*φ*}(

*y*,

*θ*) =

*ℱ*

^{−1}{

*φ̂*(

*ω*

**)}(**

*θ**y*). Likewise, we obtain the derivatives of the Radon transform:

*β*is a centered univariate B-spline of degree

^{m}*m*:

*n*-fold iteration of the finite-difference operator Δ

*(*

_{h}f*y*) = (

*f*(

*y*+

*h*/2) −

*f*(

*y*−

*h*/2))/

*h*.

**Proposition 2.** The *n*-th derivative of the Radon transform of the tensor-product B-spline is given by

*Proof.*We define

*g*(

*y*,

*θ*) =

*ℛ*

^{(n)}{

*β*}(

*y*,

*θ*). By an application of the Fourier-slice theorem

To discretize our forward model, we need to calculate the derivatives of the Radon transform of *f* at the location (*y _{j}*,

*θ*) where

_{i}*y*=

_{j}*j*Δ

*y*and

*θ*=

_{i}*i*Δ

*θ*

*ℛ*

^{(n)}

*φ*

_{k}(

*y*,

_{j}*θ*) are then determined using Proposition 2.

_{i}#### 4.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 Eq. (30)

where**c**is a vector whose entries are the B-spline coefficients

*c*

_{k},

**g**is the output vector and

**H**is the system matrix defined by:

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.

#### 4.3. 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
Eq. (32) can be derived from a single derivative of
the Radon transform, namely that of the generating function *φ*:

To improve the speed, we oversample *ℛ*^{(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

*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).

#### 4.4. 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

*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

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

*α*∈ ℝ and

**x**

*∈ ℝ*

_{k}^{2}. Figure 2(a) shows the analytic phantom that we consider as the object of interest. It was generated using 10 random radiuses

*a*, 10 random locations

_{k}**x**

*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

_{k}*π*.

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

## 5. 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:

**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

*λ*

_{1}and

*λ*

_{2}are regularization parameters and {

**Lc**}

_{k}∈ ℝ

^{2}denotes the gradient of the image at position

**. More precisely, the gradient is computed based on our continuous-domain image model using the following property.**

*k***Proposition 4.** Let *f*(** x**) =∑

_{k ∈2}

*c*

_{k}

*β*(

^{n}**−**

*x***). The gradient of**

*k**f*on the Cartesian grid is

*k*

_{1},

*k*

_{2}∈ , ${h}_{i}\left[{k}_{1},{k}_{2}\right]={\beta}^{n-1}\left({k}_{i}+\frac{1}{2}\right)-{\beta}^{n-1}\left({k}_{i}-\frac{1}{2}\right)$, and

*b*[

_{i}*k*

_{1},

*k*

_{2}] =

*β*(

^{n}*k*) for

_{i}*i*= 1, 2.

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–30] and reformulate the TV problem (42) as a linear equality-constrained problem

**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**

*α**ℒ*(

_{μ}**c**,

**u**

*,*

^{k}

*α**) is a quadratic function with respect to*

^{k}**c**whose gradient is

*ℒ*(

_{μ}**c**,

**u**

*,*

^{k}

*α**) iteratively using the conjugate gradient (CG) algorithm to solve the equation*

^{k}**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 ‖** ω** ‖

^{2}

^{n}^{−1}:

*Proof.*This follows from ${\left(\frac{\partial}{\partial y}\right)}^{*}=-\frac{\partial}{\partial y}$ and Proposition 1.

This shows that *ℛ*^{(1)*}*ℛ*^{(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 $\frac{1}{\Vert \mathit{\omega}\Vert +\mu {\Vert \mathit{\omega}\Vert}^{2}+{\lambda}_{1}}$.

A standard result for FISTA-type algorithms is that the solution of the minimization of *ℒ _{μ}*(

**c**

*,*

^{k}**u**,

*α**) with respect to*

^{k}**u**is given by the shrinkage function [31, 32]

## 6. Experimental analysis

#### 6.1. 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

*C*

_{1}and

*C*

_{2}are small constants values to avoid instability.

*μ*and

_{x}*μ*

_{x}_{̂}denote the empirical mean of image

**x**and

**x̂**in the specified window, respectively.

*σ*and

_{x}*σ*

_{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

#### 6.2. 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 as

*p*(

**c**) is the prior density of the signal. Assuming that the gradient sample values

**u**

_{k}= {

**Lc**}

_{k}are i.i.d. Laplace distributed, we have

*N*is the number of pixels in the image. Therefore ${\sigma}_{n}^{2}={10}^{-.1\text{SNR}}{\Vert \mathbf{Hc}\Vert}_{2}^{2}/N$. In practice ${\Vert \mathbf{Hc}\Vert}_{2}^{2}$ can be approximated by ${\Vert \mathbf{g}\Vert}_{2}^{2}$ and ${\sigma}_{u}^{2}$ is proportional to ${\sum}_{\mathit{k}}{\Vert {\mathbf{u}}_{\mathit{k}}\Vert}_{2}^{2}$ which is empirically estimated using ${\Vert \mathbf{g}\Vert}_{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.

#### 6.3. 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 warm-initialization 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:

*h*(

*ω*) is a lowpass filter and

_{y}*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).

## 7. 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.

## Acknowledgment

This work was funded (in part) by ERC Grant ERC-2010-AdG 267439-FUN-SP, the Swiss National Science Foundation under Grant 200020-144355 and the Center for Biomedical Imaging of the Geneva-Lausanne Universities and EPFL.

## References and links

**1. **V. Ingal and E. Beliaevskaya, “X-ray plane-wave tomography observation of the phase contrast from a non-crystalline object,” J. Phys. D: Appl. Phys **28**, 2314–2317 (1995). [CrossRef]

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

**3. **D. Chapman, S. Patel, and D. Fuhrman, “Diffraction enhanced X-ray imaging,” Phys., Med. and Bio. **42**, 2015–2025 (1997). [CrossRef]

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

**5. **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]

**6. **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 **13**, 6296–6304 (2005). [CrossRef] [PubMed]

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

**8. **K. A. Nugent, T. E. Gureyev, D. F. Cookson, D. Paganin, and Z. Barnea, “Quantitative phase imaging using hard X-rays,” Phys. Rev. Lett. **77**, 2961–2964 (1996). [CrossRef] [PubMed]

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

**10. **A. Momose, S. Kawamoto, I. Koyama, Y. Hamaishi, K. Takai, and Y. Suzuki, “Demonstration of X-ray talbot interferometry,” Jap. Jour. of Appl. Phys. **42**, L866–L868 (2003). [CrossRef]

**11. **F. Pfieffer, O. Bunk, C. Kottler, and C. David, “Tomographic reconstruction of three-dimensional objects from hard X-ray differential phase contrast projection images,” Nucl. Inst. and Meth. in Phys. Res. **580.2**, 925–928 (2007). [CrossRef]

**12. **M. Stampanoni, Z. Wang, T. Thüring, C. David, E. Roessl, M. Trippel, R. Kubik-Huch, G. Singer, M. Hohl, and N. Hauser, “The first analysis and clinical evaluation of native breast tissue using differential phase-contrast mammography,” Inves. radio. **46**, 801–806 (2011). [CrossRef]

**13. **M. Nilchian and M. Unser, “Differential phase-contrast X-ray computed tomography: From model discretization to image reconstruction,” Proc. of the Ninth IEEE Inter. Symp. on Biomed. Imag.: From Nano to Macro (ISBI’12), 90–93 (2012).

**14. **M. V. Afonso, J. M. Bioucas-Dias, and M. A. T. Figueiredo, “An augmented Lagrangian approach to the constrained optimization formulation of imaging inverse problems,” IEEE Trans. Imag. Proc. **20**, 681–695 (2011). [CrossRef]

**15. **S. Ramani and J. A. Fessler, “A splitting-based iterative algorithm for accelerated statistical X-ray CT reconstruction,” IEEE Trans. Med. Imag. **31.3**, 677–688 (2012). [CrossRef]

**16. **Z. Qi, J. Zambelli, N. Bevins, and G. Chen, “A novel method to reduce data acquisition time in differential phase contrast computed tomography using compressed sensing,” Proc. of SPIE **7258**, 4A1–8 (2009).

**17. **T. Köhler, B. Brendel, and E. Roessl, “Iterative reconstruction for differential phase contrast imaging using spherically symmetric basis functions,” Med. phys. **38**, 4542–4545 (2011). [CrossRef]

**18. **Q. Xu, E. Y. Sidky, X. Pan, M. Stampanoni, P. Modregger, and M. A. Anastasio, “Investigation of discrete imaging models and iterative image reconstruction in differential X-ray phase-contrast tomography,” Opt. Express **20**, 10724–10749 (2012). [CrossRef] [PubMed]

**19. **M. Unser, “Sampling–50 years after Shannon,” Proc. IEEE **88**, 254104–1–3 (2000). [CrossRef]

**20. **E. Y. Sidky and X. Pan, “Image reconstruction in circular cone-beam computed
tomography by constrained, total-variation minimization,” Phys. Med.
Biol. **53**, 4777–4807 (2008). [CrossRef]

**21. **A. Momose, W. Yashiro, Y. Takeda, Y. Suzuki, and T. Hattori, “Phase tomography by X-ray talbot interferometry for
biological imaging,” Jpn. J. Appl. Phys. **45**, 5254–5262 (2006). [CrossRef]

**22. **F. Pfeiffer, C. Grünzweig, O. Bunk, G. Frei, E. Lehmann, and C. David, “Neutron phase imaging and tomography,” Phys. Rev. Lett. **96**, 215505-1–4 (2006). [CrossRef]

**23. **F. Natterer, *The Mathematics of Computed Tomography* (John Wiley and sons, 1986).

**24. **H. Meijering, J. Niessen, and A. Viergever, “Quantitative evaluation of convolution-based methods for medical image interpolation,” Med. Imag. Anal. **5**, 111–126 (2001). [CrossRef]

**25. **A. Entezari, M. Nilchian, and M. Unser, “A box spline calculus for the discretization of computed tomography reconstruction problems,” IEEE Trans. Med. Imag. **31**, 1532 –1541 (2012). [CrossRef]

**26. **P. Thvenaz, T. Blu, and M. Unser, “Interpolation revisited [medical images application],” IEEE Trans. Med. Imag. **19.7**, 739–758 (2000). [CrossRef]

**27. **Y. Wang, J. Yang, W. Yin, and Y. Zhang, “A new alternating minimization algorithm for total variation image reconstruction,” SIAM Jour. on Imag. Sci. **1**, 248–272 (2008). [CrossRef]

**28. **T. Goldstein and S. Osher, “The split bregman method for l1-regularized problems,” SIAM Jour. on Imag. Sci. **2**, 323–343 (2009). [CrossRef]

**29. **M. Ng, P. Weiss, and X. Yuan, “Solving constrained total-variation image restoration and reconstruction problems via alternating direction methods,” SIAM Jour. on Sci. Comp. **32**, 2710–2736 (2010). [CrossRef]

**30. **B. Vandeghinste, B. Goossens, J. De Beenhouwer, A. Pizurica, W. Philips, S. Vandenberghe, and S. Staelens, “Split-bregman-based sparse-view CT reconstruction,” in “Fully 3D 2011 proc. ,” 431–434 (2011).

**31. **I. Daubechies, M. Defrise, and C. De Mol, “An iterative thresholding algorithm for linear inverse problems with a sparsity constraint,” Comm. Pure Appl. Math. **57**, 1413–1457 (2004). [CrossRef]

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

**33. **Z. Wang and A. Bovik, “A universal image quality index,” IEEE Sig. Proc. Lett. **9**, 81 –84 (2002). [CrossRef]

**34. **Z. Wang, A. Bovik, H. Sheikh, and E. Simoncelli, “Image quality assessment: From error visibility to structural similarity,” IEEE Trans. Imag. Proc. **13**, 600–612 (2004). [CrossRef]

**35. **S. McDonald, F. Marone, C. Hintermuller, G. Mikuljan, C. David, F. Pfeiffer, and M. Stampanoni, “Advanced phase-contrast imaging using a grating interferometer,” Sync. Rad. **16**, 562–572 (2009). [CrossRef]