Barnes-Hut Approximation for Point SetGeodesic Shooting

TL;DR

This paper introduces a Barnes-Hut based approximation method for point set geodesic shooting that significantly speeds up computations while maintaining accuracy, enabling practical applications in clinical research.

Contribution

It presents a novel Barnes-Hut approximation for point set geodesic shooting, reducing computational complexity from quadratic to near-linear, with comparable accuracy.

Findings

Achieves up to 3-fold speedup in computations.

Maintains accuracy comparable to exact methods.

Enables new clinical and practical applications.

Abstract

Geodesic shooting has been successfully applied to diffeo-morphic registration of point sets. Exact computation of the geodesicshooting between point sets, however, requiresO(N2) calculations each time step on the number of points in the point set. We proposean approximation approach based on the Barnes-Hut algorithm to speedup point set geodesic shooting. This approximation can reduce the al-gorithm complexity toO(N b+N logN). The evaluation of the proposedmethod in both simulated images and the medial temporal lobe thick-ness analysis demonstrates a comparable accuracy to the exact point set geodesic shooting while offering up to 3-fold speed up. This improvementopens up a range of clinical research studies and practical problems towhich the method can be effectively applied.

Tables

Table 1. Table 1: Process time (min) and residual error (mm) of cases (b) and (c)

	Number of points	BH-approx		Non-approx
		Time	Residual error	Time	Residual error
Circles (Fig. 2b)	1200	124.6	3.4002e-5	40.4	1.5964e-5
Flat shape (Fig. 2c)	1202	23.0	0.1558	49.5	0.1507

Table 2. Table 2: Matching accuracy and processing time.

	Time(min)	ERC DSC	BA35 DSC	BA36 DSC	PHC DSC
Non-approx (1275 pts)	$23.2\pm 4.7$	$0.92\pm 0.02$	$0.92\pm 0.02$	$0.92\pm 0.02$	$0.93\pm 0.01$
BH-approx (1275 pts)	$7.7\pm 0.7$	$0.92\pm 0.02$	$0.92\pm 0.02$	$0.92\pm 0.02$	$0.93\pm 0.01$
BH-approx (2540 pts)	$24.5\pm 2.2$	$0.93\pm 0.02$	$0.93\pm 0.02$	$0.93\pm 0.02$	$0.94\pm 0.01$

Equations

$$\{\begin{array}{r}\frac{dq}{dt}=\frac{\partial H}{\partial p}(q,p),\\ \frac{dp}{dt}=-\frac{\partial H}{\partial q}(q,p)\end{array}$$

$$p_0^{\ast }=p_0\in {\mathbb{R}}^{3L}H(q_0,p_0)+\lambda ||q(1)-X_m|{|}_2^2$$

$$v(x,t)=\sum _{i=1}^N{G}_{{\sigma }_{gs}}(||x-q_i(t)|{|}_2)p_i^{\ast }(t),x\in {\mathbb{R}}^3$$

$$v_{jc}=\sum _{k=1}^{K_j}{G}_{{\sigma }_{gs}}(||q_j^{\ast }-q_{kc}|{|}_2)p_{kc},$$

Full text

1 Introduction

In medical imaging, diffeomorphic registration is widely used to find a mapping from a biological structure to another since it ensures smoothness of the transformation and prevents folding of the manifold. In practice, the image-based computation of a diffeomorphic registration can be reduced to point matching of surfaces and curves [10] to minimize computational cost. Besides, the high-dimensional, dense time-varying velocity field in the registration can be represented as an initial momentum in a low-dimensional linear space, namely geodesic shooting [1]. This representation of the deformation field allows a convenient characterization of subject variability by simple linear statistical methods, such as principal component analysis and principal geodesic analysis [6].

Despite its favorable theoretical properties and relevance to statistical shape analysis, point set geodesic shooting has not been widely used in population studies or clinical applications because of the methods high computational complexity. Previous studies on surface matching under the diffeomorphic framework have produced valuable theoretical results and actual algorithms in applications. [15] builds a norm on the Hilbert space to represent surfaces as currents and derives a surface matching algorithm that preserves diffeomorphism. Extending this analysis under the Hilbert space, [7] proposes a matching criterion to directly compare two curves with not-necessarily matched points. Previous works have attempted on minimizing computational cost. [14] performs a Wendland kernel bundle stationary velocity field (wKB-SVF), and [12] exploits an association graph to directly compare meshes.

In this work, we present an efficient algorithm on diffeomorphic point registration via geodesic shooting suitable for practical applications. Inspired by [8] and [3], our work reduces an $O(N^{2})$ pairwise velocity computation step in geodesic shooting to two computationally light steps. The first step is an oct-tree construction step of all points by Barnes-Hut algorithm [8], with time complexity $O(NlogN)$. The second step is a discontinuous conical Gaussian kernel approximation step, which requires a $O(Nb)$ time complexity, where $b$ is the average number of points within a 3$\sigma$ neighborhood. The proposed method is evaluated with both synthetic images and T1-MRI scans in an Alzheimer’s disease(AD) medial temporal lobe (MTL) thickness study.

2 Methods

2.1 Hamiltonian Formulation of Point Set Geodesic Shooting

Let $X_{p}$ be the set of template points, $X_{m}$ the corresponding set of target points, and both sets have $N$ points. The algorithm determines a spatial transformation $\phi$ that optimally matches the points $X_{p}$ to $X_{m}$ under the diffeomorphic constraint. The objective function is a minimization of a deformation constraint term regularizing the transformation $||v||$ from time 0 to time 1 and a data attachment term minimizing the distance of the transformed points and the target. As proposed in [1], this function can be formulated as a minimization of a sum of kinetic energy and potential energy, and be solved under the Hamiltonian framework. Our notation is consistent with [1] for clarity.

Let the positions $q$ of the points be a function of time $t\in[0,1],q(t)={q_{1}(t),...,q_{N}(t)}$ and the momentum of the points be $p(t)={p_{1}(t),...,p_{N}(t)}$, $(p_{i}(t)\in\mathbb{R}^{3})$. The initial positions of the points are given by $q(0)=X_{j}$, and the initial momenta of the points $p(0)$ are the unknowns that the geodesic shooting algorithm will determine. The evolution of the system is formulated in terms of the Hamiltonian $H(p,q)=<p,K(q)p>$, where $H(p,q)$ is the kinetic energy of the system and is constant over time and $K(q)$ is a $3N\times 3N$ with $N\times N$ diagonal blocks, with the $(i,n)$-th block equal to $G_{\sigma_{gs}}(||q_{i}-q_{m}||)\textperiodcentered I_{3}$, $G_{\sigma_{gs}}$ being a Gaussian kernel and $I_{3}$ being an identity matrix. the evolution of the system is formulated in terms of the Hamiltonian:

[TABLE]

The point matching problem is formulated as an optimization of

[TABLE]

This minimization problem is discretized in time and solved using gradient-based optimization. In our implementation, we use L-BGFS [11] from VNL numeric library [9] to perform quasi-second order update.

Given an initial momenta $p^{*}_{0}$ and the corresponding point trajectories $q^{∗}(t)$ , the point transformation can be interpolated over the entire spatial domain $x\in\Omega$ to yield a smooth velocity field

[TABLE]

The velocity of each point is dependent on all other points with a Gaussian decay by equation 3: $v_{j}=\sum_{i=1,i\neq j}^{N}G_{\sigma_{gs}}(||q_{j}^{*}-q_{i}||_{2})p_{i},\forall j=1...N$. This equation is used to update each point $q_{j}$, and to calculate its exact velocity. The total computation complexity of this procedure is $O(N^{2})$ since it involves calculating the momentum of all other points weighted by the Gaussian distance kernel.

2.2 Barnes-Hut Approximation of Gaussian Kernel

Based on the observation that the Gaussian of two distant points is close to 0, a grouping algorithm is proposed to assemble a cloud of points far away from a target point as a unit and use the average information of the group to avoid redundant calculation. This approximation algorithm is composed of two steps: 1) 3D spaces are recursively sub-divided into oct-trees and 2) each point is traversed down the trees and interaction with other tree nodes is calculated. Since the same procedure applies to each time step, the time variable $t$ is dropped in the following equations for simplicity.

2.2.1 Oct-tree Representation of Space

An oct-tree of all points is constructed so that a node in a tree consisting of multiple points can be viewed as a whole to reduce calculation of all points. The oct-tree is constructed by subdividing the space recursively until each leaf node only contain a single point. First the global 3D bounding box of all points is found by $q_{\min/\max}={\min/\max}{q\in q_{i},i=1...N}q_{i}$. The root node is initialized to be $q_{min}$, $q_{max}$. Then points are added recursively to the tree and the tree structure gets updated. Suppose the $k-$th tree node $\mathcal{N}_{k}$ representing a cubic space is delimited by $q_{kmin}$, $q_{kmax}$, then each non-leave node contains eight child nodes, corresponding to the eight quadrant defined by bi-partition of the space between $q_{kmin}$ and $q_{kmax}$.

Given a tree node $\mathcal{N}_{k}$ and a point $q_{j}$, the recursion policy is devided into three scenarios: (1) $\mathcal{N}_{k}$ is a non-leaf node, then $q_{j}$ is passed to the corresponding child node and recursion preceeds (2) $\mathcal{N}_{k}$ is an empty leaf node, then $q_{j}$ is added to $\mathcal{N}_{k}$ and the recursion stops, and (3) $\mathcal{N}_{k}$ is an occupied leaf node, then $\mathcal{N}_{k}$ is subdivided to eight quadrants and $q_{j}$, $q_{i}$ are passed to the corresponding child nodes. When $q_{j}$ is passed to node $\mathcal{N}_{k}$, some node statistics are updated accordingly: (1) the node’s total momentum $p_{kc}\longleftarrow p_{kc}+p_{j}$, (2) number of points $n_{k}\longleftarrow n_{k}+1$, (3) center of point position $q_{kc}\longleftarrow(q_{kc}\cdot n_{k}+q_{j})/(q_{j}+1)$, (4) actual range of the node $q_{kmax}\longleftarrow\max(q_{kmax},q_{j})$ and $q_{kmin}\longleftarrow\min(q_{kmin},q_{j})$. This information is a compact representation of the point subset in the node.

In each step of the backward stage (gradient descent stage), the update of the point position derivative is inserted into the tree node $\mathcal{N}_{k}$ by $\alpha_{kc}\longleftarrow\alpha_{kc}+\alpha_{i}$ and momentum derivative $\beta_{kc}\longleftarrow\beta_{kc}+\beta_{i}$. The node statistics is modified in a similar way, including total point position derivative, total momentum derivative, number of points, center of point position, and range of node. The complexity of the tree building is of $O(NlogN)$ in the best balanced oct-tree case.

2.2.2 Tree-based Forward/Backward Gaussian Kernel Calculation

The forward/backward pairwise Gaussian interaction is calculated based on the spatial oct-tree. For a point $q_{j}$ and a tree node $\mathcal{N}_{k}$, the recursion policy is again split into three scenarios in both forward and backward stages. (1) If $\mathcal{N}_{k}$ is a leaf node with only one point $q_{i}$, then the exact Gaussian kernel between $q_{j}$ and $q_{i}$ is calculated. (2) If $\mathcal{N}_{k}$ contains multiple points that are all farther than a distance threshold (3$\sigma$), then a single Gaussian interaction is calculated using the summarized total momentum $p_{kc}$ and center position $q_{kc}$. (3) If $\mathcal{N}_{k}$ contains multiple points but within the distance threshold, then the node gets subdivided. The tree traversal ideally has a complexity $O(Nb+NlogN)$, where $b$ is the number of neighboring points within the distance threshold. In the experiment, we have $\sigma=2$, $k\approx 200$, $N\approx 2000$ and $b\approx 200$. Therefore, our scheme provides great speed up in theory.

All accurate calculations can now be replaced with tree node approximations. In the forward stage, the actual Gaussian kernel weighted momentum is replaced with the approximation:

[TABLE]

where $K_{j}$ is the number of tree nodes traversed from $q_{j}$. The dot product between the exact momentum is also replaced with the average momentum, i.e., $p_{ji}=p_{j}\cdot p_{i}\longrightarrow p_{jk}=p_{j}\cdot p_{kc}/n_{k}$.

In the backward stage, besides the Gaussian kernel approximation and the momentum approximation which are the same as in the forward stage, the point position derivative and momentum derivative is replaced with their tree averaged counterparts, i.e., $\alpha_{i}\longrightarrow\alpha_{kc}/n_{k}$ and $\beta_{i}\longrightarrow\beta_{kc}/n_{k}$.

2.2.3 Implementation details

The algorithm is implemented in C++ with ITK 4.11 /VTK 7.1 libraries. Our code is compiled with gcc 4.4.7 and deployed on a cluster setting with E5-2643 v3 CPU, 16GB RAM, Centos 6 without multi-threading.

3 Evaluation

3.1 Evaluation on Synthetic Data

The effectiveness of the proposed approximation is tested on (1) preservation of the diffeomorphic property (Fig. 2a, b) and (2) speed gain (Fig. 2c). In Fig. 2a, two circles (blue and green) are simultaneously registered to their counterparts (red and purple), while the two deformation fields do not collide or collapse. Synthetic case (b) registers the blue circle to the red. The shrinking-expanding artifact in the trace (green) is caused by the diffeomorphic constraint, the same as in [2]. Fig. 2c shows a registration of a long rectangular mesh (blue) to its bent counterpart (red), demonstrating a hypothetical use case where the Barnes-Hut approximation will speed up the registration significantly. Most points in this shape are far away from each other by 3$\sigma$, where Barnes-Hut approximation is useful. With a similar matching accuracy (evaluated as the sum of the square distance), our approximation achieves 2.15 folds speedup. (Table 1). On the other hand, when points are relatively concentrated, such as in case (b), the Barnes-Hut approximation is much slower than the accurate version.

3.2 Evaluation on Cortical Thickness Analysis of the Medial Temporal Lobe in Alzheimer’s Disease

To demonstrate the utility of the proposed technique in the real application, we apply the pipeline to perform thickness analysis for the medial temporal lobe (MTL) cortex in the context of early detection of AD. Since MTL is one of the earliest regions affected by neurofibilary tangle pathology (NFT), a biomarker directly linked to neuronal damage, the thickness measurement of the MTL subregions is a promising biomarker of AD.

3.2.1 Template Matching

Baseline T1-weighted MRI scans of 665 subjects from Alzheimer’s disease neuroimaging initiative (ADNI) were chosen in this study. An open-source multi-atlas segmentation pipeline proposed by Xie et al. [17] was applied to segment the MTL cortex subregions including the entorhinal cortex (ERC), Brodmann areas 35 and 36 (BA35/36) and parahippocampal cortex (PHC). The greedy registration package 111sites.google.com/view/greedyreg was applied to register the automatic segmentation to the MTL cortex template 222Publicly available at github.com/LongXie/Multitemplate-for-MTLCortex-T1 as described in [16]. The template space consists of the MTL cortex labels and a dense surface mesh of the union of ERC, BA35, BA36, and PHC labels (Fig. 3a), from which a set of $N$ points were uniformly sampled using Poisson Disk Sampling [4]. The points were then warped to the space of each subject. The alignment of the points in the subject space to the template space consists of an initial rigid registration step using the Procrustes algorithm [5], and a geodesic shooting step. The dense template mesh and the template labels are then warped to the space of the subject segmentation (Fig. 3c).

3.2.2 Cortical Thickness Measures and Statistical Analysis

To extract regional thickness, the pruned Voronoi skeleton [13] is first extracted from the target mesh and the distance between each vertex and the closest point on the skeleton is computed (Fig. 3e). The median of the thickness measures of each label is then extracted as the summary thickness measure of each label. Bilateral thickness measure of the same label is averaged. The Dice similarity coefficient (DSC) of each MTL cortex label between the warped template and the automatic segmentation is computed, as a measure of the quality of fit.

3.2.3 Discussion

The importance of faster implementation of point set geodesic shooting was motivated by the need to register data sets that are either very densely sampled or large in scale – or in many circumstances both, for example, the ADNI study of MTL. The performance of the Barnes-Hut approximation and the original geodesic point shooting are compared on the same point data set to evaluate the speed gain and point matching accuracy. The number of points is set to $N=1275$ so that the procession time of the entire dataset ($665$ subjects $\times 2$ hemispheres) is reasonable. As can be seen in Table 2, the Barnes-Hut approximation achieves approximately three-fold speedup compared to the original implementation with almost the same DSC for all the labels.

We hypothesize that both the matching accuracy and discriminative power can be improved when processing with points sampled at a higher resolution. The repeated experiment with twice the sample density, i.e., 2540 points, shows a similar processing time (24.5 min) to the 1275 points experiment (23.2 min) but a higher DSC (Table 2). This finding can potentially be attributed to a more accurate thickness measurement with a higher resolution point set sampling of the imaged anatomy.

4 Conclusion

In this work, we present an approximation scheme based on the Barnes-Hut algorithm to perform geodesic point shooting, which can have a three-fold speed up without compromising matching accuracy. The proposed method allows calculation based on more densely sampled data sets, which in turn can translate to more accurate registrations and in turn more sensitive imaging-based markers.