CT Image Sequence Restoration Based on Sparse and Low-Rank Decomposition

Blurry organ boundaries and soft tissue structures present a major challenge in biomedical image restoration. In this paper, we propose a low-rank decomposition-based method for computed tomography (CT) image sequence restoration, where the CT image sequence is decomposed into a sparse component and a low-rank component. A new point spread function of Weiner filter is employed to efficiently remove blur in the sparse component; a wiener filtering with the Gaussian PSF is used to recover the average image of the low-rank component. And then we get the recovered CT image sequence by combining the recovery low-rank image with all recovery sparse image sequence. Our method achieves restoration results with higher contrast, sharper organ boundaries and richer soft tissue structure information, compared with existing CT image restoration methods. The robustness of our method was assessed with numerical experiments using three different low-rank models: Robust Principle Component Analysis (RPCA), Linearized Alternating Direction Method with Adaptive Penalty (LADMAP) and Go Decomposition (GoDec). Experimental results demonstrated that the RPCA model was the most suitable for the small noise CT images whereas the GoDec model was the best for the large noisy CT images.


Introduction
Accurate display of anatomical detail, in the form of small structures, features and objects, is a key requirement of computed tomography (CT). These details carry important information that may be associated with distinction between normal and pathological diagnosis. However, blurring of the CT images often compromises the visibility of such fine details.
Generally, CT modality is an economical way to generate highresolution images of human anatomy. CT images show more complex details of bones structures, compared to magnetic resonance imaging (MRI) images. However, one side CT images are degraded by blurring due to many reasons, such as the imperfect resolution of the imaging system, data loss in acquisition, and acquisition noise, to name a few. These artifacts and data loss result in weak contrast and blurry organ boundaries. Furthermore, streak artifact around thick cortical bone results in degraded visibility of some soft tissue tumors present in low contrast regions. Therefore, contrast enhancement for CT images is critical to facilitate evaluating the extent of disease in soft tissue area [1]. Another, source of degraded CT image is the presence of noise, which adds further blur to organ boundaries.
Recently, many CT image deblurring methods have been studied to visualize miniature-sized features. Jiang et al proposed a deblurring method for spiral CT image based on edge signal-tonoise ratio in 2003 [2]. Their method obtained the standard deviationsof two-dimension Gaussian blur kernel by calculating the minimum value of the edge signal-to-noise ratio of the blurred CT image, and then deblurred the CT image using the EM algorithm. This deblurring method significantly improved the identification of cochlear CT details. In 2005, Wang made an improvement to Jiang's method using the Wiener filter in place of the EM algorithm [3], and accelerated the speed of deblurring algorithm without sacrificing the deblurring effect. Later, based on edge-to-noise ratio and constrained least squares, Zhou proposed a filter blind deblurring algorithm and an iterative blind deblurring algorithm respectively [4]. Both of them significantly improved the image quality, compared with Jiang's and Wang's methods. In 2010, Hussien also used the Wiener filter to restore CT soft tissue [1]. In 2012, Zohair presented a fast deblurring method for CT medical images using a novel kernel set [5].
All the deblurred methods mentioned above only deal with a single CT image without considering the context of image sequence. these methods can also deal with CT image sequences via one by one, which just likes to deblur repeat single image several times with the same operations and the information between adjacent frames are not considered. There may be timeconsuming. In this study, we will explore a new image sequence deblurred algorithm by using sparse representation and low-rank decompose model. For the reason that we always obtained CT image sequences from CT medical imaging system and the doctors judge the state of illness according to the changes of the CT image sequence, so we think that dealing with CT image sequence is more significant than dealing with a single CT image. Our method is also suitable for a single CT image, the difference is that we deblur the low-rank component working out the average low-rank image for sequence image but on own rank of a slice CT for single image. The result may slightly worse than using the average rank of image sequence if sequence image have big difference from slices. It is the fact that a single image can use information between adjacent slices. Recently, low-rank technology is used to decompose image sequence into a low-rank component and a sparse component. In low-rank model, noise and blur information is considered to exist in the sparse component of decomposition image, whereas original image information exists in the low-rank component [6]. Following a similar philosophy, CT sequence images can also be decomposed into the sparse and low-rank component, and restoration can be applied to the sparse and lowrank components, respectively. In this paper, a CT sequence image restoration algorithm is proposed based on low-rank decomposition. Three classic low-rank models have been assessed in the experiments; and comparisons were performed on different medical images.

Review
A CT medical imaging system can be modeled as a linear transformation of an incoming analog signal, which is corrupted by intrinsic measurement fluctuations or electronic noise. A blurred image g x,y ð Þis assumed to be generated by convolving the latent image f x,y ð Þ with a blur kernel h x,y ð Þwith a latent image, with additive noisen x,y ð Þ: where Ã denotes convolution operator [7][8][9][10].
The main cause of degradation in medical CT images is the imperfect resolution of the imaging system, and the degradation function of a two-dimensional medical CT image is generally approximated as an isotropic two-dimensional Gaussian Point Spread Function (PSF) [11]: whereAw0is a constant, sis the variance. Image restoration technique aims to recover an image from its observation, degraded with blur and noise. In order to restore the image, prior knowledge of the degradation and inverse filtering must be obtained. Direct inverse filtering is usually used. Unfortunately,it is invalid when noise is taken into consideration. Another, degradation process typically happens in CT medical imaging, which includes blurring caused by imperfect resolution of the imaging system and noise by random fluctuations in signal intensity. Therefore, Wiener filtering [12][13][14] is a better choice than inverse filtering since it takes noise into consideration for inverse filtering.
We can obtain the expression of the formula 1 ð Þin the frequency domain through the discrete Fourier transform: whereF u,v ð Þ,H u,v ð Þ,G u,v ð Þ and N u,v ð Þ are the discrete Fourier transform of the latent imagef x,y ð Þ, blur kernel h x,y ð Þ, blurred input imageg x,y ð Þ and noise n x,y ð Þ respectively. Wiener filtering restoration is a method that regards H u,v ð Þ as its transmission function, and generates the recovery image f Ã x,y ð Þ by minimizing the mean square error between the recovery image f Ã x,y ð Þand the latent imagef x,y ð Þ, namely f x,y ð Þ{f Ã x,y ð Þ ½ 2 ? min. The expression of the recovery image by wiener filtering is where : ½ denotes wiener filtering, H Ã u,v ð Þ is the conjugate matrix of H u,v ð Þ,p n u,v ð Þ and p f u,v ð Þ are the power spectrums of noise and latent image respectively. It is often hard to estimate the values of p n u,v ð Þandp f u,v ð Þ, and we use the following expression to approximate the wiener filtering restoration.
wherecis a constant, which is numerically takes the reciprocal value of the signal-to-noise ratio of the blur image. The recovery image can be further calculated according to equation 5 ð Þby the Fourier inverse transform:f x,y Our Method

Motivation
Generally, we obtained CT image sequences from CT medical imaging system, but the existing deblurred methods always focus on the single image. In fact, each CT image in the same sequence is similar with its adjacent frames or maybe with a little change. Another, it was shown that a low-rank matrix could be recovered from incomplete sampling of its entries by computing the matrix of minimum nuclear norm the data [15]. Inspired by the image decomposition in Ref. [6] and the video processing, CT image sequences are decomposed by sparse and low-rank here. We can find noise and blur information exist in the sparse component of decomposition image, whereas most of the original image information exists in the low-rank component. Thus, CT sequence images restoration can be operated both in sparse and low-rank components, respectively. Now, sparse and low-rank technology mainly includes three kinds of low-rank decomposition models, and each low-rank model has its emphasis. In this study, both RPCA and GoDec models will be used to deblur CT image sequences. More details can be seen from the following subsections. The frame of our work is shown in Figure 1.

Framework of sparse and low-rank matrix decomposition
Sparse and low-rank decomposition is widely used in background modeling and shadow removal; it usually focuses on image sequences due to their intrinsic sparse characteristics. Medical CT image sequence has strong correlation between adjacent frames, and it is pointed out that medical CT images are sparse [16], Inspired by the image decomposition in natural image deblurring, we applied the sparse and low-rank decomposition for CT image deblurring. In this section, we will briefly introduce the image decomposition with three low-rank models.
Model 1. X~LzS, whereX is the data matrix of the original image,Lis a low-rank matrix, Sis a sparse matrix.
By using thel 1 -norm as the proxy of sparsity and the nuclear norm as the surrogate for low-rank, as in statistics and image processing [15][16][17][18], the sparse and low-rank recovery can be accomplished by solving the following RPCA problem: k k Ã is the nuclear norm defined by the sum of singular values, lis a positive weighting parameter. This problem arises in many applications, such as image processing, web data ranking, and bioinformatic data analysis. Under surprisingly broad conditions, the RPCA problem can be exactly solved via convex optimization that minimizes a combination of the nuclear norm and the l 1 -norm. It has the advantages of high speed and precision. For the RPCA problem 6 ð Þ, we may apply the augmented Lagrange multiplier method by identifying: Then its Lagrangian function [19][20][21][22] is: whereY is the Lagrange multiplier, bis a positive scalar, S : Tdenotes the standard trace inner product and : k k 2 is the Frobenius norm. Its iterative scheme is: Formula 9 ð Þcan be solved by the following: This algorithm solves the RPCA problem through the Lagrange multiplier method, which recovers a low-rank matrix with an unknown fraction of its entries arbitrarily corrupted.
Model 2. X~XZzS, whereX is the data matrix of the original image, Zis the skinny SVD, Sis a sparse matrix.
The decomposition problem based on model 2 can be converted into the following optimization problem: Its augmented Lagrangian function [23][24][25][26][27][28] can be written as the following equation: where Y is the Lagrange multiplier, Aand B are linear mappings, f and gare convex functions,bw0is the penalty parameter, S : T denotes the inner product, : k k 2 is the Frobenius norm. Its iterative scheme is: whereA Ã is the adjoint ofA,g A w0 and g B w0are parameters. Penalty parameterbis updated according to the formula: b kz1~m in b max ,rb k ð Þ , whereb max is an upper bound of b k f g. The value ofris defined as wherer 0 §1is a constant. This is the LADMAP algorithm. It solves the low-rank representation problem through the alternating direction method, and linearizes the quadratic penalty term adaptively. Compared with the traditional alternating direction method, this method has a novel rule to update the penalty so that it converges fast, without the need to introduce auxiliary variables or to invert matrices during the operation process. It has found wide applications in computer vision and machine learning.
Model 3. X~LzSzG,rank L ð Þƒr,card S ð Þƒk, whereX is the data matrix of the original image,Lis a low-rank matrix, Sis a sparse matrix, Gis a noise matrix, ris the rank ofX , kis the candinality. During the decomposition process, the algorithm alternatingly assigns the r{rankapproximation ofX {StoLand assigns the sparse approximation with cardinalitykofX {LtoS.
The decomposition problem based on model 3 can be converted into the following optimization problem [29][30][31][32][33][34]: Formula 15 ð Þ can be converted into the following two subproblems: which is equivalent to where This is the GoDec algorithm, which is a low-rank approximation method based on bilateral random projections. The algorithm can approximate the low-rank component, the sparse component and the noise component of the input data matrix. It has the advantage of small relative error and high decomposition speed.
To sum up, RPCA and GoDec share similar motivations in that they both explore the low-rank and sparse structures in the given data matrixX . They differ, however, by their distinctive ways of handling the specific decomposition but they are intrinsically different. RPCA assumesX~LzS (Sis sparse noise) and exactly decomposesX intoLandSwithout pre-defined rank(L) and card(S), that is, RPCA offers a blind separation of low-rank data and sparse noise, however, GoDec produces approximated decomposition of a general matrix X whose exact RPCA decomposition does not exist due to the additive noise Gand pre-defined rank(L) and card(S); GoDec directly constrains the rank range of Land the cardinality range of S, while RPCA minimizes their corresponding convex polytopes (i.e., the nuclear norm of Land l 1 -norm of S); GoDec usually produces less relative error with much less CPU seconds than RPCA. And the improvement of accuracy for the model of GoDec is due to more general than that of RPCA by considering the noise component. LADMAP and RPCA are same intrinsically, but LADMAP reduces the auxiliary variables and constraints, further accelerates the convergence speed and improves the accuracy of decomposition.

Restoration using low-rank matrix decomposition on CT image sequence
In this paper, sparse and low-rank matrix decomposition was applied to the problem of medical CT image sequence restoration. For the blurred CT image, it is necessary to recover its low-rank component, because it conveys major information on the latent image. The sparse component, on the other hand, normally reflects the majority of blurring. In practice, CT images differ in their characteristics: some contain high level of noise that seriously degrades image quality and affects visualization; in others, the noise level is relatively low and has little impact on the visual quality of the image or diagnosis of the lesions. Therefore, it is necessary for us to select the appropriate low-rank model based on the image characteristics, and to perform low-rank decomposition and restoration in order to obtain the desirable recovery results.
We propose a new CT image sequence restoration model based on low-rank decomposition as follows: whereX i is the i th slice of the original CT sequence, L i is a low-rank matrix, S i is a sparse matrix, n i is a noise matrix of X i , L Lis the average of all the low-rank matrix L i ,i §1, w,u,v are recovery functions, cis 0 or 1 according to the low-rank model.
Also, there exist several methods using the information of the adjacent frames, such as the methods introduced in ref. [35,36], but they use the idea that the corrupted portions of the projection data can be substituted with the corresponding portions from an unaffected adjacent slice, and restore the image one by one. Our method can restore the image sequence by one-time processing, which can not only consider the common characters of the image sequence but also take the different features of each image into account.
In natural imaging, the image blur can be regarded as atmospheric turbulence affecting the observed scene captured by the camera due to fluctuations of the refraction index of the medium [37]. This phenomenon may also exist in the medical CT imaging process caused by the refraction of the human tissue when the X-ray go into the human body, so we propose another degradation function we call turbulence PSF in this paper, which is acted on the sparse component for wiener filtering, the expression is as follows [38]: wherel is a constant. This degradation function we provided was inspired by the blurring natural image and then we think the turbulence also exist in CT imaging based on the reference [11]. Furthermore, after CT image sequences are decomposed into sparse and low-rank components, blur information has often concentrated on the sparse component of decomposition image and there is not much difference from low-rank components. Therefore, turbulence PSF is acted on the sparse component as well as Gaussian PSF is used on the average low-rank component. In order to demonstrate the performance of provided degradation function, we did some comparison experiments that we deblur a single CT image with only turbulence PSF, only Gaussian PSF and both the two PSFs, the result shows that use both the PSFs have more advantages.
The detailed steps of our algorithm are as follows: Algorithm: CT image sequence restoration based on sparse and low-rank decomposition Input: the blurred CT image sequence Output: the latent CT image sequence Step 1: Input the medical CT image sequenceX . X i represents thei th image in the sequence X ,i §1. IfX i is acolor image, convert it to gray scale; Step 2: Use all the CT images in the sequenceX to synthesize a high dimensional matrix, each column of the high dimensional matrix represents a single CT imageX i ; Step 3: select a reasonable low-rank model to decompose the high dimensional matrix in step 2 into high dimensional sparse matrix, low-rank matrix and noise matrix according to the actual situation of the CT image sequenceX ; Step 4: convert the high dimensional sparse matrix, low-rank matrix and noise matrix into sparse image sequenceS, low-rank image sequenceLand noise image sequence n, respectively; Step 5: Using our restoration model w X i ð Þ~u( L L)zv S i ð Þ zcn i to restore the blurred CT imageX i , whereX i~Li zS i zcn i . L i ,S i and n i are thei th images of image sequenceL,Sandn,respec-,respectively. For each slice image, when Gaussian PSF acts on L L, common knowledge of sequence images is used to restore a lowrank matrix of L i whereas sparse blurring component uses own S i matrix to restore by turbulence PSF.

Evaluation of Recover
In order to objectively evaluate the performance of the image restoration algorithm, we use standard deviation, image informa- Figure 5. Shows comparisons of different methods. From top to bottom: chest CT experiment, lung CT experiment, stomach CT (1) experiment, stomach CT (2) experiment. From left to right in each row: the original image, the recovery image using Wang's method, the recovery image using Hussien's method, the recovery image using our method. It needs to say, we select the RPCA model in chest CT and lung CT, and select the GoDec model in stomach CT. We can see that using our method to restore the medical CT image sequence can obtain satisfactory results, the recovery images can show high contrast, clear detail information and clear organ boundaries which are especially shown in the components circled with red ellipses in the images. doi:10.1371/journal.pone.0072696.g005 tion entropy and image quality measurement function value [39][40][41] as evaluation indices for comparison in this paper.
The formula of standard deviation is as follows where x i,j denotes the pixel value of location i,j ð Þof the recovery image, and the size of the recovery image is M|N.The greater the standard deviation value, the more dispersed the gray level distribution, the greater the image contrast, and the more detailed information.
Information entropy is one of the key indicators to measure how much image information is: the greater the value, the larger the amount of information contained in the image, its expression is where P : ð Þdenotes the probability. Image quality measurement function value is defined as The greater the image quality measurement function value, the more uniform the gray level distribution, and the better the image quality.

Experiments and Analysis
Based on the methods we described in Section 3, we investigate the most appropriate low-rank models for various types of CT images, differing in their noise levels and properties. Each image in our CT image sequences is 512 by 512 in size, with RGB color, obtained on a Discovery CT750 HD at the Beijing Cancer Hospital. We used s~1:3 in the Gaussian PSF, andl~0:0005 in the turbulence PSF for our experiments. All methods were implemented in Matlab7. 10.0 (R2010a), running on desktop with an Intel Core i3 CPU at 2.20GHz and 2GB memory, employing a 32-bit windows 7 operating system.

A. Comparison of different low-rank models
We first use the chest CT image sequence to evaluate the contrast of different low-rank models. The noise is small on each image of the CT sequence and does not affect the basic visualization. We decompose the input chest CT image sequence into a sparse component and a low-rank component in the RPCA model, and then recover the two components using wiener filtering, we usel~0:002in this experiment, and the decomposition criterion is to keep all the low-rank component as much as possible, the procedure is shown in Figure 2.
In order to analyze the scope of application of the three different low-rank models, we applied all decomposition models to the same chest CT image, and restored the chest CT image sequence. To illustrate the effectiveness of the proposed method in recovering the CT image sequence, the 10 th original image of the chest CT image sequence and its recovery images in different models are shown in Figure 3. For RPCA and LADMAP, we usedl~0:002andl~10respectively, and for GoDec,rank~2,card~8:15e 4 .
To further compare the recovery effect of the three models, we list in Table 1 the image evaluation index values mentioned above for every image in Figure 3, as is shown.
From Table 1, we can see that the recovery images using our method have great improvement in image standard deviation, information entropy, and quality measurement function value, compared with the original chest CT image. Although the image standard deviation and quality measurement function value of GoDec model is much greater, indicating sharper contrast, the recovery image has more serious ringing effect. Therefore, for the medical CT image sequence containing a low level of noise, the RPCA model is preferable to decompose the sequence, resulting in better recovery results.
The same comparative experiments were also performed on a stomach CT image sequence, which has a much larger noise component. The results of comparison are shown in Figure 4 and the corresponding evaluation index values of every image are  shown in Table 2. For RPCA and LADMAP, we usedl~0:002 and l~5respectively, and for GoDec, rank~1,card~8:15e 4 . From Table 2, we can see the recovery images using our method have a great improvement in image standard deviation, information entropy and quality measurement function value compared with the original CT image. Although the difference among the evaluation index values of the three recovery images is small, the recovered image based on the GoDec model shows less noise than the results from the other two methods. Therefore, for the CT images with large noise, the GoDec model is more preferable to decompose the image sequence.
The four groups of comparison of evaluation index values reveal that our method have advantages in all the evaluation index values when deal with the low level noise CT image. But the information entropy values of our method is lower than some of the others in case of high level noise CT images. What we want to clear is that there is no more other effective evaluation index to address this problem so far, So we used information entropy to try it on. Overall, our method is more effective in CT image sequence restoration than the alternative methods.

Conclusion and Prospect
This paper proposed image sequence restoration for CT scan image using Wiener filtering, based on sparse and low-rank decomposition. One of the key points to this method is uniting the CT image restoration under the framework of low-rank models: we first decompose the CT image sequence into sparse component and low-rank component in the RPCA model if noise level is low, and then restore both the two components using wiener filtering; in the presence of high level noise, it is better to select the GoDec model to decompose the images into three components, namely a sparse component, a low-rank component and a noise component, and then restore the sparse component and low-rank component using wiener filtering. It takes about 1.5s for our algorithm to recover a CT image of size 512|512. Experiments demonstrated that our method generated better results, compared with the other methods. In terms of medical CT image sequence restoration, our method is highly practical, with clear organ boundary and detailed information in the recovered image. The recovered images also have good soft tissue visibility and image quality. One limitation of our method is the reliance on matrix decomposition: in particular, if the low-rank component of every image in the sequence has a large variation, then our method does not work well. In this paper, both the parameters of the blur kernel and the noise-signal ratio of wiener filtering take empirical values, in the future work, we will seek other paths to determine these parameter values automatically and to explore more accurate blur kernels for the CT images.