\section{Introduction}
Medical image registration 
consists of finding a spatial transformation that aligns two or more images. 
The intrinsic ill-posedness of registration can lead to anatomically implausible transformations
associated with unrealistic volumetric distortion of the anatomy.
For certain anatomical regions, such as the myocardium, physiologically plausible image registration requires volumetric preservation which corresponds to so-called \emph{incompressible} registration~\cite{Bistoquet2008,Mansi2011,Tobon-Gomez2013}.



Incompressible medical image registration between a pair of images $I_1, I_2$ 
can be defined as a constrained optimisation problem:
\begin{equation}
    \label{eq:reg_incompressible_continuous}
   
        \argmin_{\Phi \in \mathcal{D}\left(\Omega\right)}\,\,  
        \mathcal{L}(I_1, I_2, \Phi) + R(\Phi) \quad
        \text{s.t.} \quad
        \left[
        \forall x \in \mathcal{M},\quad \det\left(J_{\Phi}(x)\right) = 1
        \right]
   
\end{equation}
where $\Omega$ is the image domain, $\mathcal{M} \subset \Omega$ is the incompressible region, $\mathcal{L}$ an image dissimilarity measure, $R$ a regularisation term, $\mathcal{D}\left(\Omega\right)$ the group of diffeomorphic transformations from $\Omega$ onto itself, and $J_{\Phi}$ the Jacobian matrix of the transformation $\Phi$.

In practice, to solve \eqref{eq:reg_incompressible_continuous}, the transformation must be parameterised by a finite number of parameters, and the constraint $\forall x \in \mathcal{M},\, \det\left(J_{\Phi}(x)\right) = 1$ must be reduced to a finite number of equality constraints.
There are two approaches to discretise the constraint:
1) relaxing the constraint into an additive soft constraint~\cite{Aganj2015,DeCraene2012,Heyde2016,Rohlfing2003}.
2) using a specific parameterisation for the transformation~\cite{Bistoquet2008,Mang2018,Mansi2011}.
Soft constraints introduce hyperparameters that are difficult to tune reliably and cannot guarantee an incompressible transformation.
In \cite{Bistoquet2008}, the transformation is parameterised by a divergence-free displacement. Yet, this approach only provides a first order approximation of an incompressible deformation, still requiring a soft constraint.
Recently \cite{Mang2018,Mansi2011} proposed 
to parameterise the transformation by a divergence-free stationary velocity field (SVF).
The transformation is obtained via the Lie exponential mapping $\exp: v \mapsto \Phi $ that maps any divergence-free SVF $v$ into an incompressible transformation $\Phi$~\cite{Mang2018,Mansi2011}.
This reduces the non-convex constraint in \eqref{eq:reg_incompressible_continuous} to a linear constraint in the continuous domain.
However, in~\cite{Mang2018,Mansi2011} the SVF is parameterised as linear B-splines, and the linear constraint is discretised by imposing the constraint only on the points of the deformation grid.
As a result, guarantees to obtain a continuous incompressible transformation do not hold anymore.
To mitigate this issue,~\cite{Mang2018} proposed to work with images of higher resolution with a finer grid for the linear B-spline leading to the need for distributed super-computing.
Moreover,~\cite{Mang2018,Mansi2011} methods are limited to using sum of squared differences (SSD) as image similarity metric, which limits their applicability to images with similar intensity distributions.

In this paper, we propose a constrained optimisation framework for incompressible diffeomorphic registration that allows to use any smooth image similarity and regularisation penalty.
As an efficient means of discrete SVF parameterisation for this problem, we introduce multivariate divergence-conforming B-splines that have recently raised interest in computational physics~\cite{Evans2013a}. 
We demonstrate that their properties can be exploited to impose bounds on the divergence of the SVF over the entire continuous space using sparse linear constraints on its finite parameters.
Our general problem formulation allows us to solve incompressible registration using any state-of-the-art optimiser for constrained non-convex optimisation (e.g. \texttt{IPOPT}~\cite{Wachter2006}).
For evaluation, we initially apply our method for multi-modal incompressible registration of synthetically deformed brains.
We then compare our approach against the state-of-the-art results published for the STACOM'11 myocardial tracking challenge dataset~\cite{Tobon-Gomez2013} and achieve similar results while better retaining the incompressibility.

\newpage
\section{Method}
\subsubsection{Divergence-conforming B-splines}
General B-splines are very popular to parameterise deformations and velocity fields
over a continuous spatial domain $\Omega$
with a finite number of parameters
$\phi_i \in \mathds{R}^{3}$:
\begin{equation}
    \label{eq:bspline3d}
     \forall \left(x, y, z\right) \in \Omega,\quad
     v\left(x,y,z\right) = \sum_{i} B^k_{i,X}(x) B^k_{i,Y}(y) B^k_{i,Z}(z)\,\phi_{i},
\end{equation}
where 
$(x_i, y_i, z_i)$ are the knots of a regular grid of spacing $\delta x \times \delta y \times \delta z$ and 
$B^k_{i,U}$ are the
1D B-spline basis functions of order $k$ in the direction $U\in \{X, Y, Z\}$ (see supplemental material for more details).
Popular choices of the order are $k=1$ (linear B-splines)~\cite{Mang2018,Mansi2011} and $k=3$ (cubic B-splines)~\cite{Tobon-Gomez2013}.
A fundamental property of 1D B-spline basis functions of the variable $u$, on a regular grid 
of spacing $\delta u$ with $n$ knots is that:
\begin{equation}
\label{eq:deriv_bsplines}
   
    \forall i,
    \quad \dv{B^k_{i,U}}{u} = \frac{B^{k-1}_{i,U} - B^{k-1}_{i+1,U}}{\delta u}
\end{equation}
This implies that the derivative of a 1D B-spline on a regular grid is also a 1D B-spline on the same grid, albeit of lower order.
However, this property does not extend to 3D B-splines using definition \eqref{eq:bspline3d}. 
Especially, the divergence of those 3D B-splines are not B-splines,
because of the mixed order appearing with first partial derivatives.
This limitation makes it more difficult to relate properties of the velocity field $v$ to the values of the parameters
$\phi_i$.

To overcome this limitation, we propose to use a 3D divergence-conforming B-splines~\cite{Evans2013a} to parameterise $v$.
The orders of the 3D B-splines basis of each component are chosen so that the divergence of $v$ is, in a continuous sense, exactly a 1D B-spline of order $k$. 
Using the same notations as in \eqref{eq:bspline3d}, and 
using $\phi_i = (\phi^X_i,\phi^Y_i,\phi^Z_i)$,
we have:
\begin{equation}
    \forall (x,y,z) \in \Omega,\quad 
    v(x,y,z) = 
    \left(
    \begin{array}{c}
         \sum_{i} B^{k+1}_{i,X}(x) B^k_{i,Y}(y) B^k_{i,Z}(z)\,\phi^X_{i}  \\
         \sum_{i} B^k_{i,X}(x) B^{k+1}_{i,Y}(y) B^k_{i,Z}(z)\,\phi^Y_{i} \\
         \sum_{i} B^k_{i,X}(x) B^k_{i,Y}(y) B^{k+1}_{i,Z}(z)\,\phi^Z_{i}
    \end{array}
    \right)
\end{equation}
Using \eqref{eq:deriv_bsplines}, we obtain the continuous divergence $\div{v}$ for all $(x,y,z) \in \Omega$:
\begin{equation}
\label{eq:div_divergence-conforming}
    \begin{split}
        \div{v}(x,y,z) =& \sum_{i} B^{k}_{i,X}(x) B^k_{i,Y}(y) B^k_{i,Z}(z)\, \psi_i \\
        \text{s.t.}\,\, \forall i,\quad \psi_i =& \,
    \frac{\phi^X_{i} - \phi^X_{i-1}}{\delta x} +
    \frac{\phi^Y_{i} - \phi^Y_{i-1}}{\delta y} +
    \frac{\phi^Z_{i} - \phi^Z_{i-1}}{\delta z}
    \end{split}
\end{equation}
As a consequence, the following lemma states that $\div{v}$ is uniformly bounded at any point of a continuous subregion $\mathcal{M} \subset \Omega$ provided that a finite number of linear constraints are satisfied by the coefficients
$(\phi^X_i,\phi^Y_i,\phi^Z_i)$.

\begin{lemma}
\label{lemma:1}
Let $k\geq 2$, and $\mathcal{M}$ be a non-empty subset of $\Omega$. Let $\epsilon \geq 0$
and let $\mathbf{J}_{\mathcal{M}}=\{i
\,\,|\,\, 
\left(
\supp B^{k}_{i,X} \times \supp B^{k}_{i,Y} \times \supp B^{k}_{i,Z}
\right)
\cap \mathcal{M} \neq \emptyset
\}$.
If
\[
    \forall i \in \mathbf{J}_{\mathcal{M}},\quad -\epsilon \leq 
    \frac{\phi^X_{i} - \phi^X_{i-1}}{\delta x} +
    \frac{\phi^Y_{i} - \phi^Y_{i-1}}{\delta y} +
    \frac{\phi^Z_{i} - \phi^Z_{i-1}}{\delta z}
    \leq \epsilon,
\]
then at any point $m \in \mathcal{M}$, it holds that
$\,\abs{\div{v}(m)} \leq \epsilon$.
\end{lemma}

\begin{proof}
The proof follows from \eqref{eq:div_divergence-conforming} and the fact that the value of a B-spline at any point $m$ is bounded by the values of the B-spline coefficients associated to the knots that are close to $m$ (see supplementary material for more details).
\end{proof}

\subsubsection{Optimisation formulation and implementation}
In this section, we formulate the optimisation problem for diffeomorphic registration with the proposed incompressibility constraint on a subregion $\mathcal{M}$.
Using Lemma~\ref{lemma:1} and previous notations, our (symmetric) optimisation formulation of \eqref{eq:reg_incompressible_continuous} is:
\begin{equation}
    \label{eq:optim_formulation}
    \begin{split}
        \argmin_{\Theta=\{(\phi^X_i,\phi^Y_i,\phi^Z_i)\}_{i}} \,\, & \mathcal{L}(I_1 \circ \widetilde{\exp}(v(\Theta)), I_2) + \mathcal{L}(I_1, I_2 \circ \widetilde{\exp}(-v(\Theta))) + R(v(\Theta))\\
        \text{s.t. } \quad &
        \forall i \in \mathbf{J}_{\mathcal{M}},\,\,
            \frac{\phi^X_{i} - \phi^X_{i-1}}{\delta x} +
            \frac{\phi^Y_{i} - \phi^Y_{i-1}}{\delta y} +
            \frac{\phi^Z_{i} - \phi^Z_{i-1}}{\delta z}
            = 0\\
    \end{split}
\end{equation}
where 
$\widetilde{\exp}$ is an approximation
of the Lie exponential.
Thus we obtain a constrained optimisation formulation that guarantees the SVF to be \emph{exactly} divergence-free over the entire continuous subregion $\mathcal{M}$ and that can be solved with efficient state-of-the-art optimisers.
Using state-of-the-art solvers like \texttt{IPOPT}, that uses a primal-dual interior-point filter line-search method, an approximated solution of \eqref{eq:optim_formulation} that satisfies the constraints up to machine precision can be obtained.



\subsubsection{Lie exponential approximation and incompressibility}
As the Lie exponential has to be approximated in practice, having a divergence-free velocity field does not strictly guarantee that the resulting deformation will be incompressible.
We now quantify this approximation in our framework with respect to the time step $\tau$ for the Euler method.
Let $v$ denote a (SVF) solution of \eqref{eq:optim_formulation} that fulfills the divergence-free constraint up to machine precision $\epsilon_{mach}$. For any point $m \in \mathcal{M}$, the Jacobian of the first step of the Euler method, $I + \tau v$, fulfills:
\begin{equation}
    \label{eq:jac_one_step_error}
    \det\left(J_{I + \tau v}(m)\right) = 1 + \tau \div{v}(m) + \mathcal{O}\left(\tau^2\right)
    \quad \text{where} \quad \abs{\div{v}(m)} \leq \epsilon_{mach}.
\end{equation}
Then, by composing $1/\tau$ times, if the point $m$ remains inside $\mathcal{M}$ during the Euler integration, the Jacobian of $\widetilde{\exp}(v)$ satisfies (see supplementary material):
\begin{equation}
    \label{eq:jac_error}
    \det\left(J_{\widetilde{exp}(v)}(m)\right) = 1 + \mathcal{O}\left(\tau + \epsilon_{mach}\right)
\end{equation}
This approximation is independent to the spatial spacing and only depends on the time integration step $\tau$.
This is, to the best of our knowledge, a new state-of-the-art approximation for an incompressible diffeomorphic deformation parameterised by an SVF.
However, it is worth noting that \eqref{eq:jac_error} is only guaranteed
if $m$ remains inside $\mathcal{M}$ during the Euler integration.
Thus, when $\mathcal{M}$ is not equal to the entire spatial domain, this approximation may not hold for large deformations. 



\section{Evaluation on synthetic data for incompressible multi-modal registration}



\begin{table}[tb]
	\centering
	\caption{\textbf{Validation on incompressible multi-modal registration.} 
	RMSE values between predicted and ground-truth transformations are reported in mm. MAE($\abs{\text{Jacobian} - 1}$) stands for the Mean Absolute Error between the Jacobian map of the predicted transformation and a map uniformly equal to $1$ (perfect incompressibility).
	For both metrics, we reported mean value (standard deviation) over 60 registrations.
	}
	\begin{tabularx}{\textwidth}{c *{4}{Y}}
		\toprule
		 \textbf{Error measures} & Affine & Cubic B-splines & Ours\\ 
		\midrule
		RMSE (transformation)      & 2,28 (0,68) & 0.96 (0.44) & \textbf{0.90 (0.45)}\\
		MAE($\abs{\text{Jacobian} - 1}$)  & 0.0062 (0.005) & 0.054 (0.02) & \textbf{0.00079 (0.0004)}\\

	
		\bottomrule
	\end{tabularx}
	\label{tab:multimod_results}
\end{table}

We start by evaluating our method on synthetic incompressible registration in the context of multi-modal MRI data.

\subsubsection{Data generation method} We used T1, T2 and PD brain images from the IXI dataset\footnote{\url{https://brain-development.org/ixi-dataset/}}.
We generated realistic ground-truth incompressible transformations in two steps.
First, 
we non-linearly registered a pair of T1 images coming from different patients using a classical cubic B-splines SVF $v_0$.
Second, 
we generated a quasi-divergence-free SVF $v$ by projecting 
$v_0$
on the space of quasi-divergence-free SVFs.
This corresponds to solving:
\begin{equation}
    \label{eq:data_generation}
   
        \argmin_{\Theta=\{\phi_i\}_i}
        \frac{1}{2}\norm{v(\Theta) - v_0}^2 \quad
        \text{s.t.}
        \quad \forall i,
        \quad \div{v}(x_i, y_i, z_i) = 0
   
\end{equation}
where $v$ and $v_0$ are parameterised with classical cubic B-splines as in \eqref{eq:bspline3d}. We note that any potential bias towards the space of divergence-conforming B-splines is avoided in this comparison. We used \texttt{IPOPT}~\cite{Wachter2006} to solve \eqref{eq:data_generation}.

\subsubsection{Evaluation}
We generated a ground-truth quasi-incompressible SVF for 10 patients using the previous generation procedure, taking as fixed images 10 other patients.
Then for a given subject and for any pair of imaging modalities $M_1, M_2 \in \{T1, T2, PD\}$, we warped the first image $I(M_1)$ using the inverse of the ground-truth transformation. The task consists in estimating $v_{GT}$ by registering $I(M_1) \circ \widetilde{\exp}(-v_{GT})$ to $I(M_2)$.
We compared the proposed method with divergence-con\-forming B-splines SVF under incompressibility constraint to a classic registration approach based on cubic B-splines SVF similar to the one used to generate the ground-truth.
We also performed an affine registration to illustrate that the ground-truth transformation is not just affine.


\subsubsection{Implementation details}
We used \texttt{NiftyReg}~\cite{Modat2010} for the diffeomorphic registrations using cubic B-splines SVF.
For both \texttt{NiftyReg} and our incompressible implementation we used the same hyperparameters (grid size 5mm, $3$ levels of pyramid) and objective function, consisting of Normalised Mutual Information (NMI) as a similarity measure (weight $0.95$) and a bending energy regularisation term (weight $0.05$). 
The optimiser differs, as we used \texttt{IPOPT} optimiser while a Conjugate Gradient approach is used in \texttt{NiftyReg}.

\subsubsection{Results} \texttt{NiftyReg} and our approach recovered the ground-truth SVF up to a RMSE of the order of the images resolution, as shown in Table.~\ref{tab:multimod_results}. In addition, our method recovered the incompresibility with a higher accuracy.



\begin{figure}[tb!]
    \centering
    \includegraphics[width=\textwidth]{fig1_miccai19_v3_crop.pdf}
    \caption{We aim to register SSFP images at different times of the cardiac cycle.
    During the registration, the myocardium should not compress or expand.
    To impose this constraint we provide a mask of it to our algorithm (dashed in yellow).
    The floating image is warped into $\Phi$(Floating SSFP) so that it matches the reference image, while keeping the deformation of the myocardium incompressible, as shown in the log-Jacobian map.
    Red and blue represent expansion and compression respectively.
    }
    \label{fig:logJacobian}
\end{figure}

\section{Evaluation for myocardial tracking}
We evaluate the proposed method on the STACOM'11 myocardial tracking challenge~\cite{Tobon-Gomez2013}.
In particular, this allows a direct comparison of our framework with iLogDemons\cite{Mansi2011} a state-of-the-art incompressible registration method.


\subsubsection{Data}
The STACOM'11 dataset\footnote{\url{http://stacom.cardiacatlas.org/motion-tracking-challenge/}} contains 4D cine Steady State Free Precession (SSFP) of a full cardiac cycle with $30$ time points for $15$ patients.
Those SSFP come with
$12$ manually tracked landmarks.
We used only 13 of the 15 patients available because we found a shift between the images and the landmarks for two of them which we reported to the organisers.
The coordinates of the landmarks obtained by competitive methods of the challenge are available alongside the original data. This allows a direct comparison with 
our method.

\subsubsection{Implementation details} 
We used a similar implementation to the one of iLogDemons for this challenge, where we registered the first frame to all subsequent frames for each patient and used a manually delineated mask for the first frame. Yet, we used Local Normalised Cross Correlation as a similarity measure.
We used a bending energy regularisation, and a grid size of $3$mm.

\begin{figure}[bt!]
    \centering
    \begin{subfigure}[t]{0.37\textwidth}
        \centering
        \includegraphics[width=\linewidth]{Figure_boxplot_v2.png}
    \end{subfigure}
    \begin{subfigure}[t]{0.62\textwidth}
        \centering
        \includegraphics[trim={0 0 0 0.8cm},clip, width=\linewidth]{Figure_myocardium_incompressibility_v2.png} 
    \end{subfigure}
    \caption{(left) Distances to the landmarks after registering SSFP data.
    Our method achieves similar results to iLogDemons~\cite{Mansi2011} for an incompressible registration.
    UPF~\cite{DeCraene2012} achieves the highest accuracy, but it does not guarantee incompressibility.
    (right) Evolution of the Jacobian values distribution in the incompressible region (myocardium) while registering the first frame to all other time frames. The Jacobian values are uniformly close to 1 for small and moderate deformations of the myocardium. For large deformations (frames at the opposite to the cardiac cycle) 
    the mean Jacobian value is close to 1, but the dispersion of the Jacobian value distribution is large.
    This can be attributed to the use of a SVF.
    }
    \label{fig:ssfp_boxplot}
\end{figure}

\subsubsection{Results}
The evaluation of the myocardial tracking is based on the manually annotated landmarks at End Diastole and End Systole.
Using a Wilcoxon signed-rank test, we found that our results, shown in Fig.~\ref{fig:ssfp_boxplot}, are not statistically different from the results of iLogDemons~\cite{Mansi2011} in terms of landmark tracking error.
The log-Jacobian map of Fig.~\ref{fig:logJacobian} illustrates our approximation results \eqref{eq:jac_error}: for a small deformation, our registration framework can guarantee an incompressible deformation in a subregion with high accuracy.
While Fig.~\ref{fig:ssfp_boxplot} illustrates the degradation of this result for larger deformation, i.e. when registering frames at the opposite of the cardiac cycle.

\newpage

\section{Conclusion and Future Work}

\subsubsection{Limitations}
Similar to~\cite{Mang2018,Mansi2011}, the divergence-free constraint is imposed on the stationary velocity field (SVF) in a subregion $\mathcal{M}$, rather than imposing incompressibility on $\mathcal{M}$.
A voxel that is originally in $\mathcal{M}$ might exit $\mathcal{M}$ during the integration of the SVF and start being transported by velocity vectors that are not divergence-free (see supplemental material for a mathematical justification).
In this case, the deformation of this voxel is no longer incompressible.
As a result, using an SVF for incompressible registration applies only for deformations that are small enough or when the whole spatial domain is constrained to be divergence-free.
Investigating other diffeomorphic parameterisations is left for future work. 


\subsubsection{Advantages compared to previous methods}
Our method for incompressible registration relies on the parametrisation of the velocity field by a divergence conforming B-splines~\cite{Evans2013a}. 
In contrast to classical B-splines used in~\cite{Mang2018,Mansi2011}, it guarantees that the divergence of the velocity field is still a B-spline. 
This parametrisation along with constrained optimisation methods allows us to impose the velocity field to be divergence-free up to machine precision ($10^{-16}$ in our experiments). This is irrespective of the grid resolution chosen for the velocity field.
As a result, the proposed method 
is scalable to 3D images with high resolution.

We also proved an error bound for the incompressibility of the deformation for the proposed method in the case an Euler integration of the velocity field is used \eqref{eq:jac_error}.
The study of the error bound for other integrators, that may require additional interpolations (e.g. scaling-and-squaring), is left for future works.
Additionally, previous (quasi-)incompressible registration methods~\cite{Mang2018,Mansi2011} have been limited to using SSD as image similarity metric. We are proposing the first incompressible non-linear registration method that supports any smooth image similarity measure and spatial regularization.


\subsubsection*{Acknowledgments}
This project has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Sk{\l}odowska-Curie grant agreement TRABIT No 765148;
Wellcome [203148/Z/16/Z; 203145Z/16/Z; WT101957], EPSRC [NS/A000049/1; NS/A000050/1; NS/A000027/1; EP/L016478/1].
TV is supported by a Medtronic / RAEng Research Chair [RCSRF1819\textbackslash7\textbackslash34].

\bibliographystyle{splncs04.bst}
