Symplectomorphic registration with phase space regularization by entropy spectrum pathways

TL;DR

This paper introduces a flexible, robust image registration method using symplectomorphic maps and phase space regularization with entropy spectrum pathways, demonstrated on various MRI modalities with efficient processing times.

Contribution

It presents a novel registration framework combining Hamiltonian formalism, symplectomorphic maps, and ESP-based regularization, enhancing flexibility and robustness over existing methods.

Findings

Effective registration across multiple MRI modalities.

Processing times of less than a few minutes for high-resolution images.

Demonstrated robustness and flexibility in neuroimaging applications.

Abstract

The ability to register image data to a common coordinate system is a critical feature of virtually all imaging studies that require multiple subject analysis, combining single subject data from multiple modalities, or both. However, in spite of the abundance of literature on the subject and the existence of several variants of registration algorithms, their practical utility remains problematic, as commonly acknowledged even by developers of these methods because the complexity of the problem has resisted a general, flexible, and robust theoretical and computational framework. To address this issue, we present a new registration method that is similar in spirit to the current state-of-the-art technique of diffeomorphic mapping, but is more general and flexible. The method utilizes a Hamiltonian formalism and constructs registration as a sequence of symplectomorphic maps in conjunction with a novel phase space regularization based on the powerful entropy spectrum pathways (ESP) framework. The method is demonstrated on the three different magnetic resonance imaging (MRI) modalities routinely used for human neuroimaging applications by mapping between high resolution anatomical (HRA) volumes, medium resolution diffusion weighted MRI (DW-MRI) and HRA volumes, and low resolution functional MRI (fMRI) and HRA volumes. The typical processing time for high quality mapping ranges from less than a minute to several minutes on a modern multi core CPU for typical high resolution anatomical (~256x256x256 voxels) MRI volumes.

Figures

Tables

Table 1. Table 1: Comparison of accuracy and efficiency of phantom registration between three commonly used packages (ANTs, FSL and AFNI) and our symplectomorphic registration (SYM-REG) approach. Simple Root-Mean-Square Deviation (RMSD) between reference and registered images are used as an accuracy metric and time is reported to characterize an efficiency. One of the center slices of registered images is also included to provide visual insights into correspondent registration artifacts.

	ANTsa,c	ANTsb,c	FSLd	AFNIe	SYM-REG
RMSD	156.4	81.0	448.2	234.3	54.0
Time	10min 15sec	1h 21min	12min 43sec	45min 50sec	38sec

Table 2. Table 2: Comparison of accuracy and efficiency of HRA registration between ANTs SyN, AFNI and our symplectomorphic registration (SYM-REG) approach. Five different subjects and five different warpage types were used for the comparison. Two subtables include mean RMSD values and mean execution times calculated across subject (top) and across warpage types (bottom) with the best values marked in red. The average of all subjects and warpage types gives RMSD of 304.6, 348.5 and 584.0 for SYM-REG, AFNI and ANTs respectively with 976, 3414 and 1223 seconds required to finish the processing for each method.

Method	Mean RMSD / time (s) across subjects
AFNI 3dQwarp	402.2/130.8	363.8/139.4	371.6/130.0	316.5/137.8	338.1/144.8
ANTs SyN	578.2/38.2	606.0/47.0	586.2/43.2	581.1/60.0	568.5/56.2
SYM-REG	298.2/54.8	312.4/37.0	311.8/34.2	300.7/33.4	299.8/35.8
	Mean RMSD / time (s) across warpage types
AFNI 3dQwarp	162.9/134.0	373.5/137.2	482.7/135.6	404.6/130.0	368.6/146.0
ANTs SyN	140.8/50.6	705.0/47.4	602.9/53.6	781.9/51.0	689.4/42.0
SYM-REG	87.0/38.4	331.6/46.8	438.7/36.2	351.8/38.6	313.9/35.2

Table 3. Table 3: AFNI 3dQwarp RMSD and time (s) for registration of five different subjects and five different warpage types

	1	2	3	4	5
1	225.2192 / 129	161.0997 / 130	185.5213 / 123	107.5805 / 132	134.8321 / 156
2	411.3067 / 136	380.6455 / 126	380.4736 / 123	335.5967 / 157	359.5755 / 144
3	520.8890 / 119	483.1769 / 167	497.8304 / 128	452.3812 / 123	459.2277 / 141
4	445.6333 / 117	417.0812 / 132	414.8282 / 125	363.8450 / 145	381.5591 / 131
5	408.1614 / 153	376.9296 / 142	379.3344 / 151	323.2617 / 132	355.3540 / 152

Table 4. Table 4: ANTs SyN RMSD and time (s) for registration of five different subjects and five different warpage types

	1	2	3	4	5
1	144.2440 / 40	123.9576 / 46	165.7479 / 45	129.8583 / 60	140.0372 / 62
2	685.4163 / 31	741.1117 / 48	668.2548 / 44	714.8960 / 59	715.1198 / 55
3	595.7664 / 52	626.0091 / 46	611.8217 / 46	583.2656 / 63	597.7054 / 61
4	783.7732 / 34	812.4037 / 53	784.9238 / 39	774.4482 / 69	754.1112 / 60
5	681.7424 / 34	726.7434 / 42	700.0624 / 42	703.2313 / 49	635.3843 / 43

Table 5. Table 5: SYM-REG RMSD and time (s) for registration of five different subjects and five different warpage types

	1	2	3	4	5
1	81.4292 / 51	99.6207 / 50	105.7807 / 29	76.3845 / 31	71.5707 / 31
2	331.0654 / 52	339.2126 / 29	334.6728 / 52	320.4282 / 50	332.5362 / 51
3	420.8638 / 58	436.0106 / 30	429.0045 / 29	482.7290 / 28	424.9221 / 36
4	349.0831 / 50	363.1338 / 51	362.9932 / 31	325.4440 / 28	358.1071 / 33
5	308.6141 / 63	324.2186 / 25	326.5479 / 30	298.3141 / 30	311.8945 / 28

Equations

$$\mathcal{H}(\mathbf{q},\mathbf{p})=\frac{1}{2V}\int \left[{\mathbf{p}}^2+{\left(I_0(\mathbf{x})-I_1(\mathbf{q}))\right)}^2\right]d\mathbf{x}.$$

$$\frac{d\mathbf{q}}{dt}$$

$$\frac{d\mathbf{p}}{dt}$$

$$\frac{d\mathbf{q}}{dt}$$

$$\frac{d\mathbf{p}}{dt}$$

$$\frac{dJ}{dt}$$

$$\underset{0}{\overset{1}{\int }}\left|\frac{d\mathbf{x}}{ds}\right|ds=\underset{0}{\overset{1}{\int }}\sqrt{g_{ij}\frac{d{q}^i}{ds}\frac{d{q}^j}{ds}}ds,$$

$$ϵ<|J(x,t)|<{ϵ}^{-1},$$

$${\mathbf{q}}^{(n)}(\mathbf{x},0)$$

$${\mathbf{p}}^{(n)}(\mathbf{x},0)$$

$$J^{(n)}(\mathbf{x},0)$$

$$J\left(\mathbf{x},t\right)=$$

$$\int [$$

$$\displaystyle Q(\mathbf{x},\mathbf{x}^{\prime})=\begin{cases}1&\mbox{if $\mathbf{x}$ and $\mathbf{x}^{\prime}$ are connected}\\ 0&\mbox{if $\mathbf{x}$ and $\mathbf{x}^{\prime}$ are not connected}\end{cases}$$

$$\lambda \psi (\mathbf{x})=\int Q(\mathbf{x},{\mathbf{x}}^{\prime })\psi ({\mathbf{x}}^{\prime })d{\mathbf{x}}^{\prime },$$

$$\rho ({\mathbf{x}}^{\prime },\mathbf{x})=\frac{Q(\mathbf{x},{\mathbf{x}}^{\prime })\psi ({\mathbf{x}}^{\prime })}{\lambda \psi (\mathbf{x})}$$

$$\mu ({\mathbf{x}}^{\prime })=\int \rho ({\mathbf{x}}^{\prime },\mathbf{x})\mu (\mathbf{x})d\mathbf{x}$$

$${\mathcal{H}}^{nl}(\mathbf{q},\mathbf{p})=$$

$$\frac{d\mathbf{q}}{dt}$$

$$\frac{d\mathbf{p}}{dt}$$

$$\frac{dJ}{dt}$$

$${\mathcal{H}}^{nl}(\mathbf{q},\mathbf{p})=$$

$$\frac{d\mathbf{q}}{dt}$$

$$f_{lmn}^{\{0,1\}}=$$

$$R_{ln}(r)=\frac{1}{\sqrt{{\mathcal{N}}_{ln}}}{j}_l(k_{ln}r),$$

$$I_{\{0,1\}}^{NL}(r,\theta ,\varphi )=\sum _{n=1}^N\sum _{l=0}^L\sum _{m=-l}^l{\mathcal{F}}_{lmn}{f}_{lmn}^{\{0,1\}}{R}_{ln}(r)Y_l^m(\theta ,\varphi ),$$

$$I_{\{0,1\}}^N(r)$$

$$I_{\{0,1\}}^L(\theta ,\varphi )$$

$$s_r$$

$$({\theta }_r,{\varphi }_r)$$

Full text

Symplectomorphic registration with phase space regularization by entropy spectrum pathways.

Vitaly L. Galinsky1,3†, Lawrence R. Frank1,2‡

1 Center for Scientific Computation in Imaging, University of California at San Diego, La Jolla, CA 92093-0854, USA

2 Center for Functional MRI, University of California at San Diego, La Jolla, CA 92093-0677, USA

3 Electrical and Computer Engineering Department, University of California at San Diego, La Jolla, CA 92093-0407, USA

†vit@ucsd.edu

‡lfrank@ucsd.edu

Abstract

The ability to register image data to a common coordinate system is a critical feature of virtually all imaging studies that require multiple subject analysis, combining single subject data from multiple modalities, or both. However, in spite of the abundance of literature on the subject and the existence of several variants of registration algorithms, their practical utility remains problematic, as commonly acknowledged even by developers of these methods because the complexity of the problem has resisted a general, flexible, and robust theoretical and computational framework.

To address this issue, we present a new registration method that is similar in spirit to the current state-of-the-art technique of diffeomorphic mapping, but is more general and flexible. The method utilizes a Hamiltonian formalism and constructs registration as a sequence of symplectomorphic maps in conjunction with a novel phase space regularization based on the powerful entropy spectrum pathways (ESP) framework.

The main advantage of the ESP regularized symplectomorphic approach versus the standard approach of coordinates-only diffeomorphic mapping lies in use of a common metric that remains valid even with image dependent regularization. Moreover, the fusion of the Hamiltonian framework with the ESP theory goes beyond just providing an alternative spatially varying smoothing strategy - it provides an efficient and straightforward way to combine multiple modalities.

The method is demonstrated on the three different magnetic resonance imaging (MRI) modalities routinely used for human neuroimaging applications by mapping between high resolution anatomical (HRA) volumes, medium resolution diffusion weighted MRI (DW-MRI) and HRA volumes, and low resolution functional MRI (fMRI) and HRA volumes. The typical processing time for high quality mapping ranges from less than a minute to several minutes on a modern multi core CPU for typical high resolution anatomical ($\sim 256^{3}$ voxels) MRI volumes.

For validation of the framework we developed a panel of deformations expressed in analytical form that includes deformations based on known physical processes in MRI that reproduces various distortions and artifacts typically present in images collected using these different MRI modalities. Use of this panel allows us to quantify repeatability and reproducibility of our method in comparison to several available alternative approaches. The panel can be used in future studies especially for quantitative clinical validation of different registration approaches.

The registration tool will be available as a part of the QUEST suite from the UCSD Center for Scientific Computation in Imaging (CSCI).

1 Introduction

Modern imaging systems are increasingly capable of acquiring data sensitive to a wide range of physical parameters at multiple resolutions, thus offering greater sensitivity to structural and dynamical information in complex biological systems. However, these technological advancements present the increasingly important theoretical and computational challenge of how to rigorously and efficiently combine, or register, such data in order to be able to accurately detect and quantify subtle and complex system characteristics.

The ability to register image data to a common coordinate system is a critical feature of virtually all imaging studies that require quantitative statistical analysis of group populations, as well as for combining single subject modalities. Consequently, this subject has been the focus of a great deal of research. This has been a focus in computational neuroanatomy which has motivated the developed of diffeomorphic registrations methods [1, 2, 3, 4, 5] for which faster and more efficient algorithms continue to be developed[6, 7, 8, 9, 10], as well as various regularizations [11, 12] and additional enhancements such as local-global mixture, contrast changes, multichannel mapping, etc [13, 14, 15], and the use of probabilistic diffeomorphic registration methods [16, 17]. These registration advancements are important to group analyses and the development of standard atlases [18, 19, 20, 21, 22, 23, 24] which serve a critical role in the standardization of studies. The emergence of diffusion tensor imaging (DTI) methods and their variants for connectivity studies required the extension of diffeomorphic registration methods to accommodate tensor data [25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38]. These methods have had a profound effect on the success of numerous scientific studies on important clinical issues such as Alzheimer’s and traumatic brain injury [39, 40, 41, 42, 43, 44], as well as studies in other organs (cardiac, lungs, etc) [45, 46, 47, 48, 49]. Another important and even more challenging task is a multi-modal registration (i.e. registering T1 and T2 images, or T1 and DTI, etc), as the optimal choice of an appropriate objective function is unknown. Designing and evaluating a universal algorithm that can fit various applications (among subjects, multi-modal within-subject, multi-modal across subjects) is an important problem that needs to be addressed, as existing approaches do not currently posses such universality (see, e.g., [50] for a comprehensive review).

In spite of the abundance of literature and the existence of several variants of diffeomorphic algorithms their practical appeal are still rather limited (possibly due to an interplay of a variety of reasons – speed, accuracy, robustness, complexity, repeatability, etc), as commonly acknowledged even by developers of these registration methods. For example citing the developer of one of the relatively broadly used approaches – Large Deformation Diffeomorphic Metric Mapping [30, 34, 35, 37, 44, 51] – “applications of the LDDMM framework on volumetric 3D medical images still remain limited for practical reasons” [52]. Two large and thorough comparison studies (i.e. [53, 54]) also confirm that although currently available methods are in general able to perform the registration task with varying degrees of success (although some are exceedingly slow and some are not particular accurate), the practical use limitations seem to drive an interest in improvements at least in terms of speed and accuracy.

The recent review paper [55] conducted a retrospective analysis of the past two decades of the field of medical image registration since publication of the original review [56]. It is alarming again that the main conclusion of this twenty years retrospective is that in spite of all the progress in the field of registration “the two major problems mentioned in [56] – validation of registration methods and translation of these to the clinic – are major problems still, which have even been aggravated by the elaboration of registration methods.”

To address these issues we present in this paper a new method that is similar in spirit to diffeomorphic mapping, but is more general and flexible. The transformation is developed within a Hamiltonian formalism [57, 58, 59] in which not just the spatial coordinates are considered, but the entirety of phase space, which is a called a symplectomorphism. This theoretical construct enables a novel flexible, accurate, and robust computational method based on a sequence of energy shell transformations. The incorporation of phase space constraints allows us to use the same simple metric on the space of diffeomorphisms that remains valid even with image dependent regularization, something that is missing in currently available methods.

The generality of the Hamiltonian framework facilitates the inclusion of powerful prior information for spatially varying regularization in phase space using our recently developed method of entropy spectrum pathways (ESP) [61]. This is in contrast with the current state-of-the-art approaches that introduce regularization as a differential form (almost always with constant coefficients) acting on the map itself (see e.g. [62, 63, 64]), that effectively apply regularization as an additional post processing step, thus creating additional problems, especially for validation and comparison between different approaches and even different regularization techniques. Even the existence of several techniques for spatially varying smoothing strategies that have been recently proposed [65, 66] do not remediate this validation issue (and this is in addition to being of rather limited practical utility, possibly adding even more speed–accuracy–complexity issues than providing solutions). Generally speaking, the Hamiltonian approach facilitates validation of different regularizations without destroying or modifying the metric on the space of diffeomorphisms.

The importance and the main advantage of symplectomorphic approach versus coordinates only diffeomorphic mapping can be understood from the fact that use of the same common metric allows quantitative assessment of differences between registrations as well as evaluation of performance for different regularization schemes. An existence of volumetric/surface/line measures allows accurate comparison of features between subvolumes, surface areas or linear curves.

What is even more important is that this fusion of the Hamiltonian framework with the ESP theory goes beyond just providing an alternative spatially varying smoothing strategy. It provides an efficient and straightforward way to combine multiple modalities, for use in tractography, structural and functional connectivity, etc. (although the details of implementation go beyond the subject of this paper and will be reported elsewhere).

Our method also incorporates fast, accurate, and flexible spatial preconditioning using our spherical wave decomposition (SWD) [67]. The SWD approach uses fast FFT–based algorithms to expand images in spherical wave modes and therefore allows to do image resampling, scaling, rotating and filtering with the highest possible order of polynomial accuracy, but at a fraction of a time.

The method is validated on a well characterized numerical phantom and then demonstrated on a set of the “standard” neuro-MRI data acquisitions (HRA, DTI, rsFMRI) routinely collected at our UCSD Center for FMRI (CFMRI). We demonstrate the ability to accurately co-register the data volumes in computational times significantly faster and more accurately than current state-of-the-art methods. The resulting image volumes also demonstrate previously unobserved image contrasts that suggest the ability of our method to uncover more subtle and important structural features in the data.

It is well known that different MRI acquisition schemes and protocols may include a variety of incompatible distortions and artifacts due not only to variations of scanner hardware and pulse sequence designs but also due to intrinsic variations in individual subject morphology, as well as just due to simple motions. Thus validation of a registration method’s ability to disentangle the complex interplay of the acquisition details with the physical effects producing distortions within any particular individuals brain is an exceedingly non-trivial problem. Therefore, in order to facilitate a more quantitative validation of all these different conditions we developed a panel of deformations defined analytically and based on well-known physical effects present in the different MRI modalities. The deformations from the panel can be applied to images of different modalities and acquisition condition and potentially can be appropriate for quick and robust validation in clinical settings as well. This validation approach is somewhat similar to Gaussian deformations used in [68], but our panel includes deformations that can be attributed to a variety of real physical processes present in different acquisition protocols and modalities (i.e. twist, whirl, stretch, etc).

To evaluate the practical aspects of our implementation and to demonstrate the competitiveness of our approach we compared the accuracy and speed of phantom registration with several commonly used registration methods that are often reported as top performers [53] in either speed or accuracy (ANTs Diffeomorphic Demons, ANTs SyN, FSL FNIRT and AFNI 3dQWarp). While a variety of similarity metrics are available, for this paper we used a simple Root-Mean-Square Deviation (RMSD) as a metric to evaluate the accuracy of numerical phantom registration and wall–clock time (that characterizes the human perception of the passage of time from the start to the completion of a task, referred to as time afterwards) as a practical and intuitive measure of the algorithms efficiency.

In summary, this paper utilizes a Hamiltonian formalism to develop a new approach to non-linear flexible image registration. The method builds a diffeomorphic mapping as a sequence of symplectomorphic maps with each map embedded in a separate energy shell. The approach adds a novel phase space regularization based on the powerful entropy spectrum pathways framework. The framework provides a unique opportunity to tailor image details into the regularization scheme by choosing an image derived regularization kernel. A spherical wave decomposition is applied as a powerful preconditioning tool in the position domain to allow accurate and fast interpolation, resampling and estimation of fixed shape rotation and scale. The result is an efficient and versatile method capable of fast and accurate registration of a variety of volumetric images of different modalities and resolutions.

2 Symplectomorphic mapping

We introduce the Hamiltonian function $\mathcal{H}(\mathbf{q},\mathbf{p})$ on a fixed Cartesian grid $\mathbf{x}$ as

[TABLE]

Here $I_{0}$ and $I_{1}$ are two multidimensional images defined on the same fixed Cartesian grid $\mathbf{x}$, $V$ is the measure (volume) of the reference $I_{0}$ image domain ($V\equiv\int d\mathbf{x}$), and $(\mathbf{q}(\mathbf{x},t),\mathbf{p}(\mathbf{x},t))$ is a set of canonical coordinates, that define a time dependent mapping from Cartesian grid $\mathbf{x}$ to a new curvilinear grid $\mathbf{y}\equiv\mathbf{q}(\mathbf{x},t)$, such that initially at $t=0$ the grids are identical, i.e. ($\mathbf{q}(\mathbf{x},0),\mathbf{p}(\mathbf{x},0))\equiv(\mathbf{x},0$).

The Hamiltonian Eqn 1 defines a flow at each location on a fixed grid through a system of Hamilton’s equations

[TABLE]

where $\delta\mathcal{H}/\delta...$ denotes variational (or functional) derivative.

The flow defined by Eqns 2 and 3 is called a Hamiltonian flow and takes place in the space of the coordinates $(\mathbf{q},\mathbf{p})$, which is called phase space. Diffeomorphisms in this phase space are called Hamiltonian diffeomorphisms or symplectomorphisms since a phase space is a symplectic manifold. Thus symplectomorphisms preserve the symplectic structure (including the volume) of phase space. This is a very important feature that will allow the generation of a shell-like sequence of transformations suitable for volumetric measurements and quantifications.

Because the Hamiltonian function Eqn 1 and the reference image $I_{0}$ are defined on a Cartesian grid $\mathbf{x}$ we do not calculate the curvilinear gradient ${\partial I_{1}}/{\partial\mathbf{q}}$ directly. Instead we express $I_{1}(\mathbf{q})$ as a function on a Cartesian grid $I_{1}(\mathbf{q}(\mathbf{x},t))$ and use the chain rule to evaluate the curvilinear gradient through a gradient on Cartesian grid ${\partial I_{1}}/{\partial\mathbf{x}}$ and Jacobian $J\equiv{\partial\mathbf{q}}/{\partial\mathbf{x}}$ as ${\partial I_{1}}/{\partial\mathbf{x}}({\partial\mathbf{q}}/{\partial\mathbf{x}})^{-1}$.

An evolution of the Jacobian with time can be obtained by differentiating the position equation (Eqn 2) on a fixed grid, giving a closed set of equations

[TABLE]

Integrating these equations with initial conditions $\mathbf{q}(\mathbf{x},0)=\mathbf{x}$, $\mathbf{p}(\mathbf{x},0)=0$, and $J(\mathbf{x},0)=\mathds{1}$ generates a symplectomorphic transformation $\mathbf{x}\rightarrow\mathbf{q}(\mathbf{x},t)$. A new metric can be defined for the position part $\mathbf{q}$ of the canonical coordinates by introducing the metric tensor $G\equiv\{g_{ij}\}={(J^{-1})}^{T}J^{-1}$, where indices $i$ and $j$ correspond to derivatives over $q_{i}$ and $q_{j}$ components of the curvilinear coordinates $\mathbf{q}$ such that in Euclidean space $g_{ij}=\delta_{ij}$ where $\delta_{ij}$ is the Kronecker delta. The metric tensor is important for providing accurate measures of line and surface properties using the curvilinear coordinate system $\mathbf{q}$. For example, a length of a curve parameterized by ${\mathbf{x}}(s)$ with a parameter $s$ between zero and one in Cartesian space can be expressed using the metric tensor and curvilinear mapping as

[TABLE]

where repeated indices $i$ and $j$ represent summation.

To ensure that the transformation is symplectomorphic at every location on a fixed grid $\mathbf{x}$ during numerical integration we set a small constant $\epsilon$ and impose a requirement that both the Jacobian and the inverse Jacobian are bounded by this constant, i.e.

[TABLE]

For the majority of the results presented in the paper a value of $\epsilon=0.01$ was used. When the Jacobian becomes sufficiently close to zero the further integration does not make sense as it will not be able to guarantee either the symplectomorphic or diffeomorphic properties of the flow (even numerical stability of the solution can be compromised). Therefore, when the condition of Eqn 8 is violated we stop numerical integration, freeze the flow, and restart the integration (i.e., setting $t=0$) beginning at a new set of phase space coordinate $\{\mathbf{q}^{(n)}(\mathbf{x},0),\mathbf{p}^{(n)}(\mathbf{x},0)\}$ where $n$ is the number of restart times. Since the Hamiltonian is an operator that describes the “energy” of a system, we refer to these $n$ different sets of initial conditions as energy shells. Each restart of the integration therefore represents the initiation of a new energy shell.

The new initial conditions that define the energy shells are related to the stopping point of the coordinates in the previous energy shell by the following conditions:

[TABLE]

Repeating this sequence of initial conditions therefore generates a set of shell-embedded symplectomorphic transformations such that the total transformation is diffeomorphic with the Jacobian defined as a product of $J^{(n)}$

[TABLE]

It is worth noting that this updating equation for the Jacobian effectively results in an updating of the metric tensor $G={(J^{-1})}^{T}J^{-1}$ that characterizes the local geometry and assures volume preservation.

We would like to emphasize that our use of Hamiltonian framework provides a major advantage over conventional approaches in both efficiency and accuracy. For example, similar considerations for limiting the Jacobian were employed in [69] where Euler equations of viscous flow were used to describe the displacement field on a fixed grid. The introduction of fixed Eulerian reference frame resulted in frequent use of costly and inaccurate template regridding procedure that is completely avoided by our formulation.

An important practical implementation issue is that the number of shells $n$ does not have to be introduced in advance and can be determined based on overall convergence (or even devised from running time constraints). In our numerical implementation the shells were terminated as soon as $I_{1}\rightarrow I_{0}$ convergence condition

[TABLE]

was not satisfied.

3 Entropy spectrum pathways as a phase space regularization

The form of Hamiltonian function used in Eqn 1 assumes only local input from difference between $I_{0}$ and $I_{1}$ images to the flow momentum $\mathbf{p}$ at every point on the fixed grid $\mathbf{x}$. A more reasonable assumption would be an inclusion of some information relevant to the structure of $I_{0}$ and $I_{1}$ images. One possible (and by far the most straightforward) way to provide this structure based preconditioning is the entropy spectrum pathways (ESP) approach [61] that takes into account nearest neighbor coupling between adjacent grid locations.

The ESP approach starts with generating the coupling density $Q(\mathbf{x},\mathbf{x}^{\prime})$ which can be as simple and trivial as just the adjacency matrix

[TABLE]

or may in general include a strength of coupling through some kind of coupling potentials that may depend on the grid positions. The ESP approach solves the generalized eigenvalue problem

[TABLE]

finding the largest eigenvalue $\lambda$ and corresponding eigenvector $\psi(\mathbf{x})$ and then constructs the quantity

[TABLE]

calling it the transition probability density for transition between grid locations $\mathbf{x}$ and $\mathbf{x}^{\prime}$. The square of the eigenvector $\psi(\mathbf{x})$ is called the equilibrium probability $\mu(\mathbf{x})$ in the sense that it represents the stationary solution that satisfies the stationary point condition

[TABLE]

Eqn 16 can be included in Eqn 1 to take into account nonlocal effects and provide a way of regularization by defining a non-local Hamiltonian

[TABLE]

here $\delta(\mathbf{x},\mathbf{x}^{\prime})$ is Dirac delta function, $\mathbf{q}\equiv\mathbf{q}(\mathbf{x}^{\prime},t)$ and $\mathbf{p}\equiv\mathbf{p}(\mathbf{x}^{\prime},t)$. This nonlocal expression for the Hamiltonian function produces non-local Hamilton’s equations

[TABLE]

where the momentum equation (Eqn 20) is the non-local version of Eqn 5 that now includes the convolution of a local potential (gradient of squared image difference in our case) with a kernel $\rho(\mathbf{x},\mathbf{x}^{\prime})$ that depends on the coupling between grid locations.

Alternatively the non-local Hamiltonian function can be specified as

[TABLE]

providing alternative non-local form for the coordinate equation (Eqn 4) as well

[TABLE]

Assuming that the coupling density $Q(\mathbf{x},\mathbf{x}^{\prime})$ does not depend on position $\mathbf{x}$ but depends only on a difference between them (i.e. $Q(\mathbf{x},\mathbf{x}^{\prime})\equiv Q(\mathbf{x}-\mathbf{x}^{\prime})$), the ESP scheme can provide a variety of position independent regularization kernels often used as convolution filters in image registration [2]. As a trivial example, an eigenvalue problem (Eqn 15) for position independent Gaussian coupling density $Q(\mathbf{x}-\mathbf{x}^{\prime})=\exp(-(\mathbf{x}-\mathbf{x}^{\prime})^{T}S(\mathbf{x}-\mathbf{x}^{\prime})))$ in infinite $n$-dimensional domain has maximum eigenvalue $\lambda=\sqrt{\pi^{n}/\det{S}}$ and a trivial eigenvector $\psi(\mathbf{x})=\mathrm{const}$, resulting in the commonly used Gaussian regularization kernel. This simple illustration is merely meant to demonstrate that the commonly used Gaussian kernel is naturally derived from our very general procedure. In practice, more complex coupling schemes can provide more informative prior information, resulting in more robust warping schemes.

We would like to emphasize the significant advantages that ESP regularization provides. Its general formulation [61] is probabilistic in nature and provides a framework for the incorporation of available information. In the present context of image registration it naturally provides a mechanism to incorporate information from either or both of the $I_{0}$ and $I_{1}$ images. The position dependent coupling naturally creates image dependent regularization. Moreover, the ESP approach can also include any information that is not present in the images themselves but known a priori and related to images in some quantitative way can be easily included into the coupling scheme with some sort of linear or nonlinear parameterization. We have recently demonstrated this ability to incorporate multiple priors in ESP coupling in the related problem of multi-modal parameter estimation [70], where the symplectomorphic registration method of this paper was used for registration of multiple modalities. Additionally, incorporation of the ESP method into the Hamiltonian formalism provides a simple and efficient way for introduction of different image matching terms by modification of the position–based part of either local or nonlocal Hamiltonian function. This provides great flexibility for tailoring the method to specific applications.

4 Spherical waves decomposition as a position domain preconditioning

The set of Hamilton’s equations (Eqns 19, 20 and 21) used in the previous sections to generate a sequence of energy shell-embedded symplectomorphic transformations (Eqn 12) requires equal dimensionality of images $I_{0}$ and $I_{1}$. However, in many cases the images to be registered are of different spatial resolutions so that some form of interpolation is required. To provide an effective way to do position domain resampling, interpolation, filtering and estimation of best orthogonal transform in a single step we used the spherical waves decomposition (SWD) approach [67].

The SWD approach uses fast algorithms to expand both $I_{0}$ and $I_{1}$ images in spherical wave modes

[TABLE]

where $Y_{l}^{m\star}(\theta,\phi)$ are the spherical harmonics, and $R_{nl}(r)$ can be expressed through the spherical Bessel function

[TABLE]

with an appropriate choice of normalization constants $\mathcal{N}_{ln}$ and the discrete spectrum wave numbers $k_{ln}$ determined by the boundary conditions. The number of modes ($l,m=0\dots L_{max}$ and $n=1\dots N_{max}$) are determined by the highest image resolution. The details of definitions of the spherical harmonics $Y_{l}^{m}(\theta,\phi)$ and spherical Bessel Functions $j_{l}(r)$ can be found in [67]. The interpolation and resampling are then implemented as fast inverse spherical wave transform

[TABLE]

using appropriate grid locations $(r,\theta,\phi)$ and assigning $f_{lmn}$ to zeros for modes with $n>N_{max}$ or $l,m>L_{max}$. A variety of low/band/high pass filters can be used for frequency domain filter $\mathcal{F}$ following the standard image processing techniques.

The scale and the amount of rigid rotation between images can be easily and effectively estimated using the decomposition of the radial and spherical parts using the partial transforms

[TABLE]

and finding the parameters of the similarity transformation (scale $s_{r}$ and rotation angles $\theta_{r}$ and $\phi_{r}$) by solving the two (one and two dimensional) minimization problems

[TABLE]

using small number of modes ($L<L_{max}$ and $N<N_{max}$) for initial coarse search and increasing them to refine the estimate, thus avoiding being trapped in local minimums and at the same time creating computationally efficient approach. Criteria similar to the considerations about optimal order of SWD transform expressed in [67] can be used as a stopping condition for this scale refinement procedure. A similar scale refinement procedure was also used for the calculation of the symplectomorphic mapping.

5 Results and Discussion

All results shown below were coded in standard C/C++ and parallelized using POSIX threads.

5.1 Phantom data registration

To test the approach we first applied it to 3D extension of a classical “toy” example commonly used to show the performance of non-linear registration approaches – the registration of the “C” shape to the ”circle” shape. The original 200x200x200 3D “C” and ”circle” volumes are shown in Figure 1 (a) and (b). Panels (c) and (d) show 3D view and central slices for the forward and inverse maps. Panels (c) and (d) show grid lines for a subset of points selected from three groups of orthogonal planes. The overall performance of our approach seems to be very good, with nearly perfect forward map of the “circle” to the “C” shape, and only a slight signature of the original hole in the inverse map of “C” shape to the “circle”. No SWD preconditioning was used in this example and simple adjacency type matrix was used for phase space coupling.

The second stage of our evaluation procedure was the comparison with some of the existing state-of-the-art non-linear registration methods commonly used. As it is clearly beyond the scope of our paper to conduct a comprehensive evaluation using all possible variations of nonlinear deformation algorithms versus all possible variations of accuracy metrics (such as is done in evaluation studies, like e.g. [53]), we decided to report here just the simplest and the most straightforward type of metrics and have restricted our choice of registration tools to those commonly used at our institution. Following the recommendation from [53] concerning the speed and accuracy, we processed the phantom registration using Diffeomorphic Demons (which is reported as one of the fastest) and symmetric diffeomorphic image registration SyN (reported as one of the most accurate), both from the Advanced Normalization Tools (ANTs) package [71], as well as using FNIRT non-linear registration utility from the FSL [72] and 3dQwarp non-linear warping utility from AFNI[73].

As noted above, our quantitative comparison of accuracy and efficiency are based on two simple and well-understood metrics: The time is used as a practical measure of efficiency and a simple root mean square deviation between the reference and the registered image is used as a measure of accuracy. The evaluation results are summarized in table 1 and demonstrate the enhanced accuracy and efficiency of our symplectomorphic registration approach in comparison with some well established non-linear registration techniques for these two simple metrics and a well-defined standard numerical phantom.

5.2 High resolution anatomical MR data registration

5.2.1 Multi-subject registration

For our first test of the method on actual data, we addressed the most common usage of registration algorithms: to register a set of high resolution anatomical (HRA) images to a common reference image. This is a typical multi subject analysis task appearing in a variety of group studies that involve morphometry based comparison between different subjects or subject groups. We utilized HRA data collected on the 3T GE Discovery MR750 whole body system at the UCSD Functional MRI (CFMRI) using a 32 channel head coil for ten different subjects previously collected in a study to determine the effects of caffeine on the resting state brain activity [74]. However, only high resolution T1 data (all having 290x262x262 voxel resolution) were used for registration test described in the current subsection. Further details are available in [74]. Normalized T1 intensities were used as matching terms $I_{0}$ and $I_{1}$ in symplectomorphic registration.

Figure 2 shows the collage of images related to this registration test. The central plane from anatomical volume used as reference is shown in (a) panel. The same location planes for randomly selected four volumes out of ten subjects are shown in (b). Panel (d) shows the result of SWD preconditioning step equivalent to rigidly fitting each volume to reference with orthogonal transform that includes rotation and uniform scaling for the same four volumes. Panel (c) shows image obtained by averaging of SWD preconditioned volumes for all ten subjects. The next four panels show results of symplectomorphic transforms using 5 (f) and 15 (h) embedded energy shells again with correspondent all subject averages in (e) and (g) respectively. And finally, panel (i) show illustrative plots of curvilinear grids for the same randomly chosen four subjects.

Overall, as would be expected, the symplectomorphic registration shows significant improvement over rigidly fitted volumes, with additional improvement due to increase of a number of energy shells used in registration. In general, there is no obvious relationship between the number of shells and the accuracy, although practically, as the number of shells is determined by selected limits of Jacobian range ($\epsilon$) (and indirectly can be affected by a selected policy of time step adjustments), symplectomorphic registration with increased number of shells may allow to obtain better overall accuracy. The total processing time for all ten subject fitting ranges from 15 to 40 minutes based on the selected quality (this is time measured by running the registration on 12 cores Intel® CoreTM i7-4930K CPU 3.40GHz).

We note, while the quality of the fit is measured here by RMSD, the most advantageous implementation in any particular clinical or research scenario of course depends on on several parameters, including the desired quality of the fit, the type of regularization, the type of coupling in ESP step, aggressiveness of time step updating, etc. These trade-offs necessarily depend on the specific desires of the user in any particular application.

5.2.2 Synthetic deformation maps

An important practical aspect of non-linear registration methods in their application to MRI is their robustness to image distortions produced in the MRI procedure. The dependence of the spatial encoding process in MRI on magnet field gradient linearity results, in practice, in a wide range of complex non-linear distortions due to gradient non-linearities induced by both machine dependent factors (e.g. imperfect gradients) and subject morphological variations (e.g. susceptibility effects). To add to the complexity of this problem, the machine dependent variations can depend not only on the scanner vendor, but on the scanner software revision as well. And the subject morphological variations certainly differ between subjects.

Therefore, in order to provide a quantitative assessment of the symplectomorphic registration using HRA data under more realistic conditions encountered in practical applications, we took several T1 brains acquired on different hardware (Siemens and GE scanners) at different resolutions using different acquisition sequences and subjected them to artificially generated distortions designed to mimic some of the most prominent non-linear distortions common to MRI acquisitions.

We utilized 5 different high resolution datasets: (1) MNI152 T1 2mm with 91x109x91 voxels, (2) T1 MPRAGE 1.2mm with 160x200x200 voxels, (3) T1 MPRAGE 1mm with 212x240x256 voxels, (4) T1 1mm with 290x262x262 voxels, and (5) T1 1mm with 256x176x176 voxels. All subjects were resampled into MNI152 space (to provide the same accuracy and required workload for all subjects) and then were distorted several times using five different types of warpage (see Figure 3). The different types of warpage include nonlinear or differential rotation, nonlinear stretching, nonlinear compression, etc, and Figure 3 shows one of the subjects with all five warpage types applied to this subject.

The idea of generation and use of synthetic deformation maps for validation is not new. A somewhat similar procedure of synthetic Gaussian 2D deformations was used recently for registration validation in [68]. The important difference of our deformation panel is that it includes a variety of more complex deformation modes observed and frequently present in various MRI acquisition protocols and modalities.

The 25 warped volumes were then processed using the two non-linear registration methods most commonly used at our institution – ANTs SyN and AFNI 3dQwarp, in addition to SYM-REG, and compared the restored volumes with the original unwarped datasets. The results are summarized in table 2 showing mean values for the RMSD as well as mean execution times for all subjects (top) and for all warpage types (bottom). The complete processing took 976, 3414 and 1223 seconds for SYM-REG, AFNI and ANTs respectively giving 304.6, 358.5 and 584.0 for the RMSD values with SYM-REG providing both the lowest RMSD error and the lowest execution time (the default processing options were used for all packages).

5.3 Diffusion weighed MR image registration

The second important application we investigated was the registration of diffusion-weighted images (DWI) to the HRA image of the same subject. This is a critical step required for collocation of diffusion tractography based quantities with high resolution anatomical morphometry. Registration of DWI to HRA is generally a more complicated problem than registering HRA to HRA because DWI images are typically acquired with echo-planar imaging (EPI) acquisitions that are not only prone to more severe non-linear susceptibility distortions than HRA images but are invariably acquired at lower spatial resolution than the HRA for the same subject.

The data used for this example were again collected at the UCSD Center for Functional MRI (CFMRI) using 3T GE Discovery MR750 whole body system to study the effects of traumatic brain injuries (TBI). The HRA T1 volume has 168x256x256 voxel size with 1.2x0.9375x0.9375mm resolution.

A multiband DTI EPI acquisition [75] using three simultaneous slice excitations was used to acquire data with three diffusion sensitizations (at b-values b=1000/2000/3000 s/mm2) for 30, 45 and 65 different diffusion gradients (respectively) uniformly distributed over a unit sphere. Several baseline (b=0) images were also recorded. The data were reconstructed offline using the CFMRI’s multiband reconstruction routines. The DWI data has 100x100x72 voxel size with 2mm3 resolution. The normalized HRA T1 intensity was used for the reference image matching term $I_{0}$ and the equilibrium probability map (see [76]) was used as a moving image matching term $I_{1}$.

Figure 4 shows a central slice and a 3D view of the reference volume (a), DWI b=0 volume (b), DWI equilibrium probability volume that has the same resolution as b=0 volume (c), and the final symplectomorphic registration of the DWI equilibrium probability volume (d). The details about the equilibrium probability and how it is obtained can be found in [76]. The last two panels (e-f) show enlarged side by side comparison of the HRA reference and transformed DWI with the same resolution.

Two tractography examples – one that applies symplectomorphic registration to a set of tracts as a postprocessing step, and the other that includes high resolution symplectomorphic grid together with the HRA data during the tracking stage – are shown in Figure 5. The figure shows that the symplectomorphic registration method allows very accurate localization of diffusion derived tracts with the high resolution anatomical features. More details that include a SYM-REG based tractography implementation as well as multi-modal estimation in general are not relevant to this paper and are reported in [70].

We would like to mention one important consideration here. A considerable amount of work has been spent recently not just on spatial registration of diffusion imaged volumes, but also on devising techniques for local reorientation of diffusion tensors that would be consistent with the new deformed spatial grid [see e.g [29, 25, 26, 27, 30]]. These methods are both time consuming and an unnecessary intermediate step from our viewpoint. An important feature of our method is that we can directly import the diffeomorphic maps together with the high resolution data into our diffusion estimation and tractography technique GO-ESP [76], so that both the estimation of local diffusion properties and the generation of tracts are performed in the locally warped space characterized by the spatially dependent metric tensor, thereby obviating the need to use any ad hoc proceed to impose geometric consistency, which in our method is guaranteed by the symplectomorphic nature of transformation. The result is a method that provides a fast and effective way of adding a new level of details to relatively low resolution output available from diffusion weighted tractography.

5.4 Registration of functional MR images

The third and final important application we investigated was the registration of low spatial resolution functional resting state FMRI (rs-FMRI) data to HRA data from a single subject. As in DWI, FMRI data is typically acquired using EPI acquisitions that, as previously mentioned, are more prone to non-linear geometric distortions than HRA acquisitions, and are of lower spatial resolution.

The data used for this test were from the same caffeine study dataset [74] used in section 5.2. Registration of functional rs-FMRI data to anatomical images is required to establish an accurate localization of activation regions in the high resolution maps of gray matter. Only the data collected prior to caffeine administration were used. Whole brain BOLD resting-state data were acquired over thirty axial slices using an echo planar imaging (EPI) sequence (flip angle = $70\degree$, slice thickness=4mm, slice gap=1mm, FOV=24cm, TE= 30 ms,TR = 1.8 s, matrix size = $64\times 64\times 30$). Further details are available in [74]. All data were pre-processed using the standard pre-processing analysis pathway at the CFMRI (as described in [74]). Nuisance terms were removed from the resting-state BOLD time series through multiple linear regression. These nuisance regressors included: i) linear and quadratic trends, ii) six motion parameters estimated during image co-registration and their first derivatives, iii) RETROICOR (2nd order Fourier series) [77] and RVHRCOR [78] physiological noise terms calculated from the cardiac and respiratory signals, and iv) the mean BOLD signals calculated from WM and CSF regions and their first respective derivatives, where these regions were defined using partial volume thresholds of 0.99 for each tissue type and morphological erosion of two voxels in each direction to minimize partial voluming with gray matter. The normalized HRA T1 intensity was used for the reference image matching term $I_{0}$ and the temporal mean of functional activation was used as a moving image matching term $I_{1}$.

Figure 6 shows side by side comparison for 3D views of rs-FMRI (a), T1 (b) and rs-FMRI mapped to T1 (c) volumes. The processing was carried through 30 energy embedded shells and required about 5 minutes of waiting time from the start to the finish, with a subset of the final grid shown in (d). The work is currently underway to include flexible mapping grids directly to our rs-FMRI mode detection approach [79, 80].

Figure 7 shows comparison of original and registered images for several of the resting state functional modes obtained using our entropy field decomposition (EFD) technique [79, 80]. The default mode (a) and (d), the visual lateral (b) and (e), and the visual occipital (c) and (f) modes are shown for some of the subjects from Figure 2. The symplectomorphically mapped images overlayed over correspondent HRA slices (d, e, and f) show very accurate localization of functional modes in the appropriate regions of HRA volumes.

6 Conclusions

In this paper we have presented a new flexible multidimensional image registration approach that is based on the Hamiltonian formalism. The method generates a set of Hamilton’s equations capable of producing a symplectomorphic transformation for mapping between Cartesian and curvilinear grids that minimizes some predefined image difference metric. The final diffeomorphic mapping is constructed as a multiplicative sequence of symplectomorphic transforms with gradually diminishing levels of total energy, thus providing a sequence of energy embedded symplectomorphic shells. For demonstration purposes, we used both a simple local squared difference, as well as a more complicated non-local image squared difference, as a Hamilton function.

An application of the powerful and versatile ESP approach [61] to the phase space domain resulted in a non-local form of Hamilton’s equations. The non-local form represents an efficient and relatively straight forward way to introduce regularization that is capable of taking into account some image specific details or even additional knowledge based parameterizations. More generally, the Hamiltonian formalism allows easy adaptation of custom and possibly more complex forms of image difference metrics, and at the same time allows the metric in the diffeomorphism space to be kept the same, which facilitates the comparison and validation between different techniques and regularization approaches.

The resolution differences between images as well as rigid shape alignment is addressed using the preconditioning step based on the SWD technique [67]. This efficient volumetric decomposition computes a set of fast spherical harmonics and spherical Bessel transforms and is able to produce accurate interpolation, filtering and fitting of rigid shapes.

One of the major hurdles in the translation of registration methods to clinical practice [55] is the difficulty in validation of these methods in the face of the exceedingly complex interplay of the many pulse sequence details, the scanner hardware, and the anatomical variations that interact with the magnetic fields. To this end, we have developed a panel of analytically defined deformations based on known physical effects present in MR acquisitions. While these certainly do not constitute an exhaustive list of possible effects, they are a fair representation of the major distortions in the major MR image acquisition schemes. This deformation panel was used to validate our method and quantitatively demonstrate its robustness to a wide range of different types of distortions that regularly plaque MR studies. The deformations from the panel can be applied to images of different modalities and acquisition conditions and potentially can be appropriate for quick and robust validation in clinical settings. This validation approach is somewhat similar to 2D Gaussian deformations used in [68], but our panel includes deformations that can be attributed to a variety of real physical processes present in different acquisition protocols and modalities (i.e. twist, whirl, stretch, etc).

Overall the symplectomorphic registration approach is both accurate and fast and is capable of processing of a variety of volumetric images of different modalities and resolutions. In the tests reported in this paper we were able to handle all three of the major neuro-MRI modalities routinely used for human neuroimaging applications, including mapping between high-resolution anatomical volumes, medium resolution diffusion weighted volumes and high-resolution anatomicals, and low resolution functional MRI images and high-resolution anatomicals. The typical processing time for high quality mapping ranges from less than a minute to several minutes on a modern multi-core CPU for a typical high resolution anatomical MRI volumes. The speed, accuracy, and flexibility of this new method has the potential to play an important role in the quantitative assessment of neuroimaging data in a wide range of both basic research and clinical applications.

Acknowledgments

The authors thank Dr Alec Wong and Dr Tom Liu at the UCSD CFMRI for providing the resting state data and Dr Scott Sorg at the VA San Diego Health Care System for providing the diffusion weighted imaging data. LRF and VLG were supported by NSF grants DBI-1143389, DBI-1147260, EF-0850369, PHY-1201238, ACI-1440412, ACI-1550405 and NIH grant R01 MH096100.

Appendix A Individual subject and warpage accuracy and timings

Tables 3, 4 and 5 provide values of RMSD and wall execution time for all subjects and all warpage types after processing by AFNI 3dQwarp, ANTs SyN, and SYM-REG respectively.