Constrained regularized reconstruction of X-ray-DPCI tomograms with weighted-norm

: In this paper we introduce a new reconstruction algorithm for X-ray differential phase-contrast Imaging (DPCI). Our approach is based on 1) a variational formulation with a weighted data term and 2) a variable-splitting scheme that allows for fast convergence while reducing reconstruction artifacts. In order to improve the quality of the reconstruction we take advantage of higher-order total-variation regularization. In addition, the prior information on the support and positivity of the refractive index is considered, which yields signiﬁcant improvement. We test our method in two reconstruction experiments involving real data; our results demonstrate its potential for in-vivo and medical imaging.


Introduction
X-ray differential phase-contrast imaging (DPCI) is a tomographic technique that was first proposed by David et al. [1] and Momose et al. [2]. Among its advantages are its compatibility with regular laboratory X-ray sources and its high sensitivity.
The data provided by DPCI corresponds to the first derivative of the Radon transform of the refractive index of a sample. Thus, in practical applications, the common reconstruction scheme for DPCI is based on a variant of the filtered back-projection (FBP) algorithm. While FBP is a fast (non-iterative) method, it typically requires a large number of projections to achieve a good reconstruction quality [3]. This implies long exposure times which could damage the specimen.
Recently, several authors have proposed iterative techniques that exploit prior knowledge on the specimen to significantly reduce the number of required projections [4][5][6][7]. Their approaches are all based on a penalized maximum-likelihood formulation, with a standard 2 -norm datafidelity term. In this paper, we aim at further reducing the number of projections by proposing an improved iterative reconstruction algorithm for DPCI. The contributions of this paper are summarized as follows: 1. Formulation of the reconstruction as a variational problem using a weighted norm for the data term that ensures that the iteration matrix is well-conditioned and which speeds-up the algorithm considerably.
2. Use of a non-quadratic 1 regularization that consists either of total variation (TV) for piecewise-constant images or a higher-order extension that is better matched to biological specimens.
3. Inclusion of support and positivity constraints in the reconstruction algorithm.
4. Design of a novel variable-splitting scheme for the constrained optimization problem, which improves the reconstruction quality compared to [4]. As a result, the last step of every iteration is a denoising operation, which is beneficial for suppressing artifacts. Furthermore, this splitting is inherently matched to our preconditioner, leading to a fastconverging numerical scheme.
Finally, we conduct experiments on real data to evaluate the proposed method. A preliminary version of this paper has been presented at ISBI 2013 [8]. The overlap is only the use of a weighted norm (item 1), but the presented algorithm and the other refinements (items 2-4) are specific to this paper. The remainder of the paper is organized as follows: In Section 2, we review the continuousdomain formulation of the DPCI forward model and state an important relationship between the reconstruction and the measurements. In Section 3, we describe the discretization scheme of the forward model. We propose there a novel formulation of the reconstruction. In addition, the details of the algorithm are discussed. We evaluate the proposed reconstruction scheme in two experiments involving real data in Section 4. We summarize and conclude the paper in Section 5.

Continuous-domain forward model
Let f denote the 2D distribution of the refractive-index of the object. In DPCI, the measurement g is the first derivative of the Radon transform (FDRT) of this function, with where R{ f }(y, θ ) = R f (yu u u θ +tv v v θ ) dt and (u u u θ , v v v θ ) are orthonormal vectors that form an angle θ with some reference coordinate system. An important property of the continuous-domain FDRT is that it can be inverted formally using a 1D convolution operator followed by the adjoint of FDRT. This leads to Here * y denotes a 1D convolution with respect to the variable y. The frequency response of the convolution kernel h is given by h(ω) = 1/|ω|.

Model discretization
In order to discretize the forward model (1), we expand f as where β n is a tensor-product polynomial B-spline of degree n. In practice, the domain of f is bounded. Therefore, we are only interested in a finite number of expansion coefficients c k k k , which we gather in a vector c. The measurement g is only known at discrete locations and we collect the corresponding values in a vector g. Then, in the absence of noise, (1) implies the linear relationship g = Hc, where H is a matrix that represents the discretized FDRT in the B-spline basis [4].

Image reconstruction
We formulate the reconstruction as a constrained optimization problem with a generalized weighted 2 -norm data term. Specifically, we aim at finding the vector c 0 such that where g is an (M × 1) data vector, H is an (M × N) forward projection matrix, · 2 W is the weighted norm which is defined as W·, · , and C is a convex set that enforces support and positivity constraints. Since reconstruction with a limited number of projections is an ill-posed problem, we introduce regularization terms to take advantage of the prior information. The regularizations Ψ 1 (c) and Ψ 2 (c) are smooth and non-smooth, respectively. The parameters λ 1 , λ 2 ∈ R control the strength of the regularization.
For rapid convergence, our proposal is to use a weighting matrix which is the discrete counterpart of the convolution operator h in (2), with a slight modification to the frequency-domain singularity at zero. We modify it with the frequency response 1 |ω|+β , where β is an appropriate positive parameter; it is a positive-definite operator.
We solve the nonlinear regularized problem while defining an axillary variable u and using an augmented-Lagrangian (AL) scheme.
This in turn is equivalent to finding critical point of the augmented Lagrangian (AL) where α α α is a vector of Lagrange multipliers that imposes the constraint u = c. The classical AL scheme alternates between a joint minimization step and an update step, so that Moreover, we use alternating direction method of multipliers (ADMM) [9] to separate the joint minimization into the succession of simpler partial problems Since the zero frequency is in the null-space of the forward operator, we use Tikhonov regularization term Ψ 1 (u) = 1/2 u 2 .
In Step 1, c k and α α α k are fixed, therefore L µ (c k , u, α α α k ) is a quadratic function of u whose gradient is We use the conjugate-gradient (CG) method to solve this step. Since the condition number of the matrix H T WH + (µ + λ 1 )I is quite small, the corresponding iterative algorithm converges rapidly.
Step 2 of ADMM, which minimizes L µ (c, u k , α α α k ) with respect to c, is the constrained denoising problem The common expression for the regularizer is where · is a non-quadratic norm and R : R N → R (NK) is the regularization operator (e.g., gradient with K = 2 or Hessian with K = 2 × 2). For the identity regularization operator, R = I, (9) typically admits a direct threshold-based solution.
For the general case of regularization operator, we aim at solving the denoising problem, which is equivalent to the proximal map prox · ,λ ,P C (z) = argmin where P C is the convex projection that corresponds to the constraint. In order to find the solution of (11), we use Fenchel duality to rewrite the regularization term as where R T : R (NK) → R N is the adjoint of the operator R, p ∈ R (NK) and B := {p ∈ R (NK) | p * ≤ 1} with · * the dual norm. It can be shown that the solution of (11) is with We apply fast iterative shrinkage-thresholding algorithm (FISTA) [10] to solve (13). The step size is constrained by the Lipschitz constant L of ∇ ∇ ∇ f (p) that depends on the regularization operator R. The other important component is the orthogonal projection onto the set B that is specified by the chosen norm. Let us denote it by P B . Algorithm 1 describes the denoising algorithm.
Input: z, λ , τ ≤ L −1 , P B , P C Output: c (optimal solution of (11)) initialization p 0 , t 1 = 1; while stopping criterion is not satisfied do Algorithm 1: DENOISING ALGORITHM The benefits of the proposed splitting are: 1) the transformation of a complex reconstruction problem into a sequence of simpler optimizations where the constraint is applied as simple projection in each iteration of the denoising step; 2) any regularization term can be handled by knowing its corresponding denoising function; 3) the output of the algorithm is the solution of the denoising step that results in an improved quality of reconstruction.
The reconstruction method is summarized in Algorithm 2. Here, the starting point of each inner CG iteration is the outcome of the previous CG iteration called as warm initialization.
As for the regularization, we consider two distinct options 1) Our first option is the use of a total-variation (TV) regularization term to enhance the edges in the reconstructed image. Therefore, we set with Lc 1,1 = ∑ i {Lc} i 1 , where the sum is computed on all B-spline coefficients and {Lc} i ∈ R 2 is the gradient vector of the image at position i. The discrete gradient operator L : R N → R N×2 is computed according to proposition 2 in [4].
Here, the regularization operator is the discrete gradient operator and the mixed 1 − 1 norm is chosen as the potential function. As the dual norm of the 1 norm is ∞ , the dual ball is defined Input: g, H, λ 1 , λ 2 , P C Output: ( f (x x x) reconstructed image) set λ 1 , λ 2 , µ, and B-spline degree m; initialization c 0 , u 0 and α α α 0 ; while stopping criterion is not satisfied do u k+1 ← argmin u L µ (c k , u, α α α k ), using CG method with initial estimate u k ("warm initialization"); Therefore, the orthogonal projection of where [·] j is the j-th entry of the corresponding vector and y = [ y T 1 , y T 2 , ..., y T N ] T . This regularization is well-matched to piecewise-constant images.
2) Owing to the fact that biological and medical specimens consist mostly of filament-like and complicated structures, we investigate higher-order extensions of total variation. We apply Hessian Schatten-norm regularization (HS) as our second option. It can eliminate the staircase effect of TV regularization and results in piecewise-smooth variations of intensity in the reconstructed image. We set Ψ 2 (c) = Hc 1,S 1 , where H : R N → R N×2×2 is the discrete Hessian operator and Hc 1,S 1 is the mixed of 1 and nuclear norm. The norm can be computed with Hc 1,S 1 = ∑ i (σ 1,i + σ 2,i ), where σ 1,i and σ 2,i are the singular values of the Hessian matrix at position i. Therefore, the corresponding unit-norm dual ball is defined as where · S ∞ is the ∞ -norm of the singular values of the corresponding matrix (for more details, we refer the reader to [11]).

Parameter setting
The proposed algorithm has several parameters. We adjust them as follows: • parameters λ 1 and λ 2 : we use the approach proposed in [4]; λ 1 = 10 −5 and λ 2 = 10 −4 g . The experimental results suggest that this choice of parameters yields the optimal performance.
• parameter µ: this parameter affects the convergence speed of ADMM. Since the algorithm is not too sensitive to it, we use a fixed value (µ = 1).
• parameter λ : this is a parameter of the proximal map operator in Eq. 11. Since the second step of ADMM is solving Eq. 9, we have λ = λ 2 /µ.
• Lipschitz constant L: the Lipschitz constant of ∇ ∇ ∇ f (p) = −λ RP C (z − λ R T p) is approximated by the Lipschitz constant of the same operator without the convex projection P C since the projection on the convex set is firmly non-expansive. Thus, where λ max (A) is the maximum eigenvalue of the matrix A. For our regularization scheme λ max (RR T ) ≤ γ where γ = 8 for the TV regularization for two-dimensional problems, and its value is 64 for the HS regularization as computed in [11].

Experimental results
We compared the proposed algorithm to FBP and to ADMM-PCG, which appears to be the current state of the art for the reconstruction of X-ray-DPCI tomograms [4]. All experiments involved real data acquired at the TOMCAT beam line of the Swiss Light Source at the Paul Scherrer Institute in Villigen, Switzerland. The common approach for these experiments is to use a reconstruction from a large number of projections as a reference for evaluating results obtained with significantly fewer projections. In addition, the convex constraints that we apply are the positivity of the refractive index combined with the support-related constraint that the solution should be zero outside the tube that contains the specimen.
In order to identify the benefits of the proposed algorithm (CRWN), we first tested the algorithms under extreme conditions: We used only 72 projections as input, while the reference was reconstructed from 1,200 projections. For this first experiment we used a phantom that was composed of a tube and three cylinders containing liquids with different refractive indices as shown in Fig. 1(a).  The performance of different algorithms are compared in Table 1. Clearly, the new method outperforms ADMM-PCG [4]. Applying the convex constraint improves the signal-to-noise ratio (SNR) and structural similarity index measure (SSIM) [12] even further. The result of the algorithm proposed in [8] is the same as CRWN-TV without CC, but it is slower since it uses FISTA. As expected, owing to the piecewise-constant structure of the sample, TV outperforms HS regularization.
We conducted another experiment with a coronal section of a scaffold that is used for surgery. The reference image was built from 2,000 projections as depicted in Fig. 1(b). The algorithms were then benchmarked on a subset of 250 projections. Although these conditions are less severe, FBP still produces high-frequency patterns that are visible in Fig. 2(a). ADMM-PCG almost completely suppresses these artifacts, at the expense of light smoothing as shown in Fig. 2(b). Overall, CRWN yields sharper images Fig. 2(c) and 2(d), which are also reflected by the quality metrics. In addition, Hessian type regularization eliminates the staircase effect of TV which is more visible in the selected region of interest.
It is seen in Fig. 3(a) that CRWN is significantly faster at minimizing the cost functional than the standard FISTA algorithm. In addition, it appears that the convergence speed is not very sensitive to the number of inner iterations as we use warm initialization. We illustrate in Fig. 3(b) the robustness of CRWN with respect to the number of projections in terms of SNR. Owing to the poor performance of FBP in reconstructing boundaries, we compute the SNR for the specified region with dashed circle shown in Fig. 1(b).

Conclusion
We have proposed a new iterative method for the reconstruction of X-ray-DPCI tomograms, whose experimental performance exceeds that of recent algorithms. We showed that side information such as the support-related constraints and positivity of the refractive index can significantly improve the quality of reconstruction. We interpreted this side information as a convex constraint in a variational formulation. In addition, we took advantage of Hessian-type regularization to reduce the staircase effect of TV and enhance the quality of higher-order structures in the sample. Our results demonstrated its potential for in-vivo and medical imaging.