Incompressible image registration using divergence-conforming B-splines

TL;DR

This paper introduces a novel image registration method using divergence-conforming B-splines to ensure exact volumetric preservation, improving accuracy and flexibility over previous approaches.

Contribution

The paper presents a new framework for incompressible image registration with divergence-free velocity fields using divergence-conforming B-splines, enabling exact volume preservation.

Findings

Achieves exact divergence-free velocity fields at any point in the domain.

Demonstrates improved accuracy and volume preservation compared to state-of-the-art methods.

Provides theoretical insights into numerical incompressibility error.

Abstract

Anatomically plausible image registration often requires volumetric preservation. Previous approaches to incompressible image registration have exploited relaxed constraints, ad hoc optimisation methods or practically intractable computational schemes. Divergence-free velocity fields have been used to achieve incompressibility in the continuous domain, although, after discretisation, no guarantees have been provided. In this paper, we introduce stationary velocity fields (SVFs) parameterised by divergence-conforming B-splines in the context of image registration. We demonstrate that sparse linear constraints on the parameters of such divergence-conforming B-Splines SVFs lead to being exactly divergence-free at any point of the continuous spatial domain. In contrast to previous approaches, our framework can easily take advantage of modern solvers for constrained optimisation, symmetric registration approaches, arbitrary image similarity and additional regularisation terms. We study the numerical incompressibility error for the transformation in the case of an Euler integration, which gives theoretical insights on the improved accuracy error over previous methods. We evaluate the proposed framework using synthetically deformed multimodal brain images, and the STACOM11 myocardial tracking challenge. Accuracy measurements demonstrate that our method compares favourably with state-of-the-art methods whilst achieving volume preservation.

Tables

Table 1. Table 1 : Validation on incompressible multi-modal registration. RMSE values between predicted and ground-truth transformations are reported in mm. MAE( | Jacobian − 1 | Jacobian 1 \absolutevalue{\text{Jacobian}-1} ) stands for the Mean Absolute Error between the Jacobian map of the predicted transformation and a map uniformly equal to 1 1 1 (perfect incompressibility). For both metrics, we reported mean value (standard deviation) over 60 registrations.

Error measures	Affine	Cubic B-splines	Ours
RMSE (transformation)	2,28 (0,68)	0.96 (0.44)	0.90 (0.45)
MAE($\left|\text{Jacobian}-1\right|$)	0.0062 (0.005)	0.054 (0.02)	0.00079 (0.0004)

Equations

$$\underset{\Phi \in \mathcal{D}\left(\Omega \right)}{\mathrm{argmin}}\mathcal{L}(I_1,I_2,\Phi )+R(\Phi )\text{s.t.}\left[\forall x\in \mathcal{M},\det \left(J_{\Phi }(x)\right)=1\right]$$

$$\forall \left(x,y,z\right)\in \Omega ,v\left(x,y,z\right)=\sum _i{B}_{i,X}^k(x)B_{i,Y}^k(y)B_{i,Z}^k(z){\varphi }_i,$$

$$\forall i,\text{derivative}{B}_{i,U}^{k}u=\frac{B_{i,U}^{k-1}-B_{i+1,U}^{k-1}}{\delta u}$$

$$\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)$$

$$\begin{array}{rl}\text{divergence}v(x,y,z)=& \sum _i{B}_{i,X}^k(x)B_{i,Y}^k(y)B_{i,Z}^k(z){\psi }_i\\ \text{s.t.}\forall i,{\psi }_i=& \frac{{\varphi }_i^X-{\varphi }_{i-1}^X}{\delta x}+\frac{{\varphi }_i^Y-{\varphi }_{i-1}^Y}{\delta y}+\frac{{\varphi }_i^Z-{\varphi }_{i-1}^Z}{\delta z}\end{array}$$

$$\forall i\in {\mathbf{J}}_{\mathcal{M}},-ϵ\le \frac{{\varphi }_i^X-{\varphi }_{i-1}^X}{\delta x}+\frac{{\varphi }_i^Y-{\varphi }_{i-1}^Y}{\delta y}+\frac{{\varphi }_i^Z-{\varphi }_{i-1}^Z}{\delta z}\le ϵ,$$

$$\begin{array}{rl}\underset{\Theta =\{({\varphi }_i^X,{\varphi }_i^Y,{\varphi }_i^Z){\}}_i}{\mathrm{argmin}}& \mathcal{L}(I_1\circ \tilde{\exp }(v(\Theta )),I_2)+\mathcal{L}(I_1,I_2\circ \tilde{\exp }(-v(\Theta )))+R(v(\Theta ))\\ \text{s.t. }& \forall i\in {\mathbf{J}}_{\mathcal{M}},\frac{{\varphi }_i^X-{\varphi }_{i-1}^X}{\delta x}+\frac{{\varphi }_i^Y-{\varphi }_{i-1}^Y}{\delta y}+\frac{{\varphi }_i^Z-{\varphi }_{i-1}^Z}{\delta z}=0\end{array}$$

$$\det \left(J_{I+\tau v}(m)\right)=1+\tau \text{divergence}v(m)+\mathcal{O}\left({\tau }^2\right)\text{where}\text{absolutevalue}\text{divergence}v(m)\le {ϵ}_{mach}.$$

$$\det \left(J_{\widetilde{exp}(v)}(m)\right)=1+\mathcal{O}\left(\tau +{ϵ}_{mach}\right)$$

$$\underset{\Theta =\{{\varphi }_i{\}}_i}{\mathrm{argmin}}\frac{1}{2}\text{norm}{v(\Theta )-v_0}^2\text{s.t.}\forall i,\text{divergence}v(x_i,y_i,z_i)=0$$

$$\begin{array}{rl}\forall i=(i_X,i_Y,i_Z)\in \{-k,\dots ,n-1{\}}^3,& \forall \left(x,y,z\right)\in [0,1{]}^3,\\ B_{i_X,X}^k(x)=& B^k\left(\frac{x-x_{i_X}}{\delta x}-\frac{k+1}{2}\right)\\ B_{i_Y,Y}^k(y)=& B^k\left(\frac{y-y_{i_Y}}{\delta y}-\frac{k+1}{2}\right)\\ B_{i_Z,Z}^k(z)=& B^k\left(\frac{z-z_{i_Z}}{\delta z}-\frac{k+1}{2}\right)\end{array}$$

$$\begin{array}{rl}\forall t\in \text{mathds}R,& B^0\left(t\right)=1_{[-\frac{1}{2},\frac{1}{2}]}\left(t\right)\\ & B^1\left(t\right)=\left(1-|t|\right)1_{[-1,1]}\left(t\right)\\ & B^2\left(t\right)=\left(\frac{3}{4}-t^2\right)1_{[-\frac{1}{2},\frac{1}{2}]}\left(t\right)+\frac{1}{2}{\left(\frac{3}{2}-|t|\right)}^2{1}_{[\frac{1}{2},\frac{3}{2}]}\left(|t|\right)\\ & B^3\left(t\right)=\frac{1}{6}\left(4-3t^2\left(2-|t|\right)\right)1_{[0,1]}\left(|t|\right)+\frac{1}{6}{\left(2-|t|\right)}^3{1}_{[1,2]}\left(|t|\right)\end{array}$$

$$\begin{split}v(x,y,z)&=\left(\begin{array}[]{c}v^{X}(x,y,z)\\ v^{Y}(x,y,z)\\ v^{Z}(x,y,z)\end{array}\right)\\ &=\left(\begin{array}[]{c}\sum_{i_{X}=-k}^{n-1}\sum_{i_{Y},i_{Z}=-(k-1)}^{n-1}B^{k}_{i_{X},X}(x)B^{k-1}_{i_{Y},Y}(y)B^{k-1}_{i_{Z},Z}(z)\,\phi^{X}_{i_{X},i_{Y},i_{Z}}\\ \sum_{i_{Y}=-k}^{n-1}\sum_{i_{X},i_{Z}=-(k-1)}^{n-1}B^{k-1}_{i_{X},X}(x)B^{k}_{i_{Y},Y}(y)B^{k-1}_{i_{Z},Z}(z)\,\phi^{Y}_{i_{X},i_{Y},i_{Z}}\\ \sum_{i_{Z}=-k}^{n-1}\sum_{i_{X},i_{Y}=-(k-1)}^{n-1}B^{k-1}_{i_{X},X}(x)B^{k-1}_{i_{Y},Y}(y)B^{k}_{i_{Z},Z}(z)\,\phi^{Z}_{i_{X},i_{Y},i_{Z}}\end{array}\right)\end{split}$$

$$\forall k\ge 2,\forall t\in \text{mathds}R,\text{derivative}{B}^{k}t(t)=B^{k-1}\left(t+\frac{1}{2}\right)-B^{k-1}\left(t-\frac{1}{2}\right)$$

$$\begin{array}{rl}\text{derivative}{B}_{i,X}^{k}x& =\frac{B^{k-1}\left(\frac{x-x_{i_X}}{\delta x}-\frac{k+1}{2}+\frac{1}{2}\right)-B^{k-1}\left(\frac{x-x_{i_X}}{\delta x}-\frac{k+1}{2}-\frac{1}{2}\right)}{\delta x}\\ & =\frac{B^{k-1}\left(\frac{x-x_{i_X}}{\delta x}-\frac{k-1+1}{2}\right)-B^{k-1}\left(\frac{x-(x_{i_X}+\delta x)}{\delta x}-\frac{k-1+1}{2}\right)}{\delta x}\\ & =\frac{B_{i,X}^{k-1}-B_{(i_X+1,i_Y,i_Z),X}^{k-1}}{\delta x}\end{array}$$

$$\begin{array}{rl}\text{derivative}{B}_{i,Y}^{k}y=& \frac{B_{i,Y}^{k-1}-B_{(i_X,i_Y+1,i_Z),Y}^{k-1}}{\delta y}\\ \text{derivative}{B}_{i,Z}^{k}z=& \frac{B_{i,Z}^{k-1}-B_{(i_X,i_Y,i_Z+1),Z}^{k-1}}{\delta z}\end{array}$$

$$\begin{array}{r}\text{partialderivative}{v}^{X}x=\sum _{i_X=-k}^{n-1}\sum _{i_Y,i_Z=-(k-1)}^{n-1}\frac{B_{i_X,X}^{k-1}(x)-B_{i_X+1,X}^{k-1}(x)}{\delta x}{B}_{i_Y,Y}^{k-1}(y)B_{i_Z,Z}^{k-1}(z){\varphi }_{i_X,i_Y,i_Z}^X\end{array}$$

$$\begin{array}{rl}\text{partialderivative}{v}^{X}x=& \sum _{i_X=-k}^{n-1}\sum _{i_Y,i_Z=-(k-1)}^{n-1}{B}_{i_X,X}^{k-1}(x)B_{i_Y,Y}^{k-1}(y)B_{i_Z,Z}^{k-1}(z)\frac{{\varphi }_{i_X,i_Y,i_Z}^X}{\delta x}\\ & -\sum _{i_X=-(k-1)}^n\sum _{i_Y,i_Z=-(k-1)}^{n-1}{B}_{i_X,X}^{k-1}(x)B_{i_Y,Y}^{k-1}(y)B_{i_Z,Z}^{k-1}(z)\frac{{\varphi }_{i_X-1,i_Y,i_Z}^X}{\delta x}\end{array}$$

$$\begin{array}{rl}\text{partialderivative}{v}^{X}x=& \sum _{i_X,i_Y,i_Z=-(k-1)}^{n-1}{B}_{i_X,X}^{k-1}(x)B_{i_Y,Y}^{k-1}(y)B_{i_Z,Z}^{k-1}(z)\frac{{\varphi }_{i_X,i_Y,i_Z}^X-{\varphi }_{i_X-1,i_Y,i_Z}^X}{\delta x}\end{array}$$

$$\begin{array}{rl}\text{partialderivative}{v}^{Y}y=& \sum _{i_X,i_Y,i_Z=-(k-1)}^{n-1}{B}_{i_X,X}^{k-1}(x)B_{i_Y,Y}^{k-1}(y)B_{i_Z,Z}^{k-1}(z)\frac{{\varphi }_{i_X,i_Y,i_Z}^Y-{\varphi }_{i_X,i_Y-1,i_Z}^Y}{\delta y}\\ \text{partialderivative}{v}^{Z}z=& \sum _{i_X,i_Y,i_Z=-(k-1)}^{n-1}{B}_{i_X,X}^{k-1}(x)B_{i_Y,Y}^{k-1}(y)B_{i_Z,Z}^{k-1}(z)\frac{{\varphi }_{i_X,i_Y,i_Z}^Z-{\varphi }_{i_X,i_Y,i_Z-1}^Z}{\delta z}\end{array}$$

$$\begin{array}{rl}\text{divergence}v(x,y,z)=& \sum _{i_X,i_Y,i_Z=-(k-1)}^{n-1}{B}_{i_X,X}^{k-1}(x)B_{i_Y,Y}^{k-1}(y)B_{i_Z,Z}^{k-1}(z){\psi }_{i_X,i_Y,i_Z}\\ \text{s.t.}\forall i_X,i_Y,i_Z,{\psi }_{i_X,i_Y,i_Z}=& \frac{{\varphi }_{i_X,i_Y,i_Z}^X-{\varphi }_{i_X-1,i_Y,i_Z}^X}{\delta x}\\ & +\frac{{\varphi }_{i_X,i_Y,i_Z}^Y-{\varphi }_{i_X,i_Y-1,i_Z}^Y}{\delta y}\\ & +\frac{{\varphi }_{i_X,i_Y,i_Z}^Z-{\varphi }_{i_X,i_Y,i_Z-1}^Z}{\delta z}\end{array}$$

$$\begin{array}{rl}{\mathbf{J}}_{\mathcal{M}}=\{& (i_X,i_Y,i_Z)\in \{-(k-1),\dots ,n-1{\}}^3|\\ & \left(\mathrm{supp}{B}_{i_X,X}^{k-1}\times \mathrm{supp}{B}_{i_Y,Y}^{k-1}\times \mathrm{supp}{B}_{i_Z,Z}^{k-1}\right)\cap \mathcal{M}\ne \varnothing \}\end{array}$$

$$-ϵ\le \frac{{\varphi }_{i_X,i_Y,i_Z}^X-{\varphi }_{i_X-1,i_Y,i_Z}^X}{\delta x}+\frac{{\varphi }_{i_X,i_Y,i_Z}^Y-{\varphi }_{i_X,i_Y-1,i_Z}^Y}{\delta y}+\frac{{\varphi }_{i_X,i_Y,i_Z}^Z-{\varphi }_{i_X,i_Y,i_Z-1}^Z}{\delta z}\le ϵ,$$

$$\begin{array}{rl}\text{absolutevalue}\text{divergence}v(x,y,z)& =\text{absolutevalue}\sum _{i_X,i_Y,i_Z=-(k-1)}^{n-1}{B}_{i_X,X}^{k-1}(x)B_{i_Y,Y}^{k-1}(y)B_{i_Z,Z}^{k-1}(z){\psi }_{i_X,i_Y,i_Z}\\ & =\text{absolutevalue}\sum _{(i_X,i_Y,i_Z)\in {\mathbf{J}}_{\mathcal{M}}}{B}_{i_X,X}^{k-1}(x)B_{i_Y,Y}^{k-1}(y)B_{i_Z,Z}^{k-1}(z){\psi }_{i_X,i_Y,i_Z}\\ & \le \sum _{(i_X,i_Y,i_Z)\in {\mathbf{J}}_{\mathcal{M}}}{B}_{i_X,X}^{k-1}(x)B_{i_Y,Y}^{k-1}(y)B_{i_Z,Z}^{k-1}(z)\text{absolutevalue}{\psi }_{i_X,i_Y,i_Z}\\ & \le \sum _{(i_X,i_Y,i_Z)\in {\mathbf{J}}_{\mathcal{M}}}{B}_{i_X,X}^{k-1}(x)B_{i_Y,Y}^{k-1}(y)B_{i_Z,Z}^{k-1}(z)ϵ\\ & \le ϵ\end{array}$$

$$\begin{array}{rl}\tilde{\exp }& =\left(I+\frac{1}{2^K}v\right)\circ \dots \circ \left(I+\frac{1}{2^K}v\right)\\ & ={\left(I+\frac{1}{2^K}v\right)}^{2^K}\end{array}$$

$$\begin{split}\det\left(J_{I+\frac{1}{2^{K}}v}(m)\right)&=\left|\begin{array}[]{ccc}\frac{1}{2^{K}}\partialderivative{v^{X}}{x}+1&\frac{1}{2^{K}}\partialderivative{v^{X}}{y}&\frac{1}{2^{K}}\partialderivative{v^{X}}{z}\\ \frac{1}{2^{K}}\partialderivative{v^{Y}}{x}&\frac{1}{2^{K}}\partialderivative{v^{Y}}{y}+1&\frac{1}{2^{K}}\partialderivative{v^{Y}}{z}\\ \frac{1}{2^{K}}\partialderivative{v^{Z}}{x}&\frac{1}{2^{K}}\partialderivative{v^{Z}}{y}&\frac{1}{2^{K}}\partialderivative{v^{Z}}{z}+1\\ \end{array}\right|\\ &=1+\frac{1}{2^{K}}\divergence{v}(m)+\mathcal{O}\left(\left(\frac{1}{2^{K}}\right)^{2}\right)\end{split}$$

$$\{\begin{array}{rl}{m}_0& =m\\ \forall k\in \{1,\dots ,2^K\},m_k& =\left(I+\frac{1}{2^K}v\right)(m_{k-1})\\ & ={\left(I+\frac{1}{2^K}v\right)}^k(m)\end{array}$$

$$\begin{array}{rl}\det \left(J_{\widetilde{exp}(v)}(m)\right)& =\det \left(J_{{\left(I+\frac{1}{2^K}v\right)}^{2^K}}(m)\right)\\ & =\det \left(J_{\left(I+\frac{1}{2^K}v\right)\circ {\left(I+\frac{1}{2^K}v\right)}^{2^K-1}}(m)\right)\\ & =\det \left(J_{\left(I+\frac{1}{2^K}v\right)}(m_{2^K-1})\right)\det \left(J_{{\left(I+\frac{1}{2^K}v\right)}^{2^K-1}}(m)\right)\\ & =\prod _{k=0}^{2^K-1}\det \left(J_{I+\frac{1}{2^K}v}(m_k)\right)\\ & =\prod _{k=0}^{2^K-1}\left(1+\frac{1}{2^K}\text{divergence}v(m_k)+\mathcal{O}\left({\left(\frac{1}{2^K}\right)}^2\right)\right)\\ & =1+\frac{1}{2^K}\left(\sum _{k=0}^{2^K-1}\text{divergence}v(m_k)\right)+\mathcal{O}\left(\frac{1}{2^K}\right)\end{array}$$

$$\text{absolutevalue}\frac{1}{2^K}\left(\sum _{k=0}^{2^K-1}\text{divergence}v(m_k)\right)\le \frac{1}{2^K}\left(\sum _{k=0}^{2^K-1}\text{absolutevalue}\text{divergence}v(m_k)\right)\le {ϵ}_{mach}$$

$$\det \left(J_{\widetilde{exp}(v)}(m)\right)=1+\mathcal{O}\left(\tau +{ϵ}_{mach}\right)$$

Full text

1 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 incompressible registration [2, 7, 10].

Incompressible medical image registration between a pair of images $I_{1},I_{2}$ can be defined as a constrained optimisation problem:

[TABLE]

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 (1), 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 [1, 3, 5, 9].
2. using a specific parameterisation for the transformation [2, 6, 7]. Soft constraints introduce hyperparameters that are difficult to tune reliably and cannot guarantee an incompressible transformation. In [2], 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 [6, 7] 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$ [6, 7]. This reduces the non-convex constraint in (1) to a linear constraint in the continuous domain. However, in [6, 7] 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, [6] 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, [6, 7] 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 [4]. 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. IPOPT [11]). 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 [10] and achieve similar results while better retaining the incompressibility.

2 Method

2.0.1 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}$:

[TABLE]

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) [6, 7] and $k=3$ (cubic B-splines) [10]. 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:

[TABLE]

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 (2). 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 [4] 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 (2), and using $\phi_{i}=(\phi^{X}_{i},\phi^{Y}_{i},\phi^{Z}_{i})$, we have:

[TABLE]

Using (3), we obtain the continuous divergence $\divergence{v}$ for all $(x,y,z)\in\Omega$:

[TABLE]

As a consequence, the following lemma states that $\divergence{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})$.

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(\operatorname*{supp}B^{k}_{i,X}\times\operatorname*{supp}B^{k}_{i,Y}\times\operatorname*{supp}B^{k}_{i,Z}\right)\cap\mathcal{M}\neq\emptyset\}$. If

[TABLE]

then at any point $m\in\mathcal{M}$, it holds that $\,\absolutevalue{\divergence{v}(m)}\leq\epsilon$.

Proof

The proof follows from (5) 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).

2.0.2 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 1 and previous notations, our (symmetric) optimisation formulation of (1) is:

[TABLE]

where $\widetilde{\exp}$ is an approximation of the Lie exponential. Thus we obtain a constrained optimisation formulation that guarantees the SVF to be 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 IPOPT, that uses a primal-dual interior-point filter line-search method, an approximated solution of (6) that satisfies the constraints up to machine precision can be obtained.

2.0.3 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 (6) 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:

[TABLE]

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):

[TABLE]

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 (8) 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.

3 Evaluation on synthetic data for incompressible multi-modal registration

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

3.0.1 Data generation method

We used T1, T2 and PD brain images from the IXI dataset111https://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:

[TABLE]

where $v$ and $v_{0}$ are parameterised with classical cubic B-splines as in (2). We note that any potential bias towards the space of divergence-conforming B-splines is avoided in this comparison. We used IPOPT [11] to solve (9).

3.0.2 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-conforming 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.

3.0.3 Implementation details

We used NiftyReg [8] for the diffeomorphic registrations using cubic B-splines SVF. For both 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 IPOPT optimiser while a Conjugate Gradient approach is used in NiftyReg.

3.0.4 Results

NiftyReg and our approach recovered the ground-truth SVF up to a RMSE of the order of the images resolution, as shown in Table. 1. In addition, our method recovered the incompresibility with a higher accuracy.

4 Evaluation for myocardial tracking

We evaluate the proposed method on the STACOM’11 myocardial tracking challenge [10]. In particular, this allows a direct comparison of our framework with iLogDemons[7] a state-of-the-art incompressible registration method.

4.0.1 Data

The STACOM’11 dataset222http://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.

4.0.2 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.

4.0.3 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. 2, are not statistically different from the results of iLogDemons [7] in terms of landmark tracking error. The log-Jacobian map of Fig. 1 illustrates our approximation results (8): for a small deformation, our registration framework can guarantee an incompressible deformation in a subregion with high accuracy. While Fig. 2 illustrates the degradation of this result for larger deformation, i.e. when registering frames at the opposite of the cardiac cycle.

5 Conclusion and Future Work

5.0.1 Limitations

Similar to [6, 7], 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.

5.0.2 Advantages compared to previous methods

Our method for incompressible registration relies on the parametrisation of the velocity field by a divergence conforming B-splines [4]. In contrast to classical B-splines used in [6, 7], 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 (8). 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 [6, 7] 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.

Acknowledgments

This project has received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Skł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\7\34].

1 B-splines notations

To simplify the notations, let us assume that:

• •

$\Omega=[0,1]^{3}\subset\mathds{R}^{3}$ is the continuous spatial domain of the images,

• •

$k\in\mathds{N}$, the order of the B-splines basis, is less than $3$ (higher order values are rarely used in practice)

• •

$\left(\delta x,\delta y,\delta z\right)=(\frac{1}{n},\frac{1}{n},\frac{1}{n})$ is a regular spacing with $n>k$,

• •

$\Omega_{grid}^{k}$ is a regular grid of knots $\{(x_{i_{X}},y_{i_{Y}},z_{i_{Z}})=(\frac{i_{X}}{n},\frac{i_{Y}}{n},\frac{i_{Z}}{n})\}_{i_{X},i_{Y},i_{Z}=-k}^{n-1}$ for the spacing $\left(\delta x,\delta y,\delta z\right)$ on $[-\frac{k}{n},1-\frac{1}{n}]^{3}$.

We use the following notations for the B-splines basis functions of order $k$:

[TABLE]

where the $B^{k}$ for $k\in\{0,1,2,3\}$ are defined by:

[TABLE]

Using those notations, in equation (7) we define a 3D divergence-conforming SVF $v$ of order $k\geq 2$ and parameters $\left\{(\phi^{X}_{i},\phi^{Y}_{i},\phi^{Z}_{i})\right\}_{i\in\{-k,\ldots,n-1\}^{3}}\in\mathds{R}^{3(n+k)^{3}}$ as, for all $(x,y,z)\in\Omega$,

[TABLE]

Remark: Support of B-splines basis functions

The support of a function $f:\mathcal{X}\mapsto\mathds{R}$ is the set of points where $f$ is non-zero, i.e. $\operatorname*{supp}\left(f\right)=\left\{x\in\mathcal{X}\,|\,f(x)\neq 0\right\}$. For any direction $U\in\{X,Y,Z\}$, for $i\in\{-k,\ldots,n-1\}^{3}$, the support of the function $B^{k}_{i,U}$ is $]u_{i},\,u_{i}+(k+1)\delta u[$.

2 Proof of Lemma 1

In this subsection, we give more details on the proof of Lemma 1. We start by showing equation (5). The functions $B^{k}$ have the property:

[TABLE]

So using the chain rule, we obtain for all $k\geq 2$ and $i=(i_{X},i_{Y},i_{Z})\in\{1,\ldots,N\}^{3}$,

[TABLE]

Similarly, for the directions $X,Y$ we obtain:

[TABLE]

Thus, if $k\geq 2$, for all $(x,y,z)\in\Omega=[0,1]^{3}$:

[TABLE]

Separating the sum into two sums and using the change of index $i_{X}:=i_{X}+1$ in the second sum, we obtain:

[TABLE]

We have $B^{k-1}_{-k}(x)=B^{k-1}_{n}(x)=0$ for all $x\in[0,1]$ (see remark in 1.1), so we can group the two sums and we finally obtain:

[TABLE]

Similarly we obtain, for all $(x,y,z)\in\Omega=[0,1]^{3}$:

[TABLE]

As a result, for all $(x,y,z)\in\Omega=[0,1]^{3}$, the divergence of $v$ at $(x,y,z)$ is given by:

[TABLE]

We are now ready to prove Lemma 1. Let us assume that $k\geq 2$, and $\mathcal{M}$ is a non-empty subset of $\Omega$. Let $\epsilon\geq 0$ and

[TABLE]

$\mathbf{J}_{\mathcal{M}}$ contains all the indices $(i_{X},i_{Y},i_{Z})$ so that the support of the function $(x,y,z)\mapsto B^{k-1}_{i_{X},X}(x)B^{k-1}_{i_{Y},Y}(y)B^{k-1}_{i_{Z},Z}(z)$ is non-zero for at least one point of $\mathcal{M}$.

As a result, if for all $(i_{X},i_{Y},i_{Z})\in\mathbf{J}_{\mathcal{M}},$

[TABLE]

Then, for all $(x,y,z)\in\mathcal{M}$,

[TABLE]

Which concludes the proof of Lemma 1.

3 Proof for Lie exponential approximation and incompressibility

In this subsection, we give a detailed proof for the error of incompressibility with an Euler approximation of the Lie exponential (8). Let $K\in\mathds{N}$, let $\tau=\frac{1}{2^{K}}$ be the time step in the Euler integration. The Euler approximation of the Lie exponential for the time step $\tau$ is given by:

[TABLE]

where $I$ is the identity mapping.

Let $\mathcal{M}$ be a non-empty subregion of the spatial domain $\Omega$ and $m=(m_{X},m_{Y},m_{Z})\in\mathcal{M}$. The Jacobian of the first step of the Euler integration $\left(I+\frac{1}{2^{K}}v\right)$ at $m$ is given by:

[TABLE]

Let us note:

[TABLE]

Using the chain rule and the fact that the determinant of the composition of two matrices is equal to the product of their determinant, we obtain:

[TABLE]

Let us now assume that $v$ satisfies the condition of Lemma 1 for $\epsilon$ close to [math], up to machine precision (typically $\epsilon=\epsilon_{mach}=10^{-16}$). If in addition, $\forall k\in\{0,\ldots,2^{K}-1\},m_{k}\in\mathcal{M}$. Then, $\forall k\in\{0,\ldots,2^{K}-1\},\absolutevalue{\divergence{v}(m_{k})}\leq\epsilon_{mach}$, and

[TABLE]

As a result, we obtain the approximation given in (8):

[TABLE]

Remark: divergence-free SVF does not guarantee incompressibility of the transformation when $\mathcal{M}$ is local and the deformation is large

In the case at least one of the $m_{k}$ is outside of $\mathcal{M}$, inequalities (10) do not hold anymore. Although the divergence of the SVF $v$ is uniformly close to [math] on $\mathcal{M}$, some point can be deformed in a compressible manner. This happens when the deformation is large and points that were initially inside $\mathcal{M}$ end up outside of $\mathcal{M}$ during the Euler integration. In practice, deviations from an incompressible deformation in $\mathcal{M}$ can be observe for accurate divergence-free SVFs in $\mathcal{M}$, as illustrated in Fig. 2.

Fortunately, the definition of $J_{\mathcal{M}}$ implies the presence of margins around the incompressible region $\mathcal{M}$. As a consequence, inequalities (10) are verified in practice for small and moderate deformations in $\mathcal{M}$, as illustrated in Fig. 1. One can verify that thoses margins are linear in the order $k$ of the B-splines basis.

It is worth noting that this limitation is due to the use of SVFs, also used in previous works. In addition, when $\mathcal{M}$ is equal to the entire image domain $\Omega$, inequalities (10) are always satisfied for a divergence-free SVF obtained by our method.