Low-dose CT reconstruction method based on prior information of normal-dose image

Abstract

BACKGROUND:

Radiation risk from computed tomography (CT) is always an issue for patients, especially those in clinical conditions in which repeated CT scanning is required. For patients undergoing repeated CT scanning, a low-dose protocol, such as sparse scanning, is often used, and consequently, an advanced reconstruction algorithm is also needed.

OBJECTIVE:

To develop a novel algorithm used for sparse-view CT reconstruction associated with the prior image.

METHODS:

A low-dose CT reconstruction method based on prior information of normal-dose image (PI-NDI) involving a transformed model for attenuation coefficients of the object to be reconstructed and prior information application in the forward-projection process was used to reconstruct CT images from sparse-view projection data. A digital extended cardiac-torso (XCAT) ventral phantom and a diagnostic head phantom were employed to evaluate the performance of the proposed PI-NDI method. The root-mean-square error (RMSE), peak signal-to-noise ratio (PSNR) and mean percent absolute error (MPAE) of the reconstructed images were measured for quantitative evaluation of the proposed PI-NDI method.

RESULTS:

The reconstructed images with sparse-view projection data via the proposed PI-NDI method have higher quality by visual inspection than that via the compared methods. In terms of quantitative evaluations, the RMSE measured on the images reconstructed by the PI-NDI method with sparse projection data is comparable to that by MLEM-TV, PWLS-TV and PWLS-PICCS with fully sampled projection data. When the projection data are very sparse, images reconstructed by the PI-NDI method have higher PSNR values and lower MPAE values than those from the compared algorithms.

CONCLUSIONS:

This study presents a new low-dose CT reconstruction method based on prior information of normal-dose image (PI-NDI) for sparse-view CT image reconstruction. The experimental results validate that the new method has superior performance over other state-of-art methods.

Keywords

Computed tomography sparse-view prior image prior matrix image reconstruction

1 Introduction

Radiation risk brought by computed tomography (CT) is an increasing concern for patients and researchers, especially in cases where repeated CT scanning is needed because of dose cumulation [1]. Sparse-view CT scanning, implemented by setting a sparser scanning angular interval, is an effective protocol for radiation dose reduction and protecting the scanned patient from radiation injury [2, 3]. Besides radiational dose reduction, sparse sampling of image data acquisition results in higher signal-to-noise (SNR) of the projection data since the influence of electronic readout noise can be suppressed by that the radiational exposure remains high for every individual projection. However, using traditional analytical algorithms such as filtered back-projection (FBP) for sparse-view projection data reconstruction may lead to severe streak artifacts caused by data insufficiency in the resulting images. Many studies aimed at advanced algorithms for sparse-view CT reconstruction, including many novel prior image aided model-based optimization algorithms and all kinds of methods involving deep learning models for artifact elimination and denoising [4 –8]. For sufficient reduction of radiation dose, the number of projections per scan needs to be reduced further, and stronger algorithms are needed. From this aspect, the regularization methods in the mathematical model for iterative reconstruction and the usage of prior information play important roles.

Prior image associated with total variation (TV) regularization has proven to be effective for sparse-view projection data reconstruction for CT imaging [4 , 9–11]. In many CT imaging studies involving repeated X-ray scanning, such as the conservative management of asymptomatic epidural hematomas [12], the regular monitoring of the patient with cystic fibrosis [13] or the geometric registration of an image-guided surgery system [14], prior information from the previous reconstructed images can be very helpful for the current reconstruction task. For a particular patient under CT inspection, a sparse view low dose protocol can be considered if a previously obtained CT image of the same body region can be accessed from his or her medical records and the prior information from which can be introduced to the current imaging task. The prior image aided reconstruction possesses significant application prospect though there are no clinical CT imaging applications making use of previously reconstructed images. Researchers have investigated some prior image-aided algorithms, especially for this issue. Guanghong Chen et al. applied the prior image to constrain the objective image during the iterations of the CS-based algorithms in the sparse field (PICCS), from which the CT image method can be reconstructed from very sparse projections with diagnostically acceptable quality [4]. The number of projections can be further reduced with the nonconvex PICCS (NCPICCS) algorithm proposed by McCollough et al., in which a more general nonconvex sparsity measure is employed [7]. Jianhua Ma et al. proposed the PWLS-ndiTV algorithm to resolve the problem of misalignment in the PICCS method by updating the prior image according to the similarity between the estimated image and the normal-dose pre-scanning image [8]. With unregistered prior images, Pourmorteza et al. used 2D-to-3D registration to obtain a transformed prior image, obtaining the difference image between the origin prior image volume and the current anatomy [15]. All these newly proposed methods work through solution optimization in the image space. However, we note that all the algorithms mentioned above involve the usage of prior information explicitly in the image domain, which means applying the regularization or refinement directly on the image that under reconstruction, which may somewhat weaken the mapping from the acquired raw data to the reconstructed images. A different concept, involving the prior information implicitly in the process of forward projection, can be imagined by contrast to avoid direct modification on the reconstructed image data.

In this study, we propose a low-dose CT reconstruction method based on prior information of normal-dose image (PI-NDI). In our algorithm, a concept borrowed from machine learning, the kernel function, is involved and a normal-dose prior image, which represents the previous scanning, is employed for the sparse view CT reconstruction. The theory of kernel, which represents a variable conversion from its original data space to a transformed space, was originally designed for geometry nonlinearization [16]. In our algorithm, the kernel function is used for modeling the attenuation coefficients of the target tissue as a function of their corresponding feature and building a prior matrix that is involved in the forward-projection process in the reconstruction iterations, through which the prior information, namely, the pixel values of the prior image, is employed implicitly. Guobao Wang and Jinyi Qi have investigated this kernelized method for low count positron emission tomography (PET) and prove the effectiveness [17, 18]. The weakening of the data mapping from projection data to reconstructed image with direct modification of image intensities and limitation of projection number reduction are the main two limitations of the previously proposed methods for sparse CT reconstruction, however, these issues can be recognized as the two main advantages of the proposed PI-NDI method in this study since it introduce stronger prior assistance without directly weakening the data mapping. The performance of the proposed PI-NDI method is evaluated on a digital extended cardiac-torso (XCAT) phantom and a physical head phantom by measuring the quantitative indexes in terms of the reconstruction accuracy, image noise and contrast resolution for CT image evaluation.

The remainder of the paper is organized as follows. Section 2 describes the mathematical approach of the proposed PI-NDI method in detail, including the basic EM framework for image reconstruction and the method of prior information-based (PI) reconstruction. Section 3 describes the simulative and real image experiments in this work, and section 4 gives the results and some discussion on the proposed PI-NDI method. Finally, the conclusion of our study is given in section 5.

2 Materials and methods

2.1 Problem formulation

2.1.1 CT model

Without loss of generality, an X-ray CT system with a monoenergetic beam can be approximately described with a set of linear equations, which can be written in matrix form as:

$y = Ax$ (1)

where y is the projection data after system calibration and logarithm transformation, i.e., y = (y₁, y₂, …, y_M) ^T, and x represents the attenuation coefficients to be reconstructed, i.e., x = (x₁, x₂, …, x_M) ^T. ‘T’ represents matrix transposition, and M and N are the numbers of pixels of the projection data and the image, respectively. A is the system matrix corresponding to the imaging system under investigation with the specified imaging geometry conditions, which has the size of M × N. An element of A , A_i,j, describes how much the j th pixel of the image contributes to the integration of the i th ray. The objective of CT image reconstruction is to find an approximate solution of x with the measured projection data y . In our work, the system matrix A is predetermined and stored for convenient usage. In the task of CT image reconstruction from sparsely sampled projection data, the system matrix contents the information of the condition of sampling with larger angular intervals. Compared to normal dose imaging protocol, this system matrix A brings projection data with smaller size.

2.1.2 Expectation-maximization method for CT image reconstruction

In this study, the proposed PI-NDI approach is based on the maximum likelihood expectation maximization (MLEM) method [19 –23]. The MLEM method for CT image reconstruction, which explicitly include a Poisson model [19 , 22], shows

$x^{n + 1} = \frac{x^{n}}{A^{T} 1_{M}} \cdot (A^{T} \frac{y}{A x^{n}}),$ (2)

where n denotes the iteration index and 1 _M is a column vector with M elements of value 1. We note that A ^T represents the back-projection operation in equation (2).

The proposed PI-NDI method is based on the MLEM algorithm given above, for which we apply an prior matrix for the image to be reconstructed in the MLEM loop as a mask, and the information from the normal-dose prior image, alone with an appropriate kernel function, is introduced to build the prior matrix.

2.2 Proposed algorithm

2.2.1 High-dimensional image representation with kernel function

Inspired by the novel application of kernel functions, which use nonlinear operations on variables in low-dimensional space to represent the point product of variables in high-dimensional space, in machine learning [16, 18], a novel representation of CT images and reconstruction algorithms based on the application of kernel function are proposed in this study.

We can remodel a CT image by defining a mapping Γ from the feature f _j, which can be a tensor of any order, to the corresponding pixel value x_j:

$x_{j} = Γ (f_{j}) .$ (3)

Γ can be simply modeled as a linear function about the corresponding variable of f _j in the transformed space, θ ( f _j ):

$Γ (f_{j}) = w^{T} θ (f_{j}) .$ (4)

We note that the form of θ ( f _j ) is implicitly determined by a kernel function. And the weight vector, w , is also determined by all the transformed variables and a set of coefficients ξ as

$w = \sum_{k = 1}^{N} ξ_{k} θ (f_{k}) .$ (5)

Combining (3), (4) and (5) we can obtain a representation of the image pixel value by the space-transformed feature vectors θ ( f ) and the coefficient vector ξ :

$x_{j} = \sum_{k = 1}^{N} ξ_{k} θ {(f_{k})}^{T} θ (f_{j}) .$ (6)

The kernel function of f _k and f _j, ψ ( f _k , f _j ), is introduced here to replace the inner product of θ ( f _k) and θ ( f _j ) in (8), and the pixel value of the CT image is modeled by

$x_{j} = \sum_{k = 1}^{N} ξ_{k} ψ (f_{k}, f_{j}),$ (7)

or in matrix form by

$x = Ψ ξ .$ (8)

The expression of x in (8) shows that the image to be reconstructed is now represented by a masked coefficient image where ψ is the designed prior matrix. We note that in the transformed space, the variables θ is implicitly determined by the chosen kernel function ψ , and the exact mapping from f to θ is not required in our proposed method.

2.2.2 Prior information (PI) aided CT reconstruction

Inserting (8) into (1), we obtain the CT imaging model represented by the coefficient image ξ :

$y = A Ψ ξ .$ (9)

With a predetermined prior matrix Ψ , AΨ is the equivalent system matrix in the task of the reconstruction of the coefficient image ξ , and the transformed EM iteration rule becomes

$ξ^{n + 1} = \frac{ξ^{n}}{Ψ^{T} A^{T} 1_{M}} \cdot (Ψ^{T} A^{T} \frac{y}{A Ψ ξ^{n}}) .$ (10)

Once the estimated coefficient image $\tilde{ξ}$ is found, the reconstructed image $\tilde{x}$ is calculated by

$\tilde{x} = Ψ \tilde{ξ} .$ (11)

Different from other algorithms involving prior images reconstructed from previous scans, such as PWLS-PICCS and PWLS-ndiTV, which use prior information in the process of regularization directly in the image field [4, 8], the proposed PI-NDI method introduces the prior knowledge in the forward-projection model, which improves the performance of image reconstruction from sparse-view CT projection data.

2.2.3 Construction of the prior matrix

1) Selection of neighbors

As discussed in section 2.B.1, the element ψ ( f _j, f _k ) of the prior matrix Ψ is calculated based on the feature vectors f _k and f _j of their respective corresponding image pixels x_k and x_j. However, if the connection between every two pixels of the image are considered, the resulting prior matrix will be a full matrix with size N² × N². This will bring great computational cost to every reconstruction iteration and make the algorithm unfeasible. One approach to make ψ sparse is to apply the k-nearest neighbors (kNN) [24] principle to determine the n selected feature vectors for the calculation of the element of the prior matrix ψ_j,k based on f _j and f _k:

$ψ_{j, k} = {\begin{matrix} ψ (f_{j}, f_{k}), & f_{k} \in kNN cluster of f_{j} \\ 0, & otherwise \end{matrix},$ (12)

where the kNN cluster is determined based on the Euclidean distance between f _j and f _k. In this study, to involve the prior information integrally from the normal-dose previous CT image, the knearest neighbors are searched from the entire image space, which represents a classical non-local mean denoising method [25]. The number of selected neighbors, k, is also an important parameter for the PI-NDI method, which is demonstrated by a trade-off between computational cost and the extent of prior information utilization. In our work, k is set to 96, meaning that 96 neighbors are searched for every image pixel to build the prior matrix Ψ , and the discussion on how this parameter is determined is provided in section 4.C.

2) Formulation of the feature vectors

The feature tensors, which directly contain prior knowledge, significantly affect the effect and robustness of the proposed PI-NDI method. In our work, the feature tensor f _j of pixel x_j is designed as a vector consisting of the image intensities of the prior image. A pivotal parameter that needs to be discussed is the patch size n, which determines the extent to which the prior information from the spatial neighbors of the corresponding pixel in the prior image are adopted in f _j. In the case that the aiding image has a larger size than the objective image, a patch instead of a single pixel in the prior image corresponds to a pixel in the reconstructed image, and n > 1 is needed. However, in our task of CT imaging aided by a previous CT image, the prior image involved is of the same size as the object image, in which case using n > 1 may introduce some error information since the intensity of the spatial neighbors of pixel j in the prior image may not be representatively mapped to x_j via the mapping function Γ defined in (3). Thus, n = 1 is used in our work, and the feature vector contains only the corresponding pixel intensity of the prior image, namely,

$f_{j} = x_{p, j} .$ (13)

Here, f _j is then normalized with its standard deviation for the easy determination of the standard deviation σ in formula (15) for kernel function calculation by

$\tilde{f_{j}} = \frac{f_{j}}{σ_{f}},$ (14)

where σ_f is the standard deviation of f _j over all pixels.

3) Calculation of the prior matrix

The definition of the kernel function ψ ( f _j, f _k ) relating f _j and f _k represents how their inner product in the transformed space is like in the original data space. Typically, two types of kernel function can be considered in this work, one of which is the Gaussian function:

$ψ (f_{j}, f_{k}) = exp (- \frac{{| f_{j} - f_{k} |}_{2}^{2}}{2 σ^{2}}),$ (15)

where σ is the standard deviation of the Gaussian distribution that was designed previously. In this case, the mapping function θ ( f ) implicitly defined by $θ {(f_{k})}^{T} θ (f_{j}) = exp (- \frac{{| f_{j} - f_{k} |}_{2}^{2}}{2 σ^{2}})$ is of infinite dimension [16]. Another type of kernel function widely used in machine learning is the polynomial kernel, which shows

$ψ (f_{j}, f_{k}) = {({f_{j}}^{T} f_{k} + C)}^{p},$ (16)

where C and p are constant parameters.

In our study, both of the above two methods of prior matrix calculations are tested, and the Gaussian type, which involves the measurement of the distance between the two considered feature vectors, shows superior performance and is selected for the imaging tasks using the PI-NDI method.

2.3 Compared algorithms

In this study, the proposed PI-NDI method is an EM-based iterative algorithm involving a previous prior image from repeated scanning. Therefore, the reconstruction results from the traditional MLEM algorithm with TV regularization, namely, MLEM-TV; the PWLS-TV, which is also applied for prior image reconstruction; and the classical PWLS-PICCS algorithm proposed by Chen [4], are shown for comparison. These algorithms chosen for comparison are representing different type of iterative methods for CT image reconstruction. The PICCS method is an algorithm that explicitly applies the information from the prior image, since the deviation between the target image and the prior image is directly extracted and involved in the total variation regularization. On the contrary, the proposed PI-NDI represents an algorithm that implicitly apply the prior information. The MLEM-TV and PWLS-TV are iterative algorithms regularized by total variation minimization without prior image involved, of which the results are provided to see how much the proposed PI-NDI method performs better with the creatively introduced prior image.

3 Experiments and evaluation

3.1 Data acquisition

The validation and evaluation of the proposed PI-NDI method are conducted with a digital NURBS-based XCAT phantom and a physical diagnostic head phantom. In our study, we assume the strict registration between the normal-dose prior image and the current object image, and the discussion about misalignment and prior image error is given in section 4.D.

3.1.1 Digital XCAT phantom

Fig. 1a shows a digital ventral phantom used for the static CT simulation. Two regions of interest (ROIs) are selected for evaluating the performance of the reconstruction technique from the aspect of detail restoration. The phantom is composed of 512×512 pixels with a size of 0.625 mm×0.625 mm.

Fig. 1

a) Digital XCAT ventral phantom; three ROIs and a profile line are selected for performance evaluations and sectional profile display. b) Reconstructed image of the diagnostic head phantom via PWLS-TV. An ROI is selected for small object display.

3.1.2 Data acquisition from simulation

The geometry for the simulation experiments in our work is designed to represent a monoenergetic fan-beam CT scanner that acquires 360 views of projections over 360° with a circular gantry orbit. The number of detector channels is 1000, and the total length of the detector of 625 mm, which gives a size of 1000×360 for the projection data from one scanning. The distances from the X-ray source to the detector and to the origin of rotation (the SDD and SOD, respectively) are 949 mm and 541 mm, respectively. Each projection datum along an X-ray through the sectional image is calculated based on the known attenuation coefficients and intersection areas of the X-ray with the geometric shapes of the objects in the sectional image. The method for obtaining the noisy projection data is similar to that in the studies of [26 –28]. The noisy measurement of a detector bin b_i is generated based on the noise-free measurement y_i calculated from model (1), according to the statistical model of pre-logarithm projection data:

$b_{i} = Poisson (I_{0} \exp (- y_{i})) + Normal (0, σ_{e}^{2}),$ (17)

where I₀ is the incident X-ray intensity and $σ_{e}^{2}$ is the background electronic noise variance. In this work, I₀ is set to 5×10⁵ and $σ_{e}^{2}$ is set to 10 to represent the normal-dose scanning protocol. After adding the simulation noise, the noisy measurement y_i used for reconstruction is acquired by preforming a logarithm transformation on b_i. The sparse-view projection data are obtained by under-sampling the whole set of 360 views of noisy projection to 8, 15, 30, 60, 120 views evenly over 360° to represent the current low-dose scanning.

3.1.3 Diagnostic head phantom

A diagnostic head phantom (Atom Max 711 HN, CIRS Inc., VA, USA) was used for the physical experiment in this work. CT scanning was performed on our laboratory’s CT platform (Varex G-242, Varex Imaging Corporation, UT, USA) with the X-ray tube voltage and the tube current set to 120 kV and 11 mA, respectively. During the experiments, a bowtie and additional 0.5 mm copper filtration were used. An energy-resolving photon counting detector (XC-Hydra FX50, XCounter AB, Sweden) with an imaging area of 512 mm×6 mm and a native component size of 0.1 mm×0.1 mm was used. The projection data were rebinned from 5200×60 pixels to 850×10 pixels with the following two steps:

Abandon the 10 columns of data collected by both the left and right edges of the detector panel;

Rebin the dataset with a 6×6 binning mode.

The middle two layers of the rebinned 3D projection data were extracted and averaged to obtain the 2D data. The SDD was 1500 mm, and the SOD was 1000 mm. The gantry was rotated by 360° with a 1° angular projection interval, and the sparse-view data were obtained by under-sampling the whole set of 360 views of projections to 8, 15, 30, 60, 120 views evenly over 360° to represent the low-dose scanning. Figure 1b shows the image reconstructed from the normal-dose projection data with the PWLS-TV algorithm (α = 0.03, 200 iterations), which is used as the prior image in this study.

3.2 Performance evaluation

As marked in Fig. 1a, three ROIs are selected for quantitative evaluation of the proposed PI-NDI method in our simulation experiment. ROI-1 is complex, while ROI-2 and ROI-3 are uniform regions in the digital phantom. For the physical experiment, an ROI is selected for small object display.

3.2.1 Evaluation by visual inspection

All reconstructed images from both the simulation experiment and the physical experiment are shown in this paper, and the selected ROI-1s are displayed for visual evaluation in detail.

3.2.2 Evaluation of the reconstruction accuracy

The root-mean-square-error (RMSE) is calculated between the true phantoms and the reconstructed images to quantify the performance of the algorithms from the aspect of imaging accuracy. The RMSE is defined as

$RMSE (x^{recon}) = \sqrt{\frac{\sum_{j = 1}^{N} {(x_{j}^{recon} - x_{j}^{true})}^{2}}{N}},$ (18)

where x ^recon denotes the reconstructed image and x ^true denotes the ground truth phantom. N is the total number of pixels in the image or a selected ROI. For the simulation experiment, the RMSE is measured on the whole image and ROI-1.

3.2.3 Evaluation of the noise reduction

For the evaluation of performance of the proposed method and the compared methods on noise reduction, two typical quantitative metrics, including the peak signal-to-noise ratio (PSNR) and mean percent absolute error (MPAE), are used, calculated by [8, 26]:

$PSNR (x^{recon}) = 10 \times {log}_{10} (\frac{MAX {(x^{true})}^{2}}{\sum_{j = 1}^{N} {(x_{j}^{recon} - x_{j}^{true})}^{2} / N - 1}),$ (19)

and

$MPAE (x^{recon}) = \frac{100}{N} \sum_{j = 1}^{N} | \frac{x_{j}^{recon}}{x_{j}^{true}} - 1 |,$ (20)

where x ^recon, x ^true and N have the same definition as in (20). For the simulation experiment, the PSNR and MPAE are measured on the whole image and the three selected ROIs.

3.3 Implementation

The flowchart shown in Fig. 2 represents the procedure of our proposed PI-NDI method for CT image reconstruction. From the clinical aspect, the low-dose projection data come from sparse scanning. However, the sparse-view projection data are obtained from under-sampling of the fully sampled data in our work for convenience. The approach is composed essentially of three steps:

Fig. 2

The flowchart of the PI-NDI method.

Reconstruction of the prior image x _p with the fully sampled projection data using the PWLS-TV algorithm, which has the cost function

$L (x_{p}) = {(y - A x_{p})}^{T} Σ^{- 1} (y - A x_{p}) + α TV (x_{p})$ (21)

and an update rule

${x_{p}}^{n + 1} = x_{p} - β^{n} (A^{T} (\sum^{- 1} (y - A {x_{p}}^{n}))) - α \nabla TV {(x_{p})}^{n} .$ (22)

The prior image from the above step represents the image reconstructed from the pre-scanned normal-dose projection data. In formula (22), ∇ TV ( x _p) ⁿ is the gradient of TV ( x _p) after the nth updating, and βⁿ is the step size for the fidelity maintaining for the (n + 1) th step, which is adaptively calculated by [27]: $β^{n} = \frac{{(A^{T} \sum^{- 1} (y - A {x_{p}}^{n}))}^{T} (A^{T} \sum^{- 1} (y - A {x_{p}}^{n}))}{{(A A^{T} \sum^{- 1} (y - A {x_{p}}^{n}))}^{T} \sum^{- 1} (A A^{T} \sum^{- 1} (y - A {x_{p}}^{n}))} .$ (23)

Construction of the image prior matrix Ψ with the feature vectors extracted from x _p through the process described in section 2.B.3.

Reconstruction of the coefficient image ξ with the sparse-view projection data via (10) with the prior matrix Ψ and restoration of the object image x from ξ via (11).

In summary, the pseudo-code for our simulative experiments using the PI-NDI method can be described in the following chart:

The determinations of the parameters used in the proposed and compared algorithms are based on experiments. The number of iterations Iter is set to 200, which is enough for all the algorithms to converge for all the iterative methods, including the MLEM-TV and PWLS-PICCS for comparison, PWLS-TV for prior image reconstruction and PI method for sparse-data reconstruction. The TV step size α is selected to be 0.003 and 0.03 for the digital simulation experiments and for the diagnostic head phantom study, respectively. In particular, for the proposed PI-NDI method, a patch size n = 1 is used, which representing a pixel-based feature extraction as described in section 2.B.3. The number of neighbors k, which determines to what extent prior information is acquired from the pre-scanned image, is set to 96 for an acceptable trade-off between the image quality and computational cost.

The experiments in this study were conducted on MATLAB R2019b, MathWorks Inc. We employed the Image Reconstruction Toolbox (IRT) (Copyright Jeff Fessler, University of Michigan) to calculate the imaging system matrix G and conduct the simulative forward projection and the commercial FBP reconstructions.

4 Results and discussion

4.1 XCAT phantom study

4.1.1 Visual inspection

In the simulation experiments of our works, images are reconstructed with different numbers of projections, namely, 8, 15, 30, 60, 120 and 360 views, to evaluate the performance of the PI-NDI method. The influence of the number of projection views shows a difference in the proposed algorithm compared to the other algorithms, and the proposed PI-NDI method performs superiorly in the image reconstructions from ultrasparse projection data. Figure 3 shows the results from different reconstruction algorithms with variable numbers of projections. In the case of 8 views, 15 views and 30 views, only the PI-NDI method gives acceptable results, and sparse artifacts do not exist in the reconstructed images. The results become increasingly close to the ground truth as the number of projections used for reconstruction increases for all the methods compared (for the cases with 120 and 360 views, not shown here). The PWLS-PICCS method, also assisted by the prior image, give a result that is visually comparable to the PI-NDI method in the case with 60 views of projection data. Figure 6 shows the ROI-1 of the images reconstructed via MLEM-TV, PWLS-TV, PWLS-PICCS and the proposed method with 8 and 60 views of projections. The superior capability for sparse-data reconstruction of the PI-NDI method is well demonstrated since its results are closest to the ground truth.

Fig. 3

Results of the simulative experiments from MLEM-TV, PWLS-TV, PWLS-PICCS and the proposed PI-NDI method with different numbers of projections. Window of displaying is (0, 0.0139).

4.1.2 Reconstruction accuracy

Figure 5 shows the RMSE curve of the whole image and the selected ROI-1 in the reconstructed images with respect to the number of views. The advanced performance of the kernel method in CT imaging with sparse projection data is clear, and the results given by the PI-NDI method with 15 and 30 projections are of comparable quality to those given by other algorithms with fully sampled data.

Fig. 4

Residual images of the simulative experiments. The numbers on the left indicate the number of projections used for reconstruction. Window of displaying is (–0.002, 0.002).

Fig. 5

RMSE curve of the reconstructed images with different numbers of projections from the four algorithms. The RMSE is measured on the whole image and ROI-1 between the reconstructed images and the ground truth. a) The RMSE measured on the whole image; b) The RMSE measured on ROI-1.

Fig. 6

The ROI-1 in the reconstructed images of the digital XCAT phantom using 8 and 60 projections. The numbers on the left indicate the number of projections used for reconstruction. ROI-1 is marked in Fig. 1a. Display window (0, 0.0139).

Figure 7a∼b shows the sectional profile of the reconstructed images from 60 and 120 views of projections, and the pixels profiled are marked by the red line in Fig. 1a. The profiles given by the PI-NDI method are matched remarkably better with the ground truth even though all the algorithms give acceptable results with 120 projections, as shown in Fig. 3. From another aspect, the edge response of the proposed algorithm is sharper than that of the compared algorithms, which benefits the reservation of spatial resolution. According to Fig. 7c, the profiled images reconstructed by the PI-NDI methods with different numbers of projections show very small differences from each other, which indicates that the robust performance in sparse-data reconstruction of the proposed PI-NDI method.

Fig. 7

Sectional profiles of the reconstructed images. The line of profile is marked in red in Fig. 1a. a) Sectional profiles of images reconstructed via different algorithms with 60 views; b) Sectional profiles of images reconstructed via different algorithms with 120 views; c) Sectional profiles of images reconstructed via PI-NDI method with 8, 15, 30, 60 views.

4.1.3 Noise reduction measure

Table 2 shows the MPAE values calculated from the whole image and the three selected ROIs of the reconstructed images using different numbers of views of projection data. We mark the best value of every ROI with every sparse protocol. It is remarkable that for the very sparse cases, namely, from 8 views to 30 views, the PI-NDI method provides lower MPAE values compared to other algorithms, which prove the superior performance of the proposed method in sparse CT reconstruction. However, as the number of views increases, all the algorithms discussed seem giving comparable results, even PICCS and PWLS-TV methods can provide better quantitative results with 120 views and 360 views, respectively. The similar tendency with projection views increasing can be observed in the results of PSNR values shown in Table 3. In the cases with few views, the PI-NDI performs best in noise reduction. With more than 60 views of projections, the advantages of the PI-NDI method in sparse CT reconstruction is diminished.

Table 1
The pseudocode for the PI-NDI method

Algorithm Prior information of normal-dose image aided low-dose CT recon. (PI-NDI)

Input: Normal-dose sinogram y , Low-dose sparse sinogram y_s, System matrix A, A _s

Output: Reconstructed image x

Parameters: No. of iterations Iter , TV step size α , Patch size n, No. of neighbors k

1. x _p = PWLS-TV( y , A , α, Iter ); // Prior image reconstruction via PWLS-TV

2. Ψ = PriorMatrix-build( x _p, n , k ); // Build the prior matrix

3. x ₀ = 1, ξ = Ψ ^T x ₀; // Initialization

4. for i = 1: Iter do

ξ ⁿ⁺¹ = PI - update ( ξ ⁿ, Ψ , y _s, A _s) ;

endfor

5. x = Ψ ξ ^Iter ;

6. return x

Algorithm Prior information of normal-dose image aided low-dose CT recon. (PI-NDI)
Input: Normal-dose sinogram y , Low-dose sparse sinogram y_s, System matrix A, A _s
Output: Reconstructed image x
Parameters: No. of iterations Iter , TV step size α , Patch size n, No. of neighbors k
1. x _p = PWLS-TV( y , A , α, Iter );	// Prior image reconstruction via PWLS-TV
2. Ψ = PriorMatrix-build( x _p, n , k );	// Build the prior matrix
3. x ₀ = 1, ξ = Ψ ^T x ₀;	// Initialization
4. for i = 1: Iter do
ξ ⁿ⁺¹ = PI - update ( ξ ⁿ, Ψ , y _s, A _s) ;
endfor
5. x = Ψ ξ ^Iter ;
6. return x

Table 2

MPAEs measured on the whole reconstructed image and the selected ROIs (marker in Fig. 1a)

No. of views	Algorithms
		MLEM-TV	PWLS-TV	PWLS-PICCS	PI-NDI
8	whole	11.27	13.33	11.89	2.27
	ROI-1	10.08	7.74	6.36	2.75
	ROI-2	3.98	5.02	4.70	1.64
	ROI-3	78.33	128.05	125.66	5.95
15	whole	8.24	10.30	8.56	1.92
	ROI-1	7.04	7.12	5.04	2.24
	ROI-2	1.99	2.80	2.63	1.22
	ROI-3	33.55	97.06	95.00	3.45
30	whole	5.56	6.63	3.49	1.57
	ROI-1	5.93	6.61	3.42	1.95
	ROI-2	1.64	0.87	0.54	0.83
	ROI-3	11.02	12.44	9.38	3.48
60	whole	3.66	3.72	1.56	1.47
	ROI-1	4.36	4.74	2.32	1.88
	ROI-2	1.35	0.68	0.28	0.69
	ROI-3	4.33	8.84	3.02	3.30
120	whole	2.17	1.60	1.00	1.43
	ROI-1	2.72	2.51	1.70	1.84
	ROI-2	0.23	0.32	0.20	0.62
	ROI-3	1.84	1.43	0.55	3.11
360	whole	1.70	0.92	0.96	1.41
	ROI-1	2.15	1.48	1.62	1.82
	ROI-2	0.12	0.19	0.18	0.58
	ROI-3	2.13	0.39	0.42	3.24

Table 3

PSNRs measured on the whole reconstructed image and the selected ROIs (marked in Fig. 1a)

No. of views	Algorithms
		MLEM-TV	PWLS-TV	PWLS-PICCS	PI-NDI
8	whole	18.73	16.45	17.50	30.97
	ROI-1	19.68	21.23	23.13	29.84
	ROI-2	26.76	25.05	25.82	33.67
	ROI-3	2.05	– 2.18	– 2.01	22.57
15	whole	21.49	18.72	20.61	31.87
	ROI-1	21.80	21.98	24.93	31.25
	ROI-2	32.71	29.60	30.46	36.08
	ROI-3	8.88	0.22	0.42	27.27
30	whole	24.13	21.20	26.68	32.74
	ROI-1	23.61	22.91	28.11	32.04
	ROI-2	34.36	39.43	43.30	39.52
	ROI-3	16.87	15.69	18.64	27.56
60	whole	26.84	25.60	31.25	32.91
	ROI-1	25.95	25.48	30.80	31.90
	ROI-2	35.52	39.85	49.53	41.11
	ROI-3	25.99	18.75	28.45	28.17
120	whole	29.11	30.22	32.18	32.99
	ROI-1	28.87	29.53	31.94	31.99
	ROI-2	49.87	45.32	51.50	42.09
	ROI-3	34.24	34.78	43.31	28.51
360	whole	29.85	32.30	32.23	33.02
	ROI-1	30.22	32.47	32.08	32.07
	ROI-2	56.50	51.26	52.09	42.63
	ROI-3	32.93	46.33	45.59	28.28

4.2 Diagnostic head phantom study

Figure 8 shows the image reconstructed by MLEM-TV, PWLS-TV, PWLS-PICCS and PI-NDI method. For cases with sparse projection data, only the proposed method can restore an acceptable image. For the PI-NDI method, the difference between the results from different numbers of projections is reflected mainly in the level of the image noise. The remarkable distinction in smoothness is shown in Fig. 9, which shows the selected ROI of the reconstructed images. When 360 views of projection data, i.e., the fully sampled data defined in our work, are employed for reconstruction, the advantage of the proposed PI-NDI method seemed unremarkable, as the results given by all the involved algorithms are visually legible enough and comparable to each other. The advantage of the PI-NDI method is revealed if one notes that the image information contented in a sparse-data reconstructed image via PI-NDI is comparable to a full data reconstructed image via the compared methods. The remaining challenge, as shown in Fig. 9, is that the result given by the proposed method may suffer from image noise when the projection is very sparse, which may be solved by image posttreatment.

Fig. 8

Results of the diagnostic head phantom experiments from the MLEM-TV, PWLS-TV, PWLS-PICCS and PI-NDI methods with different numbers of projections. The numbers on the top indicate the number of projections used for reconstruction. Display window (0, 0.1404).

Fig. 9

The ROI-1 in the reconstructed images of the diagnostic head phantom using 8, 15 and 30 projections. The numbers on the left indicate the number of projections used for reconstruction. ROI-1 is marked in Fig. 1b. Display window (0, 0.1404).

4.3 Parameter optimization

In the proposed PI-NDI algorithm, two parameters are of great importance to be considered: the patch size n for feature extraction, which determines the length of the feature vector, and the number of neighbors k to be selected when building the image prior matrix. As discussed in section 2.B.3, only pixel-based feature extraction is needed since the prior image is of the same size as the object image matrix of 512×512, and a patch size larger than 1 may introduce error information to the reconstruction procedure. Therefore, only the selection of the number of neighbors k is discussed in this section. Theoretically speaking, a larger number of neighbors brings better results from the aspects of reconstruction accuracy and noise reduction. However, this approach naturally requires more computational cost. For the trade-off between the image quality and computational cost, an appropriate number of neighbors is needed. Figure 10 shows the reconstructed images given by the PI-NDI method with different n, namely, 8, 24, 48, 96 and 192. All the images are reconstructed from 15 views of projections. Figure 11a gives the PSNRs and MPAEs measured for all reconstructed images of the corresponding cases. The image quality is notably improved from n = 8 to n = 96 but only infinitesimally improved from n = 96 to n = 192. The computational cost, indicated by the time needed for the building of the prior matrix, shown in Fig. 11b, however, shows an explosive increase with increasing n, indicating a higher cost and lower gain when n reaches 192. Thus, in our study, the image reconstruction is implemented with a prior matrix built with 96 nearest neighbors.

Fig. 10

Reconstructed images of the XCAT ventral phantom using different numbers of neighbors k to build the prior matrix. a) Eight nearest neighbors involved; b) 24 nearest neighbors involved; c) 48 nearest neighbors involved; d) 96 nearest neighbors involved; e) 192 nearest neighbors involved.

Fig. 11

a) MPAE and PSNR measured on the reconstructed images for the cases with k = 8, 24, 48, 96 and 192. b) The consumed time for the construction of prior matrix with different k.

4.4 Influence of the quality of the prior image

In addition to the selection of the parameters, another factor that influences the performance of the proposed PI-NDI method is the quality of the prior image. The reason that the proposed method is able to reconstruct the object with ultrasparse projection data is that it absorbs the prior information strongly, which means the prior image used should be very accurate, and the misalignment between the prior image and the object image may lead to dramatic reconstruction error. A modified Shepp-Logan phantom is used to investigate the influence of the quality of the prior image on the performance of the PI-NDI method. Figure 12 shows the prior images with remarkably different quality and the corresponding reconstructed images via the PI-NDI method. The prior images are the ground truth phantom and three reconstructed images: PWLS-TV with 360 projections, FBP with 720 projections and FBP with 180 projections. It is obvious that errors in the prior image, such as noise and artifacts, are to some extent introduced to the reconstructed image. The prior image reconstructed by FBP with 180 projection contents severe noise and consequently gives a result with dramatic flecks, which may seriously affect the diagnosis from the clinical aspect.

Fig. 12

Influence of the quality of the prior image using the PI-NDI method. The images in the first row are prior images obtained by different methods (ground truth phantom, PWLS-TV, FBP with 720 views of projections and FBP with 180 views of projections), and the images in the second row are the corresponding reconstructed images from the PI-NDI method. Display window (0, 1).

4.5 Computational cost

To measure the computational cost needed for the implementation of the proposed PI-NDI method, we recorded the implementation time for the construction of the prior matrixes ψ with different numbers of neighbors, as shown in Fig. 11b, and the time needed for the iterative reconstructions using different algorithms. In our work, we used a bench top computer with Intel(R) Xeon(R) W-2102 CPU and 80GB RAM to conduct the partial parallel computation of ψ . The computational cost of the proposed PI-NDI method concentrates on the contribution of the prior matrix, which according to Fig. 11b, increases with the number of neighbors increasing. The calculation time for the iterative reconstruction with different algorithms is given in Table 4, about which 200 iterations were conducted for all the four methods. From the resultant recorded time, one can see that the proposed PI-NDI consumes more time than the compared algorithms for the sparse cases because of the involving of the prior matrix with size 252144×262144. We conclude that the reduction of computational cost is the main challenge for the clinical application of the proposed PI-NDI method, especially regarding the time needed for prior matrix construction.

Table 4
Implementation time for different methods discussed with projection data of different numbers of views

No. of views Alg. 8 15 30 60 120 360

MLEM-TV 10.59 13.09 18.21 36.96 68.73 202.84

PWLS-TV 12.77 17.25 27.09 60.72 117.57 356.40

PWLS-PICCS 17.01 21.34 31.35 64.81 121.24 358.91

PI-NDI 69.55 71.59 77.89 96.51 128.32 262.54

No. of views Alg.	8	15	30	60	120	360
MLEM-TV	10.59	13.09	18.21	36.96	68.73	202.84
PWLS-TV	12.77	17.25	27.09	60.72	117.57	356.40
PWLS-PICCS	17.01	21.34	31.35	64.81	121.24	358.91
PI-NDI	69.55	71.59	77.89	96.51	128.32	262.54

5 Conclusions

In this paper, we proposed a low-dose CT reconstruction method based on prior information of normal-dose image (PI-NDI) for CT image reconstruction from sparse-view projection data. For the pre-scanned CT image-aided reconstruction, we used pixel-based feature extraction from the prior image. We selected the number of related neighbors as 96 by successive experiments. MLEM-TV, PWLS-TV and PWLS-PICCS algorithms were compared with our PI-NDI method using a digital XCAT phantom and a diagnostic head phantom. The experimental results show that the proposed PI-NDI method is more capable for sparse-view CT reconstruction than the other methods, especially for cases with very sparse projections, from the aspect of reconstruction accuracy and noise reduction. There are several problems with the PI-NDI method, however. First, the image noise is somewhat severe, especially for the case with ultrasparse projection data. Second, the dependency on the prior image is evident, which requires a prior image with a high quality and high level of image registration. Third, the computational cost for the implementation of the PI-NDI is higher than state-of-the art methods compared in this study. Further improvement and validation of the PI-NDI method are needed for clinical application.

Conflict of interest

We declare that we do not have any commercial or associative interest that represents a conflict of interest in connection with the work submitted.

Footnotes

Acknowledgments

This work was supported by the Guangdong Special Support Program of China (2017TQ04R395), the Guangdong International Science and Technology Cooperation Project of China (2018A050506064), the National Natural Science Foundation of China (81871441), the Shenzhen International Cooperation Research Project of China (GJHZ20180928115824168), the Natural Science Foundation of Guangdong Province in China (2020A1515010733).

References

Brenner

D.J.

, Hall

E.J.

, Computed tomography—an increasing source of radiation exposure, New England Journal of Medicine 3579(22) (2007), 2277–2284.

McCollough

C.H.

, Primak

A.N.

, Braun

, et al., Strategies for reducing radiation dose in CT, Radiologic Clinics 47(1) (2009), 27–40.

Smith

, Dillon

, Gould

, Wintermark

, Radiation dose-reduction strategies for neuroradiology CT protocols, American Journal of Neuroradiology 28(9) (2007), 1628–1632.

Chen

G.H.

, Tang

, Leng

, Prior image constrained compressed sensing (PICCS): a method to accurately reconstruct dynamic CT images from highly undersampled projection data sets, Medical Physics 35(2) (2008), 660–663.

Han

, Ye

J.C.

, Framing U-Net via deep convolutional framelets: Application to sparse-view CT, IEEE Transactions on Medical Imaging 37(6) (2018), 1418–1429.

Lee

, Lee

, Kim

, et al., Deep-neural-network-based sinogram synthesis for sparse-view CT image reconstruction, IEEE Transactions on Radiation Plasma Medical Sciences 3(2) (2018), 109–119.

Ramirez-Giraldo

, Trzasko

, Leng

, et al., Nonconvex prior image constrained compressed sensing (NCPICCS): Theory and simulations on perfusion CT, Medical Physics 38(4) (2011), 2157–2167.

Zhang

, et al., Iterative image reconstruction for sparse-view CT using normal-dose image induced total variation prior, Proc SPIE 8668 (2013), 86685E.

Liu

, et al., Total variation-stokes strategy for sparse-view X-ray CT image reconstruction, IEEE Transactions on Medical Imaging 33(3) (2013), 749–763.

10.

Jiang

, Qu

, Sun

, Noise-robust Mojette reconstruction using sparse-view CT projections, Journal of X-ray Science and Technology 27(2) (2019), 237–256.

11.

Zhu

, Wahid

, Babyn

, et al., Improved compressed sensing-based algorithm for sparse-view CT image reconstruction, Computational Mathematical Methods in Medicine (2013), 185750.

12.

Sullivan

T.P.

, Jarvik

J.G.

, Cohen

W.A.

, Follow-up of conservatively managed epidural hematomas: implications for timing of repeat CT, American Journal of Neuroradiology 20(1) (1999), 107–113.

13.

De Jong

P.A.

, et al., Estimation of cancer mortality associated with repetitive computed tomography scanning, American Journal of Respiratory Critical Care Medicine 173(2) (2006), 199–203.

14.

Siewerdsen

J.H.

, et al., High-performance intraoperative cone-beam CT on a mobile C-arm: an integrated system for guidance of head and neck surgery, Proc SPIE 7261 (2009), 72610J.

15.

Pourmorteza

, Dang

, Siewerdsen

, Stayman

J.W.

, Reconstruction of difference in sequential CT studies using penalized likelihood estimation,”, Physics in Medicine Biology 61(5) (2016), 1986.

16.

Hofmann

, Schölkopf

, Smola

A.J.

, Kernel methods in machine learning,” The annals of statistics, (2008), pp. 1171–1220.

17.

Hutchcroft

, Wang

, Chen

K.T.

, Catana

, Qi

, “Anatomically-aided PET reconstruction using the kernel method,”, Physics in Medicine Biology 61(18) (2016), 6668.

18.

Wang

, Qi

, “PET image reconstruction using kernel method, IEEE Transactions on Medical Imaging 34(1) (2014), 61–71.

19.

Chuang

K.S.

, et al., A maximum likelihood expectation maximization algorithm with thresholding, Computerized Medical Imaging Graphics 29(7) (2005), 571–578.

20.

Gaitanis

, Kontaxakis

, Spyrou

, et al., PET image reconstruction: A stopping rule for the MLEM algorithm based on properties of the updating coefficients, Computerized Medical Imaging Graphics 34(2) (2010), 131–141.

21.

Lange

, Carson

, EM reconstruction algorithms for emission and transmission tomography, J Comput Assist Tomogr 8(2) (1984), 306–316.

22.

, Leahy

R.M.

, Iterative reconstruction techniques in emission computed tomography, Physics in Medicine Biology 51(15) (2006), R541.

23.

Shepp

L.A.

, Vardi

, Maximum likelihood reconstruction for emission tomography, IEEE Transactions on Medical Imaging 1(2) (1982), 113–122.

24.

Friedman

J.H.

, Bentley

J.L.

, Finkel

R.A.

, An algorithm for finding best matches in logarithmic expected time, ACM Transactions on Mathematical Software 3(3) (1977), 209–226.

25.

Buades

, Coll

, Morel

J.M.

, A non-local algorithm for image denoising, 2005 IEEE Computer Society Conference on Computer Vision and Pattern Recognition (CVPR’05) 2 (2005), 60–65.

26.

, Zhang

, Liu

, et al., A feature refinement approach for statistical interior CT reconstruction, Physics in Medicine Biology 61(14) (2016), 5311–5334.

27.

Wang

, Xu

, Shi

, A new approach for reducing beam hardening artifacts in polychromatic X-ray computed tomography using more accurate prior image, Journal of X-ray Scuience and Technology 26(4) (2018), 593–602.

28.

, Jiang

, Sun

, et al., Artifact correction in low-dose dental CT imaging using Wasserstein generative adversarial networks, Medical Physics 46(4) (2019), 1686–1696.