Patient-specific reconstruction of volumetric computed tomography images from a single projection view via deep learning

Problem formulation.

We formulate the problem of 3D image reconstruction from 2D projection(s) into a deep-learning framework. Given a sequence of 2D projections denoted as $\left\{X_1, X_2, \cdots , X_N\right\}$, where $X_i\in {ℝ}^{H_{2D}\times W_{2D}}$ for all $1<i<N$ and $N$ is the number of given 2D projections, the goal is to generate a volumetric 3D image $Y$ describing the corresponding 3D physical scene. With the sequence of 2D projections as input, the deep-learning model outputs the predicted 3D volume denoted as $Y_{pred}\in {ℝ}^{C_{3D}\times H_{3D}\times W_{3D}}$, while $Y_{truth}\in {ℝ}^{C_{3D}\times H_{3D}\times W_{3D}}$ is the ground truth 3D image as the reconstruction target. Note that network prediction $Y_{pred}$ is of the same size as ground truth image $Y_{truth}$, where each entry is a voxel-wise intensity value. Thus, the problem is formulated as finding a mapping function $F$ transforming 2D projections to volumetric images. To tackle this problem, a deep-learning model is trained to find the mapping function $F$, which uses 2D projections $\left\{X_1, X_2, \cdots , X_N\right\}$ as input and predicts the corresponding 3D image $Y_{pred}$, as expressed in equation (1).

In order to use a sequential of 2D projections as model input, we stack all the 2D projections together as a single 3D tensor. In other words, a set of 2D projections $\left\{X_1,X_2,\cdots ,X_N\right\}$ $(X_i\in {ℝ}^{H_{2D}\times W_{2D}})$ are stacked as a 3D volume $\text{Z}\in {ℝ}^{N\times H_{2D}\times W_{2D}}$, where $N$ is the number of 2D projections. In what follows, we introduce the model architecture of the deep neural network in detail.

Encoder–decoder framework.

The deep neural network is formulated into an encoder–decoder framework (Fig. 2). In the auto-encoder model⁴⁷, the encoder converts high-dimensional data into embedded representations while the decoder reconstructs high-dimensional input. In our task, instead of decoding to get the input, we developed a modified decoder to generate the corresponding volumetric images based on the codes converted by the encoder. More precisely, with a sequence of 2D projections as input, the encoder network learns the feature representation by extracting semantic information from 2D projections, such as organ position and size. In this way, the encoder network learns a transformation function $h_1$ from the 2D image domain to the feature domain. A transformation module then follows to learn the manifold mapping function $h_2$ in the feature domain to transform the feature representation across dimensions. Using the learned feature representation as the input, the decoder network is trained to generate the 3D volume. In other words, the decoder network learns a transformation function $h_3$ from the feature domain to the 3D image domain. In this way, we fit the target mapping function $F$ by decomposition: $F=h_1\circ h_2\circ h_3$. The rationale behind our network design is that both 2D projections and 3D images should share the same semantic feature representation in the feature domain, as they represent the image expressions of the same object or physical scene. Accordingly, the representation in feature space should remain invariant. In a sense, if the model can learn the transformation function between the feature space and the 2D/3D image space, it is possible to reconstruct 3D images from 2D projections. Therefore, following this encoder–decoder framework, our model is able to learn how to generate 3D images from 2D projections by leveraging from the learned representations in high-dimensional feature space.

Representation network.

Superb performance has been achieved by deep residual networks (such as ResNet)⁴⁸ in many tasks. A key step in residual learning is the identity mapping that facilitates the training process and avoids gradient vanish in back-propagation⁴⁸, which encourages residual learning of the hierarchical representation at each stage and eases the training of the deep network. Motivated by this feature, we introduce a residual-learning scheme in the representation network (Fig. 2), where the 2D convolution residual block is used to assist the deep model to learn semantic representations from 2D projections. More details about the residual learning scheme are presented in the “Ablative study and discussion” section of the Supplementary Information, and the results are summarized in Supplementary Table 2. Specifically, each 2D convolution residual block consists of a pattern of “2D convolution layer (with kernel size 4 and stride 2) → 2D batch normalization layer → rectified linear unit activation (ReLU) layer → 2D convolution layer (with kernel size 3 and stride 1) → 2D batch normalization layer → ReLU layer”. The first layer conducts 2D convolution operations using a $4\times 4$ kernel with sliding stride $2\times 2$, which down-samples the spatial size of the feature map by a factor 2. In addition, to keep the sparsity of high-dimensional feature representation, we correspondingly doubled the channel number of the feature maps by increasing the number of convolutional filters. A distribution normalization layer among the training mini-batch (batch normalization)⁴⁹ then follows before feeding the feature maps through the ReLU layer⁵⁰. Next, the second 2D convolution layer and 2D batch normalization layer are done by a kernel size of $3\times 3$ and sliding stride $1\times 1$, which keeps the spatial shape of the feature maps. Moreover, before applying the second ReLU layer, an extra shortcut path is established to add up the output of the first convolution layer to obtain the final output. By setting up the shortcut path of identity mapping, the second convolution layer is encouraged to learn the residual feature representations. In order to extract hierarchical semantic features from 2D projections, we constructed the representation network by concatenating five 2D convolution residual blocks with different number of convolutional filters. A detailed discussion of the network depth is available in the “Ablative study and discussion” section of the Supplementary Information, with some results illustrated in Supplementary Fig. 8. To be concise, we use the notation of $k\times m\times n$ to denote $k$ channels of feature maps in a spatial size of $m\times n$. In the generation network, the size of input images is denoted as $N\times 128\times 128$, where $N$ is the number of 2D projections. The change of feature map size through the network are: $N\times 128\times 128\to 256\times 64\times 64\to 512\times 32\times 32\to 1024\times 16\times 16\to 2048\times 8\times 8\to 4096\times 4\times 4$, where each $\to$ means going through a 2D convolution residual block as described above, except that batch normalization and ReLU activation are removed in the first convolution layer. Thus, the output of the representation network is a feature representation extracted from 2D projections with a size of $4096\times 4\times 4$.

Transformation module.

To bridge the representation and generation networks, a transformation module is deployed after learning the representations. As shown in Fig. 2, by taking the convolution operations with a kernel size of $1\times 1$ and ReLU activation, the 2D convolution layer learns a transformation across all 2D feature maps. Then, we reshape embedded representations from $4096\times 4\times 4$ to $2048\times 2\times 4\times 4$. In this way, we transform the feature representation across dimensions for subsequent 3D volume generation. Next, a 3D deconvolution layer with a kernel size of $1\times 1\times 1$ and sliding stride of $1\times 1\times 1$ learns a transformation among all 3D feature cubes while keeping the feature size unchanged. This transformation module bridges the 2D and 3D feature spaces. Moreover, as described in previous work⁵¹, we also remove the batch normalization in the transformation module to help the knowledge transfer through this module.

Generation network.

The generation network was built upon the 3D deconvolution block, which consists of a pattern of “3D deconvolution layer (with kernel size 4 and stride 2) → 3D batch normalization layer → ReLU layer → 3D deconvolution layer (with kernel size 3 and stride 1) → 3D batch normalization layer → ReLU layer”. Note that the “deconvolution” layer actually means the operation of “transformed convolution” or fractional stride convolution which performs up-sampling operation. The first deconvolution layer up-samples a feature map by a factor 2 with a $4\times 4\times 4$ kernel and sliding stride $2\times 2\times 2$. In order to transform from a high-dimensional feature domain to a 3D-image domain, we accordingly reduce the number of feature maps by decreasing the number of deconvolutional filters. The second deconvolution layer with a $3\times 3\times 3$ kernel and sliding stride $1\times 1\times 1$ keeps the spatial shape of feature maps. A 3D batch normalization layer and a ReLU layer follow after each deconvolution layer. Hierarchically, the 3D generation network consists of 5 concatenating deconvolution residual blocks. Following the same convention as in the representation network, we use a notation of $k\times m\times n\times p$ to denote $k$ channels of 3D feature maps with a spatial size of $m\times n\times p$. With a representation input of $2048\times 2\times 4\times 4$, the data flow of the feature maps is: $2048\times 2\times 4\times 4\to 1024\times 4\times 8\times 8\to 512\times 8\times 16\times 16\to 256\times 16\times 32\times 32\to 128\times 32\times 64\times 64\to 64\times 64\times 128\times 128$, where each $\to$ denotes a 3D residual block. At the end of the generation network, an output transformation module was constructed with a 3D convolution layer and a 2D convolution layer with a kernel size of $1$, which outputs 3D images fitting the shape of the reconstructed images. Finally, the generation network outputs the predicted 3D images of size $C_{3D}\times 128\times 128$, where $C_{3D}$ is the size of the target volumetric images along z-axis. Note that in the output transformation module, the batch normalization layer is removed, and that there is no ReLU layer after the final convolution layer.

Materials.

The approach is evaluated by using three cases of different disease sites. In the first study, a 10-phase upper abdominal 4D-CT scan of a patient for radiation therapy (RT) treatment planning is selected. To proceed, the first 6 phases are used to generate the CT-DRR pairs for model training and validation with the procedure described above. We use the anterior-posterior (AP) 2D projection as input (Fig. 2). With translation, rotation and deformation introduced to the CT volume, we obtain a total of 720 DRRs representing different scenarios of the patient anatomy for model training and 180 DRRs for validation. To ensure that the testing data are not seen in the model-training process, we generate 600 testing DRR samples independently from the remaining 4 phases of the 4D-CT. The 4D-CT images are acquired with 120 kV, 80 mA on a Positron emission tomography–computed tomography (PET-CT) simulator (Biograph mCT 128, Siemens Medical Solutions, Erlangen, Germany) together with a Varian Real-time Position Management™ (RPM) system (Varian Medical Systems, Palo Alto, CA). The 2D projection data are obtained by projecting each of the 3D CT data in the geometry of the on-board imager of TrueBeam™ system (Varian Medical System, Palo Alto, CA). In the second experiment, a lung cancer patient is chosen with two independent treatment planning 4D CT scans acquired at two different times with the same imaging parameter settings as above. Using the data-augmentation strategy described above, the first 4D-CT is used to generate training (2,400 samples) and validation (600 samples) datasets, whereas the second 4D CT was used to generate a testing dataset (200 samples). For each of the images in the training and testing datasets, the corresponding 2D projections are produced by projecting the 3D-CT volume in the geometry of the on-board imager of the TrueBeam™ system. To build a reliable model, the training and testing datasets might come from the same data distribution, but the datasets are independently sampled. The data acquisition and processing for the head-and-neck case are described in Supplementary Information.

Image preprocessing.

Data are preprocessed before feeding them into the network. First, we resize all data samples to the same size. For example, all the 2D projection images are reshaped to $128\times 128$. The volumetric images of the abdominal-CT and the lung-CT are resized to $46\times 128\times 128$ and $168\times 128\times 128$ respectively, due to their different depths in the z-axis. Each data sample is a pair of 2D projected view(s) and the corresponding 3D-CT. Similar to other deep-learning-based imaging studies¹⁵, down-sampling is introduced purely because of the memory limitation and for the purpose of computational efficiency. The formulation and algorithm are scalable to full-size images (512×512), because the component layers used in our model are also scalable to images of different sizes. At the current resolution of 128×128 (which is the same as that used in the deep-learning-based MRI reconstruction¹⁵), small motions of less than 3 mm may not be described accurately. However, we should emphasize that this resolution does not represent a fundamental limit of the deep-learning-based approach and can be improved as computational technology advances. In practice, methods such as deep-learning-based super-resolution are being actively pursued, which may be employed to improve the spatial resolution of the approach. Additionally, following the standard protocol of data pre-processing, we conduct scaling normalization for both the 2D projections and the 3D volumetric images, where pixel-wise or voxel-wise intensities are normalized to the interval $\left[0, 1\right]$. Moreover, we normalize the statistical distribution of the pixel-wise intensity values in the input 2D projections to be closer to a standard Gaussian distribution $\mathcal{N}\left(0,1\right)$. Specifically, we calculate the statistical mean and standard derivation among all the training data. When a new sample was inputted, we subtract the mean value from the input image(s) and divide the image(s) by the standard derivation to get the input 2D image(s).

Training details.

With input images $X$ containing a stacked sequence of 2D projections, we train the deep network to predict the volumetric images $Y_{pred}$, which is expected to be as close as possible to the ground-truth images $Y_{truth}$. We define the cost function as the mean squared error (MSE) between the prediction $Y_{pred}$ and the ground truth $Y_{truth}$, and the model was optimized by stochastic gradient descent iteratively. For comparison, we used the same training strategy and hyper-parameters for all experiments. We implemented the network by using the PyTorch⁵² library, and used the Adam optimizer⁵³ to minimize the loss function and to update the network parameters iteratively through back-propagation. A learning rate of 0.00002 and a mini-batch size of 1 are used because of memory limitations. At the end of each training epoch, the model is evaluated on the validation set. This strategy is commonly used to monitor the model performance and avoid overfitting the training data. In addition, the learning rate is scheduled to decay according to the validation loss. Specifically, if the validation loss remains unchanged for 10 epochs, the learning rate is reduced by a factor 2. Finally, the best checkpoint model with the smallest validation loss is saved as the final model in the experiments. We trained the network using one Nvidia Tesla V100 graphics processing unit (GPU) for 100 epochs (duration typically around 20 hours for the abdominal-CT case). During testing, the typical inference time for 3D reconstruction of one testing sample is around 0.5 seconds.

Evaluation.

To evaluate the performance of the approach, we deploy the trained model on a testing dataset, and analyze the reconstruction results using both qualitative and quantitative evaluation metrics. We use four different metrics to measure the quality of predicted 3D images: mean absolute error (MAE), root mean squared error (RMSE), structural similarity (SSIM)⁵⁴, and peak signal-to-noise-ratio (PSNR). We compute the average values across all testing samples, and they are shown in Table 1. MAE / MSE is the L1-norm / L2-norm error between $Y_{pred}$ and $Y_{truth}$. As usual, we take the square root of MSE to get RMSE. In practice, MAE and RMSE are commonly used to estimate the difference between the prediction and ground-truth images. SSIM score is calculated with a windowing approach in an image, and is used for measuring the overall similarity between two images. In general, a lower value of MAE and RMSE or a higher SSIM score indicates a better prediction closer to the ground-truth images. PSNR is defined as the ratio between the maximum signal power and the noise power that affects the image quality. PSNR is widely used to measure the quality of image reconstruction.

Comparison study.

To better benchmark the proposed method against the existing techniques, we conduct a comparative study with the published principal component analysis (PCA)-based method^55–57 and elaborate the difference and advantages of our proposed approach. The comparison is done for a special situation of 4D-CT reconstruction (abdominal CT) where the anatomical motion may be characterized by principal components. We find that the PCA and deep learning-based methods produce similar results in an ideal case when there is no inter-scan variation in patient positioning (since the results are very similar, the resultant images are not shown). However, the deep learning model outperforms the PCA method in more realistic scenarios when the patient position deviates slightly from that of the reference scan (see Supplementary Information for details).

Reporting summary.

Further information on research design is available in the Nature Research Reporting Summary linked to this article.