A Total Fractional-Order Variation Regularized Reconstruction Method for CT

+e total variation (TV) regularized reconstruction methods for computed tomography (CT) may lead to staircase effects in the reconstructed images because of using the TV regularization. +is paper develops a total fractional-order variation regularized CT reconstruction method, aiming at overcoming the weakness of the reconstruction methods based on the TV. Specifically, we propose an optimization model for CT reconstruction, including a fidelity term, a regularization term, and a constraint term. Here, the regularization is a total fractional-order variation arising from the fractional derivative of the underlying solution. To address the nondifferentiability of the resultingmodel, we introduce a fixed-point characterization for its solution through the proximity operators of the nondifferentiable functions. Based on the characterization, we further develop a fixed-point iterative scheme to solve the resulting model and provide convergence analysis of the developed scheme. Numerical experiments are presented to demonstrate that the developed method outperforms the TV regularized reconstruction method in terms of suppressing noise for CT reconstruction.


Introduction
Computed tomographic (CT) technology provides patients' anatomical information through reconstruction of measured X-ray intensities (projection data). A main research topic for CT is to improve the quality of the reconstructed images. Mathematically, CT reconstruction requires solving an ill-posed problem [1], described by the following linear system: where the measured projection data b ∈ R m is related to an unknown image x ∈ R n through the system matrix A ∈ R m×n . e conventional total variation (TV) regularized reconstruction methods may lead to staircase effects in the reconstructed images due to the use of the TV regularization.
To overcome the weakness of these methods, this paper investigates a total fractional-order variation regularized CT reconstruction method.
Regularization is necessary to CT reconstruction problem due to its ill-posedness. Many regularization methods have been proven to be effective for CT reconstruction, for instance, the TV [2,3], total generalized variation [4,5], l p (0 < p < 1) regularization [6], and physics based priors [7]. In particular, CT reconstruction methods based on the TV regularization can effectively suppress noise and preserve edges of the reconstructed images. However, using the TV regularization may lead to the so-called staircase effects for the reconstructed image.
Fractional derivative-based regularization methods are studied for overcoming the difficulty of the TV in image processing [8][9][10][11]. In particular, the authors of literatures [12,13] systematic analyzed the discretization of fractional derivative. Zhang and Chen [10] studied a fractional derivative-based regularizer and presented the fractional derivative-based total fractional-order variation (TFV) model for image restoration. As indicated in [10], the use of the fractional derivative leads to a satisfactory compromise such as no staircasing and in preserving important fine-scale features such as edges and textures. e goal of the paper is to develop a regularized CT reconstruction method with a TFV regularization for improving the quality of the reconstructed images. Specially, we propose an optimization model for CT reconstruction, including a fidelity term, a TFV regularization term, and a nonnegativity constraint term. Here, the regularization is based on the fractional derivative of the underlying solution.
e main challenge to solve the resulting model arises from its nondifferentiability. To address this challenge, we employ the proximity operators of the two nondifferentiable terms to formulate a solution of the resulting model as a system of fixed-point equations. We further develop a preconditioned fixed-point proximity scheme from the fixed-point characterization and provide convergence analysis of the developed iterative scheme. We then conduct numerical experiments to demonstrate the effectiveness of the developed method for CT reconstruction, in comparison with the TV regularized reconstruction method. Numerical results show that the developed method is superior to the competing methods in suppressing noise for CT reconstruction. e paper is organized as follows: we recall in Section 2 the definition of fractional derivative for a function and its discretization. In Section 3, we describe a TFV regularization and develop an optimization model for CT reconstruction. We introduce a fixed-point characterization for a solution of the resulting model in Section 4. Section 5 is devoted to developing an iterative scheme by the resulting characterization and analyzing its convergence. In Section 6, numerical experiments are presented to demonstrate the effectiveness of the developed method for CT reconstruction. We draw in Section 7 a conclusion for this paper.

Fractional Derivative
We describe in this section the definition of fractional derivative for a function, which is a generalization of integerorder derivative.
We now recall definitions of fractional derivative with respect to a function. ere are three popular definitions for fractional derivative [12][13][14], including the Riemann-Liouville (RL), the Grünwald-Letnikov, and the Caputo. In particular, let f denote a function on [a, b] ⊂ R, and α be a fraction satisfying 0 ≤ s − 1 < α < s for a positive integer s. e left-sided RL derivative of f is defined for and the right-sided RL derivative of f is denoted by where Γ(·) is the gamma function given by which has the basic property that Γ(s + 1) � sΓ(s). Subsequently the Riesz-RL derivative of f is given by For the other definitions, one can refer to [10,14]. We further describe the discretization of the RL derivative [12,13]. We first yield d + 1 equidistant points with the step h from the interval [a, b] through sampling; that is, t k ≔ kh (k � 0, 1, 2, . . . , d) with t 0 ≔ a and t d ≔ b. As indicated in [12], using the backward fractional difference approximation for the derivative a D α t f(t) at the points t k , k � 0, 1, . . . , d, yields the following formula: To obtain a compact approximation form of the derivative a D α t f(t) by equation (6), let and where As a result, with equations (7)-(9), accumulating formula (6) in column wise for each k ∈ 0, 1, 2, . . . , d { } may lead to the compact approximation form of the left-sided RL derivative (2) for f as where B α is the discrete analogue of the left-sided RL derivative of order α [12]. Similarly to the above analysis, we may have the approximation of the right-sided RL derivative (3) of f, denoted by with the upper triangular matrix By definition (5), the Riesz-RL derivative of order α for f can be described as a combination of the approximations (10) and (11) [13], defined by with the following matrix: Moreover, the fractional derivative matrices suitable for an image can be generated by the Kronecker product [10,15] of the above matrices and the identity matrix.

Optimization Model
We describe the total α-order variation and further develop an optimization model for CT reconstruction problem in this section. We first recall the total α-order variation. Zhang and Chen [10] studied the total α-order variation motivated by the TV [2] and the total generalized variation [4]. In par- As mentioned in [10], the total α-order variation (TV α ) is described by the following proposition.
Similarly to the pixel-based discrete form of the TV [2,15,16], we consider the pixel-based piecewise constant approximation of the anisotropic TV α due to the anisotropic structures of CT images. In particular, the resulting approximation can be characterized as a composition structure of the ℓ 1 -norm and the fractional derivative matrix suitable for an image. Compared to the first-order difference matrix involved in the discrete anisotropic TV, the fractional derivative matrix is a lower triangular matrix. is implies that using the anisotropic TV α may preserve important fine-scale features of the reconstructed image through considering a linear combination of more neighboring image intensities.
We then develop an optimization model based on the anisotropic TV α for CT reconstruction problem. e underlying model includes a fidelity term with differentiability and two nondifferentiable terms. Specially, the fidelity term is defined by a weighted least square norm. Here, let 〈·, ·〉 H ≔ 〈·, H·〉 be the weighted inner product with the corresponding norm ‖·‖ H ≔ 〈·, ·〉 1/2 H for an m × m symmetric positive definite matrix H. We further choose φ as the ℓ 1 -norm on R 2n and compute the anisotropic TV α by the composition φ(D α ·) of the ℓ 1 -norm and the fractional derivative matrix D α suitable for an image. Note that the attenuation coefficients are nonnegativity according to CT physical properties, and this is applied as the nonnegativity constraint [17] added to the cost function. Consequently, we propose optimization model with the TV α and the nonnegativity constraint for solving (1), denoted by where D α is the 2n × n matrix arising from the discrete analogue of the left-sided RL derivative of order α, μ is a regularization parameter, φ ∈ Γ 0 (R 2n ) (Γ 0 (R 2n ) denotes the space of all proper lower semicontinuous convex function mapping from R 2n to R ∪ +∞ { }, [18]), and ψ ∈ Γ 0 (R n ), given by is the indicator function of the nonnegativity constraint set denoted by R n

Fixed-Point Characterization
We introduce in this section a fixed-point characterization for a solution of model (17) through the proximity operators of two nondifferentiable functions in the model. We recall some basic definitions and notation for describing the fixed-point characterization. Denote the set of p × p symmetric positive definite matrices by S p + . For a function ϑ ∈ Γ 0 (R p ), its proximity operator with respect to the matrix J ∈ S p + , denoted by prox ϑ,J [17], is a mapping from R p to itself, defined for u ∈ R p by e subdifferential of ϑ ∈ Γ 0 (R p ) is a set-valued operator [18], defined by ere exists an equivalent relationship between the subdifferential of the function ϑ and its proximity operator with respect to J [19], described by Moreover, as mentioned in [18], define the conjugate of the function ϑ as the following form: and denote an equivalent characterization of z ϑ and z ϑ * by for u ∈ dom(ϑ) and z ∈ dom(ϑ * ). Using a way similar to [17,19], a solution of model (17) is characterized by a system of fixed-point equations in the following theorem.
Mathematical Problems in Engineering solution of model (17), then for any Q ∈ S n + and P ∈ S 2n + , there exists a vector y ∈ R 2n such that Conversely, if there exist Q ∈ S n + , .P ∈ S n + , x ∈ R n , and y ∈ R 2n satisfying (24) and (25), then x is a solution of model (17).
Proof. Suppose that x is a solution of model (17), by Fermat's rule and the Chain rule of the subdifferential, we find that for all λ > 0, the following inclusion relation holds Here, there exists y ∈ z (λμφ)(D α x) such that Following the characterization (23) for y ∈ z (λμφ)(D α x), we have D α x ∈ z (λμφ) * (y). Moreover, by the relation (21), we know that for any P ∈ S 2n + , PP − 1 D α x ∈ z (λμφ) * (y), which yields (26). For any Q ∈ S n + , multiplying the left-hand side of (27) by QQ − 1 and using (21) lead to equation (24).
Conversely, if these exist Q ∈ S n + and P ∈ S 2n + such that (x, y) ∈ R n+2n satisfies (24) and (25), then all the arguments discussed above are reversible.

□
For the development and convergence analysis of the underlying iterative scheme, we now reformulate the coupled equations (24) and (25) as the following compact formula: where the operator T: R n × R 2n ⟶ R n × R 2n is defined at a vector w ≔ (x, y) ∈ R n × R 2n as T(w) ≔ (prox λψ,Q (x), prox (λμφ) * ,P (y)) and the operator R: R n × R 2n ⟶ R n × R 2n is denoted at w by and denote with two identity matrices I n and I 2n . In particular, let Ψ(w) ≔ λψ(x) + (λμφ) * (y) and and we find that the operator T is considered as the proximity operator of the function Ψ with respect to E, given by T � prox Ψ,E . As indicated in [19], the operator T is firmly nonexpansive with respect to E, which implies that for all u, v ∈ R n+2n ,

Iterative Scheme and Its Convergence
In this section, we develop an iterative scheme based on the above fixed-point characterization and analyze its convergence. We first develop a preconditioned fixed-point proximity scheme to solve the resulting model from the characterization. We then provide convergence analysis for the developed iterative scheme. We develop an iterative scheme based on the fixed-point (28). As indicated in [15,19], the explicit Picard iteration by the fixed-point equation may not yield a convergence sequence because of the expansiveness of the matrix G. Here, N 0 ≔ 0, 1, 2, . . . { }. To overcome the difficulty, motivated by [19], we study an implicit iteration form through employing a splitting strategy to the matrix G, denoted by with We further rewrite (34) as the following preconditioned fixed-point proximity scheme e resulting iterative scheme overcomes the difficulty arising from the nondifferentiability of model (17) and involves the preconditioning matrices for developing the efficient algorithm. Moreover, it can be implemented explicitly and guaranteed to converge under proper conditions. We next analyze convergence of the developed scheme through investigating the continuity of an underlying operator and an underlying inequality motivated by [20,21]. In particular, the inequality may lead to a unique cluster point of the iterative sequence through the Picard iteration of the underlying operator, and it is easy to prove that the cluster point is a fixed point of the operator via its continuity. Some definitions are required to analyze that. As indicated in [19], we can transform the implicit iterative scheme to an explicit iteration through defining an operator. Suppose that for any v ∈ R n+2n , there exists a unique w ∈ R n+2n such that We define the operator F: R n+2n ⟶ R n+2n as is means that the implicit iterative scheme (34) can be characterized as an explicit iteration w k+1 � Fw k , k ∈ N 0 . As indicated in [20], the explicit iteration is convergent if we can prove that F is continuous and the underlying inequality holds. It is easy to prove the continuity of the operator F by [21].

Lemma 1. If F is the operator defined by (38), then F is continuous.
We now introduce a lemma that the iterative sequence by (34) satisfies an underlying inequality. To this end, we let and find that its gradient ∇ϕ is Lipschitz continuous with the constant η ≔ ‖A ⊤ HA‖ 2 . Furthermore, let W ≔ EG 1 for the above matrices E in (31) and G 1 in (35) and W ≔ W − diag(2ληI n , 0) for the above parameter λ and the constant η. e lemma is described as follows.

Lemma 2. Let U be the fixed-point set of the operator F defined by (38). If for D α and η, two matrices Q, P and the parameter λ satisfy
then for any w ∈ U, the sequence w k � (x k , y k ): k ∈ N 0 yielded by the iterative scheme (34) satisfies Proof. We first prove the symmetric positive definiteness of the involved matrices W and W for the following analysis. For W and W, the symmetric positive definiteness of W implies that of W, and W is a symmetric matrix. We need to prove that W is positive definite. Since Q − 2ληI d is a positive definite matrix, let It can be verified that for W, that is, W and diag(I n , I 2n − KK ⊤ ) are congruent. is implies that W is positive definite if and only if ‖K‖ 2 < 1. Hence, W and W are symmetric positive definite matrices.
We then verify that (41) holds as follows: Since the operator T is firmly nonexpansive with respect to E, by (32), we have that for any w � (x, y) ∈ U. By G 0 � G − G 1 and W � EG 1 , it can be verified that For the last term in (45) with R(w) � (− λQ − 1 ∇ϕ(x), 0), we let L ≔ diag(2ληI n , 0) and find that since the following inequality holds due to the convexity of ϕ and Litschitz continuity of its gradient. Multiplying (45) by 2 and combining it with (46) yield Moreover, it can be verified that Recalling W � W − L and its symmetric positive definiteness, it yields (41). Using the above lemmas, we can prove that the iterative sequence yielded by the iterative scheme (36) is convergence.
Theorem 2. Let w k � (x k , y k ) ∈ R n+2n : k ∈ N 0 be the sequence yielded by the iterative scheme (36) for any initial w 0 ∈ R n+2n and U be the fixed-point set of the operator F defined by (38). If two matrices Q, P and the parameter λ satisfy (40), then the sequence w k : k ∈ N 0 converges to a fixed point of the operator F and x k : k ∈ N 0 converges to a solution of model (17).
Proof. Since two matrices Q, P and the parameter λ satisfy (40), it can be verified that inequality (41) holds for the sequence w k : k ∈ N 0 by Lemma 2.
We first prove that the sequence w k converges a cluster point of the sequence. Using inequality (41), we find that for any w ∈ U. is means that the sequence w k is bounded. As a result, there is a subsequence w k p : p ∈ N 0 Mathematical Problems in Engineering that converges to a cluster point w * ∈ R n+2n . Summing inequality (41) for k from 0 to N − 1, we have that which implies that Hence, the sequence w k has at most one cluster point. We further prove that the cluster point is a fixed point of the operator F. Since the operator F is continuous and w k+1 � F(w k ), the cluster point w * is a fixed point of F and w * ∈ U.
As a result, the sequence w k converges to a fixed point of F, and x k converges to a solution of model (17).

Numerical Experiments
Numerical experiments are presented to show the superiority of the developed method over other competing methods for simulating low-dose CT reconstruction. We first obtain the simulated projection data from clinic CT image. We then compare the developed method to the simultaneous algebraic reconstruction technique (SART) with the nonnegative constraint and the TV regularized reconstruction method. We further conduct an experiment to demonstrate the effectiveness of the developed method with different values of order α.
We first simulate projection data from clinic CT image. Specially, we obtain noise-free parallel-beam projection data (Figure 1(b)) with 120 angles and 729 bins from clinic CT image with size 512 × 512 (Figure 1(a)) through Radon transform in MATLAB. We further have projection data with noise ( Figure 2) by adding Gaussian noise (with variance � 25, mean � 0, and standard deviation � 1) to the simulated projection data. e scatter and other image degradation factors are not simulated in the experiments.
We then apply different methods to reconstruct the images ( Figure 2) from projection data with different noise levels. e methods include the SART with the nonnegative constraint, the TV regularized reconstruction method, and the developed method. Here, the first method is to apply the first step of the iterative scheme (36) to solve model (17) without the regularization; the second method is to employ the developed iterative scheme to solve the model yielded by replacing the regularization term φ(D α x) in model (17) with the pixel-based anisotropic TV [15]; and the last one is the TFV regularized reconstruction method which applies the developed iterative scheme to solve model (17). In particular, the system matrix A is yielded by using Siddon's algorithm [22], and the matrix D α ≔ I 512 ⊗ B α B α ⊗ I 512 , where B α is the matrix with d � 511 given in (8) and I 512 ⊗ B α is the Kronecker product of matrices I 512 and B α . As indicated in [23], for the implementation of the scheme we may let P ≔ I 2n , Q ≔ β * diag(A +,1 , A +,2 , . . . , A +,n ), and H ≔ diag(1/A 1,+ , A +,j ≔ m i�1 A ij , j � 1, 2, . . . , n. Moreover, as mentioned in [15,17], we have that for x ∈ R n , and for y ∈ R 2n , We then set λ � 0.8 and β � 1 and choose μ � 0.05 for the noise-free case and μ � 0.2 for the case with noise in the experiment. Iterations stop when for the iterative sequence x k : k ∈ N 0 . e normalized mean square error (NMSE) and the peak signal-to-noise ratio (PSNR) are used to evaluate the reconstructed images by the above reconstruction methods ( Figure 3). NMSE is defined by where X is the reference image with size 512 × 512 and Z is the reconstructed image. e smaller the NMSE, the better the reconstructed image. PSNR is denoted by PSNR(X, Z) ≔ 10 log 10 max(X) 2 MSE(X, Z) , with mean squared error where max(X) is the maximum value of X. e bigger the values for PSNR, the better the quality of the reconstructed image.
We further conduct an experiment by using different values of order α for showing the effectiveness of the developed method, as compared to the TV regularized reconstruction method. In particular, with the above parameters λ and β, we choose three different regularization parameters (i.e., μ � 0.15, 0.20, 0.35) for such two reconstruction methods and set different values of order α (i.e., α ∈ 1.1, 1.2, 1.3, 1.4, 1.5, 1.6 { }) for the developed method. e above metrics, NMSE and PNSR, are applied to evaluate the reconstructed images by these reconstruction methods from projection data with noise ( Figure 4). e above experiment results show that the developed method is better than the competing methods in suppressing noise for CT reconstruction. Specially, quantitative analyses in Figure 4 shows that the TFV and the TV yield the best PSNR and NMSE when μ � 0.2, and the former is superior to the latter. Furthermore, the TFV with α � 1.2 leads to the best PSNR for the developed method when μ � 0.2 and μ � 0.15. When the regularization parameter μ is too large, the TFV outperforms the TV in the aspects of NMSE and noise suppression. is means that using the TFV regularization can overcome the staircase effects in the reconstructed images, compared with the TV regularization.

Conclusions
We developed in this paper the total fractional-order variation regularized reconstruction method for CT. e method includes an optimization model with the fractional derivative-based regularizer and a preconditioned fixedpoint proximity scheme to solve the resulting model. Numerical results showed that the developed method performs better than the competing methods in terms of NMSE and noise suppression (PSNR). In particular, the total fractional-order variation regularization is superior to the total variation regularization in suppressing noise (PSNR) for CT reconstruction.
Data Availability e involved data in this paper are CT images.