Nicolas Duchateau

Maxime Sermesant

Hervé Delingette

Nicholas Ayache

Model-based generation of large databases of cardiac images: synthesis of pathological cine MR sequences from real healthy cases

Keywords: Image synthesis and simulation, large database for machine learning, validation, myocardial infarct, biomechanical modeling, magnetic resonance

Collecting large databases of annotated medical images is crucial for the validation and testing of feature extraction, statistical analysis and machine learning algorithms. Recent advances in cardiac electromechanical modeling and image synthesis provided a framework to generate synthetic images based on realistic mesh simulations. Nonetheless, their potential to augment an existing database with large amounts of synthetic cases requires further investigation. We build upon these works and propose a revised scheme for synthesizing pathological cardiac sequences from real healthy sequences. Our new pipeline notably involves a much easier registration problem to reduce potential artifacts, and takes advantage of mesh correspondences to generate new data from a given case without additional registration. The output sequences are thoroughly examined in terms of quality and usability on a given application: the assessment of myocardial viability, via the generation of 465 synthetic cine MR sequences (15 healthy and 450 with pathological tissue viability [random location, extent and grade, up to myocardial infarct]). We demonstrate that our methodology (i) improves state-of-the-art algorithms in terms of realism and accuracy of the simulated images, and (ii) is well-suited for the generation of large databases at small computational cost.

I. INTRODUCTION

Despite decades of designing algorithms for the analysis of medical images, the collection of proper databases of annotated medical images remains crucial to validate or just correctly use such algorithms.

On the one hand, standard feature extraction methods from single images/sequences such as segmentation and registration require ground truth imaging data for validation purposes. This validation relies on manual annotations by experts, invasive tracers attached to the anatomy, or phantom data. However, these processes are cumbersome and have limited accuracy or realism. New types of ground truth data are therefore required.

On the other hand, data-based statistical analysis and machine learning algorithms are optimally exploited with large sample size. Computer vision research now provides access to datasets of 100,000+ 2D images 1 . In contrast, collecting ND, MS, HD and NA, are with the INRIA Asclepios team, Sophia Antipolis, France. ND is also with the CREATIS team, CNRS UMR 5220, Lyon, France. Address for correspondence: Nicolas Duchateau, INRIA Asclepios, 2004 route des Lucioles BP 93 06902 Sophia Antipolis Cedex, France. Tel: +33.49238.5024; Fax +33.49238.7669.

Email: nicolas.duchateau@creatis.insa-lyon.fr ; nicholas.ayache@inria.fr 1 https://en.wikipedia.org/wiki/List of datasets for machine learning research relevant sets of medical images undergoes limitations from healthcare data management (permission and storage size) [START_REF] Suinesiaputra | Big heart data: advancing health informatics through data sharing in cardiovascular imaging[END_REF] and acquisition (image quality and availability) [START_REF] Klinke | Quality assessment of cardiovascular magnetic resonance in the setting of the European CMR registry: description and validation of standardized criteria[END_REF]. In our field of interest -cardiac imaging, and especially motion and deformation analysis-current datasets with 3D sequences and coherent pathophysiological condition comprise around 1000 subjects [START_REF] Fonseca | The Cardiac Atlas Project-an imaging database for computational modeling and statistical atlases of the heart[END_REF]- [START_REF] Bai | A bi-ventricular cardiac atlas built from 1000+ high resolution MR images of healthy subjects and an analysis of shape and motion[END_REF]. The recent UK Biobank project intends to gather an unprecedented database of 100,000 cardiac sequences [START_REF] Petersen | Imaging in population science: cardiovascular magnetic resonance in 100,000 participants of UK Biobank -rationale, challenges and approaches[END_REF]. Nonetheless, it targets an outpatient or normal population, which means that the study of specific pathologies in coherent subgroups might still have limited sample size or be subject to a large variety of confounding factors [START_REF] Ganna | 5 year mortality predictors in 498,103 UK Biobank participants: a prospective population-based study[END_REF].

Thus, it is urgent to develop alternative approaches to collect (i) large datasets with (ii) coherent clinical conditions and (iii) ground truth data. For this, we propose to build upon recent initiatives for simulating and synthesizing realistic medical images [START_REF] Glatard | A virtual imaging platform for multi-modality medical image simulation[END_REF], [START_REF] Tsaftaris | Simulation and Synthesis in Medical Imaging[END_REF], and to extend the existing literature by specifically investigating database generation issues for a given clinical problem -in our case, localizing myocardial infarct.

A. Methods for simulating realistic cardiac images

Unlike modality propagation approaches, which transfer anatomical information from one modality to another [START_REF] Ye | Modality propagation: coherent synthesis of subject-specific scans with datadriven regularization[END_REF], [START_REF] Cordier | Extended modality propagation: image synthesis of pathological cases[END_REF], most cardiac synthesis methods include the creation of new data -new myocardial shape or motion, obtained from models of cardiac function. In the quest for higher realism, state-of-the-art methods base this anatomical and functional knowledge on individual patients' segmentations [START_REF] De Craene | 3D strain assessment in ultrasound (Straus): a synthetic comparison of five tracking methodologies[END_REF]- [START_REF] Zhou | Generation of realistic 4D synthetic CSPAMM tagged MR sequence for benchmarking cardiac motion tracking algorithms[END_REF], and simulated motion from realistic electromechanical models of the heart [START_REF] De Craene | 3D strain assessment in ultrasound (Straus): a synthetic comparison of five tracking methodologies[END_REF]- [START_REF] Duan | Simulation of 3D ultrasound with a realistic electro-mechanical model of the heart[END_REF].

In terms of image synthesis, a first branch of methods investigated the physics of imaging in echocardiography [START_REF] Gao | A fast convolutionbased methodology to simulate 2-D/3-D cardiac ultrasound images[END_REF], [START_REF] Shakil | Simulation in echocardiography: an ever-expanding frontier[END_REF], computed tomography [START_REF] Segars | Realistic CT simulation using the 4D XCAT phantom[END_REF] and magnetic resonance [START_REF] Kwan | MRI simulation-based evaluation of image-processing and classification methods[END_REF], [START_REF] Tobon-Gomez | Realistic simulation of cardiac magnetic resonance studies modeling anatomical variability, trabeculae, and papillary muscles[END_REF]. However, the simulated images were still far from real data as texture, artifacts, and surrounding structures were not taken into account. A significant breakthrough came from the incorporation of real image properties into the synthesized images. Highly realistic results were obtained for 3D echocardiography, with scatterers distributed as in the real sequence and displaced according to the electromechanical simulation [START_REF] Alessandrini | A pipeline for the generation of realistic 3D synthetic echocardiographic sequences: methodology and open-access database[END_REF]. This scheme ideally combines the physics of image formation and real data, but is modality-specific. Besides, generating a new sequence with different myocardial motion requires simulating again the image formation, which may be costly for the creation of large synthetic databases.

In parallel, a second branch of methods investigated the warping of existing images according to a given motion model. First attempts on real data warped a single static CT image with motion from a registered 3D echocardiographic sequence [START_REF] Li | Generation of synthetic 4D cardiac CT images by deformation from cardiac ultrasound[END_REF], and a single static MR image with a simple kinematic model [START_REF] Clarysse | Simulation based evaluation of cardiac motion estimation methods in tagged-MR image sequences[END_REF]. Substantial improvements were obtained from the use of a dynamic sequence, previously stabilized to a given instant by registration along the sequence, and realistic simulated motion (Fig. 1) [START_REF] Prakosa | Generation of synthetic but visually realistic time series of cardiac images combining a biophysical model and clinical images[END_REF]. Although this approach can be applied to any modality, it is better suited for cine MRI and CT. Indeed, undesirable warping artifacts may be observed for highly textured data as in echocardiography [START_REF] Prakosa | Generation of synthetic but visually realistic time series of cardiac images combining a biophysical model and clinical images[END_REF] and tagged-MRI [START_REF] Zhou | Generation of realistic 4D synthetic CSPAMM tagged MR sequence for benchmarking cardiac motion tracking algorithms[END_REF].

In the perspective of generating large databases of synthetic cases, we would like to build upon this approach for its portability, provided that specific aspects are improved. In particular, artifacts were reported in [START_REF] Prakosa | Generation of synthetic but visually realistic time series of cardiac images combining a biophysical model and clinical images[END_REF] when the simulated motion is too different from the original motion. Besides, inaccuracies in the registration along the sequence may propagate and alter the realism of the outcome, and bias its correspondence with the simulated motion used as ground truth. Thus, part of this paper investigates two substantial improvements to this scheme: updating the processing pipeline and therefore solving a much easier registration problem, and directly re-using mesh correspondences to generate new data from a given case without additional registration.

B. Assessment of tissue viability and infarct diagnosis

In terms of application, our objective is to evaluate the potential of this approach for computer-aided diagnosis algorithms, through the generation of large synthetic databases containing ground truth shape and motion. We would like to focus on a concrete clinical problem under active investigation: the assessment of myocardial viability from cardiac image sequences.

Coronary heart disease is actually the leading cause of death in occidental countries (1.8 million deaths per year in Europe [START_REF] Nichols | European cardiovascular disease statistics[END_REF]). Imaging data is crucial for clinical routine: it facilitates the reading of cardiac function up to shape, motion and deformation descriptors, which are directly related to changes in the myocardial viability [START_REF] Bax | Detection of viable myocardium and scar tissue[END_REF]- [START_REF] Konstam | Left ventricular remodeling in heart failure: current concepts in clinical significance and assessment[END_REF]. These descriptors are complex and in high dimension, and therefore require powerful computational tools that properly integrate them to the diagnosis process. Recent works using machine learning showed promising results for infarct diagnosis from image [START_REF] Sudarshan | Automated identification of infarcted myocardium tissue characterisation using ultrasound images: a review[END_REF] or shape data [START_REF] Suinesiaputra | Statistical shape modelling of the left ventricle: myocardial infarct classification challenge[END_REF]. In the meantime, we presented a learning method that can estimate the infarct localization [START_REF] Duchateau | Infarct localization from myocardial deformation: prediction and uncertainty quantification by regression from a low-dimensional space[END_REF]. Validation was performed on a large set of synthetic data, which consisted of meshes from realistic electromechanical simulations -but not images. Focus was on the accuracy of the localization, and the database therefore included infarcts of varying extent, location and shape.

In the present work, we would like to move forward and design a database of synthetic images. Besides, we propose to challenge the limits of diagnosis by refining the database design and incorporating different grades of tissue viability. Such cases can be obtained by altering parameters of the ... Processing pipeline used in [START_REF] Prakosa | Generation of synthetic but visually realistic time series of cardiac images combining a biophysical model and clinical images[END_REF]. The original cardiac sequence is "freezed" to the first instant via image registration (top row). Better performance is reached when initializing the registration between frames 0 and t by the composition of transformations first obtained from the registration of consecutive frames t and t + 1 (in practice, handled in the tangent space as the sum of velocity fields linking frames t and t + 1). Then, this "static" image is warped according to the simulation displacement field (bottom row). Consistency between the myocardium and the inside/outside anatomy is ensured by diffusion on this displacement field.

electromechanical simulation from healthy to locally diseased. This strategy inherits from [START_REF] Prakosa | Cardiac electrophysiological activation pattern estimation from images using a patient-specific database of synthetic image sequences[END_REF], which varied the electrical activation parameters to generate up to 180 scenarios from three given pathological image sequences. However, this design modified the parameters around existing personalized values, but did not investigate the creation of pathological cases from healthy data, and therefore the limits of detectable disease and image realism. We tackle this issue in the present paper.

C. Proposed method and contributions

We propose a revised scheme for the generation of realistic 3D cine MR cardiac sequences, which makes substantial improvements to the current state-of-the-art outlined in [START_REF] Prakosa | Generation of synthetic but visually realistic time series of cardiac images combining a biophysical model and clinical images[END_REF], [START_REF] Alessandrini | A pipeline for the generation of realistic 3D synthetic echocardiographic sequences: methodology and open-access database[END_REF]. Our contributions are both on technical and applicative aspects.

First, our revised processing only involves the registration of rather similar image sequences. This easier registration problem reduces the amount of potential artifacts. In contrast, registration along the sequence involves larger transformations and potential artifacts that could propagate along the sequence [START_REF] Prakosa | Generation of synthetic but visually realistic time series of cardiac images combining a biophysical model and clinical images[END_REF], [START_REF] Zhou | Generation of realistic 4D synthetic CSPAMM tagged MR sequence for benchmarking cardiac motion tracking algorithms[END_REF]. Besides, this strategy biases the validation of the synthesis methodology, if evaluated via tracking along the sequence by a similar algorithm [START_REF] Prakosa | Generation of synthetic but visually realistic time series of cardiac images combining a biophysical model and clinical images[END_REF], [START_REF] Alessandrini | A pipeline for the generation of realistic 3D synthetic echocardiographic sequences: methodology and open-access database[END_REF].

Then, our new pipeline takes advantage of mesh correspondences to generate variants of a given simulation. No additional registration is involved and correspondences are exact. This opens the door to the generation of large databases of synthetic cases at small computational cost, demonstrated on a given application: the assessment of myocardial viability.

Finally, we propose a thorough qualitative and quantitative evaluation of the quality and usability of the generated images, in the context of the generation of large databases. We notably investigate whether the designed database is balanced below and above the diagnosis limit of a given algorithm, within the limits of detectable disease and image realism, as illustrated in Fig. 2.

II. MATERIALS AND METHODS

A. Template sequences: real database from STACOM 2011

We used the cine MR sequences from 15 volunteers collected for the 1st Cardiac Motion Analysis (cMAC) challenge held at STACOM 2011. These volunteers had no clinical history of cardiac diseases (age 28 ± 5 years, 12(80%) male). Details about data acquisition and the subjects' characteristics can be found in the challenge paper [START_REF] Tobon-Gomez | Benchmarking framework for myocardial tracking and deformation algorithms: an open access database[END_REF]. Images had standard quality, with a spatial resolution of 1.15 or 1.25 mm in the x and y directions, and 8 mm in the z direction (11 to 16 slices, average 14 slices), and 30 frames along the cycle.

Data was already segmented at end-diastole by the semiautomatic protocol described in [START_REF] Marchesseau | Personalization of a cardiac electromechanical model using reduced order unscented Kalman filtering from regional volumes[END_REF], which involved a biventricular deformable model and manual corrections. As illustrated in Fig. 3, the available meshes were volumetric with tetrahedral elements, and included the two ventricles (23770± 5219 tetrahedra and 6633±1212 vertices for the LV, of average edge length 3.6 ± 0.2 mm, for a myocardial volume and mass of 95 ± 27 mL and 101 ± 28 g, respectively). Rule-based fibre directions were added over this geometry with an elevation angle between 70 and 70 .

B. Mesh simulations

We applied a realistic electromechanical model to each of the 15 cases to simulate the cardiac function. The model computations were updated each 2 ms of a full cardiac cycle of 1 s. Volumetric meshes were exported at uniformly sampled instants of the cycle (each 17 steps, resulting in 30 meshes along the cycle).

The model was implemented within the SOFA framework 2 , and derives from the Bestel-Clement-Sorine model for the active contraction, and the Mooney-Rivlin model for passive elasticity. It consists of energy-preserving equations during contraction and relaxation, and includes the Frank-Starling effect. We refer the reader to previous publications from our group for further details about this model, including validation issues and a full sensitivity analysis of its parameters [START_REF] Marchesseau | Personalization of a cardiac electromechanical model using reduced order unscented Kalman filtering from regional volumes[END_REF], [START_REF] Marchesseau | Preliminary specificity study of the Bestel-Clement-Sorine electromechanical model of the heart using parameter calibration from medical images[END_REF]. It was retained for its ability to reproduce the behavior of real cases, as evaluated in [START_REF] Marchesseau | Personalization of a cardiac electromechanical model using reduced order unscented Kalman filtering from regional volumes[END_REF] and more recently in [START_REF] Mollero | A multiscale cardiac model for fast personalisation and exploitation[END_REF], [START_REF] Giffard-Roisin | Non-invasive personalisation of a cardiac electrophysiology model from body surface potential mapping[END_REF]. The realism of simulated pathological behaviors was also tested through the generation of a large database of meshes with spatially varying infarcts [START_REF] Duchateau | Infarct localization from myocardial deformation: prediction and uncertainty quantification by regression from a low-dimensional space[END_REF], and single infarcted or dyssynchronous cases used as input to image synthesis algorithms [START_REF] Prakosa | Generation of synthetic but visually realistic time series of cardiac images combining a biophysical model and clinical images[END_REF], [START_REF] Alessandrini | A pipeline for the generation of realistic 3D synthetic echocardiographic sequences: methodology and open-access database[END_REF].

1) Healthy cases: The simulation parameters were adapted from the ranges in [START_REF] Alessandrini | A pipeline for the generation of realistic 3D synthetic echocardiographic sequences: methodology and open-access database[END_REF] and [START_REF] Duchateau | Infarct localization from myocardial deformation: prediction and uncertainty quantification by regression from a low-dimensional space[END_REF], and are detailed in Tab.I. Balance between the right and left ventricles was slightly updated to reach a more synchronous relaxation on these new geometries. The conductivity and the Rayleigh stiffness were respectively set to 500 mm • s 1 and 0.06 for the same reason (400 mm • s 1 and 0.1 in [START_REF] Duchateau | Infarct localization from myocardial deformation: prediction and uncertainty quantification by regression from a low-dimensional space[END_REF]). Values were also adapted over the mitral annulus to reduce its deformation and better fit the imaging data. Case-based personalization [START_REF] Mollero | A multiscale cardiac model for fast personalisation and exploitation[END_REF], [START_REF] Neumann | A self-taught artificial agent for multiphysics computational model personalization[END_REF] was not targeted here, and the same parameters were used for all cases. Nonetheless, the retained values led to realistic behaviors close to the original sequences. The personalization error was evaluated in terms of volumes and distance between surfaces. At end-systole, the distance between the manually-segmented and simulated LV endocardial surfaces 3 was 3.6 ± 1.7 mm. The agreement between the manuallysegmented and simulated LV volumes is plotted in Fig. 4a, and a representative example is represented in Fig. 5.

2) Diseased cases: Variants of each case in the database were generated by modifying the simulation parameters in a given zone, labeled as diseased.

This zone was randomly located within the LV territory irrigated by the left-anterior-descending coronary -territory of higher infarct prevalence [START_REF] Ortiz-Perez | Correspondence between the 17-segment model and coronary arterial anatomy using contrast-enhanced cardiac magnetic resonance imaging[END_REF]-and had random extent and shape. It was defined iteratively by an algorithm similar to the one used in [START_REF] Duchateau | Infarct localization from myocardial deformation: prediction and uncertainty quantification by regression from a low-dimensional space[END_REF], which was improved on two aspects: (i) the algorithm now handles mesh labels as cell attributes instead of vertex attributes for better compatibility with the SOFA mesh format; (ii) for better realism with the mechanisms of ischemia, the zone generation is now constrained to start from the endocardial surface. At each iteration, a cell is randomly selected within the region already labeled as diseased, and a spherical neighborhood of random radius between 5 and 15 mm is picked out around this point and added to the diseased region. The algorithm repeats this process a random number of times (values from 1 to 16). On our data, this led to diseased zones of 5.8 ± 3.0 mL, corresponding to 6.3 ± 3.4% of the LV myocardium.

Then, the simulation parameters with major influence on the deformation [START_REF] Marchesseau | Preliminary specificity study of the Bestel-Clement-Sorine electromechanical model of the heart using parameter calibration from medical images[END_REF] were altered in this diseased region, to mimic changes in active force and tissue elasticity. Contrary to our previous work [START_REF] Duchateau | Infarct localization from myocardial deformation: prediction and uncertainty quantification by regression from a low-dimensional space[END_REF], values were randomly modified within a given interval to simulate varied grades of tissue viability. This was achieved by setting all the selected parameters at a random portion of the healthy-infarct range in Tab.I, 100% standing for healthy and 0% for the infarcts modeled here. This range was adapted from previous publications [START_REF] Alessandrini | A pipeline for the generation of realistic 3D synthetic echocardiographic sequences: methodology and open-access database[END_REF], [START_REF] Duchateau | Infarct localization from myocardial deformation: prediction and uncertainty quantification by regression from a low-dimensional space[END_REF] through new tests on the cMAC-STACOM 2011 data, aiming at realistic motion close from the imaging data. This approach might not correspond to a fine biophysical modeling of the 3 Manual segmentation was performed on 2D slices through the Segment software (Medviso, Lund, SE). Intra-operator reproducibility was 0.58 and 0.57 mm at end-diastole and end-systole, respectively. We used a nonsymmetric distance between contours, as we compare a 3D surface ⌦ simu (the endocardium from the mesh simulation) to a stack of 2D contours (the segmented endocardium in each slice). The distance was defined as

1 K P K k=1 d(x k , ⌦ simu )
, where x k denotes the k-th point of the segmentation stack. Thus, these surfaces did not perfectly match at end-diastole: their distance was 2.3 ± 0.4 mm, and the volume difference was 11.3 ± 6.6 mL. infarct, but this would require a more in-depth study, outside the scope of this paper.

We created 30 realistic pathological variants of each of the 15 cases in the database, resulting in 450 diseased cases that complement the 15 healthy simulations described in Sec.II-B1. The characteristics of all the generated cases (extent, location, and tissue viability) are summarized in Fig. 6.

C. Image synthesis 1) Generation of a reference synthetic sequence: This part of the methodology is synthesized in Algorithm 1. It corresponds to the generation of a synthetic sequence close from the original one (the first two rows in Fig. 7) -in our case, a healthy sequence. The original sequence is cropped around the heart and upsampled in the z direction to reach an isotropic resolution (step 1, linear interpolation via the ITK image resampling filter). The first frame in the sequence is then warped according to the healthy mesh simulation (steps 2 to 4). In practice, the displacement field ' 0 (t) from the simulation between times 0 and t is only available inside the myocardium. Diffusion is therefore applied to expand this displacement field before applying it to the frame at time 0 (Fig. 8a-b). This guarantees continuity in the motion of the anatomical structures surrounding the heart, with decreasing effect according to the distance from the myocardium. As in [START_REF] Prakosa | Generation of synthetic but visually realistic time series of cardiac images combining a biophysical model and clinical images[END_REF], this was achieved via the iterative convolution with a Gaussian kernel (the ITK recursive Gaussian filter with a standard deviation of 1 voxel) followed by the reinitialization of the vector field values over the myocardium and forcing them to 0 at the image borders.

The process was iterated at most 500 times, and stopped if the difference in vectors between two iterations is lower than 0.001 voxel. Despite modifying the parameters value at the mitral annulus (Sec.II-B1), the mesh simulations may still differ from the original sequence at this level. Thus, the displacement field from the simulation is also cropped 1 cm below the mitral annulus before the diffusion step (Fig. 8c). This arbitrary cropping limit roughly corresponds to the myocardial thickness, and our tests showed that it improves the realism of the simulated images. This process provides a sequence that moves as the simulated motion, but keeps the texture of the first frame, which is still unrealistic [START_REF] Li | Generation of synthetic 4D cardiac CT images by deformation from cardiac ultrasound[END_REF], [START_REF] Clarysse | Simulation based evaluation of cardiac motion estimation methods in tagged-MR image sequences[END_REF]. Nonetheless, if the simulated motion is similar to the one in the original sequence, registering the original sequence to this warped sequence involves small displacements and should therefore be relatively easy to solve (step 5, also marked with red color in Fig. 7). We achieved this by applying the Demons algorithm with sum-of-squared-Algorithm 1: Generation of a healthy sequence. To prevent artifacts outside the heart, the vector field original!warped (t) obtained from this registration is multiplied by a diffused mask of the myocardium. This mask corresponds to a binary image of the myocardium that is expanded outside the myocardium by a diffusion process similar to the one used in Sec.II-C1. The masked vector field is then smoothed temporally to force consistency in the temporal motion (step 6, via the ITK recursive Gaussian filter applied along the temporal dimension of the 4D vector field, with a standard deviation of 2 t , where t is arbitrarily set to 1), and applied to the original sequence (step 7). Finally, the simulated sequence is downsampled to the resolution of the original sequence (step 8, ITK filter similar to the one in step 1).

2) Generation of additional image sequences from the reference synthetic sequence: Given a simulated sequence close from the original data (in our case, the healthy sequence generated in Sec.II-C1 from the healthy mesh simulation), we propose to exploit mesh correspondences to generate new sequences corresponding to the infarct configurations described in Sec.II-B2. This part of the methodology corresponds to the blue-squared part of Fig. 7, and is synthesized in Algorithm 2.

For a given infarct configuration, we first extract the vector field healthy!inf arct (t) describing the mesh correspondences at a given instant, and expand it outside the myocardium via the diffusion process of Sec.II-C1. This vector field is then composed with ⇤ original!warped (t), the one used to warp the original sequence in the healthy case, and finally applied to the original sequence.

3) Parallelization of the processes and computational time: All the mesh simulations were performed in parallel on our local cluster platform 5 . Each mesh simulation from SOFA 4 Chosen parameters: smoothing of the displacement field by 3 voxels, smoothing of the update field by 0.5 voxel, weight between the similarity and regularization terms of 0.1, 3 levels of resolution (30, 20 and 10 iterations, respectively), Gaussian smoothing of the similarity criteria by 3 voxels, and Gaussian regularization. Code v1.2, April 2016 from [START_REF] Lorenzi | LCC-Demons: a robust and accurate symmetric diffeomorphic registration algorithm[END_REF], available at: https://team.inria.fr/asclepios/software/lcclogdemons/ 5 https://wiki.inria.fr/ClustersSophia/Hardware Algorithm 2: Generation of a diseased sequence.

for i 2 N inf arct do for t 2 [0, T ] 1. Extract mesh correspondences Output: vector field healthy!inf arct (t)

2. Diffusion: D( healthy!inf arct (t))

Compose

D( healthy!inf arct (t))

⇤ original!warped (t) 4. Warp original image at time t by this field 5. Downsample sequence to the original resolution takes approximatively one hour to complete. Each registration takes at most 30 min. The diffusion is the longest process in the pipeline, and takes around 30 min. The remaining processes take less than 1 min.

For the generation of the reference sequence, all the instants can be processed independently and therefore parallelized except the temporal smoothing, which involves the whole 4D sequence. A similar parallelization was adopted for the generation of the variants of the reference sequence. Thus, simulating one or several variants of the same reference sequence takes the same amount of time -within the limits of parallel processes allowed for a given user on our cluster platformand is therefore a substantial asset for the generation of large databases of synthetic sequences.

III. EXPERIMENTS AND RESULTS

Supplementary Material6 is provided for a better assessment of the synthetic sequences quality. It includes animated views of Fig. 5 and9, and animated views of the synthetic sequences for the best, medium and worst personalized cases for the healthy and strong infarct (0% viability) simulations, which complement the illustration in Fig. 11.

A. Qualitative evaluation 1) Mimicking an ideal personalization: We first designed a toy experiment to highlight the limits of the method from [START_REF] Prakosa | Generation of synthetic but visually realistic time series of cardiac images combining a biophysical model and clinical images[END_REF]. The first frame in the sequence was warped along the cycle, as in step 4 of Algorithm 1 and was used as input to the image synthesis algorithm illustrated in Fig. 1. By definition, this toy sequence matches the mesh simulation. It therefore mimics an ideal personalization, namely a perfect match between the original sequence and the mesh simulation (Fig. 9a-b). The method from [START_REF] Prakosa | Generation of synthetic but visually realistic time series of cardiac images combining a biophysical model and clinical images[END_REF] fails at reproducing the toy sequence (Fig. 9c-f), due to the mismatch between the displacement field from the sequence registration and the one from the mesh simulation. In contrast, our method does not suffer these limitations on this toy experiment: the registration part of our method involves already matched images, as the warped sequence is used in lieu of the original sequence, and therefore original!warped (t) = 0.

2) Matching with the mesh contours: Figure 10 illustrates the representative differences between the output of our method and the one from [START_REF] Prakosa | Generation of synthetic but visually realistic time series of cardiac images combining a biophysical model and clinical images[END_REF], on the best and worst personalized cases (similar observations were made on all cases). Our method better matches the mesh contours. In contrast, for the method from [START_REF] Prakosa | Generation of synthetic but visually realistic time series of cardiac images combining a biophysical model and clinical images[END_REF] the myocardium spreads inside the LV and RV mesh blood pool, stretching is higher near the epicardium, and other anatomical zones such as the mitral annulus are insufficiently deformed. Figure 11 complements these observations by showing that the matching with the mesh contours is conserved when simulating strongly diseased tissue, despite substantial changes in the simulation (contraction is lower and mesh contours are thinner at the infarct level). The animated views [START_REF] Prakosa | Generation of synthetic but visually realistic time series of cardiac images combining a biophysical model and clinical images[END_REF] and our approach on a toy example (subject v7), demonstrating the limits of [START_REF] Prakosa | Generation of synthetic but visually realistic time series of cardiac images combining a biophysical model and clinical images[END_REF] by the generated MR image at end-systole. (a) and (b): first frame warped along the sequence (step 4 of Algorithm 1), referred to as the toy sequence. By definition, it matches the mesh simulation and mimics an ideal personalization. Thus, it also corresponds to the result of our method on this ideal case, as the warped sequence is used in lieu of the original sequence and therefore original!warped (t) = 0. (c) and (d): sequence simulated from the toy sequence using the method from [START_REF] Prakosa | Generation of synthetic but visually realistic time series of cardiac images combining a biophysical model and clinical images[END_REF], which fails at reproducing the toy sequence. (e) and (f): signed difference between the toy sequence and the result of [START_REF] Prakosa | Generation of synthetic but visually realistic time series of cardiac images combining a biophysical model and clinical images[END_REF]. An animated version of the figure is available as Supplementary Material 6 . in the Supplementary Material 6 lead to similar conclusions on the best, medium and worst personalized cases. Slight mismatch is observed along the lateral wall of the medium personalized case. This results from mis-aligned slices in the original sequence, at the infarct level where images are shown, which are not taken into account in the 3D mesh used for the simulations.

3) Examining the limits to the image realism: The warping applied to the original sequence might alter the realism of the generated synthetic images, although this depends on many factors, among which the quality of the personalization for the healthy simulation, the characteristics of the local disease, and the texture of the original sequence in the warped regions (myocardium and surroundings). Figure 13 provides a qualitative assessment of the simulated images with the highest tissue damage (8 cases with tissue viability < 10% and diseased zone volume > 10%), against the original and healthy simulated images. Warping artifacts are hardly noticeable for the healthy simulations, although changes in the myocardial thickness and shape can be noted at the apical level of some cases (v2 infarct #7 and v5 infarct #23). Due to the diffusion of the vector fields around the myocardium, the image realism is preserved in the myocardial surroundings, even for the cases with low

Our method

Prakosa et al.

Our method

Worst personalized case Fig. 10. Qualitative comparison of our method to the one from [START_REF] Prakosa | Generation of synthetic but visually realistic time series of cardiac images combining a biophysical model and clinical images[END_REF], on the best and worst personalized cases (v6 and v12, average personalization errors at end-systole: 2.1 and 8.3 mm, respectively). Our method better fits the mesh simulation, as highlighted on the blood pool and epicardium along the LV lateral wall, on the RV contours, and at the mitral annulus (white arrows). Similar observations were made on all cases. tissue viability. The main limits to the image realism actually come from the motion prescribed in the mesh simulations, with sometimes noticeable myocardial shape differences due to highly contrasted myocardial border in the original sequence (v4 infarct #2 and #18, v5 infarct #25 and v15 infarct #22).

B. Quantitative evaluation 1) Manual assessment: segmentation: We manually segmented a representative subset of the simulated images at end-systole, made of the healthy one and the ones with tissue viability closest to 50 and 0%, for each of the 15 cases in the database. Then, we compared the segmentations (stacked 2D contours) with the mesh simulations (3D surfaces, taken as ground truth). These two steps followed the protocol adopted for the original sequence, as described in Sec.II-B1. Agreement is observed on volumes (Fig. 4b), and reflects the accuracy of the synthesized images with respect to the prescribed motion. This agreement is naturally better than the match between the original sequence and the simulated end-diastole end-diastole contours (Fig. 4a), which reflected the personalization error. Similar observations can be made in terms of distance between the surfaces 3 , without substantial differences between the different grades of tissue viability: values were 2.1 ± 0.4, 2.2 ± 0.5 and 2.3 ± 0.5 mm for the healthy simulation and the tissue grades of 50 and 0%, respectively. These measurements lie in the range of differences between the 3D surface and the stack of 2D contours at end-diastole (Sec.II-B1). One common outlier was observed, which actually corresponded to a case in the cMAC-STACOM database (V11) with two mis-aligned slices at the papillary muscles level, while the 3D mesh used for the simulations does not take this into account.

2) Automated assessment: registration: We finally evaluated the sequence quality against state-of-the-art motion estimation along the sequence. We applied the Demons algorithm with sum-of-squared-differences similarity metric. We followed a two-step strategy similar to the one in Fig. 1, where pairs of consecutive frames are first registered, and the resulting transformations are composed (in practice, handled in the tangent space as the sum of velocity fields linking frames t and t + 1) and used as initialization for a second registration between frames at t = 0 and t. Parameters 4 were similar to the ones detailed in Sec.II-C1, except the smoothing of the displacement field -by 2 voxels here. Then, the first frame was propagated along the sequence according to the estimated displacement field.

The simulated contraction deficiency due to the altered tissue viability was assessed through the LV regional volume change along the sequence, on the simulated meshes and the ones propagated along the simulated images via the abovedescribed registration. The regional volume of a given region is defined as the volume formed by the LV barycenter and the endocardial surface of this region [START_REF] Marchesseau | Personalization of a cardiac electromechanical model using reduced order unscented Kalman filtering from regional volumes[END_REF]. This measure partially relates to the myocardial deformation, and is more robust to small registration errors since it averages values over a given zone, and does not involve a spatial derivative.

These measurements are plotted in Fig. 12 against the tissue viability. Differences between the healthy and diseased behaviors get more contrasted as the tissue viability decreases. They become significant for medium-grade tissue viability, although results are less contrasted for the simulated images. The latter mainly comes from the information loss when computing these values on meshes (high resolution) or images (lower resolution, in particular on the z direction), and from the accuracy of the motion estimation along the sequence: lower regional volume change is also observed for the healthy myocardium; this should not be the case, as segmentation showed that the simulated sequences agree with the prescribed motion (Sec.III-B1).

IV. DISCUSSION

We presented a new approach to synthesize cardiac image sequences from existing real sequences and a realistic electromechanical model, demonstrated on cine MR data. Our pipeline improves state-of-the-art methods by revisiting specific aspects of the processing, resulting in higher realism of the simulated images and better suitability for the generation of large databases of pathological cases. In this context, we designed a large database of 15 healthy and 450 diseased cases, and evaluated the quality and usability of the generated images.

Our method provides higher image realism and better agreement with the ground truth motion (Fig. 10). This mainly results from avoiding registration along the original sequence in our pipeline. Such a registration is a potential source of artifacts, as demonstrated on a toy example (Fig. 9). There is potentially a mixture of registration inaccuracies, and of mismatch between the displacement fields from the sequence registration and from the mesh simulation. The former might be addressed by better algorithms, although registering cine MR sequences is still highly challenging [START_REF] Tobon-Gomez | Benchmarking framework for myocardial tracking and deformation algorithms: an open access database[END_REF] due to the limited z-axis resolution and the lack of clear texture patterns attached to the myocardium, contrary to echocardiography or tagged-MRI. The latter might be improved by the personalization of the electromechanical model, using more advanced schemes [START_REF] Mollero | A multiscale cardiac model for fast personalisation and exploitation[END_REF], [START_REF] Neumann | A self-taught artificial agent for multiphysics computational model personalization[END_REF] and incorporating the match of displacement fields in the personalization metrics. In our study, the use of well-adjusted consensual values for all cases still provided simulations close to the original sequence (Fig. 4 and5). Taking advantage of this, our pipeline involved a much easier registration problem at each instant of the sequence. Better personalization might help but is not a prerequisite, as demonstrated on our data. Nonetheless, this strategy may be limited by artifacts from warping the first frame, and texture or gray-level differences with the original sequence. This registration quality was carefully checked and was not a major issue in our study.

Infarct detection and localization implies the use of dedicated algorithms [START_REF] Sudarshan | Automated identification of infarcted myocardium tissue characterisation using ultrasound images: a review[END_REF]- [START_REF] Duchateau | Infarct localization from myocardial deformation: prediction and uncertainty quantification by regression from a low-dimensional space[END_REF], whose performance depends on the quality of the features extraction (e.g. shape, motion and deformation), and therefore segmentation or registration techniques. Thus, we preferred to focus on the validation of the accuracy, quality, and usability of the simulated images assessed from the segmentation and registration point-of-view, as testing the pipeline on more advanced and applicationspecific algorithms may deviate from the original focus of the paper.

Segmentation was not used as a diagnosis tool as it does not allow to propagate the labeling of the myocardial zones (and in particular the diseased region) along the cycle, unlike registration. It only served for evaluation of the simulated images accuracy, and was therefore kept manual to prevent additional uncertainties coming from the use of an algorithm.

The registration used for the quantitative evaluation actually suffers the above-discussed limitations, and results in less contrasted differences between healthy and diseased behaviors, in contrast with the mesh simulations (these differences are also attenuated by the resolution differences between meshes and images). The optimization of a registration strategy was not the core of this paper and we preferred to test our images with a state-of-the-art algorithm. Nonetheless, the generated population still allows to see significant differences when the tissue viability decreases. As synthetic cases lie below and above the limit of the retained diagnosis method, this highlights a critical aspect for the generation of large databases (Fig. 2): better understanding the limits of the evaluated algorithms should take precedence over the race for the largest possible population.

The variety of the textures in the generated sequences was not addressed here. All the generated variants of a given sequence re-use the same texture, which is warped according to the differences between the original and the simulated motion. This may be a potential limitation for testing methods that actually make use of the texture information, such as the recent neural networks-based learning algorithms. The risk of unrealistic texture warping still exists with our method when deviating too far from the original motion, but this highly depends on the texture information of the original sequence and the prescribed motion of the simulated meshes (Fig. 13). On our dataset, this limit was substantially less problematic than in the results presented for highly textured modalities such as echocardiography [START_REF] Prakosa | Generation of synthetic but visually realistic time series of cardiac images combining a biophysical model and clinical images[END_REF] and tagged-MRI [START_REF] Zhou | Generation of realistic 4D synthetic CSPAMM tagged MR sequence for benchmarking cardiac motion tracking algorithms[END_REF]. Simulating the physics of imaging would be a robust way to tackle these issues, but at high computational cost. On the contrary, our method has the advantage of easily providing variants of a given case. We notably take advantage of existing mesh correspondences to generate new data without additional registration, at reduced computational cost and through highly parallelizable algorithms, which are key assets when targeting the design of large databases. Future improvements might incorporate texture synthesis on top of it, inspired from physicsbased [START_REF] Alessandrini | A pipeline for the generation of realistic 3D synthetic echocardiographic sequences: methodology and open-access database[END_REF] or patch synthesis [START_REF] Cordier | Extended modality propagation: image synthesis of pathological cases[END_REF] strategies.

Due to such texture modeling limits, we preferred to demonstrate the methodology on cine MR sequences. From the technical side, the interest is to support the design of computer-assisted methods that will tackle the processing of large databases with such sequences, such as the UK Biobank [START_REF] Petersen | Imaging in population science: cardiovascular magnetic resonance in 100,000 participants of UK Biobank -rationale, challenges and approaches[END_REF], eventually for the diagnosis of myocardial infarct. From the clinical side, applying such algorithms to cine MR would have potential for a first opinion diagnosis, which if positive and possible on the studied patient would lead to the recommendation of an extra LGE scan. Moreover, the proposed methodology is rather generic and could be used for echocardiography once better textures can be synthesized, for higher clinical impact.

The STACOM 2011 challenge data were already segmented in a format compatible with SOFA simulations. Focus was therefore set on the image synthesis steps coming after the segmentation of the initial frame. Atlas-based segmentation techniques could provide correspondences between the generated cases, of relevance for larger databases. Nonetheless, mesh fitting needs to be carefully checked to prevent the propagation of artifacts. Atlas-based template sequences [START_REF] Hoogendoorn | A high-resolution atlas and statistical model of the human heart from multislice CT[END_REF] might also be used, although this raises texture-related issues (variety of textures and potential warping artifacts).

Our strategy is in line with the recent initiatives for simulating and synthesizing realistic medical images [START_REF] Glatard | A virtual imaging platform for multi-modality medical image simulation[END_REF], [START_REF] Tsaftaris | Simulation and Synthesis in Medical Imaging[END_REF]. On the clinical problem chosen here for illustration purposes (assessment of the tissue viability and infarct diagnosis), expectations are high to bridge the gap between synthetic mesh simulations and real imaging data [START_REF] Duchateau | Infarct localization from myocardial deformation: prediction and uncertainty quantification by regression from a low-dimensional space[END_REF]. In this sense, we incorporated here different grades of tissue viability, and higher realism in the diseased zone generation algorithmnow initiated at the endocardial surface. Yet, aspects of higher realism of the pathological simulations (e.g. the definition of a neighboring partially affected region, and incorporating ventricular remodeling consecutive to an infarct) were outside the scope of our objectives here. Nonetheless, our method is readily applicable to these more realistic simulations, once available. Mixing synthetic and real data should be considered as a future target, as demonstrated on easier problems [START_REF] Heimann | Real-time ultrasound transducer localization in fluoroscopy images by transfer learning from synthetic training data[END_REF], provided issues of transfer learning are also taken into account.

V. CONCLUSION

We presented a new approach for the generation of synthetic cardiac image sequences from existing real sequences, which improves current methods in terms of image realism. Our pipeline, demonstrated on cine MR sequences, is highly suitable for the generation of large databases at reduced computational cost, and paves the ground for bridging the gap between synthetic mesh simulations and real imaging data [START_REF] Duchateau | Infarct localization from myocardial deformation: prediction and uncertainty quantification by regression from a low-dimensional space[END_REF].

ACKNOWLEDGMENTS

The authors acknowledge the partial support from the European Research Council (MedYMA ERC-AdG-2011-291080) and the European Union 7th Framework Programme (VP2HF FP7-2013-611823). They also thank R. Mollero and S. Giffard-Roisin (INRIA Asclepios, Sophia-Antipolis, France) for fruitful discussions about SOFA simulations, and L. Le Folgoc (Microsoft Research, Cambridge, UK) for his initial help on the implementation of [START_REF] Prakosa | Generation of synthetic but visually realistic time series of cardiac images combining a biophysical model and clinical images[END_REF].

SUPPLEMENTARY MATERIAL

Filename Description

Fig5.avi Animated view of Fig. 5 Fig9.avi Animated view of Fig. 9 Fig11 Best.avi Animated synthetic sequences for the best/medium/worst personalized cases that complement Fig. 11 Fig11 Medium.avi Fig11 Worst.avi http://www-sop.inria.fr/asclepios/docs/ImageSynthesisExamples.zip

  Fig.1. Processing pipeline used in[START_REF] Prakosa | Generation of synthetic but visually realistic time series of cardiac images combining a biophysical model and clinical images[END_REF]. The original cardiac sequence is "freezed" to the first instant via image registration (top row). Better performance is reached when initializing the registration between frames 0 and t by the composition of transformations first obtained from the registration of consecutive frames t and t + 1 (in practice, handled in the tangent space as the sum of velocity fields linking frames t and t + 1). Then, this "static" image is warped according to the simulation displacement field (bottom row). Consistency between the myocardium and the inside/outside anatomy is ensured by diffusion on this displacement field.

Fig. 2 .Fig. 3 .

 23 Fig.2. Schematic view of the tackled challenges for database generation, illustrated on the 30 variants obtained from case v2, arranged according to the disease grade (myocardial viability) and extent (diseased zone volume). Generating a synthetic database should provide a balanced amount of cases below and above the diagnosis limit of a given algorithm, within the limits of detectable disease and image realism.

Fig. 4 .

 4 Fig. 4. Comparison of LV volumes at end-systole: the 3D mesh simulation compared to the 3D segmentation from 2D slices (a) from the original sequence to evaluate the personalization error, and (b) from the synthesized healthy sequence to evaluate the accuracy of the generated images with respect to the mesh simulation. The black arrows point out a subject with representative personalization quality (v7, average personalization error in terms of distance between surfaces at end-systole: 3.0 mm).

Fig. 5 .

 5 Fig. 5. Simulated mesh superimposed to the original sequence of the representative subject v7 pointed out in Fig.4. An animated version of the figure is available as Supplementary Material6.

Fig. 6 .Fig. 7 .

 67 Fig. 6. Characteristics of the generated infarcts. Left: grade vs. volume. Right: starting location on the endocardial surface.

Fig. 8 .

 8 Fig. 8. Spread of the vector field ' 0 (t) from the electromechanical simulation (a), after diffusion to ensure continuity in the motion around the heart (b), and with cropping below the mitral annulus before the diffusion to reduce artifacts due to potential discrepancies between the original and simulated motion at this level (c).

1 .

 1 Crop and upsample sequence to isotropic resolution for t 2 [0, T ] do 2. Extract mesh local displacement: ' 0 (t) 3. Crop at the mitral annulus and diffusion: D(' 0 (t)) 4. Warp image at time 0 by D(' 0 (t))

Fig. 9 .

 9 Fig.9. Comparison between the method from[START_REF] Prakosa | Generation of synthetic but visually realistic time series of cardiac images combining a biophysical model and clinical images[END_REF] and our approach on a toy example (subject v7), demonstrating the limits of[START_REF] Prakosa | Generation of synthetic but visually realistic time series of cardiac images combining a biophysical model and clinical images[END_REF] by the generated MR image at end-systole. (a) and (b): first frame warped along the sequence (step 4 of Algorithm 1), referred to as the toy sequence. By definition, it matches the mesh simulation and mimics an ideal personalization. Thus, it also corresponds to the result of our method on this ideal case, as the warped sequence is used in lieu of the original sequence and therefore

Fig. 11 .

 11 Fig.11. Comparison of healthy and infarcted synthetic sequences along the cycle, for the best personalized case v6. The matching with the mesh contours is conserved, despite substantial changes in the simulation at the infarct level (red color), where contraction is lower and mesh contours are thinner. An animated version of the figure is available as Supplementary Material 6 , together with similar views of the medium and worst personalized cases.

Fig. 12 .

 12 Fig. 12. Regional volume change from end-diastole to end-systole computed on the simulated meshes and images. Comparison between the infarcted (red) and healthy (blue) simulations over the region labeled as infarct. Stars on the x-axis indicate statistically significant differences between the infarcted and healthy simulations (p-value < 0.01 from the Mann-Whitney U test).

Fig. 13 .

 13 Fig.[START_REF] Prakosa | Generation of synthetic but visually realistic time series of cardiac images combining a biophysical model and clinical images[END_REF]. Qualitative assessment of the simulated images with the highest tissue damage (8 cases with tissue viability < 10% and diseased zone volume > 10%), against the original and healthy simulated images. Warping artifacts are hardly noticeable for the healthy simulations. Slight myocardial shape differences can be noticed on the diseased cases with highly contrasted myocardial border in the original sequence.

TABLE I PARAMETERS

 I DEFINING THE HEALTHY-INFARCT RANGE USED FOR THE MESH SIMULATIONS (DETAILED SENSITIVITY ANALYSIS IN[START_REF] Marchesseau | Preliminary specificity study of the Bestel-Clement-Sorine electromechanical model of the heart using parameter calibration from medical images[END_REF]).

						[this paper]		[Duchateau et al. 2016]
	Category	Symbol	Description	Units Mitral annulus	RV	LV Healthy LV Infarct Healthy	Infarct
	Stiffness	c 1	Mooney-Rivlin modulus	Pa	40e4	2e4	4e4	4e4	2e4	2e4
		c 2	Mooney-Rivlin modulus	Pa	40e4	2e4	4e4	4e4	2e4	2e4
		K	Bulk modulus	Pa	3e6	3e6	3e6	3e6	3e6	3e6
	Contraction	0	maximum contraction	Pa	7.4e6	3.7e6	7.4e6	1e5	3.7e6	5e4
		k 0 katp	maximum stiffness contraction rate	Pa s 1	6e6 30	6e6 30	6e6 30	1e5 20	6e6 30	1e5 20
		krs	relaxation rate	s 1	70	70	70	80	70	80

  Preprint version accepted to appear in IEEE Transactions on Medical Imaging. Final version of this paper available at http://ieeexplore.ieee.org/ differences similarity metric4 .

	5. Register warped and original images at time t
	Output: vector field	original!warped (t)
	6. Mask this field by a diffused mask of the
	myocardium, and temporally smooth over 2 instants
	Output: vector field ⇤ original!warped (t)
	7. Warp original image at t by this field
	8. Downsample sequence to the original resolution

Preprint version accepted to appear in IEEE Transactions on Medical Imaging. Final version of this paper available at http://ieeexplore.ieee.org/

http://www-sop.inria.fr/asclepios/docs/ImageSynthesisExamples.zip


Zhenyu Wang

Jun Zheng
email: junzheng@seu.edu.cn

Yuying Wu

Nathalie Mitton

A Centrality-based RSU Deployment Approach for Vehicular Ad Hoc Networks

Keywords: centrality, deployment, roadside unit, RSU, VANET, vehicular ad hoc network I

This paper studies the RSU deployment problem in a 2-D urban or suburban road scenario of a vehicular ad hoc network (VANET). To optimize RSU deployment, we introduce the notion of centrality in a social network to RSU deployment, and use it to measure the importance of an RSU position candidate in RSU deployment. Based on the notion of centrality, we propose a centrality-based RSU deployment approach and formulate the RSU deployment problem as a linear programing problem with the objective to maximize the total centrality of all position candidates selected for RSU deployment under the constraint of a given deployment budget. To solve the formulated problem, we analogize the problem to a 0-1 Knapsack problem and thus employ a 0-1 Knapsack algorithm to solve the problem. In the analogy, the budget in the RSU deployment problem is analogous to the bag's capacity in the Knapsack problem, the cost of deploying an RSU is analogous to an item's weight, and the centrality of a position candidate is analogous to an item's value. Simulation results show that the proposed centrality-based deployment approach can effectively improve the efficiency of the RSU deployment in terms of the coverage time ratio as compared to a random deployment approach.

INTRODUCTION

Roadside Units (RSUs) are an important component of vehicular ad hoc networks (VANETs) [START_REF] Hartenstein | A tutorial survey on vehicular ad hoc networks[END_REF] [START_REF]Standard Specifications for Telecommunications and Information Exchange Between Roadside and Vehicle Systems-5 GHz Band Dedicated Short Range Communications (DSRC) Medium Access Control (MAC) and Physical Layer Specifications[END_REF]. In a VANET, RSUs are deployed at intersections or some points along a road to help improve network connectivity, data delivery, and thus network services to vehicles. Therefore, RSU deployment has a big impact on the network performance and becomes an important issue in the design of a VANET.

In general, RSU deployment is costly. To achieve a good tradeoff between deployment cost and network performance, it is expected to optimize the deployment of RSUs in a VANET. To address this challenge, considerable research work has been conducted on the optimization of RSU deployment for VANETs [START_REF] Rashidi | Placement of road side units for floating car data collection in highway scenario[END_REF][START_REF] Aslam | Optimal roadside units placement along highways[END_REF][START_REF] Tao | RSU deployment scheme with power control for highway message propagation in VANETs[END_REF][START_REF] Chi | Intersection-priority based optimal RSU allocation for VANET[END_REF][START_REF] Barrachina | Road side unit deployment: a density-based approach[END_REF][START_REF] Aslam | Optimal roadside units placement in urban areas for vehicular networks[END_REF][START_REF] Lin | Optimal roadside unit deployment in vehicle-to-infras t--ructure communications[END_REF][START_REF] Liang | Optimal placement and configuration of roadside units in vehicular networks[END_REF][START_REF] Kchiche | Centrality-based access-points deployment for vehicular networks[END_REF]. However, existing RSU deployment strategies exhibit their own limitations. For example, some of them only consider the RSU deployment in 1-D road scenarios [START_REF] Rashidi | Placement of road side units for floating car data collection in highway scenario[END_REF][START_REF] Aslam | Optimal roadside units placement along highways[END_REF][START_REF] Tao | RSU deployment scheme with power control for highway message propagation in VANETs[END_REF]; some of them only consider the same cost of deploying RSUs in different scenarios [START_REF] Rashidi | Placement of road side units for floating car data collection in highway scenario[END_REF][6] [START_REF] Kchiche | Centrality-based access-points deployment for vehicular networks[END_REF]; some of them require many parameters which may change over time, such as the density of traffic [START_REF] Aslam | Optimal roadside units placement along highways[END_REF] [START_REF] Barrachina | Road side unit deployment: a density-based approach[END_REF]. To improve the efficiency of RSU deployment, it is necessary to further study the RSU deployment problem and explore more effective RSU deployment approaches.

In this paper, we focus on the RSU deployment problem in a 2-D urban or suburban road scenario. To optimize the RSU deployment, we introduce the notion of centrality in a social network to RSU deployment, and use it to measure the importance of an RSU position candidate in RSU deployment. Based on the notion of centrality, we propose a centrality-based RSU deployment approach and formulate the RSU deployment problem as a linear programing problem with the objective to maximize the total centrality of all position candidates selected for RSU deployment under the constraint of a given deployment budget. To solve the formulated problem, we analogize the problem to a 0-1 Knapsack problem [START_REF] Sanjay | Unbounded Knapsack Problem : dynamic programming revisited[END_REF] and thus employ a 0-1 Knapsack algorithm to solve the problem. In the analogy, the budget in the RSU deployment problem is analogous to the bag's capacity in the Knapsack problem, the cost of deploying an RSU is analogous to an item's weight, and the centrality of a position candidate is analogous to an item's value. Simulation results show that the proposed centrality-based deployment approach can effectively improve the efficiency of RSU deployment in terms of the coverage time ratio as compared to a random deployment approach.

The remainder of this paper is organized as follows. Section II reviews related work on the RSU deployment problem in VANETs. Section III formulates the RSU deployment problem considered in this paper and presents the proposed centrality-based RSU deployment approach. Section IV evaluates the performance of the proposed deployment approach through simulation results. Section V concludes this paper.

II. RELATED WORK

RSU deployment has been widely studied in literature [START_REF] Rashidi | Placement of road side units for floating car data collection in highway scenario[END_REF][START_REF] Aslam | Optimal roadside units placement along highways[END_REF][START_REF] Tao | RSU deployment scheme with power control for highway message propagation in VANETs[END_REF][START_REF] Chi | Intersection-priority based optimal RSU allocation for VANET[END_REF][START_REF] Barrachina | Road side unit deployment: a density-based approach[END_REF][START_REF] Aslam | Optimal roadside units placement in urban areas for vehicular networks[END_REF][START_REF] Lin | Optimal roadside unit deployment in vehicle-to-infras t--ructure communications[END_REF][START_REF] Liang | Optimal placement and configuration of roadside units in vehicular networks[END_REF][START_REF] Kchiche | Centrality-based access-points deployment for vehicular networks[END_REF]. Some studies only consider the RSU deployment problem in 1-D road scenarios such as highways. For example, Rashidi et al. studied RSU Placement for floating car data collection in a highway scenario in [START_REF] Rashidi | Placement of road side units for floating car data collection in highway scenario[END_REF]. The trade-offs between the RSU interval and other system parameters were studied and a heuristic algorithm was proposed to calculate the interval between RSUs. In [START_REF] Aslam | Optimal roadside units placement along highways[END_REF], Aslam and Zou proposed an optimal placement of RSUs along highways with the goal of minimizing the average time taken for a vehicle to report an event of interest to a nearby RSU with a limited number of RSUs. Not only the cost but also the energy consumption of RSUs was considered in the deployment. In [START_REF] Tao | RSU deployment scheme with power control for highway message propagation in VANETs[END_REF], Tao et al. proposed an RSU deployment scheme with power control for highway message propagation in VANETs and aimed to minimize the energy consumption of RSUs and maximize the network performance.

Apart from the studies on 1-D road scenarios, there are some studies on RSU deployment in 2-D road scenarios such as urban streets. For example, Chi et al. proposed an intersection-priority based RSU placement approach in [START_REF] Chi | Intersection-priority based optimal RSU allocation for VANET[END_REF] to minimize RSU setup cost while maximizing connectivity between RSUs. In [START_REF] Barrachina | Road side unit deployment: a density-based approach[END_REF], Barrachina et al. proposed a density-based RSU deployment approach to obtain an efficient system with the lowest possible cost to alert emergency services in case of an accident. In [START_REF] Aslam | Optimal roadside units placement in urban areas for vehicular networks[END_REF], Aslam et al. studied optimal RSU placement in urban areas and also considered to minimize the average reporting time with a limited number of RSUs. In [START_REF] Lin | Optimal roadside unit deployment in vehicle-to-infras t--ructure communications[END_REF], Lin studied the optimal RSU deployment problem in vehicle-to-infrastructure communication and formulated the problem as a constrained optimization problem. In [START_REF] Liang | Optimal placement and configuration of roadside units in vehicular networks[END_REF], Liang et al. studied the optimal RSU placement and configuration problem in vehicular networks with the objective to minimize the total cost to deploy and maintain RSUs and formulated the problem as a linear programming problem. In [START_REF] Kchiche | Centrality-based access-points deployment for vehicular networks[END_REF], Kchiche and Kamoun studied centrality-based access-points deployment for vehicular networks, and considers the central positions of RSU regrading vehicle traffic. Moreover, they formulated the equidistance-based deployment problem as a coverage problem. However, the proposed heuristic algorithm is neither optimal nor scalable.

III. CENTRALITY-BASED RSU DEPLOYMENT APPROACH

In this section, we present the proposed centrality-based RSU deployment approach. We first describe the network model and then introduce the notion of centrality to RSU deployment and formulate the RSU deployment problem considered in this paper. Finally, we present a Knapsack algorithm to solve the formulated problem.

A. Network model

We consider a 2-D urban or suburban road scenario in a VANET. The urban area has an irregular road topology with unevenly distributed intersections, while the suburban area has a more regular road topology with more evenly distributed intersections, as shown in Figure 1. In the road scenario, RSUs are deployed at intersections or some points along a road. The cost of deploying an RSU depends not only on the cost of the RSU itself, but also on the cost of installing supporting facilities at a position candidate.

B. Centrality-based RSU deployment approach

Centrality is a notation originally introduced to measure the importance of a node in a social network. In a VANET, it is expected to deploy RSUs at positions where a vehicle has more opportunities to pass through. Thus, we introduce the notion of centrality to RSU deployment in a VANET and use it to measure the importance of an RSU position candidate in RSU deployment. The centrality of a node can be defined from different perspectives. Here, we introduce two centrality metrics: degree centrality and closeness centrality, which are described below.

(1) Degree Centrality

In a social network, the degree centrality of a node measures the importance of the node based on the node's degree of direct connections to other nodes. For a node k, its degree centrality is defined as

1   N d k dk  , ( 1
)
where d k is the degree of node k and N is the total number of nodes in the network.

In RSU deployment, an RSU position candidate can be analogous to a node in a social network. Moreover, an RSU position candidate is usually located at an intersection or at some point along a road. In this case, d k represents the number of RSU position candidates to each of which candidate k has a road directly connected. N is the total number of RSU position candidates in the network. Obviously, a position candidate with a larger value of dk  means that there are more roads connected to this position and accordingly there are more vehicles passing through the position. In another word, this means that a position candidate with a larger value of dk  can provide more opportunities for a vehicle to pass through. From this perspective, such a position candidate should have a higher probability to be selected for RSU deployment.

(2) Closeness Centrality

In a social network, the closeness centrality of a node measures the importance of the node based on the node's total distance to all other nodes in the network. For a node k, its closeness centrality is defined as

    N i ki ck d N 1 1  , i k  , (2)
where d ki is the geodesic distance between node k and node i, and N is the total number of nodes in the network.

In RSU deployment, d ki represents the shortest path between position candidate k and position candidate i. A larger value of ck  means a smaller total shortest path distance from candidate k to all other candidates. This means that a position candidate with a larger value of ck  allows a vehicle to enter its coverage more quickly and thus can provide better quality of service to the vehicle. From this perspective, such a position candidate should also have a higher probability to be selected for RSU deployment.

C. Selection of RSU position candidates

In an urban or suburban area, RSUs can usually be deployed at intersections to provide services to passing vehicles. However, when a road segment is long and the coverage of RSUs at intersections cannot cover the road segment, an RSU can be deployed at some point along the road segment to avoid interruptive connections and improve the quality of service to vehicles.

Assume that the coverage of an RSU is R u and the length of a road segment between two adjacent intersections is L. The selection of RSU candidates is performed according to the following rules: a) For a road segment with 2  u L R , the intersections of the road segment are selected as position candidates. The number of position candidates along the road segment is two. b) For a road segment with 2  u L R , the intersections of the road segment are selected as the position candidates. Meanwhile, a certain points along the road segment are also selected as position candidates, which are uniformly distributed along the road segment. The number of position candidates along the road segment is given by ( 2

)/2      u u L R R .

D. Problem formulation

Consider a two-dimensional urban or suburban road scenario with a set of N RSU position candidates, as shown in Figure 1. Let c k denote the deploying cost of RSU position candidate k; k  denote the centrality of position candidate k, k=1, 2, …, N, which can be calculated based on either Eq. (1) or Eq. (2); and

} 1 , 0 {  k X
denote a boolean variable indicating whether a position candidate is selected for RSU deployment, which is defined as

    selected not is candidate 0 selected is candidate 1 k k X k , N k ,..., 2 , 1  .
Moreover, assume that the RSU deployment budget for this road area is B.

Given the RSU deployment budget and the above definitions, the RSU deployment problem considered in this paper can be formulated into a linear programing problem with the objective to maximize the RSU coverage, i.e., to find a subset of RSU position candidates or determine ) ,..., 2 , 1 ( N k X k  so that the total centrality of the candidates contained in the subset is maximized, i.e.,

  N k k k X 1 max  , (5)
subject to

B X c N k k k   1 , (6)
where k  can be either dk  or ck  .

It is found that the above formulated RSU deployment problem can be analogous to the classical 0-1 Knapsack problem [START_REF] Sanjay | Unbounded Knapsack Problem : dynamic programming revisited[END_REF], which can be solved by using a Knapsack algorithm.

E. RSU deployment algorithm

Next we present the Knapsack algorithm to solve the problem.

(1) Analogous to the 0-1 Knapsack problem

The 0-1 Knapsack problem is a classical combinatorial optimization problem. Suppose that there is a bag which can hold a total weight W and there are n kinds of items. Each item has a weight w k and a value v k (k= 1, 2, …, n). Thus, the 0-1 Knapsack problem is how to find a subset of these items to put into the bag so that the total value of items in the bag is maximized. Obviously, the formulated RSU deployment problem can be analogous to the 0-1 Knapsack problem if budget B is analogous to the bag's capacity W, cost c k is analogous to each item's weight w k , and centrality k  is analogous to each item's value v k . Therefore, we can use the Knapsack algorithm to solve the formulated RSU deployment problem.

(2) Knapsack algorithm for RUS deployment

The Knapsack algorithm is a dynamic programming algorithm. According to the formulation of the RSU deployment problem, a solution to the problem is a subset of all N position candidates, whose total centrality is maximized with their total budget constrained by B. To find the solution to the problem, define m(i,b) as the total centrality of a subsolution to candidates i, i+1, . . . , N given budget b. By the nature of the optimal sub-structure of the 0-1 knapsack problem, we can establish an iteration relation of m(i,b), i.e,

            i i i i c b b i m c b c b i m b i m b i m 0 ) , 1 ( } ) , 1 ( ), , 1 ( max{ ) , (  , (7) ( , ) 0 0         N N N b c m N b b c , ( 8
) N i ,..., 2 , 1 
Obviously, the value of the optimal solution is m(1,B), which can be calculated using Eq. ( 7) through iterations, starting from m(N,b) with Eq. ( 8). After a sequence of iterations, we can obtain a table of m(i,b), which can be used to deduce the solution to the problem, i.e., a subset of position candidates that are selected for RSU deployment. The deduction process starts at m(1,B) and traces backwards where the optimal solution was obtained. If m(i,b) = m(i+1,b), candidate i is not part of the solution, and then continue tracing with m(i+1, b). Otherwise, candidate i is part of the solution, and then continue tracing with m(i+1, b-c i ). The pseudo-code of the algorithm is given below.

1 1 ( , ) ( 1,
)

[ ] 0 [ ] 1 [ ] ( , )?1:0          i for i to N do if m i B m i B then X i else X i B B c X N m N B Return X

IV. SIMULATION RESULTS

In this section, we evaluate the performance of the proposed centrality based RSU deployment approach (CDA) through simulation results. For this purpose, we developed a simulator using C++ programing language and conducted a series of experiments. Most existing RSU deployment approaches are based on some statistical parameters of a road scenario, such as traffic density, while our proposed approach can find a deployment solution without using any statistical parameters of the road scenario. Thus, it is hard to find a suitable existing RSU deployment approach to compare with in evaluation. Instead, we compare CDA with a random deployment approach (RDA). For CDA, we consider two metrics of centrality: degree centrality and closeness centrality, and use CDA-DC and CDA-CC to represent, respectively.

In the simulation experiments, we considered two road scenarios or topologies as shown in Figure 2. The number on each road segment represents the length of the road segment. A circle represents an RSU position candidate. A shaded circle represents a position candidate at an intersection with traffic light. A non-shaded circle represents a position candidate at an intersection without traffic light or at some point along a road segment. The RSU position candidates are determined based on the road topology and the selection rules given in Section III.C. Moreover, we introduce the coverage time ratio as a performance metric to evaluate the efficiency of RSU deployment, which is defined as the ratio of the amount of time a vehicle is covered by RSUs to the total amount of time the vehicle moves in the road area. Obviously, a larger value of the coverage time ratio means that a vehicle has a larger probability to be covered and thus served by RSUs during its movement in the area. Moreover, we assume that there are M vehicles moving in the road area and each vehicle moves at a constant speed in a free-flow state.

We consider different deployment cost at different position candidates. The parameters used in the simulation experiments are summarized in Table I. Figure 3 shows the coverage time ratio with CDA-DC, CDA-CC, and RDA, respectively, for different number of vehicles in the urban area. It is seen that the coverage time ratios with both CDA-DC and CDA-CC are significantly larger than that with RDA. This means RSU deployment with CDA can provide better RSU services to vehicles than that with RDA. On the other hand, it is seen that the coverage time ratio decreases with the number of vehicles increasing. This is because with the increase of the number of vehicles, more vehicles may move on road segments without the coverage of RSUs, making the coverage time ratio decrease.

Figure 4 shows the coverage time ratio with CDA-DC, CDA-CC, and RDA, respectively, for different deployment budgets in the urban area, when the number of vehicles M is 20. It is also seen that the coverage time ratios with both CDA-DC and CDA-CC are significantly larger than that with RDA. The coverage time ratio increases with the increase of the budget. Moreover, with the increase of the budget, the coverage time ratio with RDA tends to approach that with CDA. This is because in this case both CDA and RDA tend to deploy RSUs at all position candidates, which cover the whole area.

Figure 5 shows the coverage time ratio with CDA-DC, CDA-CC, and RDA, respectively, for different number of vehicles in the suburban area. Similarly, it is seen that the coverage time ratios with both CDA-DC and CDA-CC are significantly larger than that with RDA. This means RSU deployment with CDA can provide better RSU services to vehicles than that with RDA.

Figure 6 shows the coverage time ratio with CDA-DC, CDA-CC, and RDA, respectively, for different deployment budgets in the suburban area. It is also seen that the coverage time ratios with both CDA-DC and CDA-CC are significantly larger than that with RDA. The coverage time ratio increases with the increase of the budget. Moreover, with the increase of the budget, the coverage time ratio with RDA tends to approach that with CDA. V. CONCLUSIONS In this paper, we introduced the notion of centrality in a social network to RSU deployment, and proposed a centrality-based RSU deployment approach for a 2-D urban or suburban road scenario. The centrality-based RSU deployment problem was formulated as a linear programing problem with the objective to maximize the total centrality of all position candidates selected for RSU deployment under the constraint of a given deployment budget. The formulated problem is analogous to a 0-1 Knapsack problem and can be solved by employing a 0-1 Knapsack algorithm. Simulation results show that the proposed centrality-based deployment approach can effectively improve the efficiency of RSU deployment in terms of the coverage time ratio as compared to a random deployment approach.

Figure 1

 1 Figure 1 Network model

1 :

 1 based on centrality C N B Input the cost of each RSU position candidate C c c the number of RSU position candidates N the budget of deployment B Output a set of RSU deployment positions X Initia

Figure 2

 2 Figure 2 Road topology

Figure 3 Figure 4 Figure 5 Figure 6

 3456 Figure 3 Coverage time ratio vs number of vehicles (B=200000¥)

This work was supported by the National Natural Science Foundation of China under Grant No. 61372105 and the Six Talent Peaks Project in Jiangsu Province under Grant No. 2013-DZXX-010.