A Rigid Motion Artifact Reduction Method for CT Based on Blind Deconvolution

: In computed tomography (CT), artifacts due to patient rigid motion often signiﬁcantly degrade image quality. This paper suggests a method based on iterative blind deconvolution to eliminate motion artifacts. The proposed method alternately reconstructs the image and reduces motion artifacts in an iterative scheme until the di ﬀ erence measure between two successive iterations is smaller than a threshold. In this iterative process, Richardson–Lucy (RL) deconvolution with spatially adaptive total variation (SATV) regularization is inserted into the iterative process of the ordered subsets expectation maximization (OSEM) reconstruction algorithm. The proposed method is evaluated on a numerical phantom, a head phantom, and patient scan. The reconstructed images indicate that the proposed method can reduce motion artifacts and provide high-quality images. Quantitative evaluations also show the proposed method yielded an appreciable improvement on all metrics, reducing root-mean-square error (RMSE) by about 30% and increasing Pearson correlation coe ﬃ cient (CC) and mean structural similarity (MSSIM) by about 15% and 20%, respectively, compared to the RL-OSEM method. Furthermore, the proposed method only needs measured raw data and no additional measurements are needed. Compared with the previous work, it can be applied to any scanning mode and can realize six degrees of freedom motion artifact reduction, so the artifact reduction e ﬀ ect is better in clinical experiments. lower RMSE mean values. standard deviations of the proposed method are other proposed method better stability. The simulation results of slight motion show that the proposed method has better motion artifact reduction effect and


Introduction
As one of the important technologies in medical diagnosis, a CT image can achieve a high performance in detecting and measuring small lesions [1,2]. However, CT has a problem that the reconstructed image will suffer from motion artifacts if a moving object is reconstructed without motion correction, in severe cases resulting in false diagnosis [3]. To reduce motion artifacts, the main approaches are shortening the scan time [4,5], external motion monitoring techniques [6][7][8], and motion estimation and compensation methods [9]. Since the third category of methods has the advantages of neither increasing hardware cost or design difficulty nor requiring additional devices, it has been widely studied. Among these methods, some were intended for 2D parallel-beam or fan-beam geometries and needed to assume a motion model that is an approximation of the real motion [9][10][11]. In [9], a novel method based on the Helgason-Ludwig consistency condition (HLCC) for estimation of rigid motion in fan-beam geometry was presented. Once the motion has been estimated, a compensation for the motion can be performed. This method has the disadvantages of large computational complexity and poor estimation accuracy. To solve this problem, a new method based on frequency domain analysis was proposed [10]. Motion parameters can be determined by the magnitude correlation of projections in frequency domain. This method was more accurate and faster on the performance of 1.
Richardson-Lucy (RL) deconvolution with SATV regularization is brought into the ordered subset expectation maximization (OSEM) iteration.

2.
With the proposed method, image reconstruction and motion artifact reduction are completed alternately in the iteration process. 3.
The simulation results are given to verify that the proposed method can be applied to any scanning mode.
The rest of this paper is organized as follows. Section 2 presents the related algorithm and the proposed motion artifact reduction method. In Section 3, the performance of the proposed method is verified. The results are discussed in Section 4. Finally, the conclusion is given in Section 5.

Ordered Subset Expectation Maximization (OSEM) Reconstruction
In order to correct the motion during the reconstruction, an iterative reconstruction algorithm is needed. OSEM is used as the reconstruction algorithm [14,22].
where µ n j is the estimated activity in voxel j after iteration n, y i is the measured sinogram, a ij is the system matrix element representing the probability that a photon emitted from voxel j is detected in detector pixel i, and S b is one subset.

Image Degradation Model
The generation of CT images can be modeled as an image degradation system. The motion-degraded image is formed by a convolution between the ideal image and the point spread function (PSF) [23].
where g(s) is the degraded image, f (s) is the ideal image, h(s) is the impulse response of the system or PSF, s denotes the spatial coordinates, and * is the convolution in spatial domain. In practice, the PSF is not precisely known. Hence, an iterative blind deconvolution approach is needed to restore the real image without complete knowledge of the associated system PSF.

RL Deconvolution with SATV Regularization
RL can effectively solve nuclear medical imaging problems, which is derived from the maximum likelihood of a Poisson distribution. The TV model was first proposed by Rudin et al. [24] in edge preserving image restoration, it is a popular choice of being the regularization term [25]. RL algorithm with total variation regularization (RLTV) was suggested in [26,27].
where λ is the regularization parameter, div(·) refers to divergence, · represents l2-norm, and ∇ denotes the gradient operator. In Equation (4), λ plays a very important role, which controls the TV regularization strength. It should be neither too small nor too large: if λ is too small, the noise will not be well suppressed; inversely, if it is too large, the edge and detailed information in image will be blurred. Therefore, a spatially adaptive TV regularization model considering the spatial dependent property is introduced [21,28]. the term in the denominator provide an estimate of edges from an initial image f 0 at scale σ using the Gaussian kernel G σ filtered gradients, and ε is a parameter of adaptive function λ(s). The spatially adaptive function can help reduce noise without producing the staircase effect because it can self-adjust according to the smoothed gradient of the image. To summarize, the SATV regularization has two good features: (1) For edge and detail, because λ(s) is small which weakens the TV regularization, then the edge and texture will be well preserved; (2) for flat region, a large λ(s) leads to a large TV regularization strength, so the noise will be well suppressed.

Motion Artifact Reduction in Reconstruction
To reduce motion artifacts, iterative blind deconvolution algorithm is inserted into the iterative process of OSEM reconstruction. Hence, Equation (1) can be modified as In Equation (6), it is assumed that the estimation of degraded image (( f n * h n )) is still a Poisson distribution. The corresponding projection of estimated motion-degraded image are compared with the measured projection data to estimate the ideal image and PSF. The overall process of the BD-OSEM method is shown in Figure 1.
the term in the denominator provide an estimate of edges from an initial image 0 f at scale σ using the Gaussian kernel G σ filtered gradients, and ε is a parameter of adaptive function ( ) λ s .
The spatially adaptive function can help reduce noise without producing the staircase effect because it can self-adjust according to the smoothed gradient of the image. To summarize, the SATV regularization has two good features: (1) For edge and detail, because ( ) λ s is small which weakens the TV regularization, then the edge and texture will be well preserved; (2) for flat region, a large ( ) λ s leads to a large TV regularization strength, so the noise will be well suppressed.

Motion Artifact Reduction in Reconstruction
To reduce motion artifacts, iterative blind deconvolution algorithm is inserted into the iterative process of OSEM reconstruction . Hence, Equation (1) can be modified as In Equation (6), it is assumed that the estimation of degraded image ( ( ) As shown in Figure 1, px represents input image size. The result of each OSEM-iteration is used as the degraded image for the BD-iteration ( 0 ( )  As shown in Figure 1, px represents input image size. The result of each OSEM-iteration is used as the degraded image for the BD-iteration (g 0 (S) = µ n j ), while the result of the BD-iteration is returned to the OSEM-iteration as intermediate value for the next iteration. To stop the iterations, a difference measure between two successive iterations is defined [29]. If the difference is smaller than a threshold, the computation is stopped, and the last estimation is the best one. The criterion is defined as the following: Algorithms 2019, 12, 155 5 of 12 In Figure 1, if τ k ≥ r, the result of the BD-iteration will be returned to the next BD-iteration. If τ k < r and τ n ≥ r, the result of the BD-iteration will be inserted to the iteration of OSEM, which is used as intermediate value. Note that the process of BD-OSEM is stopped at two or more iterations; hence, τ k < r and τ n < r should be satisfied.

Numerical Phantom Experiments
For the numerical simulations, the modified three-dimensional version of the Shepp-Logan phantom was used to simulate the CBCT system. The additional tiny circles in the 128th slice of modified 3D Shepp-Logan head phantom are used to simulate the small tissue, as shown in Figure 2a. The size of the phantom was 256 × 256 × 256 voxels and the display window was [0 HU, 1000 HU]. The scanning time was set to t m = 1. The interval angle was 1 • , and the number of projections was 360. The distance from X-ray source to detector was 100 cm and source to object was 50 cm. Poisson noise was added to the raw simulated data before the reconstruction (assuming that 10,000 photons were detected on each detector element in the blank scan).
In Figure 1, if k r τ ≥ , the result of the BD-iteration will be returned to the next BD-iteration. If k r τ < and n r τ ≥ , the result of the BD-iteration will be inserted to the iteration of OSEM, which is used as intermediate value.
Note that the process of BD-OSEM is stopped at two or more iterations; hence, k r τ < and n r τ < should be satisfied.

Numerical Phantom Experiments
For the numerical simulations, the modified three-dimensional version of the Shepp-Logan phantom was used to simulate the CBCT system. The additional tiny circles in the 128th slice of modified 3D Shepp-Logan head phantom are used to simulate the small tissue, as shown in Figure  2a. The size of the phantom was 256 × 256 × 256 voxels and the display window was [0HU, 1000HU]. The scanning time was set to 1 m t = . The interval angle was 1°, and the number of projections was 360. The distance from X-ray source to detector was 100 cm and source to object was 50 cm. Poisson noise was added to the raw simulated data before the reconstruction (assuming that 10,000 photons were detected on each detector element in the blank scan).
(a) (b) The projection data were obtained by applying the motion to a phantom during the scan. Figure  3 shows the motion segment that was used in the simulations. The slight motion of up 4° and 7 mm for rotations and translations obtained by Kim et al. [7] (a motion segment in Figure 15 of [7]) are presented in Figure 3a. Figure 3b shows the large motion exceeding 30° and 30 mm, which is much larger than that typically expected in patients [7].
In simulations, a parameter =1 σ and the smoothness parameter 6 =10 ε − were used for computations and the iterations were stopped when the difference between two images is less than a chosen threshold For the OSEM reconstruction, four iterations with 60 subsets were applied. The deconvolution were stopped at the third iteration for the BD-iteration. The projection data were obtained by applying the motion to a phantom during the scan. Figure 3 shows the motion segment that was used in the simulations. The slight motion of up 4 • and 7 mm for rotations and translations obtained by Kim et al. [7] (a motion segment in Figure 15 of [7]) are presented in Figure 3a. Figure 3b shows the large motion exceeding 30 • and 30 mm, which is much larger than that typically expected in patients [7].
In simulations, a parameter σ = 1 and the smoothness parameter ε = 10 −6 were used for computations and the iterations were stopped when the difference between two images is less than a chosen threshold r = 10 −4 . For the OSEM reconstruction, four iterations with 60 subsets were applied. The deconvolution were stopped at the third iteration for the BD-iteration. Figure 4 shows the 128th slice of reconstructed images obtained by four methods. In this paper, "OSEM-BD" refers the method based on OSEM reconstruction algorithm with motion correction by BD-iteration, and "RL-OSEM" represents the method based on the combination of RL deconvolution and OSEM reconstruction, that is the proposed method without SATV regularization. It can be seen that the motion-corrected image obtained by the proposed method is very close to the real phantom. Since the motion information is constantly superimposed in the reconstruction process, the OSEM-BD method cannot completely eliminate the motion artifacts. In addition, the RL-OSEM method cannot balance edge preservation and noise amplification very well because it has no regularization constraints, so there are some motion artifacts and intensity oscillations in the reconstructed image. The proposed and OSEM reconstruction, that is the proposed method without SATV regularization. It can be seen that the motion-corrected image obtained by the proposed method is very close to the real phantom. Since the motion information is constantly superimposed in the reconstruction process, the OSEM-BD method cannot completely eliminate the motion artifacts. In addition, the RL-OSEM method cannot balance edge preservation and noise amplification very well because it has no regularization constraints, so there are some motion artifacts and intensity oscillations in the reconstructed image. The proposed method can overcome the problems of RL-OSEM method. Thus, the reconstructed image obtained by the proposed method is the clearest and most of the motion artifacts are eliminated.  The effects of motion artifact reduction were evaluated by visual assessment and with quantitative analysis. The root-mean-square error (RMSE), Pearson correlation coefficient (CC), and mean structural similarity (MSSIM) were chose as the metrics. RMSE was calculated in HU, while CC and MSSIM were dimensionless [7]. Table 1 provides quantitative comparisons of all reconstructed images among four and OSEM reconstruction, that is the proposed method without SATV regularization. It can be seen that the motion-corrected image obtained by the proposed method is very close to the real phantom. Since the motion information is constantly superimposed in the reconstruction process, the OSEM-BD method cannot completely eliminate the motion artifacts. In addition, the RL-OSEM method cannot balance edge preservation and noise amplification very well because it has no regularization constraints, so there are some motion artifacts and intensity oscillations in the reconstructed image. The proposed method can overcome the problems of RL-OSEM method. Thus, the reconstructed image obtained by the proposed method is the clearest and most of the motion artifacts are eliminated.   The effects of motion artifact reduction were evaluated by visual assessment and with quantitative analysis. The root-mean-square error (RMSE), Pearson correlation coefficient (CC), and mean structural similarity (MSSIM) were chose as the metrics. RMSE was calculated in HU, while CC and MSSIM were dimensionless [7]. Table 1 provides quantitative comparisons of all reconstructed images among four The effects of motion artifact reduction were evaluated by visual assessment and with quantitative analysis. The root-mean-square error (RMSE), Pearson correlation coefficient (CC), and mean structural similarity (MSSIM) were chose as the metrics. RMSE was calculated in HU, while CC and MSSIM were dimensionless [7]. Table 1 provides quantitative comparisons of all reconstructed images among four methods. The mean and standard deviation of each metric were calculated over 256 reconstructed slices. The means (RMSE, CC, and MSSIM) indicated the overall quality of the reconstructed images. Standard deviations were taken as indices of the slice-to-slice variation in the calculated metrics. As expected, images reconstructed by using the proposed algorithm provide much higher CC and MSSIM mean values and much lower RMSE mean values. All the standard deviations of the proposed method are much lower than the other three methods, which demonstrate the proposed method has better stability. The simulation results of slight motion show that the proposed method has better motion artifact reduction effect than OSEM-BD and RL-OSEM method. To further evaluate the proposed method, the projection data were obtained by applying a large motion. Figure 5 shows the 128th slice of images reconstructed from projection with large motion. The metric values for four different methods are shown in Table 2. Three metrics of the proposed method are better than the other three methods. Nevertheless, the motion-corrected image obtained by the proposed method still has residual artifacts. Because the reconstructed image is badly damaged and motion artifacts are very serious, it cannot be corrected by the proposed method. The simulation results of large motion show that the proposed method has a poor motion artifact reduction effect, but it is still better than OSEM-BD and RL-OSEM methods.

Method
Metric To further evaluate the proposed method, the projection data were obtained by applying a large motion. Figure 5 shows the 128th slice of images reconstructed from projection with large motion. The metric values for four different methods are shown in Table 2. Three metrics of the proposed method are better than the other three methods. Nevertheless, the motion-corrected image obtained by the proposed method still has residual artifacts. Because the reconstructed image is badly damaged and motion artifacts are very serious, it cannot be corrected by the proposed method. The simulation results of large motion show that the proposed method has a poor motion artifact reduction effect, but it is still better than OSEM-BD and RL-OSEM methods.

Head Phantom Experiments
In the following simulations, the proposed method was verified on head phantom. A selected slice of head phantom was shown in Figure 2b. The sinogram was acquired from head phantom with motion and noise that were the same as the Shepp-Logan phantom experiments. The head CT image was acquired from a GE Hi-Speed multi-slice CT scanner. A total of 720 views were uniformly selected over a 360 • range under 150 mA X-ray tube current and 120 kVp. The display window was [−5 HU, 75 HU]. All of the reconstructed images (512 × 512 pixels) were reconstructed by the OSEM method with three iterations and 30 subsets. BD-iterations were stopped at the sixth iteration. The other experimental parameters were set the same as those of the Shepp-Logan phantom experiments. Figure 6 shows the selected slice of reconstructed images obtained by using OSEM, OSEM-BD, RL-OSEM, and the proposed method. To view the texture of reconstructed images more clearly, the zoomed images of a region of interest (ROI) from the corresponding images of Figure 6 are shown in Figure 7. Clearly, the proposed method is superior to the other methods in image quality. slice of head phantom was shown in Figure 2b. The sinogram was acquired from head phantom with motion and noise that were the same as the Shepp-Logan phantom experiments. The head CT image was acquired from a GE Hi-Speed multi-slice CT scanner. A total of 720 views were uniformly selected over a 360° range under 150 mA X-ray tube current and 120 kVp. The display window was [-5HU, 75HU]. All of the reconstructed images (512 × 512 pixels) were reconstructed by the OSEM method with three iterations and 30 subsets. BD-iterations were stopped at the sixth iteration. The other experimental parameters were set the same as those of the Shepp-Logan phantom experiments. Figure 6 shows the selected slice of reconstructed images obtained by using OSEM, OSEM-BD, RL-OSEM, and the proposed method. To view the texture of reconstructed images more clearly, the zoomed images of a region of interest (ROI) from the corresponding images of Figure 6 are shown in Figure 7. Clearly, the proposed method is superior to the other methods in image quality.   To quantitatively evaluate the proposed method for head phantom, all slices of reconstructed images were compared by RMSE, CC, and MSSIM, as shown in Table 3. Clearly, the proposed method provides the highest CC and MSSIM mean values and the lowest RMSE mean values. All the standard deviations of the proposed method are much lower than the other three methods. The results indicate that the proposed method has a better performance on motion artifact reduction.

Patient Scan Experiments
To verify the effectiveness of the proposed method on the clinical data with realistic motion, the method was applied to clinical experiments in which motion artifacts had been observed. The scan To quantitatively evaluate the proposed method for head phantom, all slices of reconstructed images were compared by RMSE, CC, and MSSIM, as shown in Table 3. Clearly, the proposed method provides the highest CC and MSSIM mean values and the lowest RMSE mean values. All the standard deviations of the proposed method are much lower than the other three methods. The results indicate that the proposed method has a better performance on motion artifact reduction.

Patient Scan Experiments
To verify the effectiveness of the proposed method on the clinical data with realistic motion, the method was applied to clinical experiments in which motion artifacts had been observed. The scan was performed on a GE Light Speed VCT scanner. The scan parameters were tube voltage 120 kVp and tube current 200 mA. Image size was 512 × 512 pixels and the display window was [−5 HU, 75 HU]. Because the patient bed did not move with the patient during the scan, the patient bed data were removed. Three iterations with 30 subsets were applied for OSEM reconstruction, and four iterations for the BD-iteration. The other experimental parameters were set the same as those of the Shepp-Logan phantom experiments. Figure 8 shows the selected slice of reconstructed images and the repeated scan image (which is done because of the observed motion in the first scan). The uncorrected image was reconstructed with the scanner system software, as shown in Figure 8a. Table 4 shows the comparisons of RMSE, CC and MSSIM. Note that all reconstructed images were registered to the repeated scan before the calculation of these metrics. Because the positional differences are irrelevant for image quality [14,30]. The results of patient scan experiments show that the reconstructed image obtained by the proposed method can effectively eliminate motion artifacts, and the artifact reduction effect of the proposed method is better than OSEM-BD and RL-OSEM methods.

Discussion
The performance of the proposed method was verified by using a modified 3D Shepp-Logan phantom, a head phantom, and patient scan. In numerical phantom experiments, the proposed method produced higher quality images for the case of slight motion than the other three considered approaches. As illustrated in Figure 4a,b, the images included severe motion artifacts over the whole image, and thereby inner structures and the skull were noticeably distorted. In Figures 4c,d the majority of distortions were eliminated. But many intensity oscillations appeared inside the phantom (see Figure 4c). The proposed method yielded an appreciable improvement on all metrics, reducing RMSE by about 30% and increasing CC and MSSIM by about 14% and 19%, respectively, compared to the RL-OSEM method (see Table 1). However, the proposed method did not perform well in cases of large motion, such as that of Figure 3b. Since the OSEM reconstruction was corrupted severely, these motion artifacts were too severe to be corrected by the proposed method (see Figure 5d). Therefore, a limitation of the method is

Discussion
The performance of the proposed method was verified by using a modified 3D Shepp-Logan phantom, a head phantom, and patient scan. In numerical phantom experiments, the proposed method produced higher quality images for the case of slight motion than the other three considered approaches. As illustrated in Figure 4a,b, the images included severe motion artifacts over the whole image, and thereby inner structures and the skull were noticeably distorted. In Figure 4c,d the majority of distortions were eliminated. But many intensity oscillations appeared inside the phantom (see Figure 4c). The proposed method yielded an appreciable improvement on all metrics, reducing RMSE by about 30% and increasing CC and MSSIM by about 14% and 19%, respectively, compared to the RL-OSEM method (see Table 1). However, the proposed method did not perform well in cases of large motion, such as that of Figure 3b. Since the OSEM reconstruction was corrupted severely, these motion artifacts were too severe to be corrected by the proposed method (see Figure 5d). Therefore, a limitation of the method is that the proposed method usually did not perform well when the amplitude of the rotations was more than 30 • and the amplitude of the translations was more than 30 mm.
In head phantom experiments, the reconstructed images showed marked improvement when the proposed method was applied. As shown in Figure 6a,b the reconstructed images exhibited severe motion artifacts and distorted image quality. And the image reconstructed by the RL-OSEM method still had artifacts and its edge appeared thinner than the true image (see Figures 6c and 7c). All the metrics also suggested improvement motion correction performance on the proposed method. The RMSE of the proposed method was the smallest and reduced by about 27% compared to the RL-OSEM method; the CC and MSSIM of the proposed method were the greatest and increased by about 15% and 20%, respectively, compared to the RL-OSEM method (see Table 3).
Comparative results on numerical phantom and head phantom show that the OSEM-BD method can eliminate noise but not reduce motion artifacts. The BD-OSEM method can effectively reduce motion artifacts and the noise in flat regions as well as preserve the edge and detailed information.
In patient scan experiments, the reconstructed images, even with motion correction, provided severe motion artifacts (see Figure 8b,c). Most of the motion artifacts in the reconstructed image were eliminated after the motion correction by the proposed method (see Figure 8d). The RMSE of the proposed method can reduce by about 30%, compared to the RL-OSEM method; the CC and MSSIM of the proposed method can increase by about 15% and 21%, respectively, compared to the RL-OSEM method (see Table 4).
The computation will take a long time since there are a large number of CT views in clinical scans. Currently, the computation complexity of the proposed method is mainly the OSEM reconstruction for a patient scan, which improves the computation efficiency compared to the method mentioned in [14] to some extent.
The performance of the proposed method has been verified by using a patient scan. However, this experiment may not be sufficient to demonstrate the algorithm for various clinical cases. In the future, more clinical experiments should be considered.

Conclusions
In this work, a method based on iterative blind deconvolution is developed to reduce the rigid motion artifacts for CT, which only requires the measured raw data. Since the OSEM is a fast iterative algorithm and widely used in CT image reconstruction, this paper uses OSEM as the reconstruction algorithm. As well as, RL iteration algorithm can estimate a clear image only based on original data and SATV regularization can effectively preserve image edge and detail information, while removing noise effects, so this paper uses RL algorithm with SATV regularization (RLSATV) to eliminate motion artifacts. Therefore, this paper combines the OSEM algorithm with RLSATV to complete image reconstruction and motion artifact reduction in the iteration process. In the process of iteration, the RL algorithm which is regularized using SATV is inserted into the iterative process of the OSEM reconstruction technique. The iteration process does not end until the difference measure between two successive iterations is smaller than a threshold. The proposed method has been evaluated by using phantom and patient scan experiments, which can reduce motion artifacts and provide high-quality images. Quantitative analysis results show that the RMSE of the proposed method can reduce by about 30% compared to the RL-OSEM method; the CC and MSSIM of the proposed method can increase by about 15% and 20%, respectively, compared to the RL-OSEM method. Additionally, the proposed method does not have to estimate the motion and is suitable for any scanning mode, thus it can become a valuable tool for motion compensation in clinical application.
Author Contributions: Y.Z. contributed to the overall idea of the algorithms, carried out the simulations, and drafted the manuscript; L.Z. supervised the work and provided the major direction of the research. All authors have read and approved the final manuscript.