Adversarial optimization for joint registration and segmentation in prostate CT radiotherapy

TL;DR

This paper introduces an adversarial deep learning framework for joint registration and segmentation in prostate CT radiotherapy, enabling real-time contour propagation and improved accuracy over traditional methods.

Contribution

It presents a novel unsupervised 3D deep learning model that estimates deformation fields for image registration without requiring segmentation labels during testing.

Findings

Outperforms conventional elastix registration in accuracy.

Reduces computation time significantly, enabling real-time application.

Effective in propagating contours in prostate CT scans for radiotherapy.

Abstract

Joint image registration and segmentation has long been an active area of research in medical imaging. Here, we reformulate this problem in a deep learning setting using adversarial learning. We consider the case in which fixed and moving images as well as their segmentations are available for training, while segmentations are not available during testing; a common scenario in radiotherapy. The proposed framework consists of a 3D end-to-end generator network that estimates the deformation vector field (DVF) between fixed and moving images in an unsupervised fashion and applies this DVF to the moving image and its segmentation. A discriminator network is trained to evaluate how well the moving image and segmentation align with the fixed image and segmentation. The proposed network was trained and evaluated on follow-up prostate CT scans for image-guided radiotherapy, where the planning CT contours are propagated to the daily CT images using the estimated DVF. A quantitative comparison with conventional registration using \texttt{elastix} showed that the proposed method improved performance and substantially reduced computation time, thus enabling real-time contour propagation necessary for online-adaptive radiotherapy.

Tables

Table 1. Table 1: MSD (mm) values for different experiments, where † † \dagger and ‡ ‡ \ddagger represent a significant difference compared to elastix -MI and Reg-CNN, respectively.

	Prostate	Seminal vesicles	Lymph nodes	Rectum	Bladder
Evaluation	$\mu \pm \sigma$	$\mu \pm \sigma$	$\mu \pm \sigma$	$\mu \pm \sigma$	$\mu \pm \sigma$
elastix-NCC	$1.81\pm 0.7$	$2.80\pm 1.6$	$1.19\pm 0.4$	$3.79\pm 1.2$	$5.31\pm 2.6$
elastix-MI	$1.73\pm 0.7$	$2.70\pm 1.6$	$1.18\pm 0.4$	$3.68\pm 1.2$	$5.26\pm 2.6$
Reg-CNN	$1.44\pm {0.5}^{†}$	$2.09\pm {1.7}^{†}$	$1.22\pm 0.3$	$2.59\pm {1.3}^{†}$	$4.18\pm {2.6}^{†}$
JRS-CNN	$1.18\pm {0.4}^{†‡}$	$1.91\pm {1.6}^{†‡}$	$1.02\pm {0.3}^{†‡}$	$2.32\pm {1.3}^{†‡}$	$2.37\pm {2.0}^{†‡}$
Reg-GAN	$1.40\pm {0.5}^{†}$	$2.14\pm {1.7}^{†}$	$1.06\pm {0.3}^{†‡}$	$2.72\pm {1.3}^{†}$	$4.31\pm {2.8}^{†}$
JRS-GANa	$1.13\pm {0.4}^{†‡}$	$1.81\pm {1.6}^{†‡}$	$1.00\pm {0.3}^{†‡}$	$2.21\pm {1.3}^{†‡}$	$2.29\pm {2.0}^{†‡}$
JRS-GANb	$1.17\pm {0.4}^{†‡}$	$1.90\pm {1.5}^{†‡}$	$1.01\pm {0.3}^{†‡}$	$2.34\pm {1.3}^{†‡}$	$2.41\pm {2.1}^{†‡}$

Table 2. Table 2: %95HD (mm) values for different experiments, where † † \dagger and ‡ ‡ \ddagger represent a significant difference compared to elastix -MI and Reg-CNN, respectively.

	Prostate	Seminal vesicles	Lymph nodes	Rectum	Bladder
Evaluation	$\mu \pm \sigma$	$\mu \pm \sigma$	$\mu \pm \sigma$	$\mu \pm \sigma$	$\mu \pm \sigma$
elastix-NCC	$4.2\pm 1.8$	$6.1\pm 3.3$	$2.8\pm {1.0}^{‡}$	$11.0\pm 5.2$	$15.4\pm {8.4}^{‡}$
elastix-MI	$4.0\pm 1.7$	$6.0\pm 3.7$	$2.8\pm {1.0}^{‡}$	$10.9\pm 5.2$	$15.3\pm {8.3}^{‡}$
Reg-CNN	$5.3\pm 2.5$	$6.2\pm 3.5$	$4.4\pm 1.4$	$11.0\pm 6.5$	$16.6\pm 9.3$
JRS-CNN	$3.6\pm {1.5}^{†‡}$	$5.4\pm {3.4}^{†‡}$	$3.1\pm 0.9‡$	$10.3\pm {6.7}^{†‡}$	$11.6\pm {10.5}^{†‡}$
Reg-GAN	$4.3\pm 2.1‡$	$6.0\pm 3.6$	$3.4\pm 1.0‡$	$11.1\pm 6.4$	$16.2\pm 9.6‡$
JRS-GANa	$3.4\pm {1.4}^{†‡}$	$5.3\pm {3.3}^{†‡}$	$3.1\pm 0.9‡$	$10.0\pm {6.7}^{†‡}$	$11.0\pm {10.3}^{†‡}$
JRS-GANb	$3.5\pm {1.4}^{†‡}$	$5.6\pm 3.7‡$	$3.0\pm 1.0‡$	$10.5\pm {6.8}^{†‡}$	$11.4\pm {10.6}^{†‡}$

Equations

$$L_G^{GAN}$$

$$L_D^{GAN}$$

$$L_{sim}=(1-\mathrm{DSC}(S_m(\Phi (\mathbf{x})),S_f(\mathbf{x})))+(1-\mathrm{NCC}(I_m(\Phi (\mathbf{x})),I_f(\mathbf{x}))),$$

$$L_G=L_{sim}+{\lambda }_1{L}_{smooth}+{\lambda }_2{L}_G^{GAN},$$

Full text

1 Introduction

Joint image registration and segmentation (JRS) has long been an active area of research in medical imaging. Image registration and segmentation are closely related and complimentary in applications such as contour propagation, disease monitoring, and data fusion from different modalities. Image registration could be enhanced and improved using an accurate segmentation, and vice versa registration algorithms could be used to improve image segmentation.

An important application in which coupling of image registration and segmentation is crucial, is online adaptive image-guided radiotherapy. In this application, clinically approved contours are propagated from an initial planning CT scan to daily inter-fraction CT scans of the same patient. Image registration can be used to correct for anatomical variations in shape and position of the underlying organs, as well as to compensate for any misalignment in patient setup. Ideally, contours should be propagated quickly to allow immediate computation of a new dose distribution. With these propagated contours, margins can be smaller and treatment-related complications may be reduced. Thus, it is important that the daily contours are of high quality, are consistent with the planning contours, and are generated in near real-time.

In the last decade, researchers have been working on fusing image registration and segmentation. Lu et al. [1] proposed a Bayesian framework for modelling segmentation and registration such that these could alternatingly constrain each other. Yezzi et al. [2] proposed using active contours to register and segment images. Unal et al. [3], generalizing on [2], proposed to use partial differential equations without any shape prior. Most of these methods require long computation times and complex parameter tuning. Recently, the widespread adoption of deep learning techniques has led to remarkable achievements in the field of medical imaging [4]. Among these techniques are generative adversarial networks (GANs), which are defined by joint optimization of a generator and discriminator network [5]. GANs have boosted the performance of traditional networks for image segmentation [6] as well as registration [7]. Recently, Mahapatra et al. [8] proposed a GAN for joint registration and segmentation of 2D chest X-ray images. However, this method requires reference deformation vector fields (DVFs) for training. In practice, these are often unavailable and it may be more practical to perform unsupervised registration [9], i.e. training without reference DVFs.

In this paper, we introduce a fast unsupervised 3D GAN to jointly perform deformable image registration and segmentation. A generator network estimates the DVF between two images, while a discriminator network is trained simultaneously to evaluate the quality of the registration and the segmentation and propagate the feedback to the generator network. We consider the use-case in which fixed and moving images as well as their segmentations are available for training, which is a common scenario in radiation therapy. However, no segmentations are required for DVF estimation during testing. This paper has the following contributions. First, we propose an end-to-end 3D network architecture, which is trained in an adversarial manner for joint image registration and segmentation. Second, we propose a strategy to generate well-aligned pairs to train the discriminator network with. Third, we leverage PatchGAN as a local quality measure of image alignment. Fourth, the proposed network is much faster and more accurate than conventional registration methods.

We quantitatively evaluate the proposed method on a prostate CT database, which shows that the method compares favorably to elastix software [10].

2 Methods

Image registration is the transformation of a moving image $I_{m}$ to the coordinate system of a fixed image $I_{f}$. In this paper, we assume that all image pairs are affinely registered beforehand, and we focus on local non-linear deformations.

In conventional contour propagation algorithms, registration and segmentation are disjoint. First, the DVF $\Phi$ is estimated using image registration, and then $\Phi$ is used to warp the contours $S_{m}$ to the fixed coordinate space. Afterwards, during system evaluation, a similarity measure such as the Dice similarity coefficient (DSC) can be used to measure the quality of the propagated contours w.r.t. ground truth contours, but this information is not fed back to the registration algorithm. We call this an open loop system. In contrast, this paper proposes an end-to-end closed loop system to improve image registration based on feedback on the registration as well as the segmentation quality.

2.1 Adversarial Training

We propose to train a GAN containing two CNNs: a generator network that predicts the DVF $\Phi$ given $I_{f}$ and $I_{m}$, and a discriminator network that assesses the alignment of $I_{f}(\bm{x})$ and $I_{m}(\Phi(\bm{x}))$ as well as the overlap between $S_{f}(\bm{x})$ and $S_{m}(\Phi(\bm{x}))$. Hence, we assume that $S_{f}$ and $S_{m}$ are both available, but during training only. The GAN is trained using a Wasserstein objective [11], which has empirically been shown to improve training stability and convergence compared to the GAN objective in [5]. Equations (1) and (2) list the generator loss $L_{G}^{GAN}$ and the discriminator loss $L_{D}^{GAN}$ of WGAN:

[TABLE]

where $G$ and $D$ denote the generator and discriminator networks with trainable parameters and $\Phi$ is the DVF provided by $G$. In a GAN, the discriminator is trained to distinguish between real and fake samples. In this case, fake samples are the triple $(I_{f},I_{m}(\Phi),S_{m}(\Phi))$, while real samples should be well-aligned images. As we perform unsupervised registration, and assume no knowledge about the ideal alignment of two images, we synthesize such image based on the fixed image and its segmentation alone: $(I_{f},\Theta(I_{f}),S_{f})$. Hence, $\Theta$ in Equation (2) is a random combination of disturbance functions, as follows. First, to mimic imaging noise, Gaussian noise and Gaussian smoothing are added with zero mean and a standard deviation of 0.04. Second, to mimic contrast variations, we apply gamma correction with a random gamma factor in the range $[-0.4,0.4]$. Third, we mimic interpolation errors by applying a random deformation of less than 0.5 mm and resample the images using that deformation using linear interpolation.

In addition to these image-based quality measures, we include the segmentation of the deformed moving image as input to the discriminator in order to enforce DVFs that are consistent with the moving segmentation. We test two designs. The first design concatenates the segmentation as a third input channel in the discriminator, next to the fixed and moving image channels. The second design multiplies the fixed and moving image channel with the corresponding segmentation, so that the network learns to focus on the target structures and organs-at-risk instead of on the bowels and other less relevant soft tissue. These designs are named JRS-GANa and JRS-GANb, respectively.

We found that training the network using WGAN loss only, resulted in slow convergence and suboptimal registrations. Thus, a similarity loss $L_{sim}$, based on image similarity and segmentation overlap, was added to the generator:

[TABLE]

where DSC is the Dice similarity coefficient and NCC is normalized cross-correlation. Adding the DSC to $L_{sim}$ ensures that the registration improves the segmentation and vice versa. Furthermore, to ensure smooth and continuous DVFs, the bending energy penalty of the DVF, $L_{smooth}$, was added as a regularization term to the overall generator loss, which was defined as:

[TABLE]

where $\lambda_{1}$ and $\lambda_{2}$ are weights for the DVF smoothness and the generator loss.

During training of the network, for every iteration of the generator we used 100 iterations of the discriminator, for the first 25 iterations. After that we used the ratio 1:5. In each iteration, weights of the discriminator were clipped to the range $[-0.01,0.01]$ [11].

2.2 Network Architectures

Generator Network To estimate the parametric mapping function $\Phi$ between the fixed and moving images we use a 3D network similar to the U-net [12]. Figure 1 shows the network design in more detail. The input to the network is the concatenation of $I_{f}$ and $I_{m}$. The network encodes the image pairs through a set of 3$\times$3$\times$3 convolution layers followed by LeakyReLU and batch normalization layers. Strided convolutions are used in the contractive path and upsampling layers are used in the expanding path. The output size of the network is smaller than the input size in order to consider a larger field of view. A resampling network adopted from NiftyNet [13] is used to warp the images using the estimated DVF during training time so that the network can be trained end-to-end.

Discriminator Network The discriminator is responsible for assessing whether the image pairs are well-aligned or not, as well as assessing whether the segmentations overlap. Figure 1 shows the network design, which is similar to the contracting path of the generator. The discriminator network was trained using PatchGAN [16]. Hence, instead of representing the quality of the whole patch with a single number, the network could quantify the sub-patch quality locally.

.

3 Experiments and Results

3.1 Dataset, evaluation criteria and implementation details

This study includes eighteen patients who underwent intensity-modulated radiation therapy for prostate cancer in 2007 at Haukeland university hospital [14]. Each patient had a planning CT as well as 7 to 10 inter-fraction repeat CT scans. The prostate, lymph nodes, seminal vesicles, as well as the rectum and bladder were annotated. Each scan has 90 to 180 slices with a slice thickness of around 2 to 3 mm. All the slices were of size 512 $\times$ 512 with an in-plane resolution of around 0.9 mm. All the volumes were affinely registered using elastix. The volumes were resampled to isotropic voxel size of 1$\times$1$\times$1 mm. All volumes intensities were scaled to [-1, 1]. We split the dataset into 111 image pairs (from 12 patients) for training and validation and 50 image pairs (6 patients) for testing.

The quality of registration is quantified geometrically in 3D by comparing the manual delineations of the daily CT with the automatically propagated contours. We use the mean surface distance (MSD), and the 95% Hausdorff distance (HD). A Wilcoxon signed rank test at $p=0.05$ is used to compare results.

The networks were implemented using TensorFlow (version 1.13) [15] with the RMSProp optimizer using a learning rate of $10^{-5}$. The networks were trained and tested on an NVIDIA Tesla V100 GPU with 16 GB of memory. From each image pair, 1000 patches of size 96$\times$96$\times$96 voxels were sampled within the torso mask. To improve stability, the network was trained to warp the fixed patch to the moving patch and vice versa at the same training iteration. The magnitude of the three loss terms in Equation (3) was scaled by setting $\lambda_{1}=1$ and $\lambda_{2}=0.01$.

3.2 Experiments and results

Tables 1 and 2 provide quantitative results comparing the following methods. First, we include conventional iterative methods using elastix software [10] with NCC (elastix-NCC) and MI (elastix-MI) similarity measures, using the settings from [17]. Second, we evaluate two unsupervised deep learning-based methods without adversarial feedback: One uses the generator trained with the NCC loss (Reg-CNN), similar to [9]; the other uses the generator with both the NCC and DSC loss (JRS-CNN). Third, we evaluate several versions of our GAN-based approach. To study the effect of adversarial training without added segmentations, we perform an experiment named Reg-GAN. Finally, we evaluate the proposed JRS-GANa and JRS-GANb methods.

The MSD values in Table 1 show that for all organs, the GAN-based methods significantly improved over elastix. This is further shown in Figure 2. The results indicate a significant improvement when performing joint registration and segmentation instead of disjoint registration. Furthermore, the boxplot indicates that performance for JRS-GANa and JRS-GANb was very similar. Similarly, the 95% HD values in Table 2 show improvements in contour accuracy when the GAN-based method is used. Especially the organs-at-risk showed large improvements. The standard deviations of the Jacobian determinant of the estimated DVFs were $0.08\pm 0.01$ and $0.17\pm 0.04$ for elastix-MI and JRS-GANa, respectively. The average runtime for the proposed pipeline is 0.6 seconds on the GPU for a volume of size 2563 voxels, while the average runtime of elastix at 100 iterations is 13 seconds per volume on an Intel Xeon E51620 CPU using 4 cores. Figure 3 illustrates the segmentation and registration for an example case.

4 Discussion and Conclusion

In this study, we investigated the performance of an end-to-end joint registration and segmentation network for adaptive image-guided radiotherapy. Unlike conventional registration methods, our network encodes and learns the most relevant features for joint image registration and segmentation, and exploits the combined knowledge on unseen images without segmentations.

We demonstrate that including the segmentation during training boosts the system’s performance by a margin. Furthermore, adversarial feedback had a small benefit on performance, when comparing Reg-CNN with Reg-GAN. Results indicate a noticeable benefit of including segmentation masks as input to the discriminator during training. How exactly segmentation masks were embedded during training was less relevant, with only small differences observed for the seminal vesicles. This could be due to the small size and irregular nature of the seminal vesicles. A key advantage of the proposed deep learning-based contour propagation method is its runtime on new and unseen data, i.e. 0.6 s.

This work has shown that adversarial feedback can help improve registration, i.e. that a discriminator can learn a measure of image alignment. This is a promising aspect that could be further explored in future work. This will include improved GAN objectives, such as the use of gradient penalty regularization.

To conclude, we have proposed a 3D adversarial network for joint image registration and segmentation with a focus on prostate CT radiotherapy. The proposed method demonstrated the effectiveness of training the registration and segmentation jointly. Moreover, it showed a substantial reduction in the computation time making it a strong candidate for online adaptive image-guided radiotherapy of prostate cancer. Since the proposed method did not only improve accuracy for the target areas, but substantially so for the organs-at-risk, this may aid reducing treatment-induced complications.

Acknowledgements.

This study was financially supported by Varian Medical Systems and ZonMw, the Netherlands Organization for Health Research and Development, grant number 104003012. The dataset with contours were collected at Haukeland University Hospital, Bergen, Norway and were provided to us by responsible oncologist Svein Inge Helle and physicist Liv Bolstad Hysing; they are gratefully acknowledged.