11institutetext:  Department of Computing Science, University of Alberta,
Canada 22institutetext: Radiology and Diagnostic Imaging, University of
Alberta, Canada 33institutetext: Servier Virtual Cardiac Centre, Mazankowski
Alberta Heart Institute, Canada
Corresponding author. E-mail<EMAIL_ADDRESS>
Contributing authors: mnoga@ualberta.ca,punithak<EMAIL_ADDRESS>
*These authors contributed equally to this work.

# A training-free recursive multiresolution framework for diffeomorphic
deformable image registration

Ameneh Sheikhjafari  113*3* Michelle Noga  2233 Kumaradevan Punithakumar
223**3** Nilanjan Ray  1**1**

###### Abstract

Diffeomorphic deformable image registration is one of the crucial tasks in
medical image analysis, which aims to find a unique transformation while
preserving the topology and invertibility of the transformation. Deep
convolutional neural networks (CNNs) have yielded well-suited approaches for
image registration by learning the transformation priors from a large dataset.
The improvement in the performance of these methods is related to their
ability to learn information from several sample medical images that are
difficult to obtain and bias the framework to the specific domain of data. In
this paper, we propose a novel diffeomorphic training-free approach; this is
built upon the principle of an ordinary differential equation.

Our formulation yields an Euler integration type recursive scheme to estimate
the changes of spatial transformations between the fixed and the moving image
pyramids at different resolutions. The proposed architecture is simple in
design. The moving image is warped successively at each resolution and finally
aligned to the fixed image; this procedure is recursive in a way that at each
resolution, a fully convolutional network (FCN) models a progressive change of
deformation for the current warped image. The entire system is end-to-end and
optimized for each pair of images from scratch. In comparison to learning-
based methods, the proposed method neither requires a dedicated training set
nor suffers from any training bias. We evaluate our method on three cardiac
image datasets. The evaluation results demonstrate that the proposed method
achieves state-of-the-art registration accuracy while maintaining desirable
diffeomorphic properties.

###### Keywords:

Diffeomorphic registration Multiresolution deformable image registration Euler
integration Convolutional neural network.

##  1 Introduction

Deformable image registration is an important part of a variety of medical
image processing applications. For instance, deformable registration is
applied for the quantification of regional function measures and tracking
heart motion [12, 24] in cardiology. The aim of the deformable image
registration is to establish a point-wise correspondence between two images to
align them by estimating an appropriate non-linear transformation.

There is a strong need for the development of accurate and reliable nonrigid
registration algorithms for many medical imaging applications. Many deformable
registration algorithms have been proposed by various research groups and
existing algorithms could be broadly categorized into non-learning-based
(classical) methods and learning-based methods.

Commonly, non-learning-based registration techniques solve regularized
optimization problems using an image similarity/dissimilarity metric and
regularization. These methods include nonlinear elastic-based methods [6],
diffusion-based methods [15], B-splines [17, 26], demons algorithm [16],
diffeomorphic transformation-based methods [1, 15, 20]. Adapting the
conventional methods to specific problems is challenging, as they require a
precise model parameterization that can be too flexible or too restrictive.
Besides, occasionally regularization terms need to be manually adjusted for
each application to obtain accurate results. However, recent studies show that
the regularization prior can be modeled by neural networks, independent of
learning, for deformable image registration [14, 25] to remove the need for
manual adjustment of parameters.

In the past few years, learning-based approaches, notably deep CNNs, have been
successfully applied in different computer vision problems such as color
enhancement, color compensation, and image registration [28, 29, 30, 31, 32].
In supervised-learning based image registration methods, CNN is trained using
examples of medical images along with their ground truth transformations to
predict the transformations directly on test images [7, 22, 33]. Even though
the accuracy of these approaches is considerable, their performance is bounded
by the quality of the synthetic ground truth deformation field or the
annotated map [23]. Additionally, the actual ground truth of a desired neural
network output is not often available.

Unsupervised image registration methods have received a lot of attention
because they bypass the need for annotated data. In addition, the quality of
the output in unsupervised registration methods does not depend on the quality
of labels. However, they are often trained on a set of images which could
introduce bias on specific data. These methods utilize a CNN and a
differentiable similarity function to learn the dense spatial mapping between
input image pairs [3, 28]. Also, it must be noted that most of the existing
CNN-based methods often do not consider desirable diffeomorphic properties
including invertibility and inverse consistency of the transformation [19].

In this paper, we propose a novel learning-free diffeomorphic recursive
framework using an FCN capable of estimating plausible non-linear
transformations between images. The main elements of this work are:

  1. 1.

In the proposed method, an FCN models the changes in the deformation over
multiple resolutions as opposed to the deformation itself that most existing
methods do. The final deformation is estimated by a solution to an ordinary
differential equation (ODE). Thus the resulting algorithm is recursive in
nature.

  2. 2.

Following this recursion, the moving image is warped successively, enabling
the final prediction, which might consist of large displacements, to be
decomposed into smaller displacements.

  3. 3.

Within the recursive algorithmic framework, we introduce diffeomorphism, which
guarantees the inverse consistency of deformations.

  4. 4.

The proposed method is learning-free, i.e., it optimizes the parameters of the
FCN from scratch for every new pair of images and eliminates the need for a
dedicated training set.

  5. 5.

The parameters of the FCN are shared across all resolutions and the neural net
loss function values at all resolutions are optimized simultaneously.

Input: Pixel coordinates $\\{X_{t}\\}_{t=1}^{K}$ for all resolutions

Output: Deformation fields $\\{D_{t}\\}_{t=1}^{K}$ for all resolutions

-Initialize a spatial transformation

$V_{K}=f_{\theta}(X_{K},0)$

$D_{K}=EXP(V_{K})$

for  _$t\leftarrow K-1$ to $0$_ do

$D^{up}_{t+1}=UpSample(D_{t+1})$

$V_{t}=f_{\theta}(X_{t},D^{up}_{t+1})$

$V_{t}\leftarrow LowPass(V_{t})$

$D_{t}=D^{up}_{t+1}\circ EXP(V_{t})$

end for

Algorithm 1 ODE-based Multi-resolution Diffeomorphism

##  2 Methodology

Let $I_{m}\in\Omega$, $I_{f}\in\Omega$ denote the moving image and the fixed
image respectively which are both defined over d-dimensional space $\Omega$.
$D$ is the deformation field which maps $D:\Omega\xrightarrow{}\Omega$.
Commonly, deformable registration methods construct a deformation prediction
function $F$ which takes moving and fixed images as inputs and predicts a
dense deformation field that aligns $I_{m}$ to $I_{f}$. In contrast, in the
proposed approach, we construct $F$ as a function of coordinate grid (not
images) to predict the the change in the deformation field from one resolution
to the next resolution. We cascade this procedure by recursively performing
registration on the multi-resolution levels. Following this recursion, the
change of deformation field is predicted at each resolution, enabling the
final deformation field probably with large displacement to be decomposed into
cascade with small displacement.

In multi-resolution pyramids, a Gaussian pyramid of images is constructed
where the original image lies at the bottom level and subsequently higher
levels have a down-scaled Gaussian blurred version of the image. Multi-
resolution pyramids often serve to accelerate the optimization and yield
better accuracy. Using a multi-resolution recipe, two image pyramids are
built: $I_{F}^{t}$ and $I_{M}^{t}$ for $t=K,\dots,1$, where K is the maximum
level in the pyramid. Here, $I_{F}^{1}=I_{F}$ and $I_{M}^{1}=I_{M}$ are the
original fixed and moving images, respectively.

Considering resolution $t$ as the continuous variable, the change in
deformation over resolution $t$ is conceived by a FCN $f_{\theta}$:

$\frac{d}{dt}D_{t}=f_{\theta}(X_{t},D_{t}),$ (1)

where $X_{t}$, a multi-dimensional array, is the pixel coordinate grid at
resolution $t$ and $D_{t}$ is the deformation field, another multi-dimensional
array at resolution $t$. Thus, we model the change in deformation as a
function of both the coordinate grid and the deformation. Discretization of
the resolution (i,e., pyramid levels) leads to a solution by a Euler
integration with an initial value $D_{K}$:

$D_{t}=D_{t+1}+f_{\theta}(X_{t+1},D_{t+1}),~{}~{}\text{for}~{}t=K-1,K-2,...,1,$
(2)

where $t=K$ is the coarsest and $t=1$ is the finest/original resolution.
$D_{t+1}$ is the upsampled deformation field from a lower level $D_{t}$ and
$f_{\theta}(D_{t+1})$ is the changes of the deformation field at resolution
$t+1$ from resolution ${t}$.

Additionally, to have the same canonical range of pixel coordinates at every
resolution, in our implementation, we use the range $[-1,1]\times[-1,1]$ for
pixel coordinates. With this view, a multi-resolution pyramid adds more
samples in the space $[-1,1]\times[-1,1]$ as we go from lower to higher
resolutions.

A reasonable deformation field should prevent folding and be continuously
varying (diffeomorphism). By taking the composite form in the Euler
integration instead of the additive form, we introduce diffeomprhism to the
framework. Further, if we treat the output of the FCN as the velocity field,
then exponential of the velocity field will yield the multi-resolution
diffeomorphic deformation:

$D_{t}=D_{t+1}\circ
EXP(f_{\theta}(X_{t+1},D_{t+1})),~{}~{}\text{for}~{}t=K-1,K-2,...,1.$ (3)

The computational scheme for (3) appears in Algorithm 1, where we introduced
two more functionalities: upsampling and smoothing. Upsampling interpolates a
deformation from a lower to the next upper level of resolution and smoothing
of the velocity field acts as an implicit regularizer for image registration.

In this work, we use geometric optimization that uses the local canonical
coordinates [27]. Based on this method, to estimate the current deformation, a
compositional update rule is used between the previous estimate and the
exponential map $EXP$ of the displacement field. The exponential map is
efficiently calculated by using the scaling and squaring method [10] and the
composition of displacement fields. The exponentiation of the displacement
field ensures the diffeomorphism of the mapping. Therefore, our recursive
diffeomorphic image registration is obtained with Algorithm 1. Also the
schematic of the proposed recursive algorithm is shown in 1.

Figure 1: Schematic of the proposed recursive algorithm.
$PC_{K-1},PC_{K-2},...,PC_{0}$ are pixel coordinates correspond to the image
pyramids. $PC$ and $-PC$ are passed as the inputs to the network for forward
and backward directions. An initial grid with zero displacement is used as a
correction for the resolution $t=K-1$. The outputs of the FCN, forward
deformation field $FWD$ and backward deformation field $BWD$ will be smooth by
a Gaussian filter and then passed to the scaling and squaring method $SSQ$.
The final $FWD$ and $BWD$ generated at each resolution will be passed to an
upsampler and then the upsampled results will be used as a correction for the
next resolution. For each pair of images we repeat the whole process (e.g.
$800$ iterations) till we find the optimum value of $\theta$

To optimize the FCN parameters $\theta$, we use two Gaussian pyramids, one
each for the moving $I_{M}^{t}$ and the fixed $I_{F}^{t}$ images, where
$t=1,...,K$. Then, the FCN’s parameters will be optimized based on the
difference between warped moving image and fixed image at the same resolution.
In our approach, the neural network will be optimized from the scratch for
each pair of images being registered, and therefore, the trainable parameter
values of the network will at convergence will be different. Our approach is a
classical learning-free setting which would not require a dedicated training
set with annotations of any kind from experts.

The parameters $\theta$ of the FCN are found by solving the following
optimization using gradient descent (we use PyTorch’s autograd):

$\begin{split}\min_{\theta}\sum_{t=1}^{K}\sum_{x\in\Omega}&L(I_{F}^{t}(x),I_{M}^{t}(D_{t}(x)))+\gamma\left\lVert
D_{t}(x)-x\right\rVert^{2},\end{split}$ (4)

where $\Omega$ denotes the fixed and moving images domain, $I_{F}^{t}(x)$ is
the fixed image pixel value at pixel location $x$ and resolution level $t.$
$I_{M}^{t}(D_{t}(x))$ in (4) refers to the pixel value on the moving image
$I_{M}^{t}$ at location $x$ after being displaced by deformation field
$D_{t}(x).$ $L(.,.)$ is a differentiable loss function. In order to use the
same range of displacement for all levels of resolution, we normalized the
deformation field and the index at all resolutions in a canonical range
$[-1,1]$.

ES

ED

Registered ES

Flow

$D^{F}$

ED

ES

Registered ED

Flow

$D^{B}$

Figure 2: A sample result of forward and backward registrations with velocity
fields and grid deformation fields from ACDC dataset.

In most of the learning-based deformable image registration approaches, the
inverse mapping is often ignored [19]. However, the accuracy of the nonrigid
registration could be improved by computing the mappings from moving to fixed
image and vice versa as well as exploiting the invertible property of
diffeomorphic transformations.

In the proposed registration framework, we employ a FCN with two dense
velocity fields to generate two deformation fields, namely, forward $D^{F}$
and backward $D^{B}$ deformation fields, that register images $I_{M}$ to
$I_{F}$ and $I_{F}$ to $I_{M}$, respectively, using Algorithm 1. We then
enforce the consistency constraint that these two transformations are inverse
mappings of each other by adding a regularization term to the loss function.
Our FCN architecture consists of 4 convolution layers, each contains twenty-
four $5\times 5$ filters with a stride of $1$, followed by a rectified linear
unit (ReLU) activation, except for the last layer.

Note that the output of the Algorithm 1 is $D_{t}=[D_{t}^{F},D_{t}^{B}]$ in
which $D_{t}^{F}$ and $D_{t}^{B}$ are forward and backward deformation fields,
respectively, see Fig. 2 as an example. With this view the following forward-
backward optimization is used:

$\begin{split}\min_{\theta}&\sum_{t=1}^{K}\sum_{x\in\Omega}[L_{1}(I_{F}^{t}(x),I_{M}^{t}(D_{t}^{F}(x)))+L_{1}(I_{M}^{t}(x),I_{F}^{t}(D_{t}^{B}(x)))+\\\
&\alpha L_{2}(D_{t}^{B}(D_{t}^{F}(x)),x)+\alpha
L_{2}(D_{t}^{F}(D_{t}^{B}(x)),x)+\\\ &\gamma\left\lVert
D_{t}^{F}(x)-x\right\rVert^{2}+\gamma\left\lVert
D_{t}^{B}(x)-x\right\rVert^{2}]\end{split}$ (5)

We use differentiable mutual information (MI) [4] and structural similarity
index metric (SSIM) [4] as similarity metric $L_{1}=MI+SSIM$ and mean squared
error (MSE) as penalizes term to enforce the inverse consistency $L_{2}=MSE$.

Table 1: Quantitative evaluation of the results for cardiac MRI registration
on the ACDC, SCD and LA datasets. The evaluation was performed in terms of
Dice (mean$\pm$ standard deviation) and HD (mean) and the average number of
pixels with non-positive Jacobian Determinant (lower is better). Values in
bold indicate the best performance.

Methods |  | ACDC |  |  | SCD |
---|---|---|---|---|---|---
|  Dice | HD | $\det(J_{\phi})<0$ |  |  Dice | HD | $\det(J_{\phi})<0$
Unregistered |  0.65 $\pm$ 0.20 | 17.76 | N/A |  |  0.62 $\pm$ 0.15 | N/A | N/A
RTT(L2L2) |  0.76 $\pm$ 0.18 | 7.19 | 0.10 |  |  0.70 $\pm$ 0.12 | 5.36 | 0.30
RTT(L2L1) |  0.78 $\pm$ 0.18 | 6.64 | 0.25 |  |  0.70 $\pm$ 0.12 | 5.36 | 0.30
Elastix |  0.79 $\pm$ 0.18 | 11.26 | 0.32 |  |  0.79 $\pm$ 0.08 | N/A | 0.37
LCC-D* |  0.79 $\pm$ N/A | 9.21 | N/A |  |  - | - | -
VM* |  0.79 $\pm$ N/A | 8.46 | N/A |  |  - | - | -
LVM (S1)* |  0.79 $\pm$ N/A | 7.58 | N/A |  |  - | - | -
SyN* |  0.80 $\pm$ N/A | 8.24 | N/A |  |  0.81 $\pm$ 0.16 | 8.9 | 0.02
MM |  0.81 $\pm$ 0.18 | 8.59 | 0 |  |  0.72 $\pm$ 0.12 | 3.48 | 0
LVM (S3)* |  0.81 $\pm$ N/A | 6.88 | N/A |  |  - | - | -
Demons |  - | - | - |  |  0.65 $\pm$ 0.18 | 15.46 | 0.42
DIRNet** |  - | - | - |  |  0.80 $\pm$ 0.08 | N/A | N/A
Ours |  0.92 $\pm$ 0.10 | 4.76 | 0 |  |  0.85 $\pm$ 0.15 | 5.06 | 0

Methods |  | LA |  |
---|---|---|---|---
| Dice | HD | $\det(J_{\phi})<0$ |
| 2ch | 3ch | 4ch | 2ch | 3ch | 4ch |  | 2ch | 3ch | 4ch
Unregistered | 0.79 $\pm$ 0.07 | 0.78 $\pm$ 0.08 | 0.78 $\pm$ 0.09 | 7.37 | 7.70 | 8.66 |  | N/A | N/A | N/A
Elastix | 0.81 $\pm$ 0.13 | 0.86 $\pm$ 0.08 | 0.80 $\pm$ 0.11 | 7.9 | 6.99 | 8.01 |  | 0.68 | 0.5 | 0.75
RTT | 0.82 $\pm$ 0.11 | 0.85 $\pm$ 0.10 | 0.80 $\pm$ 0.10 | 7.28 | 6.82 | 7.56 |  | 0.53 | 0.45 | 0.7
MM | 0.84 $\pm$ 0.06 | 0.83 $\pm$ 0.06 | 0.83 $\pm$ 0.08 | 6.58 | 6.48 | 6.77 |  | 0.01 | 0 | 0
Demons | 0.84 $\pm$ 0.08 | 0.85 $\pm$ 0.06 | 0.82 $\pm$ 0.10 | 7.41 | 7.33 | 7.84 |  | 0.025 | 0.042 | 0
SyN | 0.87 $\pm$ 0.06 | 0.86 $\pm$ 0.13 | 0.84 $\pm$ 0.11 | 6.92 | 7.52 | 7.51 |  | 0.033 | 0 | 0.01
Ours | 0.89 $\pm$ 0.11 | 0.92 $\pm$ 0.06 | 0.86 $\pm$ 0.11 | 6.26 | 5.7 | 7.44 |  | 0 | 0 | 0

  * ∗Results are are reported in [13]. ∗∗Results are reported in [28]. N/A refers to the value that are not reported in the original paper.

##  3 Experiments

####  3.0.1 Data:

The proposed framework is evaluated on three clinical cardiac MRI datasets: 1)
Automated Cardiac Diagnosis Challenge (ACDC) dataset [5] which contains MRI
scans including 100 stack of short-axis cardiac cine MRI sequences, each
consisting of 12 to 35 temporal frames. 2) SunnyBrook cardiac dataset (SCD)
[21] which contains 45 cardiac cine MRI scans each with 20 temporal frames. 3)
Left atrium (LA) from the 2-chamber (2ch), 3-chamber (3ch) and 4-chamber (4ch)
cine long-axis MRI sequences acquired retrospectively at our center that
consists of 120 pairs of images in each chamber group. In all three datasets,
we applied the proposed framework on end-diastole (ED) and end-systole (ES)
cardiac frames from the cine sequences.

####  3.0.2 Implementation:

The proposed method is implemented in Python programming language using
Pytorch module. In neural network design, the Adam optimization with $800$
update iterations and a learning rate of $5\times 10^{-4}$ are used for all
the three datasets. The level of pyramid $K$ is considered $2$ with
$\alpha=1/(K)$, $\gamma=\frac{\lambda\times 1}{K}$, $\lambda=5$. The neural
net framework is evaluated on an NVIDIA GeForce GTX 1080 Ti GPU.

####  3.0.3 Evaluation Metrics:

The Dice metric (DM), reliability, Hausdorff distance (HD) and the determinant
of Jacobian are used as the evaluation metrics to compare the performance of
the algorithms.

##### Dice Metric

The DM [9] is a well-known segmentation based metric to measure the similarity
(overlap) between two regions, warped moving and fixed image. The DM of two
regions A and B is formulated as:

$DM(A,B)=\frac{2|A\cap B|}{A+B}$ (6)

##### Reliability

The algorithms are examined by evaluating the reliability function [2] of the
obtained DMs using (7). The complementary cumulative distribution function is
defined for each $d\in[0,1]$ as the probability of obtaining $DM$ higher than
$d$ over the entire set. The $R(d)$ measures how reliable is the algorithm in
yielding accuracy $d$.

$\begin{split}R(d)&=P_{r}(Dice>d)=\frac{\textnormal{\\# Images segmented with
DM higher than d}}{\textnormal{total number of images}}\end{split}$ (7)

##### Hausdorff Distance

The HD [11] is another well-known metric which measures the maximum deviation
between two regions’ contours. The HD between two contours $(C_{A})$ and
$C_{B}$ is formulated as:

$\textnormal{HD}(C_{A},C_{B})=\max(\max_{i}(\min_{j}(d(p^{i}_{A},p^{j}_{B}))),\max_{j}(\min_{i}(d(p^{i}_{A},p^{j}_{B}))))$
(8)

where $p^{i}_{A},p^{j}_{B}$ denote the set of all the points in $C_{A}$ and
$C_{B}$ respectively. The term $d(\cdot)$ denotes the Euclidean distance.

##### Determinant of Jacobian

We quantify and analyze deformation regularity using the determinant of the
Jacobian, $\det(J_{\phi})$. Jacobian matrix captures the local behaviours of
the deformation field. If the value of $\det(J)$ equals to $1$, the volume
remains constant after the transformation, the value smaller than $1$ shows
the local volume shrinkage and greater than $1$ shows the local volume
expansion. The negative values imply that local folding and twists have
occurred, which is physically not realizable and mathematically not inevitable
[8].

Table 2: Quantitative cardiac MRI registration results on $50\%$ patients of
ACDC, SCD and LA datasets based on different loss functions. The evaluation
was performed in terms of average Dice, HD and the forward and backward
registration time (per second).

| SCD |  |
---|---|---|---
Methods | Dice | HD | Time(s)
MSE | 0.83 | 9.1 | 15
SSIM | 0.86 | 7.9 | 20
SSIM + MI | 0.89 | 6.01 | 26

| ACDC |  |
---|---|---|---
Methods | Dice | HD | Time(s)
MSE | 0.87 | 5.73 | 15
SSIM | 0.92 | 4.31 | 21
SSIM + MI | 0.93 | 4.31 | 35

|  |  |  |  |  |  |
---|---|---|---|---|---|---|---
|  |  |  |  LA |  |  |
Methods | Dice | HD | Time (s)
|  2ch |  3ch |  4ch |  2ch |  3ch |  4ch |
MSE |  0.84 |  0.86 |  0.84 |  6.83 |  6.88 |  7.84 | 15
SSIM |  0.91 |  0.93 |  0.89 |  6.26 |  6.11 |  6.42 | 22
SSIM + MI |  0.90 |  0.93 |  0.87 |  6.89 |  5.22 |  6.40 | 35

2ch

3ch

4ch

Left Atrium

ACDC

SCD

Figure 3: Reliability versus Dice metric for different algorithms. The figures omit results for methods with no public implementation available. Table 3: Quantitative cardiac MRI registration results on the ACDC, SCD and LA based on different number of resolution. The evaluation is performed in terms of Dice (mean$\pm$ standard deviation) and HD (mean). The 2ch, 3ch and 4ch stand for the 2, 3 and 4-chamber. Values in bold indicate the best performance. | Dataset |  | Dice | HD
---|---|---|---|---
One Resolution | ACDC |  | 0.85 $\pm$ 0.20 | 6.03
SCD |  | 0.85 $\pm$ 0.17 | 7.70
LA-2ch |  | 0.87 $\pm$ 0.08 | 7.6
LA-3ch |  | 0.91 $\pm$ 0.05 | 6.59
| LA- |  | 0.88 $\pm$ 0.10 | 6.63
Two Resolutions | ACDC |  | 0.87 $\pm$ 0.17 | 6.27
SCD |  | 0.86 $\pm$ 0.16 | 7.06
LA-2ch |  | 0.89 $\pm$ 0.09 | 7.59
LA-3ch |  | 0.93 $\pm$ 0.03 | 5.47
| LA-4ch |  | 0.87 $\pm$ 0.11 | 6.67
Three Resolutions | ACDC |  | 0.83 $\pm$ 0.22 | 6.70
SCD |  | 0.85 $\pm$ 0.18 | 7.22
LA-2ch |  | 0.86 $\pm$ 0.09 | 7.03
LA-3ch |  | 0.90 $\pm$ 0.06 | 6.22
| LA-4ch |  | 0.86 $\pm$ 0.11 | 6.42
Four Resolutions | ACDC |  | 0.74 $\pm$ 0.24 | 8.96
SCD |  | 0.85 $\pm$ 0.19 | 7.50
LA-2ch |  | 0.82 $\pm$ 0.10 | 7.06
LA-3ch |  | 0.86 $\pm$ 0.07 | 6.93
| LA-4ch |  | 0.81 $\pm$ 0.10 | 7.01

###  3.1 Evaluation Results and Discussions

The performance of the proposed framework is compared with nine state-of-the-
art learning-based and classic deformable registration algorithms including
SimpleElastix (Elastix) [17], Moving Mesh (MM) [20], Real-Time Image-based
Tracker (RTT) [34], Demons [18], LCC-Demons (LCC-D) [15], Symmetric
Normalization (SYN) [1], [8], [13] and DIRNet [28].

Table 1, shows the mean and standard deviations of Dice scores, HD and the
average number of pixels with non-positive $\det(J_{\phi})$ over all subjects
ACDC, SCD, and left atrium datasets. Compared to the conventional and
learning-based methods, our framework has the best performances in terms of
Dice scores. In Fig. 4, the highest accuracy, $50th$ percentile, and the
lowest accuracy in terms of Dice scores achieved by the proposed method are
shown visually, where the red, blue, violet, gray and orange contours depict
the correspondence segmentation results by the proposed method, RRT1, RRT2,
Elastix and MM, respectively and the green contours depict the ground truth.

To quantify and analyse deformation regularity, we computed the determinant of
the Jacobian $\det(J_{\phi})$ using SimpleITK [18]. No negative values of
$\det(J)$ are observed for the proposed method as can be seen in Table 1.

We also evaluated the performance of the proposed algorithm in terms of the
reliability function $R(d)$ [2]. The corresponding $R(d)$ are plot as a
function of Dice score in Fig. 3. Our algorithm led to a higher reliability
curve on all ACDC, SCD, and left atrium datasets.

ACDC

SGD

Left-2ch

Left-3ch

Left-4ch

Figure 4: The segmentation results with the highest (first column), 50th
percentile (2nd column), lowest (3rd column) Dice score values for (a) ACDC
(1st row), (b) SGD (2nd row), and left atrium (c) 2-chamber (3rd row), (d)
3-chamber (4th row), 4-chamber (5th row) views. The green represent the ground
truth, also the red, blue, violet, gray and orange contours represent the
boundary corresponds to the registration by proposed method, RRT1, RRT2,
Elastix and MM, respectively

(a) ACDC

(b) SCD

(c) LA

ES

ED

Registered ES

Registered ED

$D^{F}$

$D^{B}$

Figure 5: Examples of pair images ES and ED from (a) ACDC, (b) SCD and (c) LA.
Each row shows forward and backward registered images and the corresponding
grid deformations. The registrations are applied to the whole image, however,
the grids are cropped to have a better resolution

We examined different loss functions, namely, MSE, SSIM and SSIM + MI, in our
framework and reported the average time to register one pair of images
(Computing both forward and backward deformation fields) as well as the
average of DM and HD over $50\%$ of patient datasets in Table 2. It can be
seen that the MSE loss function led to the best average time; however, it does
not yield the best accuracy in terms of the Dice score.

We assessed the impact of using different multi-resolution structures on the
performance, and reported the corresponding evaluations in terms of DM and HD
in Table 3. The results indicate that the high performance is obtained with
the use of two resolutions.

Examples of registered images from ACDC, SCD, and left atrium datasets with
their corresponding grid deformation are shown in Fig. 5 (a), (b), and (c),
respectively. Fig. 5 shows the original ED, ES and their correspondence
registered images and grid forward and backward deformations over images. The
algorithm is applied on whole images; however, to display the deformation
tracking on the part that has the most changes, we cropped the grid.

We have also evaluated the outputs from the forward deformation field by
ignoring the bidirectional registration to present its effects on the Dice
score and HD. The forward deformation field yields Dice score values of
$0.88\pm 0.17$ and $0.87\pm 0.14$ and HD values of $6.08$ and $5.19$ for ACDC
and SCD, respectively, while the proposed bidirectional approach yields Dice
score values of $0.92\pm 0.10$ and $0.90\pm 0.15$ and HD values of $4.76$ and
$5.06$ for ACDC and SCD, as reported in Table 1. This clearly indicates the
importance of including the bidirectional registration to obtain high
performance.

##  4 Conclusion

We present an ODE-based diffeomorphic recursive framework for multi-resolution
deformable registration using a FCN, to estimate the change velocity of
forward and backward deformation fields. We have then utilized inverse
consistency loss to further guarantee the desirable diffeomorphic properties
of the resulting solutions. The proposed framework is learning-free and does
not require a dedicated training set.

The proposed approach simultaneously estimates the forward and backward
mapping at all levels of the multi-resolution pyramid. We evaluated and
compared the method using three different MRI datasets against several state-
of-the-art traditional and learning-based registration methods. The results
demonstrate that our method outperforms both traditional and learning-based
methods.

## References

  * [1] Avants, B.B., Epstein, C.L., Grossman, M., Gee, J.C.: Symmetric diffeomorphic image registration with cross-correlation: evaluating automated labeling of elderly and neurodegenerative brain. Medical Image Analysis 12(1), 26–41 (2008)
  * [2] Ayed, I.B., Punithakumar, K., Li, S., Islam, A., Chong, J.: Left ventricle segmentation via graph cut distribution matching. In: International conference on medical image computing and computer-assisted intervention. pp. 901–909. Springer (2009)
  * [3] Balakrishnan, G., Zhao, A., Sabuncu, M.R., Guttag, J., Dalca, A.V.: An unsupervised learning model for deformable medical image registration. In: Proceedings of the IEEE conference on computer vision and pattern recognition. pp. 9252–9260 (2018)
  * [4] Belghazi, M.I., Baratin, A., Rajeswar, S., Ozair, S., Bengio, Y., Courville, A., Hjelm, R.D.: Mine: mutual information neural estimation. arXiv preprint arXiv:1801.04062 (2018)
  * [5] Bernard, O., Lalande, A., Zotti, C., Cervenansky, F., Yang, X., Heng, P.A., Cetin, I., Lekadir, K., Camara, O., Ballester, M.A.G., et al.: Deep learning techniques for automatic MRI cardiac multi-structures segmentation and diagnosis: Is the problem solved? IEEE transactions on medical imaging 37(11), 2514–2525 (2018)
  * [6] Burger, M., Modersitzki, J., Ruthotto, L.: A hyperelastic regularization energy for image registration. SIAM Journal on Scientific Computing 35(1), B132–B148 (2013)
  * [7] Cao, X., Yang, J., Zhang, J., Wang, Q., Yap, P.T., Shen, D.: Deformable image registration using a cue-aware deep regression network. IEEE Transactions on Biomedical Engineering 65(9), 1900–1911 (2018)
  * [8] Dalca, A.V., Balakrishnan, G., Guttag, J., Sabuncu, M.R.: Unsupervised learning for fast probabilistic diffeomorphic registration. In: International Conference on Medical Image Computing and Computer-Assisted Intervention. pp. 729–738. Springer (2018)
  * [9] Dice, L.R.: Measures of the amount of ecologic association between species. Ecology 26(3), 297–302 (1945)
  * [10] Higham, N.J.: The scaling and squaring method for the matrix exponential revisited. SIAM review 51(4), 747–764 (2009)
  * [11] Huttenlocher, D.P., Klanderman, G.A., Rucklidge, W.J.: Comparing images using the hausdorff distance. IEEE Transactions on pattern analysis and machine intelligence 15(9), 850–863 (1993)
  * [12] Khalil, A., Ng, S.C., Liew, Y.M., Lai, K.W.: An overview on image registration techniques for cardiac diagnosis and treatment. Cardiology research and practice 2018 (2018)
  * [13] Krebs, J., e Delingette, H., Mailhé, B., Ayache, N., Mansi, T.: Learning a probabilistic model for diffeomorphic registration. IEEE transactions on medical imaging (2019)
  * [14] Laves, M.H., Ihler, S., Ortmaier, T.: Deformable medical image registration using a randomly-initialized cnn as regularization prior. arXiv preprint arXiv:1908.00788 (2019)
  * [15] Lorenzi, M., Ayache, N., Frisoni, G.B., Pennec, X., (ADNI, A.D.N.I., et al.: LCC-Demons: a robust and accurate symmetric diffeomorphic registration algorithm. NeuroImage 81, 470–483 (2013)
  * [16] Mansi, T., Pennec, X., Sermesant, M., Delingette, H., Ayache, N.: iLogDemons: A demons-based registration algorithm for tracking incompressible elastic biological tissues. International journal of computer vision 92(1), 92–111 (2011)
  * [17] Marstal, K., Berendsen, F., Staring, M., Klein, S.: Simpleelastix: A user-friendly, multi-lingual library for medical image registration. In: Proceedings of the IEEE conference on computer vision and pattern recognition workshops. pp. 134–142 (2016)
  * [18] McCormick, M.M., Liu, X., Ibanez, L., Jomier, J., Marion, C.: ITK: enabling reproducible research and open science. Frontiers in neuroinformatics 8, 13 (2014)
  * [19] Mok, T.C., Chung, A.: Fast symmetric diffeomorphic image registration with convolutional neural networks. In: Proceedings of the IEEE/CVF conference on computer vision and pattern recognition. pp. 4644–4653 (2020)
  * [20] Punithakumar, K., Boulanger, P., Noga, M.: A GPU-accelerated deformable image registration algorithm with applications to right ventricular segmentation. IEEE Access 5, 20374–20382 (2017)
  * [21] Radau, P., Lu, Y., Connelly, K., Paul, G., Dick, A., Wright, G.: Evaluation framework for algorithms segmenting short axis cardiacMRI. The MIDAS Journal-Cardiac MR Left Ventricle Segmentation Challenge 49 (2009)
  * [22] Rohé, M.M., Datar, M., Heimann, T., Sermesant, M., Pennec, X.: SVF-Net: Learning deformable image registration using shape matching. In: International Conference on Medical Image Computing and Computer-Assisted Intervention. pp. 266–274. Springer (2017)
  * [23] Sang, Y., Xing, X., Wu, Y., Ruan, D.: Imposing implicit feasibility constraints on deformable image registration using a statistical generative model. In: Medical Imaging 2020: Image Processing. vol. 11313, p. 113132V. International Society for Optics and Photonics (2020)
  * [24] Sheikhjafari, A., Talebi, H.A., Zareinejad, M.: Robust and efficient 3d motion tracking in robotic assisted beating heart surgery. In: 2015 IEEE International Conference on Robotics and Biomimetics (ROBIO). pp. 1828–1833. IEEE (2015)
  * [25] Sheikhjafari, A., Noga, M., Punithakumar, K., Ray, N.: Unsupervised deformable image registration with fully connected generative neural network. In: Medical Imaging with Deep Learning (2018)
  * [26] Sheikhjafari, A., Talebi, H., Zareinejad, M.: 3D visual stabilization for robotic-assisted beating heart surgery using a thin-plate spline deformable model. In: 2015 3rd RSI International Conference on Robotics and Mechatronics (ICROM). pp. 743–748. IEEE (2015)
  * [27] Vercauteren, T., Pennec, X., Perchant, A., Ayache, N.: Non-parametric diffeomorphic image registration with the demons algorithm. In: International Conference on Medical Image Computing and Computer-Assisted Intervention. pp. 319–326. Springer (2007)
  * [28] de Vos, B.D., Berendsen, F.F., Viergever, M.A., Staring, M., Išgum, I.: End-to-end unsupervised deformable image registration with a convolutional neural network. In: Deep Learning in Medical Image Analysis and Multimodal Learning for Clinical Decision Support, pp. 204–212. Springer (2017)
  * [29] Xu, Y., Sun, B.: Color-compensated multi-scale exposure fusion based on physical features. Optik 223, 165494 (2020)
  * [30] Xu, Y., Sun, B., Yan, X., Hu, J., Chen, M.: Multi-focus image fusion using learning based matting with sum of the gaussian-based modified laplacian. Digital Signal Processing 106, 102821 (2020)
  * [31] Xu, Y., Yang, C., Sun, B., Yan, X., Chen, M.: A novel multi-scale fusion framework for detail-preserving low-light image enhancement. Information Sciences 548, 378–397 (2021)
  * [32] Yan, X., Liu, Y., Xu, Y., Jia, M.: Multistep forecasting for diurnal wind speed based on hybrid deep learning model with improved singular spectrum decomposition. Energy Conversion and Management 225, 113456 (2020)
  * [33] Yang, X., Kwitt, R., Styner, M., Niethammer, M.: Quicksilver: Fast predictive image registration–a deep learning approach. NeuroImage 158, 378–396 (2017)
  * [34] Zachiu, C., Papadakis, N., Ries, M., Moonen, C., de Senneville, B.D.: An improved optical flow tracking technique for real-time mr-guided beam therapies in moving organs. Physics in Medicine & Biology 60(23), 9003 (2015)