Synthesizing MR Image Contrast Enhancement Using 3D High-resolution ConvNets

Image-to-image (I2I) translation has been explored in recent years with the aim of translating an input image in a source domain to an image in a target domain. The basic idea of I2I methods is to learn a non-linear feature mapping given the input and output image pairs as training data. A large number of network architectures has been proposed to act as the non-linear mapping. For example, Long et al. [22] proposed to utilize a FCN model for image-to-segmentation. In Ronneberger et al. [8], a U-Net architecture was proposed for biomedical image segmentation, which is currently widely used in medical image translation tasks. In Chen et al. [23], the authors introduced the dilated convolution to enlarge the receptive field of neural networks. In Zhao et al. [24], a pyramid pooling module was proposed to fuse features under four different pyramid scales, which enable the model to utilize local and global context information for pixel-wise prediction. In Newell et al. [25], the stacked hourglass module was proposed to capture and consolidate information across all scales of the image for human pose estimation. Different from the traditional high-to-low and low-to-high FCN architectures, Sun et al. [26] proposed a high-resolution net, which maintains high-resolution feature maps throughout processing.

A large number of training objectives have been introduced to measure the difference between the generated image and the ground truth image in I2I translation tasks. The typical choices are $\ell_{1}$ and $\ell_{2}$ loss. In [27], a differentiable variant of $\ell_{1}$ loss, named the Charbonnier penalty function, was proposed to handle outliers. Compared with $\ell_{1}$ or $\ell_{2}$ loss, which may lead to blurry images [28], adversarial loss [28, 29, 20, 18] has become a popular choice for I2I tasks. This process trains a discriminator to distinguish generated images from ground truth images. Additionally, perceptual loss, which measures the difference in feature space, has also been widely used in I2I translation [30, 31].

Recently, an increasing number of machine learning and deep learning methods have shown potential in medical image synthesis, which estimates a desired imaging modality from other modalities or scans. For example, Li et al. applied the 3D-CNN model to predict missing PET patterns from MRI data [32]. To improve the quality of 3T MR images, Bahrami et al. collected a dataset with paired 3T and 7T images scanned from the same subjects, and proposed to reconstruct 7T-like images from 3T images [17, 33]. In Xiang et. al [34], the authors proposed to estimate standard-dose positron emission tomography (PET) images from low-dose PET and MRI images. In Huang et al. [35], the authors proposed a weakly-supervised convolutional sparse coding method to simultaneously solve the problems of super-resolution and cross-modality image synthesis. In Dar et al. [20], the authors proposed a new approach for multi-contrast MRI synthesis based on conditional generative adversarial networks, employing adversarial loss to preserve intermediate-to-high frequency details. In Han et al. [29], a GAN model is employed to synthesize rich and diverse brain MR images from existing MR images. In Nie et al. [36, 18], a 3D FCN model was trained to transform MRI to CT images using an adversarial strategy to train the FCN, which enforces the generated images to be more realistic. Finally, additional works have investigated methods for MRI to CT image synthesis [37, 38, 39], multimodal MRI synthesis [40, 41], and high-quality PET synthesis [42].

The studies that most relevant to ours can be seen in [6, 21, 43]. Specifically, in [6], the author utilized a U-Net-like model to synthesize full-dose CE-MRI from zero-dose pre-contrast MRI and the 10% low-dose postcontrast MRI. In [21], a 3D U-Net model was developed to generate CE-MRI, which utilizes 10 multiparametric MRI sequences acquired prior to GBCA application as inputs. As a result, their model demonstrates a peak signal to noise ratio (PSNR) of 22.967dB and a structural similarity index (SSIM) of 0.872 dB for the whole brain region. In [43], the author utilized the residual attention U-Net architecture to estimate CE-MRI from non-contrast T2 MRI for cerebral blood volume (CBV) mapping in mice brain.

Our dataset was acquired at UCLA on Siemens 3 Tesla MRI systems as part of standard-of-care for brain tumor patients. The protocol used was consistent with the International Standard Brain Tumor Protocol [44] and includes 3D MPRAGE T1-weighted pre- and post-contrast imaging, axial 2D T2-weighted imaging, and axial 2D diffusion-weighted imaging used in the calculation of the ADC map. A total of 426 scans from 300 brain tumor patients were included in this study. The data includes two parts, A and B. Set A consisted of 411 scans. It was used for training purposes and therefore further subdivided in 369 scans for training, and 42 scans randomly selected for validation. Set B contained 15 test samples with precise tumor masks and were used to evaluate the quantitative performance on tumor regions. Note that scans in Set B are patients from the UCLA brain tumor trial (IRB# 14-001261), and were selected randomly from the available data making sure both enhancing and non-enhancing tumor were part of the cohort. The experts with more than 10 years of experience created the tumor ROIs as part of the clinical trial reads.

Pre-contrast T1, T2 and ADC map, were utilized as input images and the contrast-enhanced T1 is utilized as the target image. Note that apparent diffusion coefficient (ADC) maps was chosen to augment T1 and T2 because it is independent of these image contrasts and may provide additional information for CE-T1 image synthesis. ADC maps were derived from standard, isotropic diffusion weighted images (DWIs) with and without diffusion weighting according to the standardized brain tumor imaging protocol (BTIP) [45]. Simply, ADC was calculated from b=0, 500, and 1000 s/mm2 by fitting the equation ADC=-1/b*ln(S(b)/S0), where b = 500 and 1,000 s/mm2, ln is the natural log, S(b) is the signal intensity for an MR image at the given b-value, and S0 is the signal intensity of the MR image without any diffusion weighting (b=0).

All the sequences of the MRI data were co-registered to match the targeted 3D contrast-enhanced T1. Bilinear interpolation was utilized to resize all the MRI data to the volume size of 192×256×192 voxels. To remove the skull, brain masks were created for the 3D contrast-enhanced T1 sequences utilizing FSL’s brain extraction tool [46]. Besides, in order to remove the side effects of the background slices, in the training and evaluation stage, we only selected the foreground slices and remove the top and bottom background slices which are less informative. Finally, all MRI scans were pre-processed by image equalization and the intensity values of the voxels within the brain region were normalized to [0,1].

Let $\mathbf{X}=[\mathbf{X}_{T1};\mathbf{X}_{T2};\mathbf{X}_{ADC}]\in\mathbb{R}^{h\times w\times d\times 3}$ denote the MRI sequences with three MRI scans, and $\mathbf{Y}\in\mathbb{R}^{h\times w\times d\times 1}$ denote the contrast-enhanced T1 sequence. To synthesize the CE-T1 sequences from the non-enhanced MRI scans, a 3D FCN model is designed to work as a non-linear mapping function $f_{\bm{\theta}}$, such that $f_{\bm{\theta}}:\mathbf{X}\rightarrow\mathbf{Y}$, where $\bm{\theta}$ is the model parameter to be learned. As shown in Fig. 2, the introduced 3D FCN generator is composed of stacked 3D convolutional layers, batch normalization layers, and non-linear activation layers. The first stage is a stem network [26] which is composed of three individual FCN branches that used to handle different input modalities. The resulting feature maps corresponding to different modalities are then fused by concatenation. Following the stem net are repeated multi-resolution subnetworks with generation multi-scale fusion stages. We start from a high-resolution subnetwork in the second stage and add high-to-low resolution subnetworks one-by-one gradually. Specifically, in the second stage, a $3\times 3\times 3$ convolution with stride 2 is used to obtain the $2\times$ downsampling feature maps. In the third stage, $4\times$ downsampling feature maps are obtained from the higher resolution feature maps. As a result, we have three different resolution feature maps for different subnetworks, which corresponds to multiple-scale information. To fuse the multi-resolution information comprehensively, multi-scale fusion stages are introduced, which ensure multi-resolution information exchange across different parallel subnetworks. Specifically, in the 3rd and 4th stages, different subnetworks aggregate the feature maps from the other parallel subnetworks. During the multi-scale fusion stages, upsample and downsample operators are utilized to match the size of the feature maps in different subnetworks. In the last stage, feature maps from different branches are fused by concatenation, and a ResBlock is utilized to obtain the final CE-T1.

As illustrated in Fig. 2, Let $F_{mn}$ denotes the feature maps generated by the $m$-th stage and the $n$-th subnetwork, then the feature maps can be calculated as: $F_{11}=\phi(T1)\oplus\phi(T2)\oplus\phi(ADC)$, $F_{12}=D_{2}(F_{11})$, $F_{21}=\phi(F_{11})$, $F_{22}=\phi(F_{12})$, $F_{23}=D_{2}(F_{22})$, $F_{31}=\phi(F_{21}\oplus U_{2}(F_{22}))$, $F_{32}=\phi(F_{22}\oplus D_{2}(F_{21}))$, $F_{33}=\phi(F_{23}\oplus D_{4}(F_{21}))$, $F_{41}=\phi(F_{31}\oplus U_{2}(F_{32})\oplus U_{4}(F_{33}))$, $F_{42}=\phi(D_{2}(F_{31})\oplus F_{32}\oplus U_{2}(F_{33}))$, $F_{43}=\phi(D_{4}(F_{31})\oplus D_{2}(F_{32})\oplus F_{33})$ and $F_{5}=F_{41}\oplus U_{2}(F_{42})\oplus U_{4}(F_{43})$. Here, $\phi(x)$ denotes the feature mapping function determined by the FCN module parameters, $\oplus$ denotes feature map concatenation along channel dimension, $D_{2}$ denotes 2x downsampling, $D_{4}$ denotes 4x downsampling, $U_{2}$ denotes 2x upsampling and $U_{4}$ denotes 4x upsampling. The FCN module consists of four ResBlocks and each ResBlock is composed of two ”Conv-BN-Relu” layers. The widths (number of feature maps) of the three parallel subnetworks are 64, 128, and 256. Detailed information regarding the model architecture can be seen in our source code, which is available at https://github.com/chenchao666/Contrast-enhanced-MRI-Synthesis.

The advantages of the presented model are four fold. First, we utilize three MRI scans as inputs and employ three individual stem nets for different modalities, which is able to preserve the modality-specific information. Second, the presented model maintains a high resolution representation throughout the processing pipeline, and contains three parallel subnetworks that can generate and process multi-scale information in parallel. Third, the repeated multi-scale fusion stages ensure better feature fusion across different scale. Fourth, 3D convolution was utilized to exploit additional information from neighboring slices.

Let $\mathbf{X}$ denote the input non-enhanced MRI sequences and $\mathbf{Y}$ denote the contrast-enhanced T1. Our goal is to learn a mapping function $f$ which can generate the CE-T1 sequences $\hat{\mathbf{Y}}=f_{\theta}(\mathbf{X})$, such that the synthetic CE-T1 is close to the ground truth $\mathbf{Y}$. The loss function utilized to train the model consists of three terms: pixel-wise MAE loss $\mathcal{L}_{MAE}$, Structural Similarity loss (SSIM) $\mathcal{L}_{SSIM}$ and a local loss $\mathcal{L}_{local}$ to focus performance on tumor regions.

$\bullet$ Pixel-wise Loss: The MAE loss and MSE loss are the most widely used pixel-wise losses for image synthesis. We found that using MSE loss resulted in blurrier images compared to MAE loss in this task. Therefore, the pixel-wise MAE loss was utilized in our model, which is given as

$\bullet$ SSIM Loss: Using pixel-wise loss alone may ignore image structures. Therefore, we also utilize SSIM loss [47], which is perceptually motivated and leads to more realistic images. The SSIM loss is defined as

where $n$ denotes the number of slices of the output 3D MRI sequences, and $\mathbf{Y}_{i}$ denotes the $i$-th slice of the ground truth CE-T1. $SSIM(\mathbf{x},\mathbf{y})$ outputs a scalar between 0 and 1, which indicates the structural similarity between images $\mathbf{x}$ and $\mathbf{y}$. The definition of the SSIM metric can be seen in 4.1.

$\bullet$ Local Loss: Tumor regions are of particular interest, but account for a very small proportion of voxels in the entire MRI seqeuences. This data imbalance problem leads to under-fitting and poor performance for the tumor regions. Therefore, we introduce a local loss to increase the contribution of the tumor regions. The local loss is defined as

where $\mathbf{M}$ is a binary mask of the tumor regions, and $\odot$ denotes voxel-wise multiplication. Since it is very expensive to assign voxel-level labels for each slice, we do not have precise tumor masks for the training samples. Fortunately, compared to the non-enhanced T1 images, the tumors and vessels are significantly enhanced in the CE-T1 images due to the utilization of GBCAs. Therefore, we can calculate a rough tumor mask by thresholding the difference of the T1 and CE-T1 images as follows,

where $\delta$ is the threshold to control the size of the mask, $\mathbf{X}_{\text{T1}}$ is the non-enhanced T1w image, and $\mathbf{Y}$ is the contrast-enhanced T1 image. Fig. 3 shows several examples of the generated binary masks. The generated masks are able to identify regions that are highly enhanced, such as tumors and vessels. These regions are used with our local loss to correct for data imbalance.

$\bullet$ Overall Loss: The overall loss function is defined as the weighted sum of all the three terms,

where $\lambda_{1}$, $\lambda_{2}$ and $\lambda_{3}$ are trade-off parameters to balance the contribution of each loss term.

The proposed network was implemented in Python with the Keras library and trained on an NVIDIA DGX system with eight NVIDIA V100 CPUs and 512G memory. The Adam optimizer was utilized for training. Since feeding the whole 3D MRI sequences into the model leads to out-of-memory (OOM) problem, we follow [18, 39] to adopt a patch-based strategy for training, and set the batch size to three. In each iteration step, three-slices of MRI sequence with size $3\times 256\times 192$ are randomly sampled from each sequence volume. Therefore, we have three input channels with size $3\times 3\times 256\times 192$ and one output with the size $3\times 256\times 192$. Model training is divided into two stages. In the first stage, we set $\lambda_{1}=1.0$, $\lambda_{2}=1.0$, and $\lambda_{3}=1.0$, and train the model for the first 40 epochs with a learning rate $lr=0.0001$. In the second stage, we alter the trade-off parameter of the local loss by setting $\lambda_{1}=0.1$, $\lambda_{2}=0.1$, $\lambda_{3}=10$ and fine-tune the model for another 10 epochs with the learning rate $lr=0.00001$. We empirically the threshold utilized to obtain the brain mask to $\delta=0.1$ throughout the experiments.

$\bullet$ Baseline Model: To evaluate the effectiveness of our proposed model, we implemented a 2D U-Net, a 3D U-Net, as well as the 2D version of our proposed network. All the baseline models are evaluated on the same training and test sets. For the 2D U-Net, we expanded the model depth to 2x of the original U-Net network. For the 3D U-Net, we expanded the model depth to 1.5x of the original model depth.

$\bullet$ Evaluation Metric: The quantitative performance of our model and the baseline models was measured using mean absolute error (MAE), peak signal-to-noise ratio (PSNR) and structural similarity index (SSIM). Given two images $\mathbf{y}\in\mathbf{Y}$ and $\hat{\mathbf{y}}\in\hat{\mathbf{Y}}$, where $\mathbf{y}$ is the ground truth image and $\hat{\mathbf{y}}$ is the predicted image. Then,

• •

  $MAE(\mathbf{Y},\hat{\mathbf{Y}})=\frac{1}{\Omega_{\mathbf{Y}}}\|\mathbf{Y}-\hat{\mathbf{Y}}\|_{1}$

• •

  $MSE(\mathbf{Y},\hat{\mathbf{Y}})=\frac{1}{\Omega_{\mathbf{Y}}}\|\mathbf{Y}-\hat{\mathbf{Y}}\|_{2}^{2}$

• •

  $PSNR(\mathbf{Y},\hat{\mathbf{Y}})=10\log_{10}(\tfrac{MAX_{I}^{2}}{MSE})$

• •

  $SSIM(\mathbf{Y},\hat{\mathbf{Y}})=\frac{1}{\Phi_{\mathbf{Y}}}\sum\limits_{\mathbf{y},\hat{\mathbf{y}}}\frac{(2\mu_{\mathbf{y}}\mu_{\hat{\mathbf{y}}}+C1)(2\sigma_{\mathbf{y}\hat{\mathbf{y}}}+C2)}{(\mu_{\mathbf{y}}^{2}+\mu_{\hat{\mathbf{y}}}^{2}+C1)(\sigma_{\mathbf{y}}^{2}+\sigma_{\hat{\mathbf{y}}}^{2}+C2)}$

where $\Omega_{\mathbf{Y}}$ is the number of voxels in $\mathbf{Y}$, $\Phi_{\mathbf{Y}}$ is the number of slices in $\mathbf{Y}$, $\mu_{\mathbf{y}}$ and $\sigma_{\mathbf{y}}$ are the mean and variance value of image $\mathbf{y}$, and $\sigma_{\mathbf{y}\hat{\mathbf{y}}}$ is the covariance between the ground truth image $\mathbf{y}$ and the predicted image $\hat{\mathbf{y}}$. $MAX_{I}$ is the maximum value of the image $\mathbf{y}$. Theoretically, lower MAE values and higher PSNR and SSIM values indicate better image generation quality. The statistical significance of experimental results was evaluated using paired t-tests.

Qualitative and Visual Assessment The results of representative test samples are shown in Fig. 4. The first two rows demonstrate a normal subject without any tumor and the remaining six rows demonstrate subjects with tumors. We compared the visual performance between our model (both 2D and 3D model) and the 2D/3D U-Net, which are widely used techniques for medical image synthesis tasks [6, 21]. As can be seen, the results reveal several interesting observations. First, both U-Net models and our proposed model generate promising visual performance for the normal patients. The vessels and high-frequency details in the generated images are very close to the ground truth images. Second, since the tumor voxels are limited relative to the entire MRI volume, models are more likely to over-fit to normal regions. As a result, the performance of the tumor regions is much worse than the performance on the normal regions. Third, compared with the U-Net model, which misses or under-estimates most of the tumors for the abnormal samples, our model often yields better visual performance for the tumor regions. Fourth, compared with the 2D model, the 3D model achieves better performance, especially for tumor regions. We believe this is because the 3D model takes advantage of the information from the neighboring slices. Fifth, while achieving promising visual performance for the whole brain MRI, our model sometimes misses or under-estimates some tumors, especially for those that are not distinct enough in the non-contrast images. Finally, Fig. 5 shows two representative MRI sequences in test set B. As observed, our model synthesized satisfactory contrast enhanced T1 images for different slices of the given 3D volume.

Quantitative Evaluation As shown in Table 1, we show the quantitative performance of our proposal and comparison methods in test set A. All metrics are computed based on the brain region. Our model significantly outperforms the U-Net model in PSNR and SSIM metrics and the 3D FCN model outperforms its 2D counterpart, which is consistent with our visual assessment. Specifically, U-Net(3D) outperforms U-Net(2D) by 0.35 points in PSRN (p-value=0.021), and Ours(3D) outperforms Our(2D) by 0.37 points (p-value=0.018), which demonstrates the effectiveness of utilizing 3D spatial information. Besides, the proposed Ours(3D) significantly outperforms the U-Net(3D) by more than one point (P-value=0.00037), we believe this is because our proposed model maintains a high resolution representation throughout the processing pipeline and contains three parallel subnetworks with multi-scale fusion stages. Compared with an existing method [21] that performs virtual contrast enhancement with deep learning, our proposal has more than five points improvement in PSNR, and has five points improvement in SSIM. Note that in [21], the author utilized 10-channel multiparametric MRI data as input while we only utilize three, which is a subset of their data. Our model outperforms [21] by a large margin even with less input data, which demonstrates the superiority of our proposed framework.

In order to evaluate the quantitative performance on the tumor region, we also collected 15 test patients with precise tumor masks in set B. The test performance in set B is shown in Table 2, with performance on both brain region and tumor region presented. The overall performance on the brain region is similar to the results on set A. Our proposed method significantly outperforms the U-Net and existing work [21]. For results on tumors, the proposed model outperforms the U-Net(3D) model by a large margin (p=0.0036), and the best performance on tumors is 21.2 in PSNR. Note that the quantitative performance on tumors is far from perfect and much worse than the performance on the whole brain region, this is because the tumor pixels are out-of-distribution and the model therefore tends to underestimate on tumor regions. It is worth noting that [21] utilizes a U-Net shape model to segment the tumor masks, which they utilize to evaluate the performance on tumors. These masks were not reviewed by a radiologist and therefore their quantitative performance on tumors may be inaccurate due to segmentation errors.

Comparison with State-of-the-Arts In addition to the U-Net structure that are used for MRI virtual contrast enhancement[6, 21], we also compare our proposal with several state-of-the-art medical image synthesis methods, including Pix2Pix[28], DECNN[37], LA-GANs[42] and MedGAN[48]. Quantitative comparison between our proposal and other state-of-the-art methods is presented in Table 4. Pipx2Pix produces the worst result, while DECNN and LA-GANs achieve similar performance and both outperform Pix2Pix. The MedGAN outperforms the previous methods by a large margin, we believe this is because MedGAN utilizes a cascade U-Blocks as generator which is deeper and better designed than the generator in [28, 37]. Compared with the MedGAN, our proposal shows impressive improvement (p-value=0.0064) on brain region, and also significantly outperforms the MedGAN on tumor region by more than 1.3 points (p-value$<$1e-4), which demonstrate the superiority of the introduced framework. Note that the performance improvement on the tumor region is much significant than the improvement on the whole brain region, this is because the comparison methods does not take into the performance on tumor regions, while we introduce a local loss to improve the performance on tumors. The significant improvement on tumors also demonstrate the effectiveness of the local loss.

Impact of the Model Architecture To demonstrate the effectiveness of the proposed framework, we compare the performance of the final model with four degraded model architectures, which are: (1) Model-A: utilize one FCN module to handle different input modalities, i.e., the three modalities (T1, T2 and ADC) are fused in the input layer; (2) Model-B: remove the 2x downsampling subnetwork and the 4x downsampling subnetwork, only utilizing the high-resolution branch as the FCN generator; (3) Model-C: remove the 1x high resolution subnetwork; (4) Model-D: remove the repeated multi-scale fusion in stage 3 and stage 4. Table 5 shows the qualitative comparison between different model architectures. As can be seen, our proposed model yields notable improvement over the comparison degraded model architectures. In particular, the final model outperforms the Model-A by more than 0.5 points in PSNR, which demonstrates that employ three individual FCN modules for different input modalities is able to preserve the modality-specific information and lead to better performance. The final model outperforms the Model-B by more than 1 point, this is because it is difficult for the model to capture the global information only using the high resolution branch. Besides, The final model also outperforms the Model-C by more than 0.8 points on average, which indicates the effectiveness of maintaining the high resolution representation. Furthermore, the final model also outperforms the degraded Model-D, which proves the effectiveness of using repeated multi-scale fusion stages.

Impact of Different Input Modalities To investigate the influence of different input modalities, we train the model with different sequences. Specifically, in addition to using all three modalities (T1, T2 and ADC) as inputs, we also train the model: (1) using only T1 as input, and (2) using both T1 and T2 modalities as inputs. Experimental results are presented in Table 3. We can conclude from the results that: (1) using T1 alone as input obtains a satisfactory performance; (2) the incorporated T2 and ADC modalities introduced additional information, which further improves performance; and (3) compared to ADC, T2 is more informative and improves performance significantly. Besides, we also present the visual performance of training the model with different input modalities in Fig. 6.

Sensitivity of the Local Loss In order to investigate the influence of the introduced local loss, we used five groups of representative trade-off parameters to train the model. The parameters $(\lambda_{1},\lambda_{3})$ were selected from $\{(1.0,1.0),(1.0,10),(1.0,30),(0.1,10),(0.01,10)\}$. Note that we can not only increase the $\lambda_{3}$ to improve the influence of the tumor regions as it will cause the gradient to be too large and the network will not converge. Therefore, we promote the influence of the local loss by increasing a ratio $r=\frac{\lambda_{3}}{\lambda_{1}}$. $\lambda_{2}$ is set to the same value as $\lambda_{1}$ throughout the experiments. Results are presented in Fig. 7 and indicate that by increasing the ratio $r=\frac{\lambda_{3}}{\lambda_{1}}$, the PSNR performance of the brain region decreased and the PSNR performance of the tumor region increased. We present two representative samples in Fig. 8 to demonstrate the influence of local loss visually. By increasing the ratio $\frac{\lambda_{3}}{\lambda_{1}}$, tumor enhancement becomes more salient while the overall image becomes more blurry, which is consistent with our quantitative results. The experimental results show that one suitable ratio can be set to $r=\frac{\lambda_{3}}{\lambda_{1}}=100$, which leads to better enhancement results for tumors and satisfactory performance for the whole brain.

Sensitivity of the threshold $\bm{\delta}$ To demonstrate the robustness of the thresholded mask, we also perform a parameter sensitivity analysis experiment on the thresholds $\delta$. Fig. 9 shows the variation of PSNR performance on brain and tumor regions when the threshold is changed from $\delta\in\{0.001,0.01,0.03,0.05,0.1,0.15,0.2,0.3\}$. As can be seen, the PSNR performance on the brain region increases as the threshold increases and gradually stabilize when $\delta>0.1$. Besides, the PSNR performance on the tumor region increases first and then decrease as $\delta$ increases, and shows a bell-shaped curve. We believe this is because when the threshold is set to $\delta<0.03$, the generated mask $\mathbf{M}$ will include more non-tumor areas, and when $\delta>0.15$, the generated mask $\mathbf{M}$ will underestimate the tumor regions. Therefore, based on the performance on brain and tumor region, the best threshold can be set to $\delta\in[0.03,0.15]$.

The proposed model works well when all three input modalities are available. However, rather than having complete three modalities, it is common to have missing modalities in clinical scenarios. To understand how the model performs when only a subset of modalities are available, we visualize the inference performance when: (b) T1 is unavailable; (c) T2 is unavailable; (d) ADC is unavailable; (e) both T2 and ADC are unavailable in Fig. 10. Note that we utilize a zero matrix to represent the unavailable input modality during inference. The results shows that the model performance will be severely degraded when one or two modalities are not available. Besides, compared with the model performance when T2 or ADC is unavailable, the model performs much worse when T1 is unavailable. To ensure that the model is able to produce satisfactory results when only a subset of inputs are available, we also train a T1$\rightarrow$T2 and a T1$\rightarrow$ADC synthesis model utilizing a similar framework. In this way, when T2 or ADC is unavailable, we can first generate the T2 and ADC and then utilize the generated data to synthesize the contrast-enhanced T1 image. Fig. 11 shows the visual performance of our trained T1$\rightarrow$T2 and T1$\rightarrow$ADC synthesis model, which demonstrates that our proposed framework also generates promising results for cross-modal image synthesis task. Fig. 10(f) illustrates the model prediction of using T1, synthetic T2 and synthetic ADC as input, which shows that using the synthetic input data, our model can generate similar results as using full three inputs.