# A Total Fractional-Order Variation Regularized Reconstruction Method for CT

Authors: Yun Chen; Baoyu Guo; Yao Lu
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

The 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. This 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 resulting model, 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.


---

## Body

## 1. 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:(1)Ax=b,where the measured projection data b∈ℝm is related to an unknown image x∈ℝn through the system matrix A∈ℝm×n. The 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], lp0<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–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.The 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. The 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.The paper is organized as follows: we recall in Section2 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.


## 2. Fractional Derivative

We describe in this section the definition of fractional derivative for a function, which is a generalization of integer-order derivative.We now recall definitions of fractional derivative with respect to a function. There are three popular definitions for fractional derivative [12–14], including the Riemann–Liouville (RL), the Grünwald–Letnikov, and the Caputo. In particular, let f denote a function on a,b⊂ℝ, and α be a fraction satisfying 0≤s−1<α<s for a positive integer s. The left-sided RL derivative of f is defined for a<t<b by(2)Datαft≔1Γn−αddts∫atfτdτt−τα−s+1,and the right-sided RL derivative of f is denoted by(3)Dtbαft≔−1sΓs−αddts∫tbfτdττ−tα−s+1,where Γ⋅ is the gamma function given by(4)Γs≔∫0∞ts−1e−tdt,which has the basic property that Γs+1=sΓs. Subsequently the Riesz-RL derivative of f is given by(5)Dabαft≔12Datαft+tDbαft.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, tk≔khk=0,1,2,…,d with t0≔a and td≔b. As indicated in [12], using the backward fractional difference approximation for the derivative Datαft at the points tk,k=0,1,…,d, yields the following formula:(6)Datkαft≈∇αftkhα=1hα∑j=0k−1jαjfk−j,k=0,1,2,…,d,where fk−j≔ftk−jh. To obtain a compact approximation form of the derivative Datαft by equation (6), let(7)F≔h−α∇αftk:k=0,1,2,…,d,f≔fk:k=0,1,2,…,d,and(8)Bα≔1hαω0α000⋯0ω1αω0α00⋯0ω2αω1αω0α0⋯0⋱⋱⋱⋱⋱⋮ωd−1α⋱ω2αω1αω0α0ωdαωd−1α⋱ω2αω1αω0α,where(9)ωjα≔−1jαj,j=0,1,2,…,d.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(10)F=Bαf,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(11)F=Cαf,with the upper triangular matrix(12)Cα≔Bα⊤.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(13)F=Sαf,with the following matrix:(14)Sα≔12hαBα+Cα.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.


## 3. 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 particular, let W1αΩ≔u∈L1Ω:uW1αΩ<+∞ be a function space embedding with the norm(15)uW1α≔∫Ωuxdx+∫Ω∇αuxdx,∇αux=∂αu∂x1,∂αu∂x2,…,∂αu∂xd⊤,Ω⊂ℝd.As mentioned in [10], the total α-order variation (TVα) is described by the following proposition.Proposition 1.
Assume thatu∈W1αΩ, then(16)TVαu=∫Ω∇αuxdx.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. This 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. The 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≔⋅,⋅H1/2 for an m×m symmetric positive definite matrix H. We further choose φ as the ℓ1-norm on ℝ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(17)min12Ax−bH2+μφDαx+ψx:x∈ℝn,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ℝ2n (Γ0ℝ2n denotes the space of all proper lower semicontinuous convex function mapping from ℝ2n to ℝ∪+∞, [18]), and ψ∈Γ0ℝn, given by(18)ψx≔0,ifx∈ℝ+n,+∞,otherwise,is the indicator function of the nonnegativity constraint set denoted by ℝ+n≔x∈ℝn:x≥0.


## 4. 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 ofp×p symmetric positive definite matrices by S+p. For a function ϑ∈Γ0ℝp, its proximity operator with respect to the matrix J∈S+p, denoted by proxϑ,J [17], is a mapping from ℝp to itself, defined for u∈ℝp by(19)proxϑ,Ju≔argmin12z−uJ2+ϑz:z∈ℝp.The subdifferential of ϑ∈Γ0ℝp is a set-valued operator [18], defined by(20)∂ϑu≔z∈ℝp:∀ω∈ℝpϑω≥ϑu+z,ω−u.There exists an equivalent relationship between the subdifferential of the function ϑ and its proximity operator with respect to J [19], described by(21)Jz∈∂ϑuif  and  only ifu=proxϑ,Ju+z.Moreover, as mentioned in [18], define the conjugate of the function ϑ as the following form:(22)ϑ∗z≔supu,z−ϑu:u∈ℝp,and denote an equivalent characterization of ∂ϑ and ∂ϑ∗ by(23)z∈∂ϑuif  and  only  ifu∈∂ϑ∗z,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.Theorem 1.
Letφ∈Γ0ℝ2n, ψ∈Γ0ℝn, H∈S+m, A∈ℝm×n, Dα∈ℝ2n×n, b∈ℝm, and μ,λ>0. If x∈ℝn is a solution of model (17), then for any Q∈S+n and P∈S+2n, there exists a vector y∈ℝ2n such that(24)x=proxλψ,Qx−λQ−1A⊤HAx−Q−1Dα⊤y+λQ−1A⊤Hb,(25)y=proxλμφ∗,Py+P−1Dαx.Conversely, if there exist Q∈S+n, .P∈S+n, x∈ℝn, and y∈ℝ2n satisfying (24) and (25), then x is a solution of model (17).Proof.
Suppose thatx 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(26)0∈λA⊤HAx−λA⊤Hb+Dα⊤∂λμφDαx+∂λψx.Here, there exists y∈∂λμφDαx such that(27)−λA⊤HAx−λA⊤Hb+Dα⊤y∈∂λψx.Following the characterization (23) for y∈∂λμφDαx, we have Dαx∈∂λμφ∗y. Moreover, by the relation (21), we know that for any P∈S+2n, PP−1Dαx∈∂λμφ∗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 existQ∈S+n and P∈S+2n such that x,y∈ℝ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:(28)w=TGw+Rw,where the operator T:ℝn×ℝ2n⟶ℝn×ℝ2n is defined at a vector w≔x,y∈ℝn×ℝ2n as Tw≔proxλψ,Qx,proxλμφ∗,Py and the operator R:ℝn×ℝ2n⟶ℝn×ℝ2n is denoted at w by(29)Rw≔λQ−1A⊤Hb−λQ−1A⊤HAx,0,and denote(30)G≔In−Q−1Dα⊤P−1DαI2n,with two identity matrices In and I2n. In particular, let Ψw≔λψx+λμφ∗y and(31)E≔diagQ,P,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∈ℝn+2n,(32)Tu−TvE2≤Tu−Tv,u−vE.


## 5. 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(33)wk+1=TGwk+Rwk,k∈ℕ0,may not yield a convergence sequence because of the expansiveness of the matrix G. Here, ℕ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(34)wk+1=TG0wk+1+G1wk+Rwk,k∈ℕ0,with(35)G0≔002P−1Dα0,G1≔In−Q−1Dα⊤−P−1DαI2n.We further rewrite (34) as the following preconditioned fixed-point proximity scheme(36)xk+1=proxλψ,Qxk−Q−1Dα⊤yk−λQ−1A⊤HAxk+λQ−1A⊤Hb,yk+1=proxλμφ∗,Pyk+P−1Dα2xk+1−xk.The 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∈ℝn+2n, there exists a unique w∈ℝn+2n such that(37)w=TG0w+G1v+Rv.We define the operator ℱ:ℝn+2n⟶ℝn+2n as(38)ℱ:v⟼w,v,wsatisfies equation37.This means that the implicit iterative scheme (34) can be characterized as an explicit iteration wk+1=ℱwk,k∈ℕ0. As indicated in [20], the explicit iteration is convergent if we can prove that ℱ is continuous and the underlying inequality holds. It is easy to prove the continuity of the operator ℱ by [21].Lemma 1.
Ifℱ is the operator defined by (38), then ℱ is continuous.
We now introduce a lemma that the iterative sequence by (34) satisfies an underlying inequality. To this end, we let(39)ϕx≔12Ax−bH2and find that its gradient ∇ϕ is Lipschitz continuous with the constant η≔A⊤HA2. Furthermore, let W≔EG1 for the above matrices E in (31) and G1 in (35) and W˜≔W−diag2ληIn,0 for the above parameter λ and the constant η. The lemma is described as follows.Lemma 2.
LetU be the fixed-point set of the operator ℱ defined by (38). If for Dα and η, two matrices Q,P and the parameter λ satisfy(40)Q−2ληIn∈S+nandP−1/2DαQ−2ληIn−1/22<1,then for any w∈U, the sequence wk=xk,yk:k∈ℕ0 yielded by the iterative scheme (34) satisfies(41)wk+1−wW2≤wk−wW2−wk−wk+1W˜2.Proof.
We first prove the symmetric positive definiteness of the involved matricesW 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ληId is a positive definite matrix, let(42)M≔Q−2ληIn−1/20KQ−2ληIn−1/2P−1/2withK≔P−1/2DαQ−2ληIn−1/2.
It can be verified that forW˜,(43)MW˜M⊤=Id00I2n−KK⊤,that is, W˜ and diagIn,I2n−KK⊤ are congruent. This implies that W˜ is positive definite if and only if K2<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 (44)wk+1−wE2≤wk+1−w,G0wk+1−w+G1wk−w+Rwk−RwE,for any w=x,y∈U. By G0=G−G1 and W=EG1, it can be verified that(45)wk+1−wW2≤wk+1−w,G1wk−wE+wk+1−w,Rwk−RwE.For the last term in (45) with Rw=−λQ−1∇ϕx,0, we let L≔diag2ληIn,0 and find that(46)2wk+1−w,Rw−RwkE≤wk−wk+1,Lwk−wk+1,since the following inequality holds(47)xk+1−x,∇ϕx−∇ϕxk=xk−xk+1,∇ϕxk−∇ϕxk+1+xk−xk+1,∇ϕxk+1+x−xk,∇ϕxk+xk+1−x,∇ϕx,≤ηxk−xk+12,due to the convexity of ϕ and Litschitz continuity of its gradient. Multiplying (45) by 2 and combining it with (46) yield(48)2wk+1−wW2≤2wk+1−w,G1wk−wE+wk−wk+1,Lwk−wk+1.Moreover, it can be verified that(49)wk+1−wW2≤wk−wW2−wk−wk+1,W−Lwk−wk+1.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.
Letwk=xk,yk∈ℝn+2n:k∈ℕ0 be the sequence yielded by the iterative scheme (36) for any initial w0∈ℝn+2n and U be the fixed-point set of the operator ℱ defined by (38). If two matrices Q,P and the parameter λ satisfy (40), then the sequence wk:k∈ℕ0 converges to a fixed point of the operator ℱ and xk:k∈ℕ0 converges to a solution of model (17).Proof.
Since two matricesQ,P and the parameter λ satisfy (40), it can be verified that inequality (41) holds for the sequence wk:k∈ℕ0 by Lemma 2.
We first prove that the sequencewk converges a cluster point of the sequence. Using inequality (41), we find that(50)wkW≤wk−wW+wW≤w0−wW+wW,for any w∈U. This means that the sequence wk is bounded. As a result, there is a subsequence wkp:p∈ℕ0 that converges to a cluster point w∗∈ℝn+2n. Summing inequality (41) for k from 0 to N−1, we have that(51)wN−wW2+∑k=0N−1wk−wk+1W˜2≤w0−wW2,which implies that(52)limk⟶+∞wk−wk+1W˜=0.Hence, the sequence wk has at most one cluster point.
We further prove that the cluster point is a fixed point of the operatorℱ. Since the operator ℱ is continuous and wk+1=ℱwk, the cluster point w∗ is a fixed point of ℱ and w∗∈U.
As a result, the sequencewk converges to a fixed point of ℱ, and xk converges to a solution of model (17).


## 6. 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 (Figure1(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. The scatter and other image degradation factors are not simulated in the experiments.Figure 1
(a) Clinic CT image; (b) the simulated projection data.
(a)
(b)Figure 2
Reconstruction images: (a) projection data without noise and with noise in the top and bottom rows, respectively; (b) results by the SART with no regularization (NR); (c) results by the TV regularized reconstruction method (TV); (d) results by the total fractional-order variation regularized reconstruction method withα=1.2 (TFV2).
(a)
(b)
(c)
(d)We then apply different methods to reconstruct the images (Figure2) from projection data with different noise levels. The 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α≔I512⊗BαBα⊗I512, where Bα is the matrix with d=511 given in (8) and I512⊗Bα is the Kronecker product of matrices I512 and Bα. As indicated in [23], for the implementation of the scheme we may let P≔I2n, Q≔β∗diagA+,1,A+,2,…,A+,n, and H≔diag1/A1,+,1/A+,2,…,1/Am,+, where Ai,+≔∑j=1nAij,i=1,2,…,m,A+,j≔∑i=1mAij,j=1,2,…,n. Moreover, as mentioned in [15, 17], we have that for x∈ℝn,(53)proxλψ,Qxj=maxxj,0,j∈ℕn,and for y∈ℝ2n,(54)proxλμφ∗,Pyi=yi−maxyi−λμ,0signyi,i∈ℕ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(55)xk+1−xk2xk+12<10−4,for the iterative sequence xk:k∈ℕ0.The 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 (Figure3). NMSE is defined by(56)NMSE≔100×∑i=1512∑j=1512Xi,j−Zi,j2∑i=1512∑j=1512Xi,j2,where X is the reference image with size 512×512 and Z is the reconstructed image. The smaller the NMSE, the better the reconstructed image. PSNR is denoted by(57)PSNRX,Z≔10log10maxX2MSEX,Z,with mean squared error(58)MSE≔15122∑i=1512∑j=1512Xi,j−Zi,j2,where maxX is the maximum value of X. The bigger the values for PSNR, the better the quality of the reconstructed image.Figure 3
Quantitative analysis for the reconstructed images from projection data with noise.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. The above metrics, NMSE and PNSR, are applied to evaluate the reconstructed images by these reconstruction methods from projection data with noise (Figure 4).Figure 4
Quantitative analyses for the reconstructed images by the TV and the TFVk from projection data with noise when we choose different regularization parametersμ1=0.15,μ2=0.20,μ3=0.35, where TFVk is the developed method with α=1.k k=1,2,…,6. (a) NMSE. (b) PSNR.
(a)
(b)The above experiment results show that the developed method is better than the competing methods in suppressing noise for CT reconstruction. Specially, quantitative analyses in Figure4 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. This means that using the TFV regularization can overcome the staircase effects in the reconstructed images, compared with the TV regularization.


## 7. Conclusions

We developed in this paper the total fractional-order variation regularized reconstruction method for CT. The method includes an optimization model with the fractional derivative-based regularizer and a preconditioned fixed-point 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.



---

*Source: 3714245-2020-02-22.xml*