Solving Low-Dose CT Reconstruction via GAN with Local Coherence

Traditional generative adversarial network [9] consists of two main modules, a generator and a discriminator. The generator $\mathcal{G}$ is a mapping from a latent-space Gaussian distribution $\mathbb{P}_{Z}$ to the synthetic sample distribution $\mathbb{P}_{X_{G}}$, which is expected to be close to the real sample distribution $\mathbb{P}_{X}$. On the other hand, the discriminator $\mathcal{D}$ aims to maximize the distance between the distributions $\mathbb{P}_{X_{G}}$ and $\mathbb{P}_{X}$. The game between the generator and discriminator actually is an adversarial process, where the overall optimization objective follows a min-max principle:

As mentioned in Section 1, optical flow can capture the temporal coherence of object movements, which plays a crucial role in many video-related tasks. More specifically, the optical flow refers to the instantaneous velocity of pixels of moving objects on consecutive frames over a short period of time [3]. The main idea relies on the practical assumptions that the brightness of the object more likely remains stable across consecutive frames, and the brightness of the pixels in a local region are usually changed consistently [10]. Based on these assumptions, the brightness of optical flow can be described by the following equation:

where $v=(v_{w},v_{h})$ represents the optical flow of the position $(w,h)$ in the image. $\nabla I=(\nabla I_{w},\nabla I_{h})$ denotes spatial gradients of image brightness, and $\nabla I_{t}$ denotes the temporal partial derivative of the corresponding region.

Following the equation (4), we consider the question that whether the optical flow idea can be applied to 3D CT reconstruction. In practice, the brightness of adjacent CT images often has very tiny difference, due to the inherent continuity and structural integrity of human body. Therefore, we introduce the “local coherence” that indicates the correlation between adjacent images of a tissue. Namely, adjacent CT images often exhibit significant similarities within a certain local range along the vertical axis of the human body. Due to the local coherence, the noticeable variations observed in CT slices within the local range often occur at the edges of organs. We can substitute the temporal partial derivative $\nabla I_{t}$ by the vertical axial partial derivative $\nabla I_{z}$ in the equation (4), where “z” indicates the index of the vertical axis. As illustrated in Figures 1, the local coherence can be captured by the optical flow between adjacent CT slices.

The proposed framework comprises three components, including a generator $\mathcal{G}$, a discriminator $\mathcal{D}$ and an optical flow estimator $\mathcal{F}$. The generator is the core component, and the flow estimator provides auxiliary warping images for the generation process.

Suppose we have a sequence of measurements $\mathbf{y_{1}},\mathbf{y_{2}},\cdots,\mathbf{y_{n}}$; for each $\mathbf{y_{i}}$, $1\leq i\leq n$, we want to reconstruct its ground truth image $\mathbf{x_{i}^{r}}$ as the equation (1). Before performing the reconstruction in the generator $\mathcal{G}$, we apply some prior knowledge in physics and run filter backward projection on the measurement $\mathbf{y_{i}}$ in equation (1) to obtain an initial recovery solution $\mathbf{s_{i}}$. Usually $\mathbf{s_{i}}$ contains significant noise comparing with the ground truth $\mathbf{x_{i}^{r}}$. Then the network has two input components, i.e., the initial backward projected image $\mathbf{s_{i}}$ that serves as an approximation of the ground truth $\mathbf{x_{i}^{r}}$, and a set of neighbor CT slices $\mathcal{N}(\mathbf{s_{i}})=\{\mathbf{s_{i-1}},\mathbf{s_{i+1}}\}$ ¹¹1If $i=1$, $\mathcal{N}(\mathbf{s_{i}})=\{\mathbf{s_{2}}\}$; if $i=n$, $\mathcal{N}(\mathbf{s_{i}})=\{\mathbf{s_{n-1}}\}$. for preserving the local coherence. The overall structure of our framework is shown in Figure 2. Below, we introduce the three key parts of our framework separately.

The optical flow $\mathcal{F}(\mathcal{N}(\mathbf{s_{i}}),\mathbf{s_{i}})$ denotes the brightness changes of pixels from $\mathcal{N}(\mathbf{s_{i}})$ to $\mathbf{s_{i}}$, where it captures their local coherence. The estimator is derived by the network architecture of FlowNet [8]. The FlowNet is an autoencoder architecture with extraction of features of two input frames to learn the corresponding flow, which is consist of $6$ (de)convolutional layers for both encoder and decoder.

The discriminator $\mathcal{D}$ assigns the label “$1$” to real standard-dose CT images and “$0$” to generated images. The goal of $\mathcal{D}$ is to maximize the separation between the distributions of real images and generated images:

where $\mathbf{x_{i}^{g}}$ is the image generated by $\mathcal{G}$ (the formal definition for $\mathbf{x_{i}^{g}}$ will be introduced below). The discriminator includes 3 residual blocks, with 4 convolutional layers in each residual block.

We use the generator $\mathcal{G}$ to reconstruct the high-quality CT image for the ground truth $\mathbf{x_{i}^{r}}$ from the low-dose image $\mathbf{s_{i}}$. The generated image is obtained by

where $\mathcal{W}(\cdot)$ is the warping operator. Before generating $\mathbf{x_{i}^{g}}$, $\mathcal{N}(\mathbf{x_{i}^{g}})$ is reconstructed from $\mathcal{N}(\mathbf{s_{i}})$ by the generator without considering local coherence. Subsequently, according to the optical flow $\mathcal{F}(\mathcal{N}(\mathbf{s_{i}}),\mathbf{s_{i}})$, we warp the reconstructed images $\mathcal{N}(\mathbf{x_{i}^{g}})$ to align with the current slice by adjusting the brightness values. The warping operator $\mathcal{W}$ utilizes bi-linear interpolation to obtain $\mathcal{W}(\mathcal{N}(\mathbf{x_{i}^{g}}))$, which enables the model to capture subtle variations in the tissue from the generated $\mathcal{N}(\mathbf{x_{i}^{g}})$; also, the warping operator can reduce the influence of artifacts for the reconstruction. Finally, $\mathbf{x_{i}^{g}}$ is generated by combining $\mathbf{s_{i}}$ and $\mathcal{W}(\mathcal{N}(\mathbf{x_{i}^{g}}))$. Since $\mathbf{x_{i}^{r}}$ is our target for reconstruction in the $i$-th batch, we consider the difference between $\mathbf{x_{i}^{g}}$ and $\mathbf{x_{i}^{r}}$ in the loss. Our generator is mainly based on the network architecture of Unet [26]. Partly inspired by the loss in [2], the optimization objective of the generator $\mathcal{G}$ comprises three items with the coefficients $\lambda_{\mathtt{pix}},\lambda_{\mathtt{adv}},\lambda_{\mathtt{per}}\in(0,1)$:

In (7), “$\mathcal{L}_{\mathtt{pixel}}$” is the loss measuring the pixel-wise mean square error of the generated image $\mathbf{x_{i}^{g}}$ with respect to the ground-truth $\mathbf{x_{i}^{r}}$. “$\mathcal{L}_{\mathtt{adv}}$” represents the adversarial loss of the discriminator $\mathcal{D}$, which is designed to minimize the distance between the generated standard-dose CT image distribution $\mathbb{P}_{X_{G}}$ and the real standard-dose CT image distribution $\mathbb{P}_{X}$. “$\mathcal{L}_{\mathtt{percept}}$” denotes the perceptual loss, which quantifies the dissimilarity between the feature maps of $\mathbf{x_{i}^{r}}$ and $\mathbf{x_{i}^{g}}$; the feature maps denote the feature representation extracted from the hidden layers in the discriminator $\mathcal{D}$ (suppose there are $t$ hidden layers):

where $\mathcal{D}_{j}(\cdot)$ refers to the feature extraction performed on the $j$-th hidden layer. Through capturing the high frequency differences in CT images, $\mathcal{L}_{\mathtt{percept}}$ can enhance the sharpness for edges and increase the contrast for the reconstructed images. $\mathcal{L}_{\mathtt{pixel}}$ and $\mathcal{L}_{\mathtt{adv}}$ are designed to recover global structure, and $\mathcal{L}_{\mathtt{percept}}$ is utilized to incorporate additional texture details into the reconstruction process.

First, our proposed approaches are evaluated on the “Mayo-Clinic low-dose CT Grand Challenge” (Mayo-Clinic) dataset of lung CT images [20]. The dataset contains 2250 two dimensional slices from 9 patients for training, and the remaining 128 slices from 1 patient are reserved for testing. The low-dose measurements are simulated by parallel-beam X-ray with 200 (or 150) uniform views, i.e., $N_{v}=200~{}(\text{or}~{}N_{v}=150)$, and 400 (or 300) detectors, i.e., $N_{d}=400~{}(\text{or}~{}N_{d}=300)$. In order to further verify the denoising ability of our approaches, we add the Gaussian noise with standard deviation $\sigma=2.0$ to the sinograms after X-ray projection in 50% of the experiments. To evaluate the generalization of our model, we also consider another dataset RIDER with non-small cell lung cancer under two CT scans [37] for testing. We randomly select 4 patients with 1827 slices from the dataset. The simulation process is identical to that of Mayo-Clinic. The proposed networks were implemented in the MindSpore framework and trained on Nvidia 3090 GPU with 100 epochs.

We consider several existing popular algorithms for comparison. (1) FBP [12]: the classical filter backward projection on low-dose sinograms. (2) FBPConvNet [11]: a direct inversion network followed by the CNN after initial FBP reconstruction. (3) LPD [1]: a deep learning method based on proximal primal-dual optimization. (4) UAR [22]: an end-to-end reconstruction method based on learning unrolled reconstruction operators and adversarial regularizers. Our proposed method is denoted by GAN-LC. We set $\lambda_{\mathtt{pix}}=1.0,\lambda_{\mathtt{adv}}=0.01~{}\text{and }~{}\lambda_{\mathtt{per}}=1.0$ for the optimization objective in equation (7) during our training process. Following most of the previous articles on 3D CT reconstruction, we evaluate the experimental performance by two metrics: the peak signal-to-noise ratio (PSNR) and the structural similarity index (SSIM) [33]. PSNR measures the pixel differences of two images, which is negatively correlated with mean square error. SSIM measures the structure similarity between two images, which is related to the variances of the input images. For both two meansures, the higher the better.

Table 1 presents the results on the Mayo-Clinic dataset, where the first row represents different parameter settings (i.e., the number of uniform views $N_{v}$, the number of detectors $N_{d}$ and the standard deviation of Gaussian noise $\sigma$) for simulating low-dose sinograms. Our proposed approach GAN-LC consistently outperforms the baselines under almost all the low-dose parameter settings. The methods FBP and UAR are very sensitive to noise; the performance of LPD is relatively stable but with low reconstruction accuracy. FBPConvNet has a similar increasing trend with our approach across different settings but has worse reconstruction quality. To evaluate the stability and generalization of our model and the baselines trained on Mayo-Clinic dataset, we also test them on the RIDER dataset. The results are shown in Table 2. Due to the bias in the datasets collected from different facilities, the performances of all the models are declined to some extents. But our proposed approach still outperforms the other models for most testing cases.

To illustrate the reconstruction performances more clearly, we also show the reconstruction results for testing images in Figure 3. We can see that our network can reconstruct the CT image with higher quality. Due to the space limit, the experimental results of different views $N_{v}$ and more visualized results are placed in our supplementary material.