Interpolation of CT Projections by Exploiting Their Self-Similarity and Smoothness

As the medical usage of computed tomography (CT) continues to grow, the radiation dose should remain at a low level to reduce the health risks. Therefore, there is an increasing need for algorithms that can reconstruct high-quality images from low-dose scans. In this regard, most of the recent studies have focused on iterative reconstruction algorithms, and little attention has been paid to restoration of the projection measurements, i.e., the sinogram. In this paper, we propose a novel sinogram interpolation algorithm. The proposed algorithm exploits the self-similarity and smoothness of the sinogram. Sinogram self-similarity is modeled in terms of the similarity of small blocks extracted from stacked projections. The smoothness is modeled via second-order total variation. Experiments with simulated and real CT data show that sinogram interpolation with the proposed algorithm leads to a substantial improvement in the quality of the reconstructed image, especially on low-dose scans. The proposed method can result in a significant reduction in the number of projection measurements. This will reduce the radiation dose and also the amount of data that need to be stored or transmitted, if the reconstruction is to be performed in a remote site.


Introduction
There is a great interest in reducing the radiation dose in computed tomography (CT). In practice, this can be done by either lowering the x-ray photon current or reducing the number of projection measurements taken. However, reconstructing a diagnostically useful image from such noisy and/or undersampled measurements is very challenging. As a result, recent years have witnessed a substantial growth in research on more sophisticated image reconstruction and processing algorithms for CT. Most of the published research have focused on iterative reconstruction algorithms [21] [14] [16] or on image-domain denoising and restoration algorithms [7] [12] [5] [25] [23] [11]. Comparatively, many fewer studies have been reported on denoising [10] [15] [13], or otherwise improving [8], the CT projection measurements (i.e., the sinogram). In this paper, we propose a sinogram interpolation algorithm for cone-beam CT. Our algorithm exploits both the smoothness and the self-similarity of the sinogram. We apply the proposed algorithm on simulated and real CT scans and compare it with an algorithm that is based on learned dictionaries.

Materials and Methods
A schematic of the cone-beam CT and a sample sinogram of a brain phantom are shown in Figure 1(a). The variation of the photon flux incident on the detectors follows a Poisson distribution. However, after the logarithm transformation, the noise follows approximately a Gaussian distribution with signaldependent variance. Specifically, let us denote the true line integral of the attenuation coefficient along the line from the x-ray source to the detector indexed with (i, j) in the k th projection with y t (i, j, k). Then the noisy measurement We assume that only a portion of the desired projections have been directly measured and the rest are to be estimated (i.e., interpolated). Due to the lack of space and in order to simplify the presentation of the algorithm, we focus on the case depicted in Figure 1(b). Specifically, we assume that half of the n θ desired projection views have been measured and the remaining half are to be estimated. In other words, the measured sinogram is y n ∈ R nu×nv×n θ /2 and we would like to estimate the interpolated sinogram y ∈ R nu×nv×n θ . As we will explain later, even though our main goal is sinogram interpolation, the proposed algorithm also has an excellent denoising effect. Therefore, we estimate the full set of n θ projection views, not only the missing n θ /2 views.
We propose to estimate y by minimizing the following cost function: The first term in J is obtained simply by maximizing the log-likelihood of the measurements. In this term, M is a binary mask matrix that removes from y those projection views that have not been measured and W is a diagonal weight matrix whose diagonal elements are inversely proportional to the measurement variance. The regularization functions R s and R h will be explained below, and λ s and λ h are regularization parameters.

Regularization in terms of sinogram self-similarity, R s
In the past decade, patch-based models have emerged as one of the most powerful models in image processing. Much of the research on these models was inspired by the success of the non-local means (NLM) denoising algorithm [2]. Let us denote the noisy image with x = x 0 + w, where x 0 and w are the true image and the noise, respectively. We also denote the i th pixel of x with x(i) and a patch centered on x(i) with x[i]. NLM estimates the value of x 0 (i) as follows: where G a denotes a Gaussian kernel with bandwidth a and N is the total number of pixels. The intuition behind this algorithm is that patches that are similar must have similar pixels at their centers. This simple idea has led to some of the best image denoising algorithms and has produced state-of-the-art results in many other image processing applications as well. Many studies have suggested solving various types of inverse problems by introducing regularization terms that are motivated by the same idea [20,24]. Because this type of self-similarity is very abundant in sinogram (even more than in natural images, as can be seen in the sample sinogram in Figure 1(a)), we suggest a similar form for R s : where we have used simplified pixel index notation I = (i, j, k) and I = (i , j , k ). Computation of y * has a few differences with the basic NLM in (2). Firstly, we work with 3D blocks instead of 2D patches. We stack the 2D projections to form a 3D image, as shown in Figure 1(b), and work with small blocks of this image. This will allow us to exploit both the correlation between adjacent pixels within a projection as well as the correlation between pixels in adjacent projection views. Secondly, for computing y * (i, j, k), we first find a small number of blocks that are very similar to y[i, j, k] and use those blocks only, instead of all blocks in the image. In (3), Ω i,j,k denotes the indices of these blocks. Thirdly, unlike (2) where patch similarities in the same image are exploited, in (3) we use a second image, as shown in Figure 2. This second image, denoted with z in Equation (3) and Figure 2 is built by stacking the projections of the normal-dose (i.e., low-noise) scan of a similar object. For instance, this can be a normal-dose prior scan of the same patient or of a different patient. Our experience shows that with a proper choice of z, this approach works much better than finding similar blocks from the same noisy scan. Therefore, computation of R s requires that for each pixel y(i, j, k) we find a set of blocks sufficiently similar to y[i, j, k] in z. We use the Generalized PatchMatch algorithm [1] for this purpose. We should also note that in finding similar blocks and in computing the block differences in (3), we only include the pixels from projections that have been directly measured. Block matching between the noisy scan to be restored y and the high-dose reference scan z.

Regularization in terms of sinogram smoothness, R h
An important characteristic of sinogram is its smoothness. We model this smoothness via the 1 norm of the Hessian of the sinogram. This is a generalization of the 1 norm of the gradient, also known as total variation (TV). TV is a good model for piecewise-constant images, but for smooth signals such as sinogram, it leads to staircase artifacts and higher-order differentials need to be considered [18]. Therefore, we suggest the following form for R h : In other words, we compute the 2D Hessians in the three orthogonal planes and add their 1 norms. We compute |∇ 2 uv y(i, j, k)| as: where Ds are second-order difference operators. Expressions for |∇ 2 vθ y(i, j, k)| and |∇ 2 θu y(i, j, k)| are very similar.

Optimization algorithm
The proposed cost function has this form: where we have dropped the pixel indices to simplify the expressions. An estimate of the interpolated (and denoised) projections is obtained as a minimizer of J(y).
To perform this minimization, we use the split Bregman iterative algorithm [6]. The split Bregman method first converts the unconstrained optimization problem of minimizing J(y) into a constrained problem: subject to: f = y, g 1 = ∇ 2 uv y, g 2 = ∇ 2 vθ y, g 3 = ∇ 2 θu y This constrained optimization problem is solved via Bregman iteration: In this algorithm, µ i are algorithm parameters and b i are auxiliary variables. The advantage of this reformulation is that the above minimization problem can be split into five smaller problems, one for each of the five variables that can be solved more easily. Minimization with respect to f has a simple closed-form solution. Minimizations with respect to g 1 , g 2 , and g 3 also have simple softthresholding solutions [18]. Minimization with respect to y is more difficult. We solve this optimization by considering each of the three Hessian terms in turn, each of which can be approximately solved using the Gauss-Seidel method. The regularization parameters λ s and λ h influence the recovered solution, whereas µ 1 and µ 2 influence the convergence speed. A study of the effect of these parameters is beyond the limits of this paper. A discussion of the role of these parameters for image restoration with second-order TV can be found in [19,18]. In our experiments we used λ h = 0.0001, µ 1 = µ 2 = 0.001, and λ s = 1.

Results and Discussion
As we mentioned above, unlike the iterative reconstruction methods that aim at reconstructing a high-quality image from undersampled projections, our goal is to estimate the "missing" projections. Therefore, in all of the experiments reported here we used the FDK algorithm for image reconstruction [4] [9]. We compare our algorithm with the dictionary-based sinogram interpolation method proposed in [17], which was shown to be better than spline interpolation.

Experiment with simulated data
We first applied our algorithm on scans simulated from a brain phantom, which we obtained from the BrainWeb database [3]. We simulated n θ projections from this phantom, for two values of n θ = 1440 and 960. For each n θ , we first reconstructed the image of the phantom from the full set of n θ projections and from n θ /2 projections; we denote these images with x n θ and x n θ /2 , respectively. We then applied the proposed algorithm and the dictionary-based interpolation algorithm to interpolate the subset of n θ /2 projections to generate n θ projections and reconstructed the image of the phantom from the interpolated projections. We will denote these images with x proposed n θ /2 and x dict.
n θ /2 . We simulated two levels of noise in the projections with different number of incident photons: N 0 = 10 6 and N 0 = 5 × 10 4 . We will refer to these simulations as low-noise and high-noise, respectively. For both simulations, we assumed the detector electronic noise to be additive Gaussian with a standard deviation of 40. As the reference scan that we need for block matching for computation of R s , we used the simulated scan of a different brain phantom from the same database.
We compared the reconstructed images with the true phantom image by computing the root-mean-square of error (RMSE) and the Structural Similarity (SSIM) index. The results of this comparison are presented in Table 1. Sinogram interpolation with the propsoed algorithm has resulted in a large improvement in the objective quality of the reconstructed image. The improvement is more substantial in the case of high-noise projections. This is because, as we mentioned above, both regularization terms R s and R h have excellent denoising effects. The proposed algorithm has also outperformed the interpolation algorithm based on learned dictionaries. The objective quality of x proposed n θ /2 is very close to x n θ in the low-noise case and better than x n θ in the high-noise case. As shown in Figure 3, sinogram interpolation with the proposed algorithm has resulted in a substantial improvement in the visual quality of the reconstructed image. Not only artifacts have been significantly reduced, noise has also been decreased substantially.  Table 1. Objective quality of the reconstructed images of the brain phantom.

Experiment with real CT data
We used two micro-CT scans of a rat for this experiment: (1) a normal-dose scan obtained at normal dose used for routine imaging, (2) a low-dose scan obtained at much reduced dose (by reducing the mAs setting to half of that in routine imaging and using additional copper filtration). Each scan consisted of 720 projections. We reconstructed a high-quality "reference" image of the rat from all 720 projections of the normal-dose scan. Then, for each of the two scans we selected subsets of n θ projections for two values of n θ = 360 and 180 and, for each n θ , performed the same analysis as described above for simulation experiment. For the block-matching in the computation of R s we used the scan of a different rat.
To assess the quality of the reconstructed images, we computed the RMSE and SSIM with respect to the reference image as well as the contrast-to-noise ratio (CNR). The results of this evaluation have been summarized in Table 2. There is a substantial improvement in image quality metrics as a result of sinogram interpolation with the proposed algorithm. The proposed algorithm has led to much better image quality than the dictionary-based sinogram interpolation. The gain in image quality is more pronounced in the low-dose case. In fact, the objective quality of x proposed n θ /2 is even better than x n θ on the low-dose scan. For a visual comparison, Figure 4 shows a slice from the reconstructed images of the rat from the low-dose scan with n θ = 240. There is a remarkable improvement in image quality as a result of sinogram interpolation with the proposed algorithm. The visual quality of x proposed n θ /2 seems to be better than both x dict. n θ /2 and x n θ . n θ = 360 n θ = 240 xn θ x n θ /2 x proposed n θ /2  Table 2. Objective quality of the reconstructed images of the rat.

Conclusion
Our experiments show that the proposed sinogram interpolation algorithm can lead to large improvements in image quality. Due to the additional denoising effect of both regularization functions used by the proposed algorithm, this improvement is more significant when applied on low-dose scans. This means that the proposed algorithm is especially well-suited for sinogram restoration in lowdose CT. Our other experiments, not reported in this paper because of space limitations, show that the proposed algorithm can also be used to effectively interpolate the sinogram in more general cases than the case considered in this paper, for example when the angular spacing of the missing projection views is non-uniform or when some detector measurements are corrupted. The experimental results suggest that our proposed method can reduce the number of projection measurements that are needed to reconstruct an image with diagnostic quality. The reduction in the number of projection measurements also means that less data will need to be transmitted in situations where the image reconstruction is performed in a different site.