Open Access
Add to BibSonomy
Abstract
Coded spectral X-ray computed tomography (CT) based on K-edge filtered illumination is a cost-effective approach to acquire both 3-dimensional structure of objects and their material composition. This approach allows sets of incomplete rays from sparse views or sparse rays with both spatial and spectral encoding to effectively reduce the inspection duration or radiation dose, which is of significance in biological imaging and medical diagnostics. However, reconstruction of spectral CT images from compressed measurements is a nonlinear and ill-posed problem. This paper proposes a material-decomposition-based approach to directly solve the reconstruction problem, without estimating the energy-binned sinograms. This approach assumes that the linear attenuation coefficient map of objects can be decomposed into a few basis materials that are separable in the spectral and space domains. The nonlinear problem is then converted to the reconstruction of the mass density maps of the basis materials. The dimensionality of the optimization variables is thus effectively reduced to overcome the ill-posedness. An alternating minimization scheme is used to solve the reconstruction with regularizations of weighted nuclear norm and total variation. Compared to the state-of-the-art reconstruction method for coded spectral CT, the proposed method can significantly improve the reconstruction quality. It is also capable of reconstructing the spectral CT images at two additional energy bins from the same set of measurements, thus providing more spectral information of the object.
© 2021 Optical Society of America under the terms of the OSA Open Access Publishing Agreement
1. Introduction
Spectral X-ray CT is an imaging modality that is used to reconstruct energy-dependent (i.e., spectral) attenuation mapping of an object, which can reveal both the inner structure of the examined object as well as its material composition [1,2]. It offers more advantages in medical diagnosis than conventional CT, which cannot resolve the material information [3,4]. In the past, several approaches have been developed to acquire spectrally resolved CT data. Fast kV switching is used to collect two measurement sets with distinct spectral attenuation characteristics by changing the voltage of the X-ray tube between successive projections [5,6]. However, the ability of spectral separation is limited because the X-ray tube voltage suffers under- and over-shooting during the switching process [7]. Dual-source CT systems are equipped with two X-ray sources with different tube voltages, which allows it to collect a pair of kV data simultaneously [8]. Nevertheless, such a system is expensive and the cross-scattered radiation between two sources can produce artifacts in CT images [9,10]. Dual-layer detectors consists of two conventional scintillation detectors, which can be used to obtain two sets of spectral CT data from a single polychromatic X-ray source [11]. The top layer and bottom layer of the detector absorb the low- and high-energy photons, respectively. However, the top layer also absorbs a portion of high-energy photons, thus leading to the inferior performance of spectral separation [12]. Photon counting detectors can acquire multiple raw spectrally resolved measurements with different truncating spectra by setting multiple spectral cut-off thresholds. Since the spectral resolution can be controlled by the cut-off thresholds, this method can achieve finer material decomposition results than others [13,14]. However, such detectors are not available currently for routine use due to the effect of pulse pile-up at high X-ray-flux rates [15,16]. Besides, the high cost of the photon counting detector is prohibitive in many applications [14].
K-edge filtering is another approach to attain spectral CT measurements. A K-edge filter is made of high atomic number material, of which the spectral attenuation has an abrupt increment at its K-edge energy [17]. Balancing the thicknesses of two filters can make their spectra almost overlap outside the range between the two K-edge energies [18]. Figure 1(a) illustrates a pair of balanced filters using Cerium (Ce) and Gadolinium (Gd) with respective thicknesses of 24 um and 17 um. Quasi-monochromatic spectrum in the range of [40.4, 50.2] KeV can be obtained by subtracting the low-Z filter spectrum (Ce) from the high-Z filter spectrum (Gd) [see Fig. 1(b)]. This method has been widely used in the fields of X-ray spectroscopy and crystallography [19]. In the past, Saito initially integrated the balanced filter pair of Erbium (Er) and Ytterbium (Yb) with conventional CT apparatus to acquire monochromatic CT data [21]. This method was further developed by Rakvongthai et al. [22], where sequential scans are performed with a set of balanced K-edge filters to obtain multiple monochromatic measurements. Standard filtered back projection (FBP) then was adopted to reconstruct the CT image of each energy bin from the monochromatic CT data [21,22].
Fig. 1. (a) Spectrum of 80kV X-ray source (black) generated with the Spektr software [20], filtered with Ce (cyan) and Gd (orange), vertical lines indicate the K-edge energy of Ce and Gd, 40.4 and 50.2 keV, respectively, where Z indicates the atomic number. (b) Quasi-monochromatic spectrum obtained by difference of the filter pair.
Although K-edge filters are cost-effective, the multiple sequential scanning process significantly increases the acquisition time and radiation dose. Inspired by spatial compressive encoded tomography [23], Cuadros et al. proposed a compressive spectral X-ray imaging (CSXI) model, which uses the encoded illumination in both spatial and spectral domains with the K-edge filters to capture the underlying information of spectral CT images in a single scan [24,25]. As shown in Fig. 2, four types of coding modes were proposed to generate various structure illuminations, including the uniform view (UV), random view (RV), uniform detector (UD), and random detector (RD) methods. In all of these modes, the X-ray source spectrum is spectrally coded by K-edge filters. By arranging the filters in different view angles and incident paths in the CT scan, the spatial information also is coded. The forward models in all of these cases follow an exponentially compressive measurement model in the spectral and spatial domains, which makes the reconstruction of CSXI data nonlinear and ill-posed. A multi-stage algorithm based on filtered intensity sinograms inpainting was subsequently developed to estimate the energy-binned sinograms from the integrating measurements [24]. The energy-binned CT images obtained from the CSXI system can be further used for material decomposition [26]. However, the inpainting process introduces errors that degrade the signal-to-noise ratio of the binned sinograms. It further impedes accurate reconstructions.
Fig. 2. Optical setup of a fan-beam spectral compressive X-ray CT with diverse encoding strategies. The different colors represent the X-ray beams filtered by different K-edge filters. (a) UV: uses filters cyclically assigned to all view angles; (b) RV: uses filters randomly assigned to all view angles; (c) UD: uses filters cyclically assigned to all projection paths; (d) RD: random coded aperture with K-edge filter is used for each of projection path.
Several algorithms have been proposed for spectral CT reconstruction [27–31]. Currently, they can be categorized into three approaches. The first approach is image-based which synthesizes the energy-binned images from the pre-reconstructed multiple polyenergetic images [2]. Since the polyenergetic image often has various artifacts, the accuracy of this method is not high [27,31]. The second is projection-based, including the aforementioned BFP reconstruction from sequential scan CT data and the multi-stage reconstruction for the CSXI system, which require first reconstructing the spectral CT sinograms [31]. However, the sinograms are always with reconstruction errors, which would impede the accuracy of subsequent reconstruction. The third method directly reconstructs the energy-combined image from the spectral CT data [29,32]. Since there are no errors introduced by the pre-reconstruction process, the direct inverse approach promises an accurate reconstruction for the fast-kV-switch-based and dual-energy-scan-based spectral CT systems [30]. Besides, material information of the object can be used as a prior to further improve the reconstruction [33,34]. Inspired by them, this paper develops a direct inverse approach for the CSXI system based on nonlinear material decomposition, which overcomes the limitations of current reconstruction algorithms for the CSXI system since it does not require the pre-estimation of the energy-binned sinograms. Given that the X-ray attenuation of an object can be approximated by the weighted sum of several basis materials’ attenuation maps [6,31,35–37], the proposed method directly estimates the weight coefficients (also known as material mass density) of the selected basis to reduce the dimensionality of the optimization variables. Joint regularization including the weighted nuclear norm (WNN) and total variation (TV) is adopted to improve the reconstruction quality by utilizing the spatial structure information of the material density maps. The regularized reconstruction problem is then solved with an alternating minimization scheme, in which the data fidelity term is minimized using the Limited-memory Broyden-Fletcher-Goldfarb-Shanno (L-BFGS) algorithm [38,39]. Finally, the energy-binned images are obtained by the weighted sum of the reconstructed density maps corresponding to the basis materials. Considering the inaccuracy introduced by the manufacturing error of the K-edge filter thickness, an empty scan is executed to predetermine the filter’s actual thickness. Simulations show that the reconstructed density images coincide with the distributions of the basis materials. The accuracy of the obtained energy-binned images is effectively improved compared with that of state-of-the-art methods including the multi-stage algorithm in [24], and the FBP reconstruction from the sequential scan measurements [22].
The organization of the rest of this paper is as follows. Section 2 describes the forward projection model of the spectral CT based on the K-edge coded illuminations. Section 3 presents the details of the proposed reconstruction algorithm. Section 4 presents the simulations to assess the proposed algorithm with diverse system settings and phantoms. Section 5 includes conclusions and discussions.
2. Forward projection model
According to the spectral version of the Beer-Lambert law, when a polychromatic X-ray beam with initial intensity $I_0(E)$ passes through an object, the attenuated intensity can be modeled as [40]:
where $D\left (E\right )$ is the spectral response function of the detector, which is proportional to the X-ray photon’s energy $E$ for the conventional energy integrating detectors. The $\mu \left (E,x\right )$ denotes the linear attenuation coefficient of the object. Given the CSXI system depicted in Fig. 2, a K-edge filter with linear attenuation coefficient $\mu _{f}(E)$ and desired thickness $t_0$ is placed in the projection path to encode the spectrum of X-ray beam. The spectral attenuation of the filter also follows the Beer-Lambert law and is formulated as $f_0(E)=\exp {[-\mu _f(E)t_0/\cos {\psi }]}$, where $\psi$ denotes the angle between the normal direction of the filter and the projection path [24]. The thickness of the filters invariably come within an unknown manufacturing tolerance. Here, an empty scan is adopted to predetermine the manufacturing thickness error $\Delta t$. The actual filtered intensity of X-ray measured without object is given by3. Proposed reconstruction method
3.1 Inverse problem
According to the forward mapping model of CSXI system described in Eq. (5), the estimation of the energy-binned attenuation images is formulated as the following optimization problem
3.2 Basis material decomposition
Based on the mixture rule [41], the linear attenuation coefficient of an object can be decomposed as the weighted sum of $M$ different basis functions that are separable in energy and space:
where $\rho _m$ represents the weight factors of $M$ linearly independent basis material attenuation functions $\tau _1, \ldots , \tau _M$. An alternative set of basis functions is material-based, where $\tau _m$ is chosen empirically according to the main components of the object under scanning [35]. In such material-based decomposition, $\tau _m$ and $\rho _m$ are referred as the mass attenuation coefficients and mass density of the material, respectively.Based on the BMD method described in Eq. (7), given a scanned object composed of $M$ known basis materials, the energy-binned linear attenuation coefficients $\mathbf {X}$ can be approximately represented by $\mathbf {X}=\mathbf {P}\mathbf {T}$, where each column $\mathbf {P}_{(:, m)}$ of $\mathbf {P}\in \mathbb {R}^{NN\times {M}}$ corresponds to the vectorized mass density image of the $m\textrm {th}$ material, and $\mathbf {T}\in \mathbb {R}^{M\times K}$ represents the energy-binned mass attenuation coefficients of $M$ basis materials at a total of $K$ energy bins. Substituting $\mathbf {X}=\mathbf {P}\mathbf {T}$ into Eq. (6), the equivalent optimization problem is formulated as
3.3 Regularized nonlinear reconstruction
For the purpose of achieving higher-quality reconstructions, we introduce the low rank and total variation (TV) regularization into Eq. (8) to simultaneously utilize the global and local prior information in the density image of the basis materials. The use of low-rank regularization is based on the assumption that a natural image usually consists of self-similar patches. Meanwhile, the image gradient often exhibits long-tailed distribution, with which the total variation regularization is thus used as another penalty term. The weighted nuclear norm (WNN) minimization is used to formulate the low rank regularization. Given a weight vector $\mathbf {w}$ with all elements being non-negative, the WNN of a 2D image $\mathbf {X}\in \mathbb {R}^{N\times {N}}$ is defined as [42]:
where $\sigma _i(\mathbf {X})$ is the $i\textrm {th}$ singular value of $\mathbf {X}$. Due to the weight vector $\mathbf {w}$ assigning the singular values with different weights, WNN enhances the capacity and flexibility of the original nuclear norm. In the image denoising application, the larger singular value is, the less it should be shrunk. A practical setting of the weights is $\mathbf {w}_i=c\sqrt {N}/(\sigma _i(\mathbf {X})+\varepsilon )$, where $c$ is a constant, and $\varepsilon$ is a small positive number to avoid dividing by 0 [42]. We use the TV regularization for per material density image $\mathbf {P}_{(:, m)}$ to simultaneously preserve edges and remove noise in flat regions [43]. For a 2D image $\mathbf {X}\in \mathbb {R}^{N\times N}$, it is defined as [44]:Update $\mathbf {P}$: The variable $\mathbf {P}$ is updated by minimizing the following large-scale nonlinear optimization problem:
The quasi-Newton direction $\mathbf {D}$ = -$\mathbf {H}_t\mathbf {g}_t$ is used as the descent direction. An exact line search on the first-order Taylor expansion of Eq. (5) is adopted to find the optimal step size $\alpha$ for each iteration. Based on the approximation $\operatorname {e}^{x}\approx {1+x}$ for $x\approx {0}$, we have
Update $\mathbf {Q}$: The subproblem of updating $\mathbf {Q}$ in Eq. (12) can be formulated as
Update $\mathbf {R}$: When $\mathbf {P}$ and $\mathbf {Z}_{\textrm {TV}}$ are fixed, we can update $\mathbf {R}$ by solving the following TV minimization problem
Update $\mathbf {Z}_*$ and $\mathbf {Z}_{\textrm {TV}}$: After the update of $\mathbf {P}$, $\mathbf {Q}$, and $\mathbf {R}$, the auxiliary variables $\mathbf {Z}_*$ and $\mathbf {Z}_{\textrm {TV}}$ can be updated as following
We alternately update the variables in Eq. (12) with the above steps until the algorithm converges, then we obtain the spectral CT image $\mathbf {X}$ from the reconstructed density image $\mathbf {P}$ by $\mathbf {X=PT}$.4. Numerical simulations
A conventional fan beam CT with a line detector array including 512 elements is used in the simulations. The distance from the source to object is 400 mm, and the distance from detector to object is 413 mm. The number of the view angles is 768 over $360^{\circ }$. The dimensionality of the object is $256\times 256$. Both the pitch size of the detector and the pixel size of the object are 2 mm. The system matrix $\mathbf {A}$ is generated using the ASTRA tomography toolbox according to the parameters mentioned above [48]. In the following simulations, different X-ray spectrum settings and phantoms are employed. The first simulation uses an X-ray source with 80 kV tube voltage, and a phantom without contrast agent. The second simulation uses an X-ray source with 120 kV tube voltage, and a phantom containing iodine contrast agent. To evaluate the performance of the material-decomposition-based method, the results of the state-of-the-art methods are also presented, including the multi-stage algorithm based on sinogram inpainting [24], and the FBP reconstruction based on sequentially scan [22]. In the simulations, we mainly demonstrate the proposed method in terms of material decomposition accuracy, reconstruction image quality, statistical analysis, and the robustness to filter thickness error.
The quality of the reconstructed image is assessed by the peak signal to noise ratio (PSNR) and the relative root mean square error (rRMSE). Suppose the reconstructed image is ${\mathbf {X}}$ with a dimension of $N\times N$, and the corresponding ground truth is $\hat {\mathbf {X}}$. The PSNR of the reconstructed image is calculated as $\textrm {PSNR}=10\log _{10}\left ({\mathbf {\hat {X}}_{\textrm {max}}^{2}/\left (\Vert \mathbf {X-\hat {X}}\Vert ^{2}_{\textrm {F}}/N^{2}\right )}\right )$, where $\mathbf {\hat {X}}_{\textrm {max}}$ denotes the maximum pixel value in $\hat {\mathbf {X}}$ [40]. In addition, the definition of rRMSE is $\operatorname {rRMSE}=\sqrt {\Vert \mathbf {X-\hat {X}}\Vert ^{2}_{\textrm {F}}/\Vert \mathbf {\hat {X}}\Vert ^{2}_{\textrm {F}}}$ [49].
4.1 Low kV setting and phantom without contrast agent
In the simulation with low kV setting, the tube voltage of X-ray source is set to be 80 kV, and the corresponding spectrum $I_0(E)$ is generated with the Spektr software [20]. A copper filter with thickness of 10 um is placed in front of the X-ray source, which is used to filter out the photons with energy lower than 20 keV, such that the source spectrum is constrained within the range from 20 keV to 80 keV. Five different filter materials [Antimony (Sb), Cerium (Ce), Gadolinium (Gd), Ytterbium (Yb), Tungsten (W)] are chosen to fabricate the K-edge filter set, such that the spectrum is divided into six energy bins, each of with has an approximate bandwidth of 10 keV. The linear attenuation coefficients of the K-edge filter are obtained from the X-ray attenuation databases provided by the National Institute of the Standards and Technology (NIST) [50]. The thicknesses of the filters are set as [61.1, 46.9, 32.5, 30.7, 10.0] um respectively, such that any two of the K-edge filters can form a Ross filter pair. Figure 3(a) shows the original spectrum and the filtered spectra. The numbers of the initial photons in the filtered spectra are $[1.17, 1.42, 1.63, 1.75, 1.75,1.78]\times 10^{6}$ respectively. Figure 3(b) shows that the quasi-monoenergetic spectra can be obtained by subtracting the consecutive spectra.
Fig. 3. (a) Original 80 kV X-ray source spectrum, and filtered with Sb, Ce, Gd, Yb, and W filters. The vertical lines indicate the K-edge energy respectively. (b) Quasi-monochromatic energy-binned spectra obtained by difference of the filter pairs, additional two spectra at intervals are indicated by chain lines.
As shown in Fig. 4(a), the scanned object has a dimension of $N=256$, which is generated from a thoracic abdominal CT scanning from the 3D-IRCADb database [51]. It is an image in DICOM format in which the gray value corresponds to the HU number. The database also provides the manual segments of the CT scan according to pixel material, which is shown in Fig. 4(b). Using the CT scan of the thorax and its segmented image, the multi-energy phantom is generated for the following simulations. Firstly, we obtain the effective attenuation coefficient $\mu$ from the HU number using the relationship given by $\textrm {HU number} = \frac {\mu - \mu _{\textrm {water}}}{\mu _{\textrm {water}}} \times 1000$, where $\mu _{\textrm {water}}$ represents the effective attenuation coefficient of water. Then, we derive the mass density by proportionally scaling $\mu$ according to the reference values of lungs’ density in NIST data [50]. Finally, we assign the effective linear attenuation coefficients for each energy bin by multiplying the material’s mass density and mass attenuation coefficient. The energy-binned data cube is then used for the K-edge encoded forward projection to generate the filtered intensity data. As mentioned in Section 1, the encoding strategies include the UV, RV, UD, and RD methods. The filtered intensity $\mathbf {y}$ is used as the parameter in the Poisson distribution to generate the noise-polluted measurement data [52]. For the CT scanning of human body, the scanned object should consist of water, soft tissue, and/or bone, etc. Figure 5 shows the mass attenuation coefficients of these basis materials. It can be found that the curves of different materials are similar except the bone and Iodine. Thus, we can combine them and get some basis material sets, such as water/bone set (WB), soft tissue/bone (SB) set, and water/iodine (WI) set. Here, the SB set is used as the basis materials to reconstruct the phantom without contrast agent contained. The parameters of regularizattion terms are set emperically, where $\mu _* = \mu _{\textrm {TV}} = 10^{-2}$, $\lambda _*=5\times 10^{-2}$, and $\lambda _{\textrm {TV}}=10^{-4}$. The choice of $\mu _*$ and $\mu _{\textrm {TV}}$ is to effectively minimize the data fidelity term. The choice of $\lambda _*$ and $\lambda _{\textrm {TV}}$ is to denoise while preserving the details in the reconstructed images. They are equal for all of our simulations.
Fig. 4. (a) 3D-IRCADb thoracic abdominal CT scan slice. (b) Segmented CT scan according to pixel material (1) heart (blood), (2) lung, (3) adipose tissue, (4) bone, (5) aorta (blood), (6) stomach (soft tissue), (7) air. (c) The estimated image of mass density in $\textrm {g}\cdot \textrm {cm}^{-3}.$
4.1.1 Performance of reconstructing the material density image
Considering that the proposed method is based on material decomposition, we first exhibit the reconstructed mass density images as shown in Fig. 6. The first column shows the material density images $\hat {\mathbf {P}}$, which are used as the ground truth. They are obtained by decomposing the original energy-binned CT images $\hat {\mathbf {X}}$ on the selected basis materials $\mathbf {T}$, i.e., solving the linear least square problem $\hat {\mathbf {P}}=\underset {\mathbf {P}}{\operatorname {argmin}}\left \lVert {\hat {\mathbf {X}}-\mathbf {PT}}\right \rVert ^{2}_{\mathrm {F}}$. The figures from the second column to the fifth column in Fig. 6 correspond to the reconstructed material mass density images for the CSXI system using the UV, RV, UD, and RD K-edge coded apertures. It is found that the results of our proposed method are consistent with the ground truth. Given that the multi-stage algorithm has similar reconstruction performance for all the encoding modes in CSXI [24], hereafter we only present its results of material decomposition using the UV encoding for comparison. For the purpose of comparison, we also present the material decomposition results from the sequential scan CT (SS-CT) with multiple balanced K-edge filter sets [22]. To be fair, the number of view angles in SS-CT is set to be $153$, i.e., $1/F$ of view angles for the CSXI system. The quasi-monoenergetic sinograms are obtained using the difference of the Ross filter pairs, where the negative values are removed. It is found that the artifacts are introduced in the reconstructed images attributed to the error of the multi-stage reconstruction algorithm (see the sixth column in Fig. 6). Artifacts in the background are also found in the SS-CT material images of soft tissue, since it uses less projection view angles.
Fig. 6. Comparison of the material decomposition reconstruction. Top to bottom: the material density images of soft tissue and bone. The first column: the ground truth of the material images. The second to fifth column: the reconstruction using proposed method for CSXI with diverse encoding modes. The sixth column: the reconstruction using method in [24]. The seventh column: the reconstruction of SS-CT [22].
Table 1 presents the PSNR and rRMSE values corresponding to the reconstructed material density images in Fig. 6. It shows that the proposed method obtains higher reconstruction quality compared to the results of CSXI [24] and SS-CT [22].
Table 1. PSNR and rRMSE Values of the Reconstructed Material Density Images
4.1.2 Performance of reconstructing the energy-binned images
Since the CSXI reconstruction results with all encoding modes are similar, in the subsequent simulations we only present the results using the UV encoding to avoid redundancy. We compare the reconstructed energy-binned images obtained by our BMD based algorithm and the existing CSXI method in [24] and the SS-CT method in [22]. The reconstructed energy-binned images are presented in Fig. 7. Artifacts in the third row indicate that the error of inpainting would degrade the image quality of CSXI. It is also observed that the reconstructed images of the proposed method are very similar to the ground truth. That means the direct reconstruction approach based on the material decomposition is feasible and profitable for the CSXI system. Table 2 presents the PSNR and rRMSE values corresponding to the energy-binned images shown in 7. Note that the linear attenuation coefficient images of the first and last energy bins can be obtained by the proposed method, although the qualities of the images are relatively low. The average reconstruction PSNR within the energy bin $k=[2,5]$ is 49.77 dB, which is 9.44 dB higher and 21.90 dB higher than the two aforementioned methods, respectively.
Fig. 7. Comparison of the energy-binned reconstruction results. From top to bottom are reconstructions obtained from: ground truth, proposed method, CSXI [24], SS-CT [22]. The display window is set to [-1000 HU, 1000 HU]. Gamma correction with the parameter $\gamma = 0.6$ is used for both the ground truth and the reconstructed images to improve the display of the lungs’ details.
Table 2. PSNR and rRMSE Values of Reconstructed Energy-Binned Attenuation Coefficients With Diverse Algorithms
Figure 8 focuses on the zoomed parts A, B, and C indicated by the red squares in Fig. 7. They are extracted from the heart(blood), lung, and bone areas, respectively. Evident noise and artifacts can be found in the zoomed parts of the image reconstructed by the methods in [24] and [22]. Compared with them, the proposed method provides more accurate image details, especially in the soft tissue like areas, such as the heart and lungs.
Fig. 8. Zoomed images corresponds to parts A, B, and C in Fig. 7. The corresponding display windows are set to [-300 HU, 200 HU], [-1000 HU, -500 HU], and [2500 HU, 3500 HU] to show more details of the reconstruction.
4.1.3 Statistical analysis
Five region-of-interests (ROIs) are respectively drawn in the adipose tissue and heart areas used for statistical analysis. As shown in Fig. 7, the yellow squares are in adipose area, and the blue squares are in the heart area. The same ROIs are put in the reconstructed images. Table 3 presents the mean and SD values of the the HU number in those ROIs. It shows that the proposed method can achieve higher reconstruction accuracy compared to the other methods in this simulation.
Table 3. Mean$\pm$SD Values of the Reconstructed Energy-Binned HU Number With Diverse Algorithms
The contrast-to-noise ratio (CNR) are calculated for each reconstruction. The CNR of heart ROIs with respect to adipose ROIs is defined as [40]
Table 4. CNR Values of the Reconstructions in Fig. 7
4.1.4 Impact of noise levels on reconstruction
It is known that the Poisson noise level introduced into the measured data is related to the number of X-ray photons. Thus, we respectively set the blank intensity as 1/1, 1/3, 1/15, and 1/75 of the initial blank intensity to validate the performance of the proposed algorithm on various noise levels. The corresponding numbers of the blank X-ray photon are $1.5\times 10^{6}$, $5\times 10^{5}, 1\times 10^{5}, \textrm {and } 2\times 10^{4}$, respectively. Table 5 provides the average PSNR and rRMSE values of the energy-binned reconstruction images using these intensity settings.
Table 6 gives the corresponding mean and SD values of the ROIs shown in Fig. 7. Compared with the images reconstructed from the high dose measurements, the reconstructed image of the proposed method with 1/15 dose still maintains good accuracy. The deterioration of the image quality is mainly caused by the variance of HU number introduced by the measurement noise. When the blank photons are further reduced to 1/75 of the original dose, the quality of reconstructed image is unacceptable since the variance of HU number is too large relative to the ground truth.
Table 5. Average PSNR and rRMSE Values of Reconstruction Using Various Blank Photons
Table 6. Mean $\pm$ SD Values of the Reconstructed Energy-Binned HU Number Using Various Blank Photons
4.1.5 Impact of the filter’s thickness error on reconstruction
This part evaluates the performance of the proposed method with respect to the inaccuracy of filter thickness. It is noted that manufacturing errors are randomly introduced for all the filters in the K-edge coded apertures. Table 7 presents the mean of the PSNR and rRMSE values of the CSXI reconstruction, where the filters are with additional $\pm 10\%$ and $\pm 20\%$ thickness errors, respectively. The results without correcting the filter thickness are also provided for comparison.
Table 7. PSNR and rRMSE Values of Reconstruction for the CSXI System With Inaccurate Filter Thickness
It is found that the performance of proposed algorithm is related to the manufacturing accuracy of the filters. Compared with the reconstruction results shown in Table 2, in the presence of manufacturing error of $\pm 10\%$ for every filter’s thickness, the average reconstruction PSNR value is degraded by 7.11 dB. For the manufacturing error of $\pm 20\%$, the degradation becomes 10.25 dB. Nevertheless, the re-correction of the filter thickness described in Section II helps to overcome this disadvantage. As presented in Table 7, after calibrating the filter’s thickness, the average PSNR values only drop by 0.65 dB and 2.93 dB, respectively. These degraded results are still better than the results of competitive methods shown in Table 2. It means that the proposed method is more robust to the inaccuracy of filter thickness.
4.1.6 Impact of the basis material set on reconstruction
The choice of basis material set is based on apriori knowledge of the components of the scanned object. Here, we study the influence of different basis material sets on the reconstruction performance. Table 8 presents the reconstruction PSNR and rRMSE values using diverse basis materials. It shows that the image quality is improved when the adipose is added into the basis material set. That is because adipose has different attenuation characteristic and thus improve the representation capacity of the basis material set (see Fig. 5). However, when the blood is added into the basis material set, the improvement of reconstruction quality is negligible. There are two main reasons for this phenomenon. On one hand, the blood has similar mass attenuation coefficients to the soft tissue, because their main constituent elements are the same. Thus, adding it to the basis material set is redundant for reconstruction. On the other hand, adding new variables will increase the dimensionality of the optimization problem, which slows down the convergence.
Table 8. PSNR and rRMSE Values of Reconstruction With Different Basis Materials Set
4.2 High kV setting and phantom with contrast agent
This section illuminates the performance of the proposed method with high voltage setting. The tube voltage of the X-ray source in this simulation is set to be 120 kV. The corresponding spectrum $I_0(E)$ is generated with the Spektr software [20]. A copper filter with thickness of 150 um is placed in front of the X-ray source to remove the photons with energy less than 20 keV. Five materials [Barium (Ba), Samarium (Sm), Thulium (Tm), Osmium(Os), Lead (Pb)] are chosen to fabricate the K-edge filters. According to their K-edge energies, the narrow energy bins are given by [20.0, 37.4], [37.4, 46.8], [46.8, 59.4], [59.4, 73.9], [73.9, 88.0], and [88.0, 120.0] keV. The linear attenuation coefficients of the K-edge filters come from the NIST X-ray attenuation databases available in [50]. The thicknesses of the filters are set as [92.6, 15.9, 7.9, 1.1, 2.0] um respectively, such that two K-edge filters can form a Ross filter pair. Figure 9(a) shows the original and filtered spectra. Subsequently, Fig. 9(b) shows the quasi-monoenergetic spectra obtained by the subtracting the consecutive spectra.
Fig. 9. (a) Original 120 kV X-ray source spectrum, and filtered with Ba, Sm, Tm, Os, and PB filters. (b) Quasi-monochromatic energy-binned spectra obtained by difference of the filter pairs, additional two spectra are indicated by chain lines.
A modified Forbild thorax phantom generated using the CONARD software is used in this simulation [53]. As shown in Fig. 10, the phantom consists of eight materials. The densities of these materials are set according to [40], and their mass attenuation coefficients can be obtained from the NIST X-ray attenuation databases available in [50]. Using this database, we assign linear attenuation coefficients in six energy bins for all the phantom pixels. The phantom is then used for forward projection with the UV-CSXI system. Poisson noises are then added to the filtered projection data.
Fig. 10. Forbild thorax phantom consists of eight materials, (1) air, (2) lung, (3) average soft tissue, (4) heart (blood), (5) $99.1\%$ blood + $0.9\%$ iodine, (6) artery (blood), (7) bone, (8) marrow.
Considering the Iodine contrast agent contained in the phantom, the basis material set of soft tissue, bone, and Iodine (SBI) is used in the reconstruction. The obtained mass density images are presented in the second column in Fig. 11. For comparison, we also present the material density images reconstructed by the methods of original CSXI [24] and SS-CT [22]. The soft tissue can be still found in the reconstructed density image of bone obtained by the method in [24], but not in the image reconstructed by the proposed method. Furthermore, it is noted that the CSXI reconstruction with method in [24] present numerous artifacts due to the estimation error introduced at the stage of the sinogram inpainting. Table 9 presents the PSNR and rRMSE values corresponding to the reconstructed material density images shown in Fig. 11. It shows that the density images reconstructed by the proposed method have higher reconstruction quality. Fig. 12 presents the reconstructed energy-binned images corresponding to the density images. The corresponding PSNR and rRMSE values are shown in Table 10. It is noted that there are some gaps between these results and the previous reconstruction results in Fig. 7 and Table 2. The PSNR values of each energy bin are lower than before. Especially, errors in the first energy bin are found near the segmented areas of bone and iodine contrast agent. Such errors are introduced by the existence of Iodine in the phantom and basis material set. As shown in Fig. 5, compared with the other materials, iodine has completely different X-ray attenuation characteristics. The K-edge absorption at the energy of 33.2 keV makes the linear attenuation coefficient no longer continuously decreasing with the increment of energy.
Fig. 11. Material decomposition reconstruction of the modified thorax phantom. Top to bottom: the material density image of soft tissue, bone, and Iodine. Left to right: the mass density images from ground truth, proposed method, CSXI [24], and SS-CT [22].
Fig. 12. Comparison of the energy-binned reconstruction results. Images from top to bottom row correspond to reconstructions obtained from ground truth, proposed method, CSXI [24], and SS-CT [22]. The display window is set to [-1000 HU, 1000 HU].
Table 9. PSNR and rRMSE Values of Material Density Images Reconstructed in Diverse Encoding Modes
Table 10. PSNR and rRMSE Values of the Reconstructed Energy-Binned Forbild Thorax Phantom
Figure 13 shows the zoomed areas indicated by the square A and square B in Fig. 12. They are drawn from the junction areas between different tissues. Note that the blood and soft tissue are distinguishable in the reconstructed images of the proposed method, but are submerged in the noise in the reconstructed images of other competitive methods.
Fig. 13. Zoomed image corresponds to parts A and B in Fig. 12. The display windows are set to [-500HU, 400HU] and [-200HU, 300HU] respectively.
4.2.1 Statistical analysis
For each tissue in the Forbild thorax phantom, we present the ROIs indicated by the blue square shown in Fig. 12. Table 11 presents the mean and SD values of the reconstructed HU number in the ROIs at the third and the fourth energy bins. It shows that the reconstructed image using the proposed method has less bias and lower variance compared to the competitive methods.
Table 11. Mean $\pm$ SD Values of the Reconstructed HU Number at the Third and Fourth Energy Bins
Given the similarity of the tissue’s X-ray attenuation, we calculate the CNR of the blood with respect to the soft tissue, as well as the CNR of the bone with respect to the iodine-mixed blood.
They are respectively denoted as $\operatorname {CNR}_{\textrm {blood}/\textrm {soft}}$ and $\operatorname {CNR}_{\textrm {bone}/\textrm {iodine}}$. Table 12 presents the CNR values of the reconstruction results. The CNR values of the proposed method outperform the sinogram inpainting based method proposed in the original CSXI system [24], and also outperform the FBP reconstruction method for the SS-CT system [22].
Table 12. Mean CNR Values of the ROIs in the Energy Binned Forbild Thorax Phantom
5. Conclusions and discussions
Image reconstruction from CSXI data is a high-dimensional ill-posed and nonlinear problem. This article develops a inverse optimization approach to directly reconstruct the energy-binned images from the CSXI measurements. This approach does not use sinogram inpainting, thus avoids the introduced error and attains significant improvement in reconstruction quality. Inspired by other forms of spectral CT, the BMD approach is used to reduce the dimensionality of variables, and to obtain two more reconstructed images in the first and last energy bins. The proposed method is evaluated in numerical studies on various simulation settings. The results show that the proposed approach yields better image quality. Statistical analysis in terms of bias, variance, and CNR values also outperforms the results of the competitive methods. Meanwhile, as a direct inverse approach, the proposed algorithm is inferior to the other methods on time consumption because the nonlinear iterative operations require additional computation overhead. Future research will explore the use of pre-conditioning to save reconstruction time.
While our proposed method provides superior results for all four coding modes, it is also necessary to compare the simplicity of implementation of the various modes. In the UV and RV modes, the type of the K-edge filter used in front of the X-ray source changes with the view angle. Furthermore, the type also changes with the projection path in the UD and RD modes. These require the coded aperture to be individually fixed to their corresponding view angles or projection paths. The K-edge coded apertures require accurate manufacturing and system calibration. Our group is currently developing them using advanced additive 3-D printing manufacturing. A description of the techniques under development will be reported in a future publication.
Future research would focus on improving the CSXI method from two aspects. First, we would explore better regularizations and constraints to improve the reconstruction quality at the first and last energy bins. Second, inspired by coded aperture optimization in compressive X-ray CT [52,54,55], we would study the optimization approach of the K-edge coded aperture to acquire the optimal measurement data from limited projections.
Funding
National Key Research and Development Program of China (2019YFB2102300, 2019YFB2102301); National Natural Science Foundation of China (61936014); Science and Technology Innovation Plan Of Shanghai Science and Technology Commission (19511103302); Fundamental Research Funds for the Central Universities; National Science Foundation (CIF 1717578); University Dissertation Award, and by a UNIDEL/W. L. Gore grant..
Disclosures
The authors declare no conflicts of interest.
Data availability
Data underlying the results presented in this paper are available in the following Refs. [20,50,51,53].
References
1. W. A. Kalender, E. Klotz, and C. Suess, “Vertebral bone mineral analysis: an integrated approach with CT,” Radiology 164(2), 419–423 (1987). [CrossRef]
2. B. J. Heismann, B. T. Schmidt, and T. Flohr, Spectral computed tomography, (SPIE, 2012), pp. 24–54.
3. A. C. Silva, B. G. Morse, A. K. Hara, R. G. Paden, N. Hongo, and W. Pavlicek, “Dual-energy (spectral) CT: applications in abdominal imaging,” RadioGraphics 31(4), 1031–1046 (2011). [CrossRef]
4. D. Cabib, R. A. Buckwald, Z. Malik, Y. Garini, N. Katzir, and D. G. Soeknsen, “Spectral bio-imaging methods for biological research, medical diagnostics and therapy,” (1998). US Patent 5, 784, 162.
5. D. Xu, D. A. Langan, X. Wu, J. D. Pack, T. M. Benson, J. E. Tkaczky, and A. M. Schmitz, “Dual energy ct via fast kVp switching spectrum estimation,” in Medical Imaging 2009: Physics of Medical Imaging, vol. 7258 (International Society for Optics and Photonics, 2009), p. 72583T.
6. M. M. Goodsitt, E. G. Christodoulou, and S. C. Larson, “Accuracies of the synthesized monochromatic CT numbers and effective atomic numbers obtained with a rapid kVp switching dual energy CT scanner,” Med. Phys. 38(4), 2222–2232 (2011). [CrossRef]
7. B. Krauss, K. L. Grant, B. Schmidt, and T. Flohr, “The importance of spectral separation: an assessment of dual-energy spectral separation for quantitative ability and dose efficiency,” Invest. Radiol. 50(2), 114–118 (2015). [CrossRef]
8. L. Yu, J. A. Christner, S. Leng, J. Wang, J. G. Fletcher, and C. H. McCollough, “Virtual monochromatic imaging in dual-source dual-energy CT: radiation dose and image quality,” Med. Phys. 38(12), 6371–6379 (2011). [CrossRef]
9. K. J. Engel, C. Herrmann, and G. Zeitler, “X-ray scattering in single-and dual-source CT,” Med. Phys. 35(1), 318–332 (2007). [CrossRef]
10. M. Petersilka, K. Stierstorfer, H. Bruder, and T. Flohr, “Strategies for scatter correction in dual source CT,” Med. Phys. 37(11), 5971–5992 (2010). [CrossRef]
11. R. Carmi, G. Naveh, and A. Altman, “Material separation with dual-layer ct,” in IEEE Nuclear Science Symposium Conference Record, 2005, vol. 4 (2005), pp. 3–1878.
12. L. Saba and J. S. Suri, Multi-detector CT imaging: abdomen, pelvis, and CAD applications, vol. 2 (CRC Press, 2013).
13. P. M. Shikhaliev, T. Xu, and S. Molloi, “Photon counting computed tomography: concept and initial results,” Med. Phys. 32(2), 427–436 (2005). [CrossRef]
14. K. Taguchi and J. S. Iwanczyk, “Vision 20/20: Single photon counting X-ray detectors in medical imaging,” Med. Phys. 40(10), 100901 (2013). [CrossRef]
15. J. Tanguay, H. K. Kim, and I. A. Cunningham, “The role of X-ray swank factor in energy-resolving photon-counting imaging,” Med. Phys. 37(12), 6205–6211 (2010). [CrossRef]
16. W. C. Barber, E. Nygard, J. S. Iwanczyk, M. Zhang, E. C. Frey, B. M. W. Tsui, J. C. Wessel, N. Malakhov, G. Wawrzyniak, N. E. Hartsough, T. Gandhi, and K. Taguchi, “Characterization of a novel photon counting detector for clinical CT: count rate, energy resolution, and noise performance,” in Medical Imaging 2009: Physics of Medical Imaging, vol. 7258 International Society for Optics and Photonics (SPIE, 2009), pp. 695–703.
17. J. A. Bearden and A. Burr, “Reevaluation of X-ray atomic energy levels,” Rev. Mod. Phys. 39(1), 125–142 (1967). [CrossRef]
18. P. Kirkpatrick, “On the theory and use of ross filters,” Rev. Sci. Instrum. 10(6), 186–191 (1939). [CrossRef]
19. W. Miao, Y. Ding, K. Sun, D. Lai, Z. Chen, H. Wang, J. Cheng, Z. Zheng, and H. Peng, “Design of the ross pairs for soft X-ray spectrometer,” in IEEE Conference Record-Abstracts. 2002 IEEE International Conference on Plasma Science (Cat. No. 02CH37340), (IEEE, 2002), p. 161.
20. J. Siewerdsen, A. Waese, D. Moseley, S. Richard, and D. Jaffray, “Spektr: A computational tool for X-ray spectral analysis and imaging system optimization,” Med. Phys. 31(11), 3057–3067 (2004). [CrossRef]
21. M. Saito, “Quasimonochromatic X-ray computed tomography by the balanced filter method using a conventional X-ray source: Quasimonochromatic X-ray CT by balanced filter method,” Med. Phys. 31(12), 3436–3443 (2004). [CrossRef]
22. Y. Rakvongthai, W. Worstell, G. El Fakhri, J. Bian, A. Lorsakul, and J. Ouyang, “Spectral CT using multiple balanced k-edge filters,” IEEE Trans. Med. Imaging 34(3), 740–747 (2015). [CrossRef]
23. D. J. Brady, A. Mrozack, K. MacCabe, and P. Llull, “Compressive tomography,” Adv. Opt. Photonics 7(4), 756–813 (2015). [CrossRef]
24. A. Cuadros, X. Ma, and G. R. Arce, “Compressive spectral X-ray tomography based on spatial and spectral coded illumination,” Opt. Express 27(8), 10745 (2019). [CrossRef]
25. A. Cuadros and G. Arce, “Pixelated k-edge coded aperture system for compressive spectral X-ray imaging,” (2020). US Patent 10, 874, 361.
26. A. P. Cuadros, X. Ma, and G. R. Arce, “Compressive X-ray material decomposition using structured illumination,” in Developments in X-Ray Tomography XII, vol. 11113 (International Society for Optics and Photonics, 2019), p. 111131H.
27. X. Zhao, P. Chen, J. Wei, and Z. Qu, “Spectral ct imaging method based on blind separation of polychromatic projections with poisson prior,” Opt. Express 28(9), 12780–12794 (2020). [CrossRef]
28. J. Liu, H. Ding, S. Molloi, X. Zhang, and H. Gao, “Ticmr: Total image constrained material reconstruction via nonlocal total variation regularization for spectral ct,” IEEE Trans. Med. Imaging 35(12), 2578–2586 (2016). [CrossRef]
29. X. Dong, T. Niu, and L. Zhu, “Combined iterative reconstruction and image-domain decomposition for dual energy ct using total-variation regularization,” Med. Phys. 41(5), 051909 (2014). [CrossRef]
30. C. Mory, B. Sixou, S. Si-Mohamed, L. Boussel, and S. Rit, “Comparison of five one-step reconstruction algorithms for spectral ct,” Phys. Med. Biol. 63(23), 235001 (2018). [CrossRef]
31. N. Ducros, J. F. P. Abascal, B. Sixou, S. Rit, and F. Peyrin, “Regularization of nonlinear decomposition of spectral X-ray projection images,” Med. Phys. 44(9), e174–e187 (2017). [CrossRef]
32. R. F. Barber, E. Y. Sidky, T. G. Schmidt, and X. Pan, “An algorithm for constrained one-step inversion of spectral ct data,” Phys. Med. Biol. 61(10), 3784–3818 (2016). [CrossRef]
33. H. Gao, H. Yu, S. Osher, and G. Wang, “Multi-energy CT based on a prior rank, intensity and sparsity model (PRISM),” Inverse Probl. 27(11), 115012 (2011). [CrossRef]
34. W. Zhao, D. Vernekohl, F. Han, B. Han, H. Peng, Y. Yang, L. Xing, and J. K. Min, “A unified material decomposition framework for quantitative dual-and triple-energy ct imaging,” Med. Phys. 45(7), 2964–2977 (2018). [CrossRef]
35. J. Weaver and A. L. Huddleston, “Attenuation coefficients of body tissues using principal-components analysis,” Med. Phys. 12(1), 40–45 (1985). [CrossRef]
36. B. Heismann, J. Leppert, and K. Stierstorfer, “Density and atomic number measurements with spectral X-ray attenuation method,” J. Appl. Phys. 94(3), 2073–2079 (2003). [CrossRef]
37. B. Heismann, “Signal-to-noise monte-carlo analysis of base material decomposed CT projections,” in 2006 IEEE Nuclear Science Symposium Conference Record, vol. 5 (IEEE, 2006), pp. 3174–3175.
38. J. Nocedal and S. Wright, Numerical optimization (Springer Science & Business Media, 2006).
39. J. Nocedal, “Updating quasi-newton matrices with limited storage,” Math. Comp. 35(151), 773 (1980). [CrossRef]
40. L. Li, Z. Chen, W. Cong, and G. Wang, “Spectral CT modeling and reconstruction with hybrid detectors in dynamic-threshold-based counting and integrating modes,” IEEE Trans. Med. Imaging 34(3), 716–728 (2015). [CrossRef]
41. J. Hubbell, “Photon cross sections, attenuation coefficients and energy absorption coefficients,” Natl. Bureau Standards Rep. NSRDS-NBS29, Wash. DC (1969).
42. S. Gu, L. Zhang, W. Zuo, and X. Feng, “Weighted nuclear norm minimization with application to image denoising,” in Proceedings of the IEEE conference on computer vision and pattern recognition, (2014), pp. 2862–2869.
43. L. I. Rudin, S. Osher, and E. Fatemi, Nonlinear total variation based noise removal algorithms, (Elsevier North-Holland, Inc., USA, 1992), pp. 259–268.
44. S. Osher, M. Burger, D. Goldfarb, J. Xu, and W. Yin, “An iterative regularization method for total variation-based image restoration,” Multiscale Model. Simul. 4(2), 460–489 (2005). [CrossRef]
45. D. P. Bertsekas and J. N. Tsitsiklis, “Parallel and distributed computation: numerical methods,” (2003).
46. M. Nikolova, “An algorithm for total variation minimization and applications,” J. Math. Imaging Vis. 20(1/2), 163–177 (2004). [CrossRef]
47. X. Yuan, “Generalized alternating projection based total variation minimization for compressive sensing,” in 2016 IEEE International Conference on Image Processing (ICIP), (IEEE, 2016), pp. 2539–2543.
48. W. van Aarle, W. J. Palenstijn, J. Cant, E. Janssens, F. Bleichrodt, A. Dabravolski, J. D. Beenhouwer, K. J. Batenburg, and J. Sijbers, “Fast and flexible x-ray tomography using the astra toolbox,” Opt. Express 24(22), 25129–25147 (2016). [CrossRef]
49. Y. Zhang, Y. Xi, Q. Yang, W. Cong, J. Zhou, and G. Wang, “Spectral CT reconstruction with image sparsity and spectral mean,” IEEE Trans. Comput. Imaging 2(4), 510–523 (2016). [CrossRef]
50. J. Hubbell and S. Seltzer, “Radiation and biomolecular physics division, MPL NIST,” National Institute of Standards and Technology (2008) [accessed 19 July 2019], https://physics.nist.gov/PhysRefData/Xcom/html/xcom1-t.html
51. S. Luc, H. Alexandre, A. Vincent, C. Arnaud, F. Jean-Baptiste, M. Johan, O. Anne-Blandine, B. Mourad, and M. Jacques, “3d image reconstruction for comparison of algorithm database: a patient-specific anatomical and medical image database,” IRCAD, Strasbourg, France, Tech. Rep (2010).
52. T. Mao, A. P. Cuadros, X. Ma, W. He, Q. Chen, and G. R. Arce, “Coded aperture optimization in X-ray tomography via sparse principal component analysis,” IEEE Trans. Comput. Imaging 6, 73–86 (2020). [CrossRef]
53. A. Maier, H. G. Hofmann, M. Berger, P. Fischer, C. Schwemmer, H. Wu, K. Müller, J. Hornegger, J.-H. Choi, C. Riess, A. Keil, and R. Fahrig, “CONRAD–A software framework for cone-beam imaging in radiology,” Med. Phys. 40(11), 111914 (2013). [CrossRef]
54. A. P. Cuadros, C. Peitsch, H. Arguello, and G. R. Arce, “Coded aperture optimization for compressive X-ray tomosynthesis,” Opt. Express 23(25), 32788–32802 (2015). [CrossRef]
55. A. P. Cuadros and G. R. Arce, “Coded aperture compressive X-ray spectral CT,” in 2017 International Conference on Sampling Theory and Applications (SampTA), (IEEE, 2017), pp. 548–551.
