Region-specific Diffeomorphic Metric Mapping

Region-specific Diffeomorphic Metric Mapping

Zhengyang Shen
UNC Chapel Hill
zyshen@cs.unc.edu &François-Xavier Vialard
LIGM, UPEM
francois-xavier.vialard@u-pem.fr &Marc Niethammer
UNC Chapel Hill
mn@cs.unc.edu
Abstract

We introduce a region-specific diffeomorphic metric mapping (RDMM) registration approach. RDMM is non-parametric, estimating spatio-temporal velocity fields which parameterize the sought-for spatial transformation. Regularization of these velocity fields is necessary. However, while existing non-parametric registration approaches, e.g., the large displacement diffeomorphic metric mapping (LDDMM) model, use a fixed spatially-invariant regularization our model advects a spatially-varying regularizer with the estimated velocity field, thereby naturally attaching a spatio-temporal regularizer to deforming objects. We explore a family of RDMM registration approaches: 1) a registration model where regions with separate regularizations are pre-defined (e.g., in an atlas space), 2) a registration model where a general spatially-varying regularizer is estimated, and 3) a registration model where the spatially-varying regularizer is obtained via an end-to-end trained deep learning (DL) model. We provide a variational derivation of RDMM, show that the model can assure diffeomorphic transformations in the continuum, and that LDDMM is a particular instance of RDMM. To evaluate RDMM performance we experiment 1) on synthetic 2D data and 2) on two 3D datasets: knee magnetic resonance images (MRIs) of the Osteoarthritis Initiative (OAI) and computed tomography images (CT) of the lung. Results show that our framework achieves state-of-the-art image registration performance, while providing additional information via a learned spatio-temoporal regularizer. Further, our deep learning approach allows for very fast RDMM and LDDMM estimations. Code is available at https://github.com/uncbiag/registration.

 

Region-specific Diffeomorphic Metric Mapping


  Zhengyang Shen UNC Chapel Hill zyshen@cs.unc.edu François-Xavier Vialard LIGM, UPEM francois-xavier.vialard@u-pem.fr Marc Niethammer UNC Chapel Hill mn@cs.unc.edu

\@float

noticebox[b]Preprint. Under review.\end@float

1 Introduction

Quantitative analysis of medical images frequently requires the estimation of spatial correspondences, i.e.image registration. For example, one may be interested in capturing knee cartilage changes over time, localized changes of brain structures, or how organs at risk move between planning and treatment for radiation treatment. Specifically, image registration seeks to estimate the spatial transformation between a source image and a target image, subject to a chosen transformation model.

Transformations can be parameterized via low-dimensional parametric models (e.g., an affine transformation), but more flexible models are required to capture subtle local deformations. Such registration models bajcsy1989multiresolution (); shen2002hammer () may have large numbers of parameters, e.g., a large number of B-spline control points rueckert1999nonrigid () or may even be non-parametric where vector fields are estimated beg2005computing (); modersitzki2004numerical (). Spatial regularity can be achieved by appropriate constraints on displacement fields haber2007image () or by parameterizing the transformation via integration of a sufficiently regular stationary or time-dependent velocity field beg2005computing (); hart2009optimal (); vercauteren2009diffeomorphic (); chen2013large (); wulff2015efficient (). Given sufficient regularization, diffeomorphic transformations can be assured in the continuum. A popular approach based on time-dependent velocity fields is LDDMM beg2005computing (); hart2009optimal (). Optimal LDDMM solutions are geodesics and minimize a geodesic distance. Consequentially, one may directly optimize over a geodesic’s initial conditions in a shooting approach vialard2012diffeomorphic ().

Most existing non-parametric image registration approaches use spatially-invariant regularizers. However, this may not be realistic. E.g., when registering inhale to exhale images of a lung one expects large deformations of the lung, but not of the surrounding tissue. Hence, a spatially-varying regularization would be more appropriate. As the regularizer encodes the deformation model this then allows anticipating different levels of deformation at different image locations.

While spatially-varying regularizers may be used in LDDMM variants schmah2013left () existing approaches do not allow for time-varying spatially-varying regularizers. However, such regularizers would be natural for large displacement as they can move with a deforming image. Hence, we propose a family of registration approaches with spatio-temporal regularizers based on advecting spatially-varying regularizers via an estimated spatio-temporal velocity field. Specifically, we extend LDDMM theory, where the original LDDMM model becomes a special case. In doing so, our entire model, including the spatio-temporal regularizer is expressed via the initial conditions of a partial differential equation. We propose three different approaches based on this model: 1) A model for which the regularizer is specified region-by-region. This would, for example, be natural when registering a labeled atlas image to a target image, as illustrated in Fig. 1. 2) A model in which the regularizer is estimated jointly with the spatio-temporal velocity fields. 3) A deep learning model which predicts the regularizer and the initial velocity field, thereby resulting in a very fast registration approach.

Contributions: 1) We propose RDMM, a new registration model for large diffeomorphic deformations with a spatio-temporal regularizer capable of following deforming objects and thereby providing a more natural representation of deformations than existing non-parametric models, such as LDDMM. 2) Via a variational formulation we derive shooting equations that allow specifying RDMM solutions entirely based on their initial conditions: an initial momentum field and an initial spatially-varying regularizer. 3) We prove that diffemorphisms can be obtained for RDMM in the continuum for sufficiently regular regularizers. 4) We explore an entire new family of registration models based on RDMM and provide optimization-based and very fast deep-learning-based approaches to estimate the initial conditions of these registration models. 5) We demonstrate the utility of our approach via experiments on synthetic data and on two 3D medical image datasets.

Figure 1: RDMM registration example. The goal is to register the dark blue area with high fidelity (i.e., allowing large local deformations), while assuring small deformations in the cyan area. Initially (t=0), a spatially-varying regularizer (fifth column) is specified in the source image space, where dark blue indicates small regularization and yellow large regularization. Specifically, regularizer values indicate effective local standard deviations of a local multi-Gaussian regularizer. Since the transformation map and the regularizer are both advected according to the estimated velocity field, the shape of the regularizer follows the shape of the deforming dark blue region and is of the same shape as the region of interest in the target space at the final time (t=1) (as can be seen in the second and the last columns). Furthermore, objects inside the dark blue region are indeed aligned well, whereas objects in the cyan region were not strongly deformed due to the larger regularization there.

2 Standard LDDMM Model

LDDMM (beg2005computing, ) is a non-parametric registration approach based on principles from fluid mechanics. It is based on the estimation of a spatio-temporal velocity field from which the sought-for spatial transformation can be computed via integration of For appropriately regularized velocity fields (dupuis1998variational, ), diffeomorphic transformations can be guaranteed. The optimization problem underlying LDDMM for images can be written as ( being the gradient)

(2.1)

Here, the goal is to register the source image to the target image in unit time. is a similarity measure between images, often sum of squared differences, normalized cross correlation, or mutual information. Furthermore, we note that , where denotes the inverse of in the target image space. The evolution of this map can be expressed as

(2.2)

where denotes the Jacobian. Equivalently, in Eq. (2.1), we directly advect the image (hart2009optimal, ; vialard2012diffeomorphic, ) via . To assure smooth transformations, LDDMM penalizes non-smooth velocity fields via the norm , where indicates the inner product and a differential operator.

At optimality of Eq. (2.1) the Euler-Lagrange equations are (div denoting the divergence)

(2.3)

Here, is the adjoint variable to , also known as the scalar momentum (hart2009optimal, ; vialard2012diffeomorphic, ). As is a smoother it is often chosen as a convolution, i.e., . Note that is the vector-valued momentum and thus . One can directly optimize over solving Eq. (2.3(beg2005computing, ) or regard Eq. (2.3) as a constraint defining a geodesic path (vialard2012diffeomorphic, ) and optimize over all such solutions subject to a penalty on the initial scalar momentum as well as the similarity measure. Alternatively, one can express (see suppl. material) these equations entirely with respect to the vector-valued momentum, , resulting in the Euler-Poincaré equation for diffeomorphisms (EPDiff) (younes2009evolutions, ):

(2.4)

which defines the evolution of the spatio-temporal velocity field based on the initial condition, , of the momentum, from which the transformation can be computed using Eq. (2.2). Both (2.3) and (2.4) can be used to implement shooting-based LDDMM (vialard2012diffeomorphic, ; singh2013vector, ). As LDDMM preserves the momentum, is constant over time and hence a shooting formulation can be written as

(2.5)

subject to the EPDiff equation (2.4) including the advection of , where .

A shooting-formulation has multiple benefits: 1) it allows for a compact representation of via its initial conditions; 2) as the optimization is w.r.t.the initial conditions, a solution is a geodesic by construction; 3) these initial conditions can be predicted via deep-learning, resulting in very fast registration algorithms which inherit the theoretical properties of LDDMM (yang2017quicksilver, ; yang2016fast, ). We therefore use this formulation as the starting point for RDMM in Sec 3.

3 Region-Specific Diffeomorphic Metric Mapping (RDMM)

In standard LDDMM approaches, the regularizer (equivalently the kernel ) is spatially invariant. While recent work introduced spatially-varying metrics in LDDMM, for stationary velocity fields, or for displacement-regularized registration (risser2013piecewise, ; schmah2013left, ; niethammer2019_cvpr, ; pace2013locally, ), all of these approaches use a temporally fixed regularizer. Hereafter, we generalize LDDMM by advecting a spatially-varying regularizer via the estimated spatio-temporal velocity field. Standard LDDMM is a special case of our model.

Following (niethammer2019_cvpr, ), we introduce a finite family of reproducing kernel Hilbert spaces (RKHS) which are defined by the pre-defined Gaussian kernels with . We use a partition of unity , on the image domain. As we want the kernel to be active on the region determined by we introduce the vector field for . On this new space of vector fields, there exists a natural RKHS structure defined by

(3.1)

whose kernel is . Thus the velocity reads (see suppl. material for derivation)

(3.2)

which can capture multi-scale aspects of deformation risser2010simultaneous (). In LDDMM, the weights are constant and pre-defined. Here, we allow spatially-dependent weights . In particular (see formulation below), we advect them via, thereby making them spatio-temporal, i.e., . In this setup, weights only need to be specified at initial time . As the Gaussian kernels are fixed convolutions can still be efficiently computed in the Fourier domain.

We prove (see suppl. material) that, for sufficiently smooth weights , the velocity field is bounded and its flow is a diffeomorphism. Following (niethammer2019_cvpr, ), to assure the smoothness of the initial weights we instead optimize over initial pre-weights, , s.t. , where is a fixed Gaussian with a small standard deviation, . In addition, we constrain to locally sum to one. The optimization problem for our RDMM model then becomes

(3.3)

subject to the constraints

(3.4)

As for LDDMM, we can compute the optimality conditions for Eq. (3.3) which we use for shooting.

Theorem 1 (Image-based RDMM optimality conditions).

With the adjoints (for ) and (for ) and the momentum the optimality conditions for (3.3) are:

(3.5)
(3.6)
(3.7)

subject to

(3.8)
Theorem 2 (Momentum-based RDMM optimality conditions).

The RDMM optimality conditions of Thm. (1) can be written entirely w.r.t.  the momentum (as defined in Thm (1)). They are:

(3.9)
(3.10)

where and subject to the constraints of Eq. (3.8) which define the relationship between the velocity and the momentum.

For spatially constant pre-weights, we recover EPDiff from the momentum-based RDMM optimality conditions. Instead of advecting via , we can alternatively advect the pre-weights directly, as . For the image-based and the momentum-based formulations, the velocity field is obtained by smoothing the momentum, .

Regularization of the Regularizer

For a well-posed optimization problem, we need to regularize the regularizer to encourage smooth solutions. We follow (niethammer2019_cvpr, ) which uses a simple optimal mass transport (OMT) penalty on the weights. Specifically, we locally penalize , where is a chosen power. As the regularization only affects the initial conditions for the pre-weights, the evolution equations for the optimality conditions (i.e., the modified EPDiff equation) do not change. Additional regularizers, such as total variation terms as proposed in (niethammer2019_cvpr, ), are possible and easy to integrate into our RDMM framework as they only affect initial conditions. For simplicity, we focus on regularizing via OMT.

Shooting Formulation

As energy is conserved (see suppl. material) the momentum-based shooting formulation becomes

(3.11)

subject to the evolution equations of Thm. 2. Similarly, the shooting formulation can use the image-based evolution equations of Thm. 1 where optimization would be over instead of .

Figure 2: Illustration of the learning framework (left) and its non-parametric registration component (right). A multi-step affine-network first predicts the affine transformation map shen2019networks () followed by an iterable non-parametric registration to estimate the final transformation map. The LDDMM component uses one network to generate the initial momentum. RDMM also uses a second network to predict the initial regularizer pre-weights. We integrate the RDMM evolution equations at low-resolution (based on the predicted initial conditions) to save memory. The final transformation map is obtained via upsampling. See suppl. material for the detailed structure of the LDDMM/RDMM units.

4 Learning Framework

The parameters for RDMM, i.e., the initial momentum and the initial pre-weights, can be obtained by numerical optimization, either over the momentum and the pre-weights or only over the momentum if the pre-weights are prescribed. Just as for the LDDMM model, such a numerical optimization is computationally costly. Consequentially, various deep-learning (DL) approaches have been proposed to instead predict displacements (balakrishnan2018unsupervised, ; cao2018deformable, ), stationary velocity (rohe2017svf, ) or momentum fields (yang2017quicksilver, ; niethammer2019_cvpr, ). Supervised yang2017quicksilver (); yang2016fast () and unsupervised jaderberg2015spatial (); de2017end (); li2017non (); balakrishnan2018unsupervised (); dalca2018unsupervised () DL registration approaches exist. All of them are fast as only a regression solution needs to be evaluated at test time and no further numerical optimization is necessary. Additionally, such DL models benefit from learning over an entire population instead of relying only on information from given image-pairs.

We use two different deep learning approaches, illustrated in Fig. 2, to predict 1) the LDDMM initial momentum only (as a special case of RDMM) and 2) the initial momentum and the pre-weights for RDMM. All these networks are trained end-to-end and are initialized by an affine transformation obtained via an affine prediction-network, which accounts for large, global displacements, whereas the LDDMM/RDMM parts account for localized deformations. Overall, we use two networks for LDDMM prediction and three networks for RDMM prediction. The first network is a multi-step affine network to predict the affine transformation following (shen2019networks, ). The second network predicts the initial momentum, . The third network predicts the initial pre-weights, , defining the regularizer.

We use low-resolution maps and map compositions for the momentum and the pre-weight networks. This reduces computational cost significantly. The final transformation map is obtained via upsampling, which is reasonable as we assume smooth transformations. We use 3D UNets (cciccek2016_3d_unet, ) for momentum and pre-weight prediction. Both the affine and the non-parametric networks can be iterated to refine the prediction results: i.e.  the input source image and the initial map are replaced with the currently warped image and the transformation map respectively for the next iteration.

Inverse Consistency: For the DL approaches we follow shen2019networks () and compute bidirectional (source to target denoted as and target to source denoted as ) registration losses and an additional symmetry loss, , where refers to the identity map. This encourages symmetric consistency.

5 Experimental Results and Setup

Datasets: To demonstrate the behavior of RDMM, we evaluate the model on three datasets: 1) a synthetic dataset for illustration, 2) a 3D computed tomography dataset (CT) of a lung, and 3) a large 3D magnetic resonance imaging (MRI) dataset of the knee from the Osteoarthritis Initiative (OAI).
The synthetic dataset consists of three types of shapes (rectangles, triangles, and ellipses). There is one foreground object in each image with two objects inside and at most five objects outside. Each source image object has a counterpart in the target image; the shift, scale, and rotations are random. We generated 40 image pairs of size for evaluation. Fig. 1 shows example synthetic images.
The lung dataset consists of 49 inspiration/expiration image pairs with lung segmentations. Each image is of size . We register from the expiration phase to the inspiration phase for all 49 pairs.
The OAI dataset consists of 176 manually labeled MRI from 88 patients (2 longitudinal scans per patient) and 22,950 unlabeled MR images from 2,444 patients. Labels are available for femoral and tibial cartilage. We divide the patients into a training group and a testing group, where 300 cross-subject pairs are evaluated with the same evaluation settings as in shen2019networks ().

Deformation models: Affine registration is performed before each LDDMM/RDMM registration.

Affine model: We implemented a multi-scale affine model solved via numerical optimization and a multi-step deep neural network to predict the affine transformation parameters.

Family of non-parametric models: We implemented both optimization and deep-learning versions of a family of non-parametric registration methods: a vector-momentum based stationary velocity field model (vSVF) (, ), LDDMM (, ), and RDMM (, ). We use the dopri5 solver using the adjoint sensitivity method (chen2018neural, ) to integrate the evolution equations in time. For solutions based on numerical optimization, we use a multi-scale strategy with L-BGFS liu1989limited () as the optimizer. For the deep learning models, we compute solutions for a low-resolution map (factor of 0.5) which is then upsampled. We use Adam kingma2014adam () for optimization.

Image similarity measure: We use multi-kernel Localized Normalized Cross Correlation (mk-LNCC) shen2019networks (). mk-LNCC computes localized normalized cross correlation (NCC) with different window sizes and combines these measures via a weighted sum.

Weight visualization: To illustrate the behavior of the RDMM model, we visualize the estimated standard deviations, i.e.the square root of the local variance .

Estimation approaches: To illustrate different aspects of our approach we perform three types of RDMM registrations: 1) registration with a pre-defined regularizer (Sec. 5.1), 2) registration with simultaneous optimization of the regularizer, via optimization of the initial momentum and pre-weights (Sec. 5.2), and 3) registration via deep learning predicting the initial momentum and regularizer pre-weights (Sec. 5.3). Detailed settings for all approaches are in the suppl. material.

5.1 Registration with a pre-defined regularizer

To illustrate the base capabilities of our models, we prescribe an initial spatially-varying regularizer in the source image space. We show experimental results for pair-wise registration of the synthetic data as well as for the 3D lung volumes.

Fig. 1 shows the registration result for an example synthetic image pair. We use small regularization in the blue area and large regularization in the surrounding area. As expected, most of the deformations occur inside the blue area as the regularizer is more permissive there. We also observe that the regularizer is indeed advected with the image. For the real lung image data we use a small regularizer inside the lung (as specified by the given lung mask) and a large regularizer in the surrounding tissue. Fig. 3 shows that most of the deformations are indeed inside the lung area while the deformation outside the lung is highly regularized as desired. We evaluate the Dice score between the warped lung and the target lung, achieving on average (over all inhalation/exhalation pairs). Fig. 3 also shows the determinant of the Jacobian of the transformation map (defined in target space): the lung region shows small values (illustrating expansion) while other region are either volume preserved (close to 1) or compressed (bigger than 1). Overall the deformations are smooth.

Figure 3: RDMM lung registration result with a pre-defined regularizer. Lung images at expiration (source) are registered to the corresponding inspiration (target) images. Last column: Inside the lung a regularizer with small standard deviation ( , ) and outside the lung with large standard deviation is used (, ). Deformations are largely inside the lung, the surrounding tissue is well regularized as expected. Columns 1 to 3 show results in image space while columns 4 to 6 refer to the results in label space. The second to last column shows the determinant of the Jacobian of the spatial transformation, .

5.2 Registration with an optimized regularizer

Figure 4: Illustration of the RDMM registration results with an optimized regularizer on the synthetic dataset. All objects are warped from the source image space to the target image space. The last two columns show the regularizer () at and respectively.
Figure 5: Average Dice scores for all objects. Left to right: SVF, LDDMM and RDMM.

In contrast to the experiments in Sec. 5.1, we jointly estimate the initial momentum and the initial regularizer pre-weights for pairwise registration. We use the synthetic data and the knee MRIs.

Fig. 5 shows that the registration warped every object in the source image to the corresponding object in the target image. The visualization of the initial regularizer shows that low regularization is assigned close to object edges making them locally deformable. The visualizations of the regularizer at the initial time point () and the final time point () show that it deforms with the image. That low regularization values are localized is sensible as the OMT regularizer prefers spatially sparse solutions (in the sense of sparsely assigning low levels of regularity). If the desired deformation model is piecewise constant our RDMM model could readily be combined with a total variation penalty as in (niethammer2019_cvpr, ). Fig. 5 compares average Dice scores for all objects for vSVF, LDDMM and RDMM separately. They all achieve high and comparable performance indicating good registration quality. But only RDMM provides additional information about local regularity.

We further evaluate RDMM on 300 images pairs from the OAI dataset. The optimization methods section of Tab. 6 compares registration performance for different optimization-based algorithms. RDMM achieves the best Dice score. While RDMM is in theory diffeomorphic, we observe some foldings, whereas no such foldings appear for LDDMM and SVF. This is likely due to inaccuracies when discretizing the evolution equations and when discretizing the determinant of the Jacobian. Further, RDMM may locally exhibit stronger deformations than LDDMM or vSVF, especially when the local regularization (via OMT) is small, making it numerically more challenging. Most of the folds appear at the boundary (and are hence likely due to boundary discretization artifacts) or due to anatomical inconsistency in the source and target images, where large deformations may be estimated.

Method OAI Dataset Dice Folds Time (s) ——— Affine Methods ———- ——– affine-NiftyReg 30.43 (12.11) 0 45 affine-opt 34.49 (18.07) 0 8 affine-net 44.58 (7.74) 0 0.20 ———– Optimization Methods ——- ———- Demonsvercauteren2009diffeomorphic (); vercauteren2008symmetric () 63.47 (9.52) 0.56 114 SyNavants2009advanced (); avants2008symmetric () 65.71 (15.01) 0 1330 NiftyReg-NMIourselin2001reconstructing (); modat2014global (); rueckert1999nonrigid (); modat2010fast () 59.65 (7.62) 0 143 NiftyReg-LNCC 67.92 (5.24) 35.19 270 vSVF-opt 67.35 (9.73) 0 79 LDDMM-opt 67.72(8.94) 0 457 RDMM-opt 68.18(8.36) 17.37 627 ———- Learning-based Methods ——- —— —- VoxelMorphdalca2018unsupervised ()(with aff) 66.08 (5.13) 3.31 0.31 vSVF-net shen2019networks () 67.59 (4.47) 0.39 0.62 LDDMM-net 67.63(4.51) 0 0.85 RDMM-net 67.94(4.40) 0.47 1.1
Figure 6: Comparison of registration methods for cross-subject registrations on the OAI dataset based on Dice scores. -opt and -net refer to optimization- and DL-based methods respectively. For all DL methods, we report performance after two-step refinement. Folds refers to the absolute value of the sum of the determinant of the Jacobian in the folding area (i.e., where the determinant of the Jacobian is negative); Time refers to the average registration time for a single image pair.
Figure 7: Illustration of RDMM registration results on the OAI dataset. "* s" and "* b" refer to the regularizer with in set to and respectively; "learn *" and "opt *" refer to a learnt regularizer and an optimized one respectively; "Jaco" refers to the absolute value of the determinant of the Jacobian; "init_w" refers to the initial weight map of the regularizer (as visualized via ). The first four columns refer to registration results in image space.

5.3 Registration with a learnt regularizer via deep learning

Finally, we evaluate our learning framework for non-parametric registration approaches on the OAI dataset. For vSVF, we follow (shen2019networks, ) to predict the momentum. We implement the same approach for LDDMM, where the vSVF inference unit is replaced by one for LDDMM (i.e., instead of advecting via a stationary velocity field we integrate EPDiff). Similarly, for RDMM, the inference unit is replaced by the evolution equation for RDMM and we use an additional network to predict the regularizer pre-weights. Fig. 6 shows the results. The non-parametric DL models achieve comparable performance to their optimization counterparts, but are much faster, while learning-based RDMM simultaneously predicts the initial regularizer which captures aspects of the knee anatomy. Fig. 7 shows the determinant of the Jacobian of the transformation map. It is overall smooth and folds are much less frequent than for the corresponding optimization approach, because the DL model penalizes transformation asymmetry. Fig. 7 also clearly illustrates the benefit of training the DL model based on a large image dataset: compared with the optimization approach (which only works on individual image pairs), the initial regularizer predicted by the deep network captures anatomically meaningful information much better: the bone (femur and tibia) and the surrounding tissue show large regularity.

6 Conclusion and Future Work

We introduced RDMM, a generalization of LDDMM registration which allows for spatially-varying regularizers advected with a deforming image. In RDMM, both the estimated velocity field and the estimated regularizer are time- and spatially-varying. We used a variational approach to derive shooting equations which generalize EPDiff and allow the parameterization of RDMM using only the initial momentum and regularizer pre-weights. We also prove that diffeomorphic transformation can be obtained for RDMM with sufficiently regular regularizers. Experiments with pre-defined, optimized, and learnt regularizers show that RDMM is flexible and its solutions can be estimated quickly via deep learning. Future work could focus on numerical aspects and explore different initial constraints, such as total-variation constraints, depending on the desired deformation model. Indeed, a promising avenue of research consists in learning regularizers which include more physical/mechanical a-priori information in the deformation model. For instance, a possible first step in this direction consists in parameterizing non-isotropic kernels to favor deformations in particular directions.
Acknowledgements Research reported in this work was supported by the National Institutes of Health (NIH) and the National Science Foundation (NSF) under award numbers NSF EECS-1711776 and NIH 1R01AR072013. The content is solely the responsibility of the authors and does not necessarily represent the official views of the NIH or the NSF. We would also like to thank Dr. Raúl San José Estépar for providing the lung data.

References

  • [1] B. B. Avants, C. L. Epstein, M. Grossman, and J. C. Gee. Symmetric diffeomorphic image registration with cross-correlation: evaluating automated labeling of elderly and neurodegenerative brain. Medical image analysis, 12(1):26–41, 2008.
  • [2] B. B. Avants, N. Tustison, and G. Song. Advanced normalization tools (ANTS). Insight j, 2:1–35, 2009.
  • [3] R. Bajcsy and S. Kovačič. Multiresolution elastic matching. CVGIP, 46(1):1–21, 1989.
  • [4] G. Balakrishnan, A. Zhao, M. R. Sabuncu, J. Guttag, and A. V. Dalca. An unsupervised learning model for deformable medical image registration. In CVPR, pages 9252–9260, 2018.
  • [5] M. F. Beg, M. I. Miller, A. Trouvé, and L. Younes. Computing large deformation metric mappings via geodesic flows of diffeomorphisms. IJCV, 61(2):139–157, 2005.
  • [6] X. Cao, J. Yang, J. Zhang, Q. Wang, P.-T. Yap, and D. Shen. Deformable image registration using a cue-aware deep regression network. IEEE Transactions on Biomedical Engineering, 65(9):1900–1911, 2018.
  • [7] T. Q. Chen, Y. Rubanova, J. Bettencourt, and D. K. Duvenaud. Neural ordinary differential equations. In Advances in Neural Information Processing Systems, pages 6571–6583, 2018.
  • [8] Z. Chen, H. Jin, Z. Lin, S. Cohen, and Y. Wu. Large displacement optical flow from nearest neighbor fields. In CVPR, pages 2443–2450, 2013.
  • [9] Ö. Çiçek, A. Abdulkadir, S. S. Lienkamp, T. Brox, and O. Ronneberger. 3D U-Net: learning dense volumetric segmentation from sparse annotation. In International conference on medical image computing and computer-assisted intervention, pages 424–432. Springer, 2016.
  • [10] A. V. Dalca, G. Balakrishnan, J. Guttag, and M. R. Sabuncu. Unsupervised learning for fast probabilistic diffeomorphic registration. arXiv preprint arXiv:1805.04605, 2018.
  • [11] B. D. de Vos, F. F. Berendsen, M. A. Viergever, M. Staring, and I. Išgum. End-to-end unsupervised deformable image registration with a convolutional neural network. In MLCDS, pages 204–212. Springer, 2017.
  • [12] J. Delon, J. Salomon, and A. Sobolevski. Local matching indicators for transport problems with concave costs. SIAM Journal on Discrete Mathematics, 26(2):801–827, 2012.
  • [13] P. Dupuis, U. Grenander, and M. I. Miller. Variational problems on flows of diffeomorphisms for image matching. Quarterly of applied mathematics, pages 587–600, 1998.
  • [14] E. Haber and J. Modersitzki. Image registration with guaranteed displacement regularity. International Journal of Computer Vision, 71(3):361–372, 2007.
  • [15] G. L. Hart, C. Zach, and M. Niethammer. An optimal control approach for deformable registration. In CVPR, pages 9–16. IEEE, 2009.
  • [16] M. Jaderberg, K. Simonyan, A. Zisserman, et al. Spatial transformer networks. In NIPS, pages 2017–2025, 2015.
  • [17] D. P. Kingma and J. Ba. Adam: A method for stochastic optimization. arXiv preprint arXiv:1412.6980, 2014.
  • [18] H. Li and Y. Fan. Non-rigid image registration using fully convolutional networks with deep self-supervision. arXiv preprint arXiv:1709.00799, 2017.
  • [19] D. C. Liu and J. Nocedal. On the limited memory BFGS method for large scale optimization. Mathematical programming, 45(1-3):503–528, 1989.
  • [20] M. Modat, D. M. Cash, P. Daga, G. P. Winston, J. S. Duncan, and S. Ourselin. Global image registration using a symmetric block-matching approach. Journal of Medical Imaging, 1(2):024003, 2014.
  • [21] M. Modat, G. R. Ridgway, Z. A. Taylor, M. Lehmann, J. Barnes, D. J. Hawkes, N. C. Fox, and S. Ourselin. Fast free-form deformation using graphics processing units. Computer methods and programs in biomedicine, 98(3):278–284, 2010.
  • [22] J. Modersitzki. Numerical methods for image registration. Oxford University Press on Demand, 2004.
  • [23] M. Niethammer, R. Kwitt, and F.-X. Vialard. Metric learning for image registration. arXiv preprint arXiv:1904.09524, 2019.
  • [24] J. Nocedal and S. Wright. Numerical optimization. Springer Science & Business Media, 2006.
  • [25] S. Ourselin, A. Roche, G. Subsol, X. Pennec, and N. Ayache. Reconstructing a 3D structure from serial histological sections. Image and vision computing, 19(1-2):25–31, 2001.
  • [26] D. F. Pace, S. R. Aylward, and M. Niethammer. A locally adaptive regularization based on anisotropic diffusion for deformable image registration of sliding organs. IEEE transactions on medical imaging, 32(11):2114–2126, 2013.
  • [27] L. Risser, F.-X. Vialard, H. Y. Baluwala, and J. A. Schnabel. Piecewise-diffeomorphic image registration: Application to the motion estimation between 3D CT lung images with sliding conditions. Medical image analysis, 17(2):182–193, 2013.
  • [28] L. Risser, F.-X. Vialard, R. Wolz, D. D. Holm, and D. Rueckert. Simultaneous fine and coarse diffeomorphic registration: application to atrophy measurement in Alzheimer’s disease. In MICCAI, pages 610–617. Springer, 2010.
  • [29] M.-M. Rohé, M. Datar, T. Heimann, M. Sermesant, and X. Pennec. SVF-Net: Learning deformable image registration using shape matching. In MICCAI, pages 266–274. Springer, 2017.
  • [30] D. Rueckert, L. I. Sonoda, C. Hayes, D. L. Hill, M. O. Leach, and D. J. Hawkes. Nonrigid registration using free-form deformations: application to breast MR images. TMI, 18(8):712–721, 1999.
  • [31] T. Schmah, L. Risser, and F.-X. Vialard. Left-invariant metrics for diffeomorphic image registration with spatially-varying regularisation. In International Conference on Medical Image Computing and Computer-Assisted Intervention, pages 203–210. Springer, 2013.
  • [32] D. Shen and C. Davatzikos. HAMMER: hierarchical attribute matching mechanism for elastic registration. TMI, 21(11):1421–1439, 2002.
  • [33] Z. Shen, X. Han, Z. Xu, and M. Niethammer. Networks for joint affine and non-parametric image registration. arXiv preprint arXiv:1903.08811, 2019.
  • [34] N. Singh, J. Hinkle, S. Joshi, and P. T. Fletcher. A vector momenta formulation of diffeomorphisms for improved geodesic regression and atlas construction. In 2013 IEEE 10th International Symposium on Biomedical Imaging, pages 1219–1222. IEEE, 2013.
  • [35] T. Vercauteren, X. Pennec, A. Perchant, and N. Ayache. Symmetric log-domain diffeomorphic registration: A demons-based approach. In MICCAI, pages 754–761. Springer, 2008.
  • [36] T. Vercauteren, X. Pennec, A. Perchant, and N. Ayache. Diffeomorphic demons: Efficient non-parametric image registration. NeuroImage, 45(1):S61–S72, 2009.
  • [37] F.-X. Vialard, L. Risser, D. Rueckert, and C. J. Cotter. Diffeomorphic 3D image registration via geodesic shooting using an efficient adjoint calculation. International Journal of Computer Vision, 97(2):229–241, 2012.
  • [38] J. Wulff and M. J. Black. Efficient sparse-to-dense optical flow estimation using a learned basis and layers. In CVPR, pages 120–130, 2015.
  • [39] X. Yang, R. Kwitt, and M. Niethammer. Fast predictive image registration. In DLMIA, pages 48–57. Springer, 2016.
  • [40] X. Yang, R. Kwitt, M. Styner, and M. Niethammer. Quicksilver: Fast predictive image registration–a deep learning approach. NeuroImage, 158:378–396, 2017.
  • [41] L. Younes, F. Arrate, and M. I. Miller. Evolutions equations in computational anatomy. NeuroImage, 45(1):S40–S50, 2009.

Appendix A Supplementary Material

This supplementary material provides additional details illustrating the proposed approach. We start by deriving the form of the smoothing kernel for the RDMM model in Sec. A.1. Based on this kernel form we can then detail the derivation of the RDMM optimality conditions in Sec. A.2. In Sec. A.3, we prove that the regularization energy is conserved over time, which allows formulating our RDMM shooting strategy based on initial conditions only. Sec. A.4 details the good theoretical behavior of our model. Sec. A.5 describes the optimization/training strategy with regard to the the initial pre-weight regularization. Sec. A.6 visualizes the inference process of the LDDMM/RDMM method. Sec. A.7 analyzes the behavior of the OMT term. Lastly, Sec. A.8 details the settings of our experiments.

a.1 Variational Derivation of the Smoothing Kernel

The derivation of our RDMM model makes use of a smoothing kernel of the form . This kernel form is a direct consequence of the definition of our variational definition of the smoothing kernel.

Recall that similar to (niethammer2019_cvpr, ) we define

(A.1)

for a given velocity field . To compute an explicit form of the norm we need to solve the constrained optimization problem of this definition. Specifically, we introduce the vector-valued Lagrange multiplier . Thus the Lagrangian, , becomes111We multiply the objective function by for convenience. This does not change the solution.

(A.2)

The variation of the Lagrangian is

(A.3)

By collecting all the terms, the optimality conditions (i.e., where the variation vanishes) are

(A.4)

Hence, we can write the norm in the following form :

(A.5)

Consequentially, the associated kernel is . Assuming the kernel can be written as a convolution, we can therefore express the velocity as:

(A.6)

a.2 Optimality Conditions

In this section we derive the RDMM optimality conditions. Both for the image-based and the momentum-based cases. Recall that the overall registration energy of RDMM can be written as:

(A.7)

under the constraints222In this section, to simplify the notation, we denote the partial derivative by only the subscript .

(A.8)

Proof of Thm. (1)
We compute the variations of the Lagrangian, to the energy (i.e., where constraints are added via Lagrangian multipliers) with respect to , , , and :

(A.9)

We use

(A.10)

According to Green’s theorem and assuming vanishes on the boundary, we get

(A.11)

Similarly, we have

(A.12)
(A.13)

Now, Eq. (A.9) reads

(A.14)

We first collect and to obtain the final condition on and the initial condition on :

(A.15)

Next, we work on . Remember, we have , where , thus

(A.16)

Note that for radially symmetric kernels (such as Gaussian kernels) ,

(A.17)

Thus, we can get

(A.18)

Substituting Eq. (A.18) into Eq. (A.16), we get

(A.19)

Next, we decompose the terms. We define the momentum, .

Now, we can collect the variation for and for and obtain the optimality conditions

(A.20)
(A.21)

Finally, we get the optimality conditions for image-based RDMMM derived from Eq. (A.7) and Eq. (A.8):

(A.22)

where and .

Proof of Thm. (2) We now derive the optimality conditions for the momentum-based formulation of RDMM. We start by taking the time derivative of the momentum and obtain

(A.23)
(A.24)

By substituting the time derivatives , , , and from Eq. (A.22) we obtain

(A.25)

Using the following two relations,

(A.26)

where denotes the Jacobian and the Hessian, we can rewrite Eq. (A.25) as

(A.27)
(A.28)
(A.29)
(A.30)

Noticing that

(A.31)

we can write

(A.32)
(A.33)

Finally, we get

(A.34)

which gives the result.

a.3 Energy Conservation

To formulate a shooting-based solution we would like to avoid integrating over time. We here show that this quantity is conserved. Hence, , which allows us to write our shooting equations only with respect to initial conditions subject to the momentum-based evolution euqations of RDMM.

Recall that the energy is preserved by the EPDiff equation since it can alternatively be written as

(A.35)

where is the adjoint of . It implies that

since . In fact, there is more than conservation of the energy, since the momentum is actually advected along the flow. Now, formula (3.9) can be shortened as and it implies that, denoting the kernel to shorten the notations,

since the first term vanishes as for the standard EPDiff equation and the two other terms cancel each other since and . Here we assumed the kernel to be symmetric in writing this equation but the result holds in general, the equations being simply modified with the transpose kernel.

a.4 Mathematical properties

In this section, we prove that given , there exists a solution solving Equations (A.7) and (A.8) at least until a time which could be less than . The notations or denote the sup norm of maps.

Theorem 3.

Let and suppose for every , . Given initial weights and time dependent vector fields , there exists a unique solution to Equations (3) until time .

Proof.

The first step consists in proving that there exists a solution locally in time. To this end, the proof follows a fixed point argument on the space