C Chnafa

S Mendez

R Moreno

F Nicoud

Using image-based CFD to investigate the intracardiac turbulence

A numerical framework designed to compute the blood flow in patientspecific human hearts is presented. The geometry of the heart cavities and associated wall motion are extracted from 4D medical images while the valves of the heart are accounted for thanks to low order geometrical models. The resulting blood flow equations are solved using a fourth-order low-dissipative finite-volume scheme and a mixed Aribtrary Lagrangian-Eulerian / Immersed Boundary framework. On top of retrieving the main fluid flow phenomena commonly observed in the left heart, the methodology allows studying the heart flow dynamics, including the turbulence characteristics and cycle-to-cycle variations.

Introduction

Heart pathologies are closely related to intracardiac hemodynamics. Recent technological innovations in imaging techniques have provided valuable opportunities for direct non-invasive in vivo assessment of hemodynamics. Blood flow velocities can be measured in vivo using phase-contrast magnetic resonance imaging (PC-MRI) or by echocardiography techniques.

PC-MRI studies have contributed to the understanding of the main hemodynamic features [START_REF] Fyrenius | Three dimensional flow in the human left atrium[END_REF][START_REF] Kilner | Asymmetric redirection of flow through the heart[END_REF][START_REF] Markl | Comprehensive 4D velocity mapping of the heart and great vessels by cardiovascular magnetic resonance[END_REF]. Although very comprehensive, the PC-MRI velocity mapping C. Chnafa I3M UMR 5149, University Montpellier II, FRANCE, e-mail: christophe.chnafa@univ-montp2.fr S. Mendez I3M, UMR 5149, University Montpellier II, FRANCE e-mail: simon.mendez@univ-montp2.fr R. Moreno Radiology Department, CHU of Rangueil, FRANCE e-mail: moreno.r@chu-toulouse.fr F. Nicoud I3M, UMR 5149, University Montpellier II, FRANCE e-mail: franck.nicoud@univ-montp2.fr 1 is not real-time. Hence, beat-to-beat variations in the flow cannot be recorded (the k-space is filled over many cardiac cycles). Moreover, PC-MRI suffers from a relatively low spatio-temporal resolution, precluding the observation of small scales and fast time-varying flow features [START_REF] Hollnagel | Comparative velocity investigations in cerebral arteries and aneurysms: 3D phase-contrast MR angiography, laser doppler velocimetry and computational fluid dynamics[END_REF].

Echocardiography techniques [START_REF] Garcia | Two-dimensional intraventricular flow mapping by digital processing conventional color-doppler echocardiography images[END_REF][START_REF] Kim | Development and validation of echo PIV[END_REF], with higher spatio-temporal resolution make an alternative to PC-MRI. However, they only give access to velocity components directed towards or away from the ultrasonic beam, while one would want to measure the full 3D flow vectors. Nevertheless, investigations have been conducted on normal and abnormal hearts and interesting potential indicators of cardiac health hemodynamics arose [START_REF] Cimino | In vivo analysis of intraventricular fluid dynamics in healthy hearts[END_REF][START_REF] Gharib | Optimal vortex formation as an index of cardiac health[END_REF][START_REF] Hong | Characterization and quantification of vortex flow in the human left ventricle by contrast echocardiography using vector particle image velocimetry[END_REF].

With the development of these cardiac imaging techniques, patient-specific geometries have been progressively used in computational fluid dynamics (CFD) [START_REF] Doenst | Fluid-dynamic modeling of the human left ventricle: Methodology and application to surgical ventricular reconstruction[END_REF][START_REF] Long | Subject-specific computational simulation of left ventricular flow based on magnetic resonance imaging[END_REF][START_REF] Merrifield | Combined CFD/MRI analysis of left ventricular flow[END_REF][START_REF] Mihalef | Patient-specific modelling of whole heart anatomy, dynamics and haemodynamics from fourdimensional cardiac CT images[END_REF][START_REF] Saber | Progress towards patient-specific computational flow modeling of the left heart via combination of magnetic resonance imaging with computational fluid dynamics[END_REF][START_REF] Schenkel | MRI-based CFD analysis of flow in a human left ventricle: Methodology and application to a healthy heart[END_REF]. Realistic heart wall movements on the basis of cine MRI or Computed Tomography (CT) scan data can be used: heart movement is prescribed from the patient-specific medical images, which can be acquired using standard clinical imaging procedures, instead of being computed. This strategy allows computation of the patient-specific hemodynamics and provides detailed insights into the cardiac flow field, providing potentially valuable clinical information. If the feasibility of this kind of approach has been shown, the results usually suffer from limited spatial resolution, partial geometry (only the left ventricle (LV) is considered in most cases) or numerical limitations (dissipative schemes). Besides, except in a few experimental works [START_REF] Falahatpisheh | High-speed particle image velocimetry to assess cardiac fluid dynamics in vitro: From performance to validation[END_REF][START_REF] Querzoli | Effect of the prosthetic mitral valve on vortex dynamics and turbulence of the left ventricular flow[END_REF], cycle-to-cycle variations in the heart flow is an issue that is rarely dealt with.

In the present work, an image-based CFD method developed to compute flows in aortas [START_REF] Midulla | Haemodynamic imaging of thoracic stent-grafts by computational fluid dynamics (CFD): presentation of a patient-specific method combining magnetic resonance imaging and numerical simulations[END_REF] is extended to compute intracardiac flows. Medical images are used to generate a moving patient-specific domain, in which the blood flow equations are solved. Heart geometry movements are generated from a 4D sequence (MRI or CT scan images) treated by an appropriate image registration algorithm [35]. A specific attention is paid to the generation of a high-quality mesh which deforms consistently with the heart motion. This allows solving the flow equations with an essentially non-dissipative scheme compatible with the description of unsteady turbulent flows by Large Eddy Simulation. In order to demonstrate the ability of the method to compute heart flows, application to a complete human left heart described by ECG-gated 3D CT scan images is presented. The resulting flow is described, emphasizing the flow characteristics usually reported in the literature. Characteristics of the cycle-to-cycle variations are also reported.

The numerical method is detailed in § 2. The characteristics of the applied case is presented in § 3 and the flow field obtained is described in § 4. Concluding remarks are given in § 5.

Methodology

The present computational method approaches couples an Arbitrary-Lagrangian Eulerian framework with an immersed boundary method in order to represent the blood flow within the moving endocardium while accounting for the heart valves motion. This section describes the flow solver and the treatment of the medical images needed to perform the computations.

Fluid problem

Blood is an incompressible, non-Newtonian fluid [START_REF] Fung | Biomechanics -Circulation[END_REF]. However, in large vessels, non-Newtonian effects are usually neglected in numerical simulations [START_REF] Mihalef | Patient-specific modelling of whole heart anatomy, dynamics and haemodynamics from fourdimensional cardiac CT images[END_REF][START_REF] Taylor | Patient-specific modeling of cardiovascular mechanics[END_REF] and constant kinematic viscosity ν is assumed in this paper. Note however that the present numerical method could be applied to non-Newtonian fluids.

The flow is thus governed by the incompressible Navier-Stokes equations (NSE), over a moving domain Ω f (t) ⊂ R 3 of boundary ∂ Ω f (t). The Arbitrary Lagrangian Eulerian (ALE) framework [START_REF] Fernandez | Cardiovascular Mathematics. Modeling and simulation of the circulatory system[END_REF] is used in order to account for the computational domain deformation over time. Introducing the pointwise computational domain velocity u g , the NSE read:

du dt + ((u -u g ) • ∇)u = - 1 ρ ∇p + ν∇ 2 u + f ∇ • u = 0,    on Ω f (t) (1)
where u is the fluid velocity, p is the pressure, ρ the density and f a force per mass unit. Note that the time derivative denotes an ALE time derivative [START_REF] Fernandez | Cardiovascular Mathematics. Modeling and simulation of the circulatory system[END_REF]. At the boundary of the computational domain, no-slip conditions are applied on walls (u = u g ) and Dirichlet conditions are applied over the inlet boundaries. These equations are implemented in the flow solver YALES2BIO (www.math.univ-montp2.fr/∼yales2bio) [START_REF] Mendez | An unstructured solver for simulations of deformable particles in flows at arbitrary reynolds numbers[END_REF], developed from the massively parallel finite-volume flow solver YALES2 [START_REF] Moureau | Design of a massively parallel CFD code for complex geometries[END_REF]. At each time step, the grid velocity u g at each node of the computational domain is calculated (see next section). A projection method is used to solve the NSE: the momentum equation is first advanced using the fourthorder Runge-Kutta time-advancement scheme. Fluxes are discretized with a 4thorder central scheme. Any forcing term coming from immersed boundaries is accounted for in this prediction step. At the end of the prediction step, the grid reaches the final position of the time step. Hence, the projection step to calculate pressure is performed over a fixed grid. A Deflated Preconditioned Conjugate Gradient algorithm is used to solve the Poisson equation [START_REF] Malandin | Optimization of the deflated conjugate gradient algorithm for the solving of elliptic equations on massively parallel machines[END_REF] involved in the projection step. Note that due to the transitional nature of the flow, Reynolds Averaged Navier-Stokes (RANS) approaches for modelling turbulence (e.g. k-ε, k-ω models, etc) are not appropriate since they rely on the assumption that turbulence is fully developed and ergodic. Even if adaptations have been proposed in order to handle transition, they essentially require the user to prescribe the transition location in advance. Instead, in the present study, the Large Eddy Simulation approach is followed, taking advantage of the low-dissipative scheme of integration used in YALES2BIO. In this view, only the smallest scales are modelled (scales smaller than the mesh size) while the evolution of the large scales is computed by solving a filtered version of NSE [START_REF] Moin | Advances in large eddy simulation methodology for complex flows[END_REF][START_REF] Pope | Turbulent flows[END_REF][START_REF] Sagaut | Large eddy simulation for incompressible flows[END_REF]. In the latter, a subgrid-scale model must be used in order to account for the effect of the unresolved scales on the dynamics of the resolved ones. This is usually done by an eddy-viscosity-based model [START_REF] Smagorinsky | General circulation experiments with the primitive equations: 1. the basic experiment[END_REF]. In the present study, where the flow is strongly confined and piloted by the wall motions, an advanced subgrid scale model able to represent the proper turbulence damping near solid walls was used [START_REF] Nicoud | Using singular values to build a subgrid-scale model for large eddy simulations[END_REF] as well as a well established formulation based on the dynamic procedure of Germano [START_REF] Germano | A dynamic subgrid-scale eddy viscosity model[END_REF].

Computational Domain

Extraction of the heart deformation

Determining the movement of the computational domain, where the NSE are solved is all but an easy task. One solution is to compute the deformations of the boundary thanks to a full electrical-fluid-structure interaction solver. Unfortunately, there are many uncertainties regarding the constitutive laws (both mechanical and electrical) of the heart muscle as well as the external constraints the heart is submitted to. Moreover, accurate heart models are still under development today [START_REF] Vadakkumpadan | Image-based estimation of ventricular fiber orientations for personalized modeling of cardiac electrophysiology[END_REF] as well as robust numerical method to solve them. A way of by-passing this issue is to extract the computational heart grid from 4D (3D + time) patient medical images. The following question must then be addressed: given several 3D images of a heart taken at different times in the heart cycle, how to extract the heart deformations from these images and how to deform a patient-specific grid accordingly? The first part of the question is actually a classical "image registration" problem.

Nowadays, there is a growing interest in the development of cardiac image registration methods [START_REF] Makela | A review of cardiac image registration methods[END_REF]. Given two cardiac images, a template and a reference one, a transformation is determined to map the template image to the reference image. The deformation field can notably provide clinical information on the myocardial contractile function. Here, the same process is used to compute the heart deformations, but instead of focusing on functional data of the heart muscle, the computed deformation is used to extract patient-specific boundary conditions for the blood flow computation.

Among many registration algorithms [START_REF] Maintz | A survey of medical image registration[END_REF][START_REF] Modersitzki | Numerical Methods for Image Registration[END_REF], voxel similarity measure techniques were preferred. This class of method operates directly on the image grey values, providing a flexible algorithm suitable for the complex heart movements.

In the present case, N 3D images I i of a patient heart are taken at different times t i , 0 ≤ t 0 ,t 1 , ...,t N-1 < T during the heart cycle of period T . One of the N images is selected as a template image. Note that this choice is somewhat arbitrary and that, without loss of generality, one can always tune the time origin so that the template corresponds to the I 0 image at time t 0 . From this image, N -1 transformations ψ i are computed such that the transformed template image becomes similar to images I i (i between 1 and N -1) successively considered as the reference image: transformation ψ i is search so that I 0 (ψ i (x)) = I i (x) for each voxel. Here, I i (x) stands for the voxel grey-level value at position x in the i-th 3D discrete image.

Each ψ i is calculated by minimizing the distance between I 0 (ψ i (x)) and I i (x), an appropriate distance measure being based on the so-called sum of squared differences. The transformation ψ i is computed through an optimisation problem reading: given two images I 0 and I i , find a mapping ψ i such that the squared intensity dif- ferences between I 0 (ψ i (x)) and I i (x) is minimized for each voxel. The number of parameters describing the transformation is too high so that the solution ψ i of this problem is not unique. Additional constraints are needed to compute the mapping ψ i .

In the present work, a constraint on ψ i is applied thanks to prior knowledge of the deformation sought [START_REF] Ashburner | Incorporating prior knowledge into image registration[END_REF]. The idea is to penalize unlikely deformations by imposing the heart deformations to be smooth. Bayesian statistics are used to obtain an a posteriori computation of the deformation field. The prior deformation probability is incorporated through the Bayes' theorem: p(Y|I) ∝ p(I|Y)p(Y), where p(I|Y) is the likelihood of observing the images data I (template I 0 and reference I i images) given the deformation parameters Y. p(Y) is the prior knowledge of the deformation translated in the a priori probability of seeing the parameters Y and p(Y|I) is the a posteriori probability of getting Y knowing the two images data I. Using this Bayesian framework, the goal is to maximise the probability p(Y|I). Knowing that a probability is related to its Gibbs form by p(Y) ∝ e -H(Y) , the problem can be seen as a minimisation of the Gibbs potential:

H(Y|I) = H(I|Y) + H(Y) + c, (2)
where c is a constant. The likelihood potential H(I|Y) of observing the images data given the deformation parameters Y is directly linked to the squared intensity difference between I 0 (ψ i (x)) and I i (x) for each voxel:

H(I|Y) = 1 2 Ω (I 0 (ψ i (x)) -I i (x)) 2 dΩ , (3)
where the integral is taken over the image volume Ω . The second term H(Y) in the right hand side of equation ( 2) is the wanted constraint, applied thanks to prior knowledge. The prior deformation knowledge put in this potential is expressed as a geometrical constraint on the mapping through the two studied images. A suitable prior probability is linked to the deformation of each voxel of the template image and to the reverse deformation [START_REF] Ashburner | High-dimensional image registration using symmetric priors[END_REF]. This potential acting as a regularization term allows a penalization based on the Jacobian of the locals deformations J. A detailed description of this term can be found in [START_REF] Ashburner | High-dimensional image registration using symmetric priors[END_REF].

This regularization term is weighted by a parameter λ linked to the belief in the amount of deformation of the heart. A high value of the λ parameter results in a high penalty on the voxels deformation, hence only small and smooth deformations are allowed. Given the high deformation of the left heart, high λ values result in partial deformation and a high residual squared difference between the template image and the reference one. In the case of small values of λ , a lower residual squared difference will be reached, but the resulting deformation can be non physical due to excessive warping. Actually, the variability of heart movements is highly different from one region to the other. The aorta needs high values of λ , whereas small values are needed for the atrium or the ventricle.

One method of increasing the likelihood of achieving a good solution without introducing a spatial variability on λ (and consequently avoiding more operatordependant work as well) is to apply successive filters to the images using a Gaussian smoothing convolution kernel G σ of width σ . High-frequency information of the image is removed thanks to this filter then gradually re-introduced as the kernel width σ becomes smaller and smaller. This iterative process is defined as the outer iterations: for each kernel width, distance between the images are minimized. This "coarse-to-fine" strategy has the effect of making the registration algorithm estimate the most global deformations during the first outer iteration, leaving out fine-scale structures. The optimum transformation for this kernel width is used to initialize the computation of the next optimum computation, which deals with finer details. This method increases the likelihood of finding the globally optimal match while avoiding the classical problem of the intensity-based method: their susceptibility to poor starting estimates.

As the regularization term, the intensity difference (first term H(I|Y) in the right hand side of the equation ( 2)) is weighted. The weight for the j-th inner iteration is defined as the inverse of the residual sum of the squared differences computed at the previous inner iteration of a given outer iteration of the algorithm and is denoted by 1 d j-1 . Because d has a high value for the first iterations, more weight is given to the regularization term, in order to get smooth deformations. As the algorithm gets close to the final solution, d theoretically tends to zero, giving less weight to the priors and letting the algorithm computing more detailed deformations.

Finally, the sought transformation ψ j,k i at the inner iteration j and the outer iteration k minimizes the function f j,k defined as:

f j,k 1 (ψ j,k i ) = 1 2d j-1 Ω ([G σ k * I 0 ](ψ j,k i (x)) -[G σ k * I i ](x)) 2 dΩ , (4)
f j,k 2 (ψ j,k i ) = g(J(ψ j,k i )), (5)
f j,k (ψ j,k i ) = f j,k 1 (ψ j,k i ) + λ f j,k 2 (ψ j,k i ), (6)
where the function g is computed from the Jacobian singular values and determinant [START_REF] Ashburner | Incorporating prior knowledge into image registration[END_REF]. Finally the N -1 deformations are computed through this iterative optimisation process. This approach was successfully applied before to large vessels as the aorta cross [START_REF] Midulla | Haemodynamic imaging of thoracic stent-grafts by computational fluid dynamics (CFD): presentation of a patient-specific method combining magnetic resonance imaging and numerical simulations[END_REF].

Patient-specific computational grid and application of the patient-specific deformation

Once the N -1 mappings ψ i are computed, a patient-specific computational grid must be extracted from the template image and warped thanks to the computed deformations. The template image corresponding to volumetric data I 0 is imported into an image processing software (ScanIP;Simpleware Ltd., Exeter, UK). Before segmentation, the image I 0 is prepared: the region of interest is isolated and the image is smoothed to erase noise inherent to the medical images protocols. The segmentation itself is done by a classical thresholding method [START_REF] Pham | Current methods in medical image segmentation[END_REF]. A suitable image intensity range which encompasses the voxel intensities of the region occupied by blood in the heart is selected. The quality of the segmentation relies on the operator's skills and knowledge of the human heart morphology and the quality of the chosen medical image. Image resolution plays a role in the potential domain simplification at this step. Trabeculae, left atrial appendage or papillary muscles can either be kept or neglected depending on the image quality and spatial resolution available. The 3D geometric reconstruction covers all the space occupied by blood in the left heart cavities. The surface of the geometric reconstruction of the heart is triangulated. Once a template 3D patient-specific surface mesh is created, a procedure to deform this surface model thanks to the images must be provided. For each couple of images (I 0 ,I i ) a suitable spatial transformation ψ i was found thanks to the method described in the previous section. These deformations ψ i are 3D deformation fields. Trilinear interpolation from these deformation fields to the template surface mesh is done. Thus, a set of N -1 successive surface meshes matching the physiological cardiac images at different times t i is produced as schematized in Fig. 1.

Position and velocity of all surface points are needed at any discrete time of the simulation, not only at the times t 0 ,t 1 , ...,t N-1 . Since all the generated surface meshes share the same topology (number and connection between nodes, number of cells), interpolation is used to compute the position of each node and the velocity by taking the time derivative of this quantity. As geometry variations are periodic, a trigonometric interpolation is used. The surface position and velocities read:

x s (t) = m ∑ i=0 [a i cos(2iπ t T ) + b i sin(2iπ t T )], u s (t) = 2iπ T m ∑ i=1 [-a i sin(2iπ t T ) + b i cos(2iπ t T )],          on ∂ Ω f (t) (7)
where T is the heart cycle period, m the number of Fourier modes (m = N-1 2 or N 2 depending on the parity of N) and a i , b i the Fourier coefficients. Surface velocity u s needed at the computation domain boundary ∂ Ω f (t) is hence not computed as a FSI problem, but entirely extracted from the medical images. In the present study, it is used to handle the cardiac chambers and their connected vessels. The template surface is imported in a commercial mesher (Gambit, ANSYS) to generate a template unstructured tetrahedral mesh. The computed boundary Fourier coefficients of Eq. ( 7) are interpolated in this template numerical domain. The com-Fig. 1 Mesh deformation procedure applied to a left human heart. The template mesh segmented from the image at time t 0 is deformed thanks to ψ i to obtain the mesh at time t i . This procedure is done for each image in the cardiac cycle in order to obtain the corresponding meshes.

putational mesh boundary now follows the shape of the patient endocardium and is updated in every step of the simulation. Motion of all internal points in the computational mesh is based on the prescribed boundary motion. At each iteration, nodal velocity u g is calculated through the computation of a Laplace equation [START_REF] Masud | An adaptive mesh rezoning scheme for moving boundary flows and fluid-structure interaction[END_REF] using the prescribed boundary motion as boundary condition for this problem.

The grid quality is monitored during the simulations. When the boundary displacement becomes too large compared to the local cell sizes, the cell quality can be highly deteriorated. This can lead to convergence problems or negative cell volumes. In this case, the numerical domain is re-meshed, Fourier coefficients are interpolated on this new discretized surface domain, as is the fluid solution at the new internal nodes.

Valves

Given their spatio-temporal resolution, MRI and CT scans fail to provide the necessary information to characterize precisely the movements of the aortic (AV) and of the mitral (MV) valves. Their geometry and movement have thus to be modelled. Accounting for the moving valves with the ALE method would be extremely complex due to grid quality issues. Here, an immersed boundary technique is used instead. A body force is imposed to drive the flow velocity to zero where the valves are located [START_REF] Chnafa | Image-based large-eddy simulation in a realistic left heart[END_REF].

The opening and closing valve time is generally small (of the order of 5% of the heart cycle [START_REF] Yoganathan | Fluid mechanics of heart valves[END_REF]). As a first modelling effort, it has been considered that valves open and close instantaneously, following the evolution of the ventricle volume. Systole lasts as long as the ventricle volume decreases. During systole, the AV is open and the MV closed. Diastole is defined as the heart period during which the ventricle volume increases: then, the AV is closed and the MV open.

Valve reconstruction starts with the definition of the valves annuli. Their geometry is reconstructed by inspecting the medical images. A number of markers are placed manually to define the valve annuli at a given time. The motion of these markers is then obtained thanks to the Fourier coefficient defining the boundary motion (Eq. 7).

As the focus is on the flow in the atrium and in the ventricle, the aortic valve is modelled very simply. Physiologically, the leaflets of the aortic valve are pushed against the vessel, offering small resistance to the flow. Hence, when the aortic valve is open, no immersed force is applied: aortic valve is only active when closed.

The mitral valve is represented by a more complex model, since its shape is expected to strongly impact the LV haemodynamics. The position of the mitral valve annulus being known over time, the mitral valve opening is defined using an ad hoc model. From visual inspection of the images, parameters as the average leaflets length, the orientation of the valve opening and an ellipse defining the opening area are imposed. Details of the procedure are described in Chnafa et al. [START_REF] Chnafa | Image-based large-eddy simulation in a realistic left heart[END_REF]. Knowing the MV leaflets position during the heart cycle, their effect on the blood flow is accounted for by using an immersed boundary method (IBM) [START_REF] Mohd-Yusof | Combined immersed-boundary/B-spline methods for simulations of flow in complex geometries[END_REF]. For this purpose, the leaflets representations are given a thickness so that a few mesh nodes are located within the valves. Then, the force f in the NS equations ( 1) is set to impose the fluid velocity to zero within the leaflets. Figure 2 shows an example of the valves models during diastole, when the MV is open.

Inlet boundary conditions

Inlet boundary conditions at the pulmonary veins must be provided. Under the present assumptions, either the aortic valve or the mitral valve is closed. Pulmonary veins thus always inject blood in a closed domain. As a consequence, blood being incompressible, the total mass flow rate entering the domain is entirely determined by the mass conservation constraint. With Q PV (t) denoting the inlet flow rate (sum of the volumetric flow rates at the four pulmonary veins) and V LA , V LV and V AO denoting respectively the LA, LV and the aortic root volumes, mass conservation yields:

Q PV (t) =          dV LA dt during systole, dV LA dt + dV LV dt during diastole. (8)
Uniform velocity profile are imposed at each inlet condition. With n j the outward normal vector and A j the area of inlet condition j ( j=1,..,4), the inflow velocity reads:

u j (t) = - Q i (t) A j ζ j n j , (9)
where ζ j determines the distribution of the flow between the four pulmonary veins.

In the absence of additional measurement (typically flow rates from PC-MRI), it is fair to assume a uniform distribution between the different pulmonary veins (ζ j =0.25 for all j.

3 Application to a patient left heart: problem formulation

Heart model

Using the framework described above, an actual CT exam of a patient treated at the University Hospital of Toulouse Rangueil (France) is used to generate a numerical domain and its deformation. The CT exam consists of N = 10 medical images along the cardiac cycle of period T = 1 s with a spatial resolution 2×2×2 mm. Due to the limited resolution of the images, the intra-cardiac geometry is simplified. As shown in Fig. 3, the numerical domain includes the LA, LV, the aortic root and four pulmonary veins. Left atrial appendum and geometrical details as the cordae tendinae or papillary muscles are omitted. The LA has a height of 5.5 cm from the MV to the upper pulmonary root and a maximum diameter of order 4 cm. The LV has a height of 8.8 cm from the MV to the apex (the lowest extremity of the LV in Fig. 3) and a maximum diameter of order 5 cm. The four pulmonary veins

Grid mesh and simulation details

A nearly isotropic grid is created from the heart model described in the previous section using the commercial software Ansys Gambit, which was selected for its ability to generate good-quality tetrahedral mesh, appropriate for non-dissipative finite-volume formulations. The spatial resolution is of order 0.8 mm, which yields grids of approximately three-million tetrahedral elements. In this study, in order to maintain the quality of the grid along the cycle, re-meshing (see section 2.2.2) was used at each acquisition instants t i , i = 1, ..., N -1.

The simulation time step is fixed by a CFL condition (CFL=0.9) consistent with the explicit time integration used in the CFD solver, which corresponds to a time step of order 10 -4 s.

Figure 4 displays the flow rates at the aortic valve (top plot), mitral valve (middle plot) and the heart inflow (bottom plot) gathering the four PV. Two verticals dotted lines separate the systolic phase and the diastolic one. For this heart, systole lasts t/T = 0.36 (from t/T 0.015 to 0.375) and diastole t/T = 0.64.

Q MV (mL.s -1 ) Q AO (mL.s -1 ) Q PV (mL.s -1 )
Fig. 4 Flow rates at the aortic valve (top plot), mitral valve (middle plot) and the total heart inflow (bottom plot) imposed at the pulmonary veins during the heart cycle. Vertical dotted lines mark the limit between the systolic phase (t/T between 0.015 and 0.375) and the diastolic phase.

Element

U max D = 2 A π Re max t m /T AV 0.96 m.s -1 2.20 cm 5300 0.16 MV 1.08 m.s -1 1.87 cm 5000 0.52 Pulm. vein 0.79 m.s -1 1.00 cm 2000 0.52 Table 1 Main flow parameters describing the simulation. The section-averaged maximum velocity is indicated as U max . Reynolds numbers are based on the diameter D of the region of interest. The instants when they reach their maximum value Re max are reported as t m /T . For valves, A is the area of the lumen when open.

The pulmonary flow rate derived from time evolution of the heart cavities is consistent with classical medical dataThe aortic valve flow rate behaves as expected: it increases during systole with a maximum flow rate of Q MV = 320 mL.s -1 , then decreases until its shutting at t/T = 0.375. The aortic flow rate stays null during diastole. The mitral flow rate is usually composed by two peaks. The E wave, or rapid filling, and a second one, the A wave corresponding to late diastole. They are separated by a phase with almost no heart motion called diastasis. The flow rate shows only one main peak in the present case: the E wave (t/T = 0.51, Q MV = 410 mL.s -1 ). The fact that the A wave is weak is symptomatic of pathologies.

The blood properties are: ρ = 1040 kg.m -3 and ν = 4 × 10 -6 m 2 s -1 . Relevant flow conditions at different locations of the heart are reported in Table 1. Inlet Reynolds number for each pulmonary vein varies from 0 to approximately 2000, based on the vein diameter. The Reynolds number at the mitral tips varies from 0 to approximately 5000 (based on the effective mitral mean diameter D = 2R a = 2 √ ab, the area of the open MV, the kinematic viscosity and the maximum flow rate). The maximum transmitral velocity U max falls into the usual measurements (around 1.0 m.s -1 ) [START_REF] Haugen | Velocity profiles in mitral blood flow based on three-dimensional freehand colour flow imaging acquired at high frame rate[END_REF]. The Reynolds number of the aortic valve is about 5300. Table 1 reports the maximum velocity magnitude and the maximum Reynolds number for different elements of the heart model: pulmonary veins, mitral valve and aortic valve. These ranges of Reynolds numbers and the pulsating nature of the inlet flow clearly indicate that this complex cyclic flow may be transitional if not fully turbulent. This justifies the use of Large-Eddy Simulation (LES) which is more suitable than other simulation approaches for unsteady and/or transitional flows. In order to achieve a quality assessment of the presented simulation, the Pope criterion [START_REF] Pope | Ten questions concerning the large-eddy simulation of turbulent flows[END_REF] is used. According to this criterion, a reliable LES should be able to resolve at least 80% of the turbulent kinetic energy. Looking at the phase where the turbulent activity is the highest (t/T =0.65), it was found that this criterion is met in 85% of the numerical domain. The last 15% are mainly located in the atrial cavity.

Results and discussion

Global description of the cardiac cycle

Fifteen cardiac cycles were simulated and phase-averaged in order to provide a clear view of the flow organization over the cardiac cycle. A detailed flow description can be found in [START_REF] Chnafa | Image-based large-eddy simulation in a realistic left heart[END_REF]; only the main flow characteristics are recalled here for the sake of completeness. Six salient instants are illustrated in Figure 5: the ventricular mid and end-systole (t/T =0.25-0.35), the beginning, peak and end of the E wave (t/T =0.45-0.55-0.65) and the end of the A wave, just before the beginning of the next ventricular systole (t/T =0.99). The velocity field is scaled by u a = ql s /V s = 0.1 ms -1 where q=7.50 × 10 -5 m 3 s -1 is the cardiac output, V s = 5.55 × 10 -5 m 3 is the end systolic volume and l s =7.40×10 -2 m is the ventricle length at the end of Fig. 5 Phase-averaged velocity field over a cutting-plane through the left heart. Velocity vector scale is not constant though the heart cycle and is indicated for each plot. Mitral valve is depicted in light grey the aortic one in dark grey.

systole. Note that due to the strong velocity variations along the cycle, the vector scale was adapted for each instant.

During systole (t/T =0.25-0.35), the mitral valve is closed (in light gray in Fig. 5), preventing backflow towards the atrium, while the aortic valve is open (in dark gray in Fig. 5). The decrease of the ventricle volume causes ejection of blood into the aorta. Note that the velocity amplitude is the highest in the ascending aortic root. The computed flow at mid-systole is also highly swirled in the atrium, as reported in vivo [START_REF] Kilner | Asymmetric redirection of flow through the heart[END_REF][START_REF] Markl | Comprehensive 4D velocity mapping of the heart and great vessels by cardiovascular magnetic resonance[END_REF]. This movement is hardly discernible in Fig. 5, because the vortical movement is mainly perpendicular to the cutting plane [START_REF] Chnafa | Image-based large-eddy simulation in a realistic left heart[END_REF]. Two recirculation zones are however visible in Fig. 5 (top center): just under the MV and within the aorta, just above the AV, in agreement with [START_REF] Kilner | Helical and retrograde secondary flow patterns in the aortic arch studied by three-directional magnetic resonance velocity mapping[END_REF]. At the end of the ventricle contraction, the aortic valve closes and the mitral valve opens: ventricle filling starts. At t/T =0.375 (not shown), the ventricle diastole starts: the LV volume increases and blood passes from the LA to the LV, forming a strong jet through the MV. The shear layer between the jet generated during the E wave and the surrounding quiescent fluid rolls-up and shapes the jet head as a vortex ring [START_REF] Le | On the three-dimensional vortical structure of early diastolic flow in a patient-specific left ventricle[END_REF]. The E wave vortex ring signature and its evolution are visible in Fig. 5 (t/T =0.45-0.55) which illustrates a mechanism already reported by several authors [START_REF] Domenichini | Three-dimensional filling flow into a model left ventricle[END_REF]. The vortex ring does not remain symmetric, as the lateral wall prevents its full development. A large recirculating cell is clearly visible in the LV at t/T =0.65 (Fig. 5, bottom center), as often described in the literature [START_REF] Hong | Characterization and quantification of vortex flow in the human left ventricle by contrast echocardiography using vector particle image velocimetry[END_REF][START_REF] Long | Subject-specific computational simulation of left ventricular flow based on magnetic resonance imaging[END_REF][START_REF] Markl | Comprehensive 4D velocity mapping of the heart and great vessels by cardiovascular magnetic resonance[END_REF]. It is characteristic of the flow in the ventricle after the E wave. Two less intense blood recirculation zones can be detected: one at the apex, which is visible during the whole diastole and an intermittent one between the aortic valve and one of the MV leaflet. These blood recirculations are also described in silico [START_REF] Doenst | Fluid-dynamic modeling of the human left ventricle: Methodology and application to surgical ventricular reconstruction[END_REF][START_REF] Mihalef | Patient-specific modelling of whole heart anatomy, dynamics and haemodynamics from fourdimensional cardiac CT images[END_REF][START_REF] Schenkel | MRI-based CFD analysis of flow in a human left ventricle: Methodology and application to a healthy heart[END_REF]. Between the E wave and the A wave, the recirculating cell core in the LV moves from the ventricle center to the septum wall. During the A wave occurring at t/T =0.99 (Fig. 5, bottom right), the blood flux passing though the MV strengthens the recirculating cell in the LV, as classically reported [START_REF] Pedrizzetti | Nature optimizes the swirling flow in the human left ventricle[END_REF]. Atrial contraction expels blood from the LA, both through the MV, as seen in the lower half part of the LA and through the pulmonary veins, as shown by the upward velocity vectors visible in the upper half part of the LA.

Cycle-to-cycle variations

The studied flow configuration is a breeding ground for weak turbulence. The transitional nature of this cyclic flow due to the highly complex evolving geometry and the unsteady inflows results in cycle-to-cycle variations. As an illustration, Figure 6 shows the time evolution of the vertical velocity at four different points within the left heart. Velocity signals corresponding to six different heart cycles are superimposed to visualize the cycle-to-cycle variations. Moderate (2u a ) to high (5u a ) variations are present at all the locations, although not uniformly distributed over the cardiac cycles. The second plot (corresponding to a probe located under the aortic valve) shows the lowest variations: almost no cycle-to-cycle variations are visible during the diastolic phase, which is expected, given the quiescent flow beneath the aortic valve during diastole. Cycle-to-cycle variations are however visible (variations around 2u a ) during the systolic phase, as expected given the typical Reynolds number. Cycle-to-cycle variations decrease during the systolic phase in the LV (see third plot), beginning from variations up to 4u a to almost null variations. This absence of fluctuations carries on even after the passage of E-wave vortex ring, visible on the signal. Variations are visible from t/T = 0.55 and amplify during the late diastole. A similar behaviour is visible in the lower part of the LV (see fourth plot): during systole, variations decrease, then rise after t/T = 0.55, reaching an amplitude of order 5u a . The largest relative fluctuations are obtained in the left atrium (top plot) with small cycle-to-cycle variations during systole (t/T between 0.015 and 0.375) but fluctuations as large as 2u a at t/T = 0.35 and 5u a at t/T = 0.6. Theses variations are related to the interaction between the four inlet flows from the four pulmonary veins.

A more quantitative assessment of the cycle-to-cycle variations is obtained by computing the kinetic energy of the velocity fluctuations. The velocity flow is decomposed in a mean (phase-averaged) velocity field u and a fluctuating part u ′ . The turbulent kinetic energy k and the mean flow kinetic energy E k are defined as:

k(t) = 1 2V (t) V (t) (u rms (x,t) 2 + v rms (x,t) 2 + w rms (x,t) 2 ) dV, (10)
E k (t) = 1 2V (t) V (t) (u(x,t) 2 + v(x,t) 2 + w(x,t) 2 ) dV, (11)
where u rms , v rms and w rms are the root-mean-square values of the velocity fluctuations in the three directions. The volume V (t) is either the volume of the ventricle or the volume of the atrium at time t.

Figure 7 shows how these energies evolve over the heart cycle, both in the ventricle (top plot) and the atrium (bottom plot). Note that 5 times k is plotted, so that the same scale is used to represent both energies. The mean flow kinetic energy evolves similarly in both cavities. During the systolic phase (t/T between 0.015 and 0.375) it increases and reaches a plateau at t/T =0.15. After a decrease just before the beginning of diastole, the mean flow kinetic energy increases again and reaches a peak a . Vertical dotted lines mark the limit between the systolic phase (t/T between 0.015 and 0.375) and the diastolic phase. at t/T =0.54, viz. 0.04T after the the peak of the E wave. The maximum value of E k corresponds to the presence of high velocities when blood flows from the atrium to the ventricle (see the mitral jet in Fig. 5, bottom left). The ventricular turbulent kinetic energy k remains low during the systolic phase thanks to the stabilizing effect of the flow acceleration, with values of less than 5% of E k . It then increases substantially, reflecting the amplification of the disturbance after the jet impingement on the lateral ventricle wall. The turbulent energetic peak is reached T /10 after the peak of mean flow kinetic energy, corresponding to the convection time of the vortex ring and the decelerating phase of the flow. The turbulence intensity k/E k in the ventricle is as high as 50 % during the k peak.

The atrial turbulent kinetic energy behaves somewhat differently. First, it increases during the whole systolic phase, because of the interaction/collision of the four inflowing jets issued from the pulmonary veins. A first peak is thus reached near the beginning of the diastolic phase and turbulence intensity is then attenuated during the flow acceleration through the atrium, as expected. The atrial turbulent kinetic energy rises again after t/T =0.5 and reaches its peak before t/T = 0.6 during the flow deceleration. On top of occurring earlier in the heart cycle, this peak is around twice less energetic than the one occurring in the ventricle. Still, it corresponds to a large turbulence intensity of approximatively 20 %. As in the ventricle, the turbulent kinetic energy then decreases until the end of the heart cycle.

Conclusion and outlook

The approach presented here allows patient-specific blood flow simulations in the heart from a series of gated 3D images. Starting from 4D medical images, the numerical domain is first extracted and the heart wall movements are then calculated thanks to a proper image registration algorithm. In order to demonstrate the ability of the method to reproduce the cardiac flow, a computation of the blood flow in a whole left heart has been conducted. Results consistent with the current knowledge in terms of left heart flow is presented. All presented features have been reported several times in the literature, both in numerical and experimental studies and by medical imaging.

Furthermore, the use of fluid numerical method well adapted to fluctuating turbulent flows enables the observation of cycle-to-cycle variations in the flow field. Such variations are expected in the present flow, due to the high Reynolds numbers encountered and the unsteadiness of the flow incoming from the pulmonary veins. The present results show that in spite of rigorously identical contraction and boundary conditions, fluid inertia makes the flow differ from one cycle to another. More precisely, cycle-to-cycle variations in the left atrium can be observed in its upper part, where the collision of the jets issuing from the pulmonary veins makes the flow particularly chaotic. Spatially averaged kinetic turbulent energy level reaches a turbulent intensity of 20 % at its peak then slowly decreases. In the left ventricle, velocity fluctuations are reported mainly during late diastole. Between the impact of the E wave jet on the lateral wall and the end of diastole, the left ventricle displays high levels of cycle-to-cycle fluctuations. Indeed, both the vortex ring impact and the E wave deceleration occur approximately at the same time, and both are features tending to generate turbulence. This translates into turbulent intensity levels as high as 50 %. This turbulent activity dissipates little by little, until the flow acceleration at early systole.

The presented method does not include or simplifies some aspects of the physiological heart. Blood is considered as a Newtonian fluid, which is an approximation commonly accepted for the heart flow. A non-Newtonian model could be included in our simulations. The spatio-temporal resolution of medical imaging imposes temporal interpolation and geometrical simplifications of the heart model. Consistently with the poor time resolution of the input medical data, a rough model of the mitral valve was used, which constitutes the main drawback of our method. Valves are instantaneously switch from their closed position to their open position and vice versa. This simple approach is justified by the fact that the opening and closing of the valves last only 5% of the whole cycle. Change of aperture area along time is not considered either. Note however that the model can be refined depending on the available data. Nevertheless, the presented approach allows to retrieve features reported in the literature and, in addition, it can provide detailed information about the complex intermittent flow in the left heart.

Fig. 2

 2 Fig. 2 Left plot: annulus markers and grid cells where a body force is added to model the MV leaflets. Diameters a(t) and b(t) of the elliptic opening are indicated. Right plot: full MV model included in the numerical domain (the AV is represented as well).

Fig. 3

 3 Fig. 3 Template computational domain extracted from a 3D medical image. The same domain is represented for four different points of view and the left ventricle (LV), left atrium (LA) and Aorta (AO) are indicated. Black line passing through the left heart indicates the position of slices used to describe the flow in section 4.

Fig. 6

 6 Fig. 6 Temporal evolutions of the scaled vertical velocity w/u a (the w direction is indicated in the figure and u a = ql s /V s ) at four different points within the left heart. Six cycles are reported to illustrate the cycle-to-cycle variations.

2 aFig. 7

 27 Fig.7Volumetric mean flow kinetic energy E k (full line) and five times the turbulent kinetic energy k (dashed line) in the left ventricle (top plot), and in the left atrium (bottom plot). The energies are nondimensionalised by u 2 a . Vertical dotted lines mark the limit between the systolic phase (t/T between 0.015 and 0.375) and the diastolic phase.

Acknowledgments

The authors thanks Prof. D. Coisne and Dr. S. Nottin for interesting discussions about the presented results and their medical/physiological concordance. Dr. V. Moureau is gratefully acknowledged for giving access to the YALES2 solver.


Rahaf Aljundi
email: rahaf.aljundi@gmail.com

Rémi Emonet

Damien Muselet
email: damien.muselet@univ-st-etienne.fr

Marc Sebban
email: marc.sebban@univ-st-etienne.fr

Landmarks-based Kernelized Subspace Alignment for Unsupervised Domain Adaptation

Domain adaptation (DA) has gained a lot of success in the recent years in computer vision to deal with situations where the learning process has to transfer knowledge from a source to a target domain. In this paper, we introduce a novel unsupervised DA approach based on both subspace alignment and selection of landmarks similarly distributed between the two domains. Those landmarks are selected so as to reduce the discrepancy between the domains and then are used to non linearly project the data in the same space where an efficient subspace alignment (in closed-form) is performed. We carry out a large experimental comparison in visual domain adaptation showing that our new method outperforms the most recent unsupervised DA approaches.

Introduction

While the standard machine learning setting assumes that training and test data come from the same statistical distribution, it turns out that many real world applications, e.g., in Computer Vision or in Natural Language Processing, challenge this assumption. Indeed, the (source) data from which the classifier is supposed to be learned often differ from the (target) data the classifier will be deployed on (see Figure 1). To deal with such situations, the learning system has to take into account the distribution shift between the two domains in order to adapt well from the source to the target. In this framework, two main categories of domain adaptation (DA) methods are available. The first one, called semi-supervised DA, can gain access to some labeled examples from the target domain (as well as labeled source data) in order to perform the adaptation process (see, e.g., [6,1,[START_REF] Mehrotra | Dictionary based sparse representation for domain adaptation[END_REF]). On the other hand, unsupervised DA assumes that no labeled information is available from the target domain. To deal with this more complex setting, DA Figure 1. Example of distribution shift between images from two datasets. First row, some bike helmets from the Amazon subset and second row, bike helmets from the webcam subset. These 2 subsets are from the Office dataset. methods (see [START_REF] Margolis | A literature review of domain adaptation with unlabeled data[END_REF] for a survey) usually assume that the distribution shift can be compensated either by a reweighting scheme applied on the source data (so as to better match the target distribution), or by a projection of both source and target data in a common (possibly latent) space. In this paper, we focus on this second category of DA methods which have the advantage of not assuming that source and target data originally lie in the same feature space.

It is worth noting that the most promising recent approaches rely on the assumption of the existence of a common embedded low dimensional manifold space that minimizes the divergence between the two distributions. In this context, subspace alignment-based DA methods have attracted a lot of interest. For example, in [START_REF] Gopalan | Domain adaptation for object recognition: An unsupervised approach[END_REF][START_REF] Gong | Geodesic flow kernel for unsupervised domain adaptation[END_REF], source and target data are projected onto intermediate (linear) subspaces that lie along the shortest geodesic path connecting the two original spaces. However, although effective, the search for such subspaces is computationally costly and subject to local perturbations. A recent work done by Fernando et al. [7] overcomes these drawbacks by optimizing a single linear mapping function that directly aligns the source and target subspaces. This new method has shown to be not only better than the state of the art but also computable in closed form. In this paper, we aim at going a We can see that no linear transformation will be able to entirely bridge the source to the target. On the other hand, we can note that, for each class, some subsets of source and target examples seem to be similarly distributed (in black). The objective is to automatically discover those so-called landmarks that will be used to map both source and target data onto a common space using a kernel.

step further by improving this approach which faces two main limitations: First, it assumes that the shift between the two distributions can be corrected by a linear transformation. We claim that this is a strong assumption that can be challenged in many real world applications. Furthermore, we show in the experimental section that performing a simple kernel PCA [2] does not solve the problem. Second, it assumes that all source and target examples are necessary to proceed to the adaptation. However, it turns out that in most of the cases only a subset of source data are distributed similarly to the target domain and vice versa. In this paper, we deal with these two issues by: • (i) Proposing a selection of landmarks extracted from both domains so as to reduce the discrepancy between the source and target distributions.

• (ii) Projecting the source and target data onto a shared space using a Gaussian kernel w.r.t. the selected landmarks. This allows us to capture the non linearity from the data in a simple way.

• (iii) Learning a linear mapping function to align the source and target subspaces. This is done by simply computing inner products between source and target eigenvectors.

As far as we know, this paper is the first attempt to select landmarks in an unsupervised way in order to reduce the discrepancy between the two domains. The intuition behind our method is summarized in Figure 2. The rest of this paper is organized as follows. Section 2 is devoted to the presentation of the related work. In Section 3, we introduce our unsupervised DA method based on a landmark selection and a subspace alignment. Section 4 is dedicated to the experimental comparison performed on computer vision datasets. It highlights the interest of our approach which significantly outperforms the state of the art DA methods.

Related Work

Domain adaptation (DA) [START_REF] Margolis | A literature review of domain adaptation with unlabeled data[END_REF] is one of the possible settings of transfer learning [21]. In this paper, we focus on the most complicated case, called unsupervised DA, where labeled data are available only from the source domain. In this context, in order to reduce the shift between the source and target distributions, one strategy consists in reweighting each source data according to its similarity with the target domain [4, [START_REF] Chu | Selective transfer machine for personalized facial action unit detection[END_REF][START_REF] Chatpatanasiri | A new kernelization framework for mahalanobis distance learning algorithms[END_REF]23,24]. Then, a classifier is learned (typically an SVM) on these weighted source instances. Another approach aims at projecting both source and target data in a common space which tends to move the two distributions closer to each other [START_REF] Pan | Domain adaptation via transfer component analysis. Neural Networks[END_REF][START_REF] Long | Transfer sparse coding for robust image representation[END_REF]3,[START_REF] Pan | Transfer learning via dimensionality reduction[END_REF]. This latter strategy has attracted a lot of interest during the past few years and gave rise to the development of subspace alignment-based DA methods. In [START_REF] Gopalan | Domain adaptation for object recognition: An unsupervised approach[END_REF], the authors learn a sequence of intermediate subspaces along the geodesic path connecting the source and target subspaces (see Figure 3). A similar work done by Gong et al. [START_REF] Gong | Geodesic flow kernel for unsupervised domain adaptation[END_REF] also leans toward the same idea by proposing a geodesic flow kernel that tries to model incremental changes between the source and target domains. However, it is clear that computing such a large number of subspaces is computationally expensive. Note that recent approaches have tried to combine both source instance reweighting and distribution projection in a common space [START_REF] Long | Transfer joint matching for unsupervised domain adaptation[END_REF]17,[START_REF] Gu | Joint feature selection and subspace learning[END_REF]. The very recent paper of Long et al. [START_REF] Long | Transfer joint matching for unsupervised domain adaptation[END_REF] is the most noticeable among them since their Transfer Joint Matching (TJM) algorithm provides the current state-of-the-art results. TJM aims at learning a new space in which the distance between the empirical expectations of source and target data is minimized while putting less importance to the source instances that are irrelevant to classify the target data. For this purpose, the authors make use of a kernel mapping of all the examples allowing by this way a non-linear transform between the two domains. It is worth noticing that, despite its efficiency, this approach suffers from several drawbacks. First, at testing time, the kernel has to be computed between each new instance and all the source and target examples, that can be very costly. Second, since the idea is to minimize the distribution difference while reweighting the source instances, the objective function is very complex to optimize and the authors resort to an alternating optimization strategy by iteratively updating one variable while the other ones are fixed. Third, because of this already complex optimization process, only the means of the source and target distributions are brought closer without accounting their standard deviations. We claim that the dispersions from the means can vary a lot between domains and should be adapted as well. This paper deals with all these issues: (i) it defines via a (cheap) greedy approach the landmarks that reduce the discrepancy between the two distributions; (ii) the mean and variance of the distribution of each landmark are taken into account in the adaptation process; (iii) the subspace alignment is performed in a closed form that only requires inner products between source and target eigenvectors (as suggested in [7]).

Since our method is directly inspired from [7], let us go into details of this one-step subspace alignment method (denoted by SA in the rest of this paper). Let us consider a labeled source dataset S and an unlabeled target dataset T lying in a D dimensional space. SA builds the source and target subspaces, respectively denoted by X S and X T , by extracting d eigenvectors (corresponding to the d largest eigenvalues, d ≤ D) from both domains thanks to two independent PCA. Since X S and X T are orthonormal, the authors show that the optimal matrix M that transforms the source subspace into the target one is M = X S X T . Note that since X S and X T are generated from the first d eigenvectors, they are intrinsically regularized. This explains why SA is extremely fast. Furthermore, the authors show that the unique hyperparameter d can be tuned by resorting to theoretical results inspired by concentration inequalities on eigenvectors.

As pointed out in the introduction, although very efficient, SA suffers from two main drawbacks. First, it allows only linear transforms between the source and target domains whereas most of the real transforms are non-linear (see Fig. 2). Second, forcing the algorithm to adapt all the source data to the target domain is a strong constraint while in most cases only a subset of source and target examples are similarly distributed.

In this paper, we cope with these two issues while preserving the theoretical guarantees and speed properties of SA. First, we automatically detect the subset of source and target examples that allow us to bring closer the two dis-tributions. Then, we kernelize the problem w.r.t. those selected landmarks before performing the subspace alignment in closed-form. Using landmarks with a local kernel allows us to annihilate the impact of the points that are far from all landmarks.

Finally, note that the use of landmarks in domain adaptation has been recently used in [START_REF] Gong | Connecting the dots with landmarks: Discriminatively learning domain-invariant features for unsupervised domain adaptation[END_REF] in a semi-supervised DA algorithm. Beyond the fact that we are working in this paper in a totally unsupervised setting, we select our landmarks from both sources while [START_REF] Gong | Connecting the dots with landmarks: Discriminatively learning domain-invariant features for unsupervised domain adaptation[END_REF] only picks them from the source domain.

Proposed Approach for Domain Adaptation

In this section, we provide a detailed description of our landmark-based subspace alignment method for domain adaption. Our method is fully unsupervised in the sense that no labels are required to perform the domain adaptation. Labeled points from the source domain are only used afterwards to learn a classifier.

Task and Overall Approach

Source S and target T points are supposed to be respectively drawn from a source distribution D S and a target distribution D T . Domain adaptation supposes that the source and target distributions are not identical but share some similarities that make it possible to adapt to the target domain what has been learned on the source domain. Said differently, if we have a set of labels L S for the source examples (all of them or only some of them), they can be used to learn a classifier that is suitable for the target domain.

Our approach combines two simple ideas: First, it projects both source and target examples in a common subspace w.r.t. some well selected landmarks. Then, it performs a subspace alignment between the two domains. After selecting landmarks among S ∪ T , all points in S and T are projected using a Gaussian kernel on the selected landmarks, leading to new representations K S and K T for the source and target points. The new representation is then used to compute a mapping using a subspace alignment approach.

Compared to [7], our two-step approach remains fast and easy to implement while improving the accuracy by capturing non-linearity. The rest of this section details each step involved in our approach from the multiscale landmark selection process to the projected subspace alignment and the classification. The complete pseudo code is given in Algorithm 2 that shows the successive steps of our approach, also illustrated in Fig. 4.

Multiscale Landmark Selection

The first step of our method is to select some points as landmarks. Intuitively, a good set of landmarks is a Figure 4. Overall workflow of our approach: First, landmarks are selected from S ∪ T so as to maximize the overlapping between the source and target distributions; Second, we get two new representations KS and KT for the source and target points using on Gaussian kernel on the selected landmarks; Then, two independent PCA are carried out before performing the subspace alignment w.r.t the d first eigenvectors; Finally, a classifier is learned (typically an SVM) from the labeled source data and deployed on the target domain.

set of points which can be used to project the source and target data onto a shared space in which their distributions are more similar. Our method chooses the landmarks among S ∪ T and does not use any label. The final output of the landmark selection procedure is a set of landmarks A = {α 1 , α 2 , ...}, A ⊂ S ∪ T . To avoid any computationally expensive iterative optimization method, we propose a direct method of asserting whether a point should be kept as a landmark or not.

Landmark selection overview Our landmark selection approach actually considers each point c from S ∪ T as a candidate landmark. Each candidate is considered independently from the others: a quality measure is computed and if this measure is above a threshold, the candidate is kept as a landmark.

To assess the quality of the candidate c, its similarity with each point p ∈ S ∪ T is first computed, using a Gaussian kernel with standard deviation s:

K(c, p) = exp -c -p 2 2s 2 (1)
The quality measure for the candidate c is computed as an overlap between the distribution of the K(c, .)-values for the source points, and the one for the target points. As a summary, a landmark is considered as a good one if the distribution of the source and the one of the target points are similar after projection using the kernel.

Multiscale analysis

The value of the kernel radius s from Eq. 1 is important as it defines the size of the neighborhood that the candidate landmark considers. Choosing the right s for a given landmark will allow us to capture the local phenomena at the proper scale and to better align the source and target distributions. Extreme values for s must be avoided as they lead to a perfect match between the distributions of source and target points: all values of K(c, .) become 0 (when s is close to 0) or 1 (when s is very big).

When calculating the quality of a candidate landmark, we actually perform a multiscale analysis: we select the best scale s to capture the local properties of the data, at the same time avoiding extreme values for s. To do so, we compute the distribution of euclidean distances between all pairs of points and try every percentile value of this distribution. With this percentile-based approach, we consider a range of values for the scale s which are all plausible. For each considered scale s, we compute an overlap measure between the source and target distributions, keeping the greatest one as the quality measure for the candidate landmark.

Distribution overlap criteria For a candidate landmark c and a scale s, we want to compute a degree of overlap between two sets of K(c, .)-values: the one for the source points KV S and the one for the target points KV T . To lower the computational cost, the two distributions are approximated as normal distributions and thus summarized by their means and standard deviations µ S , σ S , µ T , σ T . To be able to use a fixed threshold and give so semantics to it, we use a normalized overlap measure, explained below and computed as follows:

overlap(µ S , σ S ; µ T , σ T ) = N (µ S -µ T | 0, σ 2 sum ) N (0 | 0, σ 2 sum ) (2)
Where

σ 2 sum = σ 2 S + σ 2 T and N (. | 0, σ 2 sum ) is the centered 1D normal distribution.
To compute the overlap of two probability densities f and g, we can integrate their product: f (x)g(x)dx. Intuitively, the more overlap there is between the densities, the bigger their product. This overlap measure follows similar principles as the Bhattacharyya coefficient. By analytically integrating the product of two normal distributions, we obtain the numerator of Eq. 2:

N (x | µ S , σ 2 S )N (x | µ T , σ 2 T )dx = N (µ S -µ T | 0, σ 2 sum )
(3) The denominator in Eq. 2 corresponds to the maximum value of the numerator for a given σ sum (obtained when µ S = µ T ). The denominator acts as a normalization factor: it sets the overlap to 1 when the distributions are perfectly matching and gives an easier interpretation, helping in the choice of the threshold th.

Summary Algorithm 1 sums up the landmark selection process. Each point from S ∪ T is considered as a candidate landmark. For each candidate, we consider multiple scales s. If, for any of these scale, the overlap (Eq. 2) between the source and target distributions of K(c, .)-values is greater than th, the candidate is promoted as a landmark.

Algorithm 1 MLS: Multiscale landmarks selection, implementing choose landmarks in Algorithm 2. Require: Source data S, Target data T , Threshold th. Ensure: A contains the selected landmarks.

A ← {} distances ← { a -b , (a, b) ∈ (S ∪ T ) 2 } for c in S ∪ T do for s in percentiles(distances) do KV S ← {exp(-c -p 2 /2s 2 ), p ∈ S} KV T ← {exp(-c -p 2 /2s 2 ), p ∈ T } if overlap(KV S , KV T ) > th then A = A ∪ {c} end if end for end for

Kernel Projection and Subspace Alignment

Once the set of landmarks A has been selected, we use a Gaussian kernel to achieve a non-linear mapping of all the points into a common space defined by these landmarks. The subspaces from the source and the target domains are then aligned using a linear transformation.

Projection on the selected landmarks Each point p from S ∪ T is projected onto each landmark α j ∈ A using a Gaussian kernel with standard deviation σ:

K(p, α j ) = exp -p -α j 2 2σ 2 (4)
Overall, the points from both S and T are projected in a common space that has as many dimensions as there are landmarks. Following other non-linear methods we set σ to the median distance between any pair of points drawn randomly from S ∪ T . The value of σ could also be selected using some cross-validation procedure. After projection, we finally obtain new representations for the source and target, respectively denoted by K S and K T .

Subspace alignment Using the non-linearly projected set of points K S and K T , we follow a subspace alignment approach. PCA is applied on each domain separately to extract the d eigenvectors having the largest eigenvalues. Following the theoretical part of [7], we are able to determine the optimal value of d. Indeed, the authors in [7] derive a consistency theorem based on a standard concentration inequality on eigenvectors. This theorem allows them to get a bound on the deviation between 2 successive eigenvalues.

In this paper, we make use of this bound to efficiently tune d, the number of dimensions in PCA. The d eigenvectors for the source and the target domains are respectively denoted X S and X T . The points from each domain can be projected on their respective subspace as K S X S and K T X T . The goal of subspace alignment is to find a linear transformation M that best maps the source eigenvectors onto the target eigenvectors. In the end, we want to find M that minimizes the sum of the euclidean distances between the M-transformed source eigenvectors and the target eigenvectors. This minimization is equivalent to minimizing the following Frobenius norm:

F (M ) = X S M -X T 2 F . (5)
Equation 5 has a closed form solution that can be implemented efficiently as a multiplication of the eigenvectors:

M = X S X T . (6)
The alignment transformation M maps points from the source eigenspace to the target eigenspace. The transformation M can be used to bring the projected source points K S X S into the same eigenspace as the projected target points K T X T , by computing K S X S M . The complete pseudo-code of our algorithm, called LSSA (for Landmarks Selection-based Subspace Alignment), is described in Algorithm 2.

Experiments

The objective of this experimental section is twofold. First, we aim at studying the behavior of our landmark selection method, taken alone, to deal with unsupervised visual domain adaptation in comparison with other landmark selection approaches. Second, we seek to show that jointly used with a subspace alignment process, our landmark selection method yields to a significant enhancement of the adaptation process and allows us to outperform the state of [START_REF] Gong | Geodesic flow kernel for unsupervised domain adaptation[END_REF][START_REF] Gopalan | Domain adaptation for object recognition: An unsupervised approach[END_REF]7,[START_REF] Kulis | What you saw is not what you get: Domain adaptation using asymmetric kernel transforms[END_REF][START_REF] Saenko | Adapting visual category models to new domains[END_REF] for generating the source and target samples.

As indicated before, we perform the following two series of experiments.

Comparison of landmark selection methods

To achieve this task, we compare our method (MLS in Table 1) with three baselines:

• A random selection (RD). For this purpose, we randomly select 300 landmarks (150 for each domain) and repeat the task five times to get a behavior on average.

• No landmark selection (All). We use all the source and target examples as landmarks.

• Our method without the multiscale strategy (σ-LS). The same standard deviation σ is used for all candidates. σ is set (to the advantage of this baseline) to the best standard deviation allowing on average the maximum overlapping between the two distributions.

Note that the first two baselines do not perform any adaptation because they do not aim at specifically moving the distributions close to each other (the model is simply learned from the source data and deployed on the target domain). For MLS and σ-LS, we fix the overlapping rate to 0.3 to select a landmark. Because of the normalization, a threshold of 0.3 corresponds to an overlapping rate of 30%. We also compare MLS with another landmark selection method, called Connecting Dots with Landmarks (CDL), recently presented in [START_REF] Gong | Connecting the dots with landmarks: Discriminatively learning domain-invariant features for unsupervised domain adaptation[END_REF] in a semi-supervised DA setting.

Comparison with the state of the art unsupervised DA approaches The second series of experiments is devoted to the evaluation of our landmark selection method used in an unsupervised subspace alignment DA setting. We compare our approach (LSSA), as presented in Algorithm 2, with three state of the art methods:

• Geodesic Flow Kernel (GFK) [START_REF] Gong | Geodesic flow kernel for unsupervised domain adaptation[END_REF], where a sequence of intermediate subspaces is learned along the geodesic path connecting the source and target domains.

• The one step subspace alignment (SA) proposed by [7] which learns a linear transformation between the two subspaces.

• The Transfer Joint Matching (TJM) presented in [START_REF] Long | Transfer joint matching for unsupervised domain adaptation[END_REF], which is a recent work based on both feature matching and instance weighting.

In addition, we perform experiments with two baselines: the first one does not perform any adaptation (NA); the second one performs two independent KPCAs on the source and target domains and then learns the linear transformation using the algorithm SA [7] (denoted by KPCA+SA).

Analysis of the Results

Comparison of landmark selection methods From the results reported in Table 1, we can make the following comments. First, we can note that our method significantly (using a Student paired-t test) outperforms on average the other approaches (with an average accuracy of 48.1%). Among 12 DA tasks, MLS has the best behavior in terms of accuracy in 8 subproblems. Second, for two subproblems (W→D and D→W), All is better that means that keeping all the source and target examples in these two (symmetric) cases is better than trying to select landmarks. It is worth noting that the two corresponding subproblems are the simplest ones (with the highest accuracy) justifying the interest to keep all the data. Moreover, we can note that in 10 subproblems out of 12, our approach outperforms CDL which is specifically dedicated to select landmarks in a semi-supervised DA setting. Lastly, the single scale approach (using a fixed σ) does not perfom well: this illustrates how important it is to select the best radius of action for every landmark in MLS. Figure 5 gives the distribution of the landmarks selected by MLS for each DA subproblem, showing that even without class information, our approach makes a balanced selection among the classes. Comparison with the state of the art unsupervised DA approaches Table 2 reports the results of the experimental comparison between the state of the art unsupervised subspace alignment DA methods. is worth noticing that our method LSSA outperforms the approaches in 7 out of 12 DA subproblems, while TJM is better for the remaining 5 subproblems. However, on average, LSSA significantly outperforms TJM (52.6% versus 50.5%). Moreover, as pointed out by the authors, the time complexity of TJM is larger than the other approaches because it requires to solve a non trivial optimization problem, while our approach which involves a greedy strategy (for the landmark selection) as well as a closed solution (for the subspace alignment) is more efficient. We claim that the difference of accuracy between TJM and LSSA comes from the fact that the former uses a weighting scheme that mainly aims at moving closer the means of the two domains, while the latter takes into account, via the Gaussian hypothesis, not only the mean but also the standard deviation of the statistical distribution of the landmarks. From Table 2, we can also notice that LSSA significantly outperforms SA, meaning that capturing non linearity in LSSA is a big improvement over plain SA. However, the way to consider the non linearity is also key. Indeed, as shown by KPCA+SA, performing two independent KPCA before the subspace alignment leads to the worst behavior.

Conclusion

Subspace alignment-based DA methods have recently attracted a lot of interest. In this framework, the most accurate approach assumes that the shift between the source and target distributions can be corrected by a linear function. However, we experimentally show that this assumption is challenged in most of the real world applications. Thus, we argue that a non linear mapping function should be optimize to align the source and target subspaces. Furthermore, there is no reason that would justify to constrain a DA algorithm to adapt all the training source data to the target domain. In this paper, we have combined both ideas in a new unsupervised DA algorithm (LSSA) which first selects landmarks that allow us to create, by a non linear projection, a common space where the two domains are closer to each other and then performs a subspace alignment. LSSA is simple, fast and outperforms the state of the art on a visual domain adaptation task. As a future work, we plan to resort to an optimization process that jointly selects the landmarks in order to avoid redundancy between them as it might appear with the current greedy and independent selection.

Figure 2 .

 2 Figure 2. (Best shown in color) Domain adaptation problem that needs to adapt a source (S) domain to a target (T) domain w.r.t. four classes of examples (green, orange, blue and red).We can see that no linear transformation will be able to entirely bridge the source to the target. On the other hand, we can note that, for each class, some subsets of source and target examples seem to be similarly distributed (in black). The objective is to automatically discover those so-called landmarks that will be used to map both source and target data onto a common space using a kernel.

Figure 3 .

 3 Figure 3. Intuition behind the idea of learning intermediate subspaces along the geodesic path connecting the source and target subspaces.

Figure 5 .

 5 Figure 5. Distribution of the landmarks among the ten classes for each DA subproblem.

BackPack, Bike, Calculator, Headphone, Keyboard, Laptop, Monitor, Mouse, Mug, Projector.

Acknowledgement This work is supported by the ANR SOLSTICE project (ANR-13-BS02-01).

Algorithm 2 LSSA: Landmarks Selection-based Subspace Alignment and classification. Require: Source data S, Target data T , Source labels L S , Threshold th, Subspace dimension d. Ensure: L T are the predicted labels for the points in T.

A ← choose landmarks(S, T, th) σ ← median distance(S ∪ T ) K S ← project using kernel(S, A, σ) K T ← project using kernel(T, A, σ)

the art, including the two best recent methods presented in [7,[START_REF] Long | Transfer joint matching for unsupervised domain adaptation[END_REF].

Datasets and Experimental Setup

We run our experiments on standard datasets for visual domain adaptation. We use the Office dataset [START_REF] Saenko | Adapting visual category models to new domains[END_REF] that contains images acquired with webcams (denoted by W), images taken with digital SLR camera (denoted by D) and Amazon images (denoted by A). In addition, we use Cal-tech10 [START_REF] Gong | Geodesic flow kernel for unsupervised domain adaptation[END_REF], denoted by C. Each dataset provides different images from 10 classes 1 . Therefore, we can generate 12 domain adaptation subproblems from the four datasets (A, C, D, W), each time one dataset playing the role of the source S while another one is considered as the target T . We denote a DA problem by the notation S → T . The objective is to learn an SVM classifier (with a linear kernel using the SV M light implementation) from the labeled source S and to deploy it on the target T . We use the same image representation as that provided by [START_REF] Gong | Geodesic flow kernel for unsupervised domain adaptation[END_REF] for Office and Cal-tech10 datasets (SURF features encoded with a visual dictionary of 800 words). We follow the standard protocol of