# Geometry-Aware Multiscale Image Registration Via OBBTree-Based Polyaffine Log-Demons

Christof Seiler, Xavier Pennec, Mauricio Reyes Aguirre

► To cite this version:

Christof Seiler, Xavier Pennec, Mauricio Reyes Aguirre. Geometry-Aware Multiscale Image Registration Via OBBTree-Based Polyaffine Log-Demons. MICCAI 2011 - Medical Image Computing and Computer Assisted Intervention, 2011, Toronto, Canada. pp.631-638, 10.1007/978-3-642-23629-7_77. inria-00616215

HAL Id: inria-00616215

https://inria.hal.science/inria-00616215

Submitted on 28 Jun 2013

**HAL** is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire **HAL**, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Geometry-Aware Multiscale Image Registration
Via OBBTree-Based Polyaffine Log-Demons

Christof Seiler¹,², Xavier Pennec², and Mauricio Reyes¹

¹ Institute for Surgical Technology and Biomechanics, University of Bern,
Switzerland

² Asclepios Research Group, INRIA Sophia Antipolis, France

**Abstract.** Non-linear image registration is an important tool in many areas of image analysis. For instance, in morphometric studies of a population of brains, free-form deformations between images are analyzed to describe the structural anatomical variability. Such a simple deformation model is justified by the absence of an easy expressible prior about the shape changes. Applying the same algorithms used in brain imaging to orthopedic images might not be optimal due to the difference in the underlying prior on the inter-subject deformations. In particular, using an un-informed deformation prior often leads to local minima far from the expected solution. To improve robustness and promote anatomically meaningful deformations, we propose a locally affine and geometry-aware registration algorithm that automatically adapts to the data. We build upon the log-domain demons algorithm and introduce a new type of OBBTree-based regularization in the registration with a natural multiscale structure. The regularization model is composed of a hierarchy of locally affine transformations via their logarithms. Experiments on mandibles show improved accuracy and robustness when used to initialize the demons, and even similar performance by direct comparison to the demons, with a significantly lower degree of freedom. This closes the gap between polyaffine and non-rigid registration and opens new ways to statistically analyze the registration results.

# 1 Introduction

In orthopaedic research, and particularly in reconstructive trauma of mandibles, there are a number of surgical interventions such as tumor resection, fracture reconstructions, osteomyelitis or other bone loss repair, which requires implantation of a reconstructive plate. Conventional plate designs, based on shape analysis of a typically small cohort of cadaver specimens or computerized tomography (CT) derived data, fail at capturing the anatomical complexity of the mandible and its shape variation as encountered in a population. In many cases, the results are unsatisfactory and complications such as plate exposure, plate fracture, and screw failure may occur. For instance, it has been reported that suboptimal plate design leads to plate fracturing in 2.9 to 11% of implantations [4]. Additionally intra-operative bending and re-bending of a plate leads to residual stresses,
which affect the mean stress in fatigue loading. The solution is to manufacture
the plate as close to the human anatomy as possible. For this the morphometry
of the mandible should be analyzed in 3D. In the last years, methodologies based
on computational anatomy techniques have been explored for shape analysis of
the mandible [8], and for population-based orthopaedic femural plate design [7].
A common key step in these works is the need to capture the shape variability
as encountered in a population, which is typically performed through non-rigid
image registration.

Many existing non-linear registration algorithms for medical images require
target deformations to be smooth and invertible. In absence of more precise
priors, this could be a reasonable assumption, but it often leads to many local
minima far from the optimal deformation. The complexity of the shape (e.g.
condyle and coronoid process in the mandible) and the encapsulated popula-
tion shape variability makes the registration a challenging task, which calls for
improvements in non-rigid registration for orthopaedic research. Whenever ad-
ditional knowledge about the anatomy is known, we should exploit it to guide
the search of the deformation towards anatomically more meaningful ones.

We hypothesized that a hierarchical anatomically-based characterization of
the structure under study can improve the quality of non-rigid registration re-
sults. In this paper we present a novel method, which incorporates such a de-
scription of the geometry of the anatomy into the registration process. The
method builds on the work of [9], where the authors presented a polyaffine-
regularized demons algorithm, however we believe that defining the regions for
the locally affine deformations should not be left to the user. This becomes even
more evident in the case of a multiscale representation of the geometry, where
the definition and division process is not clear, not to mention the time aspect.
Thus we propose an automatic and anatomically motivated hierarchical scheme.

Our contributions are three-fold: First, we demonstrate that registration ac-
curacy and robustness can be improved by incorporating a locally affine and
anatomy-aware hierarchical scheme. Second, we demonstrate that the complex-
ity of the target deformation can be reduced, which is of great interest in statis-
tics for computational anatomy. The complexity reduction stems from the fact
that in contrast to standard multilevel schemes we introduce a higher degree of
freedom per region by allowing the region to undergo an affine transformation
as opposed to a mere translation. Furthermore, the affine transformation pro-
vides important information about the anatomy since it is implicitly defined by
the geometry of the anatomical structure. Third, we visualize that our approach
opens new ways for multiscale statistics in medical imaging.

## 2 An OBBTree-Based Polyaffine Log-Demons

Our main contribution in this paper is to show that the demons algorithm can
be enhanced in orthopedics applications with an appropriate model. The main
idea of the proposed algorithm is to introduce a multiscale regularization im-
plicitly defined via the geometry of the anatomical structure under study. In the
following sections we review three methods that will be important to formulate our registration: Log-demons, polyaffine registration and polyaffine log-demons. Then we present the new multiscale scheme and the integration into the log-demons.

**Log-Demons Registration** To setup correspondences between anatomical images, a set of images are registered to a reference. We use the recent diffeomorphic log-domain demons registration approach described in [10]. One interesting point of this registration framework is the efficient optimization in the domain of stationary velocity fields. These velocity fields can be looked at as generators for diffeomorphic deformations through the group exponential map that can be very efficiently computed using the scaling and squaring method [2]. This property explains the denomination of log-domain (or simply log-demons) registration.

**Polyaffine Registration** Polyaffine transformations were introduced in [3] to fuse locally rigid and affine transformations into a diffeomorphism. The basic idea is to consider each local affine transformation as the flow of a speed vector field obtained through the log of the affine transformation. Then, instead of averaging the affine matrices, one averages these vector fields with spatial weights describing the influence of each region. The flow of the resulting vector field automatically gives a diffeomorphic transformation. In [5] the authors introduced an efficient approximation of that framework. However, the regions are manually defined on the reference prior to the registration procedure.

**Polyaffine Log-Demons Registration** The polyaffine framework was later enhanced in [1] to work with stationary vector fields. This new formulation suggested that polyaffine transformations could be compatible in some sense with the log-demons. In [9] the authors present a marriage of the two approaches and showed a specific application to femur bones with three fixed regions, head, shaft and condyles. However manually choosing the number of regions is a difficult process which is not easy for new applications. In case of hierarchical schemes, this is an intractable task, due to the need to find regions for each level, and the fact that the number of regions usually grows at least quadratically with levels.

Now let us introduce a general formulation of [9] for $n$ regions. Let $M_i$ be the $3 \times 4$ non null components of the log of the affine transformation and $\tilde{v}(x)$ be the polyaffine velocity field model:

$$ \tilde{v}(x) = \frac{\sum_{i=1}^{n} w_i(x) \cdot M_i \cdot x}{\sum_{i=1}^{n} w_i(x)} \quad \text{with} \quad \log \left( \begin{bmatrix} A_i & t_i \\ 0 & 1 \end{bmatrix} \right) = \begin{bmatrix} M_i \\ 0 \end{bmatrix}, \qquad (1) $$

where $x, w_i, n$, are spatial position, weight for region $i$ and number of regions, respectively. Now given $v_c$, the correspondence velocity field computed by the demons (not regularized), we can solve for $M_i$ using linear least squares, i.e. minimizing $C_{\text{poly}}(M_1, \dots, M_n) = \int \lambda(x) \|v_c(x) - \tilde{v}(x)\|^2 dx$, where $\lambda$ is a binary mask indicating background voxels. This problem has an explicit solution given by $\sum_j M_i \cdot \Sigma_{ij} = B_i$, with $\Sigma_{ij} = \int \lambda(x) \cdot w_i(x) \cdot w_j(x) \cdot x \cdot x^T dx$ and
$B_i = \int \lambda(x).w_i(x).v(x).x^T dx$. To estimate $M_i$ we need to solve the system $M_i \Sigma = B$, where $\Sigma$ is symmetrical and thus diagonalizable and the minimal norm solution is given using the pseudo inverse $M = B \Sigma^{-1}$.

## OBBTree-Based Hierarchical Scheme for the Polyaffine Log-Demons

The concept of oriented bounding boxes (OBB) has been used extensively in computer graphics to speed up ray tracing and interference detection computations. In [6] the authors presented a hierarchical version and an algorithm to compute it efficiently. An OBBTree is a hierarchy of OBB's in 3D space. Let us consider first a surface (in our case an iso-surface of the CT image). The algorithm computes OBB's via principle component analysis of the vertex coordinates, which give the orientation (principle component directions) and the extent (outmost point on the principle component). A refinement to avoid bias towards densely populated patches is to sample the convex hull of the vertex coordinates and approximate the analytic surface by a linear sum of all triangle areas. There are two ways of calculating the hierarchy, bottom-up and top-down. Top-down approaches start with all vertices and subdivide the points into two groups at every subsequent hierarchical level, whereas bottom-up approaches start by assigning one box per vertex and combine vertices until one box contains all vertices. We used the top-down approach, which divides the vertices into two groups by projecting the vertex coordinates onto the principle components, and using the mean point as the group boundary. The algorithm stops once there are no more divisions possible along any component. For images, the point set could result from a random sampling weighted by the importance of points (e.g. norm of the gradient).

The weights $w_i(x)$ are defined on the OBBTree using multidimensional Gaussian functions as follows,

$$w_i(x) = \exp(-0.5(x - \bar{x}_i)^T [R_i \cdot \text{diag}(\sigma_{i1}^2, \sigma_{i2}^2, \sigma_{i3}^2) \cdot R_i^T]^{-1} (x - \bar{x}_i)) \quad (2)$$

where $\bar{x}_i$ and $R_i$ are center point and orientation of the ith OBB, and $\sigma_{i1}, \sigma_{i2}, \sigma_{i3}$ are the extent of the weights along each axis of the OBB. The parameter $\sigma_{ij}^2$ can be set by the user to enforce different smoothing behavior between regions.

To obtain an intuitive understanding of the smoothing parameter $\sigma_{ij}^2$ we introduce the concept of relative regional mass (RRM). The RRM defines the relative mass for the most prominent weight in a given region: $RRM_i = \frac{\int_{\Omega_i} w_i(x) dx}{\int_{\Omega} w_i(x) dx}$, where $\Omega_i$ and $\Omega$ are the volume given by the ellipsoid that fits inside the ith OBB and the volume of the entire image, respectively. In the special case of RRM equal to 1, there is no overlap between regions. In Fig. 2 (right), a range of possible values are shown with corresponding mean squared error and harmonic energy. In practice, we adjust $\sigma_{ij}^2$ to obtain the same RRM for all regions.

The integration of our hierarchical model into the demons (after rigid alignment) is shown in Algorithm 1. The OBBTree-based polyaffine log-demons regularizes the velocity fields during each update step. The final result is a velocity field that includes all registration steps allowing for statistics on diffeomorphisms
<figure><img src="image_2.png" /><figcaption>Fig. 1: Tree of log affine transformations for the first two levels.</figcaption></figure>

as described in [2]. In addition we obtain $M_i^l$, which is the $i$th log affine transformation at level $l$, providing us with a low-parametric representation of the anatomical structure and allowing further hierarchical statistical analysis and modeling. The hierarchical structure is shown in Fig. 1. Since we are working in the log-Euclidean framework we can subtract the previous log transform from the current level to obtain $N_i^l$, which describes the remaining transformation at that level. The importance of such a hierarchy lies in its power of decomposing features into different scales, this will be elaborated in the next section by performing a hierarchical PCA for mandibles.

**Algorithm 1 OBBTree Log-Demons (OBB-LD)**

Sequentially register OBBTree levels $l = 0, \dots, k$

- Initialize demons with previous level $v^l = v^{l-1}$ (for starting level $v^0 = 0$)

- Iterate until convergence

* Compute correspondence stationary velocity field $v_c^l$ (not regularized)

* Solve linear least square problem $M^l = B^l \cdot \Sigma^{l+}$

* Combine $M_1^l, \dots, M_n^l$ using the polyaffine model (1) to $\tilde{v}^l$

* Let $v^l = \tilde{v}^l$

## 3 Experiments on Mandible CT's

To evaluate the performance of our new method we conducted registrations on 47 CT images of mandibles. After rigid alignment of the OBB at level 0, we performed an OBBTree-based registration with 0 to 6 levels (L0 to L6), each level being initialized by the previous one. Additionally, we initialized two standard log-demons with Gaussian regularization (G0 and G6) at the first and last level (L0 and L6), respectively. Since L0 is one region, it is equivalent to a standard affine initialization.

Fig. 2 (left) shows three important results: (1) Decrease of mean squared error (MSE) with increasing amount of regions (levels). (2) G0 shows higher median, variability and more outliers than G6, which indicates the improvements in terms of robustness, when used as an initialization for standard demons. (3) L6 shows similar performance as G0 and G6, this indicates that by using $2^6 = 64$ regions, we can model a standard demons (relative degrees of freedom: three times the number of voxels). For example, given the dimensions of our mandible CT's our method uses only ≈ 0.01% degrees of freedom to reach the same MSE compared to a standard demons.
<figure><img src="image_3.png" /><figcaption>Fig. 2: (left) Boxplot of mean squared error (MSE) of intensities for each level. G0 and G6 are standard log-demons initialized with L0 and L6, respectively. L0 is one region only, i.e. standard affine initialization. The red crosses indicate outliers. (right) The relative regional mass (RRM) represents the amount of smoothing between regions, where 1 means no overlap between regions.</figcaption></figure>

<figure><img src="image_4.png" /></figure>

Fig. 2 (right) depicts MSE and harmonic energy as a function of different RRM's. The MSE converges at RRM=0.74, whereas the harmonic energy increases with increasing RRM values. For the experiments we chose a RRM of 0.74 to obtain the best MSE with highest smoothness.

Fig. 3 shows 5-level OBBTree, weight correlation structure and the first principle component for each level. The weight correlation structure is computed as $Γ_{ij} = Σ_{ii}^{-1/2} · Σ_{ij} · Σ_{jj}^{-1/2}$, and $Γ_{ij}$ is decomposed using singular value decomposition to extract the major axis of correlation. The correlation of the major axis is colored coded from cool (blue=0.4) to warm (red=1). The graph structure clearly reveals the intrinsic underlying dimensionality of the object at each scale, going from a curve to a ribbon and finally in some areas locally to a 3D volume. It would be interesting to study if this could be a robust alternative to the medial axis or surface representation.

The per level PCA can be interpreted as follows: (L0) global scaling; (L1) thickness; (L2) reorientation in the region of the masseter; (L3) relative displacement of condyles and coronoid processes; (L4) change in teeth region; (L5) change in back teeth region. This gives a visual validation of the usefulness of per level hierarchical statistical analysis, clearly distributing feature to different scales.

# 4 Conclusions

In this work we presented a geometry-aware multiscale approach for registration. We showed that our method significantly robustifies the standard demons and increases its repeatability when the full range of scales is used. When stopping at a reasonable scale (here only 64 components), it even performs similarly while
<figure><img src="image_5.png" /><figcaption>Fig. 3: (OBB) 6 levels of the OBBTree visualized on the reference mandible, starting with level 0 up to level 5. (Correlation) Major axis of correlations between regions, color coded from cold (blue=0.4) to warm (red=1.0). (PC1 front and side view) Red and white surface are -2 and +2 standard deviations from the mean.</figcaption></figure>
presenting significant reduction of degrees of freedom of the registration. We further visualized a per hierarchical level PCA, which suggested a clear division of shape features into the different levels, allowing for a better interpretation, as opposed to a PCA on the entire field with potential mixtures of features.

This new method opens a large number of potential research opportunities: How to analysis tree-like data objects? How to enforce sparsity on these objects? What is the correlation between levels? How to measure anatomical meaningfulness of deformations quantitatively (in this work a qualitative approach visualizing the modes of variation is used)? How to measure the explained variance in trees? All these questions will be addressed in future work.

**Acknowledgements** This work has been funded by the Swiss National Science Foundation.

## References

1. Arsigny, V., Commowick, O., Ayache, N., Pennec, X.: A Fast and Log-Euclidean Polyaffine Framework for Locally Linear Registration. J. Math. Imaging Vis. 33(2), 222-238 (2009)

2. Arsigny, V., Commowick, O., Pennec, X., Ayache, N.: A Log-Euclidean Framework for Statistics on Diffeomorphisms. In: Larsen, R., Nielsen, M., Sporring, J. (eds.) MICCAI 2006, LNCS, vol. 4190, pp. 924-931. Springer, Heidelberg (2006)

3. Arsigny, V., Pennec, X., Ayache, N.: Polyrigid and Polyaffine Transformations: A Novel Geometrical Tool to Deal with Non-rigid Deformations. Application to the Registration of Histological Slices. Med. Image Anal. 9(6), 507-523 (2005)

4. Coletti, Ord, R., Liu, X.: Mandibular Reconstruction and Second Generation Locking Reconstruction Plates: Outcome of 110 Patients. Int. J. Oral Max. Surg. 38(9), 960-963 (2009)

5. Commowick, O., Arsigny, V., Isambert, A., Costa, J., Dhermain, F., Bidault, F., Bondiau, P.Y., Ayache, N., Malandain, G.: An Efficient Locally Affine Framework for the Smooth Registration of Anatomical Structures. Med. Image Anal. 12(4), 427-441 (2008)

6. Gottschalk, S., Lin, Manocha, D.: OBBTree: A Hierarchical Structure for Rapid Interference Detection. In: SIGGRAPH. pp. 171-180 (1996)

7. Kozic, N., Weber, S., Büchler, P., Lutz, C., Reimers, N., González Ballester, M.A., Reyes, M.: Optimisation of Orthopaedic Implant Design Using Statistical Shape Space Analysis Based on Level Sets. Med. Image Anal. 14(3), 265-275 (2010)

8. Metzger, M.C., Vogel, M., Hohlweg-Majert, B., Mast, H., Fan, X., Rüdell, A., Schlager, S.: Anatomical Shape Analysis of the Mandible in Caucasian and Chinese for the Production of Preformed Mandible Reconstruction Plates. J. Cranio. Maxill. Surg. (2010)

9. Seiler, C., Pennec, X., Ritacco, L., Reyes, M.: Femur Specific Polyaffine Model to Regularize the Log-Domain Demons Registration. In: SPIE Medical Imaging. vol. 7962 (2011)

10. Vercauteren, T., Pennec, X., Merchant, A., Ayache, N.: Symmetric Log-Domain Diffeomorphic Registration: A Demons-Based Approach. In: Metaxas, D., Axel, L., Fichtinger, G., Székely, G. (eds.) MICCAI 2008, LNCS, vol. 5241, pp. 754-761. Springer, Heidelberg (2008)