Abstract
Purpose:
Low-dose computed tomography (LDCT) has promising potential for dose reduction in medical applications, while suffering from low image quality caused by noise. Therefore, it is in urgent need for developing new algorithms to obtain high-quality images for LDCT.
Methods:
This study tries to exploit the sparse and low-rank properties of images and proposes a new algorithm based on subspace identification. The collection of transmission data is sparsely represented by singular value decomposition and the eigen-images are then denoised by block-matching frames. Then, the projection is regularized by the correlation information under the frame of prior image compressed sensing (PICCS). With the application of a typical analytical algorithm on the processed projection, the target images are obtained. Both numerical simulations and real data verifications are carried out to test the proposed algorithm. The numerical simulations data is obtained based on real clinical scanning three-dimensional data and the real data is obtained by scanning experimental head phantom.
Results:
In simulation experiment, using new algorithm boots the means of PSNR and SSIM by 1 dB and 0.05, respectively, compared with BM3D under the Gaussian noise with variance 0.04. Meanwhile, on the real data, the proposed algorithm exhibits superiority over compared algorithms in terms of noise suppression, detail preservation and computational overhead. The means of PSNR and SSIM are improved by 1.84 dB and 0.1, respectively, compared with BM3D under the Gaussian noise with variance 0.04.
Conclusion:
This study demonstrates the feasibility and advantages of a new algorithm based on subspace identification for LDCT. It exploits the similarity among three-dimensional data to improve the image quality in a concise way and shows a promising potential on future clinical diagnosis.
Keywords
Introduction
Low-dose computed tomography (LDCT) remains an active topic [1, 2] when an ever-increasing attention has been paid to the hazards of radiation. However, when the dose gets lower, the image quality will be unavoidably degraded because of the surge in noise [3, 4]. Thus, LDCT is associated with the noise suppression in reconstructed image. To solve this problem, many investigations have focused on exploitation of prior knowledge and image similarity. Many methods have been proposed, which can be divided into three categories: sinogram restoration, image post-processing and iterative reconstruction (IR).
Among the three kinds of methods, there is considerable evidence that sinogram restoration has advantages in computational and noise-induced artifacts suppression [5]. The main idea is to restore the ideal sinogram by suppressing the noise in the transmission data before reconstruction. In the past decades, some typical methods have been proposed. For example, filter-based algorithms have been well applied such as bilateral filtering [6], adaptive filtering [7] and the non-local means filtering [8]. In 2014, Cui et al. suppressed the noise by minimizing an energy function consisting of an adaptive-weighted total variation (AWTV) [9]. In 2017, Xie et al. designed a maximum a posteriori probability (MAP) framework based on a new prior equation to restore the sinogram [10]. In 2019, Shangguan et al. proposed a sinogram restoration algorithm using the regularized Perona–Malik (P–M) equation with intuitionistic fuzzy entropy [11]. And in the recent years, deep learning methods are widely used in projection domain processing, such as GAN [12], self-supervised method [13], reinforcement Learning [14], which is discussed in the next few paragraphs. All above work plays a significant role in the research of LDCT. Nevertheless, the sensitivity of projection data may lead to spatial resolution degradation and newly generated artifacts, which limits the practical application.
Considering the difficulty in projection process, the exploitation of image post-processing is a natural idea. Generally, this means directly suppressing the noise in CT images. As it gets rid of the directly processing of projection data, the application of various natural image denoising methods is more convenient. For example, the wavelet [15] and non-local means (NLM) [16] were applied at the beginning. Further, in 2016, Michael et al. proposed the locally-consistent non-local means (LC-NLM), locally limiting spatial consistency [17]. In 2018, Hasan et al. proposed a blind source separation (BSS) based CT image method using a multiframe low-dose image sequence [18] and in 2021, they further applied this in Noise2Noise (N2N) model [19]. All above has increased the imaging accuracy and reduced the errors caused by the processing of transmission data. However, meanwhile, it is worth noting that the acquisition of noise distribution becomes difficult, which leads to over-smoothing and residual errors.
The above two kinds of methods have their own characteristics and both are sensitive to the development of LDCT. However, a critical problem is the two are inherently susceptible to noise, due to the dependence on transmission data. To overcome this, iterative reconstruction (IR) methods have been developed for decades. Different from the above two methods, IR combines the process of reconstruction and noise suppression into one step. It is seen as the solution of optimization problem which includes the basic models for imaging process and prior information [20–22]. As is known for all, much work has achieved great effect including the total variation (TV) [22], tight wavelet frames [24], dictionary learning [25] and so on. And as the development of compressive sensing (CS) [26], IR has been further improved [27–29]. Zhang et al. applied a tensor-based dictionary learning method to spectral CT and achieved better image quality [27]. Liu et al. significantly improved the image quality of LDCT by building a 3D feature dictionary [28]. Peng et al. applied convolutional sparse coding (CSC) to CT reconstruction with good results [29]. Nevertheless, the consumption of computational cost and the design of regularization terms are a bottleneck for its development.
Despite the progress in LDCT the above work has made, some problems including the methods for optimization, parameter setting and task-specific model design remain incomplete [30]. In recent years, deep learning (DL) develops rapidly, which inspires a new thinking in LDCT. DL imitates human visual neural network to design computer network to process information. And for LDCT, DL has become a hot research direction. In 2017, Chen et al. combined the auto-encoder, the deconvolution network, and shortcut connections into the residual encoder-decoder convolutional neural network (RED-CNN) and achieved a competitive performance [31]. Further, in 2019, Du et al. introduced the attention mechanism into the GAN network and demonstrated the effectiveness of the method in clinical images [32]. In 2021, Huang et al. proposed Du-GAN by combining U-net network and GAN network, which has a good enhancement effect on the edge of denoised images [33]. In 2022, Mehmet et al. exploited an approach that deep convolutional neural networks (CNNs) generate patterns easier than the noise with different loss function variants, which outperformed the other methods both qualitatively and quantitatively [34].
While reported promising results in previous works, DL methods remains to be controversy for several reasons. On the one hand, for medical CT imaging, deep learning still has some limitations [35, 36]. In 2020, Emil et al. proposed that the lack of injectivity and the high variability in images quality of DL brings a series of potential problems for image reconstruction and for medical images, real data are not easily obtained, which limits the application of deep learning [35]. On the other hand, referring to DL itself, the high demand for memory and the improvements in metrics remain to be solved [36].
Despite the surge interest in DL, there’s much potential in traditional methods to improve the image quality. As mentioned above, compared to DL, traditional methods have better interpretability and do not require a large amount of data for model training. In this paper, inspired with the application of CS [37], the purpose of this study is to propose a sinogram restoration method based on subspace identification. In the process of CT image reconstruction, we not only consider that there is similarity between the CT images themselves, but also between the data corresponding to the axially adjacent slices. Based on this hypothesis, we utilize a combination of subspace identification and block matching methods to improve imaging quality based on similarity. Specifically speaking, this algorithm sparsely represents the data corresponding to different slices based on subspace identification [37, 38]. In addition, to suppress the noise, block-matching frame is applied in eigen-images. To avoid the spatial information confusion and the loss of features due to the sparse representation, this algorithm sets regularization terms based on the correlation of data corresponding to different slices for constraint. The solution process of the model is carried out under the plug-and-play framework [39–41]. Compared with previous algorithms, the innovation of the proposed algorithm lies in the exploitation of both self-similarity [1, 2] and inter-slice similarity of CT images based on subspace identification. The algorithm has higher processing efficiency and better imaging quality through the combination of subspace identification method and block matching method. On the one hand, compared with the iterative method, it has higher imaging efficiency and lower computational cost. On the other hand, compared with the image post-processing, it reduces the over-smoothing situation. As the sinogram restoration method is prone to artifacts and reduced spatial resolution, the subspace identification method introduced by the proposed algorithm reduces unnecessary computation of redundant data in projections. At the same time, it introduces a correlation-based regularization term, which constrains the objective function and suppresses the artifacts in the imaging results. Experiments with simulated data and real data show that the proposed algorithm has certain advantages in image quality and imaging efficiency when the deviation of the spatial structure is within a certain range.
This paper is organized as follows. Section 2 describes the proposed algorithm. Section 3 shows the experiments with simulation data and real data to verify the algorithm. Finally, section 4 summarizes this paper, discusses some of the problems, and gives the direction of future research.
Method
Italic English letters (e.g., n
im
and N
b
) are used to indicate scalars. Matrices and vectors are denoted as bold capital letters and bold lowercase letters (e.g.,
Transmission data model and reconstruction task
In CT image reconstruction, the imaging process can be represented as
For low-dose CT image reconstruction, the main problem is to reduce the noise in the images. Considering that the structures of human organs change slowly in a single direction, there is a high degree of similarity among adjacent CT image slices and their corresponding projection data. This paper utilizes the similarity to suppress the noise. It integrates a certain number of consecutive slices and their corresponding projection data for processing, which can be written as
For normal-dose CT, given the transmission data and the geometry parameters, the image reconstruction can be simplified to find
In order to better solve the model, auxiliary variable
In the above process, the noise variance is estimated as the way in “Hyperspectral subspace identification”. The similarity has been utilized under the constraint of correlation. However, on the one hand, the computational cost of applying block-matching directly in a large amount of data is huge. On the other hand, the solution of the regularization term on correlation is complicated. Therefore, subspace identification based on different slice data is introduced to improve computing efficiency and method based on reference image guidance is introduced to help constrain the correlation among slices data.
Subspace identification based on low-rank and sparse representation
For CT images, the imaging results are self-similar and the data corresponding to adjacent slices has a high degree of similarity. Therefore, the data can be sparsely represented, which is beneficial for the exploitation of similarity to suppress the noise. In Equation (3),
According to Equation (8), the subspace transformation doesn’t change the statistical properties of the noise. To clarify the process, the noise is assumed to be additive zero-mean Gaussian and i.i.d. over all data. The noise in the transformation can be written as
After introducing the subspace identification method, there are two main problems in the solution process. For one thing, although the correlation constraints can keep the features of each slices independent, it’s necessary to find a correct reference. For another, model solving process is still complicated when the correlation regularization term is directly calculated.
Prior image constrained compressed sensing (PICCS) [45] makes up for the image quality problem caused by insufficient projection data by adding a priori image with complete data to the objective function. Based on the thought of PICCS, this paper introduces the reference image, which is used as the guidance to make the fidelity term and the correlation regularization term lower instead of gradient calculation in Equation (8). And that is called reference image constrained compressed sensing (RICCS). The main idea of RICCS is that the increase of the correlation regularization term is caused by the interaction of multi-channel information. As long as the single-channel processing and multi-channel processing are organically combined based on the objective function, the rapid decline of the objective function can be achieved. On this basis, the solution idea of Equation (8) is to first obtain the pre-processed data without the constraint of the correlation regularization term, and then solve the correlation regularization term through the reference image. The pre-processed data is called
Since the interaction of different slices’ data affects the correlation of P, this paper takes the projection processed without the interaction as reference. Referring to the specific acquisition of the reference images, that is processed channel by channel based on self-similarity, which can be written as
In Equation (11), the first term in the right hand is fidelity term. And the regularization term Φ (

Part of the processing of RICCS When calculating equation (12), Γ (·) sets a window of a certain size to select the pixels to be calculated. The window will first select the data in

The overall flow of the proposed algorithm. Before using FBP to get imaging results, the algorithm is mainly divided into the calculation of
The proposed method
In this section, all the data were processed under the same environment. All imaging procedures for transmission data in this section are performed on a personal computer. The computer contains an Intel(R) Core (TM) i7-9700 CPU@3.00Ghz. Meanwhile, the computer has two GPUs with Intel(R) UHD Graphics 630 and NVIDIA GeForce GT 1030.
Simulation experiment
This experiment tried to test the performance of the proposed algorithm under different noise conditions and summarize the characteristics of the proposed algorithm by comparing it with other algorithms. All the simulation experiments were conducted based on real clinical scanning three-dimensional data. The image size is 512 × 512, which was reconstructed from a sinogram of dimension 360 views. In the actual scanning process, the dose reduction will bring about a surge of noise. In the simulation experiments, we simulate LDCT by adding Gaussian noise to the transmission data [1]. The configuration was fan-beam. The distance from x-ray source to the rotation centre (SOD) and the detection centre (SDD) were 1000 mm and 1500 mm respectively. Considering that several CT images were uniformly processed, the proposed algorithm took the mean of structural similarity (MSSIM), the mean of peak signal-to-noise ratio (MPSNR) of these images as evaluation criteria.
In this part, the images were reconstructed by the following three algorithms for comparison. The first one was the conventional filtered back projection reconstruction algorithm (expressed as FBP), which represents the imaging effect of low-dose CT without targeted noise processing. The second one was the image reconstruction algorithm based on inexact alternating direction total-variation minimization (expressed as IADTVM) [46], which represents the iterative reconstruction(IR) method to deal with the noise in LDCT. To verify the effect of the proposed algorithm other than that of the state-of-art denoising algorithm, the last one is the Block-matching and 3D filtering (expressed as BM3D) algorithm. In the comparison algorithm, BM3D wasn’t applied in the subspace but in the normal projection domain to suppress noise. And then the FBP algorithm was applied to reconstruct the image. In this paper, the experiments were set up under i.i.d. Gaussian noise with variances of 0.04, 0.06, 0.08, 0.10, respectively. For the generality, three sets of the slices were randomly chosen to conduct the experiments with different algorithms. Every group includes 15 continuous slices up the axis and the closest slices in the seven groups were also separated by 50 phantom slices. In order to verify the effect of the algorithm on imaging objects with different structures, the same experiment was carried out with abdominal data. And except for the parameters of the algorithm, the settings in the abdominal data experiment are the same as those in the head data experiment.
In order to show the overall effect of the algorithm on a group of slices, Fig. 3 shows the imaging results of the first, eighth and fifteenth slices processed by different algorithms. In terms of horizontal comparison, the overall imaging effect of the proposed algorithm is relatively clear, and it has certain advantages in edge preservation and image contrast. At the same time, from the vertical point of view, In the results obtained by the proposed algorithm, although there are some differences in the processing effect of each slice, there is generally no slice with particularly poor processing effect. According to the principle of the algorithm, the higher the similarity between slices, the faster the processing speed. In order to measure the advantages of the proposed algorithm under the condition of high data similarity, this paper added processing time to the subsequent quantification.

Imaging results of three slices based on different algorithms. The first and second columns are the reconstruction results obtained by the FBP algorithm under normal dose and low dose conditions, which are signed by a, b, respectively. The third column is the reconstruction result obtained by the IADTVM algorithm under the condition of low dose, the fourth column is the result obtained by the BM3D algorithm, and the fifth column is the result obtained by the algorithm proposed in this paper. The above three columns are labeled c, d, e in order. And different rows show imaging results for different slices, which is labeled by 1, 2 and 3, respectively. The grayscale value displayed is [0, 0.04].
A representative group of each sets was chosen to be shown in detail, respectively. A portion of the images is selected as the region of interest (ROI), which is marked with red squares. Figure 4 shows the ground truth and the imaging results processed by different algorithms under additive Gaussian noise with variance 0.06. From the overall imaging results, compared with other algorithms, the proposed algorithm has surprising performance in noise suppression. And observing the region of interest, the clearer outline and detail is preserved by the proposed algorithm. Compared with other algorithms, the image edges are kept sharper by the proposed algorithm with less variation.

A representative slice of each sets with the region of interests There is a code in the upper left corner of each image in this figure, where different letters distinguish slices, and different numbers distinguish images processed by different methods. In order to show the overall processing effect of the algorithm on a set of slices, three slices were chosen to be shown and number 1 to 3 represent the imaging results of the first, eighth, and fifteenth slices, respectively. (a) stands for the results reconstructed by the FBP algorithm. (b) are the results reconstructed by IADTVM. (c) represents the result obtained by the BM3D algorithm. And the images marked with (d) are reconstructed by the proposed algorithm. In order to better reflect the characteristics of the algorithm, this experiment uses a square box to select the region of interest, enlarge and place it in the upper right corner of each image, using an arrow pointing at the ROI. And the grayscale displayed range of the head and abdominal images shown in Fig. 4 is [0,0.04].
Figure 5, Tables 2 and 3 show quantitative indicators of imaging results obtained by different algorithms under different noise levels. Figure 5 shows the noise suppression curves of different algorithms. As shown in the legend, different colored curves indicate different methods. Tables 2 and 3 shows the exact values of the indicators of different algorithms under different noise levels. Tables 2 and 3 represent the data of head and abdominal respectively. And these are the average of the three sets of experimental results. In the last row of each table, the paper shows the average processing time of each set of data to evaluate the computational cost of each algorithm. According to the results, compared with other algorithms, the proposed algorithm has certain advantages under the test of these three indicators, and as the noise level increases, this advantage becomes more and more obvious. Under different noise conditions, the MSSIM and MPSNR obtained by the proposed algorithm are higher than other algorithms, which indicates the image quality obtained by the proposed algorithm is better than the other algorithms. When the noise variance increases, the MSSIM and MPSNR values of various algorithms decrease to a certain extent.

Noise Rejection Curve of MSSIM, MPSNR. Two graphs show the variation curves of MSSIM, MPSNR with the increase of noise level.
Head data’s quantitative evaluation of reconstruction results based on different algorithms
Abdominal data’s quantitative evaluation of reconstruction results based on different algorithms
However, compared with other algorithms, the changes of the proposed algorithm are relatively small. This phenomenon shows that the proposed algorithm has strong adaptability in the case of large noise level, and it can be considered that the proposed algorithm can perform relatively well in the case of low dose. Referring to the processing time, the proposed algorithm keeps the computational cost in a low range while ensuring its image quality. What’s more, as shown in Tables 2 and 3, there are certain differences between different indicators. According to the process of CT data acquisition and the process of algorithmic imaging, there are two main reasons for this phenomenon. For one thing, due to the different composition of the head and abdomen, the abdomen contains more soft tissues, which leads to the different imaging results [3]. For another, the structural complexity of the head and abdomen is different. Compared with the structure of the head, the tissue structure of the abdomen is relatively complex. Therefore, overall, the imaging quality of abdominal data is lower than that of head data. The reason why the proposed algorithm is more efficient in processing abdominal data than head data is mainly due to the high similarity between abdominal data slices in the case of axial distance.
Combined with the observation of imaging results, it can be found that the proposed algorithm has unique advantages in noise suppression and detail preservation, especially in low dose conditions. In conclusion, the proposed algorithm has certain advantages over the conventional algorithm and has surprisingly good noise suppression performance with reasonable computational cost.
In order to verify the imaging effect of the proposed algorithm under the condition of less views, the paper conducts simulation experiments by adding noise to the transmission data with less views. In this experiment, the image size is 512 × 512, which is reconstructed from a sinogram of dimension 180 views. The scanning angle range is 360 degrees, and the number of projection acquisition frames is 360 frames. During the scanning process, scan once every interval. The rest of the scan settings were the same as that in the previous experiments. And the experiments were set up under i.i.d. Gaussian noise with variances of 0.06. Since previous experiments have partially verified the performance of the proposed algorithm, in this experiment, we only use the mean of PSNR as a quantitative indicator to measure the imaging effect. The data used in the experiment are the three sets of head and abdomen data used in the previous experiments and Table 4 shows the obtained imaging results. Observing Table 4, it can be found that under the same noise level, the proposed algorithm still has certain advantages when the angle is reduced (180 views).
Abdominal data’s quantitative evaluation of reconstruction results based on different algorithms
To assess the actual effect of the proposed algorithm, this experiment used real scan sinogram collected from the CT scans of the human head phantom. In the scanning system used in this paper, the model of the X-ray source is Hamamatsu L12161-07, and the model of the detector was Thales 4343F. The configuration was fan-beam. The scanning angle range was 360 degrees, and the number of projection acquisition frames was 360 frames. As shown in Fig. 6, in this experiment, the actual scanning experiment was carried out using the Chengdu Dose Phantom (CDP) [47] as the imaging object. The phantom is very similar to the human head, and the equivalent error between its synthetic material and the attenuation coefficient of human tissues and organs is not more than 5%. In this section, the real data refers to the data obtained by scanning the CDP model.

Actual data experimental phantom.
This experiment is divided into two parts. The purpose of the first part of the experiments is to quantitatively measure the effect of the proposed algorithm, so we simulate LDCT with simulated Gaussian noise and data obtained from actual scans in our experiments. The purpose of the second part of the experiment is to evaluate the performance of the algorithm in the actual imaging process. The experimental data is the data obtained by scanning under low current conditions.
In this first part, seven groups of slices were randomly selected for the experiment. Every group includes 15 continuous slices up the axis and the closest slices in the seven groups were also separated by 100 phantom slices at least. The specific scan parameters are set as follows. The distance from x-ray source to rotation axis (SOD) and to detector distance (SDD) are 92.2266 mm and 865.336 mm respectively. The tube voltage and current are 120 kVp and 200μA, respectively.
This experiment tried to test the performance of the proposed algorithm numerically and summarize the characteristics of the proposed algorithm by comparing it with other algorithms on real data. Like the simulation experiments, it used three algorithms of FBP, IADTVM and BM3D for comparison. Since the relationship between noise level and radiation dose in actual experiments is relatively complex, the experiments added noise based on normal dose images to simulate low dose, which accurately evaluated the performance of the algorithm. And the noise added in the actual experiment was independent and identically distributed Gaussian noises with variances 0.04.
This paper selected two representative groups of slices among the five groups. And a representative slice of each group was selected to show its difference images and region of interests in Fig. 7. In order to better demonstrate the characteristics of the algorithm, this paper used the image reconstructed by FBP under normal dose conditions as the reference image and selected the region of interest for comparison.

The detail of the imaging results based on different algorithms. The area of interest is marked with yellow squares. A code including letter and number is signed in the upper left corner of every picture in Fig. 7. The letter stands for the different methods of processing. In Fig. 7, a represents the reference images. b, c, d and e represent the imaging results processed by FBP, IADTVM, BM3D and the proposed algorithm, respectively. The reference images stand for the ground truth, which are mainly obtained in two steps. The CDP is first scanned under normal dose to get the sinogram. And then the ground truth is obtained by FBP based on the sinogram. The number represents the different manifestation of the characteristics of the proposed algorithm. Number 1 stands for the slice processed by different algorithms. And to better show the characteristics of the algorithm, difference images between the imaging results and the reference images are shown in the second row and marked with number 2. The enlarged region of interest is marked with number 3 and the representative details are pointed by an arrow. The grayscale values displayed in the first and third rows are [0, 0.2], and the grayscale values displayed in the second row are [–0.05, 0.05].
Figure 7 shows a representative slice’s detail of each group, called slice1 and slice2. According to the imaging results of the two groups of slices, the images’ overall clarity of the proposed algorithm is higher than that of the other algorithms, and the image is closer to the reference image. Observing difference images, compared with the imaging results of the other algorithms under low-dose conditions, the error in the image reconstructed by the proposed algorithm is noticeably reduced. At the same time, observing the contour of the original slice imaging results, the proposed algorithm’s images have no obvious difference, which proves the ability of edge preservation. Referring to the region of interest, the imaging results of the proposed algorithm have obvious edges and clear details. Overall, the image quality of the proposed algorithm is significantly higher than that of other algorithms, and it is also closer to the imaging results of normal doses, which illustrates the advantages in noise suppression and edge preservation of the proposed algorithm.
Based on the reference images in the experiment, this paper calculates the average value of the three objective evaluation indicators of PSNR, SSIM for a total of 75 test images in five groups of slices. To measure the computational cost, this paper documented the processing time of different algorithms for imaging a set of data. It can be seen from the Table 5 that the proposed algorithm has certain advantages in various indicators. The numerical performance is consistent with the visual perception, proving the algorithm has good performance in detail preservation and noise suppression. And observing the processing time, the computational cost of the proposed algorithm is also less than that of IADTVM and BM3D. Considering the imaging quality, the proposed algorithm’s computational overhead is reasonable. It is worth noting that the average processing time of real data is much smaller than that of simulated data, which is caused by the difference in similarity between slice data. In the case of low slice similarity, the computational cost of solving the correlation regularization term is relatively large.
Quantitative evaluation of reconstruction results based on different algorithms
The second part of this experiment aims to evaluate the performance of the proposed algorithm in practical applications. The experiment was set up to scan under low current conditions, and the specific parameters were set as follows. The distance from x-ray source to rotation axis (SOD) and to detector distance (SDD) are 711.555 mm and 1382.7 mm respectively. To better observe the results, we set three different current intensities (100μA, 50μA, 25μA) under the same voltage condition (120 kVp). In this experiment, three algorithms of FBP, IADTVM and BM3D were applied for comparison. The imaging results of two representative slices under different imaging conditions are presented in Fig. 8. According the results, the imaging results of FBP contains severe noise and the edge blurring in imaging results of IADTVM is relatively serious. Compared with the imaging results of BM3D, the proposed algorithm has certain advantages in image edge preservation and noise suppression. Generally, the experimental results are consistent with the previous ones, proving the potential of the proposed algorithm in the practical application of LDCT.

The imaging results based on different algorithms under 120kVp and different current As marked in the figure, columns 1, 2, and 3 correspond to the imaging results of FBP, IADTVM, and BM3D, respectively, and column 4 is the imaging results of the proposed algorithm. Each two rows in the figure correspond to a set of results, including the imaging results and the enlarged region of interest. The location of the region of interest is marked by a yellow square in the imaging results. The upper left corner of each two-line images is marked with a number. The number consists of two parts, the first part is used to represent the scanning current, and the second part is used to distinguish different slices. Different currents are marked with numbers 100, 50 and 25; different slices are marked with numbers 1 and 2. The imaging results based on different algorithms under 120kVp and different.
In order to verify the function of each part of the proposed algorithm, and also to determine the prospects of the similarity and correlation of different slices, this experiment was designed. It is divided into two parts for the regularization terms about similarity and correlation in equation (8). The two parts of the experiment were performed under the condition of simulation experiment. In this section, the number of slices in each batch is 15. All the experiments were performed under the condition of i.i.d Gaussian noise with variance 0.04.
The first part was aimed at exploring the role of the regularization terms about similarity in equation (8). The main idea of the experiment is to summarize the effect of the similarity regularization term based on the imaging results of projection data constrained only by the regularization terms about the similarity. This paper uses the FBP algorithm based on the solution result of Equation 8. And the imaging results are shown in Figs. 9 and 10.

The MSSIM, MPSNR and time variation curves with the number of subspaces Three sets were randomly chosen to conduct the experiment. And the data in it is the mean of the three sets’ results.

Imaging results corresponding to different numbers of subspaces. The region of interest is selected by a red square and enlarged in the upper right corner of each image pointed by a red arrow. To better distinguish these images, a code was signed in the upper left corner of each image. And the code includes a letter and a number. The letter distinguishes the images processed by different numbers of subspaces. For the convenience of comparison, (a) stands for the ground truth and (b) stands for the imaging results processed by FBP. (c) to (e) represents the imaging results with 1, 8 15 subspaces, respectively. And the number of the subspaces should be less than the number of dealt slices. To ensure the generality of the experimental results, except for the first column, the pictures shown in different columns were processed in different batches of slices. And the number in the code is used to distinguish different slices. Considering the impact of spatial structure differences, each group of slices is a random selection of 15 consecutive slices. And the grayscale displayed range of the image shown in Fig. 10 is [0 0.04].
Considering that the setting of the similarity regularization term is closely related to the number of subspaces, in order to better summarize its function, this paper changes the number of subspaces in the experiment to observe the changes of different imaging results. Figure 9 includes the curves of MSSIM value, MPSNR value and processing time as the number of subspaces increases from left to right for a set of slices. According to Fig. 10, as the number of subspaces increases, MSSIM and MPSNR increase. And when the two index reaches the maximum value, their values will no longer significantly change. As the number of subspaces increases, the image quality will improve until reaches the best level. However, when the subspaces increase, the processing time also increase, which means increase in computational cost. That is to say, the increase in computational overhead does not mean the same proportional increase in image quality.
Figure 10 shows the imaging results corresponding to different numbers of subspaces. Compared with the images in the first column and second column, there is little noise in the images processed based on equation 9. However, observing the spatial structures of the images, the structure in the third to fourth column is wrong and there are artifacts in these images. With the increase of the number of subspaces, the artifacts in the image can be reduced, which is confirmed with the numerical display results in Fig. 9. This shows when the number of participating slices is within a certain range, adjusting the number of subspaces can improve image quality and reduce the impact of inconsistency in spatial structure. However, from the details in the figure, the artifacts can’t be eliminated totally and the rest of them could also cause erroneous judgement.
The second part was aimed at the correlation-dependent regularization term. In this part, this paper compares the imaging results of the first part with the final imaging results and summarizes the role of the correlation regularization term. Figure 11 shows that of some representative slices. According to the displayed results, we can find that there are some artifacts in the results processed only by the similarity regularization term compared to the ground truth, which results from differences in the spatial structure of different slices. After adding a correlation regularization term, these artifacts were suppressed.

Imaging results under different processing conditions. The first row is the ground truth. The second row is the imaging results of the first part and the third row is the results. And the number of subspaces in the second row is 7. For ease of observation, this paper uses a square box to select the region of interest to zoom in in the upper right corner. The grayscale value displayed is [0, 0.04].
From the above two experiments, the similarity regularization term is mainly aimed at noise suppression and reducing unnecessary computational cost. The correlation regularization term is mainly aimed at the processing of artifacts. In the experiment, it can be verified that information interaction and low-rank representation of multiple slices can bring surprisingly good noise suppression effect while it could bring wrong information because of the difference among these slices. This experiment also demonstrates that the wrong information can be suppressed by the correlation regularization term.
In this paper, a projection domain processing method for low-dose CT reconstruction based on subspace identification is introduced. There are two main reasons for this work including:1) The methods based on sparse representation perform well in the natural image, which leads to the attention on the application in CT images; 2) The self-similarity of CT images and the similarity between different slices is worthy to be exploited. In this paper, the following work has been done for the aspects mentioned above. First, this paper presents a method based on subspace identification to sparsely represent the CT images, which obtain great effect. Second, batch-matching is adopted to utilize the self-similarity and the similarity between different slices. And third, a regularization term based on correlation has been proposed to avoid the confusion and loss of feature. In general, this algorithm is applied to the simulated head data, abdominal data and actual scanned head model data, which demonstrate its improvement for LDCT image quality.
Compared with the previous algorithms, the proposed algorithm mainly has the following advantages. First, the algorithm does not require a large amount of data for training, and it is more interpretable than DL. Second, as a sinogram restoration algorithm, it estimates noise more directly than other algorithms, avoiding over-smoothing. In addition, the algorithm uses the subspace method to compress the data, avoids the computational overhead of redundant data, and has high imaging efficiency. Compared with the previous sinogram recovery algorithm, reducing the processing of projection data is also beneficial to alleviate the problem of sensitivity of sinogram. In response to this problem, the proposed algorithm also introduces a correlation-based regularization term, which has a certain suppression effect on the artifacts that appear in the processing process. In general, the proposed algorithm does not require a large amount of data for training, and it can consider both imaging efficiency and imaging quality, and thus it has promising practical application prospects.
Despite the improvement in the mean of PSNR and SSIM, the proposed algorithm doesn’t perform well for every single image. Due to the difference in spatial structures, the similarity exploited by this method is only for most slices, i.e., for a certain slice, if the similarity of the slices processed together is low, the effect of the algorithm is not ideal, which causes the reduction in efficiency. Additionally, the generated images still contain some artifacts due to the sensitivity of sinogram. For the future work, on the one hand, we will consider choosing the processed pictures to improve processing efficiency. Since the existing metrics are not ideal for measuring the similarity between noisy images such as PSNR, SSIM and RMSE, we consider classifying the processed images by simply processing the noisy images. On the other hand, we will also consider improving the way the proposed algorithm search for the similar blocks. Since the effect of artifacts removal is not totally, we consider optimizing the way we obtain and utilize the prior information. And we plan to summarize the reasons why the proposed method produces artifacts in the denoising process according to the characteristics of sinograms.
Footnotes
Acknowledgments
This work was also supported by the National Key Research and Development Project of China (Grant No. 2020YFC1522002). This work was supported by the National Natural Science Foundation of China (Grant No. 62101596) and the China Postdoctoral Science Foundation (Grant No. 2019M663996).
