CT Image Reconstruction via Nonlocal Low-Rank Regularization and Data-Driven Tight Frame

: X-ray computed tomography (CT) is widely used in medical applications, where many efforts have been made for decades to eliminate artifacts caused by incomplete projection. In this paper, we propose a new CT image reconstruction model based on nonlocal low-rank regularity and data-driven tight frame (NLR-DDTF). Unlike the Spatial-Radon domain data-driven tight frame regularization, the proposed NLR-DDTF model uses an asymmetric treatment for image reconstruction and Radon domain inpainting, which combines the nonlocal low-rank approximation method for spatial domain CT image reconstruction and data-driven tight frame-based regularization for Radon domain image inpainting. An alternative direction minimization algorithm is designed to solve the proposed model. Several numerical experiments and comparisons are provided to illustrate the superior performance of the NLR-DDTF method.


Introduction
Medical imaging applications of X-ray computed tomography (CT) include cranial, chest, cardiac abdominal and pelvic imaging. The wide use of CT is due to the ability of viewing interior structures without destroying the surface of organs or subjects. However, the radiation caused by X-ray during imaging does harm to patients' health. Thus, the approximately lossless image reconstruction based on low-dose X-ray, such as reducing the number of projections, is continuously considered by scientists and doctors. The basic model of the X-ray reconstruction problem can be represented as a linear inverse problem, where P ∈ R m×n is a projection operator representing the collection of discrete line integrations at different projection angles and along different beamlets. Numerically P is realized by Siddon's algorithm [1,2]. When small number of projections are used, i.e., m < n, the matrix P has a rank deficiency. As a result, Equation (1) will have an infinite amount of solutions. There are many methods for deriving a useful result from them, such as the solutions by the filtered back-projection (FBP) algorithm [3] and algebraic reconstruction technique (ART) [4]; however, they usually contain artifacts and thus result in unreliable reconstructions.
In practice, the projection image f usually contains noise. Gaussian noise [1,[5][6][7][8][9][10][11][12][13][14] and Poisson noise [15,16] both are considered to simulate a realistic environment. Thus, the basic model (1) is modified by f = Pu + . (2) In order to suppress noise and artifacts (local irregularities) during image reconstruction, it should usually contain tiny structural or high-frequency information to maintain fidelity and effectiveness. Various differential operator-based regularization methods, such as the total variation (TV)-based regularization approach (see [6,7,15,17] and references therein), low-rank models (see [8,18]) and wavelet frame-based methods (see [9][10][11][12]), have been proposed and received acceptable results for both, in vision and computational error. A standard form of the TV-based image reconstruction model can be written as It is important to mention that the authors of [9] established a rigorous connection between a special model of wavelet frame-based approach, called the analysis-based approach, and variational models.
Two adjacent columns in the projected image f represent the information collected from adjacent projection angles. Therefore, f must be column-dependent locally. In order to exploit this kind of prior knowledge of f , the authors of [1] proposed an image inpainting method on the projection image f , following with another inverse problem for image restoration in the spatial domain. This wavelet frame-based regularity with Radon domain inpainting was shown to be useful for reconstructing high quality images from a very small number of projections, where R Λ C denotes the reconstruction on Ω\Λ, and R Λ denotes the restriction on Λ. Here, f 0 is the projection image defined on the grid Λ of size N 0 × N p , where N 0 is the total number of detectors and N p is the number of angular projections. The unknown projection image f is defined on a grid Ω ⊃ Λ of size N 0 × N p , with N p = 2N p is selected in this paper. The first three terms of Equation (4) represent the data fidelity, making sure that f is consistent with f 0 on Λ and P f ≈ f . W 1 and W 2 are two different tight wavelet frame transforms. A fast algorithm was introduced to solve the model (4) based on the split Bregman algorithm [17,19,20] and the augmented Lagrangian method [21,22]. For more details about model and algorithm, we refer to [1]. Recently, it is known that in area of image restoration, data-driven tight frames or bi-frames generally outperform the regular pre-constructed wavelet frames [5,13,14]. The author of [5] proposed the SRD-DDTF model using data-driven tight frames as their sparsity priors for both u and f , where tight frames W 1 and W 2 and corresponding coefficients v 1 and v 2 , are all treated as unknowns for enforcing the sparsity approximation to u and f . Here, 0 norm v i 0 stands for the number of non-zero elements in v i , i = 1, 2. For solving the additional W i and v i subproblems such as Equation (7), an alternative optimization algorithm was proposed in [13] based on Singular Value Decomposition (SVD) of a matrix associated with the image u (or f ) after block rearrangement. As usual, the 0 regularization of v i can be solved by applying hard-thresholding to the singular values. A short review is stated in Section 2.1, for more details and convergence analysis, we refer to [5,13,14]. It is desired that there exist lots of similar structures both in natural and medical images. This property of nonlocal self-similarity was introduced by Dong et al. in [23,24] for exploiting structured sparsity in compressive sensing of both photographic and CT/MRI images. Random sampling and pseudo-radial sub-sampling are used for under sampling. However, they did not use any inpainting method in the transform domain for reconstruction. Xie et al. [16] used the nonlocal self-similar constraint in the position emission tomography (PET) image reconstruction with Total Variation(TV) method. The nonlocal low-rank approximation and TV model were also used for reducing noise in the denoising problem of low-dose CT image [25]. Xu et al. [26] proposed a model based on nonlocal low-rank and prior images under a given wavelet framework for reducing limited-angle artifacts in CT image restoration. There are also many other applications of nonlocal low-rank methods in the field of image denoising, one may see [27] for details.
In this paper, the nonlocal patch-based method will be used for restoration of CT image u. Roughly speaking, we shall use low-rank approximation for several image block patches after block matching. The rest of this paper is organized as follows. In Section 2, we will review image restoration methods based on data-driven tight frame and nonlocal low-rank regularization. In Section 3, we shall introduce our CT image restoration model combining data-driven tight frames for f and nonlocal low-rank regularization for u. In Section 4, several experiments will be provided for validating the merits of our proposed method, and the conclusion is discussed in Section 5.

Data-Driven Tight Frames
In [1,5], the sparse approximation methods were proposed for both the CT image u and projection image f , which are based on the pre-constructed wavelet frames or datadriven tight frames. When implemented, the data-driven tight frames are constructed based on the SVD of a matrix associated with the image u (or f ) after block rearrangement.
It is well known that compared with a basis, using a frame results in a more robust representation of signals [28,29]. In discrete setting, the fast decomposition transform W and fast reconstruction transform W T , two finite matrix operators with entries from a filter {a j } m j=0 associated to a translation-invariant tight wavelet transform, can be formed by convolution operators. Actually, with convolution operator S a : 2 (Z) → 2 (Z), The tight frame property is guaranteed by WW T = I as usual, see [5]. The data-driven tight frame method takes the tight frame W as an unknown, which is also determined stepwisely by solving the following optimization problem: We also use a small-size W, thus, the phantom image f should be reshaped as the matrix F ∈ R (N 1 N 2 )×p . Here, each column is a vectorization of N 1 × N 2 patches extracted from f . The filters {a j } m j=0 used for W are selected as the columns of the matrix D ∈ R (N 1 N 2 )×m . By denoting V ∈ R (N 1 N 2 )×p as the tight frame coefficients and taking m = N 1 N 2 , Equation (7) can be rewritten as In the sequel, we follow the alternative optimization algorithm for solving Equation (8), see [5,13] for details, with X and Y obtained by taking SVD of F(V k ) T , and T λ being the hard-thresholding operator defined by otherwise. (10)

Nonlocal Low-Rank Regularization
The nonlocal low-rank regularization is used for compressed sensing [24], PET image reconstruction [16] and CT image restoration from limited-angle projection data [26]. All these problems are handled by the alternative direction multiplier method technique, which is converted into several sub-problems and solved iteratively and alternately. When implemented, the rank-minimization sub-problem is also solved by the singular value thresholding (SVT) method, and the image reconstruction sub-problem is solved by different regularization methods.
In this article, the nonlocal low-rank regularization is used for CT image restoration, for the purpose of exploiting the nonlocal property in the image. After extracting √ n 1 × √ n 1 patches at position i in CT image u, denoted by x i ∈ R n 1 , m 1 -nearest patches based on Euclidean distance are used for obtaining the matrix Since in the natural image there exist lots of nonlocal self-similar structures, the formed data X i can be split into two parts, the low-rank matrix L i and the Gaussian noise matrix W i , i.e., X i = L i + W i . Then, L i can be recovered by the low-rank regularization model, where σ 2 W denotes the variance of Gaussian noise. In order to solve Equation (11), as in [24], we consider an approximate version, where φ(L i , σ k ) = ∑ n 0 r=0 σ k ,r σ r denotes the weighted nuclear norm with σ k ,r = 1 σ k r + and n 0 = min{n 1 , m 1 }. The singular values σ r and σ k r corresponding to L i and L k i , respectively, are ordered in descending order. We adopt the weighted singular value thresholding method to solve Equation (12) efficiently and effectively, where UΣV T is the SVD of X i , (x) + = max{x, 0} and σ k =

CT Image Reconstruction Model
In this paper, we also focus on the column number of the unknown projection image satisfying N p = 2N p , with N p being the column number of f 0 . f is considered to be restored from f 0 with doubled angular sampling resolution. Unlike the symmetric model (5), we use nonlocal low-rank regularization instead of the data-driven tight frame method when dealing with CT images. Our asymmetric CT image reconstruction model based on nonlocal low-rank regularity and data-driven tight frames reads as follows, where R i u = R i,0 u, R i,1 u, · · · , R i,m 1 −1 u denotes the first m 1 closest patches under cosine similarity for every exemplar patch x i . In this section, we demonstrate that the proposed objective functional can be efficiently solved by the method of alternative minimization. After k-th iteration, the variables of step k + 1 are determined by the following sub-problems.
(1) f sub-problem If all the parameters but f are determined by k-th iteration, Equation (14) can be practically rewritten as Here, we add a special term f − f k 2 2 in order to ensure that the new f k+1 is not far from f k , which can theoretically justify the convergence of this algorithm. By differentiating the right side of Equation (15) with respect to f , we obtain the following closed-form solution where [R Λ C + κR Λ + (µ + a)] is a diagonal matrix and, hence, Equation (16) can be computed easily.
(2) W and v sub-problem In this paper, we use the data-driven tight frame method explained in [5] or Section 2, for W and v sub-problems, As discussed in the previous section, the variables f , v and W are reformulated in the practice calculation. To be more specific, we need to solve the following reformulated sub-problems: Thus, to solve sub-problems (17) and (18), we can simply compute and where T a (·) is the hard thresholding operator given in Equation (10).
(3) L i and u sub-problems Based on the augmented Lagrangian method, we formulate the L i sub-problem as, Here, σ r and σ k r are the r-th ordered singular value with respect to L i and L k i . Applying the weighted singular value thresholding operator to Equation (21) yields the algorithm where U ΣV T is the SVD of R i u k + eL k i . After L i being updated, the current estimation of CT image u can be solved by our proposed model with other variables fixed, where the term u − u k 2 2 stands for the similarity in two adjacent steps. It is a quadratic optimization problem admitting a closed-form solution, In order to derive a much faster algorithm for solving Equation (23), we also use the alternative direction multiplier method (ADMM) [30]. Thus, we obtain where z ∈ R n×n is an auxiliary variable, u ∈ R n×n is the Lagrangian multiplier, and β is a positive scalar parameter. The split iteration version of the optimization of Equation (24) can be stated as follows: Here, z k+1 admits a closed-form solution, For u k+1 , we also use the conjugate gradient method to solve a linear equation, Collecting all together, we derive the following Algorithm 1.

Experiments
In this section, we report the experimental results of the proposed asymmetric lowrank approximation and data-driven tight frame-based asymmetric CT image recovery method. The proposed joint regulation method is actually based on the SRD-DDTF model [5]. The major difference is the optimization method for updating u k . Thus, we select the same warm-started initial value u 0 as in [5], given by analysis of the model in Section 3.1. In the following, we consider three phantoms, head image with smaller details from [5], head image and brain image with much more details from [24]. All experiments are conducted on a 3.20 GHZ Intel(R) Core(TM) i7-8700cpu with 16 GB memory. The PSNR, relative error and correlation results of the reconstructed images are included in Table 1, additionally, the visual effects are reported in Figures 1-3. Both of them reveal that the nonlocal low-rank method does preserve more information from the underground truth other than noise, which contains less structural information.
The degrade projection image f 0 is constructed by the every-other-angular extraction method from synthesized data, based on the Radon transform of real data u and the Monte Carlo simulation, f 0 = R Λ (Pu + ), where is some Gaussian white noise and the projection number N p is considered to be 15,30,45,60 in all examples. We choose the standard deviation of noise as max(|Pu|)/300 in our experiments. The projection image f is calculated based on the photon numbers the detectors received. Gaussian white noise is also generated when the corresponding optical equipment counts photons, which is related to the uncertainty of photons. For an initial value of f , we directly use Pu 0 as in [5].
The patch size for u is 6 × 6 and the patch size for f is 8 × 2 to adapt to the different shape of u and f . For optimization of the W and v subproblems, we use µ = 8400λ 1 and λ 1 is selected the same as in the SRD-DDTF model (5). When splitting the image u into overlapped small patches, we extract a paragon patch in every N s pixels along both horizontal and vertical directions. Here, we use N s = 5 and N s = 3 for realizing our proposed NLR-DDTF model. For each patch of u, m 1 = 41 similar patches are selected when the low-rank approximation is executed to derive L i . In our application, we choose η = 0.8 for smaller details in the image "head1" and η = 1.6 for the many more details in images "head2" and "brain". The regularization parameter λ 2 is selected to be 1.
For quantitative accuracy evaluation, the relative error and correlation for the reconstructed u corresponding to the ground truth image u t are defined as: with u and u t being the mean of u and u t . The maximum iteration number of all our experiments is set to 1000 with the following stopping criterion: whenever the relative change between two consecutive iterations is below 10 −3 (or 10 −4 for "head1" with higher N p than 15).

Conclusions
In this work, we introduced a novel spatial-Radon domain CT image reconstruction model, which combines the spatial domain inverse problem based on the nonlocal lowrank approximation method and the Radon domain inpainting model using a data-driven tight frame-based regularization. The proposed model reconstructs high quality and high resolution CT images even when the images contain more details than the observed low-resolution projection images. We used the split Bregman algorithm for numerical simulations of the proposed model.
The success of the regularization method in CT image restoration, no matter TV-based or wavelet frame-based, lies in the property of local smoothness or structure sparsity of the ground truth image. It is shown in [5] that using data-driven tight frame as sparsity priors for Radon domain projection image preformed better than using pre-determined wavelet frame systems. The nonlocal low-rank prior can reflect the group sparsity characteristic of similar patches in CT images well. Our proposed CT image restoration method combines the advantages of the above two methods. The PSNR and subjective quality comparison shown in Table 1 clearly reveal that the proposed NLR-DDTF method performs better than just the tight frame in several cases. Compared with the SRD-DDTF method, our proposed NLR-DDFT model's shortcoming is the slightly higher time consumption. When N s = 5, the average time required for all examples is 42% longer, and the average time consumption is 120% times higher if N s = 3.
Author Contributions: Models, methodology and formal analysis, Y.S. and S.S.; conceptualization and data curation, F.X.; convergence and supervision, Y.L.; programming and experiments, Y.S. and F.X.; writing-original draft preparation, Y.S. and S.S.; writing-review and editing, Y.L., X.Y. and X.Z.; project administration, F.X., Y.L. and X.Y. All authors have read and agreed to the published version of the manuscript.