Aurélie Voisin
email: aurelie.voisin@inria.fr

Vladimir A Krylov
email: vladimir.krylov@inria.fr

Gabriele Moser
email: gabriele.moser@unige.it

Sebastiano B Serpico
email: sebastiano.serpico@unige.it

Josiane Zerubia
email: josiane.zerubia@inria.fr

Classification of multi-sensor remote sensing images using an adaptive hierarchical Markovian model

Keywords: Supervised classification, multi-sensor data, hierarchical Markov random fields, copulas, discrete wavelet transform

INTRODUCTION

Remote sensing image classification is of capital interest when investigating land-use, land-cover or areas damaged by natural disasters. The aim of this paper is to develop a novel classification method that can be applied to coregistered multi-sensor remote sensing images, an interesting typology which importance is underlined by the high diversity of the acquisition systems that are currently operating. Here, we focus on the relevant case of high resolution synthetic aperture radar amplitude (SAR) [START_REF] Oliver | Understanding Synthetic Aperture Radar images[END_REF] and optical acquisitions.

A common way to classify multi-sensor data is to combine multiple classifiers, by using, for instance, boosting [START_REF] Schapire | A brief introduction to boosting[END_REF],

A. Voisin thanks the Direction Générale de l'Armement (DGA, France) and Institut National de Recherche en Informatique et Automatique (INRIA, France) for the partial financial support of her PhD. support vector machines [START_REF] Waske | Fusion of support vector machines for classification of multisensor data[END_REF], or Dempster-Shafer theory of evidence [START_REF] Al-Ani | A new technique for combining multiple classifiers using the Dempster-Shafer theory of evidence[END_REF]. Posterior probabilities may also be combined [START_REF] Benediktsson | Consensus theoretic classification methods[END_REF]. But, most of these methods tend to degrade the acquisitions by requiring a re-sampling of the data given the fact that multi-sensor data intrinsically implies a multi-resolution aspect of such acquisitions. To our best knowledge, so far, no multi-sensor classifiers have overcome this limitation.

The Bayesian classification technique proposed in this paper is based on an explicit hierarchical Markovian context [START_REF] Graffigne | Hierarchical Markov random field models applied to image analysis: a review[END_REF] and introduces a method to jointly model and classify the different data sources instead of fusing the final classifications or the input images. Moreover, at each level of the hierarchical graph, a flexible number of the input sources (images) is allowed, to adapt the graph to the available input dataset. The method consists of two steps. The first step is the statistical pixelwise modeling of the coregistered input images: for each class and each channel in the chosen stacked-vector input dataset, the marginal probability density functions (PDFs) are estimated by finite mixtures of automatically chosen parametric families [START_REF] Krylov | Supervised high resolution dual polarization SAR image classification by finite mixtures and copulas[END_REF], typically Gaussian mixtures for optical images or generalized Gamma mixtures for SAR imagery. Next, a d-variate copula [START_REF] Nelsen | An introduction to copulas[END_REF] (d being the number of input channels) is applied to estimate multivariate joint class-conditional statistics, merging the marginal PDF estimates of the input channels. The best copula is selected among a predefined dictionary of copulas according to a chi-square goodness-of-fit test [START_REF] Krylov | Supervised high resolution dual polarization SAR image classification by finite mixtures and copulas[END_REF]. On the second step, we plug the estimated joint PDFs into a contextual model that uses a multi-scale approach via a hierarchical Markovian model [START_REF] Graffigne | Hierarchical Markov random field models applied to image analysis: a review[END_REF] based on a quad-tree structure. The hierarchical structure is quite robust to noise and preserves spatial details of high resolution acquisitions. The quad-tree allows to integrate an exact estimator of the marginal posterior mode (MPM) [START_REF] Laferte | Discrete Markov modeling and inference on the quad-tree[END_REF] that aims to determine the unknown class labels. The prior probability is iteratively updated at each level of the tree, leading to an algorithm more robust with respect to speckle noise in SAR acquisitions [START_REF] Oliver | Understanding Synthetic Aperture Radar images[END_REF] when compared to a non-updated prior [START_REF] Laferte | Discrete Markov modeling and inference on the quad-tree[END_REF][START_REF] Voisin | Multichannel hierarchical image classification using multivariate copulas[END_REF]. Multi-scale features are either the multi-resolution images, or the features extracted by discrete wavelet transforms (DWT) [START_REF] Mallat | A wavelet tour of signal processing[END_REF][START_REF] Daubechies | Orthonormal bases of compactly supported wavelets[END_REF]. A similar model has already been introduced in [START_REF] Voisin | Multichannel hierarchical image classification using multivariate copulas[END_REF] in the case of mono-resolution multi-band images, and is extended and validated here on multi-sensor and multi-resolution acquisitions that are directly integrated in the adaptive explicit hierarchical graph.

The rest of the paper is organized as follows. In Sec. 2, we develop the MPM model based on the quad-tree that we use to classify our data. In Sec. 3, we introduce the statistical multivariate copula-based model that combines the marginal PDF models of the input images at each resolution. In Sec. 4, we present an experimental validation of the developed model.

THE HIERARCHICAL CLASSIFIER MODEL

We aim to estimate a set of hidden labels X given a set of observations Y . X and Y are random processes, and X is a Markov random field (MRF) on the considered hierarchical graph. To address our classification problem, we use an explicit hierarchical graph-based model [START_REF] Graffigne | Hierarchical Markov random field models applied to image analysis: a review[END_REF][START_REF] Laferte | Discrete Markov modeling and inference on the quad-tree[END_REF] that has a specific quad-tree structure. Such a structure ensures the scale causality which allows a non-iterative algorithm to be used, thus implying a computational time decrease as compared to an iterative computation. Among the various algorithms proposed to estimate the labels on hierarchical graphs [START_REF] Graffigne | Hierarchical Markov random field models applied to image analysis: a review[END_REF], we use an exact estimator of the marginal posterior mode (MPM). A classical MPM algorithm [START_REF] Laferte | Discrete Markov modeling and inference on the quad-tree[END_REF] is generally run in two passes, referred to as bottom-up ("forward") and top-down ("backward") passes, that respectively aim to estimate the partial posterior marginals p(x s |y d(s) ), s representing a site and d(s) its descendants, and then to estimate successively the labels at each level by maximizing the posterior marginals p(x s |y). To improve the robustness with respect to the speckle noise [START_REF] Voisin | Multichannel hierarchical image classification using multivariate copulas[END_REF], we propose to truncate the top-down pass: we use the highest tree-level label estimates to update the prior, and we run a novel MPM algorithm on a smaller quad-tree (see Fig. 1).

The maximization of the posterior marginal p(x s |y) is done by employing a modified Metropolis Dynamics algorithm which has good properties for both its relative low computation time and the good precision of its results [START_REF] Berthod | Bayesian image classification using Markov random fields[END_REF].

To estimate the posterior probability, we need information about the likelihood, the prior probability and the transition probability at each site s of the quad-tree. The likelihood term p(y s |x s ) is estimated for each stacked-vector input dataset at each tree level according to the model presented in Sec. 3. The transition probabilities between the scales, that determine the hierarchical MRF by representing the causality of the statistical interactions between the different levels of the tree, are expressed in the form introduced by Bouman et al. [START_REF] Bouman | A multiscale random field model for Bayesian image segmentation[END_REF]. Specifically, for all sites s ∈ S and all scales n ∈ [0

; R -1], p(x n s = ω m |x n s -= ω k ) = θ n , if ω m = ω k 1-θn M -1 , otherwise , (1)
with the parameter θ n > 1/M and where ω m and ω k represent the classes m and k, m, k ∈ [1; M ], respectively, and M represents the number of considered classes for the final classification. This model favors an identical parent-child labeling.

Finally, the prior distributions at the current coarsest tree level are given by the estimated classification maps by resorting to a single-scale Markovian model [START_REF] Krylov | Supervised high resolution dual polarization SAR image classification by finite mixtures and copulas[END_REF]. The prior information at lower levels is determined with the knowledge of the transition probabilities. For a more comprehensive study, see [START_REF] Voisin | Multichannel hierarchical image classification using multivariate copulas[END_REF].

JOINT PDF MODEL

At each tree level n and for each considered class m we perform the joint statistics modeling of the multi-sensor input images (e.g optical color bands and/or polarization SAR image channels), or their respective DWT features. The joint PDFs are independently modeled at each scale n, hence, in this section, we consider that n is fixed to simplify the notations. At each n, the following procedure is applied: the statistics of each input band are independently modeled, and are then used to generate a joint PDF by employing the statistical instrument of copulas [START_REF] Nelsen | An introduction to copulas[END_REF].

Marginal PDF estimation

For each input image, we want to estimate the distributions of each class ω m considered for the classification, m ∈ [1; M ], given a training set. The PDF f j m (z j ) of the j th input band, j ∈ [1; d], is then modeled via finite mixtures of K j independent greylevel distributions:

f j m (z j ) = p j m (z j |ω m ) = K j i=1 P j mi p j mi (z j |θ j mi ), (2)
where z j is a greylevel, z j ∈ [0; Z -1]. P j mi are the mixing proportions such that for a given m, K j i=1 P j mi = 1 with 0 ≤ P j mi ≤ 1. θ j mi is the set of parameters of the i th PDF mixture component of the m th class. The use of finite mixtures instead of single PDFs offers the possibility to consider heterogeneous PDFs, usually reflecting the contributions of the components present in each class, for instance, different kinds of crops for the vegetation class. Moreover, the use of finite mixtures can be seen as a generalization of the determination of a single PDF, and allows to estimate both the best finite mixture model and/or the best single PDF model.

When the j th input channel is an optical acquisition, we consider that the PDF f j m (z j ) related to each class can be modeled by a finite mixture of Gaussian distributions. This choice is consistent with the fact that a Gaussian distribution is a usually accepted model for a homogeneous area in an optical remote sensing image. For each mixture component, the mean and the variance parameters are estimated by a stochastic expectation maximization (SEM) [START_REF] Celeux | Stochastic versions of the EM algorithm: an experimental study in the mixture case[END_REF] algorithm. The SAR acquisitions are known to be affected by speckle noise [START_REF] Oliver | Understanding Synthetic Aperture Radar images[END_REF], thus we propose to model each class marginal conditional PDF by a mixture of generalized Gamma distributions. When the generalized Gamma distribution can not be fitted, the PDF p j mi (z j |θ j mi ) is automatically chosen between log-normal, Weibull and Nakagami distributions, that are all commonly employed to model SAR imagery [START_REF] Krylov | Supervised high resolution dual polarization SAR image classification by finite mixtures and copulas[END_REF]. In that case, the best-fitting mixture model for each considered class is obtained by combining a density parameter estimation via the method of log-cumulants [START_REF] Tison | A new statistical model for Markovian classification of urban areas in high-resolution SAR images[END_REF] and a SEM algorithm [START_REF] Krylov | Supervised high resolution dual polarization SAR image classification by finite mixtures and copulas[END_REF].

The estimation of K j is integrated into the SEM procedure [START_REF] Krylov | Supervised high resolution dual polarization SAR image classification by finite mixtures and copulas[END_REF].

Combination of marginal distributions via multivariate copulas

A multivariate copula is a d-dimensional joint distribution defined on [0, 1] d such that marginal distributions are uniform on [0, 1]. According to Sklar's theorem [START_REF] Nelsen | An introduction to copulas[END_REF], we can obtain a unique joint d-variate PDF for any continuous random variables z 1 , ..., z d in R:

h m (z 1 , ..., z d ) = d j=1 f j m (z j ) × c m (F 1 m (z 1 ), ..., F d m (z d )), (3)
where {f 1 m , ..., f d m } are the marginal PDFs estimated for the j th input image and for the m class, defined in Sec. 3.1, and {F 1 m , ..., F d m } their corresponding cumulative distribution functions. c m is the density of the copula C m . For each m, the computation of the joint PDF h m is boiled down to the determination of the copula family C m according to Eq. ( 3).

The bivariate copulas have been studied extensively in the literature, which is not necessarily the case for multivariate copulas. We only focus in this paper on Archimedean copula families, which represent a good compromise between analytical tractability and their modeling capacity [START_REF] Nelsen | An introduction to copulas[END_REF]. We consider a predefined dictionary of three single-parameter copulas: Clayton, Ali-Mikhail-Haq (AMH) and Gumbel copula families. For each class, we determine the best fitting copula as the one reporting the highest p-value in Pearson Chi-square (PCS) goodness-of-fit test, see [START_REF] Krylov | Supervised high resolution dual polarization SAR image classification by finite mixtures and copulas[END_REF]. The null hypothesis in PCS is that the sample frequencies C m (F 1 m (u 1 ), ..., F d m (u d )), where (u 1 , ..., u d ) represent the observed data, are consistent with the theoretical frequencies C m (v 1 , ..., v d ) for each considered copula.

To estimate the copula parameters, we use the relationship between copulas and Kendall's τ [START_REF] Nelsen | An introduction to copulas[END_REF], a concordancediscordance measure that can be empirically estimated by using the training sets.

EXPERIMENTAL RESULTS

We applied the developed classification approach to a multisensor dataset that consists of two images of the port area of Port-au-Prince (Haiti): an HH-polarized COSMO-SkyMed (CSK) SAR image (StripMap acquisition mode, 2.5 m pixel spacing, geocoded, single-look) of 320 × 400 pixels and a coregistered pan-sharpened (Pan.) GeoEye acquisition of 1280 × 1600 pixels, see Fig. 2(a),(b). Initially, the GeoEye image resolution is of 0.5 meters. To fit with the dyadic decomposition imposed by the quad-tree, we re-sampled the data to obtain the 0.625 m resolution (equal to 1/4 of the SAR image's resolution). To employ the multi-resolution information, and to avoid "empty" tree levels, we integrate additional information by employing a hierarchical decomposition of the original data. For our specific dataset, the optical image is integrated at the bottom of the quad-tree, and the SAR image, at level R = 2. We, thus, apply a 2-D DWT to decompose the finest resolution image along 2 levels or more, and integrate the obtained images in the tree. At the bottom and at the first stage of the graph, we only have one observation, and at higher tree levels, we have two observations, the SAR image and the corresponding decompositions of the optical one, that are combined by using the model presented in Sec.3. We can further increase the number of levels of the tree by introducing the SAR image decomposition.

Given the size of the input image, R = 2 appears to be a reasonable choice: a higher decomposition level might remove some relevant image details (e.g., roads). A wide choice of wavelet functions exists, and according to a preliminary comparative study [START_REF] Voisin | Multichannel hierarchical image classification using multivariate copulas[END_REF], we chose Daubechies-10 wavelets and Haar wavelets to decompose the SAR and the optical images respectively. By a trial-and-error procedure, the singlescale MRF parameter (see Sec. 2) was set to 4.8 and the transition probability θ n was set to 0.8 in [START_REF] Oliver | Understanding Synthetic Aperture Radar images[END_REF], which means that a site s at scale n has a probability of 80% to belong to the same class as its ascendant s -.

Manually annotated non-overlapping training and test sets were selected in homogeneous areas (no borders were taken into account). Five classes were chosen: water, urban areas, vegetation, sand and containers. We compare classification results obtained: 1) with the proposed approach on the multisensor SAR/optical dataset, 2) with the same approach applied only to the optical image and its decomposition, and 3) with a single-scale MRF-based approach [START_REF] Krylov | Supervised high resolution dual polarization SAR image classification by finite mixtures and copulas[END_REF] applied only to the optical image. The classification results are presented in Fig. 2(c),(d) and numerically assessed via accuracy estimates (Tab. 1). The classification map obtained when using a multi-sensor input shows satisfactory results because they lead to a detailed classification (Fig. 2(c)). Numerical results highlight the improvement related to the combination of the multi-sensor images as compared to single-sensor (only Pan.) classification performances, especially in the urban areas, in which the SAR acquisition seems helpful. The presented table (Tab. 1) indicates that the single-scale MRF-based method leads to a higher numerical accuracy. But we stress that the map obtained by using the MRF-based method is severely oversmoothed (Fig. 2(d)), and this affects only marginally the numerical accuracies due to the localization of the ground truth samples within homogeneous regions. Visually, the classification obtained with the hierarchical classifier is more detailed, which renders the proposed method preferable. Yet, it may be interesting to explore a prior despeckling of the SAR image [START_REF] Parrilli | A nonlocal SAR image denoising algorithm based on LLMMSE wavelet shrinkage[END_REF] to decrease the noise effects on the final results.

CONCLUSION

The method proposed in this paper allows to deal with multisensor acquisitions for land-cover mapping purposes. The main novelty lies in the direct integration of both multiresolution and multi-sensor data in the hierarchical graph model without resorting to any external fusion as pre-or post-processing. The proposed methodological approach combines a joint copula-based statistical model for optical and SAR acquisitions, with a hierarchical MRF, leading to a statistical supervised classification approach. We have proposed a dictionary-based multivariate copulas model that enables to fuse multi-sensor acquisitions, and developed a novel MPM-based hierarchical Markov random field model that iteratively updates the prior probabilities and, thus, leads to the improved robustness of the classifier. The hierarchical MRF considered here has three advantages: it is quite robust to speckle noise, it does not require a time-expensive iterative optimization algorithm and it is sufficiently flexible to allow the integration of both multi-sensor and multi-resolution acquisitions. As a direction of future development, we intend to perform experiments with wider datasets including the novel high-resolution Pléiades ( c CNES) optical imagery.

ACKNOWLEDGMENT

Fig. 1 .

 1 Fig. 1. Proposed MPM estimation on the quad-tree on which the set of sites is hierarchically partitioned as S = S 0 S 1 ... S R , where R corresponds to the coarsest resolution (the root), 0 corresponds to the reference level (finest resolution) and S n is the collection of sites at the n th scale level (n = 0, 1, ..., R). In this figure, R = 2.

Fig. 2 .

 2 Fig. 2. (a) Optical ( c GeoEye, 2010) and (b) SAR images (CSK, c ASI, 2010), and different classification method results obtained with: (c) the proposed method applied to the multi-sensor dataset, and (d) a single-scale MRF applied to the highest resolution image. Legend: water (blue), urban areas (red), vegetation (green), sand (yellow) and containers (purple).

Table 1 .

 1 Results for the Port-au-Prince dataset: by class and overall accuracy.

			Port-au-Prince, Haiti
		water urban	vegetation sand	containers overall
	Proposed classif.	100% 75.24% 87.16%	98.89% 49.31%	82.12%
	Proposed classif. (Pan. only) 100% 67.12% 86.89%	98.83% 41.90%	78.95%
	MRF-based classif.	100% 100% 81.42%	99.94% 59.62%	88.20%

The authors would like to thank the Italian Space Agency for providing the COSMO-SkyMed (CSK R ) images (COSMO-SkyMed Product -c ASI -Agenzia Spaziale Italiana. All Rights Reserved) in the framework of the project "Development and validation of multitemporal image analysis methodologies for multirisk monitoring of critical structures and infrastructures (2010-2012)". The authors would also like to thank GeoEye Inc. and Google crisis response for providing the GeoEye image available on the website http://www.google.com/relief/haitiearthquake/geoeye.html.


Viviana Siless
email: viviana.siless@inria.fr

Joan Glaunès

J E A N -F R A N Çois Mangin

Cyril Poupon

Denis Le Bihan

Bertrand Thirion

Pierre Fillard

Joint T1 and Brain Fiber Log-Demons Registration Using Currents to Model Geometry

Keywords: Registration, neural fibers, diffeomorphism, Demons Algorithm, intensity-base registration, tensor-base registration, log-domain

We present an extension of the diffeomorphic Geometric Demons algorithm which combines the iconic registration with geometric constraints. Our algorithm works in the log-domain space, so that one can efficiently compute the deformation field of the geometry. We represent the shape of objects of interest in the space of currents which is sensitive to both location and geometric structure of objects. Currents provides a distance between geometric structures that can be defined without specifying explicit point-to-point correspondences. We demonstrate this framework by registering simultaneously T1 images and 65 fiber bundles consistently extracted in 12 subjects and compare it against non-linear T1, tensor, and multi-modal T1+ Fractional Anisotropy (FA) registration algorithms. Results show the superiority of the Log-domain Geometric Demons over their purely iconic counterparts.

Introduction

Non-linear image registration is one of the most challenging tasks in medical image analysis. For inter-individual comparison, registration should align images as well as cortical and internal structures such as sulcal lines and fibers. Non-linear registration algorithms can be categorized into iconic, geometric and hybrid.

Iconic, or image-based registration [START_REF] Thirion | Image matching as a diffusion process: an analogy with Maxwell's demons[END_REF][START_REF] Beg | Computing large deformation metric mappings via geodesic flows of diffeomorphisms[END_REF][START_REF] Rueckert | Diffeomorphic Registration Using B-Splines[END_REF] finds a voxel-wise mapping between a source and a target image. Schematically, iconic registration is mainly driven by the contours of the image, and without prior knowledge, it is difficult to coregister regions with little contrast. For instance, brain white matter appears uniformly white in T 1 images, giving no relevant information to the iconic registration, while it is composed of neural fibers connecting cortical areas. Diffusion Tensor Imaging (DTI) can be used to reveal the microscopic structure of the white matter. Tensor-based registration was recently proposed to improve white matter alignment [START_REF] Zhang | Deformable registration of diffusion tensor mr images with explicit orientation optimization[END_REF][START_REF] Yeo | Dt-refind: Diffusion tensor registration with exact finite-strain differential[END_REF]. However, misregistration may persist in regions where the tensor field appears uniform, as shown in [START_REF] Durrleman | Registration, atlas estimation and variability analysis of white matter fiber bundles modeled as currents[END_REF].

Geometric registration specifically targets the alignment of Structures of Interest (SOI), such as in [START_REF] Yeo | Spherical demons: Fast diffeomorphic landmark-free surface registration[END_REF] for cortical surfaces, or [START_REF] Durrleman | Registration, atlas estimation and variability analysis of white matter fiber bundles modeled as currents[END_REF] for fiber bundles. While those clearly improve SOI registration, they are in general not suitable for aligning other structures than those used specifically during registration.

Hybrid techniques propose to jointly consider SOI and images during registration. For instance, [START_REF] Auzias | Diffeomorphic brain registration under exhaustive sulcal constraints[END_REF][START_REF] Ha | Image Registration Driven by Combined Probabilistic and Geometric Descriptors[END_REF] used the mathematical framework of measures and currents to simultaneously register images and geometric descriptors, while [START_REF] Sotiras | Simultaneous Geometric -Iconic Registration[END_REF] proposed a Markovian solution to the same problem.

We present an hybrid registration algorithm based on the efficient framework of Demons, where we combine iconic and geometric registration. The geometry is represented in the space of currents which provides a metric sensitive to shape.

The rest of the paper is organized as follows. First, we propose a mathematically sound extension of the Geometric Demons(GD), the Log-domain Geometric Demons (LGD), that relies in the log-domain framework for computation and handles geometric constraints in the space of currents. Then, we evaluate the LGD with bundles constraints on a dataset of 12 subjects and compare them with a scalar [START_REF] Vercauteren | Symmetric Log-Domain Diffeomorphic Registration: A Demons-Based Approach[END_REF], a tensor [START_REF] Yeo | Dt-refind: Diffusion tensor registration with exact finite-strain differential[END_REF], and Ants [START_REF] Avants | A reproducible evaluation of ants similarity metric performance in brain image registration[END_REF] a multi-modal registration.

The Log-Domain Geometric Demons

The Diffeomorphic Demons

In image registration we search for a displacement field s between a fixed F and moving M image, that maps as accurately as possible corresponding structures in both images. Ideally the displacement field s minimizes the distance between the fixed and the moving image, while holding some properties such as being diffeomorphic. In the Demons framework [START_REF] Thirion | Image matching as a diffusion process: an analogy with Maxwell's demons[END_REF] a correspondence field c was introduced to make the minimization of the functional energy tractable:

E(c, s)= 1 σ 2 i Sim(F, M • c)+ 1 σ 2 x dist(s, c) 2 + 1 σ 2 T Reg(s)
where Sim is a similarity measure between images defined by the sum of square differences (SSD) and Reg a regularization term chosen to be an harmonic energy. The amount of regularization is controlled with σ T while σ i accounts for the image noise. The term dist(s, c) 2 imposes the displacement field s to be close to the correspondence field c. σ x weights the spatial uncertainty on the deformation. The energy minimization is performed by alternating minimization w.r.t. c and s.I n [START_REF] Vercauteren | Symmetric Log-Domain Diffeomorphic Registration: A Demons-Based Approach[END_REF], small deformations parametrized by a dense displacement field u are used: c ← s • exp(u), exp() being the exponential map in the Lie group sense, which ensures that the result is diffeomorphic.

Geometric Demons

To incorporate geometric constraints in the Demons framework, in [START_REF] Siless | Joint T1 and Brain Fiber Diffeomorphic Registration Using the Demons[END_REF] the definition of c was extended to carry information from both image and geometry.

Let us denote by G F (resp. G M ) the fixed (resp. moving) geometric descriptors. Therefore a new energy was defined:

E(c, s)= 1 σ 2 i Sim I (F, M • c)+Sim G (c⋆G F , G M ) + 1 σ 2 x dist(s, c) 2 + 1 σ T Reg(s), (1)
where Sim I is the image similarity measure, Sim G the geometry similarity measure, and c⋆G F denotes the action of c on the geometry. Following [START_REF] Vercauteren | Symmetric Log-Domain Diffeomorphic Registration: A Demons-Based Approach[END_REF], c was parametrized by an update field, as the additive combination of an image update field u I and a geometric update field u G .N o nintersecting domains were defined, as ideally, one should use u G only where geometric information is relevant and use u I elsewhere. In the case of fibers, the geometric domain will remain within the white matter. Thus, let Ω G be the definition domain of u G (where geometry is defined), and the definition domain of u I be Ω I = Ω -Ω G . Then we can define c =e x p ( u I + u G ). The following relationships hold:

c⋆G F =exp(u G ) ⋆ G F and M • c = M • exp(u I ).
Geometric Demons incorporates the following energy to calculate u G :

E G (s, u G )= 1 σ 2 i Sim G (s • exp(u G ) ⋆ G F , G M )+ 1 σ 2 x ΩG u G 2 , (2)
Being s the deformation field from F to M . Thus, the inverse of the s gives the geometric deformation. In section 2.3 we show an efficient approximation for obtaining the geometric deformation in the log-domain space.

Log-Domain Geometric Demons

The log-domain demons avoids the inversion of the deformation field by redefining s with the exponential map:

s =e x p ( v). Then s • exp(u)=exp( v) • exp(u)
and the Baker-Campbell-Hausdorff(BCH) formula yields log(exp(v)

• exp(u)) ≈ v + u +1/2[v, u]+1/12[v, [v, u]] + ... where [v, u]i st h eL i eb r a c k e t .
Then, the LogGeometric Demons algorithm is defined as follows:

1. Choose a starting spatial transformation s =exp(v) 2. Given s, u I , compute the update field u I as in [START_REF] Vercauteren | Symmetric Log-Domain Diffeomorphic Registration: A Demons-Based Approach[END_REF] 3. Given s, u G , compute the update field u G by minimizing Eq. ( 2)

4. Let u ← s • exp(u I + u G ) 5. v ← log(exp(v)
• exp(u)) using BCH approximation and exp(u) is efficiently computed with a few compositions, look [START_REF] Vercauteren | Symmetric Log-Domain Diffeomorphic Registration: A Demons-Based Approach[END_REF] for further details.

6. let v ← K dif f ⋆v where K dif f is a Gaussian convolution kernel 7. s =exp(v)a n ds -1 =exp(-v) 8. Go to 2.

until convergence

With this new definition we can efficiently compute the geometric deformation.

Calculation of u G in the Space of Currents. In GD the closest point distance was proposed as it does not need explicit correspondences between points and it can be a good approach for comparing single fiber bundles representatives as in [START_REF] Siless | Joint T1 and Brain Fiber Diffeomorphic Registration Using the Demons[END_REF]. In contrast, by representing geometry in the space of currents, we have a pose and shape-sensitive measure which permits to define a distance between bundles containing different number of fibers.

Let G be a set of continuous curves. We define the sequence of discretized points in G as G =(x 1 , ..., x N ), N being the number of points. We can associate to this sequence a specific measure given by the vector valued Diracs:

µ G = n-1 i=1 τ G,i δ cG ,i where c G,i =( x i + x i+1 )/2(center point) and τ G,i = x i+1 -x i (tangent vector) if x i and x i+1 belong to the same curve, otherwise τ G,i =0.
Following [START_REF] Glaunès | Large deformation diffeomorphic metric curve mapping[END_REF], let W be a reproducible kernel Hilbert space (r.k.h.s) of vector fields with kernel K W β isotropic and Gaussian of size β: vector fields in W are convolutions between any square-integrable vector fields and the convolution square root of the kernel. Then, the vector space of currents is a dense span of the set of all the vector valued Diracs currents τδ c for any τ, c ∈ R 3 .AD i r a c current may be seen as an oriented segment entirely concentrated at point c.The scalar product between two sums of vector valued Diracs expresses conveniently in terms of the kernel K W β :

µ, µ ′ = N i=1 τ i δ ci , M j=1 τ ′ j δ c ′ j = N i=1 M j=1 K W β (c i ,c ′ j )τ i .τ ′ j (3)
Having a fixed and a moving geometric descriptor G F =(x 1 , ..., x N )a n dG M = (y 1 , ..., y M ) N, M being the number of points, the distance is defined as follows:

d 2 (G F , G M )=||G F -G M || 2 W * = ||G F || 2 W * + ||G M || 2 W * -2 G F , G M W * (4)
The distance measures geometrical differences both in pose and shape. With β we define the kernel size, and points at distances much larger than β have a large distance disregarding the shape. Also, when distances are much smaller than β, they are taken as noise, and thanks to the smoothing effect of the kernel they are not taken into account. So the distance captures first misalignment and then shape dissimilarities until a noise level quantified by β is reached. Let us define the action of the correspondence field c on G as:

c⋆G = {s • exp(u G )(x i )} i∈[1,N ] ≈{s(x i )+u G (x i )} i∈[1,N ] .
Since we are dealing with discrete points, we choose to parametrize the dense update field u G by a finite set of vectors u G,i using radial basis function interpolation:

u G (x)= N i=1 h( x - x i )λ i .W h enh(x)=e -(r) 2 , λ i are calculated such that u G (x i )=u G,i ∀i.L e tu s define the matrix A such that [A] i,j = h( x i -x j )([A] i,j denotes the (i, j)entry of A), H(x) the vector such that [H(x)] i = h( x -x i )andU =[u G,1 , ..., u G,N ]. We can write: u G (x)=H(x)A -1 U .S o l v i n g∇E G (s, u G ) = 0 w.r.t. u G
narrows down to optimization for the u G,i , ∀i. After differentiation, we obtain:

u G,i = ∇ uG,i ||c⋆G F -G M || 2 W * 1+ σ 2 i σ 2 x [H(s(x i ))A -1 ] i (5)
Defining Ω G for Bundles. Since we want fibers to influence the deformation near the definition domain, we define the domain as the union of γ-radius balls B centeredateachcoordinatex i . We control the influence by varying γ and thus, dilating the domain. We define a binary map Ω γ G = N i=1 B(x i ,γ) . The domain of the image correspondence field is the complementary of Ω γ G :

Ω γ I = Ω\Ω γ G .
3J o i n t T 1 MRI and Brain Bundle Registration

Data Description

We used the NMR database of 12 healthy volunteers scanned with T 1 (256×256× 124,.9375 × .9375 × 1.2mm) and DW-MRI (128 × 128 × 60, 1.875 × 1.875 × 2mm) [START_REF] Poupon | A database dedicated to anatomo-functional study of human brain connectivity[END_REF]. 200 encoding gradients were used for the diffusion sequence. Using [START_REF] Guevara | Robust clustering of massive tractography datasets[END_REF], we obtained corresponding fiber bundles between several subjects. 100 bundles were consistently identified in all subjects. The longest 65 bundles distributed in both hemispheres were retained for the experiments. For each subject we obtained the linear transformation from B 0 to T 1 to align bundles with T 1 images.

Experiments

Two experiments were conducted. First, we exhaustively analyze the parameter γ in ω γ G defined in Sec. 2.3 to understand its effect on registration accuracy. Second, we compared the performance between Symmetric Log Domain Demons (SLDD), the Symmetric Tensor Demons (STD) and Ants. The inverted deformation field was applied to the fibers to display the registration. Each algorithm was tested on the 3-steps multi-scale approach with 15, 10 and 5 iterations at each scale (from small to large). We set the currents kernel size β by using a robust estimator of the maximum distance between bundles thresholded at 20mm.

Influence of γ.

In the first experiment, the 11 subjects were registered onto one, arbitrary chosen as the target subject. We varied γ f r o m0t o4 . 5 ,w h e r eγ is scaled by the smallest voxel size. We divided our bundles in 5 sets (13 bundles each, with bundles of ∼ 3 fibers, each of 21 points), and used jointly 4 to train, and the left one to test. The following results show the average of the 5 possible permutations of choosing the test set. We show results over training set (a), test set (b), and the image (c) for the increasing values of γ in Fig. 1. As expected, fiber matching improves as γ increases(a) at the expense of image alignment(c). Indeed, when fibers have a large influence on their neighborhood, image-driven forces are discarded, leading to poor image registration. Also, comparing (a) and (b) we note that γ =4 .5 is overfitting the fibers, misleading the overall registration. γ =3 .0 largely improves fiber alignment, while keeping a good match between images. In the sequel, a γ =3.0 will be used. In some cases image matching is slightly improved when using γ =1 .5 pointing out that geometry may indeed help image registration to avoid local minima. Comparison with SLDD, STD and ANTS. For this experiment we register all subjects to each subject for each permutation of the training set. For SLDD, we registered 11 T 1 images onto the 12 permuted target and applied the inverted deformation field to the bundles. For STD, we extracted tensors using [START_REF] Fillard | Clinical DT-MRI estimation, smoothing and fiber tracking with log-Euclidean metrics[END_REF] and registered each of them onto the 12 permuted target tensor image. Then, inverted deformation fields were applied to each subject's fibers in the DWI space. Finally, the linear transformation calculated between the target B 0 and T 1 images was used to carry fibers to the T 1 space. For Ants, we extracted the FA from the tensors obtained using [START_REF] Fillard | Clinical DT-MRI estimation, smoothing and fiber tracking with log-Euclidean metrics[END_REF], and aligned them to their T 1 image. We use the cross correlation setting with weighting equally image and FA. Then affine and nonrigid resulting transformation were applied to images, while the inverse of the non-rigid and the affine were applied to the fibers. We show the average metric over training sets (i), test sets (j) and image (k) of registering all subjects to each one with the methods mentioned above in Fig. 2.

As expected, LGD further improved fiber registration in (i) compared to the other algorithms. However, the training set contains the fibers used during registration; we explicitly optimize a metric evaluated on those fibers. Analyzing the results over the test set in (j) we see a similar performance compared to Ants, and STD, which is remarkable as STD is using information from tensors over the whole dense grid, and Ants does a cross correlation between the whole FA grid and T 1 . By contrast, in LGD the deformation field was obtained using only sparse information coming from selected fibers, which are not defined in the regions tested in (j). Therefore having a similar performance is very promising. We can also see in (k) the image registration for STD was extremely poor. We time all algorithms with an Intel Xeon 8proc. 2.53GHz, 11.8Gb and obtained: SLDD=19.61min, STD=10.75min, Ants=25.63min, and LGD=12.51min.

Discussion

We compared our algorithm against a scalar (SLDD), a tensor (STD), and a multi-modal (Ants) registration. Results show that bundle alignment was highly improved comparing to other algorithms. We get accurate results even for testing set fibers where no information was used from the support regions of the those fibers, and STD and Ants were using information from the whole grid. This shows that a small set of fibers might be sufficient for a proper registration of the white matter across subjects. Moreover, while fiber alignment is improved, the efficiency of the image alignment is maintained. When evaluating the algorithm for the different γ values, we could see that γ =1 .5 better registers missing structures than γ = 0 and than SLDD. However, we believe there is a trade-off to make between image and fiber alignment, and γ = 3 notably improved fiber alignment while still obtaining good image registration results. By using labeled bundles instead of purely tensor information, we add relevant features that were previously extracted as prior such as region connection. Nevertheless, the efficacy of trusting bundles is open to discussion, and bundles classification is an active topic in research, so we believe this information should not be discarded.

5C o n c l u s i o n

We extended the GD algorithm to the log-domain space, and combined it with currents to compare our geometric structures and respect their shape. Our algorithm perfomed similarly as other competitive algorithms on image and showed a large improvement in fiber alignment. A unique mapping for images and geometric structures is obtained, giving a consistent framework for analyzing and comparing results between voxel-based morphometry and shape of SOI.

Having disjoint domains for geometry and iconic may potentially limit the incorporation of new geometric structures such as sulcal lines or features. This choice was done to guarantee a closed-form solution when zeroing out the gradient of the criterion. Moreover, both update fields are eventually combined and smoothed out, lowering the effect of the disjoint domains. Future work will further address the remaining mathematical issues to allow overlapping domains.

Based on the efficiency of the algorithm, it would be interesting to combine it with the clustering of fibers: increased fiber registration can help clustering algorithms, which can in turn guide the registration.

Fig. 1 .

 1 Fig. 1. Influence of γ. Similarity measure average of the 5 training set permutations with varying γ. Different color curves represent the 11 subjects.

Fig. 2 .

 2 Fig. 2. Comparison of SLDD, STD, Ants and LGD Top: Target fibers overlapped with registered fibers from an arbitrary chosen registration. 29 fibers were arbitrary preselected for clarity. Corresponding fibers in subjects share colors. Bottom: Average of the metrics obtained from registering 11 subjects to the target subject.

Acknowledgements. This work was supported by the INRIA CapNeonates Equipe Associée.