Conv2Warp: An unsupervised deformable image registration with continuous convolution and warping

TL;DR

Conv2Warp introduces a novel continuous warping approach for deformable image registration that combines linear and non-linear convolutions with a learnable spline resampler, achieving more accurate and smooth alignments efficiently.

Contribution

It proposes a continuous warping method using combined convolutions and spline resampling, improving registration accuracy over traditional linear convolutional approaches.

Findings

Captures large non-linear deformations effectively

Reduces interpolation errors during registration

Achieves high accuracy with computational efficiency

Abstract

Recent successes in deep learning based deformable image registration (DIR) methods have demonstrated that complex deformation can be learnt directly from data while reducing computation time when compared to traditional methods. However, the reliance on fully linear convolutional layers imposes a uniform sampling of pixel/voxel locations which ultimately limits their performance. To address this problem, we propose a novel approach of learning a continuous warp of the source image. Here, the required deformation vector fields are obtained from a concatenated linear and non-linear convolution layers and a learnable bicubic Catmull-Rom spline resampler. This allows to compute smooth deformation field and more accurate alignment compared to using only linear convolutions and linear resampling. In addition, the continuous warping technique penalizes disagreements that are due to topological changes. Our experiments demonstrate that this approach manages to capture large non-linear deformations and minimizes the propagation of interpolation errors. While improving accuracy the method is computationally efficient. We present comparative results on a range of public 4D CT lung (POPI) and brain datasets (CUMC12, MGH10).

Tables

Table 1. Table 1: Evaluation on lung 4D CT dataset (POPI  [ 9 ] ). Dice(Jaccard) coefficients for 7 different volume pairs with larger inhalation and exhalation cyles. Grid-parameter search was done for all methods and only best performance is reported.

4D-CT	Methods
	Pre-align.	Demons	ANTS-SyN	SE	ConvNet	Conv2Warp
00-50	0.83(0.76)	0.87(0.83)	0.87(0.82)	0.89(0.85)	0.87(0.82)	0.89(0.83)
00-70	0.86(0.79)	0.90(0.85)	0.89(0.84)	0.92(0.88)	0.90(0.84)	0.91(0.85)
10-50	0.82(0.75)	0.86(0.82)	0.86(0.81)	0.88(0.85)	0.87(0.81)	0.88(0.83)
10-70	0.85(0.78)	0.88(0.83)	0.87(0.82)	0.88(0.84)	0.89(0.83)	0.90(0.84)
20-50	0.84(0.77)	0.87(0.83)	0.86(0.82)	0.91(0.86)	0.87(0.81)	0.89(0.85)
20-70	0.86(0.79)	0.89(0.84)	0.88(0.83)	0.88(0.84)	0.89(0.83)	0.90(0.84)
30-50	0.87(0.82)	0.90(0.86)	0.89(0.85)	0.91(0.87)	0.88(0.83)	0.90(0.86)
Mean	0.84(0.78)	0.88(0.83)	0.87(0.82)	0.90(0.85)	0.88(0.82)	0.90(0.84)
$\overline{t}$ (in s.)	-	23.51	216.45	416.78	2.94	2.94

Table 2. Table 2: Evaluation on T1-weighted MRI datasets. Mean Dice and Jaccard coefficients for two datasets with the state-of-the-art DIR methods and Conv2Warp (checkpoints trained on LBPA40 dataset) are shown. Only rigid registration was performed for pre-alignment prior to applying deformable registration.

Methods	Dataset				Avg.
	MGH10		CUMC12		time
	${\mu }_{dice}$	${\mu }_{jaccard}$	${\mu }_{dice}$	${\mu }_{jaccard}$	$\overline{t}$
Rigid	0.92 $\pm$ 0.011	0.86 $\pm$ 0.019	0.90 $\pm$ 0.030	0.82 $\pm$ 0.050	-
SimpleElastix (SE)	0.95 $\pm$ 0.003	0.90 $\pm$ 0.005	0.95 $\pm$ 0.002	0.91 $\pm$ 0.002	1008.2
ANTS (SyN-CC)	0.94 $\pm$ 0.007	0.90 $\pm$ 0.011	0.95 $\pm$ 0.010	0.90 $\pm$ 0.017	246.4
ConvNet	0.93 $\pm$ 0.008	0.88 $\pm$ 0.014	0.95 $\pm$0.007	0.91 $\pm$ 0.012	2.9
Conv2Warp	0.95 $\pm$ 0.007	0.90 $\pm$ 0.012	0.97$\pm$ 0.005	0.93 $\pm$ 0.008	2.9

Full text

1 Introduction

Image registration is required in many medical imaging applications - from multi-modal data fusion to inter- and intra-patient comparisons. Rigid or affine registration methods can be used to align source image to its target image by computing a simple transformation matrix, i.e., optimization of only few parameters. However, such an alignment does not model changes due to organ deformation, patient weight loss, or tumour shrinkage. In order to tackle these changes non-rigid or deformable image registration (DIR) methods are advised. Traditional DIR methods are usually mathematically complex to model and therefore require the optimization of a very larger number of parameters which makes them computationally expensive. Recently, deep learning based methods have shown tremendous success in tackling such complex problems in image registration. This is mainly due to two reasons: 1) By learning the model parameters of the neural network through optimizing a loss function directly from the data a wide range of image features are considered. 2) These methods make very efficient use of hardware acceleration as a result they require significantly less time than traditional methods.

Recently, supervised as well as unsupervised learning techniques for medical image registration have been proposed. Supervised techniques [6] rely on ground truth deformation vector fields (DVFs) computed by traditional methods. Here, a neural network is only used as a regressor for approximating traditional methods. This restricts their learning process and such training requires a large number of DVFs which are both time consuming to generate and carries the risk of mitigating unreliable deformations. On the other hand, unsupervised techniques [7, 1, 10] can learn to predict DVFs without requiring ground truth deformations. Linear convolutional neural networks (ConvNets) were used in [7, 1, 10]. Guha et al. [1] and Li et al. [7] used linear sampling techniques for upscaling the obtained DVFs. We argue that such models do not approximate complex and non-linear local deformations with a sufficient accuracy. These methods were evaluated on brain datasets where local deformations are limited. While, [10] proposed to use a B-spline as transformation model and interpolation method for the predicted DVFs and presented results on more complex cardiac cine MRI and lung datasets. However, in [10] only linear ConvNets were used. Despite of good local support, B-splines can lead to larger interpolation errors as they do not pass through data points (see Suppl. Mat. Section 1 for details).

With this paper, we introduce an unsupervised approach for deformable image registration that is capable of learning complex deformations and is able to compute smooth, accurate and plausible DVFs. We propose to: 1) Relax the fixed geometric constrain imposed by traditional convolutional filters by adding a series of deformable convolutional filters [3] combined with linear convolutional filters which allows to capture complex features. 2) Apply a learnt bicubic Catmull-Rom spline resampler to minimize error in DVF resampling, 3) Aggregate large deformations using a multiple scale warping strategy. 4) Finally, to impose further smoothness of the DVF using a $L_{2}$-norm regularization. We show that interpolation technique influence the learning phenomena. We also illustrate deformable convolutions as presented in this work through experiments. We have evaluated our model on publicly available lung CT and brain MRI datasets.

2 Method

As motivated, we now introduce the key elements to our Conv2Warp model (each highlighted in Fig. 1), and describe our constrained loss function.

Linear and Deformable Convolutional Network (LD-ConvNet)

LD-ConvNet (blue block in Fig. 1) consists of 5 linear convolutions and 3 deformable convolution layers. Each convolution layer is combined with a batch normalization and ELU activation function. After the second convolution layer an average pooling layer with downsampling factor of 0.5 is applied. Throughout we use a sequential ConvNet of (input, output) channels as follows: $(2,64)$, $(64,64)$, $(64,64)$, $(64,32)$, $(32,32)$, $(32,16)$, $(16,16)$, $(16,2)$ with kernel size $3$ and stride 1. Two images $I^{s}$ (source) and $I^{t}$ (target) are concatenated first and set as an input to the first layer of the Conv2Warp model.

Non-linear DVF Resampler (NL-DVF-R)

NL-DVF-R (green block in Fig. 1) consists of a sequential bicubic Catmull-Rom spline resampler and convolutional filter for resampling the obtained DVF from the LD-ConvNet. Since, this is integrated in our learning model, complex non-linear deformations learnt by the LD-ConvNet are guaranteed to have smooth deformation fields and the interpolation error is minimal.

Catmull-Rom spline consist of 4 basis functions with local support and are $C_{1}$ continuous and differentiable (see Suppl. Mat. Section 1). These properties makes them smoother compared to standard linear interpolation techniques. Qualitative comparison of Catmull-Rom spline with linear resampling technique is provided in Fig. 3 of Suppl. Mat.

Comparison of Catmull-Rom spline with other spline-based resamplers (including B-spline) is presented in Fig. 2. It can be observed that Catmull-Rom spline has the lowest training losses for both lung and brain datasets.

Continuous Warping Stage (CWS)

CWS (orange block in Fig. 1) warps the deformation field obtained at $l$ pyramid levels in a continuous fashion. Here $l=4$ for our 2D case with image size $256\times 256$ and $l=2$ for our 3D case with patch volume size $64\times 64\times 64$. Note that $l$ pyramids are constructed prior to feeding images into the LD-ConvNet.

The CWS employs warping of the next level source image $I^{s}_{l-1}$ with the computed DVF from previous coarse level LD-ConvNet, i.e., $\textbf{u}_{l-1}=\textbf{u}^{d\uparrow}_{l}$. This process is repeated until the final level, i.e., $l=0$. At each level the losses $\mathit{D}\left({I}^{s}_{l},{I}^{t}_{l}\right)$ between the warped image $I^{s}_{l-1}\left(\textbf{x}+\textbf{u}_{l-1}\right)$ and the target are summed and transferred to the final constrained loss function detailed in the Section 2.

Loss function

A multi-modal normalized cross-correlation (NCC) metric is used as data term $\mathit{D}$ in our loss function. We propose to use a derived sum-of-squared difference (SSD) which guarantees maximization of NCC metric [4] and is written as $\mathit{D}\left({I}^{s},{I}^{t}\right)=\tfrac{1}{2N}\sum\left(\tfrac{\mid I^{s}(\textbf{x})-\mu^{s}\mid}{{\sqrt{\sigma_{s}^{2}+\epsilon^{2}}}}-\tfrac{\mid I^{t}(\textbf{x})-\mu^{t}\mid}{\sqrt{\sigma_{t}^{2}+\epsilon^{2}}}\right)^{2},$ with mean $\mu$, standard deviation $\sigma$, the total number of pixels $N$ and $\epsilon=10^{-3}$ (used to avoid division by zero). For obtaining a smooth deformation vector field we apply an $L_{2}$-norm of the difference between the computed $\textbf{u}^{d}_{0}$ from the NL-DVF-R block and $\textbf{u}_{0}^{\uparrow}$ which is the previous upsampled input DVF to CWS block (see Fig. 1). Losses computed during $l$ pyramid levels in our CWS block $\sum_{l-1}^{1}\mathit{D}_{k}$ are aggregated to the final loss. The resulting equation for the backward propagation is: $\mathcal{L}=\mathit{D}\left({I}^{s},{I}^{t}\right)+\lambda\left[\parallel\textbf{u}_{0}^{d}-\textbf{u}_{0}^{\uparrow}\parallel^{2}\right]_{0}^{0.25}+\sum_{l-1}^{1}\mathit{D}_{k},$ where $\lambda$ is the trade-off between the regularization and the data-term. In our experiments it is set to 0.001. We restrict the regularization to the interval $[0,0.25]$ to prevent over smoothing problem during training. Omitting this constraints results in failures which are documented in Fig. 2 of Suppl. Mat. (see plot for Conv2Warp-FullReg). The loss saturates after 6000 iterations. However, regularised models with the interval continue to minimise the loss function. Our constrained loss is optimized using Adam optimizer with learning rate $10^{-3}$.

3 Experiments and Results

3.1 Datasets and training

4D CT. DIR-LAB [2] thoracic image dataset for DIR that inlcudes inspiratory and expiratory breath-hold CT image pairs were used for training. Although large permutations are possible, only breathing cycles 00-50, 10-80 and 30-90 were considered from 10 different sets (a total of 120 volumes). As a consequence we avoid an imbalance in the training dataset as most of the breathing cycles did not present strong deformations. For testing the POPI [9] lung dataset was used. The provided segmentation for air, body and lungs were used for evaluation.

MRI. The LBPA40 [8] dataset were used to train the model. We performed a random combination of these brain pairs and used nearly 300 volumes for our training. We have used MGH10 and CUMC12 [5] datasets for which masks are available for evaluation. 300 epochs with batch size of 4 were used for training performed on a 16GB NVIDIA Tesla P100 for nearly 48 hours.

3.2 Registration of 4D CT data

Conv2Warp has been evaluated by comparing it to the DIR state-of-the-methods. We also provide results where we restrict the convolution blocks in the Conv2Warp model to linear convolution only (ConvNet). For both these architectures, we will compare them only for Catmull-Rom spline based interpolation as it outperformed other spline techniques in our experiments. Seven pairs of different breathing cycles of 4D-CT data [9], each with 141 slices, where improvement of DIR methods compared to pre-alignment were significant are shown in Table 1. While the accuracy of our method measured by Dice (.90) and Jaccaard coefficient (.84) is similar to the rigorous simpleElastix (SE111http://simpleelastix.github.io), the run time is reduced by a factor of 141. Conv2Warp outperforms all other state-of-the art methods and ConvNet for almost all considered pairs. The test time on CPU is nearly 2.94 s which is multiple-folds lower than conventional DIR methods. Our model is light weight with inference time on a GPU of less than 1 s.

3.3 Registration of T1 MRI data

T1-weighted MRI scans for 9 volume pairs of MGH10 and 11 volume pairs of CUMC12 [5] where each image from different volumes were registered to the first volume data were also used for evaluation. The weights trained on LBPA [8] dataset were used in this case. In Table 2 Conv2Warp have $\mu_{dice}$ (mean dice) and $\mu_{jaccard}$ (mean Jaccard) of 0.95 and 0.90, respectively for MGH10 dataset (same as for SE), and 0.97 and 0.93, respectively for CUMC12 dataset (higher than SE). Conv2Warp has higher $\mu_{dice}$ and $\mu_{jaccard}$ for both datasets compared to ANTS (SyN-CC) and ConvNet. When compared with speed Conv2Warp is computationally the fastest among all other conventional DIR methods.

3.4 Visual validation

In Fig. 3 Conv2Warp without ($\lambda=0$) regularization and with ($\lambda=0.001$) regularization are shown in 4th and 5th columns respectively and their corresponding DVF magnitudes in 6th and 7th columns, respectively. It can be observed that the unconstrained loss function results in some unrealistic deformations (red rectangular regions in 4th column) while a more realistic deformations are visible for the constrained loss proposed in Conv2Warp (5th column). A smooth deformation can be seen in the magnitude image of the DVF $\mid\textbf{u}_{reg}\mid$ (7th column). Colour overlay images show a large improvement in the alignment of source images $I^{s}$ with the target images $I^{t}$ (3rd column) with Conv2Warp method (5th column). Fig. 4 shows the results on T1 MRI test datasets which were first rigidly aligned to MNIspace and then co-registered using Conv2Warp. It can be observed that Conv2Warp handles different magnitudes of non-linear deformations.

4 Conclusion

We have proposed a novel end-to-end convolutional neural network that consists of a sequential linear and deformable convolutions along with a learnt non-linear sampler. To handle wide range of non-linear deformations between source and target data pairs deformations are concurred by the continuous warping strategy. Our experimental results demonstrate that our proposed model outperforms most of the traditional methods and has very low computational complexity. When compared to the existing linear deep learning models such as ConvNets, Conv2Warp produces more accurate deformation fields. Additional details and evaluations are available in the supplementary material.

Acknowledgments

SA is supported by the NIHR Oxford BRC and JR by EPSRC EP/M013774/1 Seebibyte.