Keywords

1 Introduction

Radiologists often diagnose the progress of disease by comparing medical images at different temporal phases. In case of diagnosis of liver tumor such as hepatocellular carcinoma (HCC), the contrast of normal liver tissue and tumor region in contrast enhanced CT (CECT) distinctly varies before and after the infection of contrast agent. This provides radiologists an important clue to diagnose cancers and plan surgery or radiation therapy [7]. However, liver images taken at different phases are usually different in their shape due to disease progress, breathing, patient motion, etc, so image registration is important to improve accuracy of dynamic studies.

Classical image registration methods [4, 10] are usually implemented in a variational framework that solves an energy minimization problem over the space of deformations. Since the diffeomorphic image registration ensures the preservation of topology and one-to-one mapping between the source and target images, the algorithmic extensions to large deformation such as LDDMM [3] and SyN [1] have been applied to various image registration studies. However, these approaches usually require substantial time and extensive computation.

To address this issue, recent image registration techniques are often based on deep neural networks that instantaneously generate deformation fields. In supervised learning approaches [11], the ground-truths of deformation fields are required for training neural networks, which are typically generated by the traditional registration method. However, the performance of these existing supervised methods depends on the quality of the ground-truth registration fields, or they do not explicitly enforce the consistency criterion to uniquely describe the correspondences between two images.

In order to overcome the aforementioned limitations and provide topology-preserving guarantee, many unsupervised learning methods have been developed in these days. Balakrishnan et al. [2] propose 3D medical image registration algorithm using a spatial transform network. Zhang [12] presents a CNN framework that enforces an inverse-consistent constraint for the deformation fields. However, for the registration of large deformable volumes such as livers, the existing unsupervised learning methods often result in inaccurate registration due to the potential for the degeneracy of the mapping. Although Dalca et al. [5] tried to address this problem by incorporating a diffeomorphic integration layer, we found that its application of the liver registration is still limited.

In this paper, we present a novel unsupervised registration method using convolutional neural networks (CNN) with cycle-consistency [13]. We show that the cyclic constraint can be adopted for the image registration case naturally, and this cycle consistency improves topology preservation in generating fewer folding problems. Also, our network is trained with diverse source and target images in multiphase CECT acquisition so that a single neural network of our method provides deformable registration between every pairs once the network is trained. Experimental results demonstrate that the proposed method performs accurate 3D image registration for any pair of images within a few seconds in the challenging problem of 3D liver registration in multiphase CECT.

2 Proposed Method

The overall framework of our method is illustrated in Fig. 1. For the input images, A and B, in different phases, we define two registration networks as \(G_{AB}:(A, B) \rightarrow \phi _{AB}\) and \(G_{BA}:(B, A) \rightarrow \phi _{BA}\), where \(\phi _{AB}\) (resp. \(\phi _{BA}\)) denotes the 3-D deformation field from A to B (resp. B to A). We use a 3D spatial transformation layer \(\mathcal {T}\) in the networks to warp the moving image by the estimated deformation field, so that the registration networks are trained to minimize the dissimilarity between the deformed moving source image and fixed target image. Accordingly, once a pair of images are given to the registration networks, the moving images are deformed into the fixed images.

Fig. 1.
figure 1

The overall framework of the proposed method for image registration. The input images in different phases are denoted as A and B, and their phase and shape are denoted as P and S, respectively. The short-dashed line indicates the floating image and the long-dashed line denotes the fixed image.

Full size image

To guarantee the topology preservation between the deformed and fixed images, we here adopt the cycle consistency constraint between the original moving image and its re-deformed image. That is, the deformed volumes are given as the inputs to the networks again by switching their order to impose the cycle consistency. This constraint ensures that the shape of deformed images successively returns to the original shape.

2.1 Loss Function

We train the networks by solving the following loss function:

$$\begin{aligned} \mathcal {L} = \mathcal {L}_{regist}^{AB} + \mathcal {L}_{regist}^{BA} + \alpha \mathcal {L}_{cycle} + \beta \mathcal {L}_{identity}, \end{aligned}$$
(1)

where \(\mathcal {L}_{regist}\), \(\mathcal {L}_{cycle}\), and \(\mathcal {L}_{identity}\) are registration loss, cycle loss, and identity loss, respectively (see Fig. 2), and \(\alpha \) and \(\beta \) are hyper-parameters. Based on this loss function, our method is trained in an unsupervised manner without ground-truth deformation fields.

Registration Loss. The registration loss function is based on the energy function of classical variational image registration. For example, the energy function for the registration of floating image A to the target volume B is composed of two terms:

$$\begin{aligned} \mathcal {L}_{regist}^{AB}= \mathcal {L}_{sim}\left( \mathcal {T}(A, \phi ), B\right) + \mathcal {L}_{reg}({\phi }), \end{aligned}$$
(2)

where A is the moving image, and B is the fixed image. \(\mathcal {L}_{sim}\) computes image dissimilarity between the deformed image by the estimated deformation field \(\phi \) and the fixed image, and \(\mathcal {L}_{reg}\) evaluates the smoothness of the deformation field. Here, \(\mathcal {T}\) denotes the 3D spatial transformation function. In particular, we employ the cross-correlation as the similarity function to deal with the contrast change during CECT exam, and the L2-loss for regularization function. Accordingly, our registration loss function can be written as:

$$\begin{aligned} \mathcal {L}_{regist}^{AB} = -(\mathcal {T}(A, \phi _{AB}) \otimes B ) +\lambda ||\phi _{AB}||_2, \end{aligned}$$
(3)

where \(\otimes \) denotes the cross-correlation defined by

$$\begin{aligned}&x\otimes y = \frac{ |\langle x-\bar{x}, y-\bar{y} \rangle |^2}{ \Vert x-\bar{x}\Vert \Vert y-\bar{y}\Vert }, \end{aligned}$$
(4)

where \(\bar{x}\) and \(\bar{y}\) denote the mean value of x and y, respectively.

Fig. 2.
figure 2

The diagram of loss function structure in our proposed method. The short- and long-dashed lines are for floating image and fixed image, respectively.

Full size image

Cycle Loss. The cycle consistency condition is implemented by minimizing the loss function as shown in Fig. 2(a). Since an image A is first deformed to an image \(\hat{B}\), which is then deformed again by the other network to generate image \(\tilde{A}\), the cyclic consistency imposes \(A\simeq \tilde{A}\). Similarly, an image B should be successively deformed by the two networks to generate image \(\tilde{B}\). Then, the cyclic consistency is to impose that \(B\simeq \tilde{B}\).

As shown in Fig. 2(a), since the network in our registration receives both the moving image and the fixed image, the implementation of the cycle consistency loss should be given by as the vector-form of the cycle consistency condition:

$$\begin{aligned} \begin{bmatrix} A\\ B\end{bmatrix} \simeq \begin{bmatrix} \mathcal {T}(\hat{B}, \hat{\phi }_{BA})\\ \mathcal {T}(\hat{A}, \hat{\phi }_{AB}) \end{bmatrix} = \begin{bmatrix} \mathcal {T}\left( \mathcal {T}(A, \phi _{AB}), \hat{\phi }_{BA}\right) \\ \mathcal {T}\left( \mathcal {T}(B, \phi _{BA}), \hat{\phi }_{AB}\right) \end{bmatrix} \end{aligned}$$
(5)

where \( (\hat{B}, \hat{A}):= \left( \mathcal {T}(A, \phi _{AB}), \mathcal {T}(B, \phi _{BA})\right) .\) Thus, the cycle loss is computed by:

$$\begin{aligned} \mathcal {L}_{cycle} = \left\| \begin{bmatrix} \mathcal {T}(\hat{B}, \hat{\phi }_{BA})\\ \mathcal {T}(\hat{A}, \hat{\phi }_{AB}) \end{bmatrix} -\begin{bmatrix} A\\ B\end{bmatrix} \right\| _1, \end{aligned}$$
(6)

where \(||\cdot ||_1\) denotes the \(l_1\)-norm.

Identity Loss. Another important consideration for the design of loss function is that the network should not change the stationary regions of the body, i.e. the stationary regions should be the fixed points of the network. As shown in Fig. 2(b), this constraint can be implemented by imposing that the input image should not be changed when the identical images are used as the floating and reference volumes. More specifically, we use the following identity loss:

$$\begin{aligned} \mathcal {L}_{identity} = -(\mathcal {T}(A, G_{AB}(A,A)) \otimes A) - (\mathcal {T}(B, G_{BA}(B,B)) \otimes B ). \end{aligned}$$
(7)

By minimizing this identity loss (7), the cross-correlation between the deformed image and the fixed image can be maximized. Thus, the identity loss guides the stability of the deformable field estimation in stationary regions.

2.2 Network Architecture and 3D Spatial Transformation Layer

To generate a displacement vector field in width-, height-, and depth direction, we adopt VoxelMorph-1 [2] as our baseline network. Note that our model without both the cycle and identity loss is equivalent to VoxelMorph-1. This 3D network consists of encoder, decoder and their skip connections similar to U-Net [9].

The 3D spatial transformation layer [6] is to deform the moving volume with the deformation field \(\phi \). We use the spatial transformation function \(\mathcal {T}\) with tri-linear interpolation for warping the image A by \(\phi \), which can be written as:

$$\begin{aligned} \mathcal {T}(A, \phi ) = \sum \nolimits _{y \in \mathcal {N}(x+\phi (x))}^{ }A(y)\prod \nolimits _{d \in \{i,j,k\}}^{ }(1-|x_d + \phi (x_d)-y_d|), \end{aligned}$$
(8)

where x indicates the voxel index, \(\mathcal {N}(x+\phi (x))\) denotes the 8-voxel cubic neighborhood around \(x+\phi (x)\), and d is three directions in 3D image space.

3 Experiments

To verify the performance of our method, we conducted liver registration from multiphase CT images. The dataset was collected from liver cancer (HCC) patients at Asan Medical Center, Seoul, South Korea. Each scans has pathologically proven hepatic nodules and four-phase liver CT (unenhanced, arterial, portal, and 180-s delayed phases). The slice thickness was 5 mm. We did not perform pre-processing such as affine transformation except for matching the number of slices for the moving and fixed images. Here, we extracted slices only including liver by a pre-trained liver segmentation network and performed zero-padding to the above and below volumes based on the center of mass of liver.

We used 555 scans for training and 50 scans for testing. For the network training, we stacked two volumes with different phases as the input. We normalized the input intensity with the maximum value of each volume. Also, we randomly down-sampled the training data from \(512\times 512\times depth\) to \(128\times 128\times depth\) to fit in the GPU memory, while we evaluated the test data with original size of \(512\times 512\times depth\), where depth is different for each pair of input. For data augmentation, we adopted random horizontal/vertical flipping and rotation with 90\(^\circ \) for each pair of training volume. Our proposed method was implemented with pyTorch library. We applied Adam with momentum optimization algorithm to train the models with a learning rate of 0.0001, and set the batch size 1. The model was trained for 50 epochs using a NVIDIA GeForce GTX 1080 Ti GPU.

Fig. 3.
figure 3

Results of the target registration errors (TRE) of all 20 anatomical points in the deformed arterial and delayed images of 50 test data. Mean graph represents the mean TRE of the points for all subjects. D# in the x-axis indicates the patient number.

Full size image
Table 1. Tumor size differences, TRE values (mm) between the deformed arterial/delayed images and the fixed portal image, and their average time (min) to be deformed on test set.
Full size table

3.1 Registration Results

We evaluated the registration performance using the target registration error (TRE) based on the 20 anatomical points in the liver and adjacent organs at the portal phases, which are marked by radiologists. Also, we measured the tumor size that verifies the registration performance in the view point of tumor diagnosis. We compared our method to Elastix [8] that is known to have its state-of-the-art performance among the classical approaches. We also compared with VoxelMorph-1 [2]. Additionally, we performed ablation studies by excluding the cycle loss or identity loss. Apart from the loss, different ablated networks were subjected to the same training procedure for fair comparison.

Table 2. Percentage of voxels with a non-positive Jacobian determinant and normalized mean square error (NMSE) on test set.
Full size table
Fig. 4.
figure 4

Results of multiphase liver CT registration (Left) and their deformation fields (Right). The diagonal images with red-box are original images, which are deformed to other phases as indicated by each row. Specifically, the \((i,j), i\ne j\) element of the figure represents the deformed image to the i-th phase from the j-th phase original image. (Color figure online)

Full size image

Figure 3 shows the registration performance. We visualize the TRE values of the deformed arterial and delayed images with respect to each test data, and also show the average TRE values of all subjects with respect to the deformed arterial and delayed images into the fixed portal image. We can observe that the proposed method achieves significant improvement compared to VoxelMorph-1, while the error of the proposed method is slightly higher than Elastix. Also, in Table 1, we can confirm that the tumor size of deformed images from our proposed method is the most accurate for the case of delay to portal registration, and comparable in arterial to portal registration.

To demonstrate the effect of the cycle consistency, we also computed the percentage of voxels with a non-positive Jacobian determinant on the deformation fields and the normalized mean square error (NMSE) between the original moving image and re-deformed image. As shown in Table 2, we confirm that the proposed method is less prone to folding problem and improves topological preservation for liver registration.

Figure 4 illustrates an example of registration results that deforms the multiphase 3D images with the four distinct phases. Moreover, we calculated the time required for the proposed method to deform one image into the fixed image (see Table 1). Specifically, the conventional Elastix takes approximately 19.6 min for the image registration, while the proposed method takes only 10 s.

4 Conclusion

We presented an unsupervised image registration method using a cycle consistent convolutional neural network. Using two registration networks, our proposed method is trained to satisfy the cycle consistency that imposes inverse consistency between a pair of images. However, once the networks are trained, a single network can provide accurate 3D image registration with any pair of new data, so the computational complexity is same as VoxelMorph-1. Our liver registration results demonstrated that the proposed method works well for any image pairs with different contrast.