PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
Phys Med Biol. Author manuscript; available in PMC 2017 December 21.
Published in final edited form as:
PMCID: PMC5149107
NIHMSID: NIHMS832887

An Adaptive Reconstruction Algorithm for Spectral CT Regularized by a Reference Image

Abstract

The photon counting detector based spectral CT system is attracting increasing attention in the CT field. However, the spectral CT is still premature in terms of both hardware and software. To reconstruct high quality spectral images from low-dose projections, an adaptive image reconstruction algorithm is proposed that assumes a known reference image (RI). The idea is motivated by the fact that the reconstructed images from different spectral channels are highly correlated. If a high quality image of the same object is known, it can be used to improve the low-dose reconstruction of each individual channel. This is implemented by maximizing the patch-wise correlation between the object image and the RI. Extensive numerical simulations and preclinical mouse study demonstrate the feasibility and merits of the proposed algorithm. It also performs well for truncated local projections, and the surrounding area of the region-of-interest (ROI) can be more accurately reconstructed. Furthermore, a method is introduced to adaptively choose the step length, making the algorithm more feasible and easier for applications.

Keywords: Spectral CT, image reconstruction, anti-noising, adaptive step length, reference image

I. Introduction

Since Godfrey Hounsfield invented the first prototype in 1971, x-ray computed tomography (CT) has been extensively employed in many fields. Particularly, as a medical imaging modality, CT plays an important role in clinical and preclinical applications such as oncology tumor staging, acute stroke analysis, radiation therapy planning, etc. In past decades, both hardware and software of CT systems have undergone impressive developments [13]. Among all the available nondestructive medical imaging modalities, x-ray CT possesses the highest spatial and temporal resolution. However, it still cannot satisfy the requirements of many practical applications due to its inherent drawbacks. The first drawback is low contrast resolution. In reconstructed images from a traditional CT scanner, pixel values are the only elements to convey image information. Because the difference between the attenuation coefficients of soft tissue and blood/water is very small, it will results in low contrast resolution images which carry less information. In fact, the x-ray spectra for clinical/preclinical applications are polychromatic and the x-ray attenuation coefficient of an object depends on the photon energy. Because a conventional CT detector integrates photon energy throughout the whole spectrum, it doesn’t reflect the energy dependent information. This makes it difficult to distinguish different materials, especially for the case where the materials have similar attenuation coefficients in the reconstructed images. Another drawback of the conventional CT is the beam hardening artifact [4]. This is because conventional CT reconstruction theory assumes a linear integral projection model which is only valid for monochromatic x-ray. The beam hardening artifact is caused by the mismatch between the monochromatic projection model and the realistic polychromatic projections.

To overcome the aforementioned drawbacks, the dual-energy (DE) CT was proposed [59]. One of the most important merits of the DECT is that it takes advantage of spectrum information. There are three main methods to realize a DECT system. The well-known one is to use two x-ray sources with different voltages [10]. The second is to use one x-ray source combining the rapid voltage switching technique [11]. The third is to apply a dual-layer-based single-source device to the DECT [12]. Because most of the DECT also employs the energy-integrating detectors, there exists a significant spectral overlap between the two energy channels. This limits the applications of material decomposition. Although two x-ray sources have different energy spectra, the energy resolution of a DECT system is still limited. It was reported that the most prominent spectral feature of certain material is the position of its K-edge [13]. However, for those materials whose K-edge positions are adjacent, a DECT system cannot distinguish them efficiently. This implies that the DECT hasn’t fully exploited the spectrum information, and there is lots of room for improvements.

To make full use of the spectrum information, the concept of spectral CT (also known as color CT, energy sensitive CT, energy selective CT, multi-energy CT) was proposed by extending the idea of DECT along the energy dimension. Currently, the state-of-the-art photon counting detector based spectral CT systems can divide a whole x-ray spectrum into several channels [1421]. Although the photon counting detector cannot completely eliminate the overlaps between channels due to charge sharing and fluorescent x-ray escape [2225], it can significantly improve the energy resolution. By carefully setting the appropriate detector thresholds, abrupt increases of attenuation coefficients for certain materials, which are caused by the corresponding K-edges, could be observed. Consequently, more precise material decomposition results can be achieved. The photon counting detector is also characterized by high signal-noise-ratio (SNR) because the thresholds corresponding to each channel can block low-energy noise. Furthermore, it has been demonstrated that the spectral CT system has higher dose efficiency than the conventional one [26]. These merits have made the spectral CT a hot field in recent years, and inspire many research topics on different aspects of the spectral CT and its applications.

Although the spectral CT techniques have made a big progress in recent years, there still exist many technical problems that need to be conquered. As aforementioned, a photon counting detector can divide the received photons in a polychromatic spectrum into a number of channels with appropriate signal processing steps. Thus, if we assume the same amount of total photons as that for a corresponding energy-integrating detector, the amount of photons in each channel will be smaller. According to the well-accepted Poisson noise model, this implies higher noise in the projections of each channel and the reconstructed image quality will be degraded. To ensure the same image quality like the images reconstructed from projections of the energy-integrating detector, we can increase photon numbers in each individual channel. However, this will increase the total photon number and result in excessive radiation dose which is harmful to the human body. Another problem is concerned with the spectral detectors. Currently, the state-of-the-art photon counting detectors are premature, so frangibility and expensiveness are two prime weaknesses. When a detector covers the whole object, a saturated x-ray source can make the detector edge lose linearity or even damage it. In the meantime, a detector that is big enough to cover the whole imaging object is very expensive. Besides, due to weak crystal/electronics capability of dealing with clinical level X-ray flux, low radiation dose has to be applied in practice. This can seriously degrade the reconstructed image quality. All these drawbacks limit the applications of spectral CT.

For practical clinical applications, the aforementioned frangibility and expensiveness for spectral CT can be partially overcome by the interior reconstruction technique, since we are always interested in smaller internal ROIs for diagnosis. It has been proved that an interior region can be theoretically and exactly reconstructed from truncated local projections if certain prior information is available. Till now, two popular interior reconstruction theories and algorithms have been developed. One is based on the truncated Hilbert transform [2732] assuming a known sub region, and the other is based on the compressive sensing (CS) theory[33, 34] that assumes a piecewise constant/polynomial imaging model [3538]. Nevertheless, for spectral CT, one still faces the noise problem in each individual channel if we keep the same dose level and employ the conventional CT reconstruction methods. This motivates us to develop a new spectral CT reconstruction algorithm to improve the image quality in each individual channel, and apply it to spectral interior reconstruction.

By comparing the reconstructed images from different channels for the same object, we notice some common features among them e.g. the relative locations between different materials are invariable. This reveals that the reconstructed images from different channels are highly correlated. If a high quality image of the object is known, it can be used to serve as a reference image (RI) to provide prior information [3941]. By incorporating the prior information, the low-dose reconstructions of each individual channel can be significantly improved. In the high-dimensional space, both the RI and target image (TI) can be treated as vectors. One natural way to incorporate the regularization of correlation is to minimize the space angle between the RI and TI. To make sure that each modulated direction can reduce the space angle, the steepest descent method (SDM) can be applied. Furthermore, an adaptive step length (ASL) [42] is introduced to optimize the results. It makes the algorithm suitable for all kinds of images, without adjusting the parameters. The proposed algorithm has two major advantages. First, it maintains the primary structures of the TIs even when the noise is extremely strong. This feature is very useful when the dose is ultra-low. Second, it can help suppress the constant shift of intensity especially for interior reconstruction.

The rest of this paper is organized as follows. In Section II, the proposed algorithm and implementation details will be presented. In Section III, both numerical simulations and preclinical experiments will be performed to evaluate the proposed algorithm. Finally, some related issues will be discussed and conclusions will be made in Section IV.

II. Method

As we pointed out earlier, although the reconstructed images from each channel are different in spectral CT, they are highly correlated. Thus, if a high quality image is available, it can be treated as the RI to reconstruct the relevant images from low-dose projections or truncated projections. Here, the RI can be reconstructed from a set of pre-existing normal dose projections. If there are no pre-existing projections, we can always combine all the photons in different channels to synthesize a set of projections of energy-integrating detector and reconstruct a high-quality RI.

2.1. Imaging Model

Usually, a 2D digital image can be described by a W × H matrix in which each element is treated as a pixel with a constant intensity value fw,h, where the indices w =1,2 …,W and h =1,2, …, H are integers, and W and H are width and height of the image, respectively. In the context of CT reconstruction, it is convenient to rearrange a 2D digital image f = (fw,h) [set membership] RW × RH into a vector f⃑ = [f1, (...), fn, (...) fN]T [set membership] RN where N =W × H and “T” represents the transpose operation. In the vector f⃑, each element fn is equivalent to fw,h with the following relationship

n=(w1)×H+h;h=1,2,,H;w=1,2,,W.
(1)

In this paper, a 2D image will be expressed in both matrix and vector formats, when it is appropriate to do so without any confusion.

Another important component of the CT system is the projection dataset [g with right harpoon above] = [g1, (...), gm, (...), gM]T[set membership] RM, where M represents the total amount of ray integrations that equals the product of the number of projections and the number of detector elements, and gm denotes measured line integration along the mth x-ray path. In the projection procedure, the contribution of nth pixel fn to the ray path can be denoted as am,n. As a result, we have a system matrix A = (am,n) [set membership] RM × RN, and the CT system can be modelled as the following linear matrix equation:

g=Af+e,
(2)

where [e with right harpoon above] [set membership] RM is measured noise. Obviously, in Eq.(2), the matrix A represents the forward projection, and the matrix AT can represent the back-projection. There are many methods available to compute the coefficient am,n based on different projection models. For example, according to the well-known ray-driven model [4345], am,n can be calculated as the intersection length between the mth x-ray path and nth pixel area. The mth measured datum gm can be calculated by summing up the product of am,n and fn along the x-ray path.

2.2. OS-SART Algorithm

To recover the original image vector f⃑, one can make full use of the system matrix A and the measured projection dataset [g with right harpoon above] based on Eq.(2). Actually, we can minimize the noise energy to estimate f⃑. This can be implemented by minimizing the discrepancy function Δ(f):

Δ(f)=gAf2,
(3)

where ‖·‖2 represents the square of L2 norm. Generally speaking, there is no analytic solutions for the minimization of Eq.(3), and we need to employ an iterative approach. There are many iterative algorithms for image reconstruction, such as algebraic reconstruction technique (ART) [46] and simultaneous algebraic reconstruction technique (SART) [47]. Generally speaking, ART converges faster for noise-free projections. However, its result bounces around for noisy projections with a compromised image quality. In this work, SART is selected to compromise the performances and convergence. According to[47], the SART algorithm can be expressed as

fnk=fnk1+λk1a+nm=1Mam,nam+(gmAmfk1),
(4)

where am+ and a+n are respectively the sum of the mth row and the nth column of the system matrix A, Am represents the mth row of system matrix, k is the iteration index, and λk [set membership] (0, 2) is a free relaxation parameter.

To accelerate the convergence of SART, one can employ the ordered-subset (OS) technique to form an OS-SART algorithm [49]. The OS-SART can be formulated by partitioning the index set B = {1,2, …, M} into T nonempty disjoint subsets Bt={m1t,m2t,,mM(t)t} such that

B={1,2,,M}=0tT1Bt.
(5)

Then, Eq.(4) can be rewritten as

fnk=fnk1+λkmBtam,nmBtam,ngmAmfk1am+,t=k mod T.
(6)

Besides that, the OS-SART algorithm can be further accelerated by the weighting technique for the fast iterative shrinkage-thresholding algorithm (FISTA) [50], which uses the last two outputs as a feedback to refine the current output.

2.3. Regularization Function

To solve the ill-posed problems uniquely, especially for spectral CT, an additional constraint should be introduced to the object function Eq.(3). Therefore, the ill-posed image reconstruction problem can be expressed as

minfgAf2+Φ(f),
(7)

where Φ (f⃑) represents a regularization function. Eq.(7) can be solved by using the well-known alternative minimization scheme which can be implemented in two main steps. In the first step (fidelity update), the fidelity term is enforced to reconstruct an intermediate image by using Eq. (6). In the second step (regularization update), the intermediate image from the first step will be utilized as an input to minimize the regularization function. The two consecutive steps are treated as one iteration.

For the spectral CT reconstruction problem, the reconstructed images from projections of different channels are highly correlated. For a given high-quality RI [u with right harpoon above], we can define a regularization function for each individual channel image f⃑

Φ˜(f)=cos(θ)=fTufu,
(8)

where f⃑ represents the reconstructed TI and [u with right harpoon above] is the RI. Because both TI and RI are treated as vectors in high-dimensional space, Eq.(8) can be interpreted as the minus cosine of the space angle (SA) θ between the two vectors. Furthermore, if an image is treated as a random variable with a series of N measurements, Eq.(8) can be interpreted as the minus sample correlation coefficient. From Eq. (8), one can see that if images f⃑ and [u with right harpoon above] are similar, the value of Eq.(8) will approach −1. Therefore, updating the TI f⃑ ensures a smaller SA and leads to a better image. To improve the sensitivity, we modify Eq.(8) as follows

Φ(f)=(ff¯)T(uu¯)ff¯uu¯,
(9)

where f and ū are the means of f and ū, respectively. The means f and ū are also called direct current (DC) components, and (f⃑f) and ([u with right harpoon above]ū) are referred to as alternating current (AC) components. Because the main common features hidden in different channels come from the ACs, Eq.(9) should be an excellent regularization function for the correlation, and we use Eq.(9) to serve as the regularization function Φ (f⃑) in Eq.(7).

2.4. Reconstruction Algorithm

To minimize the regularization term, the SDM is employed. During the regularization update process, the updates to the TI can be described by

fk+1=fOSSARTk+αk,*pk,
(10)

where αk,* and [p with right harpoon]k represent the step length and a normalized direction vector to update the estimated image in the kth iteration. The direction vector [p with right harpoon]k can be computed as

p=Φ(f)Φ(f).
(11)

Let [f with tilde] = f⃑f, and ũ = [u with right harpoon above] − ū. Each component of the gradient direction [nabla]Φ(f⃑) can be determined as

fn(Φ(f))=u˜nf˜u˜f˜nf˜Tu˜u˜f˜f˜2u˜2,
(12)

where [f with tilde]n and ũn are the nth elements of vectors [f with tilde] and ũ.

In this work, the step length in Eq. (10) plays an important role to control the convergence. In practice, it is very difficult to choose the optimal parameters for an iterative algorithm, and the parameters are usually empirically selected based on extensive experiments. To improve the computing efficiency, we propose to employ one adaptive step length method. For this purpose, the Wolfe conditions based method is chosen [42]. It is a typical inexact line search condition described as follows

Φ(fk+αk,*pk)Φ(fk)+c1αk,*ΦkTpk,
(13)

Φ(fk+αk,*pk)Tpkc2ΦkTpk,
(14)

where c1 and c2 are constants satisfying 0 < c1 < c2 < 1, and [nabla] Φk is the short form of [nabla]Φ(f⃑k). Formulas (13) and (14) actually offer two constraints, which are named Armijo condition and curvature condition, respectively. In formula (13), the left side should be considered as a function with respect to the variable α such that ϕ(α) = Φ(f⃑ + α[p with right harpoon]), and the right side is a monotonously decreasing function, which can be denoted by [var phi](α). Because the range of c1 is from 0 to 1, for certain small positive value of α, [var phi](α) must be greater than ϕ(α). Therefore, condition (13) guarantees that the acceptable step length α must decrease ϕ(α) in a definite extent. However, the Armijo condition is not sufficient to ensure reasonable progress, because it always accepts all the sufficiently small step lengths. To exclude those too short step lengths, the curvature condition has to be incorporated. In formula (14), the left side presents the gradient ϕ′(α), and right side is the product between the gradient ϕ′(0) and a constant factor c2. Because ϕ′(α) should not be smaller than c2 timesϕ′(0), it is not hard to understand why the curvature condition can help to retain sufficiently large step lengths. To rule out large step lengths, whose slope is positive, the curvature condition could be revised as

|Φ(fk+αk,*pk)Tpk|c2|ΦkTpk|.
(15)

Then, formulas (13) and (15) can be used to construct strong Wolfe conditions. To implement the strong Wolfe conditions, we first need to set an initial step length αk,0 for each inner iteration. Actually, this initial step length is always the same for all outer iterations. The pseudo-codes for the steepest decent method are summarized in Table 1.

Table 1
Pseudo-codes for gradient descent method with adaptive step lengths.

To keep the image details and suppress the effect of direct current component, a patch-wise processing method is adopted. We set the patch size (PS) as D = d × d. One patch works like a window, which can traverse the image pixel by pixel. When the ith pair of corresponding patches are extracted from RI and TI at the same position, they will be processed via the aforementioned strategy. In this way, each tiny structure in the RI could be used to guide the update of the reconstructed image. After a pair of image patches has been processed, the normalized resulting patch will be multiplied by the module value corresponding to the one in TI. This can eliminate the deviation of AC on magnitude. Then, the DC component needs to be added back. Because the patches have overlaps, each pixel may have multiple values in different patches. For the same pixel, all the values are uniformly weighted to generate the mean as the final result of the regularization term. The aforementioned procedure is referred to as adaptive space angle (AdSA) method.

An important issue for an iterative algorithm is how to terminate the iteration. For the OS-based algorithm, T iterations form a complete period. Therefore, we make a judgment for the termination of the program every T iterations. We first compute the module value of the difference dk [triangle, equals]f⃑kf⃑k−1‖ between the results from two successive iterations, and then the average value dτ can be calculated within τth complete period. Consequently, the convergent ratio of the program can be observed by comparing the adjacent dτand dτ−1 If the difference of dτ and dτ−1 is small enough, the algorithm will be terminated. In this paper, the threshold of the variation is set to be 1/2000 of the averaged attenuation coefficient of TI. More details of the algorithm are shown in Table 2.

Table 2
Pseudo-codes for the proposed algorithm.

III. Results

The proposed algorithm is implemented according to the pseudo-codes summarized in Tables 1 and and2.2. To demonstrate the feasibility and evaluate the proposed algorithm for spectral CT, extensive numerical simulations and preclinical experiments are performed. In the numerical simulation studies, the classical filtered backprojection (FBP) is first used to reconstruct images from noise-free projections as benchmarks, which are compared with the results from reconstruction by the proposed method for both global and locally truncated projections. In the preclinical experiment, a set of realistic micro-CT projections of a mouse are applied to evaluate the performance of the proposed method. Material decomposition is performed for further performance analysis based on the method proposed in [51].

3.1 Numerical Simulations

A 2D analytic phantom (see Figure 1) containing 17 disks is constructed. The biggest disk serves as background to mimic the soft tissue. The rest of 16 disks are used to mimic six different materials and different spatial resolutions and their parameters are listed in Table 3. To analytically generate projections of the 2D numerical phantom, we assume fan-beam geometry with an equidistant detector, and the detector element size is 0.08 mm. 512 and 256 detector elements are assumed to simulate a long and a short detector to collect global and truncated local projections, respectively. The distance between the x-ray source focal spot and the origin (the center of the phantom) is assumed to be 5cm. We also assume that the detector is perpendicular to and symmetric with the line passing the origin and the x-ray focal spot, and the distance between the focal spot and the detector is 10cm. 720 projections are uniformly acquired over a range of a full scan. To simulate a polychromatic x-ray source (GE Maxiray 125), 951 uniform samples (from 25KeV~120KeV) were extracted. There are four k-edges among the material attenuation coefficients around 33, 37.4, 50.2 and 80.7KeV for Iodine, Barium, Gadolinium and Gold, respectively. To optimize the channels for the best image quality, it is better to make the total photons in each channel as equal as possible, and the difference of attenuation coefficients between different channels as great as possible. According to the aforementioned principles, eight spectral channels are set to be 25~32KeV, 32~37KeV, 37~43KeV, 43~50KeV, 50~58KeV, 58~65KeV, 65~80KeV, and 80~120KeV. According to the spectrum, the percentage of incident photons at each energy level is determined. After the projections at all energy levels are collected, exponential, weighted summation and logarithmic operations are performed to produce linearized projections for the polychromatic x-ray spectrum. To further evaluate the performance of the proposed algorithm, different photon numbers are assumed. Here, we mainly use 2 × 104 and 105 incident photons to simulate two different dose levels.

Figure 1
Sketch map of the 2D phantom, in which different colors represent different materials.
Table 3
Parameters for 17 disks of the 2D phantom.

3.1.1 Global Reconstruction

The reconstructed images are set as 256 × 256 matrix to cover a 2.0 × cm2 region which is sufficient to cover the largest disk of the phantom. The coefficients c1 and c2 are respectively set to 0.0001 and 0.01, and the PS is set as 8 × 8. For qualitative and quantitative comparison, the FBP and a soft-threshold filtering based TV (StTV) minimization [36] method are also implemented.

Figure 2 shows the global reconstruction results from noisy projections. Poisson noise is superimposed to the raw projections by assuming 2 × 104 incident photons for each spectral detector element. It can be seen that the FBP images have heavy noise due to the ultra-low photon number. The four smallest disks cannot be visually identified, especially in the 8th channel. The restored information is very limited. The StTV results are in the same boat: about half of disks in the 8th channel disappear into the background. The images also become seriously blurred with the increase of the photon energy. However, the AdSA reconstructs almost all the details. The material edges are clear, and even the smallest disk can be recognized in some channels. Although the AdSA cannot reconstruct perfect results from the noisy data (e.g. some bright-spot-like artifacts can be seen on the disks), the AdSA results are still the best among the three reconstruction methods. From the material decomposition results, one can see that FBP fails to restore part of the spectrum information for the simulated phantom. Colors used to represent different materials are scattered all around the result. Particularly, the 1st and 6th disks are recognized as similar color. For the StTV, the material edges always reflect some inconsistent colors, which are mainly caused by noisy data. However, the AdSA gives a much better result that is consistent with the ground truth.

Fig. 2
Spectral reconstruction results from noisy projections assuming 2× 104 incident photons per detector element. From top to bottom rows, the images are reconstructed by the StTV, AdSA and FBP, respectively. The images in the first three columns, ...

In order to quantitatively evaluate the performance of the proposed method, the root-mean-square error (RMSE) and structural similarity (SSIM) are used as assessment metrics. The range of RMSE is from 0 to infinite. The smaller the RMSE value is, the better the image quality is. The range of SSIM index is from 0 to 1, where 0 represents totally irrelevant and 1 means exactly the same structure [52]. To suppress the fluctuation in flat region, the noise-free FBP reconstruction is filtered by the soft-threshold filtering method to serve as the reference image for RMSE computing. The RMSEs of all the aforementioned results are summarized in Table 4. Obviously, the performance of FBP is the worst for noisy projections. For the noise-free projections, the RMSEs of StTV and AdSA are almost the same except for the 1st and 3rd channel. This is because the numerical phantom assumes the piecewise constant image model which shows priority by using TV minimization method, and the FBP results may have some noise in the flat regions due to discrete errors. However, the AdSA outperforms the StTV for noisy projections. The greater the photon energy is, the more obvious the advantage of the AdSA is. The results in Table 5 are the corresponding values of SSIMs, from which similar conclusions like of the RMSEs can be made. It reveals that the AdSA has strong capability to suppress noise.

Table 4
Quantitative evaluation results of global reconstruction in term of RMSE (unit: cm−1).
Table 5
Quantitative evaluation results of global reconstruction in term of SSIM.

Furthermore, to evaluate the spectral accuracy of reconstructed results, relative biases of different material regions are calculated (Figure 3). Here, the smoothed noise-free FBP results are set to be the reference and the selected regions (100 × 100 pixels) are marked in figure 2. One can see that the AdSA offers the best results. The relative biases of FBP are obvious and much greater than its other two competitors. Although the performance of StTV is similar to that of AdSA in several channels, the relative biases of the StTV fluctuate heavily from channel to channel, which are mainly caused by the noisy data. On the other hand, the performance of AdSA is steadier for all channels and materials. The relative biases for the AdSA are always limited to within 10−3, which is negligible in practice. One can see that the proposed algorithm ensures the spectral accuracy of reconstructions, even under noisy conditions.

Fig. 3
Relative biases of different material regions, which are marked with red squares in figure 2. From left to right, the plots are for the soft tissue, 12.4%Ca + 87.6%water and 1.6%gold + 98.4%water, respectively.

3.1.2 Interior Reconstruction

In interior reconstruction experiments, only the central 256 detector elements are used to collect projections. The reconstructed image size and the corresponding region are the same as the global reconstruction. For AdSA, the RIs are the corresponding global images reconstructed by the StTV.

The aforementioned experiments for global reconstruction are repeated for the interior reconstruction with locally truncated noise-free and noisy projections (105 and 2 × 104 incident photons). Figure 4 are the results from noise-free truncated projections where the ROIs are marked by circles. Although the ROIs are similar between the StTV and AdSA images, the AdSA significantly outperforms StTV for the regions outside of the ROIs. The material decomposition of AdSA looks similar to the global one, while StTV only conserves the color map inside the ROIs and part of the outside region. With the regularization of the RI, the AdSA can make full use of the truncated projections to obtain better image quality not only inside the ROIs but also outside them. This is exactly the color diffusion phenomenon, which we previously reported [38]. To further demonstrate the difference between these two methods, Figure 5 shows the corresponding profiles along the line indicated in Figure 4. In Figure 5, we find that although the curves of StTV have similar shape with standard ones, there are some errors in the material edge areas.

Fig. 4
Spectral reconstruction results from noise-free truncated projections. The top and bottom rows are for the images reconstructed by the StTV and AdSA, respectively. The images in the first three columns, from left to right, are for the 1st, 4th and 8th ...
Fig. 5
Profiles along the line in the top-left of figure 4. From left to right, the plots are for the 1st, 4th and 8th channels in Figure 4, respectively.

Representative images reconstructed from projections with 2 × 104 incident photons per detector element are shown in Figure 6. Although the image quality of these two algorithms become worse, the images reconstructed by AdSA are still very close to the global reconstructions, and the edges in the AdSA reconstructions are much clearer than its counterparts in the StTV. For the material decomposition results, the AdSA shows a consistent performance. Figure 7 shows the corresponding profiles. Although strong noises are superimposed into the raw projections, the AdSA can still recover high quality ROI images with fine details especially for the high frequency components. However, the StTV appears to have certain constant shift in the whole ROIs. This demonstrates that the AdSA is robust with respect to noise.

Fig. 6
Same as Figure 4, but from noisy projections assuming 2×104 incident photons per detector element.
Fig. 7
Same as Figure 5, but for the images in figure 6.

The RMSE (see Table 6) and SSIM (see Table 7) are also used to quantitatively evaluate the AdSA and StTV for interior reconstruction. Because only the ROI can be theoretically and exactly reconstructed from truncated projections, we calculate the indices only using the pixels inside the ROIs. From Table 6, one can see that AdSA significantly outperforms StTV, especially for the case with heavy noise. This is due to the additional prior information from RI. The SSIM values in Table 7 reveal that the AdSA is not sensitive to noise. When the photon energy increases, the corresponding SSIM values increase to close 1. On the other hand, the performance of StTV is worse and worse, no matter what the photon energy is increased or the incident photon number is reduced.

Table 6
Quantitative evaluation results of interior reconstruction in term of RMSE (unit: cm−1).
Table 7
Quantitative evaluation results of interior reconstruction in term of SSIM.

3.2 Preclinical Mouse Study

To demonstrate the application of the proposed algorithm, realistic preclinical projections are adopted. The projections are provided by the medipix all resolution system (MARS) team in New Zealand, and they were acquired using a MARS micro-CT with a Medipix 3 CdTe detector. It is equipped with a network of charge summing circuits, that communicate between adjacent 2 × 2 pixels, for detecting coincidences and reconstructing charges[22], which significantly helps to reduce the effects such as charge sharing and K-fluorescence especially when the detector size is very small. In this work, the detector row contains 1024 elements corresponding to a unit length of 55 µm. This covers a field of view (FOV) of 34.89 mm in diameter. The mouse was euthanized, and injected with 0.2 mL of 15 nm Aurovist II GNP (Nanoparticles; Yaphank, NY) into the tail vein. A Source-Ray SB-120-400 (produced by Source-Ray Inc., New York), which has a minimum focal spot of 75 µm, was applied to illuminate the mouse with 120 kVp and 175 mA. The detector and the source spot were positioned on the opposite directions with a distance of 255mm between them, and the center of the mouse was 158mm away from the source spot. Although a flat panel detector was employed with cone-beam geometry, we only extracted the central row in terms of fan-beam geometry. 371 projections of 13 energy bins were uniformly acquired over 360 degrees. In order to improve the quality of projection data, two adjacent detector elements were merged to reduce the data noise. As a result, 512 samples were obtained per projection. After that, a wavelet and Fourier based algorithm [53] was applied to get rid of ring artifacts from the data. Because the first channel was allotted with the lowest threshold and contains less noise, it can be treated as the gray scale projection. For the mouse, we set the size of the reconstructed image to 512 × 512 pixels covering a region of 18.41 × 18.41 mm2. The gray scale channel was first reconstructed by the FBP to obtain a preliminary image. Then the reconstructed image was de-noised by the BM3D [54] to obtain the RI. The rest of the 12 channels were used for reconstruction and comparison.

For the global experiment, both the FBP and StTV are implemented to compare with the proposed AdSA. The reconstructed images of the 2nd and 13th channels are shown in Figure 8. Two ellipses and one red arrow indicate representative features in the reconstructed images. Although all the three methods successfully reconstruct the crevice in channel 2 which is indicated by the arrow due to a relatively low noise level, for the same structure in channel 13, only the one reconstructed by the AdSA can be visually identified. Some tiny bony structures show the similar phenomenon. The most representative is the one at the far right of the upper ellipse, whose hole can hardly be seen in the FBP and StTV reconstructions. This kind of blurry result is derived not only from heavy noises but also from the high frequency structure. Although the main structure of the mouse can be recognized by all three methods in the material decomposition results, the FBP result is too noisy, and the AdSA reconstructs more clear edges (e.g. the narrow central crevice, which is not clear for StTV).

Fig. 8
Spectral reconstruction results of the mouse study from global projections. From top to bottom rows, the images are reconstructed by the StTV, AdSA and FBP, respectively. The 1st and 2nd column images are respectively for the 2nd and 13th channels, and ...

The central 128 samples (256 detector elements) of the combined 512 version in the previous global reconstruction are manually extracted to simulate truncated local projections. The corresponding local reconstruction results are shown in Figure 9, where the ROIs are marked by circles. The long and narrow crevice, the most obvious feature at the center of the image, is successfully reconstructed by the AdSA. However, it cannot be visually recognized in the images reconstructed by the StTV. For the rest of the structural features, few are clearly restored by the StTV regardless of whether they are inside or outside the ROI. The color image of AdSA shows a comparable performance to that of the global reconstruction, and most of the main features are conserved and successfully identified. This further confirms AdSA as a robust interior reconstruction algorithm.

Fig. 9
Spectral reconstruction results of the mouse study from truncated projections. The top and bottom images are reconstructed by the StTV and AdSA, respectively. The 1st and 2nd column images are for the 2nd and 13th channels, respectively, and the display ...

IV. Discussion and Conclusion

From the aforementioned numerical simulations and preclinical application, one can see that the advantages of spectral CT over the conventional CT are apparent. In the numerical phantom, the materials in the 2nd – 6th disks have similar attenuation curves which imply a very challenging task to distinguish them using the conventional CT. However, by employing the photon counting detector, these materials can be reconstructed in different channels and mapped into different colors (see Figure 2). which is much easier to distinguish. Therefore, it is no wonder that spectral CT has attracted a great deal of attention in recent years, and is a major direction of CT development. Consequently, the proposed algorithm has arrived in the right time to promote spectral CT.

In the implementation of the proposed AdSA, there are several parameters that need to be selected. The most primary one is the PS. In our experiments, the PS is set as 8 × 8 pixels. Actually, extensive numerical simulations have been conducted to optimize the PS. We found that, when the PS is small, the AdSA can accurately capture the small-sized features which correspond to high frequency information. However, because most of the noises are also in high frequency, it is hard for the AdSA to get rid of noises with a small PS. On the other hand, when the PS is large, an apparent DC component will be introduced even though the AdSA can effectively suppress noises. For example, when the PS is 64 × 64, the resulting images are almost the same as the RI for all the channels. This is because most of the DCs are kept in the large patches, although they have been subtracted by the mean values. Another weakness brought by large PS is the exponential growth of the computational cost, which is a serious issue. In summary, the PS has to be selected to compromise the aforementioned problems, resulting in the size 8×8 for practical applications.

Besides the patch size, the selection of the RI also significantly affects the performance of the AdSA. In the numerical simulations, the StTV is selected to reconstruct the RI from projections of all the photon energies, because the phantom is piecewise constant, and it precisely matches the image model of TV minimization. In the preclinical experiments, an algorithm called BM3D is used to suppress noise. In fact, finding a good RI is one of the key issues for the AdSA, and the dependency of RI is also a weakness of the proposed method. An additional experiment is conducted to support this point. When the RI is selected as the noise-free reconstruction to regularize the noisy reconstructions with 2 × 104 incident photons, the reconstructed image quality is much better than its counterpart using the RI reconstructed from the projection with 2 × 104 incident photons. Nevertheless, no matter what image is selected as the RI, the spectral images reconstructed by the AdSA are much better than the ones reconstructed merely from the data in one channel. This is because the noise level in any channel of the spectral CT is much greater than the noise level in whole energy integrating data.

Furthermore, the convergence of the proposed algorithm is highly concerned. This can be explained by the physical meaning. Because the reconstructed images from all the channels are highly correlated, the most intuitionistic relevant information comes from their edges and textures. These features are relatively invariable even when the attenuation coefficients of different materials vary in their own patterns. If the selected reference image contain some fine features, by using the constraint Eq.(9), the features will be preserved well in each channel. Actually, the magnitudes and direct current components of image patches are only determined by the fidelity term. Thus, the convergence of Eq.(9) is guaranteed.

The aforementioned numerical simulations and preclinical application not only have validated the feasibility, but also have demonstrated several merits of the proposed algorithm. In comparison to the FBP and StTV, the AdSA has a strong capability for de-noising, which can be seen by comparing the numerical simulations without noise and with 2 × 104 incident photons. When the amount of incident photons is reduced to an ultra-low level, the AdSA still could reconstruct images with very similar image quality to the noise-free case. The AdSA is good at dealing with low contrast images. Both the quantitative RMSE and SSIM measurements have shown that the higher the photon energy is, the greater the advantage of the AdSA is over the competitors. Because higher photon energies make lower contrast images, this implies that the AdSA works better for low contrast images. When the patch size is appropriately selected, high frequency information can also be restored well. This point can be confirmed by both the reconstructed smallest disk in the numerical simulations and the long and narrow crevice in the preclinical mouse study. With the truncated data, the AdSA can almost accurately reconstruct the whole spectral image if a global RI is available. This color diffusion phenomenon is attributed to the global prior information provided by the RI and the limited-angle projections for the surrounding regions of the ROI. Once an RI has been selected, no more parameters need to be selected. This is an important merit for the proposed iterative algorithm. In this work, the SLs of patches need to be determined thousands of times in one iteration. The introduced ASL successfully skips the complex procedure, and two relevant parameters (c1 and c2) are always fixed. This demonstrates the AdSA is more convenient and flexible for all conditions.

The key of the proposed AdSA method is to employ the RI to improve the spectral CT image quality. Actually, the method of making use of RI to improve the image quality for spectral CT has been employed in several literatures, but in different forms. In [38], Xu et al. adopted the gray scale image reconstructed from the global projections as an initial guess of the reconstructed images in a different channel. However, we treat the gray scale image as an RI to improve the image quality throughout the whole reconstruction procedure. In [55], Zhang et al. proposed an algorithm named TV-SM-LM. While the TV-SM-LM calculates the L2-norm value of the differences between different channels, the proposed AdSA measures the similarities between TIs and the RIs. While the RI of the AdSA is fixed, it is dynamic in the TV-SM-LM. Therefore, the proposed AdSA is different from the existing methods.

In conclusion, an algorithm that incorporates the information of an RI is proposed for spectral image reconstruction. Both the extensive numerical simulations and preclinical application demonstrate that the proposed algorithm has strong capability for anti-noising. It is also robust and stable for interior reconstruction. This can help to advance spectral CT by improving the accuracy of material decomposition from low-dose projections. Furthermore, the strategy to adjust the SL makes the proposed algorithm more competitive.

Acknowledgments

This work was supported in part by NIH/NIBIB U01 grant EB017140 and R21 grant EB019074, in part by NSF CAREER award No. 1540898 and NSF grant No. 1619550. The authors would like to thank the MARS team in New Zealand for providing the mice datasets. The authors are also grateful to Mr. Ibrahim Mkusa at UMass Lowell for his proof-reading.

References

1. Kalender WA. X-ray computed tomography. Physics in medicine and biology. 2006;51:R29. [PubMed]
2. Wang G, Yu H, De Man B. An outlook on x-ray CT research and development. Medical physics. 2008;35:1051–1064. [PubMed]
3. Shefer E, Altman A, Behling R, Goshen R, Gregorian L, Roterman Y. State of the art of CT detectors and sources: a literature review. Current Radiology Reports. 2013;1:76–91.
4. Brooks RA, Di Chiro G. Beam hardening in x-ray reconstructive tomography. Physics in medicine and biology. 1976;21:390. [PubMed]
5. Alvarez RE, Macovski A. Energy-selective reconstructions in x-ray computerised tomography. Physics in medicine and biology. 1976;21:733. [PubMed]
6. Macovski A, Alvarez R, Chan J-H, Stonestrom J, Zatz L. Energy dependent reconstruction in X-ray computerized tomography. Computers in biology and medicine. 1976;6:325–336. [PubMed]
7. Alvarez RE, Macovski A. X-ray spectral decomposition imaging system. ed: Google Patents. 1977
8. Marshall WH, Jr, Alvarez RE, Macovski A. Initial results with prereconstruction dual-energy computed tomography (PREDECT) Radiology. 1981;140:421–430. [PubMed]
9. Stonestrom JP, Alvarez RE, Macovski A. A framework for spectral artifact corrections in X-ray CT. IEEE Transactions on Biomedical Engineering. 1981;2:128–141. [PubMed]
10. Graser A, Johnson TR, Chandarana H, Macari M. Dual energy CT: preliminary observations and potential clinical applications in the abdomen. European radiology. 2009;19:13–23. [PubMed]
11. Zou Y, Silver MD. Analysis of fast kV-switching in dual energy CT using a pre-reconstruction decomposition technique. Medical Imaging. 2008:691313-691313-12.
12. Carmi R, Naveh G, Altman A. Material separation with dual-layer CT; Nuclear Science Symposium Conference Record, 2005 IEEE; 2005. p. 3.
13. Roessl E, Proksa R. K-edge imaging in x-ray computed tomography using multi-bin photon counting detectors. Physics in medicine and biology. 2007;52:4679. [PubMed]
14. Jakůbek J. Semiconductor pixel detectors and their applications in life sciences. Journal of Instrumentation. 2009;4:P03013.
15. Fredenberg E, Lundqvist M, Cederström B, Åslund M, Danielsson M. Energy resolution of a photon-counting silicon strip detector. Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment. 2010;613:156–162.
16. Fredenberg E, Hemmendorff M, Cederström B, Åslund M, Danielsson M. Contrast-enhanced spectral mammography with a photon-counting detector. Medical physics. 2010;37:2017–2029. [PubMed]
17. Kappler S, Hannemann T, Kraft E, Kreisler B, Niederloehner D, Stierstorfer K. First results from a hybrid prototype CT scanner for exploring benefits of quantum-counting in clinical CT. SPIE Medical Imaging. 2012:83130X-83130X-11.
18. Kappler S, Glasser F, Janssen S, Kraft E, Reinwand M. A research prototype system for quantum-counting clinical CT. SPIE Medical Imaging. 2010:76221Z-76221Z-6.
19. Ballabriga R, Campbell M, Heijne E, Llopart X, Tlustos L. The Medipix3 prototype, a pixel readout chip working in single photon counting mode with improved spectrometric performance; Nuclear Science Symposium Conference Record, 2006. IEEE; 2006. pp. 3557–3561.
20. Steadman R, Herrmann C, Mülhens O, Maeding DG, Colley J, Firlit T. ChromAIX: a high-rate energy-resolving photon-counting ASIC for spectal computed tomography. SPIE Medical Imaging. 2010:762220-762220-8.
21. Schlomka J, Roessl E, Dorscheid R, Dill S, Martens G, Istel T. Experimental feasibility of multi-energy photon-counting K-edge imaging in pre-clinical computed tomography. Physics in medicine and biology. 2008;53:4031. [PubMed]
22. Taguchi K, Iwanczyk JS. Vision 20/20: Single photon counting x-ray detectors in medical imaging. Medical physics. 2013;40:100901. [PubMed]
23. Aillon EGD, Gentet MC, Montémont G, Rustique J, Verger L. Simulation and experimental results on monolithic CdZnTe gamma-ray detectors. Nuclear Science, IEEE Transactions on. 2005;52:3096–3102.
24. Guerra P, Santos A, Darambara D. Development of a simplified simulation model for performance characterization of a pixellated CdZnTe multimodality imaging system. Physics in medicine and biology. 2008;53:1099. [PubMed]
25. Xu C, Danielsson M, Bornefalk H. Evaluation of energy loss and charge sharing in cadmium telluride detectors for photon-counting computed tomography. Nuclear Science, IEEE Transactions on. 2011;58:614–625.
26. Tapiovaara MJ, Wagner R. SNR and DQE analysis of broad spectrum x-ray imaging. Physics in Medicine and Biology. 1985;30:519.
27. Ye Y, Yu H, Wei Y, Wang G. A general local reconstruction approach based on a truncated Hilbert transform. Journal of Biomedical Imaging. 2007;2007:2–2. [PMC free article] [PubMed]
28. Courdurier M, Noo F, Defrise M, Kudo H. Solving the interior problem of computed tomography using a priori knowledge. Inverse problems. 2008;24:065001. [PMC free article] [PubMed]
29. Kudo H, Courdurier M, Noo F, Defrise M. Tiny a priori knowledge solves the interior problem in computed tomography. Physics in medicine and biology. 2008;53:2207. [PubMed]
30. Ye Y, Yu H, Wang G. Exact interior reconstruction with cone-beam CT. International journal of biomedical imaging. 2008;2007 [PMC free article] [PubMed]
31. Ye Y, Yu H, Wang G. Exact interior reconstruction from truncated limited-angle projection data. Journal of Biomedical Imaging. 2008;2008:5. [PMC free article] [PubMed]
32. Yu H, Ye Y, Wang G. Interior reconstruction using the truncated Hilbert transform via singular value decomposition. Journal of X-ray science and technology. 2008;16:243. [PMC free article] [PubMed]
33. Yu H, Yang J, Jiang M, Wang G. Supplemental analysis on compressed sensing based interior tomography. Physics in medicine and biology. 2009;54:N425. [PMC free article] [PubMed]
34. Yang J, Yu H, Jiang M, Wang G. High-order total variation minimization for interior tomography. Inverse problems. 2010;26:035013. [PMC free article] [PubMed]
35. Yu H, Wang G. Compressed sensing based interior tomography. Physics in medicine and biology. 2009;54:2791. [PMC free article] [PubMed]
36. Yu H, Wang G. A soft-threshold filtering approach for reconstruction from a limited number of projections. Physics in medicine and biology. 2010;55:3905. [PMC free article] [PubMed]
37. Xu Q, Mou X, Wang G, Sieren J, Hoffman E, Yu H. Statistical interior tomography. Medical Imaging, IEEE Transactions on. 2011;30:1116–1128. [PMC free article] [PubMed]
38. Xu Q, Yu H, Bennett J, He P, Zainon R, Doesburg R. Image reconstruction for hybrid true-color micro-CT. Biomedical Engineering, IEEE Transactions on. 2012;59:1711–1719. [PMC free article] [PubMed]
39. Xi Y, Chen Y, Tang R, Sun J, Zhao J. United iterative reconstruction for spectral computed tomography. Medical Imaging, IEEE Transactions on. 2015;34:769–778. [PubMed]
40. Leng S, Yu L, Wang J, Fletcher JG, Mistretta CA, McCollough CH. Noise reduction in spectral CT: reducing dose and breaking the trade-off between image noise and energy bin selection. Medical physics. 2011;38:4946–4957. [PubMed]
41. Yu Z, Leng S, Li Z, Ritman EL, McCollough CH. Fully three-dimensional image reconstruction in radiology and nuclear medicine. Newport, RI: 2015. Spectral PICCS reconstruction for photon-counting CT; pp. 198–201.
42. Wright SJ, Nocedal J. Numerical optimization. Vol. 2. New York: Springer; 1999.
43. Joseph PM. An improved algorithm for reprojecting rays through pixel images. Medical Imaging, IEEE Transactions on. 1982;1:192–196. [PubMed]
44. Zhuang W, Gopal S, Hebert T. Numerical evaluation of methods for computing tomographic projections. Nuclear Science, IEEE Transactions on. 1994;41:1660–1665.
45. Siddon RL. Fast calculation of the exact radiological path for a three - dimensional CT array. Medical physics. 1985;12:252–255. [PubMed]
46. Gordon R, Bender R, Herman GT. Algebraic reconstruction techniques (ART) for three-dimensional electron microscopy and X-ray photography. Journal of theoretical Biology. 1970;29:471–481. [PubMed]
47. Andersen A, Kak C. A SUPERIOR IMPLEMENTATION OF THE ART ALGORITHM [PubMed]
48. Elbakri IA, Fessler JA. Statistical image reconstruction for polyenergetic X-ray computed tomography. IEEE transactions on medical imaging. 2002;21:89–99. [PubMed]
49. Wang G, Jiang M. Ordered-subset simultaneous algebraic reconstruction techniques (OS-SART) Journal of X-ray Science and Technology. 2004;12:169–178.
50. Beck A, Teboulle M. A fast iterative shrinkage-thresholding algorithm for linear inverse problems. SIAM journal on imaging sciences. 2009;2:183–202.
51. Granton P, Pollmann S, Ford N, Drangova M, Holdsworth D. Implementation of dual-and triple-energy cone-beam micro-CT for postreconstruction material decomposition. Medical physics. 2008;35:5030–5042. [PubMed]
52. Wang Z, Bovik AC, Sheikh HR, Simoncelli EP. Image quality assessment: from error visibility to structural similarity. Image Processing, IEEE Transactions on. 2004;13:600–612. [PubMed]
53. Münch B, Trtik P, Marone F, Stampanoni M. Stripe and ring artifact removal with combined wavelet—Fourier filtering. Optics express. 2009;17:8567–8591. [PubMed]
54. Dabov K, Foi A, Katkovnik V, Egiazarian K. Image denoising by sparse 3-D transform-domain collaborative filtering. Image Processing, IEEE Transactions on. 2007;16:2080–2095. [PubMed]
55. Zhang Y, Xi Y, Yang Q, Cong W, Zhou J, Wang G. Spectral CT reconstruction using image sparsity and spectral correlation; Biomedical Imaging (ISBI), 2015 IEEE 12th International Symposium on; 2015. pp. 1600–1603.