Recursive Cascaded Networks for Unsupervised Medical Image Registration

TL;DR

This paper introduces recursive cascaded networks for unsupervised deformable medical image registration, improving alignment accuracy by successively refining the deformation through a simple, end-to-end trainable architecture.

Contribution

The proposed recursive cascaded network architecture is novel in enabling progressive deformation learning for image registration, applicable to any base network and trained end-to-end.

Findings

Achieves significant performance improvements over state-of-the-art methods.

Performance improves with more cascades, with no observed limit.

Effective on 3D medical images for deformable registration.

Abstract

We present recursive cascaded networks, a general architecture that enables learning deep cascades, for deformable image registration. The proposed architecture is simple in design and can be built on any base network. The moving image is warped successively by each cascade and finally aligned to the fixed image; this procedure is recursive in a way that every cascade learns to perform a progressive deformation for the current warped image. The entire system is end-to-end and jointly trained in an unsupervised manner. In addition, enabled by the recursive architecture, one cascade can be iteratively applied for multiple times during testing, which approaches a better fit between each of the image pairs. We evaluate our method on 3D medical images, where deformable registration is most commonly applied. We demonstrate that recursive cascaded networks achieve consistent, significant gains and outperform state-of-the-art methods. The performance reveals an increasing trend as long as more cascades are trained, while the limit is not observed. Code is available at https://github.com/microsoft/Recursive-Cascaded-Networks.

Tables

Table 1. Table 1: Comparison among traditional methods (ANTs SyN and Elastix B-spline), our baseline networks (VoxelMorph and VTN), and our proposed recursive cascaded networks with and without shared-weight cascading. r 𝑟 r × \times n 𝑛 n -cascade means that every deformable cascade is repetitively applied for r 𝑟 r times during testing, using our proposed shared-weight cascading method. For liver datasets (SLIVER, LiTS, and LSPIG), the Dice score measures the overlap of liver segmentations, and Lm. Dist. means an average distance among 6 annotated landmarks. Avg. Dice means an average Dice score among all 56 segmented anatomical structures for the brain dataset LPBA. Standard deviations across instances are in parentheses.

Method	SLIVER		LiTS	LSPIG	LPBA	Time (sec)
	Dice	Lm. Dist.	Dice	Dice	Avg. Dice	GPU	CPU
ANTs SyN [4, 5]	0.895 (0.037)	12.2 (5.7)	0.862 (0.055)	0.825 (0.063)	0.708 (0.015)	-	748
Elastix B-spline [35, 48]	0.910 (0.038)	12.6 (6.6)	0.863 (0.059)	0.825 (0.059)	0.675 (0.013)	-	115
VoxelMorph1  [9]	0.883 (0.034)	14.0 (4.6)	0.831 (0.061)	0.715 (0.090)	0.685 (0.017)	0.20	17
VoxelMorph (reimplem.)2	0.913 (0.025)	13.1 (4.7)	0.870 (0.048)	0.833 (0.057)	0.688 (0.015)	0.15	14
5-cascade VoxelMorph	0.944 (0.017)	12.4 (4.9)	0.903 (0.055)	0.849 (0.062)	0.708 (0.015)	0.41	69
3$\times$5-cascade VoxelMorph	0.950 (0.014)	11.9 (4.9)	0.905 (0.065)	0.842 (0.066)	0.715 (0.014)	1.09	201
VTN (ADDD)3  [37]	0.942 (0.020)	12.0 (4.9)	0.897 (0.049)	0.846 (0.064)	0.701 (0.014)	0.13	26
10-cascade VTN	0.953 (0.014)	10.8 (4.9)	0.909 (0.060)	0.855 (0.060)	0.716 (0.013)	0.25	87
2$\times$10-cascade VTN	0.956 (0.012)	10.2 (4.7)	0.908 (0.070)	0.849 (0.063)	0.719 (0.012)	0.42	179

Table 2. Table 2: Comparison among different number of recursive cascades. n 𝑛 n -cascade means n 𝑛 n recursive cascades of the base network, excluding the affine cascade. Standard deviations across instances are in parentheses.

Architecture	SLIVER		LiTS	LSPIG	LPBA	Time (sec)
	Dice	Lm. Dist.	Dice	Dice	Avg. Dice	GPU	CPU
Affine only	0.794 (0.042)	14.8 (4.7)	0.754 (0.059)	0.727 (0.054)	0.628 (0.017)	0.08	0.4
1-cascade VoxelMorph	0.913 (0.025)	13.1 (4.7)	0.867 (0.050)	0.833 (0.057)	0.688 (0.015)	0.15	14
2-cascade VoxelMorph	0.933 (0.021)	12.8 (4.8)	0.888 (0.048)	0.845 (0.057)	0.699 (0.014)	0.21	27
3-cascade VoxelMorph	0.940 (0.018)	12.6 (5.0)	0.897 (0.049)	0.849 (0.060)	0.706 (0.014)	0.28	40
4-cascade VoxelMorph	0.943 (0.017)	12.5 (5.1)	0.900 (0.052)	0.851 (0.058)	0.707 (0.014)	0.35	54
5-cascade VoxelMorph	0.944 (0.017)	12.4 (4.9)	0.903 (0.055)	0.849 (0.062)	0.708 (0.015)	0.41	69
1-cascade VTN	0.914 (0.025)	13.0 (4.8)	0.870 (0.048)	0.833 (0.054)	0.686 (0.014)	0.10	10
2-cascade VTN	0.935 (0.020)	12.2 (4.7)	0.891 (0.045)	0.843 (0.061)	0.697 (0.014)	0.12	18
3-cascade VTN	0.943 (0.018)	11.8 (4.7)	0.900 (0.045)	0.850 (0.060)	0.703 (0.014)	0.13	26
4-cascade VTN	0.948 (0.016)	11.6 (4.8)	0.906 (0.047)	0.852 (0.063)	0.708 (0.014)	0.15	35
5-cascade VTN	0.949 (0.015)	11.5 (4.8)	0.908 (0.051)	0.853 (0.064)	0.709 (0.014)	0.17	47
6-cascade VTN	0.951 (0.015)	11.3 (4.9)	0.910 (0.050)	0.852 (0.064)	0.712 (0.014)	0.18	57
7-cascade VTN	0.951 (0.015)	11.2 (4.9)	0.908 (0.055)	0.852 (0.061)	0.712 (0.013)	0.20	65
8-cascade VTN	0.952 (0.014)	11.1 (4.7)	0.910 (0.056)	0.854 (0.059)	0.714 (0.013)	0.22	75
9-cascade VTN	0.953 (0.014)	10.9 (4.7)	0.910 (0.059)	0.851 (0.064)	0.716 (0.013)	0.23	90
10-cascade VTN	0.953 (0.014)	10.8 (4.9)	0.909 (0.060)	0.855 (0.060)	0.716 (0.013)	0.25	87

Table 3. Table 3: Results of recursive cascaded networks built on the base network VTN, with different times (1 × \times , 2 × \times , 3 × \times , or more) of shared-weight cascades. Similarity is measured by the correlation coefficient between the warped moving image and the fixed image.

Cascade	SLIVER			LiTS		LSPIG		LPBA
	Dice	Lm. Dist.	Similarity	Dice	Similarity	Dice	Similarity	Avg. Dice	Similarity
1$\times$1	0.914 (0.025)	13.0 (4.8)	0.7458 (0.0396)	0.870 (0.048)	0.7386 (0.0468)	0.833 (0.054)	0.7527 (0.0515)	0.686 (0.014)	0.9814 (0.0021)
2$\times$1	0.932 (0.020)	12.6 (5.0)	0.8108 (0.0289)	0.886 (0.048)	0.8045 (0.0376)	0.840 (0.057)	0.8162 (0.0392)	0.694 (0.014)	0.9845 (0.0016)
3$\times$1	0.937 (0.019)	12.5 (5.1)	0.8333 (0.0248)	0.888 (0.050)	0.8272 (0.0336)	0.839 (0.057)	0.8369 (0.0338)	0.695 (0.013)	0.9854 (0.0014)
4$\times$1	0.938 (0.018)	12.5 (5.2)	0.8444 (0.0227)	0.887 (0.053)	0.8381 (0.0314)	0.837 (0.057)	0.8467 (0.0305)	0.692 (0.013)	0.9857 (0.0011)
5$\times$1	0.939 (0.018)	12.5 (5.2)	0.8510 (0.0214)	0.886 (0.056)	0.8446 (0.0300)	0.835 (0.058)	0.8518 (0.0289)	0.686 (0.013)	0.9857 (0.0010)
1$\times$2	0.935 (0.020)	12.2 (4.7)	0.8270 (0.0297)	0.891 (0.045)	0.8209 (0.0367)	0.843 (0.061)	0.8435 (0.0369)	0.697 (0.014)	0.9854 (0.0017)
2$\times$2	0.947 (0.017)	11.6 (4.8)	0.8779 (0.0198)	0.900 (0.049)	0.8715 (0.0282)	0.847 (0.063)	0.8919 (0.0243)	0.701 (0.014)	0.9885 (0.0011)
3$\times$2	0.948 (0.016)	11.5 (4.8)	0.8930 (0.0171)	0.900 (0.054)	0.8865 (0.0254)	0.845 (0.063)	0.9039 (0.0211)	0.697 (0.014)	0.9895 (0.0008)
1$\times$3	0.943 (0.018)	11.8 (4.7)	0.8584 (0.0245)	0.900 (0.045)	0.8535 (0.0318)	0.850 (0.060)	0.8774 (0.0282)	0.703 (0.014)	0.9876 (0.0014)
2$\times$3	0.951 (0.015)	11.2 (4.8)	0.8977 (0.0168)	0.905 (0.052)	0.8927 (0.0246)	0.852 (0.061)	0.9102 (0.0210)	0.710 (0.014)	0.9904 (0.0009)
3$\times$3	0.951 (0.015)	11.1 (4.9)	0.9088 (0.0146)	0.904 (0.058)	0.9037 (0.0225)	0.850 (0.062)	0.9189 (0.0188)	0.711 (0.014)	0.9916 (0.0007)
1$\times$5	0.949 (0.015)	11.5 (4.8)	0.8926 (0.0186)	0.908 (0.051)	0.8893 (0.0254)	0.853 (0.063)	0.9088 (0.0223)	0.709 (0.014)	0.9894 (0.0010)
2$\times$5	0.954 (0.013)	10.8 (4.9)	0.9215 (0.0131)	0.908 (0.061)	0.9184 (0.0198)	0.851 (0.063)	0.9334 (0.0164)	0.715 (0.013)	0.9921 (0.0006)
3$\times$5	0.954 (0.013)	10.6 (5.0)	0.9308 (0.0115)	0.906 (0.067)	0.9278 (0.0182)	0.845 (0.065)	0.9406 (0.0145)	0.715 (0.013)	0.9930 (0.0005)
1$\times$10	0.953 (0.014)	10.8 (4.9)	0.9163 (0.0145)	0.909 (0.060)	0.9129 (0.0211)	0.855 (0.059)	0.9290 (0.0174)	0.716 (0.013)	0.9918 (0.0008)
2$\times$10	0.956 (0.012)	10.2 (4.7)	0.9384 (0.0106)	0.908 (0.070)	0.9355 (0.0171)	0.849 (0.062)	0.9471 (0.0132)	0.719 (0.012)	0.9942 (0.0005)
3$\times$10	0.956 (0.012)	10.2 (4.7)	0.9461 (0.0094)	0.905 (0.076)	0.9434 (0.0158)	0.841 (0.068)	0.9534 (0.0112)	0.717 (0.012)	0.9951 (0.0004)

Table 4. Table 4: Comparison against other variants of VoxelMorph (VM), including VM x2 (doubling the feature counts of every convolutional layer), VM-double (doubling the number of convolutional layers at each level), and VM xx2 (doubling the encoder-decoder architecture cascade-like).

Architecture	SLIVER		LiTS	LSPIG	LPBA
	Dice	Lm. Dist.	Dice	Dice	Avg. Dice
VoxelMorph	0.913 (0.025)	13.1 (4.7)	0.867 (0.050)	0.833 (0.057)	0.688 (0.015)
VM x2	0.922 (0.024)	13.0 (4.9)	0.879 (0.047)	0.839 (0.058)	0.691 (0.015)
VM-double	0.919 (0.025)	12.9 (4.9)	0.877 (0.048)	0.833 (0.059)	0.689 (0.015)
VM xx2	0.925 (0.023)	12.8 (4.9)	0.881 (0.047)	0.843 (0.057)	0.693 (0.014)
2$\times$1-cascade VM	0.930 (0.021)	12.8 (4.8)	0.883 (0.051)	0.840 (0.060)	0.697 (0.014)
2-cascade VM	0.933 (0.021)	12.8 (4.8)	0.888 (0.048)	0.845 (0.057)	0.699 (0.014)

Table 5. Table 5: Experiments on different atlases in LPBA.

Method	Avg. Dice
	Atlas1	Atlas2	Atlas3
ANTs SyN	0.708 (0.015)	0.717 (0.011)	0.707 (0.015)
Elastix B-spline	0.675 (0.013)	0.684 (0.011)	0.670 (0.013)
VoxelMorph	0.688 (0.015)	0.694 (0.010)	0.678 (0.015)
5-cascade VoxelMorph	0.708 (0.015)	0.714 (0.011)	0.702 (0.014)
3$\times$5-cascade VoxelMorph	0.715 (0.014)	0.721 (0.012)	0.713 (0.013)
VTN (ADDD)	0.701 (0.014)	0.709 (0.011)	0.695 (0.015)
10-cascade VTN	0.716 (0.013)	0.723 (0.010)	0.712 (0.013)
2$\times$10-cascade VTN	0.719 (0.012)	0.725 (0.011)	0.716 (0.013)

Equations

$$I_m^{\prime }=\varphi \circ I_m.$$

$$F(I_m,I_f)=\varphi \circ F_1(\varphi \circ I_m,I_f),$$

$${\varphi }_k=f_k(I_m^{(k-1)},I_f).$$

$$F(I_m,I_f)={\varphi }_n\circ \cdots \circ {\varphi }_1,$$

$$I_m^{(n)}=F(I_m,I_f)\circ I_m.$$

$$\mathrm{Dice}(A,B)=2\cdot \frac{|A\cap B|}{|A|+|B|}.$$

Full text

1 Introduction

Deformable image registration has been studied in plenty of works and raised great importance. The non-linear correspondence between a pair of images is established by predicting a deformation field under the smoothness constraint. Among traditional algorithms, an iterative approach is commonly suggested [2, 3, 4, 7, 10, 18, 27, 52], where they formulate each iteration as a progressive optimization problem.

Image registration has drawn growing interests in terms of deep learning techniques. A closely related area is optical flow estimation, which is essentially a 2D image registration problem but the flow fields are discontinuous across objects and the tracking is mainly about motion with rare color difference. Occlusions and folding areas requiring a guess are inevitable in optical flow estimation (but not expected in deformable image registration). Automatically generated datasets (e.g., Flying Chairs [24], Flying Things 3D [41]) are of great help for supervising convolutional neural networks (CNNs) in such settings [24, 29, 30, 54, 55]. Some studies also try to stack multiple networks. They assign different tasks and inputs to each cascade in a non-recursive way and train them one by one [30, 45], but their performance approaches a limit with only a few (no more than 3) cascades. On the other hand, cascading may not help much when dealing with discontinuity and occlusions. Thus by intuition, we suggest that cascaded networks with a recursive architecture fits the setting of deformable registration.

Learning-based methods are also suggested as an approach in deformable image registration. Unlike optical flow estimation, intersubject registration with vague correspondence of image intensity is usually demanded. Some initial works rely on the dense ground-truth flows obtained by either traditional algorithms [14, 56] or simulating intrasubject deformations [36, 53], but their performance is restricted due to the limited quality of training data.

Unsupervised learning methods with comparable performance to traditional algorithms have been presented recently [8, 9, 19, 20, 37, 38]. They only require a similarity measurement between the warped moving image and the fixed image, while the gradients can backpropagate through the differentiable warping operation (a.k.a. spatial transformer [32]). However, most proposed networks are enforced to make a straightforward prediction, which proves to be a burden when handling complicated deformations especially with large displacements. DLIR [19] and VTN [37] also stack their networks, though both limited to a small number of cascades. DLIR trains each cascade one by one, i.e., after fixing the weights of previous cascades. VTN jointly trains the cascades, while all successively warped images are measured by the similarity compared to the fixed image. Neither training method allows intermediate cascades to progressively register a pair of images. Those non-cooperative cascades learn their own objectives regardless of the existence of others, and thus further improvement can hardly be achieved even if more cascades are conducted. They may realize that network cascading possibly solves this problem, but there is no effective way of training deep network cascades for progressive alignments.

Therefore, we propose the recursive cascade architecture, which encourages the unsupervised training of an unlimited number of cascades that can be built on existing base networks, for advancing the state of the art. The difference between our architecture and existing cascading methods is that each of our cascades commonly takes the current warped image and the fixed image as inputs (in contrast to [30, 45]) and the similarity is only measured on the final warped image (in contrast to [19, 37]), enabling all cascades to learn progressive alignments cooperatively. Figure 1 shows an example of applying the proposed architecture built on 10 deformable cascades of the base network VTN.

Conceptually, we formulate the registration problem as determining a parameterized flow prediction function, which outputs a dense flow field based on the input of an image pair. This function can be recursively defined on the warped moving image with essentially the same functionality. Instead of training the function in a straightforward way, the final prediction can be considered a composition of recursively predicted flow fields, while each cascade only needs to learn a simple alignment of small displacement that can be refined by deeper recursion. Figure 2 verifies our conception. Our method also enables the use of shared-weight cascades, which potentially achieves performance gains without introducing more parameters.

To summarize, we present a deep recursive cascade architecture for deformable image registration, which facilitates the unsupervised end-to-end learning and achieves consistent gains independently of the base network; shared-weight cascading technique with direct test-time improvement is developed as well. We conduct extensive experiments based on diverse evaluation metrics (segmentations and landmarks) and multiple datasets across image types (liver CT scans and brain MRIs).

2 Related Work

Cascade approaches have been involved in a variety of domains of computer vision, e.g., cascaded pose regression progressively refines a pose estimation learned from supervised training data [23], cascaded classifiers speed up the process of object detection [25].

Deep learning also benefits from cascade architectures. For example, deep deformation network [57] cascades two stages and predicts a deformation for landmark localization. Other applications include object detection [13], semantic segmentation [17], and image super-resolution [16]. There are also several works specified to medical images, e.g., 3D image reconstruction for MRIs [6, 49], liver segmentation [46] and mitosis detection [15]. Note that shallow, non-recursive network cascades are usually proposed in those works.

In respect of registration, traditional algorithms iteratively optimize some energy functions in common [2, 3, 4, 7, 10, 18, 27, 52]. Those methods are also recursive in general, i.e., similarly functioned alignments with respect to the current warped images are performed during iterations. Iterative Closest Point is an iterative, recursive approach for registering point clouds [12, 58], where the closest pairs of points are matched at each iteration and a rigid transformation that minimizes the difference is solved. In deformable image registration, most traditional algorithms basically works like this but in a much more complex way. Standard symmetric normalization (SyN) [4] maximizes the cross-correlation within the space of diffeomorphic maps during iterations. Optimizing free-form deformations using B-spline [48] is another standard approach.

Learning-based methods are presented recently. Supervised methods entail much effort on the labeled data that can hardly meet the realistic demands, resulting in the limited performance [14, 56, 36, 53]. Unsupervised methods are proposed to solve this problem. Several initial works shows the possibility of unsupervised learning [19, 20, 38, 50], among which DLIR [20] performs on par with the B-spline method implemented in SimpleElastix [40] (a multi-language extension of Elastix [35], which is selected as one of our baseline methods). VoxelMorph [8] and VTN [37] achieve better performance by predicting a dense flow field using deconvolutional layers [44], whereas DLIR only predicts a sparse displacement grid interpolated by a third order B-spline kernel. VoxelMorph only evaluates their method on brain MRI datasets [8, 9], but shown deficiency on other datasets such as liver CT scans by later work [37]. Additionally, VTN proposes an initial convolutional network which performs an affine transformation before predicting deformation fields, leading to a truly end-to-end framework by substituting the traditional affine stage.

State-of-the-art VTN and VoxelMorph are selected as our base networks, and the suggested affine network is also integrated as our top-level cascade. To our knowledge, none of those work realizes that training deeper cascades advances the performance for deformable image registration.

3 Recursive Cascaded Networks

Let $I_{m},I_{f}$ denote the moving image and the fixed image respectively, both defined over $d$-dimensional space $\Omega$. A flow field is a mapping $\phi:\Omega\to\Omega$. For deformable image registration, a reasonable flow field should be continuously varying and prevented from folding. The task is to construct a flow prediction function $F$ which takes $I_{m},I_{f}$ as inputs and predicts a dense flow field that aligns $I_{m}$ to $I_{f}$.

We cascade this procedure by recursively performing registration on the warped image. The warped image $I_{m}^{\prime}$ is exactly the composition of the flow field and the moving image, namely

[TABLE]

Conceptually,

[TABLE]

where $F_{1}$ may be the same as $F$, but in general a different flow prediction function. This recursion can be infinitely applied in theory.

Following this recursion, the moving image is warped successively, enabling the final prediction (probably with large displacement) to be decomposed into cascaded, progressive refinements (with small displacements). One cascade is basically a flow prediction function ($f_{k}$), and the $k$-th cascade predicts a flow field of

[TABLE]

$I_{m}^{(k)}$ denotes the moving image warped by the first $k$ cascades. Figure 3 depicts the proposed architecture. Assuming for $n$ cascades in total, the final output is a composition of all predicted flow fields, i.e.,

[TABLE]

and the final warped image is constructed by

[TABLE]

3.1 Subnetworks

Each $f_{k}$ is implemented as a convolutional neural network in this paper. Every network is designed to predict a deformable flow field on itself based on the input warped image and the fixed image. $f_{1},\dots f_{n}$ can be different in network architecture, but surely using a common base network is well-designed enough for convenience. Those cascades may learn different network parameters on each, since one cascade is allowed to learn a part of measurements or perform some type of alignment specifically. Note that the images input to the networks are discretized and so are the output flow fields, thus we treat them by multilinear interpolation (or simply trilinear interpolation for 3D images), and out-of-bound indices by nearest-point interpolation [37].

An architecture similar to the U-Net [31, 47] is widely used for deformable registration networks, such as VTN [37] and VoxelMorph [8]. Such network consists of encoders followed by decoders with skip connections. The encoders help to extract features, while the decoders perform upsampling and refinement, ending with a dense prediction.

For medical images, it is usually the case that two scans can be roughly aligned by an initial rigid (or affine) transformation. VoxelMorph [8] assumes that input images are pre-affined by an external tool, whereas VTN [37] integrates an efficient affine registration network which outperforms the traditional stage. As a result, we also embed the affine registration network as our top-level cascade, which behaves just like a normal one except that it is only allowed to predict an affine transformation rather than general flow fields.

3.2 Unsupervised End-to-End Learning

We suggest that all cascades can be jointly trained by merely measuring the similarity between $I_{m}^{(n)}$ and $I_{f}$ together with regularization losses. Enabled by the differentiable composition operator (i.e., warping operation), recursive cascaded networks can learn to perform progressive alignments cooperatively without supervision. To our knowledge, no previous work achieves good performance by stacking more than 3 deformable registration networks, partly because they train them one by one [19] (then the performance can hardly improve) or they measure the similarity on each of the warped images [37] (then the networks can hardly learn progressive alignments).

Regularization losses are basically the smooth terms of $\phi_{1},\dots,\phi_{n}$, and thus are necessary. Every predicted flow field is penalized by an L2 variation loss as done in [8, 37]. The affine cascade works with its own regularization losses introduced in VTN [37].

3.3 Shared-Weight Cascading

One cascade can be repetitively applied during recursion. I.e., multiple cascades can be shared with the same parameters, and that is called shared-weight cascading.

After an $n$-cascade network is trained, we can still possibly apply additional shared-weight cascades during testing. For example, we may replicate all cascades as an indivisible whole by the end of $I_{m}^{(n)}$, i.e., totally $2n$ cascades are associated with flow prediction functions $f_{1},\dots,f_{n},f_{1},\dots,f_{n}$ respectively. We develop a better approach by immediately inserting one or more shared-weight cascades after each, i.e., totally r$$\times$$n cascades are constructed by substituting each $f_{k}$ by $r$ times of that. This approach will be proved to be effective later in the experiments.

Shared-weight cascading during testing is an option when the quality of output flow fields can be improved by further refinement. However, we note that this technique does not always get positive gains and may lead to over deformation. Recursive cascades only ensure an increasing similarity between the warped moving image and the fixed image, but the aggregate flow field becomes less natural if the images are too perfectly matched.

The reason we do not use shared-weight cascading in training is that shared-weight cascades consume extra GPU memory as large as non-shared-weight cascades during gradient backpropagation in the platform we use (Tensorflow [1]). The number of cascades to train is constrained by the GPU memory, but they would perform better with the allowance of learning different parameters when the dataset is large enough to avoid overfitting.

4 Experiments

4.1 Experimental Settings

We build our recursive cascaded networks mainly based on the network architecture of VTN [37], which is a state-of-the-art method for deformable image registration. Note that VTN already stacks a few cascades of their deformable subnetworks, and a single cascade is being used as our base network. Up to 10-cascade VTN (excluding the affine cascade) is jointly trained using our proposed method. To show the generalizability of our architecture, we also choose VoxelMorph [9] as another base network. We train up to 5-cascade VoxelMorph, because each cascade of VoxelMorph consumes more resources.

We evaluate our method on two types of 3D medical images: liver CT scans and brain MRI scans. For liver CT scans, we train and test recursive cascaded networks for pairwise, subject-to-subject registration, which stands for a general purpose of allowing the fixed image to be arbitrary. For brain MRI scans, we follow the experimental setup of VoxelMorph [8], where each moving image is registered to a fixed atlas, called atlas-based registration. Both settings are common in medical image registration.

Implementation.

Inherited from the implementation of VTN [37] using Tensorflow 1.4 [1] built with a custom warping operation, the correlation coefficient is used as the similarity measurement, while the ratios of regularization losses are kept the same as theirs. We train our model using a batch size of $4$, on 4 cards of 12G NVIDIA TITAN Xp GPU. The training stage runs for $10^{5}$ iterations with the Adam optimizer [33]. The learning rate is initially $10^{-4}$ and halved after $6\times 10^{4}$ steps and again after $8\times 10^{4}$ steps.

Baseline Methods.

VTN [37] and VoxelMorph [8] are state-of-the-art learning-based methods. We cascade their base networks and also compare with the original systems. Besides, we also compare against SyN [4] (integrated in ANTs [5] together with the affine stage) and B-spline [48] (integrated in Elastix [35] together with the affine stage), which are shown to be the top-performing traditional methods for deformable image registration [8, 34, 37]. We run ANTs SyN and Elastix B-spline with the parameters recommended in VTN [37].

Evaluation Metrics.

We quantify the performance by the Dice score [22] based on the segmentation of some anatomical structure, between the warped moving image and the fixed image, as done in [8, 19]. The Dice score of two regions $A,B$ is formulated as

[TABLE]

Perfectly overlapped regions come with a Dice score of $1$. The Dice score explicitly measures the coincidence between two regions and thereby reflects the quality of registration. If multiple anatomical structures are annotated, we compute the Dice score with respect to each and take an average.

In addition, landmark annotations are available in some datasets and can be utilized as an auxiliary metric. We compute the average distance between the landmarks of the fixed image and the warped landmarks of the moving image, also introduced in VTN [37].

4.2 Datasets

For liver CT scans, we use the following datasets:

• •

MSD [42]. This dataset contains various types of medical images for segmenting different target objects. CT scans of liver tumours (70 scans excluding LiTS), hepatic vessels (443 scans), and pancreas tumours (420 scans) are selected since liver is likely to be included.

• •

BFH (introduced in VTN [37]), 92 scans.

• •

SLIVER [28], 20 scans with liver segmentation ground truth. Additionally, 6 anatomical keypoints selected as landmarks are annotated by 3 expert doctors, and we take their average as ground truth.

• •

LiTS [39], 131 scans with liver segmentation ground truth.

• •

LSPIG (Liver Segmentation of Pigs, provided by the First Affiliated Hospital of Harbin Medical University), containing 17 pairs of CT scans from pigs, along with liver segmentation ground truth. Each pair comes from one pig with (perioperative) and without (preoperative) 13 mm Hg pneumoperitoneum pressure.

Unsupervised methods are trained on the combination of MSD and BFH with $1025^{2}$ ($1025=70+443+420+92$) image pairs in total. SLIVER ($20\times 19$ image pairs) and LiTS ($131\times 130$ image pairs) are used for regular evaluation, while LSPIG is regarded as a challenging dataset which entails generalizability. Only 34 intrasubject image pairs in LSPIG, each of which comes from a same pig (preoperative to perioperative, or vice versa), are evaluated.

For brain MRI scans, we use the following datasets:

• •

ADNI [43], 66 scans.

• •

ABIDE [21], 1287 scans.

• •

ADHD [11], 949 scans.

• •

LPBA (LONI Probabilistic Brain Atlas) [51]. This dataset contains 40 scans, each of which comes with segmentation ground truth of 56 anatomical structures.

ADNI, ABIDE, ADHD are used for training, and LPBA for testing. All 56 anatomical structures are evaluated by an average Dice score. For atlas-based registration, the first scan in LPBA is fixed as the atlas in our experiments, which is shown to be without loss of generality later in the atlas analysis.

We carry out standard preprocessing steps referring to VTN [37] and VoxelMorph [8]. Raw scans are resampled into $128\times 128\times 128$ voxels after cropping unnecessary area around the target object. For liver CT scans, a simple threshold-based algorithm is applied to find a rough liver bounding box for cropping. For brain MRI scans, skulls are first removed using FreeSurfer [26]. The volumes are visualized for quality control so that seldom badly processed images are manually removed. (An overview of the evaluation datasets is provided in the supplementary material.)

4.3 Results

Table 1 summarizes our overall performance compared with state-of-the-art methods. Running times are approximately the same across datasets, so we test them on SLIVER, with an NVIDIA TITAN Xp GPU and an Intel Xeon E5-2690 v4 CPU. No GPU implementation of ANTs or Elastix has been found, nor in previous works [5, 8, 19, 35, 37]. Figure 4 visualizes those methods on an example in the brain dataset LPBA. (See the supplementary material for more examples.)

As shown in Table 1, recursive cascaded networks outperform the existing methods in all our datasets with significant gains. More importantly, the proposed architecture is independent of the base network, not limited to VTN [37] and VoxelMorph [8]. Although the number of cascades causes linear increments to the running times, a 10-cascade VTN still runs in a comparable (GPU) time to the baseline networks, showing the efficiency of our architecture.

Number of Cascades.

Table 2 presents the results with respect to different number of recursive cascades, choosing either VTN or VoxelMorph as our base network. As shown in the table, recursive cascaded networks achieve consistent performance gains independently of the base network. Our 3-cascade VTN (in Table 2) already outperforms VTN (ADDD) (in Table 1) although they have similar network architectures, mainly because our intermediate cascades learn progressive alignments better with only the similarity loss drawn on the final warped image. Figure 5 plots our results for better illustrating the increasing trend. Note that our architecture requires a linear time increment, but cascading a small-size base network like VTN is quite efficient.

Shared-Weight Cascading.

Deeper cascades can be directly constructed using weight sharing. As we suggest, an r$$\times$$n-cascade network successively repeats each of the jointly trained $n$ cascades for $r$ times during testing. A linear time increment is also required. This technique ensures an increasing similarity between the warped moving image and the fixed image, but we note that it does not always get positive performance gains.

Table 3 presents the results of shared-weight cascaded networks, together with the image similarity (correlation coefficient is used in this paper). The image similarity is always increasing as we expect. Shallower cascaded networks benefit more from this technique relatively to the deeper ones, since the images are still not well-registered (with relatively low similarity, as shown in the table). Less expected results on LiTS and LSPIG datasets may imply that this additional technique has a limited generalizability.

Note that shared-weight cascades generally perform worse than their jointly trained counterparts. More than 3 times of shared-weight cascades are very likely to deteriorate the quality (which partly coincides with previous studies), further proving the end-to-end learning to be vital.

Cascades vs. Channels vs. Depth.

VoxelMorph (VM) [9] suggests that the number of channels in the convolutional layers can be doubled for a better performance. We compare this variant (VM x2) against the jointly trained 2-cascade VM as well as a shared-weight 2$\times$1-cascade VM, shown in Table 4. VM x2 performs better than the original one as they suggest, but worse than both of our cascade methods. On the other hand, the number of parameters in VM x2 is 4 times as large as that in VoxelMorph (as well as 2$\times$1-cascade VM), and 2 times as large as that in 2-cascade VM.

However, one may wonder that whether simply deeper networks would do the trick. To this end, we construct VM-double by doubling the number of convolutional layers at each U-net level, and also an encoder-decoder-encoder-decoder architecture denoted VM xx2, which looks similar to a 2-cascade VM except the explicit warping. They have approximately the same amount of parameters compared to the 2-cascade VM, but are outperformed by a considerable margin. This experiment implies that our improvements are essentially based on the proposed recursive cascade architecture rather than simply introducing more parameters.

Atlas Analysis.

The performance for atlas-based registration may vary depending on the chosen atlas. As a comparison, we retrain the models on two more (the second and the third) atlases in the LPBA dataset, shown in Table 5. These results indicate that our performance is consistent and robust to the choice of atlas.

5 Discussion

Recursive cascaded networks are quite simple to implement, and also easy to train. We do not tune the ratios of losses when training more cascades, nor the training schedule, showing the robustness of our architecture. If more resources are available or a distributed learning platform is being used, we expect that the performance can be further improved by deeper cascades, and also, training or fine-tuning shared-weight cascades would be an alternative choice. A light-weight base network is also worth an exploration.

A possible limitation of this work would be on the smoothness of the composed field. Theoretically, recursive cascaded networks preserve the image topology as long as every subfield does. However, folding area is common in currently proposed methods and may be amplified during recursion, which brings challenges especially for the use of weight sharing techniques. This problem can be reduced by taking a careful look on the regularization terms, or designing a base network that guarantees invertibility.

6 Conclusion

We present a deep recursive cascade architecture and evaluate its performance in deformable medical image registration. Experiments based on diverse evaluation metrics demonstrate that this architecture achieves significant gains over state-of-the-art methods on both liver and brain datasets. With the superiority of good performance, the general applicability of the unsupervised method, and being independent of the base network, we expect that the proposed architecture can potentially be extended to all deformable image registration tasks.