# Prior Image Guided Undersampled Dual Energy Reconstruction with Piecewise Polynomial Function Constraint

Authors: Dufan Wu; Li Zhang; Liang Li; Le Shen; Yuxiang Xing
Journal: Computational and Mathematical Methods in Medicine
 (2013)

---

## 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.


---

## Body

## 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 imagexL and the high-energy image xH can be reconstructed independently from the two sets of projections, the low-energy projections pL and the high-energy projections pH. Although there are various techniques for DECT reconstruction, for example, prereconstruction [2], postreconstruction [3], and iterative reconstruction [4], we concentrate on reconstructing xL and xH 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, xL and xH 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 xH while keeping xL fully sampled. During the reconstruction of xH, the structural information extracted from the well-reconstructed xL 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 xL to reduce the number of variables during the reconstruction of xH [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 xL to xH 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 imagexL is used to sparsify the undersampled high-energy image xH. 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 xL to the target image xH. An l1-norm constraint is applied on the distance from xH 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 Section2, 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:(1)min∥xH-φ(xL,xH)∥1s.t.∥AxH-pH∥2≤εxH≥0,
where φ(xL,xH) is the piecewise polynomial fitting (PPF) which approximates xH by xL. The images are represented by the row vectors xH and xL. A is the system matrix, and pH is the high-energy projections. After a high-quality xL is reconstructed from the fully sampled low-energy projections, it is used as the prior knowledge to iteratively solve xH form the undersampled pH 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 sparsifyxH by xL is to find a transform φ(·) which makes ∥xH-φ(xL,xH)∥0 small. We will show that piecewise polynomial fitting is a good way to approximate xH by xL, which actually makes the difference between them sparse.Image approximation by polynomial fitting is to solve the following equation:(2)min∥xH-f(xL,xH)∥2,
where
(3)f(xL,xH)=∑k=0p-1akxLk,
where xk is defined as the element by element power of x and x0 is an all-ones vector. (p-1) is the order of the polynomial and a0 to ap-1 are the corresponding coefficients which are determined by both xH and xL.Equation (2) can be written into matrix form, which is
(4)a=argmina∥VTa-xH∥2,
where
(5)V=(11⋯1(xL)1(xL)2⋯(xL)n⋮⋮⋱⋮(xL)1p-1(xL)2p-1⋯(xL)np-1)
is an order p Vandermonde matrix. (xL)i means the ith element of xL. To further reduce the scale of (4), we adopted the following method for least squares solution:
(6)Ha=VxH,
where
(7)Hij:=(VVT)ij=∑k=1n(xL)ki+j-2
is an order p 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
(8)f(xL,xH)=VTa=VTH-1VxH=P(xL)xH.The approximation has some good properties. Firstly, low-order polynomials have the property of smoothness, which maps similar values inxL to similar values in f(xL,xH). 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 xL and xH, an polynomial of order 4 is able to precisely map xL to xH. 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 O(np)+O(p2) and O(n)+O(p2). Since the polynomial order p is small, the method is both time and memory saving.In our experiments, the orderp 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.Figure1 shows the PPF results on noisy phantom with different patch sizes and offsets. Noisy xL is used to fit the noisy xH. Obvious block-like artifacts can be observed when the patch size and offset are (8,8) and (16,16). 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).Figure 1
The PPF results on noisy phantom with different patch sizes and offsets. From the top row to the bottom row, the corresponding offsets between adjacent patches are 4, 8, and 16. From the leftmost column to the rightmost column, the corresponding patch sizes are 8, 16, and 32. The grey scale value window is[0.185,0.192].When approximatingxH by xL with PPF, for each pair of patches, the high-energy image patch xH(s) is approximated by the low-energy image patch xL(s) according to (8). Then the entire image is generated by aggregating the patches. Define the patch selection matrix E(s) which satisfies
(9)xH(s)=E(s)xH.Then putting the patch back to its original position can be achieved by the transposed matrixE(s)T. The aggregation can be expressed as
(10)φ(xL,xH)=∑sE(s)TGP(xL(s))E(s)xH∑sE(s)TGE(s)·1=(∑sE(s)TW(s)P(xL(s))E(s))xH=BxH.In (10), the division is element by element. G is a Gaussian kernel which is used for smoothness during aggregation. W(s) is the normalized weighting for each patch. It is a diagonal matrix and
(11)Wii(s)=1(E(s)∑lE(l)TGE(l)·1)iGii,
which can be calculated by dividing the Gaussian kernel G by the aggregated weightings.The transpose ofB is
(12)BT=∑sE(s)TPT(xL(s))W(s)TE(s).SinceP=VTH-1V=VT(VVT)-1V is a symmetric matrix, PT=P holds. Furthermore, W(s) is diagonal, so it also holds that W(s)T=W(s). The transposed matrix can be written as
(13)BT=∑sE(s)TP(xL(s))W(s)E(s).The transposed operator will be of future use and we will take a look at its properties here. The PPF matrixB is realized by first fitting the patches and then weighted aggregation, but the transposed PPF matrix BT 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.Algorithm 1:Image approximation by PPF.
(1)Let φ=0,w=0;
(2)For each patch pair xL(s) and xH(s)
(3)if (||xL(s)||2≥δ) Calculate φ(s)=P(xL(s))xH(s) by (8);
(4)else φ(s)=0;
(5)φ=φ+E(s)TGφ(s),w=w+E(s)TG·1, aggregate the weighted images;
(6)End for
(7)Output φ/w;Algorithm 2:Transposed PPF.
(1)Let φ=0,  w=0;
(2)Calculate aggregated weighting w=∑sE(s)TG·1;
(3)For each patch pair xL(s) and xH(s)
(4)Calculate W(s)=G/(E(s)w);
(5)if (||W(s)xL(s)||2≥δ)φ(s)=P(W(s)xL(s))W(s)xH(s);
(6)else φ(s)=0;
(7)φ=φ+E(s)Tφ(s);
(8)End for
(9)Output φ;δ is a small value in case the values in xL(s) are almost zeros and the polynomial fitting algorithm fails. It is set at 1×10-4 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
(14)g(xH)=∥xH-φ(xL,xH)∥1=∥(I-B)xH∥1.Using the chain derivative rule, one can easily derive the gradient of the objective function:(15)∇g(xH)=(I-BT)sgn((I-B)xH),
where sgn(·) is the element by element sign function:
(16)sgni(x)={-1xi≤-δ2,xiδ2-δ2<xi≤δ2,1xi>δ2,
where δ2 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 estimatexH and PPF is used to approximate the current image x^H by xL. The approximation result φ(xL,x^H) is used as the initial value of the iterations.The sparsity of the PPFC transform(I-B)xH 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.Figure 2
The approximation result of PPF.xL is used to approximate xH. The top row is fully sampled. The bottom row is sparsely sampled by 15 projections. The left most column is the reconstruction result by FBP. The middle column is the PPF approximation results. The right most column is the difference between the reconstructed images and the approximation results, which is the PPFC transform. The grey scale window for the first and second columns is [0.15,0.25]. For the third column, the window is [0,0.01].


## 2.1. An Overview of the Method

The reconstruction formula of piecewise polynomial function constrained (PPFC) method is as follows:(1)min∥xH-φ(xL,xH)∥1s.t.∥AxH-pH∥2≤εxH≥0,
where φ(xL,xH) is the piecewise polynomial fitting (PPF) which approximates xH by xL. The images are represented by the row vectors xH and xL. A is the system matrix, and pH is the high-energy projections. After a high-quality xL is reconstructed from the fully sampled low-energy projections, it is used as the prior knowledge to iteratively solve xH form the undersampled pH 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 sparsifyxH by xL is to find a transform φ(·) which makes ∥xH-φ(xL,xH)∥0 small. We will show that piecewise polynomial fitting is a good way to approximate xH by xL, which actually makes the difference between them sparse.Image approximation by polynomial fitting is to solve the following equation:(2)min∥xH-f(xL,xH)∥2,
where
(3)f(xL,xH)=∑k=0p-1akxLk,
where xk is defined as the element by element power of x and x0 is an all-ones vector. (p-1) is the order of the polynomial and a0 to ap-1 are the corresponding coefficients which are determined by both xH and xL.Equation (2) can be written into matrix form, which is
(4)a=argmina∥VTa-xH∥2,
where
(5)V=(11⋯1(xL)1(xL)2⋯(xL)n⋮⋮⋱⋮(xL)1p-1(xL)2p-1⋯(xL)np-1)
is an order p Vandermonde matrix. (xL)i means the ith element of xL. To further reduce the scale of (4), we adopted the following method for least squares solution:
(6)Ha=VxH,
where
(7)Hij:=(VVT)ij=∑k=1n(xL)ki+j-2
is an order p 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
(8)f(xL,xH)=VTa=VTH-1VxH=P(xL)xH.The approximation has some good properties. Firstly, low-order polynomials have the property of smoothness, which maps similar values inxL to similar values in f(xL,xH). 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 xL and xH, an polynomial of order 4 is able to precisely map xL to xH. 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 O(np)+O(p2) and O(n)+O(p2). Since the polynomial order p is small, the method is both time and memory saving.In our experiments, the orderp 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.Figure1 shows the PPF results on noisy phantom with different patch sizes and offsets. Noisy xL is used to fit the noisy xH. Obvious block-like artifacts can be observed when the patch size and offset are (8,8) and (16,16). 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).Figure 1
The PPF results on noisy phantom with different patch sizes and offsets. From the top row to the bottom row, the corresponding offsets between adjacent patches are 4, 8, and 16. From the leftmost column to the rightmost column, the corresponding patch sizes are 8, 16, and 32. The grey scale value window is[0.185,0.192].When approximatingxH by xL with PPF, for each pair of patches, the high-energy image patch xH(s) is approximated by the low-energy image patch xL(s) according to (8). Then the entire image is generated by aggregating the patches. Define the patch selection matrix E(s) which satisfies
(9)xH(s)=E(s)xH.Then putting the patch back to its original position can be achieved by the transposed matrixE(s)T. The aggregation can be expressed as
(10)φ(xL,xH)=∑sE(s)TGP(xL(s))E(s)xH∑sE(s)TGE(s)·1=(∑sE(s)TW(s)P(xL(s))E(s))xH=BxH.In (10), the division is element by element. G is a Gaussian kernel which is used for smoothness during aggregation. W(s) is the normalized weighting for each patch. It is a diagonal matrix and
(11)Wii(s)=1(E(s)∑lE(l)TGE(l)·1)iGii,
which can be calculated by dividing the Gaussian kernel G by the aggregated weightings.The transpose ofB is
(12)BT=∑sE(s)TPT(xL(s))W(s)TE(s).SinceP=VTH-1V=VT(VVT)-1V is a symmetric matrix, PT=P holds. Furthermore, W(s) is diagonal, so it also holds that W(s)T=W(s). The transposed matrix can be written as
(13)BT=∑sE(s)TP(xL(s))W(s)E(s).The transposed operator will be of future use and we will take a look at its properties here. The PPF matrixB is realized by first fitting the patches and then weighted aggregation, but the transposed PPF matrix BT 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.Algorithm 1:Image approximation by PPF.
(1)Let φ=0,w=0;
(2)For each patch pair xL(s) and xH(s)
(3)if (||xL(s)||2≥δ) Calculate φ(s)=P(xL(s))xH(s) by (8);
(4)else φ(s)=0;
(5)φ=φ+E(s)TGφ(s),w=w+E(s)TG·1, aggregate the weighted images;
(6)End for
(7)Output φ/w;Algorithm 2:Transposed PPF.
(1)Let φ=0,  w=0;
(2)Calculate aggregated weighting w=∑sE(s)TG·1;
(3)For each patch pair xL(s) and xH(s)
(4)Calculate W(s)=G/(E(s)w);
(5)if (||W(s)xL(s)||2≥δ)φ(s)=P(W(s)xL(s))W(s)xH(s);
(6)else φ(s)=0;
(7)φ=φ+E(s)Tφ(s);
(8)End for
(9)Output φ;δ is a small value in case the values in xL(s) are almost zeros and the polynomial fitting algorithm fails. It is set at 1×10-4 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
(14)g(xH)=∥xH-φ(xL,xH)∥1=∥(I-B)xH∥1.Using the chain derivative rule, one can easily derive the gradient of the objective function:(15)∇g(xH)=(I-BT)sgn((I-B)xH),
where sgn(·) is the element by element sign function:
(16)sgni(x)={-1xi≤-δ2,xiδ2-δ2<xi≤δ2,1xi>δ2,
where δ2 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 estimatexH and PPF is used to approximate the current image x^H by xL. The approximation result φ(xL,x^H) is used as the initial value of the iterations.The sparsity of the PPFC transform(I-B)xH 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.Figure 2
The approximation result of PPF.xL is used to approximate xH. The top row is fully sampled. The bottom row is sparsely sampled by 15 projections. The left most column is the reconstruction result by FBP. The middle column is the PPF approximation results. The right most column is the difference between the reconstructed images and the approximation results, which is the PPFC transform. The grey scale window for the first and second columns is [0.15,0.25]. For the third column, the window is [0,0.01].


## 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 imagexL is reconstructed by FBP and the high-energy image xH 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 (∥x^H∥2/∥xL∥2)xL. The simulation phantoms are shown in Figure 3.The phantoms used for simulation: (a) the cylinder phantom; the base cylinder is PMMA and the inner cylinders from the biggest to the smallest are made of bone, water, and 1%, 2%, 5%, and 10% NaI solutions and (b) the dental phantom.
(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 projectionspH and pL. Then the low-energy image xL is reconstructed from 900 projections by FBP. The region outside the ROI of FBP is set at zero. For the high-energy image xH, 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 2×105 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.Table 1
The scan geometry for the cylinder phantom.
Parameter
Value
Source to center distance
75 cm
Source to detector distance
105 cm
Total width of detector
9.6 cm
Number of detector bins
512
Reconstruction resolution
256
×
256
Reconstruction size
6.4 cm× 6.4 cmTable 2
The quantitative estimations of the reconstruction results of the cylinder phantom.
Method
RMSE
Iteration number
Time (s)
Noiseless PPFC
2.9
×
1
0
-
4
10
31.2
Noiseless PICCS
4.8
×
1
0
-
4
100
34.3
Noisy PPFC
2.9
×
1
0
-
4
10
31.2
Noisy PICCS
8.7
×
1
0
-
4
100
34.4
The RMSE is calculated by comparing to the noiseless result of FBP.
The iterations numbers do not include the 10 times pre-iterations and are enough for convergence.
The code is implemented by C++ on a machine with Intel Core i5 processor and 3 GB memory.Figure 4
The reconstruction results of cylinder phantoms. The top row is the results from noiseless projections and the bottom row is the results from noisy projections. From left to right: the first column is the fully sampledxL reconstructed by FBP, the second column is the fully sampled xH reconstructed by FBP, the third column is the reconstruction results of undersampled xH by PICCS, and the fourth column is the results from PPFC, the grey scale windows for the low-energy images and high-energy images are [0.1,0.6] and [0.15,0.25].In the results of PICCS, the cylinder made of 1% NaI solution is blurred. The reason is that in the reference imagexL, the attenuation factor of the 1% NaI solution is larger than the background PMMA cylinder, but in the target image xH, 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 imagexL is reconstructed by FBP from 720 projections. xH is downsampled to 15 projections and reconstructed by PPFC and PICCS. Noises with 1×105 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.Table 3
The scan geometry for the dental phantom.
Parameter
Value
Source to center distance
50 cm
Source to detector distance
80 cm
Total width of detector
30 cm
Number of detector bins
960
Reconstruction resolution
512
×
512
Reconstruction size
19.2 cm× 19.2 cmTable 4
The quantitative estimations of the reconstruction results of the dental phantom.
Method
RMSE
Iteration number
Time (s)
PPFC
1.1
×
1
0
-
3
20
238
PICCS
2.5
×
1
0
-
3
200
288
The RMSE is compared to the noiseless results of FBP.Figure 5
The reconstruction results of the dental phantom. From left to right: the first one is thexL reconstructed from 720 projections by FBP, the second one is the xH reconstructed from 70 projections by FBP, the third one is the results by PICCS from 15 projections, and the last one is the results of PPFC. The grey scale windows for the low and high-energy images are [0, 0.85] and [0, 0.65]. All the results are reconstructed from the noisy projections.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.


## 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 imagexL is reconstructed by FBP and the high-energy image xH 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 (∥x^H∥2/∥xL∥2)xL. The simulation phantoms are shown in Figure 3.The phantoms used for simulation: (a) the cylinder phantom; the base cylinder is PMMA and the inner cylinders from the biggest to the smallest are made of bone, water, and 1%, 2%, 5%, and 10% NaI solutions and (b) the dental phantom.
(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 projectionspH and pL. Then the low-energy image xL is reconstructed from 900 projections by FBP. The region outside the ROI of FBP is set at zero. For the high-energy image xH, 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 2×105 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.Table 1
The scan geometry for the cylinder phantom.
Parameter
Value
Source to center distance
75 cm
Source to detector distance
105 cm
Total width of detector
9.6 cm
Number of detector bins
512
Reconstruction resolution
256
×
256
Reconstruction size
6.4 cm× 6.4 cmTable 2
The quantitative estimations of the reconstruction results of the cylinder phantom.
Method
RMSE
Iteration number
Time (s)
Noiseless PPFC
2.9
×
1
0
-
4
10
31.2
Noiseless PICCS
4.8
×
1
0
-
4
100
34.3
Noisy PPFC
2.9
×
1
0
-
4
10
31.2
Noisy PICCS
8.7
×
1
0
-
4
100
34.4
The RMSE is calculated by comparing to the noiseless result of FBP.
The iterations numbers do not include the 10 times pre-iterations and are enough for convergence.
The code is implemented by C++ on a machine with Intel Core i5 processor and 3 GB memory.Figure 4
The reconstruction results of cylinder phantoms. The top row is the results from noiseless projections and the bottom row is the results from noisy projections. From left to right: the first column is the fully sampledxL reconstructed by FBP, the second column is the fully sampled xH reconstructed by FBP, the third column is the reconstruction results of undersampled xH by PICCS, and the fourth column is the results from PPFC, the grey scale windows for the low-energy images and high-energy images are [0.1,0.6] and [0.15,0.25].In the results of PICCS, the cylinder made of 1% NaI solution is blurred. The reason is that in the reference imagexL, the attenuation factor of the 1% NaI solution is larger than the background PMMA cylinder, but in the target image xH, 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 imagexL is reconstructed by FBP from 720 projections. xH is downsampled to 15 projections and reconstructed by PPFC and PICCS. Noises with 1×105 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.Table 3
The scan geometry for the dental phantom.
Parameter
Value
Source to center distance
50 cm
Source to detector distance
80 cm
Total width of detector
30 cm
Number of detector bins
960
Reconstruction resolution
512
×
512
Reconstruction size
19.2 cm× 19.2 cmTable 4
The quantitative estimations of the reconstruction results of the dental phantom.
Method
RMSE
Iteration number
Time (s)
PPFC
1.1
×
1
0
-
3
20
238
PICCS
2.5
×
1
0
-
3
200
288
The RMSE is compared to the noiseless results of FBP.Figure 5
The reconstruction results of the dental phantom. From left to right: the first one is thexL reconstructed from 720 projections by FBP, the second one is the xH reconstructed from 70 projections by FBP, the third one is the results by PICCS from 15 projections, and the last one is the results of PPFC. The grey scale windows for the low and high-energy images are [0, 0.85] and [0, 0.65]. All the results are reconstructed from the noisy projections.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.



---

*Source: 437917-2013-10-27.xml*