Monochromatic image reconstruction via machine learning

In medical CT, the x-ray source generally emits polychromatic x-ray photons, and the x-ray detector array operates in the energy-integrating mode to acquire energy intensity. The x-ray linear attenuation coefficient depends on both material composition of the imaged object and the x-ray photon energy. If an x-ray beam passes through an object, the x-ray transmitted beam intensity $I\left( l \right)$ is accurately described by the non-linear integral model based on Beer–Lambert law [1, 8]:

$$\begin{equation}I\left( l \right) = \int {S\left( E \right)D} \left( E \right)\exp \left( { - \int\limits_l {\mu \left( {r,E} \right){\text{d}}r} } \right){\text{d}}E,\end{equation} \tag{ 1 }$$

where $S\left( E \right)$ is the energy spectrum of the x-ray source, $D\left( E \right)$ is the detection efficiency, and $\mu \left( {r,E} \right)$ is the energy-dependent linear attenuation coefficient at an energy $E$ and a spatial position $r$ along a linear path $l$ through the object. Using the integral mean value theorem, there is an energy level ${\varepsilon _l}$ for the x-ray transmission path such that the following formula holds:

$$\begin{equation}\left\{ \begin{gathered} I\left( l \right) = {I_0}\left( l \right)\exp \left( { - \int\limits_l {\mu \left( {r,{\varepsilon _l}} \right){\text{d}}r} } \right) \hfill \\ {I_0}\left( l \right) = \int\limits_{{E_{\min }}}^{{E_{\max }}} {S\left( E \right)D\left( E \right){\text{d}}E} \hfill \\ \end{gathered} \right..\end{equation} \tag{ 2 }$$

Equation (2) is equivalent to

$$\begin{equation}\int\limits_l {\mu \left( {r,{\varepsilon _l}} \right){\text{d}}r} = \log \left[ {\frac{{{I_0}\left( l \right)}}{{I\left( l \right)}}} \right]{} ,\end{equation} \tag{ 3 }$$

where ${I_0}\left( l \right)$ is the x-ray intensity along the path $l$ without any object in the field of view. Equation (3) indicates that the line integral along a transmission path $l$ relies on energy levels ${\varepsilon _l}$. Because lower energy photons are more attenuated by tissues than higher energy photons, the x-ray spectral distribution on different x-ray transmission paths may be different. Thus, the line integral along different transmission paths may correspond to different energy levels in equation (3). Hence, the x-ray physical model described in equation (3) is different from the Radon transform without energy dependence. The mismatch of the physical model may induce beam-hardening artifacts in the image reconstruction [2].

The purpose of this research is to establish a transform relation between polychromatic CT images and monochromatic sinograms at a pre-determined energy level $\varepsilon$ in the detectable energy range, using a ML technique. A polychromatic attenuation image $\mu \left( r \right)$ is reconstructed from log-transformed raw data collected at all x-ray transmission paths and viewing angles. The image $\mu \left( r \right)$ contains abundant prior information on the object structure and x-ray attenuation information. A practical method is to map attenuation distribution on the x-ray transmission path $l$ in the polychromatic image $\mu \left( r \right)$ to corresponding monochromatic projection datum at the specific energy level $\varepsilon$, denoted by ${p_{{\text{mon}}}}\left( {l,\varepsilon } \right)$. This transform relation can be described mathematically by an optimization model as follows [25]:

$$\begin{equation}\mathscr{M} = {\text{arg}}\min \sum\limits_{l\; \in \;{\text{all}}\;{\text{path}}\;{\text{set}}} {\left\| {\mathscr{M}\left( {\mu \left( r \right),\;r \in l} \right) - {p_{{\text{mon}}}}\left( {l,\varepsilon } \right)} \right\|} ,\end{equation} \tag{ 4 }$$

where $\mathscr{M}$ denotes a transformation relation, $\mu \left( r \right)$ is the polychromatic attenuation image reconstructed from a single-spectrum CT, $l$ is an x-ray transmission path indexed by the location variable $r$, and ${p_{{\text{mon}}}}\left( {l,\varepsilon } \right)$ is the monochromatic projection data at a pre-specified energy level $\varepsilon$, which can be calculated from the corresponding labeled monochromatic CT image at the energy level $\varepsilon$ using the standard ray-tracing method. Specifically, we assume that all x-ray transmission paths have the same length, equal to the diameter of the field of view. The transmission path is divided into a fixed number of segments, which is denoted as n. Typically, n is chosen larger than $\sqrt 2$× the input image size for sufficient sampling. The attenuation coefficient of each segment is calculated from an input polychromatic CT image using bi-linear interpolation, as shown in figure 1. The involved computation is fairly straightforward, identical to the ray-tracing process in a typical iterative reconstruction, such as algebraic reconstruction technique (ART) or simultaneous algebraic reconstruction technique (SART). As input data, the attenuation distribution along the x-ray transmission path is calculated by multiplying the attenuation coefficient of each segment by the length of the segment.

Zoom In Zoom Out Reset image size

A well-defined fully connected neural network is capable of learning any function [30]. To implement the optimization model in equation (4), we establish a one-dimensional fully connected deep network, allowing a much lower computational cost and much less memory requirement than higher dimensional networks. The neural network consists of seven layers, including an input layer, five hidden layers, an output layer, and a shortcut connection, as shown in figure 2. The first hidden layer contains 1024 neurons, the second hidden layer contains 512 neurons, the third hidden layer has 256 neurons, the fourth hidden layer contains 128 neurons, and the fifth hidden layer has 64 neurons. For every hidden layer, neurons receive the weighted combinations of neuron values on the previous layer and perform corresponding sigmoid activations. Inspired by the idea of the two-dimensional ResNet [31], the shortcut connection technique is adapted to stack the line integral of the input polychromatic image along the transmission path to the neuron on the output layer to implement a residual mapping. For implementation of the shortcut connection, a summation of the attenuation distribution along the transmission path on the input polychromatic CT image can reproduce an approximate energy-integrating projection value $\left\{ {\log \left[ {{I_0}\left( l \right)/I\left( l \right)} \right]} \right\}$. The linear integral model is approximate in CT imaging, and usually leads to inaccurate quantification of reconstructed images. The proposed VM imaging approach is to correct the errors/mismatches of the line integral model, which can be effectively implemented by the residual mapping scheme. Prior works show that modeling residual mapping is easier than the original mapping [31]. The shortcut connection helps train neural networks efficiently. The input of the network is attenuation distribution along an x-ray transmission path $l$ on the polychromatic CT image $\mu \left( r \right)$. The x-ray transmission path is divided into 1024 segments to describe the corresponding attenuation coefficient distribution. The output layer only contains a single neuron yielding a weighted linear combination of all neuron units on the last hidden layer, and outputs an x-ray monochromatic projection along the corresponding x-ray transmission path. The loss function uses the l₁ norm to evaluate the difference between the predicted monochromatic projection value and the monochromatic projection value along the same transmission path for the label image.

Zoom In Zoom Out Reset image size

A well-trained fully connected deep network would be a non-linear transform $\mathscr{M}$ that maps the polychromatic CT images to the monochromatic projection data at energy level $\varepsilon$. Once the monochromatic projections are obtained through ML, a monochromatic image $\mu \left( {r,\varepsilon } \right)$ can be reconstructed using a standard image reconstruction algorithm such as a filtered backprojection (FBP) or iterative algorithm (ART or SART) based on following line integral model:

$$\begin{equation}\int\limits_l {\mu \left( {r,\varepsilon } \right){\text{d}}r} = {p_{{\text{mon}}}}\left( {l,\varepsilon } \right).\end{equation} \tag{ 5 }$$

Equation (5) is an accurate x-ray imaging model at energy level $\varepsilon$. Therefore, the monochromatic CT image can be reconstructed based on equation (5), achieving accurate quantification of attenuation images and overcoming beam-hardening artifacts.

We performed an ablation study on the number of neurons per layer, the number of layers in the network, activation functions, and so on. The proposed network architecture has five hidden layers, where the number of neurons in each layer is half the number of neurons in the previous layer. This network was selected based on our extensive simulation and experimental results, some of which are described below. To demonstrate the performance of the proposed network, we also trained a multilayer perceptron (MLP) with the same number of neurons in each of the five hidden layers, and evaluated its convergence and accuracy. We found that our proposed network showed an excellent convergence and accuracy, comparable to that of the MLP, while the proposed network only requires a half the memory and half the computational cost. Especially, compared with the network architecture without a shortcut connection, the convergence speed of the network architecture with a shortcut connection was increased by more than three times. Moreover, we applied different activation functions for our neural network. The experiments show that the sigmoid is an appropriate activation function for this application. For example, if we use the ReLU activation function for all the hidden layers, the average relative error of the reconstructed monochromatic projection would exceed 20%, which is ten times higher than that of the proposed network architecture. The proposed network represents an excellent balance between accuracy, efficiency, and robustness.

A dataset of VM CT images produced from a GE Discovery CT750HD dual-energy scanner at Ruijin Hospital in Shanghai was used for the network training, validation, and testing. The scanning was performed with fast kVp switching between 80 kVp and 140 kVp with tube current of 260 mAs. The gantry rotation time was 0.5 s. The dataset includes 274 VM CT images at 50 keV, 60 keV, 65 keV, 70 keV, 80 keV, 90 keV, 100 keV, and 110 keV, respectively. Based on Beer–Lambert law, we first synthesized polychromatic projection data from these multi-energy images using an x-ray source spectral distribution at 120 kVp/30 mA generated by public software [32]. Then, 274 polychromatic CT images were reconstructed from polychromatic projection data using the FBP algorithm. The 274 polychromatic CT images and corresponding 274 VM images at 80 keV formed the training dataset I, while 274 polychromatic CT images and corresponding 274 VM images at 110 keV formed the training dataset II. The dataset was divided into training dataset of 200 images, validation data of 50 images, and testing data of 24 images. A total of 105 million (200 images × 720 views × 729 rays) data pairs as training data were extracted along x-ray line integral paths from 200 polychromatic CT images and corresponding 200 VM images. For the one-dimensional model, 105 million data pairs were sufficient for the network training with supervised learning.

The training procedure was programmed in Python and Tensorflow on a PC with a NVIDIA Titan XP GPU with 12 GB memory. The network training was conducted using the ADAM optimization algorithm with β₁= 0.9, β₂ = 0.999, epsilon = 1.0 × 10⁻⁷, and weight decay = 0.0. The learning rate was set to 10⁻³. Weight and bias parameters of the network were initialized randomly. Data were randomly sampled in the training dataset, maximizing the probability of finding the global minimum. Based on dataset I, the network is trained to generate monochromatic sinogram at 80 keV, while the network is trained using dataset II to generate monochromatic sinogram at 110 keV. The networks were trained in 1000 epochs within 12 h. Figure 3(a) shows the loss function during the network training using dataset I. Figure 3(b) shows the loss function during the network training using dataset II. The training of the fully connected deep network showed an excellent convergence behavior.

Zoom In Zoom Out Reset image size

For the testing of the trained network, 24 polychromatic images at 120 kVp were input to the trained neural network models I and II to generate monochromatic sinogram at 80 keV and 110 keV, respectively. CT image were randomly selected in testing dataset as examples to present the quality of monochromatic imaging. The labeled monochromatic sinogram data were obtained from the labeled monochromatic CT image using the ray-tracing method. Figure 4 presents the comparison between the estimated monochromatic sinograms and the corresponding (ground-truth) labeled monochromatic sinograms at 80 keV and 110 keV respectively. Figure 5 presents the comparison between the polychromatic x-ray projection at 120 kVp, the estimated monochromatic projection, and the corresponding labeled monochromatic projection at the horizontal views. The trained neural network delivered high-quality monochromatic projection data in the testing phase, with an average relative error of less than 2%. With VM images reconstructed from DECT as the reference, we used the popular signal-to-noise ratio (PSNR) and structural similarity (SSIM) indices to evaluate the estimated monochromatic sinogram against the labeled monochromatic sinogram. The average PSNR measures are 37.54 dB and 36.89 dB for monochromatic sinograms at 80 keV and 110 keV respectively, and the average SSIM values are 0.9941 and 0.9902 for monochromatic sinograms at 80 keV and 110 keV respectively.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Furthermore, the VM CT images at 80 keV and 110 keV were reconstructed from the estimated monochromatic sinogram at 80 keV and 110 keV using the FBP algorithm. Figure 6 presents the comparison between the estimated VM image and the label VM image at 80 keV and 110 keV. We calculated PSNR and SSIM indices to evaluate the estimated VM images at 80 keV and 110 keV. The average PSNR was 69.79 dB and 63.34 dB for the VM images at 80 keV and 110 keV, respectively, while the average SSIM was 0.9999 and 0.997 for the VM images at 80 keV and 110 keV, respectively. The proposed ML-based method well preserved structural information especially texture features and gave superior image quality.

Zoom In Zoom Out Reset image size

We obtained a second series of clinical DECT dataset to further evaluate the performance of the deep learning-based VM imaging method. We now used the originally reconstructed kVp images as the input to the network, instead of the synthesized kVp images used in the previous section. DECT data of eight patients (3182 slices in total) were collected on a SOMATOM Definition Flash DECT scanner (Siemens Healthineers, Forchheim, Germany) at Ruijin Hospital in Shanghai, China. The DECT scanner worked in the dual-source scanning mode, which operated at 100 kVp/210 mAs and 140 kVp/162 mAs with a wedge filter and a flat filter respectively. The scanning was set for exposure time 0.5 s. CT images were reconstructed using the FBP algorithm. The dataset includes 3182 polychromatic CT images at 140 kVp, being associated with 3182 VM CT images at 80 keV and 3182 VM CT images at 110 keV. The dataset was split into training, validation and testing sets, which respectively came from five, two and one patients. In other words, 2195 polychromatic CT images at 140 kVp and the corresponding 2195 monochromatic images at 80 keV formed training dataset I, while 2195 polychromatic CT images at 140 kVp and the corresponding 2195 monochromatic images at 110 keV formed training dataset II. In total, 1152 million (2195 images × 720 views × 729 rays) data pairs were extracted along x-ray line integral paths from the 2195 polychromatic CT images at 140 kVp and corresponding 2195 monochromatic images in the datasets I and II for training the network models I and II respectively.

Polychromatic CT images at 140 kVp for the test datasets were input to the trained neural network models I and II to generate monochromatic sinograms at 80 keV and 110 keV, respectively. Figure 7 presents the comparison between the estimated monochromatic sinograms and the corresponding label monochromatic sinograms at 80 keV and 110 keV respectively. Figure 8 presents the comparison between the polychromatic x-ray projection at 140 kVp, the estimated monochromatic projection, and corresponding label monochromatic projection in the horizontal views. We calculated PSNR and SSIM to evaluate the estimated monochromatic sinograms with label monochromatic sinograms as reference. The average PSNR was 33.40 dB and 33.27 dB for the estimated monochromatic sinograms at 80 keV and 110 keV, respectively, while average SSIM was 0.9995 and 0.9978 for the estimated monochromatic sinograms at 80 keV and 110 keV, respectively.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Again, VM images at 80 keV and 110 keV were reconstructed from the estimated monochromatic sinogram data at 80 keV and 110 keV respectively using the FBP algorithm. Figure 9 shows the comparison between the estimated VM image and the label VM images at 80 keV and 110 keV, respectively. The proposed estimation method well preserved structural information, in particular texture features, yielding a superior image quality. We calculated the PSNR and SSIM measures to evaluate the reconstructed VM images against the label VM images. The average PSNR was 48.92 dB and 48.73 dB for the reconstructed VM images at 80 keV and 110 keV respectively, while the average SSIM was 0.9940 and 0.9933 for the reconstructed VM images at 80 keV and 110 keV, respectively. Comparing image quality in terms of SSIM and PSNR in sections 3.1 and 3.2, the quality of the VM image reconstructed from real kVp images is slightly compromised relative to the quality of the VM image reconstructed from synthesized kVp images, due to non-ideal noise and spectral distortion in clinical data.

Zoom In Zoom Out Reset image size