Abstract
Dual energy CT has the ability to give more information about the test object by reconstructing the attenuation factors under different energies. These images under different energies share identical structures but different attenuation factors. By referring to the fully sampled low-energy image, we show that it is possible to greatly reduce the sampling rate of the high-energy image in order to lower dose. To compensate the attenuation factor difference between the two modalities, we use piecewise polynomial fitting to fit the low-energy image to the high-energy image. During the reconstruction, the result is constrained by its distance to the fitted image, and the structural information thus can be preserved. An ASD-POCS-based optimization schedule is proposed to solve the problem, and numerical simulations are taken to verify the algorithm.
1. Introduction
Computed tomography has become an important nondestructive detection method in medicine, industry, and security. Typically CT scans the object by a single energy to reconstruct the attenuation factors in order to evaluate the density distributions inside the test object. However, some materials’ attenuation factors are close and hard to distinguish, which brings trouble for diagnosis. Since the attenuation factors are different under different X-ray energies, DECT [1] has been brought about to enhance the material distinguishing ability in CT. Furthermore, atomic numbers, electron densities, or specific material equivalent densities can also be reconstructed from DECT for better visualization.
In DECT, the test object is scanned under different energies while keeping the object fixed. Thus, two different images, the low-energy image and the high-energy image can be reconstructed independently from the two sets of projections, the low-energy projections and the high-energy projections . Although there are various techniques for DECT reconstruction, for example, prereconstruction [2], postreconstruction [3], and iterative reconstruction [4], we concentrate on reconstructing and in this paper to demonstrate the mathematics of our method.
Dose has been concerned more and more recently with the increasing public awareness of the possible risks brought about by the radiation of CT scans. One of the most efficient ways to reduce dose is to reduce the sampling number. According to the concept of compressed sensing (CS) [5, 6], when the sampling numbers are reduced beneath the conventional required sampling rate, one can still accurately recover the signal by incorporating prior knowledge to the reconstruction. In the instance of DECT, and share identical structures because they are taken from the same object at almost the same time. One of the straightforward strategies for dose reduction is to undersample while keeping fully sampled. During the reconstruction of , the structural information extracted from the well-reconstructed can be made utility to improve the reconstruction quality.
Although the two images share the same structure, their attenuation factors are different under the two different energies, which leads to the grey scale difference in the reconstructed images. Similar situations can be found in multimodality imaging, where the grey scale values of the images are far different from each other. Bowsher prior has been used for MR/SPECT imaging to improve the reconstruction quality of SPECT [7, 8]. For DECT, Liu et al. used image segmentation on the fully sampled to reduce the number of variables during the reconstruction of [9]. Recently, some methods invoking CS have been proposed for DECT. Szczykutowicz and Chen applied PICCS to a slow kVp switching acquisition scheme and achieved good results [10, 11]. Xing and Zheng applied ART-TV on the ratio image of to for sparser presentation of the image [12].
It has been shown that invoking reference images in CS-based reconstruction is able to improve the reconstruction quality, but grey scale value compensation remains a challenge in DECT. For example, PICCS algorithm requires the reference image and the target image to be as close as possible, but the attenuation factor may differ widely under different energies. The ratio image method, on the other hand, requires that the changes in the grey scale values are proportional so that there are less edges in the ratio image and its gradient image is sparser. However, the change of the attenuation factor under different energies is unpredictable and the conditions required for the above methods may be violated sometimes.
Here, we propose a novel CS-based method for undersampled DECT reconstruction. The well-reconstructed low-energy image is used to sparsify the undersampled high-energy image . To compensate for the grey value scale difference, both images are divided into patches and polynomial fitting is used on each pair of patches to fit the reference image to the target image . An l1-norm constraint is applied on the distance from to the fitted image, and ASD-POCS [13] is adopted for the optimization. Since the piecewise polynomial fitting is able to almost precisely approximate the target image under most occasions, promising results can be achieved when the sampling rate is greatly reduced.
The method is very much motivated by one of our previous works [14], feature constrained compressed sensing (FCCS). In the FCCS method, linear space extracted from training image set is used as the prior knowledge rather than a single image. During the reconstruction, the distance between the target image and the prior space is used as the constraint for the result. In this paper, the prior space is generated from a single reference image by taking its powers. The constraint on the distance between the target image and the polynomial space is achieved by piecewise polynomial fitting.
The paper is organized as follows. In Section 2, the mathematical principles of the algorithm are presented. Numerical simulation results are shown in Section 3. Section 4 is conclusions and discussions.
2. Methodology
2.1. An Overview of the Method
The reconstruction formula of piecewise polynomial function constrained (PPFC) method is as follows: where is the piecewise polynomial fitting (PPF) which approximates by . The images are represented by the row vectors and . is the system matrix, and is the high-energy projections. After a high-quality is reconstructed from the fully sampled low-energy projections, it is used as the prior knowledge to iteratively solve form the undersampled according to (1) by ASD-POCS algorithm.
2.2. Image Approximation by Polynomial Fitting
Before introducing PPF, we will firstly show the least squares polynomial fitting method and some of its properties. In DECT, one of the ways to sparsify by is to find a transform which makes small. We will show that piecewise polynomial fitting is a good way to approximate by , which actually makes the difference between them sparse.
Image approximation by polynomial fitting is to solve the following equation: where where is defined as the element by element power of and is an all-ones vector. is the order of the polynomial and to are the corresponding coefficients which are determined by both and .
Equation (2) can be written into matrix form, which is where is an order Vandermonde matrix. means the th element of . To further reduce the scale of (4), we adopted the following method for least squares solution: where is an order Hankel matrix. Although (6) is not the most numerically stable solution for least squares problems, it is enough for low-order polynomial fitting problems. Efficient algorithms such as the Gohberg-Kailath-Koltracht (GKK) algorithm [15] can be used to solve (6). The approximation result can be expressed as
The approximation has some good properties. Firstly, low-order polynomials have the property of smoothness, which maps similar values in to similar values in . Thus, the fitting result will not be seriously influenced by noises. Secondly, if there are not too many kinds of values presented in the images, low-order polynomials can always guarantee a good approximation. For example, if there are only 5 pairs of values presented in and , an polynomial of order 4 is able to precisely map to . If the values distribute around a few points, the approximation error will also remain small. Furthermore, low-order polynomial fitting is numerical stable. High-order least squares polynomial fitting requires high computational precision, and the result will be sensitive to rounding error. To solve this problem, we divide the image into small patches so that there are not too many kinds of grey values in each patch, and low-order polynomial fitting is applicable. Last but not least, the least squares polynomial fitting problem can be efficiently solved by the GKK algorithm, whose time and space complexities are and . Since the polynomial order is small, the method is both time and memory saving.
In our experiments, the order is set to 4, which is the minimum value required for fitting in our case. A small order makes computation with single floating point data type feasible. A higher order is applicable, but it means heavier computational load. It may be preferred when the test object is more complex; for example, the object contains more kinds of materials in a small region.
2.3. Piecewise Polynomial Fitting for the Entire Image
As we have stated, low-order polynomial fitting requires that the grey values of the image distribute around only a few points. Thus, the size of the image for fitting should be small for accuracy. This can be achieved by extracting patches from the original image and applying polynomial fitting on the sub images. On the other hand, the patches should be big enough to hold the structural information with them. There should also be adequate overlapping areas between adjacent patches or the difference between patches will lead to obvious block-like artifacts in the approximated image. In our experiments, the patch size is selected as 16 by 16, and the offset of adjacent patches is 4 by 4, which leads to an overlapping size of 12 by 12.
Figure 1 shows the PPF results on noisy phantom with different patch sizes and offsets. Noisy is used to fit the noisy . Obvious block-like artifacts can be observed when the patch size and offset are and . This is because in each small patch, the noise distributions are not identical, and without proper smoothing between adjacent patches, block-like artifacts emerge. Furthermore, the patch size should be larger than the offset; otherwise, part of the image will be missing in the result (see the result at the bottom left corner of Figure 1).
When approximating by with PPF, for each pair of patches, the high-energy image patch is approximated by the low-energy image patch according to (8). Then the entire image is generated by aggregating the patches. Define the patch selection matrix which satisfies
Then putting the patch back to its original position can be achieved by the transposed matrix . The aggregation can be expressed as
In (10), the division is element by element. is a Gaussian kernel which is used for smoothness during aggregation. is the normalized weighting for each patch. It is a diagonal matrix and which can be calculated by dividing the Gaussian kernel by the aggregated weightings.
The transpose of is
Since is a symmetric matrix, holds. Furthermore, is diagonal, so it also holds that . The transposed matrix can be written as
The transposed operator will be of future use and we will take a look at its properties here. The PPF matrix is realized by first fitting the patches and then weighted aggregation, but the transposed PPF matrix is realized by first weighting the patches, then fitting, and then unweighted aggregation. The algorithms for PPF and its transpose are shown in Algorithms 1 and 2.
| ||||||||||||||
| ||||||||||||||||||
is a small value in case the values in are almost zeros and the polynomial fitting algorithm fails. It is set at in our experiment.
2.4. Optimization by ASD-POCS
The formula for optimization is shown in (1). When adapting ASD-POCS algorithm to the problem, the only change to make is the first derivative of the objective function. The objective function is
Using the chain derivative rule, one can easily derive the gradient of the objective function: where is the element by element sign function: where is a small regularization factor near zero. In our experiments, this value is set at 0.01.
The initial value for the iterations is crucial for the convergence speed. For PPFC, 10 times SART is used to estimate and PPF is used to approximate the current image by . The approximation result is used as the initial value of the iterations.
The sparsity of the PPFC transform is shown in Figure 2. It can be seen that the PPFC transform on the well-reconstructed image is sparse, while the zero entities are much more in the image with artifacts.
3. Simulations
3.1. Simulation Method
Multienergy projections are used for the experiments. The spectrum of the X-ray source is generated by the Monte Carlo method. Three different spectrums are used for the experiments: a 90 kVp spectrum, a 120 kVp spectrum, and a beam-hardened 160 kVp spectrum. There are two different phantoms for testing, a cylinder phantom and a realistic dental phantom. The cylinder phantom is forward projected by the 90 kVp spectrum and the 160 kVp spectrum. The energy used for the dental phantom is 90 kVp and 120 kVp. The low-energy image is reconstructed by FBP and the high-energy image is reconstructed by the proposed PPFC method with only 15 projections being uniformly sampled on 360 degrees. PICCS is also realized under the ASD-POCS framework as a reference method and its initial value is set at the energy normalized image . The simulation phantoms are shown in Figure 3.
(a)
(b)
3.2. Cylinder Phantom Experiments
The cylinder phantom is first forward projected by the 90 kVp and 160 kVp spectrums with fan beam geometry to get the high and low projections and . Then the low-energy image is reconstructed from 900 projections by FBP. The region outside the ROI of FBP is set at zero. For the high-energy image , only 15 projections uniformly sampled across 360 degrees are used for reconstruction. Then both PPFC and PICCS are employed to reconstruct the phantom. Furthermore, the algorithms are also tested against noises. The noises are added to the projections by setting the photon number of the X-ray source at per detector bin per projection. The scan geometry is shown in Table 1, and the simulation results are shown in Figure 4 and Table 2.
In the results of PICCS, the cylinder made of 1% NaI solution is blurred. The reason is that in the reference image , the attenuation factor of the 1% NaI solution is larger than the background PMMA cylinder, but in the target image , the attenuation factor of the 1% NaI solution becomes smaller than the PMMA’s attenuation factor. The way PICCS uses the prior image by subtraction actually reduces the sparsity on the cylinder of 1% NaI solution. As its consequence, PICCS is not able to recover the cylinder of 1% NaI solution well while other cylinders are well preserved.
As for the results of PPFC, all the cylinders including the 1% NaI solution cylinder are well recovered. Furthermore, the results of the noisy projections have not degraded much comparing to the noiseless results. Thus, the experiments indicate that PPFC is compensating grey scale values difference and noise stable.
3.3. Realistic Phantom Experiments
The dental phantom is forward projected by the 90 kVp and 120 kVp spectrums with fan beam geometry. The low-energy image is reconstructed by FBP from 720 projections. is downsampled to 15 projections and reconstructed by PPFC and PICCS. Noises with initial photons are also added to the projections. The scan geometry is shown in Table 3, and the corresponding results are shown in Figure 5 and Table 4.
The experiments on the realistic phantom show that PPFC is able to reconstruct objects with complicated structures as well as the simple objects. It also shows an advantage over PICCS on the aspect of RMSE.
4. Conclusion and Discussions
In this paper, we propose a CS-based method for undersampled DECT reconstruction with piecewise polynomial function constraint. The low-energy image is reconstructed from fully sampled projections and the high-energy image is reconstructed from severely corrupted samples with the well-reconstructed low-energy image as the reference. The proposed piecewise polynomial fitting method has good ability to compensate for the grey scale value difference between the high- and low-energy images. Under most conditions, the target image can be well approximated by the reference image using the PPF, which ensures the sparsity of the PPFC transform. The simulation results show that our method is both accurate and stable.
The drawback of the algorithm is that the piecewise polynomial fitting is still not efficient enough. However, the fitting of each patch is independent and the algorithm can be further accelerated by parallel computation. Furthermore, the algorithm has the potential to reconstruct the decomposition coefficients images in DECT, whose values are far from the values in the low-energy image. Applying the method to a dual-effect or dual-material decomposition reconstruction is of future concerns.
Acknowledgments
This work was partly supported by the grants from NNSFC 10905030 and the Beijing Natural Science Foundation (research on key techniques of medical cone-beam CT reconstruction from little data based on compressed sensing theory).