David Coeurjolly
email: david.coeurjolly@liris.cnrs.fr

Isabelle Sivignon
email: isabelle.sivignon@gipsa-lab.grenoble-inp.fr

Efficient Distance Transformation for Path-based Metrics

In many applications, separable algorithms have demonstrated their efficiency to perform high performance volumetric processing of shape, such as distance transformation or medial axis extraction. In the literature, several authors have discussed about conditions on the metric to be considered in a separable approach. In this article, we present generic separable algorithms to efficiently compute Voronoi maps and distance transformations for a large class of metrics. Focusing on path-based norms (chamfer masks, neighborhood sequences), we propose efficient algorithms to compute such volumetric transformation in dimension n. We describe a new O(n • N n • log N • (n + log f )) algorithm for shapes in a N n domain for chamfer norms with a rational ball of f facets (compared to O( f n 2 • N n ) with previous approaches). Last we further investigate a more elaborate algorithm with the same worst-case complexity, but reaching a complexity of O(n • N n • log f • (n + log f )) experimentally, under assumption of regularity distribution of the mask vectors.

Introduction

Volumetric analysis of digital shapes is crucial in many geometry processing applications, for instance to measure distances between two points in Z n , or to measure the width of a shape or the proximity between two shapes. Since early works on digital geometry, distance transformation has been widely investigated (e.g. [START_REF] Rosenfeld | Distance functions on digital pictures[END_REF]). Given a finite input shape X ⊂ Z n , the distance transformation labels each point in X with the distance to its closest point in Z n \ X. Labeling each point by the closest background point leads to Voronoi maps (e.g. the restriction to Z n of Voronoi diagrams from computational geometry [START_REF] De Berg | Computational Geometry[END_REF]). Distance transformation (or distance field) is a key tool in many applications such as shape modeling, shape matching, geometry processing, motion planing, object tracking. . . (see [START_REF] Fabbri | 2D euclidean distance transform algorithms: A comparative survey[END_REF] or [START_REF] Jones | 3D distance fields: a survey of techniques and applications[END_REF] for surveys of tech-niques and applications). In this article, we propose a theoretical analysis to speed up the distance transformation computation for a large class of metric, allowing efficient and parallel implementations.

As the distance transform is parametrized by a distance function, many authors have addressed this distance transformation problem with trade-offs between algorithmic performances and the accuracy of the digital distance function with respect to the Euclidean one. Hence, authors have considered: distances based on chamfer masks [START_REF] Rosenfeld | Distance functions on digital pictures[END_REF][START_REF] Borgefors | Distance transformations in digital images[END_REF][START_REF] Fouard | 3-D chamfer distances and norms in anisotropic grids[END_REF] or sequences of chamfer masks [START_REF] Rosenfeld | Sequential operations in digital picture processing[END_REF][START_REF] Mukherjee | On approximating euclidean metrics by digital distances in 2D and 3D[END_REF][START_REF] Strand | Distance Functions and Image Processing on Point-Lattices With Focus on the 3D Face-and Body-centered Cubic Grids[END_REF]Normand et al., 2013a); vector displacement based Euclidean distance [START_REF] Danielsson | Euclidean distance mapping[END_REF][START_REF] Ragnemalm | The Euclidean Distance Transform[END_REF]; Voronoi diagram based Euclidean distance [START_REF] Breu | Linear time euclidean distance transform algorithms[END_REF][START_REF] Maurer | A linear time algorithm for computing exact euclidean distance transforms of binary images in arbitrary dimensions[END_REF] or square of the Euclidean distance [START_REF] Hirata | A unified linear-time algorithm for computing distance maps[END_REF][START_REF] Meijster | A general algorithm for computing distance transforms in linear time[END_REF]. For the Euclidean metric, separable volumetric computations have demonstrated to be very efficient with the design of optimal O(n • N n ) time algorithms for shapes in [1, N] n domains, optimal multithread/GPU implementation or extensions to toric domains (please refer to Coeur-jolly (2012) for a discussion).

Path-based approaches (e.g. chamfer mask orweighted-neighborhood sequences) approximate the Euclidean distance as the length of shortest paths defined from sequences of displacement vectors on the grid (from a finite set of possible moves). Aside distance information, path-based approaches provide an explicit notion of discrete path that is not accessible for the Euclidean norm. Furthermore, the discrete and combinatorial nature of the distance function has been used to define efficient algorithms to extract discrete medial axis [START_REF] Borgefors | Efficient shape representation by minimizing the set of centres of maximal discs/sphere[END_REF][START_REF] Remy | Medial axis for chamfer distances: computing look-up tables and neighbourhoods in 2D or 3D[END_REF][START_REF] Saha | A survey on skeletonization algorithms and their applications[END_REF] as local maxima of the distance map. Normand et al. (2013b[START_REF] Normand | A streaming distance transform algorithm for neighborhood-sequence distances[END_REF] further exploits the combinatorial structure of path-based distances to compute distance transformation in on-the-fly streaming context. In terms of distance transform computation, two main techniques exist. The first one considers a weighted graph formulation of the problem and Dijkstra-like algorithms on weighted graphs to compute distances. If m denotes the size of the chamfer mask, computational cost could be in O(m • N n ) using a cyclic bucket data structure as suggested by [START_REF] Verwer | An efficient uniform cost algorithm applied to distance transforms[END_REF]. Another approach consists in a raster scan of the domain: first the chamfer mask is decomposed into disjoint sub-masks; then the domain grid points are scanned in a given order (consistent with the sub-mask construction) and a local computation is performed before being propagated [START_REF] Rosenfeld | Sequential operations in digital picture processing[END_REF][START_REF] Borgefors | Distance transformations in digital images[END_REF]. Scanning the domain several times (one per submask) leads to the distance transformation values. Again, we end up with a O(m • N n ) computational cost. Besides specific applications which use the anisotropic nature of the chamfer mask, rotational dependency is usually enforced by increasing the mask size m (its number of vectors, see below) leading to expensive computational costs.

Contributions The goal of this work is to demonstrate that the linear factor in the mask size can be lowered down to a logarithmic one in any dimension for path-based metrics. This is achieved by first detailing and analyzing the separable distance transformation algorithm and briefly recalling the preliminary analysis of [START_REF] Coeurjolly | 2D subquadratic separable distance transformation for path-based norms, in: Discrete Geometry for Computer Imagery[END_REF] for the 2D case, before extending it to higher dimensional distance transformation problems. More precisely, we describe efficient and parallel algorithms in arbitrary dimension n to compute error-free distance transformation and Voronoi map for chamfer norms and other path-based metrics. Overall computational costs are summarized in Table 1 (see 3.1 for the predicate definitions).

The article is organized as follows: First, we recall basic definitions and properties of path-based norms (Section 2). In Section 3 we clarify the separable ndimensional Voronoi map extraction. Section 4 is a short discussion about the complexity of this algorithm for L p metrics. Then Section 5 is dedicated to the design of a fast implementation of the separable algorithm for path-based metrics. In Section 6, we present and analyse the proposed n-dimensional algorithm for path-based metrics.

Preliminaries

Metric space and distance transformation

A metric space (E, F, d) is a set E together with a metric d : E × E → F on the set. When E is equal to Z n and d is an integer-valued metric, also called digital metric, i.e. d : Z n × Z n → Z, we say that (Z n , Z, d) is a digital metric space. A digital shape is a finite subset of Z n .
Definition 1 (Voronoi Map and Distance Transformation) For a digital shape X ⊂ Z n , the Voronoi map V X associated with a digital metric space (Z n , Z, d) is the map

X → P(Z n \X) such that V X (a) = arg min b∈Z n \X {d(a, b)}. The distance transformation DT X is a map X → Z such that DT X (a) = d(a, b) for b ∈ V x (a).
The Voronoi map V X corresponds to the intersection between the continuous Voronoi diagram for the metric d of the pointset Z n \X and the lattice Z n . Note that V X (a) may contain several equidistant points to a in Z n \ X. In the following, we consider a restricted Voronoi map, denoted Π X , such that Π X (a) = b with b ∈ V X (a) choosen arbitrarily. Π X are not unique but provide the same distance map DT X . In the following, we focus on the restricted Voronoi map computation and we may omit the word restricted for the sake of clarity. Interested readers may refer to [START_REF] Couprie | Discrete bisector function and Euclidean skeleton in 2D and 3D[END_REF][START_REF] Hesselink | A linear-time algorithm for euclidean feature transform sets[END_REF] for separable algorithms to compute the complete Voronoi map of Definition 1.

Defining digital metrics spaces, notably in the context of digital image processing, has been the object of many works for the past fourty years. Note that (weighted, with

O(n) O(n) Θ(n • N n ) Hirata (1996) L ∞ O(n) O(n) Θ(n • N n ) Meijster et al. (2000) L 1 O(n) O(n) Θ(n • N n ) Meijster et al. (2000) L p (exact pred.) O(n • log p) O(n • log p • log N) O(n 2 • N n • log p • log N) Lemma 1 L p (inexact pred.) O(n) O(n • log N) O(n 2 • N n • log N) Lemma 1 2D Chamfer norm O(log m) O(log 2 m) O(log 2 m • N 2 ) Coeurjolly (2014) 2D Neig. seq. norm O(log m) O(log 2 m) O(log 2 m • N 2 ) Normand et al. (2013a) with Coeurjolly (2014) nD Chamfer norm O(n + log f ) O((n + log f ) • log N) O(n • N n • log N • (n + log f )) Lemma 4 w i ≥ 0) L p metrics d L p (a, b) = n ∑ k=1 w k |a k -b k | p 1 p , (1)
define metric spaces (Z n , R, d L p ) which are not digital.

However, rounding up the distance function, (Z n , Z, d L p ) is a digital metric space [START_REF] Klette | Digital Geometry: Geometric Methods for Digital Picture Analysis[END_REF]). However, it is not precise enough in many situations, and other approaches have been designed. Among them, the family of path-based metrics (chamfer norms, -weighted-neighbourhood sequences) aim at defining digital metrics induced by norms. In the following and for the sake of simplicity, we focus on chamfer norms but similar results can be obtained for more generic path-based metrics such as neighborhood sequences. Elements to support this claim are provided further in the following section.

Chamfer norms

Definition 2 (Chamfer Mask) A weighted vector is a pair ( v, w) with v ∈ Z n and w ∈ N\{0}. A chamfer mask M is a central-symmetric set of weighted vectors with no null vectors and containing at least a basis of Z n .

In most situations, vectors of a chamfer mask exhibit axial symmetries. As examples, see Figure 1(a) and (c), where a subset of vectors (together with their weights) defining chamfer masks by symmetries (called generators) are depicted.

From a chamfer mask, we can define a path between two points a and b as a sequence of k points {c i } such that

c 0 = a, c k-1 = b and c i c i+1 = v i ∈ M for i ∈ {0 . . . k -2}.
The length of this path is thus the sum of all weights associated with the vectors v i (i.e. ∑ w i ). As M contains a basis of Z n , such path between a and b always exists and we can define the chamfer distance between two points a and b in Z n as the length of the shortest path between a and b. Since weights are positive integers (see Def. 2), distance values are scaled by the weight of the first vector ((1, 0 . . . , 0) T by convention). Hence, using masks defined in Fig. 1,1 3

• d M 3-4 (a, b) and 1 5 • d M 5-7-11 (a, b) are approximations of d L 2 (a, b).
For all positive weights, a chamfer mask defines a metric. In many shape processing applications, we usually consider a subset of chamfer masks, the chamfer norms, with weights such that the induced metrics have convex unit balls, and thus leading to homogeneous distance functions. Chamfer norms can be characterized by a set of linear constraints on the mask weights [START_REF] Borgefors | Distance transformations in digital images[END_REF]; [START_REF] Strand | Distance Functions and Image Processing on Point-Lattices With Focus on the 3D Face-and Body-centered Cubic Grids[END_REF]. In the following, we define the size of a chamfer norm simply by the number of vectors of its associated chamfer mask.

Many authors have proposed algorithmic and/or analytic approaches to construct chamfer norms approximating the Euclidean metric. Following [START_REF] Thiel | Géométrie des distances de chanfrein[END_REF] and [START_REF] Fouard | 3-D chamfer distances and norms in anisotropic grids[END_REF], we briefly recall here the classical construction of chamfer norms from Farey set since it will be the base of the study proposed in Section 6.

The Farey sequence F n m of dimension n and order m is defined as follows :

F n m = {( x 2 x 1 , . . . , x n x 1 ), gcd i∈1..n (x i ) = 1, 0 ≤ x n ≤ x n-1 ≤ • • • ≤ x 1 ≤ m}. Then a Farey sequence F n m is in bijection with all the points (x 1 , . . . , x n ) in Z n , 0 ≤ x n ≤ • • • ≤ x 1 ≤
m visible from the origin1 . The vectors v k of a chamfer norm in dimension n can be defined using a subset of a particularF n m : the weights w k are set so that the rational ball B R (see Definition 3 below) is convex. By construction, such chamfer masks have axis symmetric unit balls and thus define chamfer norms.

Distance computation for chamfer norms

To evaluate distances between two digital points for a given chamfer norm, direct formulations have been proposed with a simple geometrical interpretation [START_REF] Thiel | Géométrie des distances de chanfrein[END_REF][START_REF] Normand | Medial axis lookup table and test neighborhood computation for 3d chamfer norms[END_REF], using the so-called rational ball. Definition 3 (Rational ball, minimal H-representation) Given a Chamfer norm M, the rational ball associated with M is the polytope

B R = conv v k w k ; ( v k , w k ) ∈ M . ( 2
)
where conv denotes the convex hull of a set of points.

Rational balls for some 2D and 3D chamfer norms are illustrated in Figure 1. As any convex polytope, the rational ball B R can also be described as the intersection of f linear constraints in dimension n, f being the number of (n -1)-facets of B R . This is the H-representation of the polytope which can be written in a matrix form:

B R = x ∈ R n ; x = m ∑ k=1 α k v k w k , α k ≥ 0 , m ∑ 1 α k = 1 = {x ∈ R n ; Ax ≤ y} ,
where A is a f × n matrix and y a vector of n values, so-called the H-coefficients (Ziegler, 2012). The Hrepresentation of a polytope P is with minimal param- [START_REF] Normand | Medial axis lookup table and test neighborhood computation for 3d chamfer norms[END_REF]. 2 In other words,

eter if P = {x ∈ Z n ; Ax ≤ y } with A being such that ∀k ∈ [1 . . . f ], ∃x ∈ P A k x = y k
A is the minimal H-representation of B R if each linear hyperplane of the H-representation of B R contains at least one point in B R ∩ Z n .
From [START_REF] Normand | Medial axis lookup table and test neighborhood computation for 3d chamfer norms[END_REF], an important result for distance computation can be summarized as follows:

Proposition 1 (Direct Distance Computation) Given a chamfer norm M and (A, y) its minimal parameter Hrepresentation, then for any a ∈ Z n , the chamfer distance of the point a from the origin is

d M (O, a) = max 1≤k≤ f {A k a T } .
(3)

Coming back to general path-based digital metrics, (weighted) neighborhood sequences have been proposed to have a better approximation of the Euclidean metric [START_REF] Rosenfeld | Sequential operations in digital picture processing[END_REF][START_REF] Mukherjee | On approximating euclidean metrics by digital distances in 2D and 3D[END_REF][START_REF] Strand | Distance Functions and Image Processing on Point-Lattices With Focus on the 3D Face-and Body-centered Cubic Grids[END_REF]Normand et al., 2013a). The main idea is to combine sequences of elementary chamfer norms. A key result has been demonstrated by Normand et al. (2013a) stating that for such distance functions, a minimal parameter polytope representation exists and that distances can be obtained from an expression similar to (3):

d(O, a) = max 1≤k≤ f {γ k (A k a T )} , (4)
γ k : N → N being some integer sequence characterizing the neighborhood sequence metric. As we will see in the next sections, direct distance computation is key to design an efficient distance transformation algorithm. Similarity of Equations 3 and 4 makes the algorithms presented in the following sections for chamfer norms easily generalizable to neighborhood sequences.

To conclude this preliminary section, algorithms efficiency is characterized by their asymptotic behavior using the O(•) and Θ(•) notations 3 as a function of the dimension, the number of vectors defining the chamfer norm and the domain size.

2 A k being the k th row of A.

3 In a computational model where arithmetic operations and scalar comparisons are constant time with:

f (x) = O(g(x)) ⇔ ∃C, x 0 ∈ R + , ∀x > x 0 , | f (x)| ≤ C • g(x), and f (x) = Θ(g(x)) ⇔ ∃C,C , x 0 ∈ R + , ∀x > x 0 ,C • g(x) ≤ | f (x)| ≤ C • g(x).

Separable distance transformation

3.1. Voronoi map from separable approach and metric conditions

Several authors have described optimal in time and separable techniques to compute error-free Voronoi maps or distance transformations for L 2 and L p metrics [START_REF] Breu | Linear time euclidean distance transform algorithms[END_REF][START_REF] Hirata | A unified linear-time algorithm for computing distance maps[END_REF][START_REF] Meijster | A general algorithm for computing distance transforms in linear time[END_REF][START_REF] Maurer | A linear time algorithm for computing exact euclidean distance transforms of binary images in arbitrary dimensions[END_REF]. Separability means that computations are performed dimension by dimension. In the following, we consider the Voronoi Map approach as defined by [START_REF] Breu | Linear time euclidean distance transform algorithms[END_REF]. Given a digital set X defined on an hyperrectangular domain

[1..N 1 ] × . . . × [1..N n ], let us first de- fine the image I X : [1..N 1 ]×. . .×[1..N n ] → {0, 1} such that I X (a) = 1 for a ∈ [1..N 1 ]×. . .×[1..N n ] iff a ∈ X (I X (a) = 0 otherwise).
The separable algorithm that computes the Voronoi Map for I X is defined in Algorithm 1 and works on the image spans for each dimension. An image span S along dimension q is a vector of N q points with same coordinates except at their q th one. The q th coordinate of a point a ∈ Z n is denoted by a q . A given span S in dimension q is denoted by {s i } i=1...N q . In Algorithm 1 the Voronoi map is first initialized by processing each span of the input image along the first dimension in order to create independent 1D Voronoi maps for the metric (lines 5 -6). Then, for each further dimension q, the partial Voronoi map Π X is updated using one dimensional independent processes on each span along the q th dimension (line 8). Algorithm 2 describes the function VORONOIMAPSPAN. This function is the core of the separable algorithm as it defines the 1D processes to perform on each row, column and higher dimensional image span. In this process, metric information are embedded in the following key predicates (see Fig. 2. HIDDENBY(a, b, c, S): given a 1D image span S parallel to the q th coordinate axis, and three points a, b, c ∈ Z n such that a q < b q < c q , this predicates returns true if there is no s ∈ S such that d(b, s) < d(a, s) and d(b, s) < d(c, s) .

(5)

Algorithm 1: VORONOIMAP(BINARY MAP I X )

1 Π X = empty image, same size as I X ; 2 for q in {1 . . . n} do

3 for (x 1 , ..x q-1 , x q+1 , ..x n ) in [1..N 1 ] × ..[1..N q-1 ] × [1..N q+1 ].. × [1..N n ] do 4 S = {s i } i∈[1.
.Nq] where s i = (x 1 ..x q-1 , i, x q+1 ..x n ); // all the coordinates are fixed in S except the q th one

5 if q == 1 then // Π X is initialized span by span 6 Π X = Π X ∪ VORONOIMAPSPAN(I X , q, S); 7 else // Π X is updated along span S 8 Π X = VORONOIMAPSPAN(Π X , q, S);
2: VORONOIMAPSPAN(MAP M X , DIMEN- SION q , 1D SPAN S)
Data: q is an integer in {1 . . . n}; S is a 1D span along dimension q, with points {s 1 , . . . , s Nq } sorted by their q th coordinate; M X is either a binary map if q = 1 or a partial Voronoi Map.

Result: Partial Voronoi map Π X updated along S. 1 if q == 1 ; // Special case for the first dimension 2 then 3 Π X = empty image, same size as M X ; 24 return Π X Remark. Note that by construction, for a given span S along dimension q, points a, b, c given as parameters to the HIDDENBY predicate necessarily verify a q b q c q . Indeed, these points are defined as the (partial) Voronoi map images of three points s i , s j , s k (i j k) of S, therefore having their q th coordinate equal to i, j, and k respectively (see lines 11, 12, 15 of Algorithm 2, and Figure 3 for an illustration).

4 k = 0; 5 foreach point s in S do 6 if M X (s) == 0 then // if s ∈ Z n \X 7 L S [k] = s; // L S = list of the sites visible on S 8 k + +; 9 else 10 Π X = M X ; 11 L S [0] = M X (s 1 ); 12 L S [1] = M X (s 2 ); 13 k = 2 , l = 3; // Update the list L S 14 while l ≤ N q do 15 w = M X (s l ); 16 while k ≥ 2 and HIDDENBY(L S [k -1], L S [k], w, S) do // L S [k] is no longer visible, unstack
For L 1 , L 2 and L ∞ metrics, CLOSEST and HIDDENBY predicates can be computed in O(n) in dimension n [START_REF] Breu | Linear time euclidean distance transform algorithms[END_REF][START_REF] Maurer | A linear time algorithm for computing exact euclidean distance transforms of binary images in arbitrary dimensions[END_REF]. Hence, Algorithm 2 is in O(n • N q ) for the dimension q, leading to an overall computational time for the Voronoi Map (Algorihtm 1) and Distance Transformation computations in Θ(n 2 • N n ) (if we assume that ∀q ∈ [1 . . . n], N q = N). Note that for L p metrics, we can derive a Θ(n • N n ) algorithm as suggested in [START_REF] Hirata | A unified linear-time algorithm for computing distance maps[END_REF]; [START_REF] Meijster | A general algorithm for computing distance transforms in linear time[END_REF] using the following observation: when evaluating the CLOSEST predicates in line 21 of Algorithm 2, we compare distances along the 1-D span of dimension q. If we store the partial power p of the distance to the closest site a for each grid point y for previous dimensions (i.e the sum ∑ q-1 i=1 (a iy i ) p ), such distance comparisons can be obtained in O(1). Similar argments can be used for the HIDDENBY predicates of line 16, leading to the overall computational cost in Θ(n • N n ).

s j M X (s j ) s i M X (s i ) S
Figure 3: For two points s i and s j (in purple) on a span S along dimension q (in red), the partial Voronoi map images M X (s i ) and M X (s j ) (in black) respectively have i and j as q th coordinate. [START_REF] Hirata | A unified linear-time algorithm for computing distance maps[END_REF] or [START_REF] Maurer | A linear time algorithm for computing exact euclidean distance transforms of binary images in arbitrary dimensions[END_REF] discussed about conditions on the metric d to ensure that Algorithm 2 is correct. The key property can be informally described as follows: given two points a, b ∈ Z n such that a q < b q and a straight line l along the q th direction and if we denote by v l (a) (resp. v l (b)) the intersection between the Voronoi cell of a (resp. b) and l, then v l (a) and v l (b) are simply connected Euclidean segments and v l (a) appears before v l (b) on l (so called monotonicity property by [START_REF] Maurer | A linear time algorithm for computing exact euclidean distance transforms of binary images in arbitrary dimensions[END_REF] and is related to quadrangle inequality by [START_REF] Hirata | A unified linear-time algorithm for computing distance maps[END_REF]. These contributions are summed up in Definition 4 and Proposition 2.

Definition 4 (Axis symmetric ball norm) A metric d induced by a norm whose unit ball is symmetric with respect to grid axes is called axis symmetric ball norm.

Proposition 2 (Metric conditions [START_REF] Hirata | A unified linear-time algorithm for computing distance maps[END_REF]) Algorithm 1 exactly computes the Voronoi Map Π X of a binary input image I X for any axis symmetric ball norm.

Proposition 2 implies that most chamfer norms and neighborhood sequence based norms can also be considered in separable Algorithm 1 (see Fig. 4). However, note that Algorithm 2, and as a by-product Algorithm 1, are exact only if the distance comparison predicate is exact, i.e. if we can compare distances, through the CLOSEST and HIDDENBY predicates, without error.

Furthermore, computational efficiency of the algorithm requires the design of efficient algorithmic tools to implement these predicates, and this the purpose of the next sections.

Generic predicates and complexity analysis for axis symmetric ball norms

We first detail the overall computational cost of Algorithms 2 and 1. We assume in the following that ∀q ∈ [1 . . . n], N q = N.

Lemma 1 [START_REF] Maurer | A linear time algorithm for computing exact euclidean distance transforms of binary images in arbitrary dimensions[END_REF][START_REF] Coeurjolly | 2D subquadratic separable distance transformation for path-based norms, in: Discrete Geometry for Computer Imagery[END_REF]) Let (Z n , F, d) be a metric space induced by a norm with axis symmetric unit ball. If C denotes the computational cost of CLOSEST predicate and H is the computational cost of the HIDDENBY predicate, then Algorithm 2 is in O(N

• (C + H)), leading to a complexity of O(n • N n • (C + H)) for Algorithm 1.
For a given axis symmetric ball norm d, generic Algorithms 3, 4 and 5 were defined in [START_REF] Coeurjolly | 2D subquadratic separable distance transformation for path-based norms, in: Discrete Geometry for Computer Imagery[END_REF]. Note that these algorithms are valid for any dimension n. The computational cost of the CLOSEST predicate is simply the one of a distance evaluation. As a first approach, Algorithms 4 and 5 show that the HIDDENBY predicate can be obtained by a binary search on the 1D image span S to localize the abscissa of Voronoi edges of sites {a, b} and {b, c}.

Algorithm 3: Generic CLOSEST(a, b, c ∈ Z n ). 1 return d(a, b) < d(a, c);
The complexity H of Algorithm 5 can be expressed as a function of the complexity C of Algorithm 3, leading to the general result below:

Lemma 2 (Coeurjolly 2014) Let M be a chamfer norm with axis symmetric unit ball in dimension n whose rational ball has f facets, Algorithm 1 can be implemented with a computational complexity of O(n

• N n • C • log N),
where N n is the size of the image.

Distance transformation for L p metrics

As a direct consequence of Lemma 1, we briefly derive computational costs for L p metrics. For such metrics, as discussed in Section 3.1, the CLOSEST and HIDDENBY Algorithm 4: Generic VORONOIEDGE(a, b, s i , s j ∈ Z n ) with i < j, a q < b q .

1 if ( ji = 1) then 2 if i = 1 and CLOSEST(s i , b, a) then 3 return -∞; 4 if i = N q and CLOSEST(s i , a, b) then 5 return ∞; 6 return i; Algorithm 5: Generic HIDDENBY(a, b, c ∈ Z n ; S in the q th direction) with a q < b q < c q .

7 mid = i + ( j -i)/2; 8 if CLOSEST(s mid , a, b) then //
1 v ab = VORONOIEDGE (a, b, s 1 , s Nq ); 2 v bc = VORONOIEDGE (b, c, s 1 , s Nq ); 3 return (v ab > v bc );
predicates are in O(n) for p = {1, 2, ∞} with exact integer only computations [START_REF] Maurer | A linear time algorithm for computing exact euclidean distance transforms of binary images in arbitrary dimensions[END_REF][START_REF] Meijster | A general algorithm for computing distance transforms in linear time[END_REF]. We thus have distance transformation algorithms in Θ(n 2 • N n ). Let us now show that Algorithm 3 and 5 lead to a faster algorithm for any p ≥ 1.

For p ∈ R, p ≥ 1, we can use approximations of the evaluation of distances on IEEE 754 double and then consider the Generic HIDDENBY predicate in O(n • log N) (Alg. 5). As predicates being based on floating point computations, numerical issues may occur but we have an O(n 2 • N n • log N) distance transformation algorithm (L p inexact predicates in Table 1). If p ∈ Z, p ≥ 3, we use exact integer number based computations of distances storing sum of power p quantities (which can be computed in O(n • log p) thanks to exponentiation by squaring). The HIDDENBY predicate is also based on Algorithm 5, leading to an O(n 2 • N n • log p • log N) distance transformation algorithm (L p exact predicates in Table 1).

Distance transformation in higher dimension for chamfer norms

In this section, we consider digital metrics given by chamfer norms and propose an efficient algorithm to compute the separable distance transformation for such metrics in nD. In [START_REF] Coeurjolly | 2D subquadratic separable distance transformation for path-based norms, in: Discrete Geometry for Computer Imagery[END_REF], the structure of the rational ball of a chamfer mask in dimension 2 was used to obtain an O(log m) algorithm for the CLOSEST predicate, and an O(log 2 m) one for the HIDDENBY predicate, leading to an overall O(log 2 m • N 2 ) algorithm for the separable distance transformation (see Table 1). The following sections extend these results to higher dimensions.

Definitions and general principle

Let us consider a general chamfer norm in arbitrary dimension n with m weighted chamfer vectors. As explained in Section 2, these vectors define a rational ball B R (see Definition 3), the center of which can be any point p. We define a wedge (p, f k ) as the conical hull of p (called the apex) and the vertices of a given facet

f k of B R (k ∈ {1.. f } if B R
has f facets). Thus, to each wedge is associated a row A k of the matrix of the minimal Hrepresentation of B R . Note that A k can also be seen as anon-unitary -normal vector to the facet f k , as quoted by [START_REF] Normand | Medial axis lookup table and test neighborhood computation for 3d chamfer norms[END_REF]. In the following, given a point a and a point p, we denote by (a, F a (p)) the wedge of apex a containing point p, F a (p) being one facet of B R (see Figure 5).

Using similar notations, [START_REF] Thiel | Géométrie des distances de chanfrein[END_REF] and [START_REF] Strand | Distance Functions and Image Processing on Point-Lattices With Focus on the 3D Face-and Body-centered Cubic Grids[END_REF] demonstrated that distance evaluation between a point a and a point p can be obtained in two steps: first, compute the wedge (a,

F a (p) = f k ); then d M (a, p) = A k .(p -a) T . ( 6
)
Thus, implementing the CLOSEST predicate comes down to computing the wedge a given point belongs to. In 2D, it was shown in Coeurjolly ( 2014) that a binary search over the chamfer vectors was enough. The nD case is discussed in section 5.2.

Let us now see how to optimize the HIDDENBY predicate, which comes down to optimizing the VORONOIEDGE function. Given two points a and b (a q < b q ) and a 1D image span S along the q th dimension, we have to find the point e of S (e ∈ Z n ) with abscissa e q such that all the points of S of abscissa lower than e q are in the Voronoi cell of a while all the points with a greater abscissa are in the Voronoi cell of b. To compute e q let us define l(S) as the one-dimensional flat that contains all the points of S and consider the Euclidean point ξ ∈ l(S) satisfying:

d M (a, ξ ) = d M (b, ξ ) . ( 7
)
In other words, we are looking for the point in l(S) which is equidistant to a and b. To compute ξ , we first suppose that we know the two wedges (a, F a (ξ ) = f k ) and (b, F b (ξ ) = f l ) (see Fig. 6-(b)). In this situation, ξ is the solution of

A k • (ξ -a) T = A l • (ξ -b) T . ( 8
)
Note that since ξ ∈ l(S), we have one linear equation with only one unknown, ξ q . As a consequence, if we know the two wedges point ξ belongs to, we have its q th coordinate ξ q in O(1). As A is the minimal representation of B R , it has rational components and thus ξ has integer coordinates execpt ξ q which is rational. Finally, the integer point e ∈ S is given such that e q = ξ q . The next section is dedicated to the CLOSEST predicate while in section 5.3 we detail how to efficiently compute the wedges (a, F a (ξ ) = f k ) and (b, F b (ξ ) = f l ), summarized in Algorithm 7 (see also Fig. 6-(b)).

CLOSEST predicate and first results

To begin with, let us discuss about the combinatorics of the chamfer norm rational ball. If m denotes the number of weighted vectors of M, its rational ball B R has O(m n 2 ) i-facets (for all 0 ≤ i ≤ n) (de [START_REF] De Berg | Computational Geometry[END_REF]. If f denotes the number of (n -1)-facets of B R , then we have:

Lemma 3 Let M be a chamfer norm whose rational ball B R has f (n -1)-facets in dimension n, then distance computation and thus CLOSEST predicate are in (amor-

tized) O(n + log f ) with O f n 2 (log f ) n 2 -δ
space and preprocessing time 4 . 4 δ is an arbitrarily small positive constant.

The proof is given in Appendix A.

Algorithm 4 being valid in any dimension, we can merge this result with Lemma 2 to straightforwardly obtain the result below:

Lemma 4 Let M be a chamfer norm whose rational ball

B R has f (n -1)-facets in dimension n, separable exact Voronoi Map Π X can be obtained in O(n • N n • log N • (n + log f )), thanks to a preprocessing in O f n 2 (log f ) n 2 -δ .
However, we show below that we can still expect faster VORONOIEDGE function even in higher dimension.

Improved HIDDENBY predicate

In dimension 2, it was shown in Coeurjolly ( 2014) that it was possible to reduce the complexity from a logarithmic factor on the size N of the image to a logarithmic factor on the size m of the mask using binary search over chamfer vectors. This process cannot be extended straighforwardly in higher dimensions since chamfer vectors cannot be ordered to perform a binary search anymore. However, it is interesting to notice that, whatever the dimension n, vectors from a given point a to any point of a span S lie in the smallest affine subspace containing a and the one-dimensional flat l(S).

l(S) s 1 s N a (a, F a (s N )) F a (s 1 )
Figure 5: Vectors s ia lie in a 2-flat defined by l(S) and a, in light green. The distance d M between a point on S and a is computed via a ray shooting that returns the n-1-facet of B R traversed by the ray s i -a : wedge (a, F a (s N )) and facet F a (s 1 ) are depicted in light red.

In the general case where a does not lie on S, this is actually always a 2-flat, denoted by P, and the intersection of this 2-flat with the rational ball B R is a polygon (see Fig. 5). The vertices of this polygon could be used to define a set of vectors on which a binary search could be performed as in the 2D case [START_REF] Coeurjolly | 2D subquadratic separable distance transformation for path-based norms, in: Discrete Geometry for Computer Imagery[END_REF]. However, since a 2-flat is actually the intersection of n -2 hyperplanes, computing this polygon comes down to intersecting a n-polytope with n -2 hyperplanes. By duality, each of these operations is equivalent to a convex hull computation, with a complexity of O( f n/2 ) [START_REF] Chazelle | An Optimal Convex Hull Algorithm in Any Fixed Dimension[END_REF][START_REF] Bajaj | Splitting a complex of convex polytopes in any dimension[END_REF]. As a consequence, in order to make the approach efficient, we must avoid to compute explicitly this intersection. This is achieved by rewriting the VORONOIEDGEWEDGE function as presented in Algorithm 6. As before, let e ∈ Z n be the point on S such that all the points of S of abscissa lower than e q are in the Voronoi of a while all the points with a greater abscissa are in the Voronoi cell of b. Given a, the goal of this function is to find the wedge of B R (of apex a) e belongs to.

The algorithm computes two points s i and s j such that s i belongs to the Voronoi cell of a (if a q < b q , b otherwise), s j to the Voronoi cell of b (if a q < b q , a otherwise) and either ji = 1 or F a (s i ) = F a (s j ). Similarly to generic Algorithm 4, this is done by performing a binary search over the points of S, with a key difference on exit conditions: now, the algorithm does not wait until point e is found, but exits as soon as the two points s i and s j belong to the same wedge (line 1). Indeed, by convexity of a wedge, this implies that any point on S between s i and s jand in particular e -also belongs to the same wedge. Correctness of the algorithm is ensured by maintaining two invariants : (i) s i is lower than s j on span S (s i q < s j q ) ; (ii) if a q < b q , s i is in a's Voronoi cell, s j in b's Voronoi cell, and conversely if b q < a q .

Figure 6 illustrates the first step of the binary search in (a), and the situation at the end of the search in (b) (projection in plane P).

It remains now to use Algorithm 6 to compute point e. Algorithm 7 implements the VORONOIEDGEND function as the nD counterpart of the VORONOIEDGE function of [START_REF] Coeurjolly | 2D subquadratic separable distance transformation for path-based norms, in: Discrete Geometry for Computer Imagery[END_REF]. First, lines 1 to 11 are dedicated to checking whether the bisector of a and b crosses span S or not. If it does not, algorithm exits with an error code (lines 5 or 11). Otherwise, wedges (a, F a (e)) (b, F b (e)) are computed lines 12-13 calling the VORONOIEDGEWED-GEND function (see Figure 6(b) for an illustration).

Algorithm 6: VORONOIEDGEWEDGEND(a, b ∈ Z n ; i, j ∈ Z, i < j; span S along the q th direction; faces Proposition 3 Let M be a chamfer norm in dimension n whose rational ball B R has f (n -1)-facets. Let W be the computational time complexity of the VORONOIEDGEWEDGEND function. Then, the separable exact Voronoi Map can be obtained in (amortized)

f k i = F a (s i ), f k j = F a (s j ) ) 1 if f k i = f k j or j -i = 1 then 2 return f k i ; 3 else 4 mid = i + ( j -i)/2; //
O(n • N n • (n + log f + W )) with a O f n 2 (log f ) n 2 -δ
space and preprocessing time. More precisely, the worst-case complexity W being O((n

+ log f ) • log N), this leads to a global (amortized) complexity of O(n • N n • log N • (n + log f )) (same preprocessing).
The proof is given in Appendix B.

Note that in the worst-case, this approach does not improve the result presented in Lemma 4 (using the generic VORONOIEDGE of Algorithm 4). However, in Section 6, we give some experimental insights on a finer analysis of the complexity W under distribution hypothesis.

S a

F a (s 1 ) s N s 1 F a (s N ) s mid (a) S a F a (s i ) = F a (s j ) = fk s j s i b s i s j F b (s i ) = F b (s j ) = fl e ? (b)

Experimental analysis

Insights on the complexity in dimension n

The complexity W of Algorithm 7 depends on the number of recursion steps done until points s i and s j are in the same wedge. Wedges being defined by the n -1 dimensional faces of B R , this complexity depends on the distribution of the chamfer vectors defining B R .

Let us denote by P the intersection between the (n -1)-faces of B R and the 2-flat P (see the red polygon on Fig. 5). Note that P goes through the center of B R . If we assume that the vectors defining B R are uniformely distributed on the unit sphere S n and that the faces of P are also uniformely distributed on B R ∩ P, then we can expect that W = O((n + log f ) • log |P|). Even if studying precisely these questions is out of scope of this work, in the following we give insights on both the relevance of these assumptions and the behaviour of |P| in the context of chamfer norms. // Check that the bisector of a and b crosses span

S 1 Compute f k 1 = F a (s 1 ), f l 1 = F b (s 1 ); 2 d(a, s 1 ) = A k 1 • (s 1 -a) T ; d(b, s 1 ) = A l 1 • (s 1 -b) T ;
3 if (a q < b q and d(b, s 1 ) < d(a, s 1 )) ; // or (b q < a q and d(a, s 1 ) < d(b, s 1 )) 4 then 5 Bisector does not cross S. return -1.

6 else 7 Compute f k N = F a (s N ), f l N = F b (s N ); 8 d(a, s N ) = A k N • (s N -a) T ; d(b, s N ) = A l N • (s N -b) T ; 9 if (a q < b q and d(a, s N ) < d(b, s N )) ; // or (b q < a q and d(b, s N ) < d(a, s N )) 10 then 11
Bisector does not cross S. return -1.

// Compute e 12 f k =VORONOIEDGEWEDGEND(a, b, 1, N, S, f k 1 , f k N ); 13 f l =VORONOIEDGEWEDGEND(b, a, 1, N, S, f l 1 , f l N ); 14 Compute abscissa ξ q of the point ξ ∈ S such that A k • (ξ -a) T = A l • (ξ -b) T ; 15 return ξ q ; 6.1.1

. Some observations on the distribution hypothesis

To study the distribution of chamfer vectors, we consider chamfer masks where vectors are defined from a subset of Farey sequences, as presented in Section 2 (see also [START_REF] Thiel | Géométrie des distances de chanfrein[END_REF] and [START_REF] Fouard | 3-D chamfer distances and norms in anisotropic grids[END_REF]). Studying the distribution of such sets of vectors is a field of research in itself, and we simply mention below several results relevant to our context.

First, it is well-known from [START_REF] Marklof | Fine-Scale Statistics for the Multidimensional Farey Sequence[END_REF]; [START_REF] Marklof | Visibility and directions in quasicrystals[END_REF] that n-dimensional lattice points visible from the origin have a constant density in R n . Moreover, [START_REF] Boca | Distribution of lattice points visible from the origin[END_REF] studied in the 2D case the distribution of the angles of straight lines from the origin through visible points. More precisely, they study the proportion of differences between consecutive angles which are larger than the average: they show that this proportion is smaller than what is expected for a random distribution, and give an explicit formulation of the repartition function. Similar results in higher dimension remain an open question.

These results tend to support the hypothesis of a uniform distribution of the vectors of B R , but the question of the distribution of the faces of the polygon P has not been investigated to our knowledge.

Experimental behaviour of |P|

In this part, we investigate the number of faces of P when B R is a rational ball defined from Farey Sequences. The results are presented in Figure 7 and we detail below how the rational balls are generated, how the 2-flats P are selected, and how the intersection between B R and P is performed.

In the four subfigures of Figure 7, rational balls are defined from Farey sequences:

• In (b-d), the vectors of B R are all normalized vectors of a Farey sequence of order m (the higher the order, the greater the number of vertices -and (n -1)-faces -of B R ). The order of the Farey sequences ranges from 1 to 10 in (b-c), from 1 to 6 in (d);

• in (a), B R is computed thanks to the algorithm presented by [START_REF] Fouard | 3-D chamfer distances and norms in anisotropic grids[END_REF]. 5 Given a (odd) mask size m, and a maximal error ε, the algorithm computes a subset of vectors of F m-1 2 and weights such that the rational ball B R is convex and the error with respect to the optimal theoretical error expected (wrt the Euclidean distance) for this mask size is below ε.

Once the sets of vectors defined, we use Qhull [START_REF] Barber | The quickhull algorithm for convex hulls[END_REF] to compute both the rational ball itself and its intersection with a 2-flat P that goes through the center of B R . This intersection is performed by randomly picking the coefficients of n -2 (n -1)-hyperplanes containing the center of B R , and iteratively adding each (n -1)hyperplane. The vertices of P are the points lying on all (n -1)-hyperplanes. 6For each rational ball, a certain number of cuts is performed: from 1000 in dimension 3 to only 6 in dimension 5 for rational balls obtained from Farey sequences of order 5 and 6 (due to precision issues in Qhull). 95% confidence intervals are depicted for each point (i.e. for each rational ball) as error bars, but most of the time too small to be visible on the graphs. Note that this remark suggests that the size of |P| does not depend on the position of P, thus supporting the uniform distribution hypothesis discussed in the previous section.

Analysing these results, we see that |P| seems to behave as f α , with α < 0 and decreasing when the dimension increases. This suggests that, in practice, the complexity W of Algorithm 7 is expected to be O((n + log f ) • log f ). Similarly to dimension 2 [START_REF] Coeurjolly | 2D subquadratic separable distance transformation for path-based norms, in: Discrete Geometry for Computer Imagery[END_REF], this approach is expected to lower down the worst case complexity of the computation of the distance transformation for chamfer norms in dimension n from a logarithmic factor on the size N of the image, to a logarithmic factor on the size f of the rational ball.

Distance transformation in dimension 2

We evaluate the performance of the separable approach to compute restricted Voronoi diagrams and distance transformation for chamfer norms in dimension 2. First, we observe that using Algorithm 1 with the nD VORONOIEDGE (Alg. 7), we obtain an overall complexity in O(log m • log N • N 2 ) which is close to the O(log 2 m • N 2 ) complexity of the ad-hoc 2D version of the problem [START_REF] Coeurjolly | 2D subquadratic separable distance transformation for path-based norms, in: Discrete Geometry for Computer Imagery[END_REF]. These complexities have to be compared with the O(m • N 2 ) complexity of the classical raster scan approach for chamfer norms [START_REF] Borgefors | Distance transformations in digital images[END_REF]. In Fig. 8-(a), we first illustrate some restricted Voronoi map results on small domains.

In Fig. 8-(bc), we have considered a 2D domain 2048 2 with 2048 random sites. First, we observe that fixing N, the log 2 m term is clearly visible in the computational cost of the Voronoi map (single thread curve). Bumps in the single thread curve may be due to memory cache issues. Please note that if we consider classical chamfer norm DT from raster scan (and sub-masks), the computational cost is in O(m • N 2 ) and thus has a linear behavior (green curves in Fig. 8-(bc)). Since we have a separable algorithm, we can trivially implement it in a multi-thread environment (here using OpenMP). Hence, on a bi-processor and quad-core (hyper-threading) Intel(R) Xeon(R) cpu (16 threads can run in parallel), we observe a speed-up by a factor 10 (orange curves in Fig. 8-(bc)). In Fig. 9, we present a shape processing experiment: when considering a chamfer norm mask with m = 100, both the raster scan and our separable approaches produce the same distance transformation (Fig. 9-(c)) but the raster scan approach has been obtained in 520ms whereas our separable approach only requires 25ms (note that since the domain is small -128 × 127-, the OpenMP scheduler has only considered a single core for the separable Voronoi map computation).

Implementation of all separable algorithms are publicly available in the DGtal library ( dgt ).

Conclusion and Discussion

In this article, we have proposed generic algorithms to efficiently solve the restricted Voronoi map and distance transformation problems for a large class of metrics in any dimension. Focusing on chamfer norms, geometrical interpretation of this generic approach allows us to design an algorithm with logarithmic factors in the chamfer mask size compared to a linear one for previous approaches. Thanks to separability, parallel implementation of the distance transformation leads to efficient distance computation for path based metrics.

For the L 2 metric, (additively) weighted Voronoi maps, also known as power maps, can be used to solve the reverse distance transformation and medial axis extraction problem using similar separable techniques [START_REF] Coeurjolly | Optimal separable algorithms to compute the reverse euclidean distance transformation and discrete medial axis in arbitrary dimension[END_REF]. A challenging future work would be to extend these results for path-based norms such as chamfer norms.

Figure 1 :

 1 Figure1: Chamfer masks and rational balls: in dimension 2, generator vectors for the mask M 3-4 (a), its rational ball (b). Generator vectors for M 5-7-11 (c) and its rational ball (d). In dimension 3, rational ball of a chamfer mask obtained using generator vectors (x, y, z) ∈ [[-3, 3]] 3 and weights computed following[START_REF] Fouard | 3-D chamfer distances and norms in anisotropic grids[END_REF].

  2): 1. CLOSEST(a, b, c): given three points a, b, c ∈ Z n this predicate returns true if d(a, b) < d(a, c);

  9 return Π X In other words, HIDDENBY returns true if and only if the Voronoi cells of sites a and c hide the Voronoi cell of b along S.

Algorithm

  ; l + +; 19 L S [k] = w; 20 foreach point s in S by increasing q th coordinate do 21 while (k < |L S |) and CLOSEST(s, L S [k + 1], L S [k]) do // s is closer to L S [k + 1], look further 22 k + + ; 23 Π X [s] = L S [k];

Figure 2 :

 2 Figure 2: Geometrical predicates for Voronoi map construction Coeurjolly (2014): HIDDENBY(a, b, c, S) returns true in (a) and false in (b) (straight segments correspond to Voronoi diagram edges). (c) illustrates the CLOSEST(a, b, c) predicate for a ∈ S.

Figure 4 :

 4 Figure 4: Distance transformation from a single source for different metrics satisfying Definition 4 and thus Proposition 2: (from left to right) L 1 , L 2 , L 4 , L 80 , M 3-4 and M 5-7-11 .

56

  Check whether s mid is closest to a or to b // O(n + log f ) evaluation of distances w.r.t. a and bComputef k = F a (s mid ); d(a, s mid ) = A k • (s mida) T ; Compute f l = F b (s mid ); d(b, s mid ) = A l • (s midb) T ; 7 if d(a, s mid ) < d(b, s mid ) then 8 if a q < b q then 9 return VORONOIEDGEWED-GEND(a, b, mid, j, S, f k , f k j ) 10 else 11 return VORONOIEDGEWED-GEND(a, b, i, mid, S, f k i , f k ) b, mid, j, S, f k , f k j )

Figure 6 :

 6 Figure 6: View of the 2-flat P: (a) Binary search initialization to compute the wedge (a, F a (e)). (b) After completion of two binary searches, both wedges (a, f k = F a (e)) and (a, f l = F b (e)) are known.

Algorithm 7 :

 7 VORONOIEDGEND(a, b ∈ Z n , span S).

Table 1 :

 1 Computational cost summary for separable Voronoi map computation on N n domains (m being the size of the chamfer norm and f the number of row in a H-representation of the mask, see below).

	Metric	CLOSEST	HIDDENBY	Sep. Voronoi Map	Reference
	L 2

  s mid closer to a

	9	return VORONOIEDGE(a, b, s mid , s j )
	10 else
		// s mid closer to b

11

return VORONOIEDGE(a, b, s i , s mid )

A point p ∈ Z n is visible from the origin in Z n if there is no point of Z n on (Op) between O and p.

Code is available on the TC18 website www.tc18.org/code_ data_set/code.php

Python code used to generate Farey sequences and to compute these graphs is available on http://www.gipsa-lab.fr/ ~isabelle. sivignon/recherches_en.html.

Acknowledgments

This project was supported in part by French ANR under Grant No. ANR-15-CE40-0006 (CoMeDiC).

Appendix A. Proof of Lemma 3

Similarly to the 2D case, the distance d M (O, a) for a ∈ Z n is given by first solving a ray-shooting problem on convex polytopes which consists in first computing the (n -1)-facet of B R pierced by the ray (O, a) (see Fig. 5). Once the facet is obtained, the associated A k row is used to evaluate d M (O, a) = A k • a T in O(n). Following [START_REF] Matousek | On ray shooting in convex polytopes[END_REF] Theorem 10, such a ray-shooting query on convex polytopes can be solved in

. In the case when the ray hits a facet of dimension strictly lower than n -1, the algorithm returns one of the adjacents (n -1)-facets. Propositions 3 and 4 from [START_REF] Normand | Medial axis lookup table and test neighborhood computation for 3d chamfer norms[END_REF] ensure that the choice of any (n -1)facet leads to the same distance evaluation. Please note also that the preprocessing time is roughly equivalent to the convex hull computation in higher dimension which is in O( f n 2 ). Hence, preprocessing for ray-shooting can be done while computing the rational ball B R using Eq. (2).

Appendix B. Proof of Proposition 3

Following Lemma 1, the generic separable algorithms computes the Voronoi map in O(n

space and preprocessing time. Remains to evaluate H, i.e. the complexity of the VORONOIEDGEND function. In Algorithm 7, the first eleven lines are in O(C) since only distance computations are involved. Lines 12 and 13 are calls to the VORONOIEDGEWEDGEND function, with a complexity in O(W ). In the worst case, we have W = O((n + log f ) • log N) thanks to the test ji = 1 on line 1 of Algorithm 6. Last, the system to solve in line 14 has only one unkwown ξ q since ξ belongs to the onedimensional span S, with a complexity of O(1). Figure 7: Number of faces of P with respect to the number of (n -1)-faces of B R in different settings. In (a), B R is a rational ball as computed by [START_REF] Fouard | 3-D chamfer distances and norms in anisotropic grids[END_REF] in dimension 3. In (b-c), B R is defined from a Farey sequence of given order and dimension, taking all the fractions : (b) dimension 3, for orders between 1 and 10, (c) dimension 4 for order between 1 and 10, (d) dimension 5 for orders between 1 and 6. Each point is the mean of a certain number of random cuts (1000 in dimension 3, 500 in dimension 4, 400 in dimension 5 for orders up to 4, and 6 in dimension 5 for orders 5 and 6.).

M 5-7-11 M 5-7-11 , and the L p metrics L 2 and L 80 . The digital shape X considered is the whole [0, 256] 2 domain except two random points in the first two columns, and 10 random points in the last two ones.(b) Experimental evaluation of the subquadratic algorithm when increasing the mask size on a [0, 2048] 2 image and following the chamfer norm construction of [START_REF] Fouard | 3-D chamfer distances and norms in anisotropic grids[END_REF] (zoom in (c)). We compare the efficiency of Algorithm 1 in single thread and multi-thread settings, with the classical raster scan approach of [START_REF] Borgefors | Distance transformations in digital images[END_REF]


D Andrault
email: denis.andrault@uca.fr

G Morard

G Garbarino

M Mezouar

M A Bouhifd

T Kawamoto

Melting behavior of SiO2 up to 120 GPa

Keywords: SiO2-silica, melting diagram, melt structure, lower mantle pressures

The structure of liquid silicates is commonly described as a statistical mixture of various atomic entities with relative abundances that can vary with pressure, temperature and composition.

Unfortunately, this view remains largely theoretical due to scarce experimental reports on the silicate melt structure, in particular under pressure. We performed X-ray diffraction of the SiO2 end-member to probe the melting curve up to ~120 GPa and 7000 K, and the melt structure up to ~80 GPa. We confirm the steep increase of the melting curve above ~14 GPa when stishovite becomes stable over coesite in subsolidus conditions, with a slope of about 80 K/GPa. Then, around 45 GPa and 5400 K, the melting curve flattens significantly, an effect most likely reflecting the densification of the SiO2 melt structure. The signal of diffuse X-ray scattering is compatible with a change of the Si coordination number from 4 to 6 along the melting curve, in agreement with previous works reporting a similar evolution during the cold compression of SiO2-glass. Because of the limited pressure range (within 10 to 20 GPa) in which the melting curve changes its slope, we speculate a difficult coexistence of tetrahedral SiO4 and octahedral SiO6 units in SiO2 melt at high pressures.

I. INTRODUCTION

Although SiO2 is rarely present as a pure phase in mantle rocks, it is a dominant component of terrestrial planets. Because SiO2 is highly refractory, its melting curve significantly affects the melting diagram of silicates present in the crust and the mantle. At subsolidus conditions, SiO2 undergoes a major transition from coesite to stishovite, at a pressure of ~14 GPa for a temperature of 3000 K [START_REF] Zhang | Melting and subsolidus relations of SiO2 at 9-14 GPa[END_REF]. The transition involves a change of the Si coordination number from tetrahedral "SiO4" to octahedral "SiO6" units, accompanied with a density change of ~30% [START_REF] Akaogi | Thermodynamic properties of stishovite by low-temperature heat capacity measurements and the coesite-stishovite transition boundary[END_REF].

This transition occurs in all silicates with a transition pressure that depends on composition. The structural changes are major and induce a divergence of physical and chemical properties between the shallower and deeper parts of planetary mantles.

The change of the Si coordination number is also expected to occur in silicate melts [START_REF] Sanloup | Structural change in molten basalt at deep mantle conditions[END_REF]. For SiO2, the experimental investigation on the melt is strongly hampered by a steep increase of the melting curve at pressures higher than the coesite to stishovite transition: The melting temperature was reported to increase from 3000 to 3900 K between ~14 GPa and ~22 GPa [START_REF] Shen | Measurement of melting temperatures of some minerals under lower mantle conditions[END_REF].

Such a steep Clapeyron slope (dT/dP=V/S) denotes a large volume of fusion between stishovite and a low density melt (LDM) containing Si in 4-fold coordination. In contrast, the Clapeyron slope of the melting curve is almost flat below ~14 GPa, which implies a negligible density contrast between coesite and LDM [START_REF] Zhang | Melting and subsolidus relations of SiO2 at 9-14 GPa[END_REF]. The brutal change of slope of the melting curve at ~14 GPa shows that LDM is maintained at pressures above the Si coordination change in the solid state.

The pressure-induced Si coordination change in amorphous SiO2 was investigated experimentally at room temperature [START_REF] Benmore | Atomistic modeling of multiple amorphousamorphous transitions in SiO2 and GeO2 glasses at megabar pressures[END_REF][START_REF] Lin | Electronic bonding transition in compressed SiO2 glass Phys Rev B 75[END_REF][START_REF] Murakami | Spectroscopic Evidence for Ultrahigh-Pressure Polymorphism in SiO2 Glass[END_REF][START_REF] Sato | High-pressure structural transformation of SiO2 glass up to 100 GPa[END_REF] and theoretically [START_REF] Meade | High-pressure X-ray diffraction of SiO2 glass[END_REF][START_REF] Stixrude | Structure and freezing of MgSiO3 liquid in the Earth's lower mantle[END_REF][START_REF] Takada | Molecular dynamics study of liquid silica under high pressure[END_REF][START_REF] Usui | Ab Initio Two-Phase Molecular Dynamics on the Melting Curve of SiO2[END_REF]. In both cases, the transformation to an octahedral-based SiO2 structure is clearly demonstrated. At high temperatures, the formation of a high-density melt (HDM), with a similar structure to the high-density glass potentially, would induce a smaller density contrast between stishovite and HDM, compared to stishovite and LDM. This would induce the flattening of the Clapeyron slope of the melting curve with increasing pressure. This line of reasoning agrees with the report of a moderate increase of the melting temperature from 5000 K to 5800 K with increasing pressure from 40 to 120 GPa [START_REF] Usui | Ab Initio Two-Phase Molecular Dynamics on the Melting Curve of SiO2[END_REF]. Still, large inconsistencies remain. For example, the melting temperature was reported from 4400 K [START_REF] Lin | Electronic bonding transition in compressed SiO2 glass Phys Rev B 75[END_REF] to 5250 K (Usui and Tsuchiya 2010) at a pressures around 50 GPa. In addition, neither the change in Clapeyron slope of the SiO2 melting curve at increasing pressure, nor the change in Si coordination number in the melt have been properly documented.

II. MATERIALS AND METHODS

Experimental design

We first refined the melting curve of SiO2 up to ~120 GPa using the laser-heated diamond anvil cell (LH-DAC) coupled with in situ X-ray diffraction. Samples consisted of small shards of pure SiO2 glass, typically of 20-50 µm across and 10-15 µm thick. The glass was loaded in the pressure chamber between two KCl layers of similar thickness. The use of a KCl makes the X-ray analysis more difficult, due to intense KCl peaks; however, a critical advantage is to minimize the axial thermal gradients within the sample. In addition, soft KCl layers favor homogenous sample pressure at high temperature. The amount of water trapped in the pressure chamber was minimized by performing all sample-loading operations under N2-flux in a glove bag. Samples were heated alternatively using a CO2 laser on one side or two fiber-lasers on both sides of the DAC. We found that SiO2 absorbs relatively well the 1 µm radiation of fiber-lasers after the sample is already heated to high temperatures using the CO2. In all cases, laser spots were more than 20 µm in diameter.

In a second series of experiments, we probed the structure of SiO2 melt at high pressures. We did not use a pressure medium in this case to prevent artifacts that could arise from its diffuse scattering, if the pressure medium eventually melts. Yet, a small grain of KCl was loaded in the gasket hole for the measurement of the nominal pressure. The sample geometry consisted in two thin disks of SiO2 quartz single crystal in direct contact with the diamonds, between which a little bit of Si-metal powder was distributed. The role of Si-metal is to couple with fiber-lasers in order to initiate the sample heating before the annealed SiO2 could eventually absorb the laser radiation by itself. We preferred the addition of Si to Pt powder that would have also provided an internal pressure calibrant, because Pt melts at lower temperature than SiO2 at high pressures; its signal of diffuse scattering would have hampered the study of SiO2-melt. The presence of a minor excess of Si in SiO2 melt should not affect significantly the melt structure, nor the mechanism of the pressure-induced Si coordination change.

Upon laser heating, SiO2 eventually melts at mid-distance between the two diamond culets. Then, further increase of the laser power induces a thicker layer of molten SiO2. In such a sample experiencing large axial temperature gradient, some SiO2 always remain solid at the contact with cold diamonds. More importantly, the temperature gradient within the molten part of the sample should be minor, due to turbulent convection in the high-pressure melt.

Determination of pressure and temperature in the laser heated diamond anvil cell

We used the equation of state of KCl [START_REF] Dewaele | High-pressure high-temperature equation of state of KCl and KBr[END_REF] to determine the nominal pressure at 300K before and after laser heating. Then, the sample pressure at high temperature was determined using two different cross-calibrated methods: (i) the P-V-T equation of state of stishovite [START_REF] Wang | P-V-T equation of state of stishovite up to mid-lower mantle conditions[END_REF]) and (ii) the nominal pressure corrected by XThKT to account for the increase of pressure in the laser spot (XTh is an adjustable parameter,  and K are thermal expansion and bulk modulus, respectively, of coesite [START_REF] Akaogi | Thermodynamic properties of stishovite by low-temperature heat capacity measurements and the coesite-stishovite transition boundary[END_REF]) at nominal pressures below 14 GPa and stishovite [START_REF] Wang | P-V-T equation of state of stishovite up to mid-lower mantle conditions[END_REF] for pressures above). We used samples loaded in KCl and heated below the melting curve to calibrate the value of XTh. For this sample geometry, the two methods give pressures within +/-5 GPa for XTh values increasing linearly from 35% to 90%, when temperature is increased from 4000 to 7000 K, respectively. At the melting temperature, the second method is preferred, because the remaining fraction of stishovite grains may encounter a temperature lower than the SiO2-melt, because they are located closer to the diamonds. The second method was also used to determine the pressure for DAC loadings without KCl pressure medium.

Temperatures up to 7000 K were determined based on spectroradiometric measurements (Fig. 1).

To make these measurements accurate, a critical parameter of the optical set-up is the use of reflective Schwarzfield-type telescopes and the absence of any refractive device along the optical path, in order to prevent chromatic aberration [START_REF] Schultz | Double-sided laser heating system for in situ high pressure-high temperature monochromatic x-ray diffraction at the ESRF High Pressure[END_REF]. Temperature measurements are performed with a precision of ±100 K and a reproducibility better than ±30 K. When not using KCl insulating layers in the DAC, the temperature uncertainty could be twice larger, due to steeper axial temperature gradient.

We note that the measured temperature corresponds to the hottest part of the sample, where thermal emission is the most intense.

The sample behavior was monitored in situ using an X-ray beam much smaller than the hot-spot on the sample (see below). The fluorescence (in the visible range of wavelength) of the irradiated KCl pressure medium, or irradiated Re-gasket, was used to adjust the spectrometer position, in order to perform the temperature measurement at the position of the X-ray beam. The spectrometer alignments were carefully checked before and after each melting experiments.

X-ray diffraction measurements

At the ID27 beamline of the ESRF (Grenoble, France), we used a monochromatic X-ray beam of 0.3738 Å wavelength focused by two Kirkpatrick-Baez mirrors onto a spot of less than 2x2 μm 2 . Diffraction patterns were recorded using a MAR-160 CCD detector and integrated using DIOPTAS [START_REF] Prescher | DIOPTAS: a program for reduction of two-dimensional X-ray diffraction data and data exploration[END_REF]. Acquisition times extended to 10 sec and 300 sec for solid and molten SiO2 samples, respectively. For experiments dedicated to the determination of the melting curve, volumes and relative abundances of SiO2 and KCl were refined using the XRDUA code [START_REF] De Nolf | XRDUA: crystalline phase distribution maps by two-dimensional scanning and tomographic (micro) X-ray powder diffraction[END_REF]. For experiments dedicated to probe the structure of SiO2-melt, we inserted Soller slits between the DAC and the CCD detector. These slits filter out a major fraction of the Compton diffusion coming from diamonds [START_REF] Weck | Use of a multichannel collimator for structural investigation of low-Z dense liquids in a diamond anvil cell: Validation on fluid H2 up to 5 GPa[END_REF], which results in a significant improvement of the signal of X-ray diffuse scattering coming from the melt. We first recorded the X-ray diffuse scattering of SiO2 glass under compression up to ~50 GPa, before we performed measurements at high temperatures in SiO2 melt up to 80 GPa. In both cases, we also recorded the diffuse background coming from diamond anvils (Compton diffusion). For the melt, small diffraction peaks of not-molten stishovite were eventually visible in the Q region between 1 and 8 Å -1 . The peaks were carefully fitted by Gaussian functions and subtracted. Then, after subtraction of the diamond Compton diffusion, the scattering signal coming from the sample was normalized to obtain the scattering factor S(Q). Finally, the wavelength position of the first sharp diffraction peak (FSDP) was extracted from the S(Q) using a Gaussian function coupled with a polynomial background.

III. RESULTS AND DISCUSSION

A steep SiO2 melting curve up to 45 GPa

We first refine the melting curve of SiO2 up to ~120 GPa using the LH-DAC coupled with in situ X-ray diffraction. With increasing temperature, the glass starting material first crystallizes in coesite below 14 GPa, stishovite or in the CaCl2-form above ~60 GPa [START_REF] Andrault | Pressure-induced landau-type transition in stishovite[END_REF]. The presence of thin and continuous diffraction rings are evidences for fine crystalline SiO2 powders under a moderate temperature gradient within the X-ray beam (Fig. 2A). At a sample temperature interpreted hereafter as the melting point of SiO2, several changes occur simultaneously (such as reported in [START_REF] Andrault | Melting of subducted basalt at the core-mantle boundary[END_REF])): (i) Additional increase of the laser power does not induce an increase of the sample temperature (discussed in (Geballe and Jeanloz 2012)); (ii) texture of diffraction rings becomes spottier; (iii) position of diffraction spots on the image plate changes rapidly at constant temperature;

(iv) integrated intensity of the SiO2 diffraction peaks decreases severely relative to the KCl pressure medium (Fig. 2); (v) upon quenching, intense stishovite peaks reappear, which certifies the presence of non-crystalline (i.e. molten) SiO2 within the X-ray spot at high temperature. The reason why all diffraction peaks of SiO2 do not disappear abruptly at the melting temperature is the unavoidable axial temperature gradient between the two diamonds. Also, the SiO2 melt may absorb the laser radiation less than solid SiO2, resulting in a relatively lower heating efficiency when the sample is partially molten. This can explain the temperature plateau despite an increase of the laser power.

Our measurements confirm the steep increase of the SiO2 melting temperature above ~14 GPa (Fig. 3, Table 1), in agreement with previous experimental [START_REF] Shen | Measurement of melting temperatures of some minerals under lower mantle conditions[END_REF] and theoretical [START_REF] Belonoshko | Molecular-dynamics of stishovite melting Geochim[END_REF][START_REF] Usui | Ab Initio Two-Phase Molecular Dynamics on the Melting Curve of SiO2[END_REF] studies. Within experimental uncertainties, the melting curve appears linear between ~15 and ~45 GPa. The Clapeyron slope of ~80 K/GPa remains ~4 times smaller than that reported at ~300 K/GPa for the coesite to stishovite transition [START_REF] Akaogi | Thermodynamic properties of stishovite by low-temperature heat capacity measurements and the coesite-stishovite transition boundary[END_REF]. Basic thermodynamic relations apply around the triple point where coesite, stishovite and the LDM of SiO2 coexist (Table 2). As a consequence, we estimate the volume of fusion of stishovite (2xVMelting/(VLDM+VSti)) at ~28% from the combination of the volume change of ~28% at the subsolidus transition [START_REF] Akaogi | Thermodynamic properties of stishovite by low-temperature heat capacity measurements and the coesite-stishovite transition boundary[END_REF]) and a negligible volume of fusion of coesite derived from the quasi-horizontal melting line between 9 to 14 GPa [START_REF] Zhang | Melting and subsolidus relations of SiO2 at 9-14 GPa[END_REF]). On the other hand, the entropy of melting of stishovite (Sm) is estimated to ~60 J/molK, a value 4 times larger than the entropy change at the coesite-stishovite transition [START_REF] Akaogi | Thermodynamic properties of stishovite by low-temperature heat capacity measurements and the coesite-stishovite transition boundary[END_REF].

A relatively flat SiO2 melting curve from 45 to 90 GPa

The response of our samples to laser heating was different below and above a pressure of ~45 GPa.

Between ~15 and ~45 GPa, after the onset of sample melting between ~3000 and ~5400 K (Fig. 3), further increase of the laser power always yield to the same pressure-temperature conditions of ~45 GPa and ~5400 K. This behavior could result from the very large volume of fusion (~28%, see above).

Increasing the amount of melt by LH within a nearly constant sample volume in the DAC would logically result in a significant increase of the sample pressure. Above ~45 GPa, this effect is not observed anymore; full sample melting could be achieved within the X-ray spot at each given experimental pressure.

Between ~45 to ~90 GPa, the SiO2 melting curve presents a Clapeyron slope of ~13 K/GPa, which is ~6 times lower than that observed at lower pressures (Fig. 3). The difference cannot result from subsolidus properties, because the stability field of stishovite extends continuously between ~15 and ~90 GPa below the melting line. Instead, it can be related to the evolution of the structure of the SiO2 melt. A smaller Clapeyron slope (dT/dP) above 45 GPa denotes either a smaller volume of fusion (Vm) and/or a higher entropy of fusion (Sm). A higher Sm above 45 GPa is improbable, because we expect large structural similarities between stishovite and HDM, as they both contain most of the Si in a 6-fold coordination [START_REF] Benmore | Atomistic modeling of multiple amorphousamorphous transitions in SiO2 and GeO2 glasses at megabar pressures[END_REF][START_REF] Lin | Electronic bonding transition in compressed SiO2 glass Phys Rev B 75[END_REF][START_REF] Meade | High-pressure X-ray diffraction of SiO2 glass[END_REF][START_REF] San | Polyamorphism of liquid silica under compression based on five order-parameters and two-state model: a completed and unified description High[END_REF][START_REF] Sato | High-pressure structural transformation of SiO2 glass up to 100 GPa[END_REF]. Instead, their structural similarities suggest a lower Sm. Therefore, the flattening of the melting slope above ~45 GPa is likely to result from a decrease of the volume of fusion due to the presence of a HDM of SiO2 above the melting curve. Now, if we consider a fictive triple point between LDM, HDM and stishovite at ~45 GPa and 5400 K, we can use Clapeyron relations (V=0, S=0 around triple points, and dT/dP=V/S) to retrieve the thermodynamical parameters of these three phases. Assuming a negligible entropy difference between LDM and HDM (SLDM-HDM), the stishovite volume of fusion (Vm) is estimated to ~5.3% above ~45 GPa (see details of the calculation below Table 2). Vm becomes 6.6 % or 4.0 % for hypothetical SLDM-HDM values of e.g. -15 or 15 J/mol.K, respectively. Such range of Vm values is comparable to the density contrast (reported at ~3 %) between stishovite and a glass compressed above 45 GPa at room temperature [START_REF] Petitgirard | SiO2 Glass Density to Lower-Mantle Pressures[END_REF]). Then, we can estimate a volume difference of ~25% between LDM and HDM melts. This value is marginally affected (+/-2%) by changing SLDM-HDM in a range of possible values. Such volume change is comparable to the 28% reported between coesite and stishovite [START_REF] Akaogi | Thermodynamic properties of stishovite by low-temperature heat capacity measurements and the coesite-stishovite transition boundary[END_REF]. We finally estimate the SiO2 LDM and HDM densities at ~3.26 10 3 and ~4.17 10 3 kg/m 3 , respectively, at 45 GPa and 5400 K (Table 2). Such LDM density is consistent with a previous theoretical report considering less than 10% Si in 6-fold coordination in the melt [START_REF] Takada | Molecular dynamics study of liquid silica under high pressure[END_REF].

Melting of the CaCl2-form of SiO2 above 90 GPa

Above ~90 GPa, the slope of the SiO2 melting curve increases slightly to ~35 K/GPa to yield a melting temperature of ~7000 K at ~120 GPa. The progressive increase of the Clapeyron slope could be related to the 2 nd order phase transformation from stishovite to the CaCl2-type polymorph (Fig. 3), even though the location of this phase boundary in the P-T diagram remains controversial (see [START_REF] Fischer | Equations of state and phase boundary for stishovite and CaCl2-type SiO2[END_REF]) and references therein). In fact, our X-ray measurements confirm that the subsolidus phase transformation occurs between ~80 and ~100 GPa at ~5400 K (Fig. 4). The change of slope of the melting curve is unlikely to come from a difference of volume between the two SiO2 polymorphs, because it was reported to be insignificant up to ~120 GPa, at least at 300 K [START_REF] Andrault | Equation of state of stishovite to lower mantle pressures[END_REF]). In contrast, the higher entropy of the CaCl2-form, compared to stishovite, implies a lower entropy of melting (Sm) above the transition pressure and, consequently, a steeper melting curve (because dT/dP=Vm/Sm). In addition, the strain energy becomes significant with increasing pressure above the transition. The CaCl2 distortion of the stishovite-based lattice observed at ~107 GPa and 5400 K corresponds to a strain energy of ~2 kJ/mole [START_REF] Andrault | Equation of state of stishovite to lower mantle pressures[END_REF]. This value is significant compared to the enthalpy of fusion of quartz of 9.4 kJ/mole [START_REF] Akaogi | Thermodynamic properties of stishovite by low-temperature heat capacity measurements and the coesite-stishovite transition boundary[END_REF]). Thus, the stabilization energy associated to the CaCl2-distortion can modify the Clapeyron slope of the melting curve.

Comparison with previous works

Previous LH-DAC experiments suggested an even steeper Clapeyron slope than ours above ~14 GPa [START_REF] Shen | Measurement of melting temperatures of some minerals under lower mantle conditions[END_REF]. However, after correction of the nominal pressures of the early experiments regarding the pressure increase due to laser heating in the DAC, the two data sets are in fact quite similar to each other. On the other hand, our melting curve plot in relatively good agreement with theoretical calculations, except in a few aspects (Fig. 3). The experimental melting slope remains steep, with little pressure-induced flattening, up to pressures higher than reported in theoretical studies [START_REF] Belonoshko | Molecular-dynamics of stishovite melting Geochim[END_REF][START_REF] Lin | Electronic bonding transition in compressed SiO2 glass Phys Rev B 75[END_REF][START_REF] Usui | Ab Initio Two-Phase Molecular Dynamics on the Melting Curve of SiO2[END_REF]. When it eventually flattens around 45 GPa and up to 90 GPa, the SiO2 experimental melting curve is found parallel, but at 200-300 K higher temperatures, compared to the most recent calculation [START_REF] Usui | Ab Initio Two-Phase Molecular Dynamics on the Melting Curve of SiO2[END_REF].

The significant increase of melting temperature observed experimentally above ~90 GPa is not reported in any theoretical studies yet.

We note an important discrepancy, of e.g. ~1000K for the melting temperature at 45 GPa, among the previous theoretical studies [START_REF] Belonoshko | Molecular-dynamics of stishovite melting Geochim[END_REF][START_REF] Lin | Electronic bonding transition in compressed SiO2 glass Phys Rev B 75[END_REF][START_REF] Usui | Ab Initio Two-Phase Molecular Dynamics on the Melting Curve of SiO2[END_REF]. It could come from the difficulty to model properly the local structure in SiO2 melt, in particular the strength of the Si-based structural units in LDM and HDM, respectively. Then, compared to shock compressions, our measurements plot at the upper limit of the Hugoniot plots of (solid) quartz and fused quartz (see [START_REF] Akins | Dynamic compression of SiO2: A new interpretation[END_REF] and references therein). Using this technique, one generally identifies the melting temperature when melting is truly achieved at higher impact energies, rather than in the solid at the highest temperatures, before the solid melts. However, this remains a matter of interpretation and the range of temperature uncertainties is large.

Analysis of the SiO2-melt structure factors

Polyamorphism in liquid SiO2 has been extensively discussed in the past by theoretical approaches (e.g. [START_REF] Benmore | Atomistic modeling of multiple amorphousamorphous transitions in SiO2 and GeO2 glasses at megabar pressures[END_REF][START_REF] Lin | Electronic bonding transition in compressed SiO2 glass Phys Rev B 75[END_REF][START_REF] San | Polyamorphism of liquid silica under compression based on five order-parameters and two-state model: a completed and unified description High[END_REF][START_REF] Takada | Molecular dynamics study of liquid silica under high pressure[END_REF])). To address experimentally the structural evolution of liquid SiO2 under high pressures, we performed LH-DAC experiments coupled with the Soller slits system installed at the ID27 beamline [START_REF] Weck | Use of a multichannel collimator for structural investigation of low-Z dense liquids in a diamond anvil cell: Validation on fluid H2 up to 5 GPa[END_REF].

This maximizes the signal of the SiO2 melt over the noise from the sample environment. No KCl pressure medium was used for these experiments. First, we collected diffraction patterns for the coldcompression of SiO2 glass up to ~45 GPa, in order to validate our experimental measurements of the structure factors S(Q). Our results are found in good agreement with previous studies (e.g. [START_REF] Prescher | Beyond sixfold coordinated Si in SiO2 glass at ultrahigh pressures[END_REF]). Spectra present a first sharp diffraction peak (FSDP) around 2 Å -1 , a second, less intense, around 5 Å -1 and a new contribution growing with increasing the pressure between 3 and 3.5 Å -1 (Fig. 5). These experimental features have been correlated to the gradual compaction of the structure of amorphous SiO2 [START_REF] Benmore | Atomistic modeling of multiple amorphousamorphous transitions in SiO2 and GeO2 glasses at megabar pressures[END_REF][START_REF] Prescher | Beyond sixfold coordinated Si in SiO2 glass at ultrahigh pressures[END_REF][START_REF] Sato | High-pressure structural transformation of SiO2 glass up to 100 GPa[END_REF].

We then recorded diffuse patterns of SiO2 melt between ~25 and ~75 GPa, which we normalized to obtain S(Q) (Fig. 6). In comparison to the glass, the same FSDP position is found at higher pressure in the melt (Figs. 7 and8). In addition, the new diffraction peak starts growing between 40 and 54 GPa, and 18.5 and 30.5 GPa, for the melt and the glass, respectively (Fig. 5). In the melt, the shift of the FSDP occurs prior to the appearance of the new diffraction peak, a similar behavior as in the glass [START_REF] Benmore | Atomistic modeling of multiple amorphousamorphous transitions in SiO2 and GeO2 glasses at megabar pressures[END_REF][START_REF] Prescher | Beyond sixfold coordinated Si in SiO2 glass at ultrahigh pressures[END_REF][START_REF] Sato | High-pressure structural transformation of SiO2 glass up to 100 GPa[END_REF]. This effect has been associated to the first regime of the compression where the structural changes are dominated by the compaction of the oxygen network [START_REF] Wu | Structure and Properties of Dense Silica Glass Scientific Reports 2[END_REF], before Si changes its coordination number at higher pressures.

At pressures higher than ~54 GPa, the additional contribution between 3 and 3.5 Å -1 appears clearly in the signal of SiO2-melt diffuse scattering. At this point, the position of the FSDP is found above 2.2 Å - 1 and, at 76 GPa, the melt presents the same FSDP position as the glass at ~45 GPa (Fig. 8). All experimental features are similar to those observed in the glass upon change of Si coordination number from 4 to 6 [START_REF] Benmore | Atomistic modeling of multiple amorphousamorphous transitions in SiO2 and GeO2 glasses at megabar pressures[END_REF][START_REF] Inamura | Transformations in the intermediate-range structure of SiO2 glass under high pressure and temperature[END_REF][START_REF] Meade | High-pressure X-ray diffraction of SiO2 glass[END_REF][START_REF] Prescher | Beyond sixfold coordinated Si in SiO2 glass at ultrahigh pressures[END_REF][START_REF] Sato | High-pressure structural transformation of SiO2 glass up to 100 GPa[END_REF]).

It appears clear that the structural changes occur at higher pressures in the liquid than in the glass.

Regarding the extremely high temperatures considered (between 5250 K and 6200 K, Fig. 6), this could be an effect of thermal pressure; the sequence of structural compaction may take place for similar SiO2 volumes in the glass and in the melt. However, higher temperature requires higher pressure to achieve a similar SiO2 volume. This hypothesis is supported by the large similarities between diffraction patterns recorded in the liquid and in the glass.

IV. CONCLUSION

The compaction mechanisms of liquid and amorphous SiO2 are complex and involve various atomic processes. Our experimental measurements suggest a change of the SiO2 melt structure in a relatively narrow range of pressures, within 10-20 GPa, around 45 GPa. Such evolution is compatible with a previous work investigating the Si coordination change from 4 to 6, e.g. in a SiO2-glass using transverse acoustic wave velocities [START_REF] Murakami | Spectroscopic Evidence for Ultrahigh-Pressure Polymorphism in SiO2 Glass[END_REF]. In our experiments, the change in the SiO2 melt structure is recognized based on (i) a change in the Clapeyron slope of the SiO2 melting curve and (ii) the occurrence of the new diffraction peak between 3 and 3.5 Å -1 . The relative sharpness of the structure change in the melt suggests a dominant role of the covalent character of Si-O bonds in the melt [START_REF] Lin | Electronic bonding transition in compressed SiO2 glass Phys Rev B 75[END_REF], which could induce a significant difference in energy between the tetrahedral "SiO4" and octahedal "SiO6" electronic configurations as a function of pressure and temperature. The two different melt structures could be relatively incompatible with each other, in agreement with previous reports [START_REF] Benmore | Atomistic modeling of multiple amorphousamorphous transitions in SiO2 and GeO2 glasses at megabar pressures[END_REF][START_REF] Lin | Electronic bonding transition in compressed SiO2 glass Phys Rev B 75[END_REF][START_REF] Takada | Molecular dynamics study of liquid silica under high pressure[END_REF]).

FIGURES AND TABLES

Figure 1: Typical temperature measurements. Thermal emission spectra were recorded at the center of the laser spot for sample temperatures of e.g. 5010 K (left) and 6140 K (right), using the online optical system installed at the ID27 beamline. The raw spectra (upper frames) are corrected by the response of the optical system (itself calibrated with a W-lamp at 2600 K) and fitted using the Wien law (middle frames). The so-called "two-color" methods (lower frames) shows a temperature variation of less than 10% of the mean temperature within a range of wavelength from ~670 to ~870 nm, which establishes the good quality of the thermal emission spectra. Temperature measurements are performed with a precision of ±100 K and a reproducibility better than ±30 K. We report temperatures of solid (blue) and liquid (red) SiO2 together with a smooth trend through the data set (grey band). Our temperature and pressure uncertainties are 200 K and 5 GPa, respectively. We also plot previous melting curve determinations based on early experimental works [START_REF] Shen | Measurement of melting temperatures of some minerals under lower mantle conditions[END_REF][START_REF] Zhang | Melting and subsolidus relations of SiO2 at 9-14 GPa[END_REF]) (light blue), theoretical calculations (purple [START_REF] Usui | Ab Initio Two-Phase Molecular Dynamics on the Melting Curve of SiO2[END_REF], blue [START_REF] Belonoshko | Molecular-dynamics of stishovite melting Geochim[END_REF] and green [START_REF] Lin | Electronic bonding transition in compressed SiO2 glass Phys Rev B 75[END_REF]) and shock wave experiments (orange [START_REF] Akins | Dynamic compression of SiO2: A new interpretation[END_REF]). For the latter, we report Hugoniot plots in solid and molten states for starting materials composed of quartz and fused-quartz. GPa: The stishovite model fits well the upper diffraction pattern recorded at ~75 GPa and ~4500 K. In contrast, the diffraction pattern recorded at ~102 GPa and ~5400 K is not well fitted with the stishovite model (middle frame), but instead with the CaCl2-form model (lower pattern). At these conditions, the CaCl2-form presents (a, b, c) unit cell parameters of (3.831(6), 3.967(5), 2.543(3)). Black, red and blue profiles are experimental data, background and Rietveld fit, respectively. Spectra recorded at high temperature in the melt at pressures of 38 GPa (orange), 40 GPa (bleu), 54 GPa (grey) and 76 GPa (yellow) are bracketed between spectra recorded at 300 K in the glass compressed to ~18.5 GPa (green) and ~49.5 GPa (dark blue). The contributions observed for the melt below and above ~45 GPa are comparable to those of the glass at ~18.5 GPa and ~49.5 GPa, respectively. It indicates an evolution of the SiO2 melt structure between 40 and 54 GPa comparable to that reported during the compression of the glass (Fig. 5). Table 1. Melting temperatures: Sample temperatures are directly provided by radiospectrometry (Fig. 1). Sample pressures are derived from cross correlation between two methods: the P-V-T equation of state of stishovite and an empirical model of thermal pressure (see text). Experimental uncertainties are within 5 GPa and 100 K.

Table 2. Clapeyron relations: Volumes and entropies relations calculated around (i) the true triple point where coesite, stishovite and LDM coexist and (ii) the fictive triple point where LDM, HDM and stishovite would coexist (see text). Thermodynamical parameters come from [START_REF] Akaogi | Thermodynamic properties of stishovite by low-temperature heat capacity measurements and the coesite-stishovite transition boundary[END_REF][START_REF] Wang | P-V-T equation of state of stishovite up to mid-lower mantle conditions[END_REF][START_REF] Zhang | Melting and subsolidus relations of SiO2 at 9-14 GPa[END_REF]. ) [START_REF] Meade | High-pressure X-ray diffraction of SiO2 glass[END_REF][START_REF] Prescher | Beyond sixfold coordinated Si in SiO2 glass at ultrahigh pressures[END_REF][START_REF] Inamura | Transformations in the intermediate-range structure of SiO2 glass under high pressure and temperature[END_REF][START_REF] Benmore | Atomistic modeling of multiple amorphousamorphous transitions in SiO2 and GeO2 glasses at megabar pressures[END_REF][START_REF] Sato | High-pressure structural transformation of SiO2 glass up to 100 GPa[END_REF] Our data; Quench Glass Our data; Glass compress. ) triple points e, derived from the flat melting curve of coesite between 9 and 14 GPa [START_REF] Zhang | Melting and subsolidus relations of SiO2 at 9-14 GPa[END_REF] f, calculated using DV=VHDM-VLDM g, calculated at the considered P-T conditions using the EoS of [START_REF] Wang | P-V-T equation of state of stishovite up to mid-lower mantle conditions[END_REF] h, calculated using Vcoes=Vsti+DV(C-S) i, calculated using VLDM=Vsti+DV(LDM-S) j, calculated using VHDM=Vsti+DV(HDM-S) k, DV of the transition normalized to the mean value between the two phases considered

Figure

  Figure 2A: X-ray evidences for melting of SiO2 at ~40 GPa. Time elapses from bottom to top. Temperature uncertainty is +/-100 K. After crystallization of stishovite, continuous rings remain clearly visible up to more than 5000 K. At 5100 K the intensity of SiO2 diffraction lines is already significantly decreased. At sample temperatures of 5600 and 5400 K, the only remaining diffraction lines come from the KCl pressure medium. Diffraction lines of stishovite reappear when the laser power is decreased to 5050 K, evidencing the SiO2 recrystallization.

Figure 2B :

 2B Figure 2B: X-ray evidence for the melting of SiO2 at 105 GPa. Time elapses from top to bottom. For temperatures lower than 6200 K, the set of diffraction lines evidences the presence of the CaCl2-form of SiO2 (see Fig.4) in addition to the KCl pressure medium. The intensity of the SiO2 peaks decreases severely at 6300-6400 K before they almost disappear at 6500 K. The stishovite to KCl ratio of diffraction peaks intensity decreases by a factor of 10 around between 6200 and 6500 K.

Figure 3 :

 3 Figure 3: Melting curve of SiO2 up to ~120 GPa. We report temperatures of solid (blue) and

Figure 4 :

 4 Figure 4: Subsolidus transition from Stishovite to the CaCl2-form between 75 and 102

Figure 5 :

 5 Figure 5: Structure factor of the compressed glass: As previously noticed (e.g. (Prescher et al. 2017) and references therein), a new contribution around 3-3.5 Å -1 grows above ~20 GPa.

Figure 6 :

 6 Figure 6: Structure factors S(Q) of SiO2 recorded in the melt between ~38 and ~76 GPa.

Figure 7 :

 7 Figure 7: Diffuse scattering signals of melt and quenched glass. The quenched glass(orange, recorded at 8 GPa and 300 K) presents a FSDP position significantly sharper and at slightly higher Q position compared to the corresponding melt (blue, recorded at 26 GPa and ~5150 K). This suggests less densification of the local structure of SiO2 melt at high temperature, compared to the glass. We did not successfully produced a glass at experimental pressures exceeding 45 GPa, because of full sample crystallization.

Figure 8 :

 8 Figure 8: Pressure dependence of the first sharp diffraction peak. Our measurements performed in melts (Red) and in glasses either quenched from high temperature (Blue) or during cold-compression (Green) are compared with previous reports. The shift of the FSDP to high wavenumbers occurs at lower pressures in the glass, compared to the melt.

  Fig.1

  LDM and HDM stand for Coesite, Stishovite, low-density and high-density SiO2 melts *, this study a, from Akaogi et al. (2011) b, calculated based on the Clapeyron relation dT/dP=DV/DS c, value assumed to be -15, 0 or 15 J/molK (see text) d, derived from Clapeyron relations around the (C, S, LDM) or (C, LDM, HDM

	Fig. 8
	Transition	Press. Temp. dT/dP	DS	DV	Vstish	Vcoes	VLDM	VHDM	DV/V
	considered	GPa	K	K/GPa	mol/K	cm3/mol	cm3/mol	cm3/mol	cm3/mol	cm3/mol
	C => S	15	3023	312 a -15.28 a -4.76736 b	14.6916 g	19.459 h			-0.279 k
	C => LDM	15	3023		74.9 d	0 e		19.459 h
	S => LDM	15	3023	80 *	59.59 b	4.76736 g	14.6916 g		19.459 i		0.2792 k
	S => HDM	45	5400	12.5 *	74.59 d	0.9324 b	13.6434 g			14.5758 j 0.0661 k
	LDM=>HDM	45	5400		-15 c -3.83496 f			18.4108 i	14.5758 j	-0.233 k
	S => HDM	45	5400	12.5 *	59.59 d	0.7449 b	13.6434 g			14.3883 j 0.0531 k
	LDM=>HDM	45	5400		0 c -4.02246 f			18.4108 i	14.3883 j	-0.245 k
	S => HDM	45	5400	12.5 *	44.59 d	0.5574 b	13.6434 g			14.2008 j	0.04 k
	LDM=>HDM	45	5400		15 c -4.20996 f			18.4108 i	14.2008 j	-0.258 k
	C, S,
								Our data; Melt

ACKNOWLEDGEMENTS

We thank anonymous reviewers for helpful comments. This research was financed by the French Government Laboratory of Excellence initiative n°ANR-10-LABX-0006, the Région Auvergne and the European Regional Development Fund. This is Laboratory of Excellence ClerVolc contribution N°XX.