Zhiqiang Wang
email: zhiqiang.wang@univ-nantes.fr

Catherine Da Cunha

Mathieu Ritou

Benoît Furet

on Changeable, Agile, Reconfigurable and Virtual Production Comparison of K-means and GMM methods for contextual clustering in HSM

Keywords: HSM, Industry4.0, clustering, K-means, GMM

High speed machining (HSM) is widely used for the manufacturing of aircraft structures, turbine blades, etc. It greatly increases the efficiency and automation for the machining. However, in HSM, operators cannot detect incidents when they manage several machines of a production cell. Robust monitoring systems are required to protect the machine tool and the high value added parts. In the global context of the Industry 4.0, abundant digital data is available in a modern manufacturing company and could be used to turn the machines-tools smarter and to support the decision making of the operational management. One of the first step of data mining approach is the accurate selection of relevant. To do so, the raw data need to be classified into different contextual clusters. This paper compares two different methods of the unsupervised classification of machining context: K-means and GMM (Gaussian Mixture Model). It was found that GMM method can classify correctly the machining context, whereas K-means is not suitable.

Introduction

Compared with the traditional machining, high speed machining (HSM) has greatly increased the cutting speeds. However, in HSM, the operator cannot detect incidents when he manages several machines of a production cell. Moreover, he has not enough time to rapidly stop the machine. Consequently, monitoring systems are required to protect the machine tool and the high value added parts. Especially in aeronautic parts manufacturing, those aeronautic pieces are high added value parts. To guarantee the quality of the surface of the parts, and to avoid the tool broken and prevent the great chatter. ( [START_REF] Godreau | Extraction des connaissances à partir des données de la surveillance de l'usinage[END_REF] presents the link between the quality of surface and chatter), a smart machining monitoring system need to be set up. [START_REF] Zhou | Review of tool condition monitoring methods in milling processes[END_REF] presents a review of the state-of-the-art methods employed for conducting tool condition monitoring (TCM) in milling processes including sensors, feature extraction, and monitoring models.

Abundant digital data is available in a modern manufacturing company and could be used to turn the machinestools smarter and to support the decision making of the operational management [START_REF] Emeric | Un Système à Base de Connaissances pour la gestion de données d'usinage dans une perspective de fouille de données[END_REF]. Data and data analytics are considered ever more important in manufacturing. [START_REF] Lenz | Holistic approach to machine tool data analytics[END_REF] discusses the current situation exhibited in machine tool data analytics in manufacturing, and proposes a holistic approach to machine tool data analytics in order to tackle some of the identified shortcomings of current practices. During data analytics, the objective is to set up the adaptive criterions (KPI) and calculate the threshold for every criterion. Some algorithms will be used to exploit these criterions and thresholds. Data mining has been widely used in manufacturing. [START_REF] Da Cunha | Data mining for improvement of product quality[END_REF] Uses data mining to improve the products quality in assembly production. [START_REF] Martínez-Arellano | In-process Tool Wear Prediction System Based on Machine Learning Techniques and Force Analysis[END_REF] presents an in-process tool wear prediction method based on deep learning CNN (Convolutional Neural Network). [START_REF] Denkena | Self-optimizing cutting process using learning process models[END_REF] presents a method that uses a support vector machine (SVM) as a machine learning approach to model the obtained process data. There are some other algorithms in machine learning which can be used in manufacturing data mining, ex. K-means [START_REF] Julien | Fouille de données / Data mining Université Lumière Lyon[END_REF]. Also, there are some Statistics algorithms interesting for our database, such as Gaussian mixture models GMM [START_REF] Reynolds | Gaussian mixture models[END_REF]. Next chapter will introduce the methods for the data analytics and how it will be used to the monitoring system for the machine-tool.

Proposed approach

LS2N has developed Emmatools during the UsinAE project (Fig. 1). It collects the data from the CNC (Computer Numerical Control) and from additional sensors like accelerometers, every 0.1s. This project finally led to the creation of criteria (SBBPFO or SBN) to analyse wear variations at different time scales [START_REF] De Castelbajac | Monitoring of distributed defects on HSM spindle bearings[END_REF]. Then, Godreau et al. [START_REF] Godreau | Continuous improvement of HSM process by data mining[END_REF] proposed 3 criterions to detect defects during machining: chatter; tool breakage; and collision. Abundant digital data have been collected by Emmatools. As cited in [START_REF] Lemberger | Big Data et Machine Learning-Manuel du data scientist-2e éd[END_REF], data clustering is very important at the beginning of data mining. There are some criterions have been calculated in different clusters. For example, criterion of chatter is calculated in cluster when tool is cutting materials; while criterion of tool breakage is calculated when tool is out of materials [START_REF] Godreau | Continuous improvement of HSM process by data mining[END_REF]. A better selection of data is better for a more relevant calculation of criterions (KPI). For example, the tool cuts materials only when the machine-tool is not stopped (i.e. feed rate (Vf) of spindle is not null). Also, when tool is cutting materials, the tool must be rotating. Eventually, the machine-tool can (de)accelerate. The varying Vf involves varying cutting conditions. Therefore, classifying the raw data by the feed rate Vf (m/min) into three clusters is necessary: the machine-tool can be stopped or moving at constant or varying speed.

The ). And the objective is to find the threshold T ΔVf which can define which data belongs to machine-tool moving at constant speed and which data belongs to machine-tool moving at varying speed. Moreover, the data whose Vf are less than the threshold T ΔVf will be considered as Vf null (machine-tool stops) in industry. Which means, firstly, the data whose spindle feed rate Vf less than T ΔVf is to be found, and is labelled 'Vf null' (machine-tool stops). And then, the data whose Vf bigger than T ΔVf , is labelled 'Vf not null'. And in the cluster 'Vf not null', the data whose ΔVf less than threshold T ΔVf is to be found, and is labelled 'machine-tool moves at constant speed'. Also in the cluster 'Vf not null', the data whose ΔVf bigger than threshold T ΔVf is labelled 'machine-tool moves at varying speed'. See the Fig. 2.

K-means

To classify the machine-tool context information into 3 clusters, an unsupervised machine learning method named K-means may be helpful. K-means clustering aims at partitioning n observations into k clusters in which each observation belongs to the cluster with the nearest mean, serving as a prototype of the cluster [START_REF] Julien | Fouille de données / Data mining Université Lumière Lyon[END_REF]. There are three clusters: machine-tool stops; machine-tool moves at constant or varying speed. The data during one day from an aircraft company was taken for example. Firstly, the data with Vf equal to 0 m/min were removed from the selection. And the K-means was applied for the ΔVf (m/s²) for this given day. The threshold of machine-tool moving at constant or varying speed locates on the boundary between two clusters: T ΔVf =1 (m/s²). The performance of this classification can be tested by manual mining of feed rate Vf classification in 3 clusters, by detecting the errors of classification. Fig. 3 illustrates the results of contextual classification of Vf, during a period of 30s. The Vf is labeled with different colors according to their clusters. Machine-tool stops are presented by green diamond; Machine-tool moving at constant speed are presented by red asterisk while machine-tool moving at varying speed are presented by blue plus. It is found manually that, during the Vf acceleration, there are many red asterisk (cluster machine-tool moving at constant speed) in this acceleration period. Which means the machine-tool is accelerating while the K-means suppose that these points are in the cluster machine-tool moving at constant speed. And there are many other classification errors in the machine-tool moving at constant speed (cycle purple). It means the threshold found by K-means T ΔVf = 1 (m/s²) is too high. In fact, there is a limitation when using K-means. K-means assumes that it deals with spherical clusters and that each cluster has roughly equal numbers of observations. For this case, the distribution of ΔVf is not spherical. And there are much more population of ΔVf around zero, compared to other population. Which means K-means will partition the data where machine-tool moving at varying speed into the cluster machine-tool moving at constant speed. So other data clustering methods will be tested.

GMM classification

Gaussian mixture models (GMM) [START_REF] Reynolds | Gaussian mixture models[END_REF] are often used for data clustering. Usually, fitted GMMs cluster by assigning query data points to the multivariate normal components that maximize the component posterior probability given the data. Its advantage is showing how to fit a GMM to data, cluster using the fitted model, and estimate component posterior probabilities. In this case, the same data as in chapter 3 is to be used for example. The Vf=0 m/min was taken away too, and the same column ΔVf is created also. Then GMM was to be used to classify ΔVf into 2 clusters: machine-tool moves at constant speed, machine-tool moves at varying speed. Here, the density of of ΔVf during this day is drawn. There are 2 great populations of ΔVf: ΔVf around zero (Y1) which is fine and high; ΔVf distributed along the whole day (Y2) which is wide and low. Therefore, 2 Gaussian distributions are set up to model the true data. The sum of the probability of these 2 Gaussians along with the ΔVf (Y3) is drawn to be compared with our true density of probability of ΔVf (Histogram of ΔVf). The result is in Fig. 4: the abscissa X is ΔVf, the abscissa Y is the density of probability in scale linear. It is found that the sum of these 2 Gaussians modeled well our true data. Next step, the population of Y1 (fine and high) is supposed to follow the distribution Gaussian. To find the threshold of the machine-tool moving at constant speed or varying speed, the σ, 2 σ, 3 σ and 4 σ will be tested by drawing the Vf along with the time, and zoom a period of 30 seconds to analyze.

Results and discussion

As the ΔVf within Y1 follows a normal distribution, the threshold T ΔVf can be defined to be equal to σ, 2σ, 3σ or 4σ. (While the mean of normal distribution μ = 0, and σ is the standard deviation of ΔVf within Y1). Because [-2σ, 2σ] represents the 95.45%, and [-3σ, 3σ] represents the 99.73% of all the data within Y1. For this day, the threshold T ΔVf =2σ = 0.0033 m/s² and 3σ = 0.005 m/s² as well as 4σ = 0.0067 m/s². To verify which threshold is the most suitable for this case, the feed rate Vf can be classified into 3 clusters by using these 3 thresholds. The same period 30s is always taken to be analyzed. Firstly, the classification by threshold T ∆N =1σ= 0.00167 m/s² is tested. The result is in the Fig. 5 (a): there is no wrong classification during acceleration which is better than the K-means method. However, just at the beginning of the acceleration, there is a blue plus point while it should be green. It is because the 1σis too small. And then the threshold 2σis tested in the Fig. 5 (b): there is no wrong classification before and during acceleration which is better than 1σ. And then the threshold 3σis tested in the Fig. 5 (c): classification by threshold 3σ is not better than that of 2σ. To confirm this threshold, the threshold 4σis tested in the Fig. 5 (d): there is a classification error at the end of this 30s, during the varying speed Vf. It is because its ΔVf=0.0059<4σ, 4σsuppose that this point belongs to machine-tool moving at constant speed. 4σ is too big for threshold. Therefore, the threshold T ∆N =2σ=0.0033 m/s² can be set as the threshold for this day's data to classify machine-tool stops and machine-tool moving at constant or varying speed. See the Fig. 5.

Conclusion

Two methods of unsupervised machine learning, K-means and GMM, have been tested for contextual classification, based on the machine-tool feedrate, into 3 clusters: machine stopped; constant or varying move. The results are as below:

 K-means is not suitable because the data is not spherical data (the population of ΔVf around zero is too large). However, GMM can classify the raw data well.  According to GMM classification and statistical analyzes, the distribution of ΔVf in the Y1 population follows a normal distribution, the threshold can be defined as T ΔVf = 2σ=0.0033m/s².  The data whose Vf is less than T ΔVf is labeled as 'machine-tool stopes'; the data whose Vf is greater than 2σ and its ΔVf is less than 2σ is labeled as 'machine-tool moves at constant speed'; the data whose Vf and ΔVf are greter than 2σ is labeled as 'machine-tool moves at varying speed'.  The threshold was chosen at 2σ by manual mining, through the verification of classification results.  As the raw data has been classified into 3 clusters, the new KPIs can be calculated in each cluster in the future. Such as, the tool will cut materials linear in the cluster 'machine-tool moves at constant speed' while the tool will cut materials in bending surface in the cluster 'machine-tool moves at varying speed'.

Fig. 1 .

 1 Fig. 1. Emmatools for process monitoring and data collection.

Fig. 2

 2 represents the flow chart how to classify the raw data. A criterion for classification is set up: It is equal to the half of the difference between the Vf the time after and the Vf the time before (∆𝑉𝑓 = 𝑉𝑓 𝑛+1 -𝑉𝑓 𝑛-1 2

Fig. 2 .

 2 Fig. 2. Process of classification according the Feed rate Vf and ΔVf

Fig. 3 .

 3 Fig. 3. Vf-T according to K-means classification during one day (zoom 30 seconds)

Fig. 4 .

 4 Fig. 4. Density of probability of ΔVf for one day modeled by GMM (2 Gaussians).

Fig. 5 .

 5 Fig. 5. (a) Vf-T according to 1σ classification of ΔVf during 30 seconds; (b) Vf-T according to 2σ classification of ΔVf during 30 seconds; (c) Vf-T according to 3σ classification of ΔVf during 30 seconds; (d) Vf-T according to 4σ classification of ΔVf during 30 seconds;

Acknowledgements

This research was funded by the French National Research Agency (ANR). We gratefully acknowledge for the financial support within the ANR SmarEmma project (ANR-16-CE10-0005).


Pierre Biasutti

Jean-Franc ¸ois Aujol

Mathieu Brédif

Aurélie Bugeau

Fast Image and LiDAR alignment based on 3D rendering in sensor topology

Mobile Mapping Systems are now commonly used in large urban acquisition campaigns. They are often equiped with LiDAR sensors and optical cameras, providing very large multimodal datasets. The fusion of both modalities serves different purposes such as point cloud colorization, geometry enhancement or object detection. However, this fusion task cannot be done directly as both modalities are only coarsely registered. This paper presents a fully automatic approach for LiDAR projection and optical image registration refinement based on LiDAR point cloud 3D renderings. First, a coarse 3D mesh is generated from the LiDAR point cloud using the sensor topology. Then, the mesh is rendered in the image domain. After that, a variational approach is used to align the rendering with the optical image. This method achieves high quality results while performing in very low computational time.

Results on real data demonstrate the efficiency of the model for aligning LiDAR projections and optical images.

Introduction

Over the past decades, the interest in urban acquisition systems has been growing continuously, especially Mobile Mapping Systems (MMS). These systems are vehicles equipped with many sensors that produce different modalities in order to acquire all the details of a scene. These systems are used on wide acquisition campaigns in cities, roads, highways, resulting in the production of very large -multimodal -datasets. Among all available sensors often met on such systems, 3D LiDAR sensors joint with optical cameras enable a geometrical acquisition of the scene as well as the acquisition of textures and colors. However, due to the complexity of such acquisition systems, the calibration from one sensor to the other does not generally meet pixel accuracy. This can be caused by the instability of the sensors throughout a mobile acquisition, where the calibration slowly deteriorates while the systems are being operated. Therefore, the different modalities are slightly misaligned which can compromise further processing requiring data fusion, e.g. point cloud colorization or multimodal object detection which may be required for autonomous driving applications or to filter out non-permanent objects for cartographic purposes. It is possible to interactively reduce this misalignment, but most of the time this solution is not suitable as the datasets are typically composed of thousands of examples. The automatic alignement of LiDAR data to optical image is therefore a crucial issue.

The problem of LiDAR to image alignment rises several issues. First of all, the comparison between the two modalities can only be done if they share common attributes (colors or reflectances). However, in many systems, each sensor solely acquires a specific aspect of the scene. Moreover, optical sensors and LiDAR sensors are located at different positions on the MMS. This implies that the different sensors do not acquire the scene from the same point of view, resulting in visual ambiguities. The correlation Email address: pierre.biasutti@labri.fr (Pierre Biasutti) between both modalities is therefore irrelevant for some parts of each data.

The paper contribution is two-fold, as illustrated Figure 1: first, we propose a very fast approach for mesh generation from 3D LiDAR data using sensor topology. The second contribution of the paper is an extension of a variational multimodal registration method (Sutour et al.) to the problem of LiDAR to image alignment with higher degree of freedom.

The paper is organized as follows: first, an overview of the related works is presented. Then, the mesh generation from the sensor topology is explained. After that, the extension of the variational method is detailed. Finally, evaluation and results are shown and a conclusion is drawn.

Previous works

Multi-modal registration has been a subject of interest over the past decades. In this section, previous works on multi-modal registration as well as previous works on LiDAR to optical image are introduced.

Multi-modal image registration

In computer vision, registration methods often consist in the detection and the matching of corresponding features from two different modalities to recover the 2D transformation that provides the best alignment between the two input images. It often assumes that the displacement between both image is small, otherwise there would not be any good alignment between the two images because of the differences in perspective. Feature points are extracted using common methods (SIFT (Lowe) or SURF [START_REF] Bay | SURF: Speeded Up Robust Features[END_REF]), or more specified adaptations (Mikolajczyk and Schmid; Rublee et al.). These features are then matched using the RANSAC algorithm (Fischler and Bolles) to estimate the optimal transformation, as it can be seen in many biomedical imaging works [START_REF] Allaire | Full orientation invariance and improved feature selectivity of 3D SIFT with application to medical image analysis[END_REF][START_REF] Paganelli | Scale Invariant Feature Transform as feature tracking method in 4D imaging: a feasibility study[END_REF][START_REF] Toews | Feature-based alignment of volumetric multi-modal images[END_REF]. However, these methods rely on strong similarities between each modalities which can be limited in a multimodal context. This problem can also be solved using variational approaches. In this case, the optimal alignment can be defined as the maximum of a given metric, typically Mutual Information (Viola and Wells III) or Cross-correlation [START_REF] Roshni | Using mutual information and cross correlation as metrics for registration of images[END_REF], which aim at finding correlations between two distributions of intensities. These methods perform well as long as there exists a bijection between both modalities (e.g. between CT and MR images) which is not the case between LiDAR range measurements and optical images. Another approach presented in (Sutour et al.) aligns the gradients of both modalities, thus being agnostic to any correlation between the modalities. However, this method only estimates translation and scaling without rotation. Moreover, it implies that both modalities are defined in the same domain and that gradients can be computed on both modalities, which is typically not possible when dealing with 3D points.

LiDAR to optical registration

The problem of LiDAR to optical registration can be divided into three main kinds of approaches: 2D feature-based methods, 3D-based method and statistical methods.

2D feature-based methods aim at establishing correspondences between feature points of the optical image and the point cloud projected in the optical image domain. In (Moussa et al.), the authors propose a method that uses ASIFT features (Morel and Yu) to match a colorized point cloud with an optical image. Aberrant correspondences are then filtered out using RANSAC (Fischler and Bolles). The final 3D pose is estimated by solving a Perspective-n-Point problem (Lepetit et al.) in which the 2D coordinates of feature points in the optical image is associated with the 3D locations of the corresponding feature points in the point cloud. González et al. propose a method for estimating the location of an optical image relatively to a 3D colorized point cloud of the same scene. The image is first enhanced to increase its contrasts. Then, the projection of the point cloud is manually resized in order to fit the optical image as well as possible. After that, correspondences are estimated by averaging cross-correlation and least square metrics. Finally, the 3D pose is retrieved using RANSAC. This method assumes that the original image and the point cloud are acquired at very close location otherwise the distortion brought by the resizing method would affect the correspondence finding step. Although 2D feature-based methods provided straight forward ways to estimate the optimal alignment between optical image and point cloud, they typically rely on shared information between the two modalities. This can be a major drawback on LiDAR systems are primarily designed to collect range rather than spectral measurements.

3D-based methods offer to align the 3D LiDAR point cloud with the 3D reconstruction of a set of optical images. Corsini et al. propose a two-step method for 3D-based point cloud to image alignment. First, a 3D sparse point cloud is reconstructed from a set of input optical images by using Structure From Motion (SFM) algorithm. The SFM algorithm is designed to find 2D correspondences in images of an input set of images and to regress the 3D pose of each image as well as the 3D position of each feature point, producing a sparse point cloud. After that, the 4-points congruent set (Aiger et al.) algorithm is used to align the sparse 3D point cloud with the 3D LiDAR point cloud. Later, Abayowa et al. propose a similar method for aligning a 3D Li-DAR point cloud with a set of aerial optical images. A dense 3D point cloud model is built from the set of optical images using the dense 3D reconstruction method described by Furukawa and Ponce. Then, the pose of the dense point cloud is recovered by using Iterative Closest Point (ICP) (Besl and McKay) algorithm in order to minimize the distance error between the dense point cloud and the LiDAR point cloud. Although these methods achieve high quality results, they require a set of input optical images instead of a single image. Moreover, 3D registration methods are largely sensitive to missing data that often appear in real urban LiDAR data.

Statistical methods for point cloud to image registration try to define metrics that can be used to measure similarities between the two input modalities. Most of the time, the metric is computed in the 2D image domain. The work described in Miled et al. proposes to align the sparse projection of a LiDAR point cloud with an optical image by comparing both modalities using Mutual Information (MI). In Miled et al., this metric is used to find the dependency between the colors carried by the optical images and the reflectances brought by the LiDAR point cloud. The pose between the image and the point cloud is computed using a variational model that maximizes the MI metric between the two modalities. This method achieves very convincing results. However it strongly relies on the quality of the reflectances aquired by the LiDAR sensor. In practical use, only very few high quality LiDAR sensors can reach such levels of accuracy. Most common sensors acquire reflectance with high level of noise. Moreover, the reflectance is only relevant in certain scenarios and cannot be used on wet surfaces or highly reflective surfaces for example. To overcome the problem of using reflectance, a method for the registration of a raw LiDAR point cloud with a single image is proposed in [START_REF] Castorena | Autocalibration of lidar and optical cameras via edge alignment[END_REF]. There, the authors propose to fuse and align the modalities at the same time by computing a dense image from the projection of the point cloud and by aligning edges of both modalities. This method can only perform well if the acquisition center of both modalities are very close. Otherwise a lot of ambiguities can arise from the LiDAR projection in the image domain which often leads to large errors in the calibration estimation. Later, (Guislain et al.) proposed a method that aims at aligning only visible points of the LiDAR point cloud with the optical image. To do so, they first estimate the visible points given the optical image point of view using [START_REF] Rubinstein | Improved seam carving for video retargeting[END_REF]. The remaining points are used to produce a dense image of reflectances by performing bilinear inpainting. This dense reflectance image is aligned with the optical image using a metric that is less sensitive to missing data than Mutual Information. In the case where the reflectance is not available, they offer to compute the same metric on a dense normals map of visible points. This method achieves very good results when the visibility estimation performs well. This is the case when each different objects of the 3D scene are well separated. However, in the case of urban scenes, the amount of missing data as well as the heterogeneity of the shapes and object is very challenging for visibility estimation methods as shown in [START_REF] Biasutti | Visibility Estimation in Point Clouds with Variable Density[END_REF]. Therefore, the quality of the results on real urban data often lacks of accuracy.

In this paper, we only consider the case of fine alignment between a LiDAR point cloud and an optical image. Thus, we estimate that the deformation induced by the perspective between the two modalities can be ignored as the original alignment is close to the optimal one. Therefore, it is sufficient to only estimate a 2D transformation between the optical image and the point cloud in the image domain, which largely simplifies the problem, but limits its usage to close initializations. In the next section, we propose a novel method for LiDAR point cloud and optical image alignment that uses the topology of the LiDAR sensor to generate a dense image without any visibility ambiguities. This dense image is later aligned with the optical image using a variational model.

Methodology

In this section, we present each step of the proposed framework for point cloud to image registration. The proposed framework is highlighted Figure 1: first, an image is created by rendering the triangulation based on the sensor topology of the point cloud. Then, this rendering is aligned with the optical image using a variational approach to align the gradients of both modalities.

Fast mesh reconstruction in sensor topology

The first step of the proposed framework consists in the reconstruction of the mesh of the point cloud. The problem of mesh reconstruction consists in linking points of a point cloud with triangles in order to approximate the surface of the objects in the scene. Surface reconstruction is traditionally done by smoothness approaches [START_REF] Lipman | Parameterization-free projection for geometry reconstruction[END_REF][START_REF] Xiong | Robust surface reconstruction via dictionary learning[END_REF], primitive approximation [START_REF] Schnabel | Completion and reconstruction with primitive shapes[END_REF][START_REF] Lafarge | Surface reconstruction through point set structuring[END_REF] or global regularity approaches (Li et al., 2011a,b;[START_REF] Monszpart | RAPter: rebuilding man-made scenes with regular arrangements of planes[END_REF]. However, these methods are often computationally expensive. Moreover, they often require strong assumptions on the homogeneity of the point cloud, which is not suitable in the case of LiDAR acquisitions. To overcome these problems, we propose a very fast approach for mesh reconstruction that exploits sensor topology to instantly create a raw mesh from the point cloud. Note that more precise meshes can be reconstructed using the analogue method proposed in [START_REF] Guinard | Sensor-topology based simplicial complex reconstruction from mobile laser scanning IV[END_REF] but with a substantive impact on the computational time. However this work focuses on the efficiency and the performance of the final alignment between LiDAR point cloud and optical image. Thus, the use of the method proposed in [START_REF] Guinard | Sensor-topology based simplicial complex reconstruction from mobile laser scanning IV[END_REF] is out of the scope of this work, although it would be interesting to test.

Modern LiDAR sensors often acquire 3D points following a defined pattern from which we can build a dense image [START_REF] Biasutti | Range-Image: Incorporating sensor topology for LiDAR point cloud processing[END_REF]. Indeed each point is defined by two angles and a depth, (θ, φ, d) respectively, with steps of (∆θ, ∆φ) between two consecutive positions. Each point p of the LiDAR point cloud can be mapped to the coordinates (x, y) with x = θ ∆θ , y = φ ∆φ of a 2D map, hereinafter referred to as u. An example of a point For each pixel (x, y) of u, 2 triangles ul , bt are created as follows: ul = {u(x, y), u(x + 1, y), u(x, y + 1)} br = {u(x + 1, y), u(x + 1, y + 1), u(x, y + 1)} This principle is illustrated on Figure 3. After that, triangles are filtered by discarding the ones that have at least one edge that is longer than a certain threshold t, typically t = 1.0m. This step prevents separate objects from being connected together which enhances the overall quality of the mesh. An example of reconstructed mesh is showed Figure 2(c). Finally, the mesh is being rendered from the optical camera location, with the same intrinsic parameters. This produces a dense image d u of the point cloud. As the mesh is not textured, d u is filled by the values of the zbuffer of the rendering (i.e. the depth of each pixel). Figure 4 displays an example of a sparse projection of the point cloud (b) in the image domain of (a) compared to texture-less rendering (c) and depth rendering (d). We can see that the renderings are largely denser than the sparse projection, resulting in the appearance of strong depth gradients.

Depth to optical image alignment

As mentioned in Section 2, the alignment between a LiDAR point cloud P and an optical image I is non-trivial as both modalities do not share any common attribute. The mesh rendering d u provides strong depth gradients in the image domain. These gradients correspond to object contours which can also be met in the optical image. Although strong depth gradients can occur without appearing in the optical image, and vice-versa, it is reasonable to assume that most depth gradients also appear in the optical image in real data. Therefore, aligning P and I in the domain of I can be simplified as the alignment between the gradients of d u and I. However, this assertion is only true if the initialization of the alignment between d u and I is relatively close. Indeed, the perspective induced by the 3D rendering introduces deformations that are proportional to the depth of the scene. Thus, if the initialization is too far from the optimal alignment, the alignment between the gradients of d u and I is not possible.

The method described in Sutour et al. offers to align gradients of two modalities expressed in the same image domain. To that extent, they define a variational model in which gradient alignment between images u 1 and u 2 is done by maximizing the following criterion:

C(T ) = Ω ∇u 1 (T t x ,t y ,z (X)) • ∇u 2 (X) dX, T t x ,t y ,z (X) =           1 + z 0 t x 0 1 + z t y 0 0 1           X
where Ω is the domain of definition of I and T z,t x ,t y represents a 2D affine transform with 3 degrees of freedom: vertical and horizontal translation t x , t y as well as zooming z. In the case of LiDAR point cloud to optical image alignment, rotation should also be considered in the transform as we cannot assume that the rotation between both sensors is always null. Therefore, we propose to extend the model presented in Sutour et al. in order to estimate rotation as well as translation and zooming. We define Tz,t x ,t y ,θ the 4 degrees of freedom (t x , t y translation, z zoom and θ rotation) transformation matrix such that:

Tt x ,t y ,z,θ = T t x ,t y ,z           cos θ -sin θ 0 sin θ cos θ 0 0 0 1           =           s cos θ -s sin θ t x s sin θ s cos θ t y 0 0 1          
with s = 1 + z to simplify notations. Similarly to Sutour et al., the gradients of d u and I are aligned by maximizing the following criterion:

C( T ) = Ω ∇d u ( Tt x ,t y ,z,θ (X)) • ∇I(X) dX.
Using this formulation, an explicit optimization scheme is built to maximize the proposed criterion at each iteration n, by performing a gradient ascent on each parameters of the transformation Tz,t x ,t y ,θ : where the partial derivatives of C( Tt x ,t y ,z,θ ) are defined as follows for each iteration:

                 t n+1 x = t n x + λ 1 ∂C ∂t x ( Tt x ,t y ,z,θ ) t n+1 y = t n y + λ 2 ∂C ∂t y ( Tt x ,t y ,z,θ ) z n+1 = z n + λ 3 ∂C ∂z ( Tt x ,t y ,z,θ ) θ n+1 = θ n + λ 4 ∂C ∂θ ( Tt x ,t y ,z,θ )
∂C ∂t x ( Ttx,ty,z,θ ) = Ω σ∇ 2 du (X) 1 0 • ∇I(X)dX, ∂C ∂t y ( Ttx,ty,z,θ ) = Ω σ∇ 2 du (X) 0 1 • ∇I(X)dX, ∂C ∂z ( Ttx,ty,z,θ ) = Ω σ∇ 2 du (X) x cos θ + y sin θ -x sin θ + y cos θ • ∇I(X)dX, ∂C ∂θ ( Ttx,ty,z,θ ) = Ω σ∇ 2 du (X) -x • s sin θ -y • s cos θ x • s cos θ -y • s sin θ • ∇I(X)dX
having du (X) = d u ( Tt x ,t y ,z,θ (X)) and σ = sign(∇d u (X) • ∇I(X)).

The functional we aim at optimizing is not convex. Therefore, it is highly subject to local maxima. However we consider that the alignment we seek to perform only concerns data provided by calibrated Mobile Mapping Systems. Therefore, the provided alignment of the LiDAR point clouds and the optical images is assumed to be close to the optimal alignment, as discussed here after in Section 4.1. For the gradient ascent scheme, we set λ 1 = λ 2 = 10 -3 to be larger than λ 3 = λ 4 = 10 -5 as the translation expressed in pixel is likely to be larger than the rotation or the zooming factor. We set the maximum number of iterations to 200. However, most of our experiments have shown that the method converges in less than 30 iterations on the data presented Section 4.

Finally, we propose to improve the gradient ascent scheme by refining the search steps at each iteration. The search step λ n

x at iteration n is then defined as follows:

λ n x = λ n-1 x if C n ( T ) > ρ C n-1 ( T ) λ n-1
x /2 otherwise with C n ( T ) the energy at iteration n, ρ = 0.99. This improvement prevents the algorithm from being directly stuck in a local maxima, and provides better results in practice as demonstrated in Section 4.1.

Experiments and results

We conclude this paper by presenting different results obtained using the proposed framework. The proposed pipeline is evaluated on the RobotCar dataset (Maddern et al.) which provides images of resolution 1280 × 960px as well as point clouds composed of millions of points. We demonstrate the efficiency of the proposed method through a quantitative and qualitative analysis.

Quantitative analysis

The calibration of the RobotCar dataset does not provide a perfect alignment between LiDAR point clouds and optical images. We propose to manually align mesh renderings with optical images to create ground truths. We found out that the original data alignment compared to the ground truth alignment presents a Mean Absolute Error (MAE) of about 19px for translation, 0.9 degree for rotation and 0.01 for zooming. We propose to apply comparable transformations on the manually align renders to generate evaluation data. The transformations are generated by randomly and uniformly shifting the renders between -20 and 20 pixels on both x and y axis, rotating the renders between -1 and 1 degree and zooming by a uniform random factor between 0.95 and 1.05.

We compare our method with and without the refinement of the search steps to the method proposed in (Sutour et al.), as this method presents the baseline of gradient alignment without the estimation of the rotation. Moreover, we also compare our method to an exhaustive search of the maximum of the Mutual Information (Viola and Wells III) as done in recent multimodal alignment methods, such as (Miled et al.). We compute the MAE between each estimated parameter (t x , t y , z, θ) and the ground truth. The results of this experiment are summarized in Table 1. We can see that our method achieves very fine alignment of LiDAR point cloud and optical image. The method with refinement of search steps provides finer results than each other method. The use of the functional defined in (Sutour et al.) as well as the extension presented in this paper outperforms the exhaustive search with Mutual Information metric.

Moreover, we can see that extending the original functional by adding the regression of rotation improves the results not only in the estimation of the rotation, but also in the overall alignment. This is due to the fact that limiting the transformations to translation and scaling prevents the algorithm from finding the optimal alignment. Therefore, the baseline algorithm finds another local maxima which does not align well both modalities. This shows the importance of predicting the rotation as well as the baseline parameters of the transformation. Finally, the refinement of the search steps prevents the variational model from being stuck in local maxima, which makes it more robust to largely shifted initialization while keeping the same computational cost.

Qualitative analysis

We conclude our experiments with a qualitative analysis. Figure 5 presents the results of LiDAR point cloud to optical image alignment using our method. The first row shows the original alignment, the second row shows the results of the alignments using our method, with closeup looks at the original alignments and our results on the last two rows respectively. On each image, the strong gradients of the depth renderings are represented by green lines on the optical images.

The results presented in Figure 5 highlight that our model succeeds in aligning gradients of both modalities, producing a very good 2D registration between LiDAR point clouds and optical images. From initialization with shifted alignments (shown in the first row), our method produces results where both modalities are seamlessly aligned (second row). In particular, the last row of Figure 5 shows some areas where the variational model perfectly matches the renders and the optical images on structures that display strong gradients such as roof lines or windows. Moreover, the method only requires to match a small amount of gradients in order to correctly align both modalities. This property makes it more robust to outliers as some gradients of the depth rendering do not correspond to any gradient in the optical image, and vice-versa, as discussed previously in Section 3.2. Finally, our method is able to produce good alignment even when initialized with large shifts between both modalities. This is specially visible in the last column where we can see that in the original alignment, the optical image is shifted from the gradients of the depth render. Despite this initialization, our method succeeds in producing a very fine alignment of the two modalities as it can be seen on the lowest line.

Conclusion

In this paper, we have proposed a novel framework for LiDAR point clouds to optical images alignment. The first step of this framework offers to reconstruct the mesh from the point cloud by exploiting the topology of the sensor. After that, the mesh is rendered with the same pose as the optical image. Finally, the depth gradients of the rendered LiDAR mesh and the color gradients of the optical image are aligned using a modified variational approach from Sutour et al.. The qualitative and quantitative results demonstrate that the framework succeeds in very fine alignment between both modalities.

Although the overall results of the proposed method are satisfying, it depends on the initialization of the alignment. In the future, we would like to extend the variational model to perform the gradient alignment from coarse to fine scale to make it less dependant on the initialization. We also would like to compare our current results to the one obtained on renderings done with the mesh produced by the method presented in [START_REF] Guinard | Sensor-topology based simplicial complex reconstruction from mobile laser scanning IV[END_REF].

Acknowledgement

This project has also received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No 777826.

Figure 1 :

 1 Figure 1: Scheme of the proposed framework.

Figure 2 :

 2 Figure 2: Mesh reconstruction scheme. (a) is the input point cloud, (b) the point cloud as seen in sensor topology and (c) the reconstructed mesh.

Figure 3 :Figure 4 :

 34 Figure 3: Triangle construction from image in sensor topology.

Figure 5 :

 5 Figure5: Example of alignments produced by our method. The green lines correspond to the strong gradients of the depth render. The first row shows the original alignment between the optical images and the mesh gradients. The second row shows the alignment produced by our method. A closeup look at details of original alignment and our results is showed on the last two rows respectively.

Table 1 :

 1 MAE of each method compared to the manually aligned data for each parameter on 50 randomly generated transformations.

	Method	Mean Absolute Error t x t y z	θ
	Mutual Information	16.3 11.9	0.05	0.46
	Sutour et al. (baseline)	2.91 6.76 0.006 0.57
	baseline + rotation	2.96 6.29 0.004 0.04
	baseline + rotation + refined 1.93 3.31 0.005 0.03


Xi Deng

Pierre Boivin

Feng Xiao

A new formulation for two-wave Riemann solver accurate at contact interfaces

This study proposes a new formulation for Harten, Lax and van Leer (HLL) type Riemann solver which is capable of solving contact discontinuities accurately but with robustness for strong shock. It is well known that the original HLL, which has incomplete wave structures, is too dissipative to capture contact discontinuities accurately. On the other side, contact-capturing approximate Riemann solvers such as Harten, Lax and van Leer with Contact (HLLC) usually suffer from spurious solutions, also called carbuncle phenomenon, for strong shock. In this work a new accurate and robust HLL-type formulation, so-called HLL-BVD (HLL Riemann solver with BVD) is proposed by modifying the original HLL with BVD (boundary variation diminshing) algorithm. Instead of explicitly recovering the complete wave structures like the way of HLLC, the proposed method restores the missing contact with a jump-like function. The capability of solving contact discontinuities is further improved by minimizing the inherent dissipation term in HLL. Without modifying the original incomplete wave structures of HLL, the robustness for strong shock has been reserved. Thus the proposed method is free from shock instability problem. The accuracy and robustness of the new method are demonstrated through solving several one-and two-dimensional tests. Results indicate that the new formulation based on two-wave HLL-type Riemann solver is not only capable of capturing contact waves more accurately than the original HLL or HLLC but, most importantly, is free form carbuncle instability for strong shock.

I. INTRODUCTION

The modern finite volume method (FVM) for hyperbolic systems is developed based on Godunov scheme 1 which assumes that (1) the flow variables within each computing elements are constant and (2) numerical fluxes are calculated from the solution of the Riemann problem. Owing to its robustness and flexibility, Godunov approach has been further developed for decades, and is popularly applied to solve Euler equations for high speed compressible flows. One branch in the development is the construction of approximate Riemann solvers, or in another word, upwind fluxes. Different from the exact Riemann solver, the approximate algorithms are noniterative and hence less expensive. Over the years, several upwind fluxes have been proposed, which can be generally categorized into flux vector splitting (FVS) and flux difference splitting (FDS) methods [START_REF] Leveque | Finite volume methods for hyperbolic problems[END_REF] .

The FVS methods split flux vectors into upstream and downstream travelling informations. These methods are generally less complicated and more easily compatible with implicit methods. However, they are not effective in capturing the discontinuities represented by linear waves. Thus they are prone to intensively smear out contact surfaces and shear waves. Schemes such as the Steger-Warming 3 and the FVS scheme of van Leer 4 belong to a) Electronic mail: deng.xi98@gmail.com; Also at Department of Mechanical Engineering, Tokyo Institute of Technology b) Electronic mail: pierre.boivin@univ-amu.fr c) Electronic mail: xiao.f.aa@m.titech.ac.jp this approach.

On the other hand, the FDS methods are constructed from approximate solution of local Riemann problem between two adjacent states. The FDS methods can be further divided into two groups as complete or incomplete Riemann solvers. The complete solvers have the same wave structures as the exact Riemann solver. For Euler equations, the wave structures comprise a contact discontinuity and two nonlinear waves which corresponds to either a shock or a rarefaction wave. By virtue of comprising all wave components, the complete FDS methods are capable of capturing linear as well as non-linear waves. This attractive quality makes complete Riemann solvers such as the Roe solver 5 , the Osher solver [START_REF] Osher | Upwind difference schemes for hyperbolic systems of conservation laws[END_REF] and Harten, Lax and van Leer with Contact (HLLC) solver widely used. Incomplete Riemann solvers, however, employ coarser wave structures which do not include the contact discontinuity. Thus, they generally bring excessive dissipation across contact. The incomplete solvers include the HLL [START_REF] Harten | On upstream differencing and godunov-type schemes for hyperbolic conservation laws[END_REF] , the Harten Lax van Leer Einfeldt (HLLE) [START_REF] Einfeldt | On godunov-type methods for gas dynamics[END_REF] and the Rusanov solver 9 . Among these approximated Riemann solvers introduced above, HLLC method, which is less dissipative, is considered to be one of the most preferred FDS method because it can capture contact very accurately. Unfortunately, it may encounter a problem of numerical instability known as the carbuncle phenomenon. In fact, it has been known that Riemann solvers, which can resolve the contact and shear waves exactly, usually suffer from numerical instability in the vicinity of strong shock waves [START_REF] Kim | Robust hllc riemann solver with weighted average flux scheme for strong shock[END_REF] . On the other side, HLL, which is more dissipative due to its incomplete wave structures, is very robust.

It has an entropy satisfaction property, exactly resolves isolated shock and preserves positivity [START_REF] Harten | On upstream differencing and godunov-type schemes for hyperbolic conservation laws[END_REF][START_REF] Toro | Riemann solvers and numerical methods for fluid dynamics: a practical introduction[END_REF][START_REF] Einfeldt | On godunovtype methods near low densities[END_REF] although it is notorious that it cannot resolve contact discontinuities exactly.

Consequently, there have been several attempts to construct Riemann solvers which can overcome the problem of either shock instability or exactly resolving contact discontinuities. Among these attempts, the most popular methodology is a hybrid method proposed in [START_REF] Quirk | A contribution to the great riemann solver debate[END_REF] . It is suggested that using a dissipative scheme in shock region while a less dissipative one elsewhere may tackle this problem. From then, several successful approaches [START_REF] Kim | Robust hllc riemann solver with weighted average flux scheme for strong shock[END_REF][START_REF] Huang | Cures for numerical shock instability in hllc solver[END_REF][START_REF] Wu | A hybrid numerical method to cure numerical shock instability[END_REF][START_REF] Nishikawa | Very simple, carbuncle-free, boundary-layer-resolving, rotated-hybrid riemann solvers[END_REF] have been proposed by relying on hybridization. However, the choice and design of switching parameters required in the hybrid methods usually have some arbitrariness. This arbitrariness in parameters may lead to failures for complex flow cases [START_REF] Shen | A robust hllc-type riemann solver for strong shock[END_REF] . Moreover, the hybridization of complementary solvers may break the strict Rankine-Hugoniot jump conditions, which will result in smearing contact discontinuities [START_REF] Wu | A hybrid numerical method to cure numerical shock instability[END_REF] . Realizing the drawbacks of hybrid methods, other ideas have been tried to modify the HLL-type approximate Riemann solver. For example, in the work of [START_REF] Mandal | Robust hll-type riemann solver capable of resolving contact discontinuity[END_REF] Euler equations are split into convective and pressure parts before formulating the HLL in order to recover the contact in HLL. Through modifying the dissipation term in HLL flux, the ability to resolve the contact is further enhanced. In 17 , new HLLtype and HLLC-type Riemann solvers are constructed by modifying the shear viscosity of the original HLL and HLLC methods respectively. Although these modifications give good results in two dimensional problems, the performance has not been evaluated in one dimensional cases. Recently, several methods have been proposed to reduce shock instability for complete wave Riemann solvers through controlling numerical dissipation. In the work [START_REF] Rodionov | Artificial viscosity in godunov-type schemes to cure the carbuncle phenomenon[END_REF] , artificial viscosity has been introduced in the control equations to cure the carbuncle phenomenon. Shockstable items or shear viscosity are introduced into numerical flux in [START_REF] Chen | Affordable shock-stable item for godunov-type schemes against carbuncle phenomenon[END_REF][START_REF] Chen | A new robust carbuncle-free roe scheme for strong shock[END_REF] to stabilize strong shock. A method named normal reconstruction procedure is proposed in [START_REF] Chen | General procedure for riemann solver to eliminate carbuncle and shock instability[END_REF] to add dissipation to transverse velocity which tends to cause propagation of shock anomalies. In [START_REF] Simon | A cure for numerical shock instability in hllc riemann solver using antidiffusion control[END_REF][START_REF] Simon | A simple cure for numerical shock instability in the hllc riemann solver[END_REF] , shock stability is realized by effective controlling dissipation term. Based on the fact that different reconstruction processes have different dissipation properties, the work [START_REF] Jiang | Effective technique to improve shock anomalies and heating prediction for hypersonic flows[END_REF] alleviates shock anomalies through combing with reconstruction schemes.

Despite these efforts, there is still no consensus on how to construct an accurate and robust approximate Riemann solver. In the present work, a simple novel method, so-called HLL-BVD (HLL Riemann solver with boundary variation diminishing algorithm), is proposed to solve contact discontinuities accurately without suffering from numerical shock instability problems. By realizing that the dissipation term in HLL flux leads to the failure of accurately resolving contact, the current method makes an attempt to minimize this dissipation term through introduction of a jump-like function. In line with [START_REF] Sun | Boundary variation diminishing (bvd) reconstruction: A new approach to improve godunov schemes[END_REF] , which uses Tangent of Hyperbola for INterface Capturing (THINC) function [START_REF] Xiao | Revisit to the thinc scheme: a simple algebraic vof algorithm[END_REF][START_REF] Xiao | A simple algebraic interface capturing scheme using hyperbolic tangent function[END_REF] with boundary variation diminishing (BVD) algorithm to effectively control numerical dissipation, the current work employs a jump-like THINC function to represent the information of contact discontinuities and derives a consistent BVD algorithm to minimize the dissipation term. This leads to significant improvements in the resolution of contact discontinuities. Since the new method restores the contact with an implicit algebraic way, the incomplete wave structure of HLL is preserved. Thus, unlike other methods which retrieve contact by modifying the wave structures of HLL, the robustness of HLL for strong shocks can be kept. Examined through several benchmark tests, the proposed HLL-BVD approach is shown to better represent contact interfaces and be free from shock instability problems.

The outline of the paper is as follows. In Section 2, after a brief review of HLL-type Riemann solvers, details of our HLL-BVD method are presented. In Section 3, numerical results of benchmark tests are presented and systematically compared with conventional HLL and HLLC methods. Some concluding remarks are given in Section 4.

II. NUMERICAL METHODS

A. Governing equations

The one-dimensional Euler equations can be expressed as

∂U ∂t + ∂F(U) ∂x = 0 (1)
under conservative form, where the physical variables vector U and flux vectors F(U) read

U =    ρ ρu ρv E    and F(U) =    ρu ρu 2 + p ρuv u(E + p)    .
Following the conventional notation, ρ denotes the density, p the pressure and E the total energy per unit volume. u is the fluid velocity in the x direction and v is the transverse velocity. System (1) is closed through introduction of the perfect gas equation of state

p = (γ -1)(E - 1 2 ρ(u 2 + v 2 )), ( 2
)
where γ is the specific heat ratio. Within the finite volume framework, System (1) in semi-discrete form reads

∂ Ūi ∂t = - 1 ∆x (F i+ 1 2 -F i-1 2 ), (3)
where the volume-integrated average value Ūi in cell I i is defined as

Ūi ≈ 1 ∆x ˆxi+1/2 x i-1/2 U(x, t) dx. (4)
The numerical fluxes F i± 1 2 at cell boundaries can be computed by a Riemann solver. In the following, let the subscripts L and R identify the left and right values of a given quantity on either side of the cell interface.

B. HLL Riemann solvers

Within the HLL framework, each cell boundary is regarded as a two-wave system. The system, illustrated in Fig. 1, then consists of two waves travelling at (S L , S R ).

In that context, writing the conservation equations (1) in weak form

˛U.dx -F.dt = 0, (5)
and assuming U to be constant on each branch of the integration path illustrated in Fig. 1 directly leads to

U * = S R U R -S L U L + F L -F R S R -S L , (6)
where the * superscript denotes the middle state. The corresponding flux can be found as

F * = S R F L -S L F R + S L S R (U R -U L ) S R -S L . (7)
The final flux F i+ 1 2 , required to solve (3), then read

F hll i+ 1 2 =    F L if 0 ≤ S L , F * if S L ≤ 0 ≤ S R , F R if 0 ≥ S R . (8)
Following [START_REF] Davis | Simplified second-order godunov-type methods[END_REF] , the right and left wave velocities are approximated as

S L = min(u L -c L , u R -c R ), S R = max(u L +c L , u R +c R ), (9)
where c L/R is the sound speed on either side of the interface.

As with most shock capturing schemes, robustness is assured by inclusion of diffusion at the shock. This can be easily seen by separating the flux (7) into two contributions

F hll = 1 2 (F L + F R ) + ε hll , (10)
where the first term corresponds to the classical central difference flux, and ε hll can be interpreted as a numerical diffusion term

ε hll = S R + S L 2(S R -S L ) (F L -F R ) ε 1 - S L S R S R -S L (U L -U R ) ε 2 . ( 11
) x t S R S L U R , F R U L , F L U * , F * FIG.
1. The wave structure in HLL Riemann solvers. In red, the integration contour leading to Eq. (6) .

When the flow presents density discontinuities, the density (and therefore U) is not constant between the S L and S R waves. This leads to excessive numerical diffusion and the failure of HLL Riemann solver to accurately solve contact discontinuities [START_REF] Mandal | Robust hll-type riemann solver capable of resolving contact discontinuity[END_REF] .

C. Towards the reduction of numerical diffusion

There are several strategies to cope with these deficiencies. One strategy, proposed in 18 is to assume local isentropy, as to replace the density jump terms in the second term of the HLL diffusion ε hll (11)

ε 2 = S L S R S R -S L    ρ L -ρ R (ρu) L -(ρu) R (ρv) L -(ρv) R E L -E R    , (12)
into pressure jump terms

ε isen 2 = S L S R c2 (S R -S L )     p L -p R (pu) L -(pu) R (pv) L -(pv) R c2 γ-1 (p L -p R ) + 1 2 ((p(u 2 + v 2 )) L -(p(u 2 + v 2 )) R )     , (13)
through introduction of an interface sound speed c = cL+cR 2

. This method can improve the resolution across contact waves and keep shock stable. However, the isentropic hypothesis may not always hold.

A second strategy is the introduction of an antidiffusive term in (11). In that direction, the HLLEM solver 8 is a possible alternative, with the diffusive term now reading

ε hllem = S R + S L 2(S R -S L ) (F L -F R ) - S L S R S R -S L (U L -U R ) + φ S L S R S R -S L M (U)(U L -U R ), (14)
where the scalar φ ∈ [0, 1] is a free parameter, and M (U) arises from the eigen-structure of Jacobi matrix of hyperbolic system [START_REF] Balsara | Self-adjusting, positivity preserving high order schemes for hydrodynamics and magnetohydrodynamics[END_REF][START_REF] Dumbser | A new efficient formulation of the hllem riemann solver for general conservative and nonconservative hyperbolic systems[END_REF] .

Clearly, the HLLEM Riemann solver reduces to the HLL solver for φ = 0. Just as the HLLC solver, the HLLEM solver is flexible and able to resolve contact waves with accuracy. Unfortunately, HLLC and HLLEM solvers share similar shock instability problems, believed to be caused by perturbation transverse to the shock waves [START_REF] Flandrin | An improved approximate riemann solver for hypersonic bidimensional flows[END_REF] . These comments are compatible with the observation [START_REF] Pandolfi | Numerical instabilities in upwind methods: analysis and cures for the "carbuncle" phenomenon[END_REF] that Riemann solvers explicitly dealing with the contact surface display clear evidence of carbuncle phenomenon.

This motivates the derivation of an implicit algebraic method to capture the missing contact waves in HLL Riemann solver without modifying the original two-wave system.

D. HLL-BVD: A new formulation for HLL accurate at contact

As stated in the previous section, the failure of HLL Riemann solver to resolve contact waves is caused by the numerical dissipation term ε hll (11). Assuming constant flow variables in each cell, in particular, leads to excessive numerical dissipation due to the density jump |ρ L -ρ R |

In this work, we then suggest a procedure to minimize that contribution to the numerical dissipation ε hll . Inspired by [START_REF] Sun | Boundary variation diminishing (bvd) reconstruction: A new approach to improve godunov schemes[END_REF][START_REF] Deng | Limiter-free discontinuity-capturing scheme for compressible gas dynamics with reactive fronts[END_REF][START_REF] Deng | High fidelity discontinuity-resolving reconstruction for compressible multiphase flows with moving interfaces[END_REF] , we use the THINC function [START_REF] Xiao | Revisit to the thinc scheme: a simple algebraic vof algorithm[END_REF][START_REF] Xiao | A simple algebraic interface capturing scheme using hyperbolic tangent function[END_REF] to reconstruct density jumps

ρi (x) = ρmin + ∆ρ 2 1 + θ tanh β x -x i-1/2 x i+1/2 -x i-1/2 -xi , (15)
where ρmin = min(ρ i-1 , ρi+1 ), ∆ρ = |ρ i+1 -ρi-1 |, and θ = sgn(ρ i+1 -ρi-1 ). The jump "thickness" is controlled by parameter β, set to 1.6 throughout the validation section. It is noted that from our numerical experiments, the parameter β valued from 1.4 to 2.0 is also able to give good or acceptable results. A smaller β will results in more diffusive solutions while a larger one tends to introduce extra anti-diffusion effect. xi represents the position of jumps, computed as to satisfy ρi = 1 ∆x ´xi+1/2

x i-1/2 ρi (x) dx. Having reconstructed the THINC function in every cell, it becomes straightforward to retrieve the corresponding density values ρ T L/R on either side of the face (note the T superscript for THINC reconstruction).

On the other hand, we have at our disposal the interface values ρ P L/R , coming from the Polynomial reconstruction, which may vary depending on the numerical scheme. In the present work, we use MUSCL (Monotone Upstream-centered Schemes for Conservation Law) scheme [START_REF] Van Leer | Towards the ultimate conservative difference scheme. v. a second-order sequel to godunov's method[END_REF] as the polynomial reconstruction.

In order to minimize the dissipation term, a BVD algorithm which is originally proposed in [START_REF] Sun | Boundary variation diminishing (bvd) reconstruction: A new approach to improve godunov schemes[END_REF] is devised to minimize the density jumps on each cell boundaries. Final values for ρ L/R are picked from ρ P L/R and ρ T L/R values, as to minimize the density jump. In other words,

|ρ L -ρ R | = min(|ρ T L -ρ T R |, |ρ P L -ρ P R |, |ρ P L -ρ T R |, |ρ T L -ρ P R |), (16)
leading to a minimal density jump in the dissipative term in ε hll (11), now reading

ε new = S R + S L 2(S R -S L ) (F L -F R ) - S L S R S R -S L      ρ L -ρ R ρ L u L -ρ R u R ρ L v L -ρ R v R ( p L γ -1 + 1 2 ρ L (u 2 L + v 2 L )) -( p R γ -1 + 1 2 ρ R (u 2 R + v 2 R ))      . ( 17
)
To summarize, the numerical diffusion term of the new formulation for HLL-type Riemann solver is calculated by following steps:

1. Prepare the density values ρ P L/R at cell interfaces by using classical Godnunov scheme or high order schemes;

2. Retrieve the discontinuous density field in the intermediate state with THINC function of (15) to get ρ T L/R ;

3. Determine the intermediate density value ρ L and ρ R from ρ T L/R and ρ P L/R with BVD algorithm in (16) to minimize the density jump |ρ L -ρ R | at cell boundaries;

4. The final flux is then computed through (10) and (17).

Since the contact waves are restored through an algebraic way without modifying the wave structures explicitly like HLLC or HLLEM, the numerical tests in following section will show that the new HLL-type solver can accurately solve sharp contact whilst preserving shock stability.

III. NUMERICAL RESULTS

This section presents a series of numerical tests conducted to evaluate the performance of the proposed new HLL-BVD Riemann solver and systematic comparison with the HLL and HLLC solvers are presented. Throughout this section, the specific heat ratio is set as γ = 1.4. The CFL number is set to 0.4 unless specifically noted.

A. Isolated contact

Let us first consider contact wave propagation at low Mach numbers. A one-dimensional domain is initialized as

(ρ 0 , u 0 , p 0 ) = (1.4, 0.1, 1.0) 0 ≤ x ≤ 0.5 (1.0, 0.1, 1.0) otherwise , (18)
corresponding to a contact discontinuity moving at M a = 0.1. The computation is carried out with 100 uniform cells and with time integrated up to t = 2.0. Figure 2.a presents the density profile as obtained through integration with the classical HLL and HLLC solvers, as well as the present HLL-BVD solver. It can be seen that the moving contact wave resolved by the original HLL Riemann solver is the most diffusive. Although HLLC can resolve the moving contact wave better than HLL, it still poorly represents the contact wave as a sharp discontinuity. However, the proposed HLL-BVD has the ability to capture the sharp jump of density field across the contact.

The test shows that the proposed HLL-BVD is superior to original HLL and HLLC regarding to represent sharp contact discontinuities.

To further examine the performance of the proposed formulation, the same test case is carried out at a very low Mach Number, M a = 10 -5 . The results, reported in Fig. 2.b, exhibit the failure of the HLL solver, consistently with previous reports [START_REF] Guillard | On the behaviour of upwind schemes in the low mach number limit[END_REF][START_REF] Thornber | On entropy generation and dissipation of kinetic energy in highresolution shock-capturing schemes[END_REF][START_REF] Rieper | On the dissipation mechanism of upwind-schemes in the low mach number regime: A comparison between roe and hll[END_REF] . Our HLL-BVD proposal, however, successfully capture the contact discontinuity, even at very low Mach number.

B. Sod's problem

Our second test case is the Sod problem, with initial conditions set as [START_REF] Sod | A survey of several finite difference methods for systems of nonlinear hyperbolic conservation laws[END_REF] (ρ 0 , u 0 , p 0 ) = (1, 0, 1) 0 ≤ x ≤ 0.5 (0.125, 0, 0.1) otherwise .

This initial condition will result in a right running shock, a right running contact wave and a left running rarefaction wave. With the same number of cells as in previous case, the computation lasts until time t = 0.25. The density and pressure profiles, displayed in Fig 3, are very similar whether the HLL, HLLC or HLL-BVD solver is used, with the exception of the contact discontinuity, slightly thinner in the HLL-BVD case.

C. Lax's problem

The Lax problem is employed here to check the performance of the proposed numerical scheme to capture relatively strong shock [START_REF] Shu | Efficient implementation of essentially nonoscillatory shock-capturing schemes[END_REF] . The initial condition is (ρ 0 , u 0 , p 0 ) = (0.445, 0.698, 3.528) 0 ≤ x ≤ 0.5 (0.5, 0, 0.571) otherwise .

(20) A relatively strong right running shock, a right running contact discontinuity and a left running expansion wave will evolve with this initial condition. Numerical solutions for density and pressure field are presented in Fig. 4. Again, solutions for expansion fan and shock are almost identical for all three types of HLL Riemann solvers. Again, HLL-BVD solves the contact discontinuity with slightly less numerical diffusion.

Strong Lax problem

In this test, a numerical test involving a strong shock is considered here to examine the ability of numerical methods to capture contact discontinuities as well as strong shock waves. The initial condition is set as following:

(ρ 0 , u 0 , p 0 ) = (1.0, -19.59745, 1000.0) 0 ≤ x ≤ 0.8 (1.0, -19.59745, 0.01) otherwise .

(21) The solution consists of a strong shock, a stationary contact and a left running expansion wave. The numerical results at t = 0.012 from different Riemann solvers are presented in Fig. 5 for density field and Fig. 6 for pressure field. For this case, it seems that the HLLC Riemann solver can resolve stationary contact with the least amount of numerical diffusion. However, it produces the largest overshoot values across the shock waves as shown in zoomed region of Fig. 5 and Fig. 6, which may lead unstable solutions for strong shock. However, the proposed one has the ability to represent stationary contact with much less numerical diffusion than the original HLL but also capture strong shock waves without overshooting.

D. Slow moving shock

To further investigate the abilities of different HLLtype Riemann solvers to capture shock waves stably, the slowly moving shock problem is considered here. The computational domain is still one unit long, but now consists of 1000 elements. The initial condition of slowly right moving shock is set as follows (ρ 0 , u 0 , p 0 ) = (1.0, -3.44, 1.0) 0 ≤ x ≤ 0.5 (3.86, -0.81, 10.33) otherwise .

(

) 22
The post-shock condition is set at the left boundary while the pre-shock condition is set at the right side. The numerical solutions are obtained at t = 1.5 and shown in Fig. 7. Among the three solvers, HLLC leads to the maximum oscillation amplitude in the post-shock state, while HLL-BVD clearly leads to the most stable solution. As shown in two-dimensional tests later, the unstable shock waves from HLLC may result in unacceptable solutions.

E. 2D Riemann problems

Let us consider two-dimensional problem, starting with the two-dimensional Riemann problem proposed and extensively studied in [START_REF] Schulz-Rinne | Classification of the riemann problem for two-dimensional gas dynamics[END_REF][START_REF] Kurganov | Solution of two-dimensional riemann problems for gas dynamics without riemann problem solvers[END_REF] . It has been reported in [START_REF] Dumbser | A posteriori subcell limiting of the discontinuous galerkin finite element method for hyperbolic conservation laws[END_REF] that sharp capturing of discontinuities as well as shear waves is critical to capture the small-scale structures created by the Kelvin-Helmholtz instability.

The computational domain is [-0.5, 0.5] × [-0. Case 1:

(ρ 0 , u 0 , v 0 , p 0 ) =      0.75, 0.5, 1.0)
x ≤ 0.0, y ≥ 0.0 (1.0, -0.75, 0.5, 1.0) x < 0.0, y < 0.0 (1.0, 0.75, -0.5, 1.0) x > 0.0, y > 0.0 (3.0, -0.75, -0.5, 1.0) x > 0.0, y < 0.0

. ( 23
) Case 2: (ρ 0 , u 0 , v 0 , p 0 ) =     
(2.0, -0.75, 0.5, 1.0) x ≤ 0.0, y ≥ 0.0 (1.0, 0.75, 0.5, 1.0) x < 0.0, y < 0.0 (1.0, -0.75, -0.5, 1.0) x > 0.0, y > 0.0 (3.0, 0.75, -0.5, 1.0) x > 0.0, y < 0.0

. ( 24
)
Uniform grids of 600 × 600 are employed in our calculation. The numerical results calculated from different HLL-type Riemann solvers are presented in Fig. 89. These results clearly indicate the superiority of our HLL-BVD solver over the classical HLL and HLLC options in capturing contact and shear waves. The solutions produced by the original HLL Riemann solver are the most diffusive, which will lead to failure of resolving Kelvin-Helmholtz instability. The tests in this subsection demonstrate that the proposed HLL-BVD can work well in two dimensional problems whilst significantly improving the capability to solve contact discontinuities and associated fine flow structures.

F. Odd-even perturbation problem

It is well known that Riemann solvers which can accurately capture contact discontinuities usually suffer from spurious solutions, called carbuncle instabilities, in the vicinity of strong shock. Solving strong shock is especially important for problems involving reaction and heat release [START_REF] Deng | High resolution multi-moment finite volume method for supersonic combustion on unstructured grids[END_REF][START_REF] Deng | New accurate and efficient method for stiff detonation capturing[END_REF][START_REF] Yang | Initiation structure of oblique detonation waves behind conical shocks[END_REF][START_REF] Zhang | Effects of activation energy on the instability of oblique detonation surfaces with a one-step chemistry model[END_REF] . In this case, the odd-even decoupling test originally proposed in [START_REF] Quirk | A contribution to the great riemann solver debate[END_REF] is employed to examine the ability of the proposed HLL-BVD to obtain the stability of strong shocks. The computational domain consists of 51 × 26 elements. Initially, a steady normal shock of Mach=20 is set in the domain with initial conditions for left (i ≤ 39) and right (i ≥ 41) sides following Rankine-Hugoniot conditions. An intermediate state is prescribed at i = 40 where the internal shock conditions are given as ρ M = λρ L + (1 -λ)ρ R with other variables calculated based on ρ M . A perturbation of shock along y-direction is produced by assigning different λ as

λ = 0.2 if j is odd 0.4 if j is even . ( 25
)
The initial density contour of odd-even perturbation is presented in the first plot of Fig. 10. The CFL=0.5 is used. We examine the solutions from different Riemann solvers. The contour plots after 1 × 10 4 time steps and 4 × 10 4 time steps are presented in Fig. 10. The density distributions along x-direction in the middle of domain are also plotted in the Fig. 11. From these results, oddeven decoupling is clearly visible with HLLC, with which original steady shock suffers from instabilities and dies out after 4×10 4 time steps. Whereas, HLL and proposed one are free from carbuncle instability. The test shows that the proposed formulation is robust than HLLC when solving strong shocks.

G. Sedov blast wave problem

A spherically symmetric explosion, known as Sedov blast wave problem, is considered here. The computational domain is [-1.1, 1.1] × [-1.1, 1.1] with 160 grid elements in each direction. The initial conditions are given by (ρ 0 , u 0 , v 0 , p 0 ) = (1.0, 0.0, 0.0, 10 -6 ) except in the central cell where an initial delta-function energy source was set in terms of pressure as

p c = (γ -1)ρ 0 ǫ 0 V c ( 26
)
in which V c denotes the volume of the cell containing the origin energy source and ǫ 0 is the total amount of released energy. In our numerical tests, the released energy is ǫ 0 = 0.244816. The simulations are carried out up to time t = 1.0 at which the shock front reaches a radius of R = 1. Fig. 12 shows the numerical results obtained by the HLL, HLLC and the proposed method.

As one can see, the results from the HLLC exhibit carbuncle instability across the shock front. Whereas, both HLL and the proposed method are capable of solving stable shock front. To make a quantitative comparison, we plot numerical solutions with different solvers and exact solution [START_REF] Sedov | Similarity and dimensional methods in mechanics[END_REF][START_REF] Peng | A direct ale multi-moment finite volume scheme for the compressible euler equations[END_REF][START_REF] Kumar | A skewed kernel approach for the simulation of shocks using sph[END_REF] in Fig. 13. It is obvious that HLLC produces the largest numerical errors due to carbuncle issue.

H. Double Mach reflection

To further examine that the proposed HLL-type Riemann solver is free from carbuncle phenomenon, a M a = 10 hypersonic propagating planar shock reflected by 30 • ramp is simulated in this case. It is well known that the kinked Mach stem may occur due to carbuncle instability from Riemann solvers. This problem is proposed in [START_REF] Woodward | The numerical simulation of twodimensional fluid flow with strong shocks[END_REF] . The computational domain is [0, 3.2] × [0, 1]. A right-moving Mach 10 shock is imposed with 60 • angle relative to x-axis. At the right boundary of the computational domain, the boundary condition is given by setting all gradients to be zero. The solution is computed up to t = 0.2 with different Riemann solvers. The results are displayed in Fig. 14. It is clearly visible that kinked Mach stem occurs in case of HLLC, which gives rise to spurious triple point and polluted region after shock. Whereas the spurious kinked Mach stem is completely eliminated with HLL or the proposed one. Moreover, the roll-up region around slip line is solved with less numerical diffusion by the new formulation compared with the original HLL.

I. Rayleigh-Taylor instability

Let us conclude the validation section with a simulation of the Rayleigh-Taylor instability problem [START_REF] Xu | Anti-diffusive flux corrections for high order finite difference weno schemes[END_REF][START_REF] Kumar | A multimass correction for multicomponent fluid flow simulation using smoothed particle hydrodynamics[END_REF] . Initially, a single contact discontinuity separating two fluids with different densities is set at y = 0.0. A small perturbation is given to the velocity. Due to the density difference, the heavy fluid will accelerate into the light one. The Rayleigh-Taylor instability will develop as the spike of the heavy fluid falls into the light one and the bubble of the light fluid rises into the ambient heavy fluid. The computation domain is [0, 0.25] × [0, 1]. The initial condition is given by (ρ 0 , u 0 , v 0 , p 0 ) = (2, 0, -τ cos(8πx), 1 + 2y) 0 ≤ y ≤ 0.5 (1, 0, -τ cos(8πx), 1.5 + y) otherwise .

(

) 27
where τ = 0.025 γp/ρ, and the specific heat is chosen to be γ = 5/3 in this case. To model the effect of density difference, a source term ρ is added to the right hand of the momentum equations in the y-direction and a source term ρv is added to the energy equation. The computation is carried out up to the t = 1.95s with 60 × 240 cell elements. As studied in [START_REF] Deng | Limiter-free discontinuity-capturing scheme for compressible gas dynamics with reactive fronts[END_REF][START_REF] Xu | Anti-diffusive flux corrections for high order finite difference weno schemes[END_REF] , low-dissipation methods are essential to resolve the sharp contact and associated small scale structures due to Rayleigh-Taylor instability. We plot the density contour calculated by different Riemann solvers in the left panel of Fig. 15.

Clearly HLL is the most diffusive solver and can not resolve sharp contact. Both HLLC and proposed HLL-BVD methods produce high resolution results. The proposed HLL-BVD method obtains sharper material interfaces compared with HLLC. To evaluate the numerical solution quantitatively, we plot the temporal evolution of the spike and bubble tip in the right panel of Fig. 15.

As studied in [START_REF] Kumar | A multimass correction for multicomponent fluid flow simulation using smoothed particle hydrodynamics[END_REF] , numerical solvers with excessive errors are not able to correctly predict the temporal evolution of instability. In our results, we find that HLL underestimates the evolution of bubble and spike characteristics due to excessive numerical dissipation errors. And the HLL and HLL-BVD produce very similar prediction.

IV. CONCLUDING REMARKS

This work proposes a novel formulation called HLL-BVD which is simultaneously capable of solving contact discontinuities accurately and of being free from shock instability. The missing contact in original HLL is restored with a jump-like THINC function. The resolution of contact wave is further improved by minimizing the dissipation term in HLL. Without modifying the waves structure explicitly, the robustness of original HLL has been preserved, as shown through one-dimensional and two-dimensional benchmark test cases. The new formulation significantly improves the resolution of contact wave and, most importantly, is free from carbuncle instability for strong shock. Future work will investigate whether the reduced dissipation across the contact interface also improves the behavior of the HLL solver for the description of viscous flow, often penalized by excessive dissipation across the shear wave [START_REF] Nishikawa | Very simple, carbuncle-free, boundary-layer-resolving, rotated-hybrid riemann solvers[END_REF][START_REF] Chen | Affordable shock-stable item for godunov-type schemes against carbuncle phenomenon[END_REF] .

- 5 FIG. 2 .

 52 FIG.2. Density profiles for an isolated right-moving contact wave travelling at Mach=0.1 (a) and Mach=10 -5 (b), as obtained with the HLL (solid), HLLC (dashed) and present (dot-dashed) solvers. The reference solution is represented by the thick line.

FIG. 3 .

 3 FIG.3. Density (a) and pressure (b) profiles for Sod problem with 100 cell elements, as obtained with the HLL (solid), HLLC (dashed) and present (dot-dashed) solvers. The reference solution is represented by the thick line.

FIG. 4 .

 4 FIG.4. Density (a) and pressure (b) profiles for Lax problem with 100 cell elements, as obtained with the HLL (solid), HLLC (dashed) and present (dot-dashed) solvers. The reference solution is represented by the thick line.

FIG. 5 .

 5 FIG.5. Global (a) and magnified (b) density profiles for strong Lax problem, as obtained with the HLL (solid), HLLC (dashed) and present (dot-dashed) solvers. The reference solution is represented by the thick line.

FIG. 6 .

 6 FIG.6. Global (a) and magnified (b) pressure profiles for strong Lax problem, as obtained with the HLL (solid), HLLC (dashed) and present (dot-dashed) solvers. The reference solution is represented by the thick line.

FIG. 7 .

 7 FIG.7. Density (a) and pressure (b) profiles for the slow moving shock problem, as obtained with the HLL (solid), HLLC (dashed) and present (dot-dashed) solvers. The reference solution is represented by the thick line.

FIG. 8 .

 8 FIG. 8. Density contour for 2D Riemann problems with Case 1 initial condition. Solutions using uniform grids of 600 × 600 from different formulation for HLL-type Riemann solvers are presented.

FIG. 9 .

 9 FIG. 9. Density contour for 2D Riemann problems with Case 2 initial condition. Solutions using uniform grids of 600 × 600 from different HLL-type Riemann solvers are presented.

FIG. 10 .FIG. 11 .

 1011 FIG. 10. Density contours for the odd-even perturbation problem as obtained with the HLL (left), HLLC (middle) and HLL-BVD (right) solvers. Initial condition (top), after 10000 (middle) and 40000 time steps.

FIG. 12 . 6 FIG. 13 .

 12613 FIG.12. Density contour calculated by different Riemann solvers until t = 1.0 for the Sedov blast wave problem. It is noted that in the x-and y-direction, some numerical oscillations caused by the carbuncle instability can be observed for HLLC.

FIG. 14 . 8 FIG. 15 .

 14815 FIG. 14. Density contours for double Mach reflection at t = 0.2 with different Riemann solvers.

  5, 0.5]. Two different configurations are studied in the present work. The initial conditions are set as follows,

	3.88									10.4
	3.87									10.35
	3.86									10.3
	3.85									10.25
	3.84									10.2
	3.83									10.15
	3.82									10.1
	3.81									10.05
	3.8									10
	0.59	0.6	0.61	0.62	0.63	0.64	0.65	0.66	0.67	0.59	0.6	0.61	0.62	0.63	0.64	0.65	0.66	0.67
				(a) Density							(b) Pressure

ACKNOWLEDGMENT

This work was supported in part by the fund from JSPS (Japan Society for the Promotion of Science) under Grant Nos. 17K18838 and 18H01366. We also acknowledge support from Labex MEC (ANR-10-LABX-0092) and the A*MIDEX project (ANR-11-IDEX-0001-02), funded by the "Investissements d'Avenir", for funding Xi Deng at M2P2.