Demian Wassermann

James Ross

George Washko

William M Wells

Raul San

Jose -Estepar

Deformable Registration of Feature-Endowed Point Sets Based on Tensor Fields

The main contribution of this work is a framework to register anatomical structures characterized as a point set where each point has an associated symmetric matrix. These matrices can represent problem-dependent characteristics of the registered structure. For example, in airways, matrices can represent the orientation and thickness of the structure. Our framework relies on a dense tensor field representation which we implement sparsely as a kernel mixture of tensor fields. We equip the space of tensor fields with a norm that serves as a similarity measure. To calculate the optimal transformation between two structures we minimize this measure using an analytical gradient for the similarity measure and the deformation field, which we restrict to be a diffeomorphism. We illustrate the value of our tensor field model by comparing our results with scalar and vector field based models. Finally, we evaluate our registration algorithm on synthetic data sets and validate our approach on manually annotated airway trees.

Introduction

Point-set based representations arise in a wide variety of medical imaging applications. Examples include the extraction of structures like blood vessels and airways [START_REF] Aylward | Initialization, noise, singularities, and scale in height ridge traversal for tubular object centerline extraction[END_REF][START_REF] Kindlmann | Sampling and Visualizing Creases with Scale-Space Particles[END_REF]. The ability to register two different point-sets representing the same anatomical structure is critical to enable populationbased studies. It is also important for tracking longitudinal characteristics of the structure of interest. Non-rigid pointset registration algorithms exist (e.g. [START_REF] Tsin | A Correlation-Based Approach to Robust Point Set Registration[END_REF][START_REF] Jian | Robust Point Set Registration Using Gaussian Mixture Models[END_REF]); however, these algorithms represent structures as a collection of points in R 3 neglecting valuable information regarding the shape of the structure. This representation was recently improved by the currents model [START_REF] Durrleman | Statistical Models of Sets of Curves and Surfaces based on Currents[END_REF][START_REF] Vaillant | Statistics on diffeomorphisms via tangent space representations[END_REF] enriching the information of each point with a vector. Currents has proved to be a extremely useful in registration situations that involve orientable surfaces, where the vector is the surface normal; and curves, where it is the tangent to the line. Despite these successes, the Currents model has shortcomings in cases where a vector space is not powerful enough to represent the shape in-formation of the structure of interest. An example of this are airways, in which it is desirable to represent the direction of the structure at a point as well as its thickness.

The main contribution of this work is a framework to perform registration of anatomical structures characterized as a point set where each point has an associated symmetric matrix. These matrices may represent problem-dependent characteristics of the registered structure. Examples of these are airways in which the matrices represent the orientation and thickness of the near-tubular structure The proposed registration framework relies on a dense tensor field representation, which we represent sparsely as a kernel mixture of tensor fields. We equip the tensor space with a norm that provides a similarity measure. By optimizing this measure, or matching criterion, we calculate the optimal transformation between two structures.

Borrowing tools from differential geometry and matrix calculus we derive an analytical gradient for the matching criterion and thus, the velocity field. Using our closed-form gradient and the transform field representation of one parameter subgroups of diffeomorphisms [START_REF] Ashburner | A fast diffeomorphic image registration algorithm[END_REF][START_REF] Arsigny | A Log-Euclidean Framework for Statistics on Diffeomorphisms[END_REF], the resulting registration algorithm yields a diffeomorphic transform. To illustrate the value of the tensor representation, we compare our results with the scalar-based [START_REF] Jian | Robust Point Set Registration Using Gaussian Mixture Models[END_REF] and vector based [START_REF] Vaillant | Statistics on diffeomorphisms via tangent space representations[END_REF] representation methods. We evaluate our registration algorithm on synthetic data sets and illustrate the utility of our approach on manually annotated airway trees. All algorithms presented in this paper can be downloaded at http://github.com/demianw/pyMedImReg-public

Methods

Feature Field Representation and Distance

Let us represent an anatomical structure as a point set where each point is endowed with a feature. We will call these enriched points particles, and represent a set as P = {(p i , F i )} where p i ∈ Ω is the spatial location of the particle and F i is a feature representing characteristics like direction or volume. If the features belong to a space F admitting linear combination, we can represent P as a feature field w:

Ω →F w(x) = i k(x, p i )F i (1)
where K : Ω × Ω →R is a kernel function. Moreover, if F is equipped with an inner product •, • F , this induces an inner product operation between feature fields:

w 1 , w 2 W = Ω w 1 (x), w 2 (x) F dx = Ω ij k(x, p 1 i )F 1 i , k(x, p 2 j )F 2 j F dx = ij K(p 1 i , p 2 j ) F 1 i , F 2 j F , with K(x, y) = Ω k(η, x)k(η, y)dη.
(

) 2
This inner product operation in eq. ( 2) induces a distance between objects in the feature field space,

w 1 -w 2 2 W = w 1 2 W + w 2 2 W -2 w 1 , w 2 W , w 2 W = w, w W (3)
This distance measures the similarity between the two feature fields globally, without assuming point correspondence between the anatomical structures. Particular cases of this are the Gaussian mixture model (GMM) and the discrete setting of the currents model [START_REF] Jian | Robust Point Set Registration Using Gaussian Mixture Models[END_REF][START_REF] Durrleman | Statistical Models of Sets of Curves and Surfaces based on Currents[END_REF]. In GMM, the features F are a scalar, usually 1, while in currents they are 3D vectors and •, • F is the inner product in R 3 . In both cases, the kernel K is a Gaussian function. In this work, we extend this representation to include the feature space of symmetric positive semidefinite (SPS) matrices, then w(x) becomes a tensor field. This feature representation is powerful enough to consider directionality information of each particle as well as orientation and volumetric information. Let us define the feature space F ⊂ Sym 3 ≥0 with the inner product F, G F = trace(FG T ) which is the Frobenius inner product. These tensors have a graphical representation as ellipsoids where the axes are aligned with the eigenvectors and have the length of the eigenvalues. Our choice of inner product measures the similarity between two symmetric matrices taking in account volume, calculated as the determinant, and angle. It enables us to use the rich field of matrix calculus to derive our registration algorithm. We show examples of these representations in fig. 1 where we illustrate how the tensor field represents orientation, inter-point distance and width while the vector field or currents approach represents inter-point distance and direction only. Our choice of inner product measures the similarity between two symmetric matrices taking in account determinant and angle. It enables us to use the rich field of matrix calculus to derive our registration algorithm. We show examples of these representations in fig. 1 where we illustrate how the tensor field represents orientation, inter-point distance and width while the vector field, or currents, model represents only inter-point distance and direction.

Tensor Field Deformation

To implement a registration algorithm, we need to be able to transform tensor fields by a smooth map. Let φ : R 3 → R 3 be a smooth map, transforming w(x) by φ can be formulated through a pullback operation φ * :

φ * w = w(φ(x)) = i k(φ(x), p i ) φ -1 * F i . (4)
In the case where we apply the map φ only to discrete points p i , then

φ * w = w(φ(x)) = i k(x, ψ(p i ))ψ * F i ( 5
)
where ψ is φ -1 with respect to the kernel function K and ψ * F i is the action of ψ on F i . In the case that the feature F i is a vector f ∈ R 3 , the action ψ on f is defined as

D p i ψ T f , with D p i ψ the Jacobian of ψ(x) at p i [9]
. However, in our case F i are symmetric positive semidefinite (SPS) matrices F. Then, we derive the action ψ by decomposing the SPS matrix as F = 3 j=1 f j [f j ] T , and applying the vector action

ψ * F = j=1...3 (D T p i ψf j )(D T p i ψf j ) T = D T p i ψFD p i ψ. (6)
In this work, we picked, in agreement with [START_REF] Durrleman | Statistical Models of Sets of Curves and Surfaces based on Currents[END_REF], k(x, y) = r( xy 2 ) with r a radial kernel.

Tensor Field Registration

Having a representative space for our anatomical structures, we now derive a registration algorithm. Given two feature fields w 1 and w 2 , we look for the transform φ minimizing:

E(φ) = w 1 -φ * w 2 2 W + Reg(φ) (7)
where the first term of E(φ), or data attachment term, quantifies the similarity between the feature fields and the second one is a regularization term on the transform. To minimize eq. ( 7) we use a gradient-based algorithm. Hence, assuming that φ(x) = φ(x; θ 1 , . . . , θ M ) is a parametric transform, we need to obtain the gradient of E(φ) with respect to each parameter. As we show in appendix A, using properties of the inner product space of tensor fields we obtain the gradient of the data attachment term with respect to the transform parameter θ: In the case that our features are tensors, we use eq. ( 5) and matrix calculus properties (see appendix A) to derive

∂ θ φ * w 2 -w 1 2 W = 2 ( φ * w 2 ,∂ θ φ * w 2 W -w 1 ,∂ θ φ * w 2 W )
∂ θ φ * w 2 = i ∇k(x, ψ(p i ))∂ θ T ψ(p i ) J T i F i J i + i k(x, ψ(p i )) ∂ θ J T i F i J i + J T i F i (∂ θ J i ) , (8)
with

J i = D p i ψ, ψ = φ -1 .
Diffeomorphic deformation case: If we restrict φ to be a diffeomorphism parameterized by a stationary velocity field [START_REF] Ashburner | A fast diffeomorphic image registration algorithm[END_REF][START_REF] Arsigny | A Log-Euclidean Framework for Statistics on Diffeomorphisms[END_REF] we can define a velocity field u : R 3 →R 3 , and characterize the deformation by the differential equation dφ/dt = u(φ (t) ) subject to φ (0) = Id and t ∈ [0, 1]. To generate the transform φ, we start with φ (0) = Id, and integrate it over unit time to obtain φ = φ (1) . We implement this integration by discretizing the unit time interval into R equal size steps:

φ (1/R) (x) = x+ u(x) R , φ ((i+1)/R) = φ (1/R) •φ (i/R) (9) with i = 2 . . . R.
Larger values of R yield a more accurate approximation to φ. From eq. ( 9) we deduce that the deformation field approximation is uniquely determined by the velocity field u and the number of steps R. If we take a parametrical representation of u, with parameters Θ = {θ i }, the derivative of φ with respect to a parameter θ is:

∂ θ φ (1/R) (x) = ∂ θ u(x) R , ∂ θ φ ((i+1)/R) = ∂ θ φ (1/R) •φ (i/R) +∂ θ φ (i/R) D x φ (1/R)
As our particle-sets are not bound to a specific grid, we represent the velocity field sparsely, through the interpolating function u(x) = j exp(-p ix 2 /σ 2 )α i , where σ controls its smoothness and the parameters of the transform are Θ = {α i , p i ∈ R 3 } i≤N . This representation has advantages of being implementable by mixture of finite support functions and an easy calculation of its derivatives with respect to the parameters Θ and position x. It is worth noting that our tensor field approach is independent of the particular choice of transform provided that it is a diffeomorphism. How to model and parameterize such transforms is a wide area [START_REF] Holden | A Review of Geometric Transformations for Nonrigid Body Registration[END_REF] which we do not address in this work as it is centered on the choice of metric between anatomical structures.

Regarding regularization on eq. ( 7), as we are using a stationary velocity field, we chose a space-only regularizer Reg(φ) = Lu(x) 2 dx with the linear difference operator L = ∆ which ensures the smoothness of the velocity field.

Finally, having described in detail the terms of eq. ( 7) and the gradient of the data attachment term, we minimize eq. ( 7) through a quasi-newton approach, the L-BFGS-B method [START_REF] José | Remark on "algorithm 778: L-BFGS-B: Fortran subroutines for large-scale bound constrained optimization[END_REF]. This gradient-based optimization method uses an approximation of the inverse Hessian of E(φ) to speed up convergence while keeping the memory requirements fixed during the optimization process.

Airway datasets

In the case of our airway datasets, we sample the airway centerline with scale space particles [START_REF] Kindlmann | Sampling and Visualizing Creases with Scale-Space Particles[END_REF]. The particles are equally distributed along the centerline of the airways at a distance s and have an ellipsoidal shape derived from the Hessian matrix at the scale-space location of the airway point. From this particles, define the tensor field: where one eigenvalue of F i is set to the inter-particle distance s/2 with its eigenvector aligned with the extracted airway direction, and the other two have the lumen size at p i as their magnitude. In the limiting where the ellipsoids are infinitely thin (s→0) they become discs modeling the tubular structure of the airways.

w(x) := i exp -(x -p i ) 2 /σ 2

Experiments

Synthetic Experiments Our first synthetic experiment is to show the efficacy of the Frobenius norm to register tensor fields. We generated two synthetic linear transformations and registered them back using the algorithm in section 2.3, only here we take the linear transform φ(x) = Ax + t and no regularization term. We show the results of this experiment in fig. 2 where the red tensor field is being registered onto the green one.

On a second synthetic experiment we tested the efficacy of the deformable registration algorithm and compared the results of using a tensor field against the currents and GMM models. The deformation field parameterization and optimization was kept equal for all approaches in order to test the performance of the choice of representation. To obtain a dataset with a gold standard we produced 100 random deformations of an airway dataset. The initial dataset is inscribed in a box of 190 × 125 × 250mm. We generated smooth random deformations by producing Gaussian ran-Gen 1 2 3 4 5 Position Overlap with Manual Labeling Tens 91% ± 2 90% ± 6 73% ± 5 62% ± 3 61% ± 3 Cur 90% ± 2 91% ± 7 73% ± 2 61% ± 3 63% ± 3 Volume Overlap with Manual Labeling Tens 87% ± 5 75% ± 3 72% ± 5 62% ± 4 60% ± 5 Cur 42% ± 8 33% ± 10 38% ± 7 31% ± 9 28% ± 15

Table 1: Comparison of the Tensor and Current methods in with manually labeled datasets. The performance regarding overlap of labelled regions is comparable between the Tensors and Currents methods. Regarding the volume of the particles our methods outperforms Currents.

dom vector fields with standard deviation γ taking values of 0.1; 1.0; and 10.0mm and convolving them with a Gaussian kernel. Then, we applied each transformation to an airway image transforming the points and the particle shapes according to section 2.2. Finally, we registered each set to the original template and quantified the registration error by taking the mean squared error (MSE) of the position of each registered particle and the original position in the template; its original volume; and the angle with respect to its original direction. In fig. 3 we show: on the left panel ten realizations of randomly deformed airways; center: the original template in white, a deformed case in green and its resulting registration to the template in red; right: the quantitative results for our registration comparison. Our method proved to perform better then GMM and currents: our results show to have a lower MSE in terms of the distance to the original point than the GMM and the current approaches for all noise levels. Human Data Airway Experiments We analyzed 20 volumetric chest CT datasets acquired at full inspiration. The scans were performed either with GE scanner and Siemens scanners. In-plane pixel spacing ranged from 0.54mm to 0.85mm across all scans. The CT were down-sampled with a x4 ratio and the airway tree extracted using scale-space particles [START_REF] Kindlmann | Sampling and Visualizing Creases with Scale-Space Particles[END_REF]. We registered all datasets affinely and then deformably. We show the results of this process in fig. 4 where we gave a different colored to each subject. Were the high level of agreement can be noticed. To validate our results, an expert labeled the airways generations from 1 to 5 (1 being the trachea and 5 the smaller bronchi) [START_REF] Diaz | Airway count and emphysema assessed[END_REF]. Then, we quantified the classification agreement of the overlapping particles, on a per-generation basis, we show these results in table 1.

We show the results of our registration algorithm in the top right image of fig. 4. Our proposed method provides a sound mathematical framework to incorporate different information types in anatomical structure registration such as thickness or volume of the structure at every point.

Conclusion

In this work we present a method for registering anatomical structures represented as point-sets endowed with a shape feature represented these structures as tensor fields. Our method computes a diffeomorphic transform between pairs of particle sets through a minimization process using an analytical gradient.

At first sight, this work is similar to registration work in diffusion tensor imaging (DTI) [START_REF] Daniel C Alexander | Spatial transformations of diffusion tensor magnetic resonance images[END_REF][START_REF] Zhang | Deformable registration of diffusion tensor MR images with explicit orientation optimization[END_REF][START_REF] Cao | Large deformation diffeomorphic metric mapping of vector fields[END_REF]. However, 3 main issues distinguish our application field and approach from DTI: 1) DTI tensors are Symmetric Positive Definite (SPD) while we are concerned with SPS ones. The SPD constraint requires DTI algorithms to keep tensors on the boundary of the SPS set. Hence the metric and transformation is different to ours; 2) The DTI deformation is different than ours: Most DTI deformable registration techniques transform each tensor through functions of rotations, i.e. without changing their eigenvalues. A main reason for this is the lack of knowledge on tissue microstructure changes (e.g. myelin density) driving eigenvalue changes across subjects; 3) We deal with a sparse set in space while DTI deals with an evenly sampled grid. The equivalence between these two problems is not simple. To convert the grid to a sparse sample we would need techniques like compressed sensing, the converse would require to pick a grid resolution and sample the continuous representation of the sparse set. All in all, both application domains deal with tensor fields, however they are far from equivalent in mathematical formalism and application domain.

Finally, we tested our algorithm in airway trees, synthetic and human. Our method outperformed the GMM and currents approaches. We also showed the performance of our method in the case of inter-subject registration. We validated these results using expert annotated airway datasets and comparing the overlap of the annotated regions after registration. Our results show an average overlap of higher than 70%.

A. Gradient of the Data Attachment Term

In this appendix we show the analytic gradient of the data attachment, eq. ( 7), term with respect to the transform parameter θ where eq. (11c) depends on the existence of the chain and product rules for derivation in F. When F is the space of square matrices with F, G F = trace(FG T ). Then, it can be shown that ∂ F F, G F = G and that if γ : F → R, then ∂ F γ(F) = γ (F), F F [START_REF] Peter | Linear Algebra and Its Applications[END_REF]. These two properties are the ones allowing use to derive eq. ( 5) obtaining the tensor field derivative eq. [START_REF] Arsigny | A Log-Euclidean Framework for Statistics on Diffeomorphisms[END_REF]

Figure 1 :

 1 Figure 1: Comparison of airway tree models. As shown in figs b and c, lower zoom in box, the tensor field represents changes in the bronchi thickness while the currents approach does not.

≥0

  F i , w : R 3 → Sym 3

Figure 2 :

 2 Figure 2: Registration of tensor fields under synthetic deformations using the algorithm in section 2.3. The metric • W effectively registers tensor fields under affine transforms.

Figure 3 :

 3 Figure 3: Synthetic airways. Left: ten realizations of random airways with deformation noise γ = 10. Center: we register the green airway tree onto the white one. The result is shown in red and the deformation field as arrows. The registration quality is such that the red is hard to tell from the white. Right: quantitative comparison of on the synthetic registration of 300 cases. It is noticeable that our method has a lower mean distance and volume difference with the original position than GMM and currents; while retaining a comparable performance in terms of angle

Figure 4 :

 4 Figure 4: Human data. Left: Airway tree with four generations colored: 1: Red; 2: Yellow; 3: Green; 4: Cyan; 5: Blue. Center: Airways of 10 different subjects affinely registered. Right: The same airways registered with our deformable tensor field approach

∂ θ d 2 = ∂ θ φ * θ w 2 -w 1 2 F∂ θ w 1 2 F + φ * θ w 2 2 F -φ * θ w 2 , w 1 F

 2222221 (w 1 , φ * θ w 2 ) = ∂ θ φ * θ w 2dxdue to the cosine rule in F:= dx (11b) deriving w.r.t. θ: = 2 φ * θ w 2 ,∂ θ φ * θ w 2 Fw 1 ,∂ θ φ * θ w 2 F dx (11c)using eq. (2):= 2 φ * θ w 2 , ∂ θ φ * θ w 2 W -2 w 1 , ∂ θ φ * θ w 2 W (11d)

Acknowledgements

This work was supported by NIH grants 1R01HL116931-02,1R01HL116473-02,K25HL104085-04,P41EB015898