# A Generalized Framework for Analytic Regularization of Uniform Cubic
B-spline Displacement Fields

Keyur D. Shah1, James A. Shackleford1, Nagarajan Kandasamy1, Gregory C. Sharp2
1 Electrical and Computer Engineering Department, Drexel University,
Philadelphia, PA 19104, USA  2 Department of Radiation Oncology, Massachusetts
General Hospital, Boston, MA 02114, USA<EMAIL_ADDRESS>
###### Abstract

Image registration is an inherently ill-posed problem that lacks the
constraints needed for a unique mapping between voxels of the two images being
registered. As such, one must regularize the registration to achieve
physically meaningful transforms. The regularization penalty is usually a
function of derivatives of the displacement-vector field, and can be
calculated either analytically or numerically. The numerical approach,
however, is computationally expensive depending on the image size, and
therefore a computationally efficient analytical framework has been developed.
Using cubic B-splines as the registration transform, we develop a generalized
mathematical framework that supports five distinct regularizers: diffusion,
curvature, linear elastic, third-order, and total displacement. We validate
our approach by comparing each with its numerical counterpart in terms of
accuracy. We also provide benchmarking results showing that the analytic
solutions run significantly faster — up to two orders of magnitude — than
finite differencing based numerical implementations.

Keywords: deformable registration, B-spline registration, regularization,
analytic regularization

##  1 Introduction

The goal of image registration is to find a geometric transformation between
corresponding image data that brings them into a common coordinate frame. By
fusing multiple images, a physician gains a more complete understanding of
patient anatomy. The images can be acquired using similar or different imaging
modalities — for example, CT or MRI — and may represent different stages of
growth or disease. A registration is called _rigid_ if the motion or change is
limited to global rotations and translations, and _deformable_ when the
registration includes complex local variations. Deformable registration is
preferred over rigid when locally precise alignment is needed; for example, in
image-guided surgery (Hartkens et al., 2003) and image-guided radiotherapy
(Zhang et al., 2007).

Given a three-dimensional fixed image $F$ with voxel coordinates
$\boldsymbol{x}=(x_{1},x_{2},x_{3})$ and voxel intensity
$f=F(\boldsymbol{x})$, and moving image $M$ with voxel coordinates
$\boldsymbol{x^{\prime}}=(x_{1}^{\prime},x_{2}^{\prime},x_{3}^{\prime})$ and
voxel intensity $m=M(\boldsymbol{x^{\prime}})$ representing the same
underlying anatomy as $F$ within the image overlap domain
$\boldsymbol{\Omega}$, the two images $F$ and $M$ are said to be registered
when the cost function

$C=\sum_{\boldsymbol{T}(\boldsymbol{x})\in\boldsymbol{\Omega}}{\boldsymbol{\Psi}\left(f,m\right)}$
(1)

is minimized according to a similarity metric $\boldsymbol{\Psi}$ under the
coordinate mapping
$\boldsymbol{T}(\boldsymbol{x})=\boldsymbol{x}+\boldsymbol{\nu}$. Here
$\boldsymbol{\nu}$ is the dense displacement field defined for every voxel
$\boldsymbol{x}\in\boldsymbol{\Omega}$, which is assumed capable of providing
a good diffeomorphism from $F$ to $M$. A diffeomorphism is a globally one-to-
one smooth and continuous mapping with derivatives that are invertible
(Ashburner, 2007). Deformable image registration is an inherently ill-posed
problem and so the unconstrained formulation in (1) can lead to physically
unrealistic transforms such as the one shown in Figure 1. Moving image, Figure
1 (b) is registered to the fixed image, Figure 1 (a). However, the resulting
warped image Figure 1 (c) exhibits areas of irregular compression and
expansion that are marked using red circles.

Figure 1: Example of intra-subject registration demonstrating the need for
regularization.

It is desirable to confine the solution space to prevent physically
unrealistic transforms. This can be done by introducing a penalty term which
regularizes the transformation. We can modify (1) to include the
regularization term as

$C=\sum_{\boldsymbol{T}(\boldsymbol{x})\in\boldsymbol{\Omega}}{\boldsymbol{\Psi}\left(f,m\right)}+S(\boldsymbol{\nu}),$
(2)

where the smoothness $S(\boldsymbol{\nu})$ is added to $\boldsymbol{\Psi}$ to
drive $\boldsymbol{T}$ to a physically meaningful coordinate mapping. Several
formulations for the regularization term $S$ have been developed in the
literature. Rueckert et al. (1999) penalize high bending energy in thin-plate
spines to achieve smoother local deformations. This technique is implemented
within the Medical Image Registration Toolkit (MIRTK). Rohlfing et al. (2003)
penalize local deviations from a unity Jacobian determinate to preserve
incompressibility of volume regions since soft tissue is incompressible for
small deformations. Linear elastic energy is minimized by Miller et al. (1993)
to ensure that the deformation field generated by the registration is
physically smooth. Chun and Fessler (2009) develop a regularizer based on
sufficient conditions to enforce local invertibility. Cahill et al. (2009)
generalize the Demon’s algorithm (Thirion, 1998) to allow image-driven locally
adaptive regularization. Sorzano et al. (2005) propose a regularizer based on
the curl and divergence of the underlying displacement field as it is the
measure of the true roughness of the deformation. Use of Fourier methods to
solve partial differential equations of various standard regularizers is
studied by Cahill et al. (2007). Burger et al. (2013) develop a regularizer
based on hyper-elasticity in the context of a mass-preserving registration
problem. Using a Demons’ framework, Tustison and Avants (2013) use the
directly manipulated free-form deformation as a regularizer for the resulting
displacement field to provide biologically plausible solutions. Vishnevskiy et
al. (2016) use an isotropic version of the total variation regularization to
correctly represent non-smooth displacement fields, that occur at sliding
interfaces in the thorax and abdomen in image time-series during respiration.
Schmidt-Richberg et al. (2012) and Delmon et al. (2013) proposed a direction-
dependent regularization to estimate slipping organ motion. Miura et al.
(2017) use the anatomically constrained deformation algorithm (ANACONDA) of
RayStation, which penalizes invertibility of the displacement field and any
large shape deviations to the region of interest. Fu et al. (2018) develop an
adaptive direction dependent regularization technique using a Gaussian
isotropic filter and a bilateral filter in order to preserve sliding motion.
Ghaffari and Fatemizadeh (2017) proposed a rank-regularized sum of squared
differences similarity measure in order to overcome the challenge of spatially
varying intensity distortion. Mang and Biros (2016) constrain the divergence
of the velocity field to control the compressibility of the displacement field
for a 2D case.

Returning to the regularized cost function in (2), the penalty term
$S(\boldsymbol{\nu})$ can be calculated either numerically or analytically.
The numerical approach, while simple and flexible, is computationally
expensive depending on the image size, whereas the analytical methods requires
solving a series of logical steps, but is computationaly inexpensive to solve.
Moreover, the numerical methods provide the approximate solution and the
analytical method gives the exact solution and therefore analytic methods are
often preferred. Shackleford et al. (2012) have previously developed an
analytic method to calculate the bending energy of a vector field that is
parameterized via uniform cubic B-spline basis function and report significant
speedup compared to the numerical counterpart. Along similar lines, Shusharina
and Sharp (2012) present an analytic method to regularize the radial basis
function.

While the above prior work has developed analytic methods for specific types
of regularizers, the novelty of this paper lies in the development of a
_generalized mathematical framework_ for this problem. Using cubic B-splines
as the deformable registration transform, our framework accommodates five
distinct types of regularization: diffusion (Thirion, 1998), curvature
(Modersitzki, 2004), linear elastic (Broit, 1981), third-order (Lellmann et
al., 2013), and total displacement. In addition to exhibiting improved
computational efficiency with respect to numerical approaches, a key advantage
of our approach is the ability to seamlessly combine multiple regularizers to
realize custom or domain-specific smoothness constraints as dictated by the
underlying registration problems at no additional cost in performance.

The paper is organized as follows. Section 2 discusses relevant theory and the
regularizers of interest. Section 3 develops the framework and its analytical
solution; which is validated and benchmarked in Section 4. We conclude the
paper in Section 5.

##  2 Background

Here we describe the formulation of the diffusion, curvature, linear elastic,
third-order, and total displacement regularizers for a three-dimensional
displacement field parameterized by a uniform cubic B-spline basis function.

For volumetric or 3D registration, the displacement field at any given voxel
is determined by the $4^{3}=64$ control points in the immediate vicinity of
the voxel. We use the term _tile_ to denote the set of voxels which receives
local support from the same set of 64 control points. The tile forms the
backbone of an analytic expression for the continuous displacement field
$\boldsymbol{\nu}$. B-spline interpolation is performed for each vector within
a tile using the 64 control-point coefficients that provide local support for
the operation. The B-spline interpolation yielding the first component of the
displacement vector for a voxel located at $\boldsymbol{x}$ is

$\nu_{1}(\boldsymbol{x})=\sum_{l=0}^{3}\sum_{m=0}^{3}\sum_{n=0}^{3}\beta_{l}(u_{1})\beta_{m}(u_{2})\beta_{n}(u_{3})p_{1,l,m,n}$
(3)

where $p_{1}$ is one of the 64 B-spline coefficient used to interpolate the
$\nu_{1}$ component of the displacement vector $\boldsymbol{\nu}$ for the
voxel located at $\boldsymbol{x}$. The $\beta_{l}$, $\beta_{m}$, and
$\beta_{n}$ terms represent the uniform cubic B-spline basis functions in the
$x_{1}$, $x_{2}$, and $x_{3}$ directions, respectively, and $u_{1}$, $u_{2}$,
and $u_{3}$ represent the normalized local coordinates of the voxel within its
tile (Shackleford et al., 2010). The uniform cubic B-spline basis function
$\beta_{l}$ along the $x_{1}$ direction is given by

$\displaystyle\beta_{l}(u_{1})=\begin{dcases}\frac{(1-u_{1})^{3}}{6}&l=0\\\
\frac{3{u_{1}}^{3}-6{u_{1}}^{2}+4}{6}&l=1\\\
\frac{-3{u_{1}}^{3}+3{u_{1}}^{2}+3u_{1}+1}{6}&l=2\\\
\frac{{u_{1}}^{3}}{6}&l=3\end{dcases}$ (4)

and similarly for $\beta_{m}$ and $\beta_{n}$ in the $x_{2}$ and $x_{3}$
directions respectively. The displacement vector components $\nu_{2}$ and
$\nu_{3}$ in the remaining two directions are calculated similarly. The
optimizer updates the B-spline coefficients during each iteration until an
optimal registration between the moving and fixed images is achieved. The
regularization penalty term is a function of the displacement vector field
_i.e._ a function of the B-spline coefficients. Thus, the regularization
penalty term for the entire deformation may be expressed as a sum of the
regularization penalty terms over all tiles. Therefore, our approach is to
develop an operator that computes the penalty term for a tile as a function of
its B-spline control points.

We now describe the five regularizers of interest. The functions in (5)
describe the first, second, and third-order partial derivatives of the
displacement field, which are the building blocks of these regularizers:

$\displaystyle f_{1}(\nu_{i},x_{j})$
$\displaystyle=\frac{\partial{\nu_{i}}}{\partial{x}_{j}}$ $\displaystyle
f_{2}(\nu_{i},x_{j},x_{k})$
$\displaystyle=\frac{\partial^{2}\nu_{i}}{\partial{x}_{j}\partial{x}_{k}}$
$\displaystyle f_{3}(\nu_{i},x_{j},x_{k})$
$\displaystyle=\frac{\partial{\nu_{i}}}{\partial{x}_{j}}\frac{\partial{\nu_{i}}}{\partial{x}_{k}}$
$\displaystyle f_{4}(\nu_{i},x_{j},x_{k},x_{q})$
$\displaystyle=\frac{\partial^{3}\nu_{i}}{\partial{x}_{j}\partial{x}_{k}\partial{x}_{q}}$
(5)

Diffusion Regularizer: Thirion originally introduced the demons algorithm as a
diffusion process (Thirion, 1998) and in later work Modersitzki coined the
term diffusion regularizer (Modersitzki, 2004) since this gradient-based
partial differential equation was viewed as a generalized diffusion equation.
The penalty term is given by

$\displaystyle
S_{1}=\int_{\boldsymbol{\Omega}}{\sum_{i,j=1}^{3}f_{1}(\nu_{i},x_{j})^{2}}d\boldsymbol{x},$
(6)

where $\boldsymbol{\Omega}$ denotes the image region over which the
regularization penalty is to be calculated.

Curvature Regularizer: A curvature regularizer uses second-order derivative
terms and aims to reduce the bending energy of the displacement field
(Modersitzki, 2004). The penalty term is specified as

$\displaystyle
S_{2}=\int_{\boldsymbol{\Omega}}\sum_{i,j,k=1}^{3}f_{2}(\nu_{i},x_{j},x_{k})^{2}d\boldsymbol{x}.$
(7)

Linear Elastic Regularizer: The penalty term is given by

$\displaystyle
S_{3}=\int_{\boldsymbol{\Omega}}{\sum_{i,j=1}^{3}f_{1}(\nu_{i},x_{j})^{2}}+{\sum_{i,j,k=1,j\neq
k}^{3}f_{3}(\nu_{i},x_{j},x_{k})}d\boldsymbol{x},$ (8)

which aims to strike a balance between the global registration achieved via
affine mapping versus the more local elastic registration (Broit, 1981).

Third-order Regularizer: The penalty term is specified in terms of third-order
derivative terms as

$\displaystyle
S_{4}=\int_{\boldsymbol{\Omega}}\sum_{i,j,k,o=1}^{3}f_{4}(\nu_{i},x_{j},x_{k},x_{o})^{2}d\boldsymbol{x}.$
(9)

Total displacement Regularizer: The penalty term is specified in terms of the
magnitude of the displacement vector field at each voxel as

$\displaystyle
S_{5}=\int_{\boldsymbol{\Omega}}\sum_{i=1}^{3}{\nu_{i}}^{2}d\boldsymbol{x}.$
(10)

Combining the above-described regularizers, the total smoothness penalty $S$
for the unified framework can be expressed as

$\displaystyle
S=\mu_{1}S_{1}+\mu_{2}S_{2}+\mu_{3}S_{3}+\mu_{4}S_{4}+\mu_{5}S_{5},$ (11)

where $\mu_{1}$, $\mu_{2}$, $\mu_{3}$, $\mu_{4}$, and $\mu_{5}$ are the
weights corresponding to the diffusion, curvature, linear-elastic, third-
order, and total displacement regularizers, respectively. In the simplest
case, one regularizer can be chosen over the others; for example, setting
$\mu_{2}$, $\mu_{3}$, $\mu_{4}$, and $\mu_{5}$ to zero chooses the diffusion
regularizer. Other regularizers can be selected similarly. Several other
regularizers can be derived using this framework as long as the regularizer
penalty terms consists of product of two partial derivative terms of the
displacement filed from (5) or one of the components of the displacement
field.

##  3 Development of Analytical Algorithm

The analytic algorithm developed in this section comprises the following three
major steps: (1) taking the first, second, or third-order derivative of the
displacement vector field; (2) squaring or multiplying the derivative terms;
and (3) integrating the products of the derivative terms over a tile. The
various derivative terms that constitute (11) are recast as simple matrix
operations that can be efficiently performed using Basic Linear Algebra
Subprograms (BLAS). The algorithmic steps are described in greater detail
below.

When represented sparsely via the uniform cubic B-spline basis, the
displacement field $\boldsymbol{\nu}$ is parameterized by the set of B-spline
basis coefficients $\boldsymbol{p_{1}},\boldsymbol{p_{2}},\boldsymbol{p_{3}}$,
where

$\boldsymbol{p_{1}}=\begin{bmatrix}p_{1,0,0,0}\\\ \vdots\\\ p_{1,I,J,K}\\\
\end{bmatrix},\boldsymbol{p_{2}}=\begin{bmatrix}p_{2,0,0,0}\\\ \vdots\\\
p_{2,I,J,K}\\\ \end{bmatrix},\boldsymbol{p_{3}}=\begin{bmatrix}p_{3,0,0,0}\\\
\vdots\\\ p_{3,I,J,K}\\\ \end{bmatrix}$ (12)

are the control points that define the displacement field within a single
tile. For cubic B-splines, $I=J=K=2$. The tile has dimensions of
$\boldsymbol{r}=r_{1}\times r_{2}\times r_{3}$ mm3. To express the first
component of the displacement field $\nu_{1}$ as a function of
$\boldsymbol{p_{1}}$, we start with the matrix $\boldsymbol{B}$ containing the
coefficients for the cubic B-spline basis function as described in (4) and
matrix $\boldsymbol{R_{1}}$ which controls for tile size as

$\displaystyle\boldsymbol{B}=\frac{1}{6}\begin{bmatrix}1&-3&3&-1\\\
4&0&-6&3\\\ 1&3&3&-3\\\
0&0&0&1\end{bmatrix},\;\boldsymbol{R_{1}}=\begin{bmatrix}1&0&0&0\\\
0&\frac{1}{r_{1}}&0&0\\\ 0&0&\frac{1}{r_{1}^{2}}&0\\\
0&0&0&\frac{1}{r_{1}^{3}}\end{bmatrix}.$ (13)

We also generate the matrix $\boldsymbol{\Delta}^{(\delta)}$ which is defined
for $\delta\in$ [0,3] as

$\displaystyle\boldsymbol{\Delta}^{(0)}=\begin{bmatrix}1&0&0&0\\\ 0&1&0&0\\\
0&0&1&0\\\
0&0&0&1\end{bmatrix},\;\boldsymbol{\Delta}^{(1)}=\begin{bmatrix}0&0&0&0\\\
1&0&0&0\\\ 0&2&0&0\\\ 0&0&3&0\end{bmatrix},\;$
$\displaystyle\boldsymbol{\Delta}^{(2)}=\begin{bmatrix}0&0&0&0\\\ 0&0&0&0\\\
2&0&0&0\\\
0&6&0&0\end{bmatrix},\;\boldsymbol{\Delta}^{(3)}=\begin{bmatrix}0&0&0&0\\\
0&0&0&0\\\ 0&0&0&0\\\ 6&0&0&0\end{bmatrix}.\;$ (14)

Matrices $\boldsymbol{Q_{1}}$, $\boldsymbol{Q_{2}}$, and $\boldsymbol{Q_{3}}$
can now be calculated as

$\displaystyle\boldsymbol{Q_{1}}^{(\delta)}=\boldsymbol{B}\boldsymbol{R_{1}}\boldsymbol{\Delta}^{(\delta)},\;\boldsymbol{Q_{2}}^{(\delta)}=\boldsymbol{B}\boldsymbol{R_{2}}\boldsymbol{\Delta}^{(\delta)},\;\mathrm{and}\;\boldsymbol{Q_{3}}^{(\delta)}=\boldsymbol{B}\boldsymbol{R_{3}}\boldsymbol{\Delta}^{(\delta)},$
(15)

to provide a convenient method for obtaining the first, second and third-order
derivatives, $\boldsymbol{\nu^{\prime}}$, $\boldsymbol{\nu^{\prime\prime}}$
and $\boldsymbol{\nu^{\prime\prime\prime}}$, respectively, with respect to the
Euclidean basis as required by the calculation of the smoothness penalty. For
example, $\boldsymbol{Q_{1}}^{(0)}$ is used when $\boldsymbol{\nu}$ is needed
as is, and $\boldsymbol{Q_{1}}^{(1)}$, $\boldsymbol{Q_{1}}^{(2)}$ and
$\boldsymbol{Q_{1}}^{(3)}$ are used when the first, second and third
derivatives, respectively, of $\boldsymbol{\nu}$ are needed in the
calculations. These matrices also map the domain of the B-spline basis
function to lie in the interval [0,1]. Now, the vector $\nu_{1}$ may now be
expressed at a point $\boldsymbol{x}={x_{1},x_{2},x_{3}}$ using the 64
B-spline coefficients supporting $\boldsymbol{x}$ as the tensor product:

$\displaystyle\nu_{1}=\sum_{i,j,k=0}^{3}{p}_{i,j,k}\sum_{a=0}^{3}\boldsymbol{Q_{1}}^{(0)}\boldsymbol{x_{1}}\sum_{b=0}^{3}\boldsymbol{Q_{2}}^{(0)}\boldsymbol{x_{2}}\sum_{c=0}^{3}\boldsymbol{Q_{3}}^{(0)}\boldsymbol{x_{3}}$
(16)

where the four rows of $\boldsymbol{Q_{1}}^{(\delta)}$ can be separated into
vectors as

$\boldsymbol{Q_{1}}^{(\delta)}=\begin{bmatrix}q_{1,0}^{\boldsymbol{T}}\\\
q_{1,1}^{\boldsymbol{T}}\\\ q_{1,2}^{\boldsymbol{T}}\\\
q_{1,3}^{\boldsymbol{T}}\end{bmatrix}^{(\delta)}$ (17)

We define

$\boldsymbol{x_{1}}=[1\quad x_{1}\quad x_{1}^{2}\quad
x_{1}^{3}]^{\boldsymbol{T}}$ (18)

Cartesian basis vectors $\boldsymbol{x_{2}}$ and $\boldsymbol{x_{3}}$ are
defined similarly. Finally, vectors $\nu_{2}$ and $\nu_{3}$ are defined in a
similar fashion to $\nu_{1}$ to complete the B-spline interpolation operation.

Now that we have the basic mathematical framework in place, the next steps of
the algorithm are designed to calculate the various squared and product terms
listed in (11) — specifically, to square or to multiply the
$\sum\boldsymbol{Q_{1}}^{(\delta)}\sum\boldsymbol{Q_{2}}^{(\delta)}\sum\boldsymbol{Q_{3}}^{(\delta)}$
sub-term in (16).

A $4\times 4$ matrix is calculated as the outer product

$\boldsymbol{\Xi_{1,a,b}}\ =q_{1,a}^{\delta_{i}}\otimes q_{1,b}^{\delta_{j}},$
(19)

where $a,b\in\\{0,1,2,3\\}$ and $i,j\in\\{1,2,3\\}$. Here, $\delta_{i}$ and
$\delta_{j}$ are the same when taking the square of the derivative term
whereas they are different when taking the product of two distinct terms —
which is the case for the linear elastic regularizer. The $\delta$ terms
indicate the variable on which the partial derivative of the displacement
field is obtained. For example, if $i$ is 1, the partial derivative is taken
with respect to the first component, and so on. The matrices
$\boldsymbol{\Xi_{2,a,b}}$ and $\boldsymbol{\Xi_{3,a,b}}$ for the second and
third components can be calculated similarly. For the ease of readability, the
$\delta$ terms are not included in the follow-up equations.

Taking the outer product of the Cartesian basis vector $\boldsymbol{x_{1}}$ as
defined in (18) with itself results in a $4\times 4$ matrix

$\boldsymbol{X_{1}}=\boldsymbol{x_{1}}\otimes\boldsymbol{x_{1}}=\begin{bmatrix}1&x_{1}&x_{1}^{2}&x_{1}^{3}\\\
x_{1}&x_{1}^{2}&x_{1}^{3}&x_{1}^{4}\\\
x_{1}^{2}&x_{1}^{3}&x_{1}^{4}&x_{1}^{5}\\\
x_{1}^{3}&x_{1}^{4}&x_{1}^{5}&x_{1}^{6}\end{bmatrix},$ (20)

and $\boldsymbol{X_{2}}$ and $\boldsymbol{X_{3}}$ can be obtained similarly.
Taking the Hadamard product of (19) and (20) results in $4\times 4$ matrix

$\psi_{1,a,b}\ =\boldsymbol{\Xi_{1,a,b}}\ \odot\boldsymbol{X_{1}}.$ (21)

The elements of $\psi_{1,a,b}$ can be combined into a matrix
$\boldsymbol{\Gamma_{1}}$, where each element $\boldsymbol{\Gamma_{1}}(a,b)$
is formed as

$\boldsymbol{\Gamma_{1}}(a,b)=\sum_{i=1}^{4}\sum_{j=1}^{4}\psi_{1,a,b}(i,j).$
(22)

Since there are 16 combinations for $a$ and $b$, $\boldsymbol{\Gamma_{1}}$ is
a $4\times 4$ matrix; $\boldsymbol{\Gamma_{2}}$ and $\boldsymbol{\Gamma_{3}}$
are obtained similarly, allowing for the desired composite matrix operator

$\boldsymbol{\Gamma}=\boldsymbol{\Gamma_{1}}\otimes\boldsymbol{\Gamma_{2}}\otimes\boldsymbol{\Gamma_{3}},$
(23)

to calculate the smoothness metric over a tile, for the specified choice of
$\delta$’s, as

$\int_{0,0,0}^{r_{1},r_{2},r_{3}}\sum\boldsymbol{p_{i}}^{\boldsymbol{T}}\boldsymbol{\Gamma}\boldsymbol{p_{j}}d\boldsymbol{x}.$
(24)

Since B-spline coefficients are constants, we can rewrite the above expression
as

$\sum\boldsymbol{p_{i}}^{\boldsymbol{T}}\boldsymbol{V}\boldsymbol{p_{j}},$
(25)

where
$\boldsymbol{V}=\int_{0,0,0}^{r_{1},r_{2},r_{3}}\boldsymbol{\Gamma}d\boldsymbol{x}$
and $i,j\in{1,2,3}$. The smoothness penalty $S$ for the entire volume can be
calculated as the sum of all its constituent tiles.

This analytic implementation is extremely memory efficient. The 32
$\boldsymbol{V}$ matrices each having dimensions of $64\times 64$
corresponding to every term in the five regularizers are calculated beforehand
and reused in the entire optimization process. During the optimization step
only $\boldsymbol{p^{T}}\boldsymbol{V}\boldsymbol{p}$ is calculated for each
tile according to (25). The overall memory requirement is 512 kB
($32\times\boldsymbol{V}=32\times 64\times 64\times 4\ kB)$.

Illustrative Examples: Given the above-described formalism, the generalized
equation for squaring or multiplying the various derivative terms can be
written in terms of the B-spline basis coefficients and the composite matrix
operator $\boldsymbol{\Gamma}$ as

${\left(\frac{\partial^{n}\nu_{i}}{\partial{x_{i}}}\frac{\partial^{n}\nu_{j}}{\partial{x_{j}}}\right)}=\boldsymbol{p_{i}}^{\boldsymbol{T}}\boldsymbol{\Gamma}^{n\times(\boldsymbol{\delta_{i}},\boldsymbol{\delta_{j}})}\boldsymbol{p_{j}},$
(26)

where $n$ determines the order of the partial derivative and
$\boldsymbol{\delta_{i}}=(\delta_{i1},\delta_{j1},\delta_{k1})$ and
$\boldsymbol{\delta_{j}}=(\delta_{i2},\delta_{j2},\delta_{k2})$. Individual
terms of $\boldsymbol{\delta_{i}}$ depend on the variable on which the partial
derivative of the displacement field is obtained. For example, $\delta_{i1}=1$
if the partial derivative is taken with respect to the first component. For
the squared terms $\boldsymbol{\delta_{i}}$ and $\boldsymbol{\delta_{j}}$ are
the same. The terms corresponding to the diffusion, curvature, linear-elastic,
third-order, and total displacement regularizers are described by the
following equations using the appropriate $\delta$ values discussed earlier in
(19) and (26).

$\displaystyle{\left(\frac{\partial\nu_{2}}{\partial{x_{1}}}\right)}^{2}$
$\displaystyle=\boldsymbol{p_{2}}^{\boldsymbol{T}}\boldsymbol{\Gamma}^{(100,100)}\boldsymbol{p_{2}}$
$\displaystyle{\left(\frac{\partial^{2}\nu_{1}}{\partial{x_{1}}\partial{x_{3}}}\right)}^{2}$
$\displaystyle=\boldsymbol{p_{1}}^{\boldsymbol{T}}\boldsymbol{\Gamma}^{(101,101)}\boldsymbol{p_{1}}$
$\displaystyle{\left(\frac{\partial\nu_{1}}{\partial{x_{1}}}\frac{\partial\nu_{3}}{\partial{x_{3}}}\right)}$
$\displaystyle=\boldsymbol{p_{1}}^{\boldsymbol{T}}\boldsymbol{\Gamma}^{(100,001)}\boldsymbol{p_{3}}$
(27)
$\displaystyle{\left(\frac{\partial^{3}\nu_{3}}{\partial{x_{1}}^{3}}\right)}^{2}$
$\displaystyle=\boldsymbol{p_{3}}^{\boldsymbol{T}}\boldsymbol{\Gamma}^{(300,300)}\boldsymbol{p_{3}}$
$\displaystyle{\left(\nu_{3}\right)}^{2}$
$\displaystyle=\boldsymbol{p_{3}}^{\boldsymbol{T}}\boldsymbol{\Gamma}^{(000,000)}\boldsymbol{p_{3}}$

Using similar reasoning, the linear elastic regularization penalty would be
computed as the sum of twelve vector-matrix-vector products of the form
$\boldsymbol{p^{T}}\boldsymbol{V}\boldsymbol{p}$.

##  4 Experiments and Results

We quantify performance of the developed methods in terms of: accuracy and
speedup. The DIR-Lab dataset used in our experiments consists of 4D-CT images
from ten patients who were treated for malignancies in the esophagus or lung
(Castillo et al., 2009). The CT images were masked to only include the lungs,
trachea and bronchi. The maximum inhalation phase (fixed image) was registered
with the maximum exhalation phase (moving image) for intra-subject
registrations. Five of the ten studies have volumes of $512\times 512\times
128$ voxels with physical separation of $0.92\times 0.92\times 2.5$ mm, and
the remaining volumes have $256\times 256\times 94$ voxels. In each case, 300
landmark points were placed within the lung by a medical expert. Figure 2
shows a coronal slice for an intra-subject registration example with a
landmark point overlaid. The moving landmarks are warped using the underlying
transformation $\mathbf{T}$ to get the warped landmarks. Registration accuracy
is measured as the average separation between the fixed and warped landmark
points.

Figure 2: Example of intra-subject registration showing a coronal slice of the
(a) fixed, (b) moving and (c) warped images, with the corresponding selected
landmark point overlaid.

The regularizers developed in this paper have been implemented within
Plastimatch, an open source software for image computation with focus on high-
performance volumetric registration of medical images. Plastimatch is
distributed under a BSD-style license and can be downloaded from
www.plastimatch.org. The fixed and moving images are registered using a three-
stage pyramidal registration approach. Grid spacing for the first and second
stage was kept constant at 100 mm and 80 mm, respectively. Grid spacing for
the third stage was varied from a coarse value of 60 mm to a finer value of 10
mm. Registrations are performed using the mean-squared error similarity
metric, penalized by $S_{n}$ with weight $\mu_{n}$. Referring to (11), we
choose a single regularization strategy over others by setting the remaining
weights to zero. For example, setting $\mu_{1}$, $\mu_{3}$ , $\mu_{4}$, and
$\mu_{5}$ to zero chooses the curvature regularizer. The B-spline coefficients
$\boldsymbol{P}$ describing the transform $\boldsymbol{T}$ are optimized via
the L-BFGS-B optimizer using an analytically computed cost function and
gradient.

Accuracy is measured using the _mean landmark separation_ (_MLS_) between the
inhale landmarks and warped exhale landmarks (or vice versa) after
registration, and smoothness is measured using the _minimum Jacobian
determinant_ of the resulting displacement vector field over the entire
volume. Experiments were performed to assess both quantities as a function of
control-point spacing as well as $\mu_{n}$.

Table 1 shows the relative difference (in %) between the MLS achieved by the
analytic and numeric implementations. The maximum relative difference of 7.4%
occurs for the third-order regularizer. This is because voxels located along
the image boundaries are not used to calculate the smoothness penalty in case
of the finite differencing numeric solutions. Figure 3 shows accuracy and
smoothness results for the curvature, linear elastic, and third-order
regularizers. For a given control-point spacing, a desirable $\mu_{n}$ is one
which produces the best compromise between a small MLS and smooth
$\boldsymbol{T}$. A smooth $\boldsymbol{T}$ can be inferred by a positive
minimum Jacobian determinant. A smooth nonrigid transformation ensures that
the warped image after registration is free from unrealistic compression
and/or expansion artifacts (Chun and Fessler, 2009). Considering the curvature
regularizer, for example, the MLS achieved when $\mu_{2}=10^{-3}$ is very
similar without application of any regularization penalty, and smaller weights
need not be explored. Conversely, the MLS achieved when $\mu_{2}=10^{3}$ is
nearly equal to the MLS prior to registration and changes little when
$\mu_{2}$ is increased further. Thus, weights larger than $10^{3}$ are not
shown. Lower and upper bounds for the weights of the remaining four
regularizers are determined similarly. Additionally looking at the minimum
Jacobian determinant heatmaps from Figure 3 (b), the increased need of
regularization at smaller control point spacing can be inferred. The minimum
Jacobian determinant of the resulting displacement vector field without any
regularization penalty is negative for control point spacing of 10 mm for all
the three regularizers shown. This is due to the higher degrees of freedom
available to the displacement vector field at such smaller control point
spacing.

Table 1: Relative difference (in %) between MLS achieved by the proposed analytic schemes versus numerical finite differencing of the displacement field. Volume size | Grid spacing | Diffusion | Curvature | Elastic | Third-order | Tot. displacement
---|---|---|---|---|---|---
$256\times 256\times 94$ | $20\times 20\times 20$ | 1.5 | 1.9 | 1.3 | 0.9 | 4.4
$256\times 256\times 94$ | $30\times 30\times 30$ | 0.1 | 3.1 | 0.2 | 3.9 | 4.1
$512\times 512\times 128$ | $20\times 20\times 20$ | 2.3 | 2.3 | 3.0 | 7.4 | 1.9
$512\times 512\times 128$ | $30\times 30\times 30$ | 5.9 | 2.0 | 3.9 | 5.7 | 2.9

(a)

(b)

Figure 3: Mean landmark separation (a) and and minimum Jacobian determinant of the transform $\boldsymbol{T}$ (b) as a function of B-spline control-point spacing and the weight $\mu_{n}$, over the 10 thoracic cases using curvature (left), linear elastic (middle), and third-order (right) regularizer for the masked images.  Table 2: Lowest MLS value achieved for each of the ten cases and the corresponding registration configuration. Case # | MLS (mm) | Regularizer ($\mu$) | Grid Spacing (mm)
---|---|---|---
1 | 1.15 | Diffusion ($10^{-2}$) | 30
2 | 1.28 | Diffusion ($10^{-3}$) | 10
3 | 1.52 | Linear Elastic ($10^{-8}$) | 10
4 | 1.56 | Curvature ($10^{-1}$) | 20
5 | 1.80 | Curvature ($10^{-2}$) | 10
6 | 1.56 | Curvature ($10^{-1}$) | 10
7 | 1.74 | Curvature ($10^{-2}$) | 10
8 | 1.49 | Curvature ($10^{-2}$) | 10
9 | 1.50 | Third-order ($1$) | 10
10 | 1.42 | Curvature ($10^{-2}$) | 10

Table 2 summarizes the results for each of the ten cases by providing the
least MLS achieved and the registration configuration used. Curvature
regularizer works best for most of the cases, but other registration problems
may be best optimized with a different regularization choice.

(a) Image overlay

(b) $\boldsymbol{T}$ with no regularization

(c) $\boldsymbol{T}$ for $\mu_{2}=10^{-3}$

(d) $\boldsymbol{T}$ for $\mu_{2}=10^{-2}$

Figure 4: Coronal slice showing: (a) the fixed and the moving CT images
overlaid; (b) displacement field generated without regularization; (c)–(d)
displacement field generated by the curvature regularizer for different values
of $\mu_{2}$. Displacement field shown in black indicates displacement less
than 1 mm, green indicates displacement of around 2 mm, yellow indicates
displacement of around 5 mm, and red indicates displacement greater than 10
mm,

The effect of the penalty factor for the curvature regularizer $\mu_{2}$ on
the transformation $\boldsymbol{T}$ is qualitatively shown in Figure 4 (c-d).
The inhaled and exhaled thoracic volumes are registered using a B-spline
control-point spacing of $10\times 10\times 10$ mm while varying $\mu_{n}$,
and the resulting transform exhibits increased smoothness upon application of
the penalty term. Taking the curvature regularizer for example, we see in
Figure 3 (a)-left that for a control point spacing of 10 mm, the MLS is
comparable when $\mu_{2}$ is either $10^{-3}$ or $10^{-2}$. Inspecting the
minimum Jacobian determinant in Figure 3 (b)-left, we see that the minimum
Jacobian determinant is negative when $\mu_{2}$ is $10^{-3}$, which indicates
a non-smooth displacement field and is positive when $\mu_{2}$ is $10^{-2}$,
which indicates a smooth displacement field. This makes $\mu_{2}=10^{-2}$ a
better choice for the regularization penalty term. This can also be verified
by comparing Figures 4 (c) and 4 (d). Quantitatively, if the registration
accuracy in terms of MLS is the same, the smooth displacement field with a
positive minimum Jacobian determinant is preferred. Similar effect is observed
for the other four regularizers.

The main advantage of an analytic regularization method is significantly
reduced computational time compared to the numerical approach. Computation
times incurred by the analytic solutions depend only on the number of tiles
defined by the B-spline control-point spacing and not on the number of voxels,
drastically reducing the complexity. The composite V matrices can be
calculated once for each control point spacing, and reused throughout the
optimization.

Table 3: Execution times for the analytic methods and corresponding speedup
over numerical approach.

Volume Size | Grid Spacing | Analytic (s) | Diffusion | Curvature | Elastic | Third-order | Tot. displacement
---|---|---|---|---|---|---|---
$256\times 256\times 94$ | $20\times 20\times 20$ | 0.43 | 4.3x | 13.8x | 4.3x | 29.6x | 0.8x
$256\times 256\times 94$ | $30\times 30\times 30$ | 0.14 | 13.6x | 42.7x | 13.3x | 91.6x | 2.6x
$512\times 512\times 128$ | $20\times 20\times 20$ | 1.63 | 6.3x | 19.8x | 6.3x | 41.9x | 1.5x
$512\times 512\times 128$ | $30\times 30\times 30$ | 0.67 | 15.6x | 48.2x | 14.7x | 100.4x | 2.9x

Table 3 lists the execution times incurred by the analytic regularization and
corresponding speedup compared to numerical solutions based on finite
differencing. The benchmarking results reported here were performed using a
machine equipped with dual Intel Octo-core Xeon processors clocked at 2.4 GHz
and 512 GB of main memory. Examining the third column of the table, note that
execution time for a given volume size and grid spacing is the same for any of
the analytic regularizers. This is because there is one composite matrix V per
partial derivative for a total of 32 such matrices in the generalized
framework, all of which are computed irrespective of the regularizer to be
used. The regularizer of interest is later selected by setting a non-zero
weight to the corresponding term in (11). Columns four through eight of the
table list speedup over numerical implementations of the various regularizers.
We consider single-threaded implementations here. Speedup depends on three
factors: volume size, grid spacing, and complexity of the numerical solution.
For a volume size of $512\times 512\times 128$ voxels with grid spacing of
$20\times 20\times 20$ mm. The B-spline grid has 7744 tiles; with grid spacing
$30\times 30\times 30$ mm it has 3177 tiles. In theory, execution time is
linear with the number of tiles.

Referring back to Table 3, third-order and curvature regularizers exhibit
higher speedup than the other regularizers due to the high complexity of their
numerical solution. For example, for a grid-spacing of 20 mm, the analytic
form of the curvature regularizer executes in 450 ms whereas the numerical
form requires 6.2 seconds. On the other hand, the analytic form of the total
displacement regularizer executes in 450 ms but the numeric form only takes
380 ms because regularization involves simply summing the squared vector-field
magnitude at each voxel as per (10).

$0$$2000$$4000$$6000$$8000$$10000$$10^{-2}$$10^{-1}$$10^{0}$Number of
tilesExecution time (seconds)Analytic (1 thread)Analytic (2 threads)Analytic
(4 threads)Analytic (8 threads)Analytic (16 threads)

(a)

$2^{0}$$2^{1}$$2^{2}$$2^{3}$$2^{4}$$10^{-2}$$10^{-1}$$10^{0}$Number of
threadsExecution time (seconds)Number of tiles $=8\times 8\times 8$Number of
tiles $=10\times 10\times 10$Number of tiles $=22\times 2\times 22$

(b)

Figure 5: Execution times for OpenMP implementations as a function of number
of tiles (a) and function of thread count for different numbers of tiles (b).

Finally, the analytical formulation developed in the previous section is
easily parallelized. Returning to (25), this expression can be calculated for
each tile in parallel to calculate partial sums which can then be reduced to a
single value to obtain the smoothness penalty for the entire volume. We
compare the execution time incurred by the single-threaded implementation
against a multi-threaded one using OpenMP. Execution time for a $512\times
512\times 512$ image is shown in Figure 5 as as a function of the number of
number of tiles in the image as well as the number of threads (from one to
sixteen) for the linear-elastic regularizer. Experiments are repeated twenty
times to avoid any discrepancy in the timings caused by a cold cache and the
Figure 5 shows average execution times. The effect of varying the control-
point grid spacing on the execution time of the analytic implementation of the
linear-elastic regularizer for different number of threads can be seen in
Figure 5 (a). Notice the nearly linear speed-up of about $16x$ when the number
of threads is increased to sixteen, which is the number of cores available on
the system used to produce these benchmarks. Figure 5 (b) shows variation in
execution time as a function of thread count, for specific grid-spacing
settings. Execution time decreases when the control-point grid spacing is
increased since this reduces the number of tiles in the overall image. Since
the underlying implementation of all five regularizer is similar, so no
comparison of performance is provided.

##  5 Conclusions

We have developed a fast and general framework which supports five unique
regularizers to calculate the smoothness penalty using an analytical approach
— specifically, by deriving composite matrix operators that operate on a set
of 64 B-spline control points to calculate the regularization penalty within a
given region of support. In terms of accuracy, the maximum relative difference
between the analytical and numerical solutions was $7.4\%$. Furthermore, the
analytic solutions run up to two orders of magnitude faster than finite
differencing based numerical solutions. Fast analytical methods such as these
provide effective regularization without imposing a computational burden to
the deformable image registration pipeline.

This material is based upon work supported by the National Science Foundation
under Grant Nos. 1553436, 1642345 and 1642380 and the National Institutes of
Health under NCI R01CA229178.

## References

  * Ashburner (2007) Ashburner, J., 2007. A fast diffeomorphic image registration algorithm. Neuroimage 38 (1), 95–113.
  * Broit (1981) Broit, C., 1981. Optimal registration of deformed images. Ph.D. thesis, University of Pennsylvania, Philadelphia, PA, USA.
  * Burger et al. (2013) Burger, M., Modersitzki, J., Ruthotto, L., Jan. 2013. A hyperelastic regularization energy for image registration. SIAM Journal on Scientific Computing 35 (1), B132–B148.
  * Cahill et al. (2007) Cahill, N. D., Noble, J. A., Hawkes, D. J., Jun. 2007. Fourier methods for nonparametric image registration. In: IEEE Conf. Computer Vision & Pattern Recognition. pp. 1–8.
  * Cahill et al. (2009) Cahill, N. D., Noble, J. A., Hawkes, D. J., Sep. 2009. A demons algorithm for image registration with locally adaptive regularization. In: Int’l Conf. Medical Image Computing & Computer-Assisted Intervention. pp. 574–581.
  * Castillo et al. (2009) Castillo, E., Castillo, R., Martinez, J., Shenoy, M., Guerrero, T., Jan. 2009. Four-dimensional deformable image registration using trajectory modeling. Physics in Medicine & Biology 55 (1), 305–327.
  * Chun and Fessler (2009) Chun, S. Y., Fessler, J. A., Feb. 2009. A simple regularizer for B-spline nonrigid image registration that encourages local invertibility. IEEE Journal Selected Topics Signal Proc. 3 (1), 159–169.
  * Delmon et al. (2013) Delmon, V., Rit, S., Pinho, R., Sarrut, D., 2013. Registration of sliding objects using direction dependent b-splines decomposition. Physics in Medicine & Biology 58 (5), 1303.
  * Fu et al. (2018) Fu, Y., Liu, S., Li, H. H., Li, H., Yang, D., 2018. An adaptive motion regularization technique to support sliding motion in deformable image registration. Medical physics 45 (2), 735–747.
  * Ghaffari and Fatemizadeh (2017) Ghaffari, A., Fatemizadeh, E., 2017. Image registration based on low rank matrix: Rank-regularized SSD. IEEE transactions on medical imaging 37 (1), 138–150.
  * Hartkens et al. (2003) Hartkens, T., Hill, D. L., Castellano-Smith, A. D., Hawkes, D. J., Maurer, C., Martin, A. J., Hall, W. A., Liu, H., Truwit, C. L., 2003. Measurement and analysis of brain deformation during neurosurgery. IEEE Transactions on Medical Imaging 22 (1), 82–92.
  * Lellmann et al. (2013) Lellmann, J., Morel, J.-M., Schönlieb, C.-B., Jun. 2013. Anisotropic third-order regularization for sparse digital elevation models. In: International Conference on Scale Space and Variational Methods in Computer Vision. Leibnitz, Austria, pp. 161–173.
  * Mang and Biros (2016) Mang, A., Biros, G., 2016. Constrained ${H^{1}}$-regularization schemes for diffeomorphic image registration. SIAM journal on imaging sciences 9 (3), 1154–1194.
  * Miller et al. (1993) Miller, M. I., Christensen, G. E., Amit, Y., Grenander, U., Dec. 1993. Mathematical textbook of deformable neuroanatomies. Proc. National Academy of Sciences 90 (24), 11944–11948.
  * Miura et al. (2017) Miura, H., Ozawa, S., Nakao, M., Furukawa, K., Doi, Y., Kawabata, H., Kenjou, M., Nagata, Y., 2017. Impact of deformable image registration accuracy on thoracic images with different regularization weight parameter settings. Physica Medica 42, 108–111.
  * Modersitzki (2004) Modersitzki, J., 2004. Numerical methods for image registration. Oxford University Press on Demand, New York, NY.
  * Rohlfing et al. (2003) Rohlfing, T., Maurer, C. R., Bluemke, D. A., Jacobs, M. A., Jun. 2003. Volume-preserving nonrigid registration of MR breast images using free-form deformation with an incompressibility constraint. IEEE Transactions on Medical Imaging 22 (6), 730–741.
  * Rueckert et al. (1999) Rueckert, D., Sonoda, L. I., Hayes, C., Hill, D. L., Leach, M. O., Hawkes, D. J., Aug. 1999. Nonrigid registration using free-form deformations: Application to breast MR images. IEEE Transactions on Medical Imaging 18 (8), 712–721.
  * Schmidt-Richberg et al. (2012) Schmidt-Richberg, A., Werner, R., Handels, H., Ehrhardt, J., 2012. Estimation of slipping organ motion by registration with direction-dependent regularization. Medical image analysis 16 (1), 150–159.
  * Shackleford et al. (2010) Shackleford, J. A., Kandasamy, N., Sharp, G., Oct. 2010. On developing B-spline registration algorithms for multi-core processors. Physics in Medicine & Biology 55 (21), 6329–6351.
  * Shackleford et al. (2012) Shackleford, J. A., Yang, Q., Lourenço, A. M., Shusharina, N., Kandasamy, N., Sharp, G. C., Oct. 2012. Analytic regularization of uniform cubic B-spline deformation fields. In: International Conference on Medical Image Computing and Computer-Assisted Intervention. Nice, France, pp. 122–129.
  * Shusharina and Sharp (2012) Shusharina, N., Sharp, G., Mar. 2012. Analytic regularization for landmark-based image registration. Physics in Medicine & Biology 57 (6), 1477–1498.
  * Sorzano et al. (2005) Sorzano, C. O. S., Thévenaz, P., Unser, M., Apr. 2005. Elastic registration of biological images using vector-spline regularization. IEEE Transactions on Biomedical Engineering 52 (4), 652–663.
  * Thirion (1998) Thirion, J.-P., Sep. 1998. Image matching as a diffusion process: An analogy with Maxwell’s demons. Medical Image Analysis 2 (3), 243–260.
  * Tustison and Avants (2013) Tustison, N., Avants, B., Dec. 2013. Explicit B-spline regularization in diffeomorphic image registration. Frontiers in neuroinformatics 7 (39), 1–13.
  * Vishnevskiy et al. (2016) Vishnevskiy, V., Gass, T., Szekely, G., Tanner, C., Goksel, O., 2016. Isotropic total variation regularization of displacements in parametric image registration. IEEE transactions on medical imaging 36 (2), 385–395.
  * Zhang et al. (2007) Zhang, T., Chi, Y., Meldolesi, E., Yan, D., Jun. 2007. Automatic delineation of on-line head-and-neck computed tomography images: Toward on-line adaptive radiotherapy. International Journal of Radiation Oncology* Biology* Physics 68 (2), 522–530.