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

Keyur D. Shah¹, James A. Shackleford¹, Nagarajan Kandasamy¹, Gregory C. Sharp² ¹ Electrical and Computer Engineering Department, Drexel University, Philadelphia, PA 19104, USA ² Department of Radiation Oncology, Massachusetts General Hospital, Boston, MA 02114, USA [email protected]

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.

Refer to caption — 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}$ mm³. 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.

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

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.

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

(a)

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