Nicolas Ray
email: nicolas.ray@inria.fr

Dmitry Sokolov
email: dmitry.sokolov@univ-lorraine.fr

Bruno Lévy

Practical 3D Frame Field Generation

Keywords: smooth frame fields, remeshing Concepts:, Theory of computation → Computational geometry;, Computing methodologies → Mesh geometry models;

Figure 1: Our algorithm produces smooth frame fields in volumes. Frames (a) are represented by spherical harmonic functions (b), attached to each vertex of a tetrahedral mesh. Streamlines and singularities of the field are shown in yellow and red, respectively.

Introduction

Given a volume of interest, a 3D frame field can be defined at each point of the volume as a set of 6 unit vectors, globally invariant by rotations of π/2 around any of its vectors. One of the main motivations for generating such direction fields is to use them to steer the placement of objects within the volume of interest, as in hex-dominant meshing. Thus, hexahedral meshing is decomposed into two steps: (1) smooth frame field design and (2) hexahedral partitioning of the domain aligned with the frame field.

We make the remark that previous frame field generation methods address the 2D and 3D cases with very different strategies:

• In 2D, frame field design can be restated as vector field design (a simpler problem), thanks to the introduction of the "representation vector". It replaces a frame (a set of 4 vectors invariant by rotations of kπ/2) with a 2D vector, and defines the distance between two frames as the distance between the corresponding 2D vectors. It leads to a simple optimization problem that can easily support constraints and data fitting terms.

• In 3D, it is (unfortunately) not trivial to extend the idea of a "representation vector". Instead, [START_REF] Huang | Boundary aligned smooth 3d cross-frame field[END_REF] propose to represent frames by functions defined on the unit sphere (refer to Fig. 1). They derive a definition of the field smoothness from this representation. The existing approaches optimize it in two steps: (1) initialize the field by optimizing spherical harmonics coefficients (Huang et al.) or propagate the boundary constraints greedily [START_REF] Kowalski | A PDE based approach to multi-domain partitioning and quadrilateral meshing[END_REF]) (2) further smooth this result by minimizing a non-convex objective function, using L-BFGS [START_REF] Liu | On the limited memory bfgs method for large scale optimization[END_REF]]. In this smoothing step, the frames are represented by Euler angles.

We experimented with Huang's approach in 2D ( §2): each frame is represented by a sine function. The corresponding 2D optimization problem means finding for each frame its coefficients in the Fourier basis ( §2.2), that will minimize the L 2 norm between adjacent frames ( §2.3), subject to unit function norm constraint ( §2.4).

We found that the instantiation of [START_REF] Huang | Boundary aligned smooth 3d cross-frame field[END_REF] in 2D reproduces exactly the same optimization problem as in the current 2D state-of-the-art ( §2.5).

Therefore, using the notations introduced for the 2D case, we can describe the 3D version in §3 and extend the 2D optimization mechanism to 3D in §3.5. In the end, our 3D version is very similar to Huang's initialization, except that it exactly enforces boundary constraints.

• The first difficulty is to find the expression of the boundary conditions. On the boundary the frames are free to rotate around the surface normal ( §3.4) and it is difficult to enforce the alignment constraint while preserving this rotational degree of freedom. Previous work [START_REF] Huang | Boundary aligned smooth 3d cross-frame field[END_REF]] only partially enforces the alignment condition, resulting in a poor initialization of the optimization procedure as evaluated in §4.2.

• The second difficulty is that a frame is represented in a 9D Fourier basis, but the set of admissible functions (those that correspond to a frame) has dimension 3 (all possible rotations of the reference frame). Thus, in 2D one simply needs to normalize the representation vector, whereas in 3D we need to find the nearest point on the 3D manifold of admissible vectors embedded in 9D space.

In the present paper, we discuss some evaluation criteria related to our main envisioned application (hexahedral dominant meshing), and we show that our algorithm outperforms the state of the art methods [START_REF] Huang | Boundary aligned smooth 3d cross-frame field[END_REF][START_REF] Kowalski | A PDE based approach to multi-domain partitioning and quadrilateral meshing[END_REF]. Furthermore, we show that minor modifications of existing algorithms make them competitive with ours, and that our initialization alone can be sufficient in many cases.

Previous works

Frame field design on surfaces

Optimizing a frame field inside a 2D shape is very similar to optimizing a direction field on a 3D surface. The differences between a 2D shape and a 3D surface come from the curvature of the surface in the 3D case (angle defect, hard constraints, curvature fitting term).

Direction field design on surfaces was introduced by [START_REF] Hertzmann | Illustrating smooth surfaces[END_REF] for orienting strokes in non-photorealistic rendering. In their method, directions are represented by an angle rotation θ per vertex, and the smoothing is performed by a non-linear solver (BFGS). In the reference cited above, the "representation vector"

was not yet introduced, and the optimization mechanism they use is very similar to the smoothing step used in [START_REF] Huang | Boundary aligned smooth 3d cross-frame field[END_REF][START_REF] Kowalski | A PDE based approach to multi-domain partitioning and quadrilateral meshing[END_REF] for the 3D frame fields. Solving with a representation vector v = (r cos(θ), r sin(θ)) for each θ was suggested in [START_REF] Ray | Periodic global parameterization[END_REF] for faster results, and improved later for better control over the field topology [START_REF] Ray | Geometry-aware direction field processing[END_REF]]. Based on representation vectors, [Palacios and Zhang 2007a] proposed to control the field topology by local operations. For direction fields without constraints, fixing the norm of the entire solution instead of the norm of each representative vector [Knöppel et al. 2013] allows to find optimal direction fields by solving an eigenvector problem.

Directly optimizing for the angle θ makes it possible to perfectly control the field topology, but this is obtained at the expense of solving a mixed-integer system [START_REF] Ray | Nsymmetry direction field design[END_REF][START_REF] Bommes | Mixedinteger quadrangulation[END_REF]. In [Palacios and Zhang 2007a], representation vectors are introduced together with their duality with N th order traceless symmetric tensors. This relation is very useful to unify 2D and 3D frame fields.

Frame field design in volumes

In the computer graphics domain, [START_REF] Huang | Boundary aligned smooth 3d cross-frame field[END_REF] bring the functional representation of 3D frames: they represent frames by a function defined on a sphere and they store its decomposition in the spherical harmonics basis.

Their initialization step first creates a smooth spherical harmonic field. Then, for each sample, they compute the 3D frame that is best aligned with the spherical harmonic. This initial solution is finally improved by applying an optimized rotation to each frame. These rotations are defined by Euler angles and obtained by minimizing the field smoothness with L-BFGS. Boundary alignment is enforced by a penalty term. et al. [2012] propose an alternative initialization method. They design a 2D frame field on the volume boundary, convert it to a 3D frame field by adding the surface normal and its opposite vector, and then propagate it inside the volume. The resulting field is finally smoothed by optimizing a rotation for each frame, similarly to what is done in [START_REF] Huang | Boundary aligned smooth 3d cross-frame field[END_REF]], but with a better optimization scheme. They also optimize the singularity graph of the field by local combinatorial operations, as done by [START_REF] Jiang | Frame field singularity correction for automatic hexahedralization[END_REF].

Li

Our method optimizes for the same objective function as the references cited above, but with an initialization step that may be sufficient for remeshing applications. Our smoothing (same as [START_REF] Kowalski | A PDE based approach to multi-domain partitioning and quadrilateral meshing[END_REF]] but with a different sampling) is only required for surfaces of revolution. The performance of our method is compared with previous work in §4.

Functional representation in 2D

In this section, we explain how 2D algorithms like [START_REF] Ray | Periodic global parameterization[END_REF][START_REF] Kowalski | A PDE based approach to multi-domain partitioning and quadrilateral meshing[END_REF]] can be reformulated with a functional representation for each frame (instead of a "representation vector"). The goal of this section is to introduce the concepts and notations in the simple 2D setting in order to facilitate understanding the 3D setting (later in Section 3). Note that in the 3D setting there will be no Gauss curvature. For this reason, we discuss here only the flat 2D setting (field defined on a subset of IR 2 instead of a 3D surface).

Definitions for the 2D setting

Given a 2D triangulated shape, frame field design in 2D consists in finding a smooth frame field aligned with the boundary of the shape. We formulate it as minimizing the field curvature, based on the following definitions:

• A frame is a set of 4 unit vectors f = {f k }, k ∈ [0, 3] that is
invariant by a rotation of π/2 (Figure 2). It can be represented by the angle θ such that ∀k, f k = (cos(θ + kπ/2), sin(θ + kπ/2));

• a frame field is a frame per vertex of the triangulation;

• the boundary constraint: a normal on the boundary must match one of member vectors of the corresponding boundary frame;

• the rotation angle between two frames is the angle ∆θ of the rotation that transforms one frame into the other. This angle being defined modulo π/2, we choose the ∆θ with minimum absolute value;

• a triangle is said to be singular if the sum of the rotation angles over its three edges is not equal to 0.

Representing a frame field by angles is a natural and simple idea, but it makes it difficult to optimize the field curvature [START_REF] Bommes | Mixedinteger quadrangulation[END_REF][START_REF] Ray | Nsymmetry direction field design[END_REF]]. More importantly, this representation cannot be trivially generalized to 3D frame fields. For these reasons, we propose to study and fully characterize an alternative representation, where each frame field is represented by the graph of a function with the same symmetries as the frame field.

Functional approach: frame representation

In what follows, the trivial frame f that is aligned with the coordinate axes is called the reference frame. It is defined by f = {(1, 0), (0, 1), (-1, 0), (0, -1)}. We represent the reference frame by the function

F (α) = cos(4α) with α ∈ [0, 2π[ (Figure 3-left).
One can notice that the graph of this function F (α) is similar to the frame f it represents1 . In particular, it exhibits the same symmetries and π/2 rotation invariance as f .

Any frame f can be obtained through a rotation of f by a given angle θ. The functional counterpart is to rotate the graph of the function F , namely any frame f can be represented by a function

F (α) = F (α -θ) = cos(4(α -θ)) with α ∈ [0, 2π[ (Figure 3- right).
Since a rotation of the graph does not change the frequencies of the underlying function, any frame function F (α) can be exactly represented as a weighted sum of Fourier basis functions cos(4α) and sin(4α). Indeed, if F (α) represents a rotation of the reference frame by angle θ, then F (α) = cos(4α-4θ) = cos(4θ) cos(4α)+ sin(4θ) sin(4α).

Let us now introduce some notations: let B = (cos(4α), sin(4α)) denotes the row vector of the orthogonal basis functions and a = (a0, a1) = (cos(4θ), sin(4θ)) denotes the column vector of coordinates of a function F in the basis B. With these notations, any frame function F can be represented by its coordinates a in the basis B, or F = Ba.

A coefficient vector a is feasible if and only if there exists θ such that a = (cos(4θ), sin(4θ)) . Geometrically, a is constrained to lie on a curve parameterized by θ. This curve represents, in coefficient space, all possible rotations of the reference frame. In 2D, this curve is the unit circle, thus the feasibility constraint on a is simply a a = 1.

At this point, we can observe that the coefficient vector a exactly corresponds to the representation vector used in the direction field literature, but it comes with a different viewpoint that will have an importance when moving to the 3D setting in Section 3. Before moving to the 3D setting, we explain how the objective function that expresses the field smoothness can be defined, as well as the constraints it should satisfy.

Functional approach: objective function

The objective function that we minimize is defined as a sum over each edge of the squared differences between the frames at the edges extremities. In our formulation, the difference between two frames (at vertices i and j) is not the rotation angle, but the L 2 distance 2π 0 (F j (α) -F i (α)) 2 dα, which can be simplified as:

E = ij 2π 0 (F j (α) -F i (α)) 2 dα = ij 2π 0 (Ba j -Ba i ) 2 dα = ij (a j -a i ) 2π 0 B Bdα (a j -a i ) = π ij a j -a i 2 (1)
In Equation 1above, the last step is justified because the function basis B is orthogonal (< cos(4α), sin(4α) > = 0). One can also notice that all feasible functions are of norm √ π.

Functional approach: constraints and boundary conditions

As discussed in §2.2, the first constraint is that the variables a i must be feasible (i.e. there should exist a frame represented by a i ).

As for boundary conditions, frames on boundary vertices i must have one member vector equal to the normal direction. If θ i denotes the normal direction, the frame functions on boundary vertices are constrained by two equations:

a i 0 = cos(4θ i ) ; a i 1 = sin(4θ i )

Implementation

At this point, our optimization problem consists in minimizing the objective function E (eq. ( 1)) subject to linear equality constraints on boundary vertices (eq. ( 2)) and quadratic equality constraints a i a i = 1 for each vertex i. We first minimize the quadratic energy without the constraint a i a i = 1, then choose the nearest feasible solution by normalizing a i .

Here, we relax the feasibility constraints, so we need to minimize the quadratic function E subject to linear boundary constraints.

To do that, we simply replace the linear constraints with a strong penalty term in the objective function, leading to a new quadratic function to optimize. This penalty method is very simple and sufficient in practice.

In more details, the new quadratic function is expressed in the form AX -b 2 where A is a matrix, X our variable vector (X2i = a i 0 and X2i+1 = a i 1 ) and b is a vector. The system of equations AX -b = 0 is constructed line-by-line:

• initial objective function E: for each edge ij, we add two equations (eq. ( 1)):

√ π(X2i -X2j) = 0 √ π(X2i+1 -X2j+1) = 0
• boundary constraints: for each boundary vertex i, we add two equations (eq. ( 2)):

CX2i = C cos 4θ i CX2i+1 = C sin 4θ i ,
where the constant C is set as C = 100 in all our experimental results reported here.

From A and b, we just need to solve the linear system A AX = A b to obtain a minimizer of AX -b 2 . Then from X we can obtain an initialization of a i by projecting the corresponding vectors onto the set of feasible coefficients:

a i ← (X2i, X2i+1) / (X2i, X2i+1) .
To solve the linear system, we use the OpenNL library [Lévy 2001], which automatically constructs A AX = A b from the set of linear equations and then solves it by the Jacobi-preconditioned conjugate gradient method. We obtain in the end an algorithm that is very similar to the one in [Knöppel et al. 2013] or to the initialization in [START_REF] Ray | Periodic global parameterization[END_REF][START_REF] Kowalski | A PDE based approach to multi-domain partitioning and quadrilateral meshing[END_REF]]. Smoothness can be further improved by a nonlinear solver initialized with this initial guess, as done in [START_REF] Ray | Periodic global parameterization[END_REF][START_REF] Kowalski | A PDE based approach to multi-domain partitioning and quadrilateral meshing[END_REF] 3 Optimization of 3D frame fields

The representation vector used by previous work was the key to efficiently optimize direction fields in 2D, however, generalizing it to the 3D setting is non-trivial. We now explain how the functional representation viewpoint of the 2D problem ( §2) can be naturally generalized to the 3D setting. We present below the functional representation of 3D direction fields, the constraints they need to satisfy and derive an efficient optimization algorithm.

Definitions for the 3D setting

Given a tetrahedral mesh, our goal is now to create a smooth frame field that is aligned with the boundary of the mesh. We first introduce the following notions:

• The reference frame f is the set of 6 unit vectors aligned with coordinate axes (Figure 4);

• a frame is the reference frame rotated by a 3 × 3 orthonormal matrix R: f = R f ;

• a frame field associates a frame to each vertex of the tetrahedral mesh;

• the boundary constraint ensures that the frame of a boundary vertex has one of its member vectors equal to the normal at the boundary;

• a tetrahedron is said to be singular if any of its triangles is singular. The triangle ijk is singular if and only if R ij × R jk × R ki = Id, where R ij denotes the rotation matrix that brings the frame f i to the frame f j .

• the stable direction of a singular tetrahedron's facet is the set of

2 unit vectors { v, -v} such that R ij × R jk × R ki v = v.
Having a stable direction is a necessary condition to produce hexahedral meshes that [START_REF] Kowalski | A PDE based approach to multi-domain partitioning and quadrilateral meshing[END_REF][START_REF] Jiang | Frame field singularity correction for automatic hexahedralization[END_REF]] are enforcing (on the dual mesh). It is also assumed in the definition of singularities [START_REF] Persson | 23rd international meshing roundtable (imr23) block-structured hexahedral meshes for cad models using 3d frame fields[END_REF]] using local projections.

Frame representation

Like [START_REF] Huang | Boundary aligned smooth 3d cross-frame field[END_REF], we represent a frame by a polynomial, whose restriction to the unit sphere exhibits the 24 symmetries of a cube. The lowest degree of such a polynomial is 4, thus our reference frame (up to a scaling factor and a constant summand) is the polynomial x 4 + y 4 + z 4 restricted to the unit sphere, refer to Figure 4 for an illustration.

In the 2D setting, we decomposed the reference frame polynomial onto the Fourier basis, and our variables were the coordinates in this basis (they correspond to the representation vector). We apply the same technique in 3D: we decompose the function x 4 + y 4 + z 4 onto the basis of spherical harmonics. Note that it is composed of a single component of frequency 4 (we removed the constant summand and normalized the function). Thus, the reference frame f is represented by the function F defined by:

F = 7 12 Y4,0 + 5 12 Y4,4, (2)
where Y l,m is the real-valued spherical harmonic of degree l and order m (see the supplemental material for more details). These harmonics are sometimes called tesseral [Ferrers 1877, p. 74]. The list of harmonics of degree 4 can be found in [Görller-Walrand and Binnemans 1996, p. 239]). The function F is defined over the entire space R 3 → R, but here we are only interested in its restriction to the unit sphere S 2 → R.

An arbitrary frame f can be obtained as a rotation of f by a matrix R. Thus, it is represented by the function F (P ) = F (RP ), where P is a point of the unit sphere (Figure 4).

The family of functions Y l,m forms a function basis over the unit sphere. Clearly, applying a rotation to a spherical harmonic of degree l does not change its degree (it produces another spherical harmonic of degree l). This invariance is quite natural: for example, translating a 1D function does not change its frequencies, thus the degree of functions in the Fourier basis is not changed upon a translation. We represent the reference frame F by a sum of two spherical harmonics of degree 4 (Equation ( 2)), so any frame function F can be exactly represented in the basis B = (Y4,-4, Y4,-3, . . . , Y4,4).

Let us express the reference frame function in this basis: F = Bã with ã = 0, 0, 0, 0, 7 12 , 0, 0, 0,

5 12
, here B is a row vector and ã is a column vector, so F is a scalar function defined on the unit sphere. Any other frame f = R f can be represented as F = Ba, where a = RB ã with RB being a 9 × 9 rotation matrix acting on the coefficients space. The formulas and technical details behind these rotation matrices are also provided in the supplemental material.

A feasible coefficient vector a corresponds to the reference vector ã transformed by a rotation. In other words, it is a vector that can be written as a = RB ã where RB is a 9D rotation matrix that can be derived from a 3D rotation (determined by 3 parameters). Thus, a is geometrically constrained to lie on a manifold of dimension 3 embedded in the 9D coefficient space.

At this point, one can consider the coefficient vector a as an extension of the representative vector used in the direction field literature.

It also corresponds to the representation introduced in [START_REF] Huang | Boundary aligned smooth 3d cross-frame field[END_REF].

Objective function

As in the 2D setting, the objective function is defined as the sum over each edge ij of the squared difference between frames located at the edges extremities. In our formulation, the difference between two frames is the L2 norm of the difference between the corresponding functions, defined as:

E = ij S 2 (F j (α) -F i (α)) 2 dα
Since the function basis B is orthonormal, E can be simplified as:

E = ij a j -a i 2 (3)

Constraints and Boundary Conditions

There are two types of constraints that need to be satisfied by the coefficient vectors: first, each coefficient vector a i needs to be feasible and second, frames on boundary vertices must have one vector aligned with the normal at the volume boundary. We already mentioned the first constraint at the end of §3.2. It is enforced by our optimizer in a dedicated "projection" step (3D counterpart of the normalization of a in 2D case). More details about it are given in the supplemental material.

We now focus on the boundary conditions, which have more particular cases than in 2D, depending whether the considered boundary vertex is smooth, on a crease or on a corner.

Smooth vertex

We assume first that there is only one normal associated with the vertex, it can be computed as the average normal vector of incident boundary triangles.

+ c0(0, 0, 0, 0, 0, 0, 0, 0, 1) + c1(1, 0, 0, 0, 0, 0, 0, 0, 0)

With this equation, all frames with one of their vector that corresponds to z can be represented by the 2D vector c. As in the 2D case, this equation comes with a norm constraint: c 2 0 + c 2 1 = 5 12 . The variable c defines the rotation of the frame around the surface normal. In other words, it defines a 2D frame field. Optimizing this 2D frame field using c as variables is exactly what we did in 2D by introducing the coefficient vector a. Our 3D solution restricted to the object boundary is therefore almost 2 equivalent to our 2D solution.

Case 2: the normal is not equal to the z axis.

In this (more general) case, we rotate the constraint. If we want the vector n to be preserved, we first compute a rotation that brings the z axis to n. From this rotation, we compute the corresponding 9D rotation matrix RB, and derive the constraints: a = 7 12 RB(0, 0, 0, 0, 1, 0, 0, 0, 0)

+ c0RB(0, 0, 0, 0, 0, 0, 0, 0, 1) + c1RB(1, 0, 0, 0, 0, 0, 0, 0, 0)

This expression of the normal constraint gives us a set of 9 linear equations per boundary vertex. It introduces two new variables c0 and c1, and a quadratic constraint c c = 5/12.

Note As in the 2D case, the boundary constraint has a simpler expression [START_REF] Huang | Boundary aligned smooth 3d cross-frame field[END_REF]]: a RB(0, 0, 0, 0, 1, 0, 0, 0, 0) =

Non-smooth vertices

Frames of vertices located on non-smooth vertices (creases, hard edges and corners) have to conform to more than one normal. These vertices have multiple normal constraints, we pick two normals that are almost orthogonal, perturb them (by rotations around their cross product vector) to make them orthogonal, and compute the rotation that brings x to the first normal, and y to the second normal. We compute the corresponding coefficient-space rotation RB and fix the frame coefficient vector a i to RB ã.

Implementation

We now have to minimize our objective function (eq. ( 3)) with linear equality constraints on boundary vertices (eq. ( 5)), quadratic equality constraints c i • c i = 1 on boundary vertices, and feasibility constraints for a i . Algorithm 1 describes our optimization process.

First of all, as in the 2D case (in §2.5), our initialization step is formulated as a least squares problem (minimize AX -b 2 ), constructed without the feasibility constraint on a i (line 1).

Then, the minimizer of the least-squares problem is projected onto the manifold of valid frames (lines 2-4). The projection f i ← closest frame(a i ) onto the set of feasible coefficient vectors is no longer a simple normalization of representation vectors, as it was in 2D. Instead we perform, for each vertex, a gradient descent. In more details, starting with a seed frame, we rotate it gradually in order to minimize the distance between the current frame function and the function to be projected. The gradient is evaluated by calculating the variation of the L 2 norm induced by infinitesimal rotation matrices with Euler's angles. We do not have a formal proof, but we conjecture that there is a single minimum of the L 2 norm.

Finally, starting with the initial field, we smooth it using L-BFGS (lines 5-7). Since at this step we already have valid frames, we use Euler angles (rotation with respect to the reference frame) as variables in the L-BFGS optimization procedure. On the boundary of the volume each frame is represented by one angle representing a rotation around the surface normal. Note that each frame can be represented by 48 triplets of equivalent Euler angles. In our implementation we choose the triplet that maximizes the distance to the nearest gimbal lock. An extensive description of all steps is provided in the supplemental material.

Algorithm 1: Frame field optimization Input: A tetrahedral mesh M with nv vertices including n l vertices with normal constraint Output: A frame f i for each vertex i 1 {a i } nv i=1 ← arg min E({a i }) subject to normal (5), but not feasibility constraints ; // solve a linear system 2 foreach i < nv do 3 f i ← closest frame(a i ); // project onto feasible manifold 4 end 5 repeat 6 L-BFGS iteration on E({f i }) subject to feasibility and normal constraints 7 until ∆E/E < 10 -5 or time exceeds;

Results

First, we observe that minor modifications of the methods introduced in previous work can significantly improve the result. We obtained much better results than Huang et al. by sampling the field on vertices instead of tetrahedra facets. It does not avoid degenerate cases, as shown if Figure 9, but better captures the global shape of the object. The running time also significantly decreases. So, this section does not only compare our algorithm to previous work, but also evaluates benefits of our initialization alone and the impact of sampling on vertices (instead of tetrahedra or tetrahedra facets).

The main challenge here is to find a fair way to evaluate the field quality: a natural idea would be to evaluate and compare the values of the objective function (eq.3). However, we noticed that in the global value of the objective function, singular regions of the field have a dominating influence. As a consequence, a lower value does not mean that the corresponding field is better suited to hexdominant meshing. We illustrate this remark in Figure 5, which shows an example where a field with higher objective function produces a considerably lower proportion of non-hexahedral elements.

A better option would be to compare the quality of the mesh produced from these fields, but the state-of-the art algorithms [Nieser et al. 2011;[START_REF] Kowalski | A PDE based approach to multi-domain partitioning and quadrilateral meshing[END_REF][START_REF] Jiang | Frame field singularity correction for automatic hexahedralization[END_REF] are not robust enough for testing a large collection of objects. Moreover, it is difficult to separate the impact of the meshing scheme from the underlying field quality. We propose an alternative solution that estimates how the field's singularity graph is suitable for hex remeshing, and compare it in Table 1. The next section describes the protocol that we used. In addition, in the supplemental material, we show the hexahedra produced by CubeCover applied to the original mesh minus the singular tetrahedra. It gives an intuition about what could be expected from a hex-dominant meshing algorithm.

After a brief description of our experimental protocol, we present the results and compare the methods.

Experimental protocol

We cannot predict whether it is possible to produce a nice hex mesh from a given frame field. However, we know a necessary condition: [Nieser et al. 2011, §2.2] have proven that a singularity of a volume parameterization is always similar to a 2D singularity extruded along the third coordinate. In our setting, it means that the singularities must admit a stable direction (refer to § 3.1 for the definitions of singularities and stable directions) and it must be equal to this third coordinate. Note that only the existence of a stable direction was enforced by local operations in [START_REF] Kowalski | A PDE based approach to multi-domain partitioning and quadrilateral meshing[END_REF]].

According to this observation, we derive the following measure of field quality: first, we extract the singularity graph i.e. for each singular triangle, we create a segment linking the barycenters of its adjacent tetrahedra, then we smooth the geometry of this graph to remove the high-frequency noise that is due to the discretization. Next, for each segment, we determine the stable direction of its corresponding singular triangle. If it exists and its deviation from the segment direction is smaller than π/4, then we consider that it can be used for hex-meshing. Otherwise this part of the singularity graph is considered to be incompatible with hex meshing (see Figure 6 for an illustration). Column 1 is the model identifier used in the supplemental material, column 2 gives the model size (in thousands of tetrahedra), columns 3,4,5,6,7,8 report the total length of the singularity graph that is incompatible with hex remeshing (refer to §4.1), columns 9,10,11,12,13,14 report running times in seconds. We stop L-BFGS iterations when ∆E/E < 10 -5 .

evaluate the impact of the sampling, we report the statistics of previous works with the initial sampling (on tetrahedra) and with our sampling (on vertices). Note that for Huang et al.'s algorithm, our sampling also introduces non-smooth vertices ( §3.4), which completely fix the drifting problem encountered on hard edges. The comparison is made on a database of canonical geometrical shapes and industrial models (see supplemental materials for images). It includes a selection of failure cases for each algorithm: objects of revolution(p1, p6, p10, p17, p27), object including long cylindrical shapes (p12, 27), degenerated case for Huang's boundary condition (p4, p5), and "traps" for the front propagation initializer (p16, p17, p20, p21, p23, p26, p27, p29, p34).

Observations

Our initializer is very fast, easy to implement, and produces fields that are often equivalent to other methods. There is one exception: for objects of revolution, our initializer fails to pick a solution in the infinite set of equivalent solutions. Our method fixes these failure cases with the non-linear solver. On other models, the non-linear solver improves the quality a little bit.

Previous works are slower and have failure cases (see below). However, if we improve them by sampling on vertices instead of sampling on tetrahedra, the nonlinear solver is often able to produce fields with a quality that is similar to ours (see Figure 6). It just takes more time to converge from a weak initialization (see Figure 7).

For dense input tet meshes, even when sampled on vertices, L-BFGS is unable to reach a good configuration from a weak initialization. Such cases are illustrated in Figure 8 with failure cases of the initialization for each algorithm.

Huang et al.'s original method suffers from the boundary condition that is only partially enforced during the initialization step. It can lead to dramatic consequences, as illustrated on some simple shapes (Fig. 9), where a constant field perfectly satisfies their incomplete boundary condition whereas it does not match the bound- ary. One can notice that, in such cases, the harmonic function is also very far from being valid.

For less specific objects, results presented in their article are all "polycube like". We preferred to test their algorithm with a higher penalty term associated to the boundary condition: it produces more interesting fields, but also introduces new failure cases as illustrated in Figure 10.

The improved version of their method samples the field on vertices and introduces locked frames that avoid drifting on hard edges. It Even when sampled on vertices, the nonlinear solver fails to find a good singularity graph from a poor-quality initial field. We used the following initializations: A: improved Huang's method applied to a tetrahedron model (865K tets), B: improved Li's method applied to a higher resolution version of object p26 (2639K tets), C: our initialization on an object of revolution -ellipsoid (2713K tets), D: Li's method on the same ellipsoid. is therefore very similar to our algorithm.

Fields generated by Li et al.'s original algorithm match the boundary, however they are likely to have a singularity graph that is incompatible with hexahedral remeshing. When their initialization is too much distorted across the medial axis of the object, the nonlinear smoother is not always able to find a good quality singularity graph (Fig. 11). The modified version sampled on vertices is better but the lower quality of the initialization makes the whole algorithm slower than ours (Fig. 7).

Conclusion

In this paper, we draw a bridge between existing 2D and 3D frame field generation algorithms, and use it to extend some 2D frame field design algorithms to 3D. As a practical consequence, the ob- tained algorithm computes a better initial solution than previous work in 3D, and can be easily extended with additional constraints and/or internal boundary conditions. The initial solution computed by our algorithm (i.e. without the smoothing step) can be sufficient for some applications and has a virtue of being simple to implement: it requires only a sparse linear system solver.

The smoothness of this initial field can be further improved by a modification of [START_REF] Kowalski | A PDE based approach to multi-domain partitioning and quadrilateral meshing[END_REF]] that samples the field on vertices. The discretization on vertices that we propose makes L-BFGS much more efficient than in the original version of [START_REF] Kowalski | A PDE based approach to multi-domain partitioning and quadrilateral meshing[END_REF] where the field was sampled on tetrahedra. With this approach we are able to generate (on a laptop) 3D frame fields on models with up to several millions tetrahedra in a matter of minutes.

To ease the reproducibility of our results, extensive explanations with all the derivations and algorithms are provided in the supplemental material.

Figure 2 :FFigure 3 :

 23 Figure 2: A 2D frame is a set f of 4 vectors f0, f1, f2, f3 invariant under rotation by π/2. Its angle representation is the rotation θ between the global axis x and f0. F (α) f F (α) = F (αθ)

Figure 4 :

 4 Figure 4: An arbitrary frame f corresponds to the reference frame f rotated by a 3 × 3 matrix R. The graphs of the corresponding functions F and F are also rotated by R, and their coefficients vectors verify a = RB ã where RB is a 9 × 9 rotation matrix (derivation detailed in the supplemental material).

Case 1 :

 1 the normal is equal to the z axis.Let us first consider the case where the fixed vector (the surface normal) corresponds to the z axis. If we rotate F around z by an angle θ, we obtain a = RB ã with RB being a rotation around z. The simple structure of RB together with the zero coefficients of ã gives the equation: the construction of coefficient vectors in the 2D case, we can get rid of the angle θ by replacing it by a vector c = (c0, c1)

Figure 5 :

 5 Figure 5: The frame field initialized by Li et al.'s front propagation (left) has a better smoothness energy than ours (right). However it yields a lower proportion of hexahedra. Therefore, the quality of the singularity graph is more important than the smoothness energy.

Figure 6 :

 6 Figure 6: Visualization of results for model p31 in Table 1. Top: hexahedra extracted from CubeCover. Bottom: singularity graphs, with the hex-incompatible parts colored in red. A: our algorithm; B: Huang's algorithm; C: Li's algorithm; D: our initialization alone; E: our modified version of Huang's algorithm; F: our modified version of Li's algorithm.

Figure 7 :

 7 Figure 7: The upper row shows the singularity graphs obtained by Li et al. sampled on vertices after 1 minute (left), after 3 minutes (middle) and at convergence after 20 minutes (right). The lower row shows our result at the same time steps.

Figure 8 :

 8 Figure8: Importance of the initialization for high-resolution input meshes. Even when sampled on vertices, the nonlinear solver fails to find a good singularity graph from a poor-quality initial field. We used the following initializations: A: improved Huang's method applied to a tetrahedron model (865K tets), B: improved Li's method applied to a higher resolution version of object p26 (2639K tets), C: our initialization on an object of revolution -ellipsoid (2713K tets), D: Li's method on the same ellipsoid.

Figure 9 :

 9 Figure 9: Typical failure cases of Huang's method are objects where a constant field of SH functions respects all boundary conditions of the model.

Figure 10 :

 10 Figure 10: Impact of penalty term for Huang's normal constraint: high value (left) avoids "polycube like" fields (right), but introduces degenerate cases (sphere).

Figure 11 :

 11 Figure 11: Typical failure cases of Li's method are object with high field discontinuities across the medial axis.

Table 1

 1 compares the total length of hex-incompatible singularities and the timings for our algorithm, Huang et al.'s algorithm and Li et al.'s algorithm. We also evaluate the potential of using our initialization alone, without the subsequent nonlinear smoothing steps, as a much simpler and faster alternative to the overall algorithm. To

	Model	Size			Incompat. singularities						Time (s)
			our method	Huang et al.	Li et al.	our method	Huang et al.	Li et al.
		(K tets)	full init only original modified original modified full init only original modified original modified
	p1	6	0	2	0	0	0	0	1	0	4	33	1	2
	p2	7	1	0	58	1	3	1	0	0	9	22	1	0
	p3	11	0	0	3	0	4	0	0	0	10	25	2	0
	p4	12	0	0	0	0	1	0	1	0	33	66	2	1
	p5	17	2	2	117	2	3	2	1	0	33	174	3	2
	p6	17	0	0	0	0	0	0	1	0	20	31	2	2
	p7	18	8	16	21	10	12	6	2	0	58	198	4	3
	p8	21	0	0	1	0	0	0	1	0	45	34	3	1
	p9	22	0	5	11	0	9	0	2	0	21	72	8	3
	p10	27	17	293	154	12	47	19	9	0	27	278	52	6
	p11	34	71	125	136	78	61	85	2	1	64	248	28	4
	p12	40	21	13	50	24	34	32	5	1	62	123	38	11
	p13	46	85	170	200	93	145	101	7	1	71	246	72	9
	p14	61	0	1	28	0	1	0	5	1	295	419	29	7
	p15	69	89	113	69	107	36	94	3	1	184	280	54	4
	p16	112	53	101	127	53	103	50	19	3	513	623	50	18
	p17	121	25	326	154	18	110	33	30	1	183	606	278	57
	p18	130	61	104	162	68	128	72	12

Table 1 :

 1 Evaluation of our algorithm (with and without the smoothing step), and Huang et al.'s algorithm (sampled on faces and sampled on vertices) with a high penalty given to the boundary condition, and Li et al.'s algorithm (sampled on tetrahedra and sampled on vertices).

Note that representing frames by functions is not completely new: when expressed in Cartesian coordinates, our reference frame function F is the polynomial 4(x 4 + y 4 ) -3 restricted to the unit circle, thus this formulation is equivalent to traceless symmetric 4 th order tensors (in other words, degree 4 polynomials) manipulated in[Palacios and

Zhang 2007b].

12 that is valid only if all a i are feasible, which is not ensured by their method. As a consequence, it cannot be safely used during the initialization step. We discuss the practical consequences later, and show a comparison in Figure9.

The boundary has curvature that is not present in the (flat) 2D setting.


Chen Chen

Adrien Ugon

Xun Zhang

Amara Amara

Patrick Garda

Jean-Gabriel Ganascia

Carole Philippe

Andrea Pinna

Personalized Sleep Staging System by combining Symbolic Fusion and Feedback System Control

In this paper, a personalized sleep staging system is proposed by combining symbolic intelligence and feedback system control technique. Symbolic Fusion is dedicated to mimic decision-making process of clinical sleep staging. It starts from extraction of digital parameters from raw polysomnography (PSG) signals and goes up-to high level symbolic interpretation of feature parameters. At last, decision is generated using rules inspired by international guidelines in sleep medicine and applied to feature parameters. Meanwhile, Feedback System Control (FSC) is designed to provide optimal thresholds for Symbolic Fusion in transforming digital parameters into feature parameters which can take individual variability into consideration. Proposed system can be also integrated with portable PSG system to realize personal sleep evaluation or home care application in the future.

Introduction

Sleep is an indispensable part in daily life which contributes to self-repairing and self-recovering. However, sleep disorders are affecting more and more people. It becomes a significant cause of morbidity and mortality [START_REF] Jamil | The prevalence, cost implications, and management of sleep disorders : an overview[END_REF].

Sleep Staging, as a fundamental step to diagnose and treat sleep disorders, is usually visually scored by an expert according to American Academy of Sleep Medicine (AASM) manual [START_REF] Iber | The AASM manual for the scoring of sleep and associated events : rules, terminology and technical specifications[END_REF] in analysis of an overnight PSG recordings of electroencephalography (EEG), electrooculography (EOG) and electromyography (EMG). However, it is a time-consuming task and inter-rater variability exists. 82% inter-rater reliability is reported for visual sleep staging performed by different experts [START_REF] Danker-Hopfe | Interrater reliability for sleep scoring according to the rechtschaffen & kales and the new aasm standard[END_REF]. Due to the limitations of clinical sleep staging, automatic sleep staging methods have gained a wide spread attention of researches.

In the past few decades, different automatic sleep staging methods have been proposed. Without the requirement of complex clinical knowledge, numerical methods received much attention involving Decision Tree, Artificial Neural Network and Support Vector Machine. However, it suffers from several severe limitations: 1) Numerical classification is often thought as an independent process ; whereas sleep staging is a time dependent classification problem. The stage scored at one position can be influenced by the stage assigned to previous or/and to next position. 2) Smaller number of signals has been involved in order to reduce the dimension of the feature sets. However, to increase the accuracy and reliability of the sleep staging, there is a need to include more signals in the analysis [START_REF] Clete | Practice parameters for the indications for polysomnography and related procedures : an update for 2005[END_REF].

A Symbolic Fusion method was proposed in [START_REF] Ugon | Fusion symbolique et données polysomnographiques[END_REF] to overcome the limitations of numerical methods. It used AASM manual as guidances to imitate clinical sleep staging process and allows providing enhanced and complementary decision by combining EEG, EOG and EMG. Symbolic fusion starts from extraction of digital parameters from raw PSG signals and goes up-to high level symbolic interpretation of different features. However, this approach is limited by the manual setting-up of thresholds in transforming digital parameters to symbolic interpretation, which require manual efforts.

In this paper, a personalized system combining FSC and symbolic fusion is proposed. FSC is adopted for each subject to avoid inter-subject variability and provide optimal thresholds to symbolic fusion.

Personalized Sleep Staging System

In 2003, an automatic sleep staging system was proposed which adopted the three-level architecture defined by B.Dasarathy [START_REF] Belur | Sensor fusion potential exploitation-innovative architectures and illustrative applications[END_REF]. As shown in Figure 1a, it consists of data fusion, feature fusion and decision fusion. In data fusion, digital parameters were extracted from PSG recordings to maximize the useful information and to minimize noise and artifacts. In feature fusion, digital parameters were transformed into feature parameters. It simplifies the interpretation of digital parameters, and also performs normalization, reduction and matching of digital parameters. In decision fusion, inference method was used to fulfill sleep staging on the basis of feature parameters.

FIGURE 1. Sleep Staging System

However, visual interpretation was performed to define thresholds in [START_REF] Ugon | Fusion symbolique et données polysomnographiques[END_REF]. These thresholds are used in transforming digital parameters into feature parameters to differentiate the boundary of different feature parameters. E.g. According to AASM manual, chin EMG tone is low in stage Rapid Eye Movement (REM). The boundary for low or high feature of EMG is determined by a visual setting of threshold for digital parameter. Total 15 thresholds are defined for translating nine digital parameters into 24 feature parameters. It is a time-consuming process and individual variability also exists.

To overcome the limitations in [START_REF] Ugon | Fusion symbolique et données polysomnographiques[END_REF], a personalized sleep staging is proposed, as shown in Figure 1b. A Thresholds Selection part using Feedback System Control technique is involved. A generic FSC platform was proposed by Chih-Ming Ho [START_REF] Ho | Control of complex systems : Fsc.x technology enabled personalized medicine[END_REF] to find optimal input parameters combination for guiding the bio-system toward to a desired state in 2012. With robustness and rapid searching ability, FSC has been demonstrated to be very effective in combinatorial parameters selection problems.

In this paper, a specific FSC for sleep staging system is proposed. As shown in Figure 2, it is dedicated to provide optimal automatic thresholds selection for releasing the burden of visual interpretation and avoiding individual variability. Digital parameters extracted from PSG are used as the inputs of Thresholds Selection. Meanwhile, initial thresholds combinations are selected by Stochastic Searching Algorithm (SSA). Based on the initially selected thresholds combinations and digital parameters, sleep staging is performed by combining feature fusion and decision fusion. To evaluate the accuracy of sleep staging, a Assessment of Sleep Staging (ASS) part is adopted. With the assessment result of sleep staging, SSA will generate a new thresholds combination for the next loop. This process is repeated until optimal thresholds combination is found which can provide good assessment in ASS. With this Thresholds Selection, an optimal threshold combination can be generated automatically for each subject.

Conclusion

In this paper, a personalized sleep staging system is proposed by combining Symbolic Fusion and Feedback System Control technique. With the feasibility of FSC technique and the elasticity of symbolic intelligence, proposed system can be a reliable computer-assisted tool for assisting clinical sleep analysis and can be integrated with other medical systems for remote sleep monitoring.

FIGURE 2 .

 2 FIGURE 2. Thresholds Selection using Feedback System Control


A Jardin

D Mazon

M O'mullane

M Chernyshova

T Czarski

K Malinowski

G Kasprowicz

A Wojenski

K Pozniak

P Malard

C Bourdelle

Synthetic SXR diagnostic using GEM detectors on WEST: development in the prospect of tungsten monitoring

Introduction

WEST (Tungsten (W) Environment in Steady-state Tokamak) will start operating in 2016 as a test bed for ITER divertor components in long pulse operation. In this context, radiative cooling of heavy impurities like tungsten (W) sputtered from Plasma Facing Components (PFC) into the plasma core is a critical issue for the plasma performances and stability. Indeed, even small W concentrations such as 3. Finally, work is ongoing to apply the synthetic diagnostic on WEST cases with W impurities.

A simple model of GEM response to SXR radiation

The GEM diagnostic is structured as follows.

Photoionization occurring for each detected X-ray photon in the GEM gas mixture Ar-CO 2 (70-30%) produces an electron cloud in the conversion chamber.

This electron cloud drifts and diffuses in the gas toward the anode due to electric fields 𝐸 𝐼,𝑇,𝐶 ⃗⃗⃗⃗⃗⃗⃗⃗⃗ of a few kV/cm, see Fig. The electric field ~ 50-100 kV/cm in GEM holes of r 1 -r 0 ~ 50µm is high enough, see Fig. 4, to trigger electron avalanches by collision with gas atoms:

where 𝐴 ̅ > 1 is the amplification factor, α is the first Townsend coefficient and 𝐺 ̅ ~ 10 2 -10 4 is the average gain. The associated relative standard deviation is determined by: 𝑛(𝑥, 𝑦, 𝑧, 𝑡) = ( 1

( 𝜎 𝐺 𝐺 ̅ ) 2 = ( 𝜎 𝑁 𝑁 ̅ ) 2 + 1 𝑁 ̅ ( 𝜎 𝐴 𝐴 ̅ ) 2 = 𝐹 + 𝑓 𝑁 ̅ = (𝐹 + 𝑓). 𝑤 ℎ𝑣 (
√4𝜋𝐷 𝐿 𝑡 ) 2 ( 1
√4𝜋𝐷 𝑇 𝑡 ) 𝑒𝑥𝑝 (- 𝑥 2 + 𝑦 2 4𝐷 𝑇 𝑡 - (𝑧 -𝑣 𝑍 𝑡) 2 4𝐷 𝐿 𝑡 ) (3)
The that the deflection of electron trajectories (so called the Lorentz angle) due to the ExB drift should not exceed ~5° in the worst case, well below the deflection limit of 1 pixel (~10°).

Furthermore, an additional µ-metal shielding will decrease the magnetic field below 0.1 T to avoid any perturbative effect on the GEM.

Tests with a Fe 55 radioactive source

The synthetic diagnostic is applied to a radioactive Fe 55 source (5.9 keV) with a 2D square GEM. Fig. 5 shows the results of the simulation, including the charge spatial distribution on the GEM anode (left) and the Fe 55 spectrum acquired by the GEM (right) with its Argon Two Fe 55 sources will be implemented in WEST SXR cameras in order to allow an in situ spectral calibration of the GEM before and after every plasma discharge.

Perspectives for tungsten monitoring on WEST

WEST will be the first tokamak to use the GEM technology as SXR diagnostic for plasma tomography, using a Minimum Fisher algorithm [START_REF] Jardin | proc. ICFDT[END_REF]. Preliminary ongoing studies on the standard WEST Physics basis scenario with heating power of 12 MW and I p = 0.6 MA [START_REF] Bourdelle | [END_REF] are encouraging in terms of incoming photon flux (10 5 -10 7 ph/strip/s), tomographic reconstruction capabilities and observation of W poloidal asymmetries. Furthermore, W line radiation around ~2 keV might be disentangled from the continuum background plasma radiation ≳ 4 keV by choosing appropriately the GEM energy bands.

  10 -5 increases by 20% the minimum triple product nTτ E required to make the thermonuclear burn possible [1], and can sometimes lead to radiative collapse. Thus reliable tools are needed to monitor W density and avoid its accumulation in the plasma core, where its radiation is dominant in the Soft X-ray (SXR) range 0.1 keV -15 keV with complex contributions from line transition, radiative recombination and Bremsstrahlung emission. The SXR diagnostic of WEST will be equipped with two new GEM (Gas Electron Multiplier) based poloidal cameras [2], [3] allowing to perform 2D tomographic reconstructions in a poloidal cross-section [4], see Fig.1, and with spectral resolution in tunable energy bands, in the contrast with former Tore Supra silicon barrier diodes working in current mode. Thus once the GEM response to plasma emissivity is characterized thanks to a synthetic diagnostic, it will offer new possibilities [5] to disentangle the different SXR contributions in the prospect of W monitoring in harsh fusion environments

Figure 1 .

 1 Figure 1. Simulated SXR emissivity of WEST in a poloidal cross-section and edge Lines of Sight (in red) of the two GEM cameras.

  photoionization of an Argon atom, see[7] and Fig.3, releasing a free electron of energy E ≃ hv -E b , where E b is the binding energy of the electron. This energetic electron creates an electron cloud by collision with neighboring gas atoms, resulting in 𝐸 𝑤 ⁄ electron-ion pairs, where the mean ionisation energy w(Ar-CO 2 ) ≈ 28 eV. The residual energy E b in the photoionized Argon atom can be released by Auger emission such that the total average number of electron-ion pairs 𝑁 ̅ ≃ ℎ𝑣 𝑤 ⁄ ~ 10 2 . For hv ≳ 3.2 keV there is a probability of 0.14 to have instead emission of a 2.9 keV photon by Argon K fluorescence [8] that can escape the detecting volume. This missing energy leads to a parent Argon escape peak in the GEM signal.

Figure 2 .Figure 3 .

 23 Figure 2. Layout of the GEM principle

2 )

 2 where F (~ 0.19 for Argon) is the Fano factor [9] and f ~ 0.6 is the electron avalanche standard deviation. This corresponds to a loss of spectral resolution of ~ 5-10% per GEM foil.The spatial electron distribution n(x,y,z,t) is determined by solving the 3D advection-diffusion equation [10], in presence of electric field 𝐸 ⃗ = 𝐸 𝑧 𝑢 𝑧 ⃗⃗⃗⃗ , and applied on the electron cloud.

  drift velocity v z and diffusion coefficients D L,T of the electron cloud are determined using the program Magboltz [6], see Fig. 4. Magboltz solves the Boltzmann equations to study electron avalanche and transport in a given gas mixture under influence of electric and magnetic fields. The magnetic field created by WEST magnetic coils will be up to 0.3 T at the GEM location and might interfere with the electron cloud transport. However, Magboltz simulations have shown [11]

Figure 4 .

 4 Figure 4. Magboltz parametrization of electron avalanche and transport coefficients at B = 0 T

Figure 5 .

 5 Figure 5. 2D GEM simulated response to a Fe 55 radioactive source (5.9keV) with the charge spatial distribution (left) and the charge value distribution (right).