A Total Fractional-Order Variation Regularized Reconstruction Method for CT

Chen, Yun; Guo, Baoyu; Lu, Yao

doi:https://doi.org/10.1155/2020/3714245

Mathematical Problems in Engineering

On this page

Abstract Introduction Conclusions Data Availability Conflicts of Interest Acknowledgments References Copyright Related Articles

Research Article | Open Access

Volume 2020 | Article ID 3714245 | https://doi.org/10.1155/2020/3714245

A Total Fractional-Order Variation Regularized Reconstruction Method for CT

Yun Chen,^1,2Baoyu Guo,^1,2and Yao Lu^2,3,4

Academic Editor: Ines Tejado Balsera

Received11 Nov 2019

Accepted18 Jan 2020

Published22 Feb 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.

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:where the measured projection data is related to an unknown image through the system matrix . 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], 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 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.

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 , and α be a fraction satisfying for a positive integer s. The left-sided RL derivative of f is defined for byand the right-sided RL derivative of f is denoted bywhere is the gamma function given bywhich has the basic property that . Subsequently the Riesz-RL derivative of f is given byFor the other definitions, one can refer to [10, 14].

We further describe the discretization of the RL derivative [12, 13]. We first yield equidistant points with the step h from the interval through sampling; that is, with and . As indicated in [12], using the backward fractional difference approximation for the derivative at the points , yields the following formula:where . To obtain a compact approximation form of the derivative by equation (6), letandwhereAs a result, with equations (7)–(9), accumulating formula (6) in column wise for each may lead to the compact approximation form of the left-sided RL derivative (2) for f aswhere 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 bywith the upper triangular matrixBy definition (5), the Riesz-RL derivative of order α for f can be described as a combination of the approximations (10) and (11) [13], defined bywith 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.

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 be a function space embedding with the norm

As mentioned in [10], the total α-order variation (TV^α) is described by the following proposition.

Proposition 1. Assume that , then

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 -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 be the weighted inner product with the corresponding norm for an symmetric positive definite matrix H. We further choose φ as the -norm on and compute the anisotropic TV^α by the composition of the -norm and the fractional derivative matrix 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 bywhere is the matrix arising from the discrete analogue of the left-sided RL derivative of order α, μ is a regularization parameter, ( denotes the space of all proper lower semicontinuous convex function mapping from to , [18]), and , given byis the indicator function of the nonnegativity constraint set denoted by .

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 of symmetric positive definite matrices by . For a function , its proximity operator with respect to the matrix , denoted by [17], is a mapping from to itself, defined for byThe subdifferential of is a set-valued operator [18], defined byThere exists an equivalent relationship between the subdifferential of the function ϑ and its proximity operator with respect to J [19], described byMoreover, as mentioned in [18], define the conjugate of the function ϑ as the following form:and denote an equivalent characterization of and byfor and . 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 , , , , , , and . If is a solution of model (17), then for any and , there exists a vector such thatConversely, if there exist , ., , and satisfying (24) and (25), then is a solution of model (17).

Proof. Suppose that is a solution of model (17), by Fermat’s rule and the Chain rule of the subdifferential, we find that for all , the following inclusion relation holdsHere, there exists such thatFollowing the characterization (23) for , we have . Moreover, by the relation (21), we know that for any , , which yields (26). For any , multiplying the left-hand side of (27) by and using (21) lead to equation (24).
Conversely, if these exist and such that 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 is defined at a vector as and the operator is denoted at byand denotewith two identity matrices and . In particular, let andand we find that the operator is considered as the proximity operator of the function with respect to E, given by . As indicated in [19], the operator is firmly nonexpansive with respect to E, which implies that for all ,

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 equationmay not yield a convergence sequence because of the expansiveness of the matrix G. Here, . To overcome the difficulty, motivated by [19], we study an implicit iteration form through employing a splitting strategy to the matrix G, denoted bywithWe further rewrite (34) as the following preconditioned fixed-point proximity schemeThe 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 , there exists a unique such thatWe define the operator asThis means that the implicit iterative scheme (34) can be characterized as an explicit iteration . 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 letand find that its gradient ∇ϕ is Lipschitz continuous with the constant . Furthermore, let for the above matrices E in (31) and in (35) and for the above parameter λ and the constant η. The lemma is described as follows.

Lemma 2. Let be the fixed-point set of the operator defined by (38). If for and η, two matrices and the parameter λ satisfythen for any , the sequence yielded by the iterative scheme (34) satisfies

Proof. We first prove the symmetric positive definiteness of the involved matrices W and for the following analysis. For W and , the symmetric positive definiteness of implies that of W, and is a symmetric matrix. We need to prove that is positive definite. Since is a positive definite matrix, letIt can be verified that for ,that is, and are congruent. This implies that is positive definite if and only if . Hence, W and are symmetric positive definite matrices.
We then verify that (41) holds as follows: Since the operator is firmly nonexpansive with respect to E, by (32), we have that for any . By and , it can be verified thatFor the last term in (45) with , we let and find thatsince the following inequality holdsdue to the convexity of ϕ and Litschitz continuity of its gradient. Multiplying (45) by 2 and combining it with (46) yieldMoreover, it can be verified thatRecalling 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 be the sequence yielded by the iterative scheme (36) for any initial and U be the fixed-point set of the operator defined by (38). If two matrices and the parameter λ satisfy (40), then the sequence converges to a fixed point of the operator and converges to a solution of model (17).

Proof. Since two matrices and the parameter λ satisfy (40), it can be verified that inequality (41) holds for the sequence by Lemma 2.
We first prove that the sequence converges a cluster point of the sequence. Using inequality (41), we find thatfor any . This means that the sequence is bounded. As a result, there is a subsequence that converges to a cluster point . Summing inequality (41) for k from 0 to , we have thatwhich implies thatHence, the sequence 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 , the cluster point is a fixed point of and .
As a result, the sequence converges to a fixed point of , and 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 (Figure 1(b)) with 120 angles and 729 bins from clinic CT image with size (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.

(a)

(b)

(a)

(b)

(c)

(d)

We then apply different methods to reconstruct the images (Figure 2) 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 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 , where is the matrix with given in (8) and is the Kronecker product of matrices and . As indicated in [23], for the implementation of the scheme we may let , , and where . Moreover, as mentioned in [15, 17], we have that for ,and for ,

We then set and and choose for the noise-free case and for the case with noise in the experiment. Iterations stop whenfor the iterative sequence .

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 (Figure 3). NMSE is defined bywhere X is the reference image with size and Z is the reconstructed image. The smaller the NMSE, the better the reconstructed image. PSNR is denoted bywith mean squared errorwhere is the maximum value of X. The 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., ) for such two reconstruction methods and set different values of order α (i.e., ) 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).

(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 Figure 4 shows that the TFV and the TV yield the best PSNR and NMSE when , and the former is superior to the latter. Furthermore, the TFV with leads to the best PSNR for the developed method when and . 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.

Data Availability

The involved data in this paper are CT images.

Conflicts of Interest

The authors declare that there are no conflicts of interest.

Acknowledgments

This work was supported in part by the China Department of Science and Technology under Key grants 2016YFB0200602 and 2018YFC1704206, Natural Science Foundation of China under grants 11401601 and 81971691, Science and Technology Innovative Project of Guangdong Province under grants 2016B030307003, 2015B010110003, and 2015B020233008, Project of Guangzhou Innovation Leader Team (201809010012), Guangdong Provincial Science and Technology under Key grant 2017B020210001, Guangzhou Science and Technology Creative under Key grant 201604020003, and Construction Project of Shanghai Key Laboratory of Molecular Imaging (18DZ2260400).

References

L. Epstein, Introduction to the Mathematics of Medical Imaging, Society for Industrial and Applied Mathematics, Philadelphia, PA, USA, 2007.
L. I. Rudin, S. Osher, and E. Fatemi, “Nonlinear total variation based noise removal algorithms,” Physica D: Nonlinear Phenomena, vol. 60, no. 1–4, pp. 259–268, 1992.
View at: Publisher Site | Google Scholar
A. Sawatzky, “(Nonlocal) total variation in medical imaging,” Westfälische Wilhelms-Universität Münster, Münster, Germany, 2011, Ph.D. thesis.
View at: Google Scholar
K. Bredies, K. Kunisch, and T. Pock, “Total generalized variation,” SIAM Journal on Imaging Sciences, vol. 3, no. 3, pp. 492–526, 2010.
View at: Publisher Site | Google Scholar
S. Niu, J. Ma, J. Huang et al., “Second order total generalized variation for low-dose computed tomography image reconstruction,” in Proceedings of the IEEE International Symposium on Biomedical Imaging, pp. 173–176, Beijing, China, April 2014.
View at: Google Scholar
C. Miao and H. Yu, “A general-thresholding solution for ℓp (0 < p < 1) regularized CT reconstruction,” IEEE Transactions on Image Processing, vol. 24, no. 12, pp. 5455–5468, 2015.
View at: Publisher Site | Google Scholar
M. Sadiq, J. Simmons, and C. Bouman, “Model based image reconstruction with physics based priors,” in Proceedings of the IEEE International Conference on Image Processing, pp. 3176–3179, Phoenix, AZ, USA, September 2016.
View at: Google Scholar
Y. Chen, B. Guo, and X. Ma, “Image processing based on regularization with fractional calculus,” Mathematica Numerica Sinica, vol. 39, no. 4, pp. 393–406, 2017.
View at: Google Scholar
Y. Pu, J. Zhou, and X. Yuan, “Fractional differential mask: a fractional differential-based approach for multiscale texture enhancement,” IEEE Transactions on Image Processing, vol. 19, no. 2, pp. 491–511, 2010.
View at: Google Scholar
J. Zhang and K. Chen, “A total fractional-order variation model for image restoration with nonhomogeneous boundary conditions and its numerical solution,” SIAM Journal on Imaging Sciences, vol. 8, no. 4, pp. 2487–2518, 2015.
View at: Publisher Site | Google Scholar
J. Zhang and Z. Wei, “A class of fractional-order multi-scale variational models and alternating projection algorithm for image denoising,” Applied Mathematical Modelling, vol. 35, no. 5, pp. 2516–2528, 2011.
View at: Publisher Site | Google Scholar
I. Podlubny, “Matrix approach to discrete fractional calculus,” Fractional Calculus and Applied Analysis, vol. 3, no. 4, pp. 359–386, 2000.
View at: Google Scholar
I. Podlubny, A. Chechkin, T. Skovranek, Y. Chen, and B. M. Vinagre Jara, “Matrix approach to discrete fractional calculus ii: partial fractional differential equations,” Journal of Computational Physics, vol. 228, no. 8, pp. 3137–3153, 2009.
View at: Publisher Site | Google Scholar
S. Das, Functional Fractional Calculus, Springer Science & Business Media, Berlin, Germany, 2011.
C. A. Micchelli, L. Shen, and Y. Xu, “Proximity algorithms for image models: denoising,” Inverse Problems, vol. 27, no. 4, pp. 045009–045038, 2011.
View at: Publisher Site | Google Scholar
R. Chan, T. Chan, and A. Yip, “Numerical methods and applications in total variation image restoration,” in In Handbook of Mathematical Methods in Imaging, Springer-Verlag, New York, NY, USA, 2011.
View at: Google Scholar
A. Krol, S. Li, L. Shen, and Y. Xu, “Preconditioned alternating projection algorithms for maximum a posteriori ECT reconstruction,” Inverse Problems, vol. 28, no. 11, Article ID 115005, 2012.
View at: Publisher Site | Google Scholar
H. Bauschke and P. Combettes, Convex Analysis and Monotone Operator Theory in Hilbert Spaces, Springer, New York, NY, USA, 2011.
Q. Li, L. Shen, Y. Xu, and N. Zhang, “Multi-step fixed-point proximity algorithms for solving a class of optimization problems arising from image processing,” Advances in Computational Mathematics, vol. 41, no. 2, pp. 387–422, 2015.
View at: Publisher Site | Google Scholar
Y. Chen, J. Huang, S. Li, and Y. Xu, “A content-adaptive unstructured grid based integral equation method with the TV regularization for SPECT reconstruction,” Inverse Problems & Imaging, vol. 14, no. 1, pp. 27–52, 2020.
View at: Publisher Site | Google Scholar
Z. Wu, S. Li, X. Zeng, Y. Xu, and A. Krol, “Reducing staircasing artifacts in SPECT reconstruction by an infimal convolution regularization,” Journal of Computational Mathematics, vol. 34, no. 6, pp. 626–647, 2016.
View at: Publisher Site | Google Scholar
R. Siddon, “Fast calculation of the exact radiological path for a three-dimensional CT array,” Medical Physics, vol. 27, no. 12, 1985.
View at: Google Scholar
Y. Chen, Y. Lu, X. Ma, and Y. Xu, “Regularized CT reconstruction on unstructured grid,” Proceedings of the Medical Imaging 2016: Physics of Medical Imaging, vol. 97834G, pp. 1–6.
View at: Google Scholar

Copyright

Copyright © 2020 Yun Chen et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

PDF Download Citation

Download other formats

Order printed copies

Views

300

Downloads

738

Citations