Christof Seiler

Xavier Pennec

Mauricio Reyes

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

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.

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 [START_REF] Coletti | Mandibular Reconstruction and Second Generation Locking Reconstruction Plates: Outcome of 110 Patients[END_REF]. 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 [START_REF] Metzger | Anatomical Shape Analysis of the Mandible in Caucasian and Chinese for the Production of Preformed Mandible Reconstruction Plates[END_REF], and for population-based orthopaedic femural plate design [START_REF] Kozic | Optimisation of Orthopaedic Implant Design Using Statistical Shape Space Analysis Based on Level Sets[END_REF]. 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 population shape variability makes the registration a challenging task, which calls for improvements in non-rigid registration for orthopeadic research. Whenever additional 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 results. In this paper we present a novel method, which incorporates such a description of the geometry of the anatomy into the registration process. The method builds on the work of [START_REF] Seiler | Femur Specific Polyaffine Model to Regularize the Log-Domain Demons Registration[END_REF], where the authors presented a polyaffineregularized 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 accuracy and robustness can be improved by incorporating a locally affine and anatomy-aware hierarchical scheme. Second, we demonstrate that the complexity of the target deformation can be reduced, which is of great interest in statistics 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 provides 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.

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 implicitly 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 logdemons.

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 [START_REF] Vercauteren | Symmetric Log-Domain Diffeomorphic Registration: A Demons-Based Approach[END_REF]. 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 [START_REF] Arsigny | A Log-Euclidean Framework for Statistics on Diffeomorphisms[END_REF]. This property explains the denomination of log-domain (or simply log-demons) registration.

Polyaffine Registration Polyaffine transformations were introduced in [START_REF] Arsigny | Polyrigid and Polyaffine Transformations: A Novel Geometrical Tool to Deal with Non-rigid Deformations Application to the Registration of Histological Slices[END_REF] 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 [START_REF] Commowick | An Efficient Locally Affine Framework for the Smooth Registration of Anatomical Structures[END_REF] 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 [START_REF] Arsigny | A Fast and Log-Euclidean Polyaffine Framework for Locally Linear Registration[END_REF] to work with stationary vector fields. This new formulation suggested that polyaffine transformations could be compatible in some sense with the log-demons. In [START_REF] Seiler | Femur Specific Polyaffine Model to Regularize the Log-Domain Demons Registration[END_REF] 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 [START_REF] Seiler | Femur Specific Polyaffine Model to Regularize the Log-Domain Demons Registration[END_REF] for n regions. Let M i be the 3 × 4 non null components of the log of the affine transformation and ṽ(x) be the polyaffine velocity field model:

ṽ(x) = n i=1 w i (x).M i .x n i=1 w i (x) with log A i t i 0 1 = M i 0 , (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 poly (M 1 , . . . , M n ) = λ(x) v c (x) -ṽ(x) 2 dx, where λ is a binary mask indicating background voxels. This problem has an explicit solution given by j M i .Σ ij = B i , with Σ ij = λ(x).w i (x).w j (x).x.x T dx and

B i = λ(x).w i (x).v(x).
x T dx. To estimate M i we need to solve the system M.Σ = B, where Σ is symmetrical and thus diagonalizable and the minimal norm solution is given using the pseudo inverse M = B.Σ + .

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 [START_REF] Gottschalk | OBBTree: A Hierarchical Structure for Rapid Interference Detection[END_REF] 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 -xi ) T .[R i . diag(σ 2 i1 , σ 2 i2 , σ 2 i3 ).R T i ] -1 .(x -xi ) , (2)
where xi and R i are center point and orientation of the ith OBB, and σ i1 , σ i2 , σ i3 are the extent of the weights along each axis of the OBB. The parameter σ 2 ij can be set by the user to enforce different smoothing behavior between regions.

To obtain an intuitive understanding of the smoothing parameter σ 2 ij 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 = Ω i wi(x)dx Ω wi(x)dx , where Ω i and Ω 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 σ 2 ij 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

N 0 1 = M 0 1 N 1 1 = M 1 1 -M 0 1 N 1 2 = M 1 2 -M 0 1
Fig. 1: Tree of log affine transformations for the first two levels.

as described in [START_REF] Arsigny | A Log-Euclidean Framework for Statistics on Diffeomorphisms[END_REF]. In addition we obtain M l i , which is the ith 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 l i , 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, . . . , 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 l c (not regularized) • Solve linear least square problem M l = B l .Σ l + • Combine M l 1 , . . . , M l n using the polyaffine model (1) to ṽl • Let v l = ṽl

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. 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 = Σ -1/2 ii .Σ ij .Σ -1/2 jj
, 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.

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

  demons (L0 to L6), Gaussian demons (G0 and G6)Distribution of registration error per levelmean squared error (MSE)

Fig. 2 :

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

Fig. 3 :

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

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


Christof Seiler

Xavier Pennec

Lucas Ritacco

Mauricio Antonio

Reyes Aguirre

Lucas

Mauricio Reyes

Femur Specific Polyaffine Model to Regularize the Log-Domain Demons Registration

Keywords: Registration, Regularization, Log-Domain Demons, Polyaffine, Femur, Allotransplantation

  L'archive ouverte pluridisciplinaire

INTRODUCTION

Tumor excision is the primary treatment of aggressive or recurrent benign bone tumors and malignant bone sarcomas. This requires an invasive surgical intervention that entails a residual bone defect, which can be reconstructed with a fresh frozen bone allotransplantation. In orthopedics, recent improvements in diagnostic and therapeutic techniques have produced an increase of the patient survival as well as a reduction of complication rate. [START_REF] Muscolo | Use of Distal Femoral Osteoarticular Allografts in Limb Salvage Surgery[END_REF] It has been shown that the selection of bone allograft in terms of shape and size is crucial to prevent changes in joint mobility and load distribution, which can lead to joint fractures and early joint degeneration. [START_REF] Enneking | Retrieved Human Allografts: A Clinicopathological Study[END_REF] However, current selection approaches are very time-consuming, mostly based on manual measurements performed directly on the bones or on three-dimensional models reconstructed from images.

Ideally, an automatic allograft selection system could perform shape comparisons between the patient's anatomy and the complete bone bank. In this case, it would be necessary to allocate extra time to the clinical work flow in order to segment the patient's anatomy and register it to the existing bone bank.

Practically, to avoid preprocessing, hence reducing preoperative planning time, the authors in Ritacco et al. [START_REF] Ritacco | Three-Dimensional Morphometric Analysis of the Distal Femur: A Validity Method for Allograft Selection Using a Virtual Bone Bank[END_REF] proposed a simplified protocol for bone size characterization. The so-called ABC protocol is based on well-known anatomical landmarks which can be extracted from the CT image and used in a selection criteria. The protocol defines 6 landmarks Further author information: (Send correspondence to C. Seiler) C. Seiler: E-mail: christof.seiler@istb.unibe.ch on the distal femur: Transepicondyle distance (A), anterior-posterior distance in the medial condyle (B) and in the lateral condyle (C). One possible automation of this protocol requires two steps. First, CT images of all donor bones are segmented and registered to a preselected reference. The ABC measurements can now be performed on the reference and correspondences in the other images can be extracted and stored. Second, the surgeon identifies manually the ABC landmarks on the patient CT image and uses the previously stored ABC measures to select the best match.

In order to automatize the procedure, we first tried to register non-linearly the CT image to the candidate bones using a demons registration approach. [START_REF] Vercauteren | Symmetric Log-Domain Diffeomorphic Registration: A Demons-Based Approach[END_REF] However, we observed that the direct use of this method failed to capture torsion and rotation in head and condyle regions of the femur bone. We propose in this paper to model the expected deformations in the specific case of femur bone registration considering torsion and rotation in head and condyles, and to use this femur specific polyaffine model as an adapted regularization prior for the demons registration in replacement of the uninformed Gaussian smoothing.

We build on the concept of polyaffne registration, introduced in Arsigny et al. [START_REF] Arsigny | Polyrigid and Polyaffine Transformations: A Novel Geometrical Tool to Deal with Non-rigid Deformations Application to the Registration of Histological Slices[END_REF] and fully developed into an efficient and practical framework in Commowick et al. [START_REF] Commowick | An Efficient Locally Affine Framework for the Smooth Registration of Anatomical Structures[END_REF] We initialize the demons with a polyaffine registration of three regions (head, shaft and condyles) using block matching [START_REF] Ourselin | Block Matching: A General Framework to Improve Robustness of Rigid Registration of Medical Images[END_REF] and regularize the subsequent iterations with a polyaffine model integrated into the symmetric demons algorithm. As a last step, we relax the polyaffine constraint to allow more general local deformations regularized using a Gaussian smoothing.

The pipeline was developed using open source libraries: VTK, ITK and Qt (UI framework). In the following, we develop the methodology and show comparison results between multiresolution demons with Gaussian regularization and our approach considering demons with the new polyaffine regularization term.

A POLYAFFINE LOG-DOMAIN DEMONS FOR THE FEMUR

Our main contribution in this paper is to show that the demons algorithm can be improved in specific cases with an appropriate model. The main idea of the proposed algorithm is to introduce a polyaffine regularization of the velocity field. The polyaffine model divides the femur in three regions to capture torsion and rotation in head and condyles. In the following sections we review the two methods that we later combine: Log-domain and polyaffine registration. We show how to add our specific femur model during the regularization step of demons through an estimate of the polyaffine parameters given a velocity field. This estimate replaces the standard Gaussian smoothing.

Log-Domain Registration

To setup correspondences between anatomical images, a set of images are registered to a reference. We use the novel symmetric diffeomorphic demons registration approach described in Vercauteren et al. [START_REF] Vercauteren | Symmetric Log-Domain Diffeomorphic Registration: A Demons-Based Approach[END_REF] What is new in this registration framework is the efficient optimization in the log-domain. As a consequence, the results of the registration are 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. 8

Polyaffine Registration

In Commowick et al. [START_REF] Commowick | An Efficient Locally Affine Framework for the Smooth Registration of Anatomical Structures[END_REF] the authors introduced an iterative framework to fuse regional affine transformation to a global transformation. In a first step, an affine registration is performed for each region, second, the affine transformation are regularized, and finally a global transformation is calculated through a weighted average of all regions. The averaging of affine transformations can be done on the matrix logarithms to ensure invertibility of the final transformation. The regions are manually defined on the reference prior to the registration procedure.

Femur Polyaffine Model

For the polyaffine femur model we define three regions on the reference image. The middle section of the femur is defined as shaft region. The shaft region is centered in the middle of the reference femur and covers half of the volume of the image. The distal and the proximal end are condyle and head region, respectively. All images to be registered are anisotropically scaled to fit the size of the reference femur before starting with the polyaffine registration.

Let Mi be the 3 × 4 non null components of the log of the affine transformation log Ai ti

0 1 = Mi 0 , ( 1
)
where i = 1 is an affine transform for condyles, i = 2 a rigid transform for the shaft and i = 3 an affine transform for the head region. We use the following weight function to connect the transformations:

wi(x) = f e T z .(x -xi) ||x3 -x1|| , ( 2
)
where f (•) is a Gaussian function of the spatial position x ∈ R 3 with μ = 0 and σ = 0.125. The Gaussian parameters are defined heuristically with the aim to find the best compromise between smooth transitions from one region to the next and considering each region individually with the least amount of influence by others. The center positions of each region x1, x2 and x3 are set to the center of the transverse plane. Along the transversal axis the positions are set to achieve an equal weighting at the boundary of two adjacent regions. The weights are normalized to wi(x) = 1. Thus the velocity field predicted by our model is v poly (x) = w1(x).M1.x + w2(x).M2.x + w3(x).M3.x.

(3)

Estimation of Femur Polyaffine Model

Here we introduces the estimation of our polyaffine model which will be needed in the next section for the demons regularization. We show how to find a closed form solution on logarithms of affine transformations. Given a velocity field, the log affine parameters of (3) can be estimated by solving a linear least square problem. The least square error with respect to the velocity field observed is given by:

C(M1, M2, M3) = d(x)||w1(x).M1.x + w2(x).M2.x + w3(x).M3.x -v(x)|| 2 dx, ( 4
)
where d(x) is a binary variable set to 0 in the background. Using the Frobenius inner product which generalizes the dot product to matrices, ||W || 2 = Tr(W.W T ), we obtain the following directional derivative:

∂W C = lim →0 C(M1 + W1, M2 + W2, M3 + W3) -C(M1, M2, M3) = d(

x). Tr[(w1(x).W1.x + w2(x).W2.x + w3(x).W3.x)(w1(x).M1.x + w2(x).M2.x + w3(x).M3.xv(x)) T ] = Tr[W1. d(x).w1(x).x.(w1(x).M1.x + w2(x).M2.x + w3(x).M3.xv(x))

T ] + same for W2 and W3.

At the optimum, the directional derivative should be null in each direction W , i.e. ∂W C = 0 for all matrix W. As we have Tr(W.A) = 0 for all W iff A = 0, we end up with the system

j Mi.Σij = Bi, (5)
with

Σij = d(x).wi(x).wj(x).x.x T dx, Bi = d(x).wi(x).v(x).x T dx.

To estimate M1, M2 and M3 we need to solve the following system:

M1 M2 M3 ⎡ ⎣ Σ11 Σ12 Σ13 Σ12 Σ22 Σ23 Σ13 Σ23 Σ33 ⎤ ⎦ = B1 B2 B3 , ( 6
)
where Σ is symmetrical and thus diagonalizable and the solution is given by:

M = B.Σ -1 . ( 7
)
This system can be solved very efficiently since the matrix Σ is only of dimension 12 × 12. Σ does not change for different velocity field observations as long as the weighting functions are fixed, hence only B needs to be recomputed. This fact will be helpful during the iterative regularization of demons which is described in the next section.

Algorithm 1 Relaxed Femur Polyaffine Log-domain Demons (RFPLD)

• Affine registration of condyles, shaft and head region

• Combine affine transformations to one velocity field v using (3)

• Initialize demons with v

• Iterate until convergence:

-Compute the demons forces u

-Let v ← Z X (v, u) (Baker-Campbell-Hausdorff approximation 4 )
-For polyaffine regularization let v ← solve ( 6) for v and combine M 1 , M 2 and M 3 using (3)

• Iterate until convergence:

-Compute the demons forces u

-Let v ← Z X (v, u) -For diffusion-like regularization let v ← K diff v (Gaussian smoothing)

Regularization with a Femur Polyaffine Model

Since the demons algorithm decouples the registration problem in force computations using optical flow equations and a smoothing step it is straight forward to introduce a new regularization by replacing the latter. Our new polyaffine regularized demons estimates the polyaffine femur model at each step of the demons registration, see Algorithm 1, and thus connects the non-linear and image intensity based computation of forces with the femur model.

VALIDATION WITH ABC GROUND TRUTH AND VISUALIZATION OF TORSION AND ROTATION

We compare three different registration methods: First, a standard multiresolution demons with Gaussian smoothing for the regularization, referenced hereafter with Standard Log-domain Demons (SLD). Second, the newly developed polyaffine regularized demons, referenced hereafter with Femur Polyaffine Log-domain Demons (FPLD). Third, the second registration method with a subsequent relaxation of the regularization term to Gaussian smoothing as described in Algorithm 1, referenced hereafter with Relaxed Femur Polyaffine Log-domain Demons (RFPLD). We conduct a task specific validation of the newly developed regularization, testing the performance of all three registration methods on ABC measurements for bone tumor replacement. The ground truth of the measurements was obtained by manually selecting the landmarks on 3D reconstructed surface models. Two medical doctors performed these measurements independently and one expert repeated the measurements after a few days to yield the interobserver variability.

The results in Fig. 1 show that for the SSD only a small difference is visible for SLD and RFPLD, whereas FPLD is clearly higher. The smoothness decreases slightly in our new approach. The newly developed method RFPLD provides smaller ABC measurement errors compared to the standard demons. In two of three cases the FPLD shows smaller ABC measurement errors than the old approach. The accuracy of the A and B measurements (and in a smaller measure for C) are of the same order of magnitude than the intra and inter-rater errors for our relaxed polyaffine model.

To analyze the obtained results we perform a principle component analysis (PCA) on the velocity fields obtained by the three registration methods. In Fig. 2 the variance along the first three PCs of the six landmarks are visualized with curves. Each curve corresponds to the flow of a landmark point by integrating the velocity field starting at the mean (t = 0) in both direction with end points t = 3σ and t = -3σ. These curves are also called streamlines and used in fluid dynamics. Since the exponential map and streamlines of a flow are equivalent, it is an intrinsic visualization tool to understand the spatial deformations.

Similarly, to show the torsion and rotation variance in the femur, streamlines were calculated for the femoral head (Fig. 3). Two different views were chosen to best visualize the variability patterns. It can clearly be seen that the first PC describes rotations along the transversal axis, the second PC a change of the caput collum diaphysis angle and the third PC a S-shaped deformation pattern with a stretching of the femoral head.

DISCUSSION AND CONCLUSIONS

In this work we introduced a femur specific polyaffine model to replace the regularization step of the demons algorithm.

The goal was to improve the registration in terms of anatomical correspondences to enable automatic ABC measurements for allograft selection.

Comparisons to the standard demons clearly indicate the power of the new method to capture the anatomical variability. This was visualized in Fig. 2 with the much higher variances at landmark positions. The validation in Fig. 1 shows that our method provides higher or equal SSD but better ABC correspondences suggesting that the standard demons finds a local minimum. Instead our method guides the demons to a minimum that is anatomically more meaningful in terms of correspondences. The ABC measurements give us a ground truth that is clinically relevant. The accuracy of the measurements are of the same order of magnitude than the intra and inter-rater errors for our relaxed polyaffne model, which shows that we are clearly reaching the performances of the human raters. This indicates that it is important to include our intuitive understanding of the anatomy and shows that in specific cases the demons can be improved by incorporating a simple a priori model like our femur specific polyaffine model.

Future work includes changing the resampling of the images from 1.25 mm to a smaller voxel spacing, to further improve the accuracy thus potentially being more accurate than the human raters. The demons provides a straight forward and fast vehicle for the connection of models with medical images since it divides force computations and the regularization step. This opens up new ways of model-based registration with the demons algorithm for other anatomical structures.

Figure 1 .Figure 2 .

 12 Figure 1. Boxplots comparing the three registration methods: SLD (•1), FPLD (•2) and RFPLD (•3). (•4) and (•5) are interobserver and intraobserver errors. (S•), (E•), (A•), (B•) and (C•) are SSD, harmonic energy and the three ABC distances.

Figure 3 .

 3 Figure 3. Visualization of the variance of the velocity fields through streamlines. The first three PCs of RFPLD is shown for randomly picked points which are integrated from t = -3σ to t = 3σ: (a) PC1, (b) PC2 and (c) PC3. The color describes the length of the trajectory.