 # On Smooth 3D Frame Field Design

We analyze actual methods that generate smooth frame fields both in 2D and in 3D. We formalize the 2D problem by representing frames as functions (as it was done in 3D), and show that the derived optimization problem is the one that previous work obtain via "representation vectors." We show (in 2D) why this non linear optimization problem is easier to solve than directly minimizing the rotation angle of the field, and observe that the 2D algorithm is able to find good fields. Now, the 2D and the 3D optimization problems are derived from the same formulation (based on representing frames by functions). Their energies share some similarities from an optimization point of view (smoothness, local minima, bounds of partial derivatives, etc.), so we applied the 2D resolution mechanism to the 3D problem. Our evaluation of all existing 3D methods suggests to initialize the field by this new algorithm, but possibly use another method for further smoothing.

## Authors

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

In computer graphics, a frame field can be defined on a surface () or inside a volume (). For each point of the domain, it defines a set of (resp. ) unit vectors invariant by rotations of around the surface normal (resp. around any of its member vector). The main motivation to study these fields is to split the quad and hexahedral remeshing problems into two steps: (1) the design of a smooth frame field, (2) and the partitioning of the domain by quads or hexes aligned with the frame field. Our objective is to unify the formulation of the and frame field design problems, and to use it to extend an efficient algorithm to the case.

In most cases, frame field design is formalized as the following optimization problem: find the smoothest frame field subject to some constraints. What makes them different from each others is obviously the dimension of the frames ( or ), but also the definition of the field smoothness, the expression of the constraints, and the optimization method. Interestingly, the case and the case are addressed by very different strategies:

• In , the frame field design problem can be restated as a vector field design problem thanks to the introduction of the “representation vector”. In local polar coordinates, each vector of a frame has the same angle modulo , if we multiply it by we obtain a unique representation vector (modulo ). It is easy to derive optimization algorithms acting on the representation vectors. For simplicity reasons, we limit ourselves to planar frame fields and use the algorithm proposed by Kowalski et al. [Kowalski et al. (2012)] as reference.

• In , it is not possible to extend the idea of “representation vector”. Instead, Huang et al. [Huang et al. (2011)] propose to represent frames by functions defined on the sphere, refer to figure 1 for an illustration. A definition of the field smoothness is derived from this representation and optimized in a two step procedure: (1) initialization based on optimization of spherical harmonics coefficients in [Huang et al. (2011)] or front propagation of boundary constraints in [Li et al. (2012)], followed by (2) smoothing iterations performed by L-BFGS on Euler angles representation of frames. Figure 1: Orthonormal frames (close-up (a)) are represented by spherical harmonic functions (close-up (b)), attached to each vertex of a tetrahedral mesh. Streamlines and singularities of the field are shown in yellow and red, respectively.

Thus our goal is to better understand how and problems are related to each other and to extend [Kowalski et al. (2012)] to . We first show that [Kowalski et al. (2012)] can be derived with the formalization approach inherited from the case, and then we extend it to by the same logical flow. Busy readers interested in only reproduction of the method can skip to implementation section §3.5, the only required tool is a linear solver, all calculations are given explicitely.

### The 2d algorithm with frames represented by functions

Solutions developed for are very different from solutions because the “representation vector” trick does not extend nicely into . To unify both problems, we propose to go in the other direction §2: we apply the functional frame representation to the case. By doing so, we found another way to introduce the “representation vectors”: they appear as coefficient vector of the function decomposed in the Fourier function basis §2.2. Following the logical flow introduced for the

case, we derive an estimation of the field smoothness §

2.3 and formulate the corresponding optimization problem §2.6. We end up with exactly the formulation of methods based on the “representation vectors”.

This common formulation of the problem is not the simplest one: the difference between adjacent frames is not evaluated by the rotation between them as in [Bommes et al. (2013), Ray et al. (2008)], but approximated by the Euclidian distance in a plane. We observe the impact of this approximation on the objective function §2.7, and show how it simplifies the optimization problem §2.8.

The resolution mechanisms of the optimization problem are strongly inspired by the geometric intuition of the representation vector field. To make abstraction of this intuition, we re-explain the algorithm from [Kowalski et al. (2012)] with the notations/vocabulary introduced by the functional representation of frames.

### Extension of the 2d algorithm to 3d

Now that we can describe the algorithm without referring to the representation vector, it is possible to extend it to . Using the notations introduced in the case, we describe the version in §3 and extend the optimization mechanism to §3.5.

A first difficulty was to find the expression of the boundary conditions because the boundary frames are free to rotate around the surface normal §3.4. Incomplete enforcing of this condition by Huang et al. [Huang et al. (2011)] results in a poor initialization of the optimization procedure as evaluated in §3.6.1.

The second difficulty is that the frame is represented in a function basis, but the set of functions that corresponds to a frame has dimension (the rotation that brings the axis aligned frame to the current frame). The extension of the normalization of the representation vector in [Kowalski et al. (2012)] becomes: find the nearest point on the manifold of the set of functions representing a frame.

Our extension of [Kowalski et al. (2012)] nicely completes the tool set of frame field design algorithms §3.6. Our initial field has lower energy, often a better topology, and is more robust to surface noise. The optimization step is easy to implement from the initialization step, is competitive when the initial topology is good enough, but does not outperform the current state of the art.

### Previous works

The orientation of objects is commonly represented by symmetric tensors in physics to model the orientation distribution of fibers. For example, Moakher

[Moakher (2009)] introduced the notion of cubic orientation distribution functions. However, in computer graphics, more compact representations are often preferred: “representation vectors” in and vectors of spherical harmonics coefficients in .

#### On surface

The optimization of a frame field inside a shape is very similar to the optimization of direction fields on surfaces 111The differences between these two problems (angle defect, hard constraints, curvature fitting term, etc.) have an impact on the optimization algorithm, as detailed in the supplemental material.

The problem of direction field design on surfaces was introduced by [Hertzmann and Zorin (2000)] for orienting strokes in non photo realistic rendering. Directions are represented by an angle rotation per vertex, and the smoothing is performed by a non linear solver (BFGS). The ”representation vector” was not introduced yet, and the optimization mechanism was similar to actual approach of frame fields smoothing. Solving with a representation vector for each was used in [Ray et al. (2006)] for faster results, and improved later for better control over the field topology [Ray et al. (2009)]. Based on the similar representative vectors, [Palacios and Zhang (2007a)] propose to control the field topology by local operations. For direction fields without constraints, [Knöppel et al. (2013)]

proved that the norm of the representative vector does not affects the result, leading to finding optimal direction fields by solving an eigenvector problem.

Directly working with angle allows to perfectly control the field topology, but at the expense of solving a mixed integer system [Ray et al. (2008), Bommes et al. (2009)]. In [Palacios and Zhang (2007a)], representative vectors are introduced with its duality with order traceless symmetric tensors. This relation is very useful to unify and frame fields.

#### Inside a volume

The pioneer work of [Huang et al. (2011)] discovered that each frame can be represented by a spherical harmonic. Their initialization step defines a smooth spherical harmonic field then, for each sample, defines the frame that better aligns with the spherical harmonic. This initial solution is then improved by rotating each frame. These rotations are defined by Euler angles and obtained by minimizing the field smoothness with L-BFGS and enforcing the boundary alignment by a penalty term.

Li et al. Li:2012 propose an alternative initialization method. They optimize a frame field on the volume boundary, convert it to a frame field by adding the surface normal and its opposite, then propagate it inside the volume. The resulting field is then smoothed by optimizing a rotation for each frame, as in [Huang et al. (2011)], but with an improved optimization scheme. They also optimize the singularity graph of the field by local combinatorial operations, as done by [Jiang et al. (2014)].

Our extension of [Kowalski et al. (2012)] is optimizing for the same objective function, but with a very different solution mechanism. It performances are compared with previous work in §3.6.

## 2 Functional representation of frames in 2d

This section introduces how to optimize frame fields using a functional representation for each frame. While we do not claim any technical contribution in this section, we think that it is important to reformulate existing concepts using the functional representation, because case inherits exactly the same difficulties and the intuition we gained in helps to motivate the choices made for fields. We derive an energy and the boundary conditions from this new representation. The resulting optimization problem is exactly the one usually solved by direction field algorithms based on the standard “representation vector”. We show on a toy and a real examples that this optimization problem is much easier than directly minimizing the frame rotation. We then present the algorithm [Kowalski et al. (2012)] that we extend to in the following section.

### 2.1 Problem settings

Given a shape, frame field design in consists of 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 222The problem is directly presented in discrete settings. Interesting results on its (non trivial) continuous counterpart are given in supplemental materials.:

• A frame is a set of unit vectors that is invariant by a rotation of (Figure 2). It can be represented by the angle such that .

• A frame field is a frame per vertex of a shape triangulation.

• The boundary constraint: a frame located on a boundary vertex must have one of its member vectors equal to the normal on the boundary.

• The rotation angle between two frames is the angle of the rotation that transforms one frame into the other. This angle being defined modulo , we choose the with minimum absolute value.

• The curvature of a frame field is the sum over each edge of the squared rotation angle between frames that are defined at the edge extremities.

• A triangle is said to be singular 333Frame field topology is further discussed in supplemental material. if the sum of the rotation angles over his boundary is not equal to . 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.

Representing frames by angles is simple, but it makes the field curvature hard to optimize [Bommes et al. (2009), Ray et al. (2008)], and it does not extend nicely in . For these reasons, we propose an alternative representation based on functions, and use this curvature based formalization as a reference.

### 2.2 Functional approach: frame representation

The frame aligned with coordinate axes is called the reference frame . Instead of using the angle approach, we represent it by the function with (Figure 3–left), that exhibits the same rotation invariance as the frame.

Any other frame can be obtained by a rotation of by angle . The functional counterpart is to rotate the graph of the function , namely any frame can be represented by a function with (Figure 3–right).

A compact representation of these functions is given by the trigonometric relation : we see that a frame function can be represented by a vector of coefficients in Fourier basis i.e. .

A coefficient vector is feasible if and only if there exists such that . Geometrically, is constrained to live on a curve parameterized by . This curve represents, in coefficient space, all possible rotations of the reference frame. In , this curve is the unit circle, so the constraint on is simply : .

At this point, we can observe that the coefficient vector is exactly the representative vector used in the direction field literature. We can also notice that expressed in Cartesian coordinates, our reference frame function is the polynomial restricted to the unit circle, thus it is also equivalent to the traceless symmetric order tensors manipulated in [Palacios and Zhang (2007b)]. Figure 3: Parametric plot of the reference frame representation ~F(α) (left) and an arbitrary frame F(α) (right). The plot is made using x(α)=(1+F(α))cos(α) and y(α)=(1+F(α))sin(α) for α∈[0,2π[. It is easy to see that corresponding frames are aligned with maxima of the representation functions.

### 2.3 Functional approach: objective function

We have defined the field curvature as the sum over each edge of the squared difference between frames at the edges extremities. In our formalization, the difference between two frames (at vertices and ) is no longer the rotation angle, but the norm of the difference between the corresponding functions : . It leads to the new objective function:

 E = ∑ij∫2π0(Fj(α)−Fi(α))2dα = ∑ij∫2π0(Baj−Bai)2dα = ∑ij(aj−ai)⊤(∫2π0B⊤Bdα)(aj−ai)

As the function basis is orthogonal, and all functions are of norm , so the expression simplifies to:

 E=π∑ij∥aj−ai∥2 (1)

### 2.4 Functional approach: constraints

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

The second constraint is that frames of boundary vertices must have one member vector equals to the normal direction. If denotes the normal direction, the frame can be directly fixed by satisfying two equations:

 ai0 = cos(4θi) (2) ai1 = sin(4θi)

However, as we already have the feasibility constraint , enforcing only one equation has the same effect:

 ai0cos(4θi)+ai1sin(4θi)=1. (3)

### 2.5 Toy example

It is natural to ask the question: “Does minimizing our energy minimizes the field curvature as well?”

Two frames and are represented by and , both located on the unit circle. The field curvature measures the circle arc length between them, whereas our norm is the chord length between and .

From a practical point of view, we want to produce smooth fields, so most edges will have low curvature. In this case the objective function is almost proportional to the field curvature. If, however, two adjacent frames are not similar (e.g. they are close to singularities), then the function is smoother than the field curvature, making it easier for the optimization algorithm to move singularities.

Let us illustrate our intuition on a very simple interpolation example: a chain of four vertices having its extremities locked. The toy problem is therefore to find two frames interpolating the extremity frames. Figure 4: Top row: interpolation problem with two locked extremity frames. Bottom row, left: field curvature plot. Bottom row, right: our objective function E. The plots are made as functions of the rotations θ1,θ2 of interpolated frames. Both functions share same local minima P0 and P1.

If we represent two free frames by their angle and , we can plot and compare the field curvature versus our objective function (Figure 4). The field curvature is not smooth (it is piecewise quadratic) and we can observe that there are two local minima. Our objective function is smooth, and has exactly the same minima on this example. Note that it could also have a smaller number of minima e.g. if the constrained frames are more similar. Figure 5: Two minima for the toy problem shown in Figure 4. P0 turns the frames counterclockwise while P1 turns clockwise. P0 minimizes energy E and has better field curvature.

Figure 5 shows the frames corresponding to minima and : they differ by the sense of rotation. The point minimizes objective function and has better field curvature. In next two sections we expose current state of the art approach to minimization of the objective function.

### 2.6 Implementation

We have to minimize our objective function (eq. (1)) with linear equality constraints on boundary vertices (eq. (2) or eq. (3)) and quadratic equality constraints for each vertex .

We use [Kowalski et al. (2012)]’s algorithm to solve this problem. It finds an initial solution by relaxing the quadratic feasibility constraints on and finding the nearest feasible solution. Then it performs a number of smoothing iterations to ameliorate the solution. In order to respect the feasibility, the quadratic constraints are linearized at each smoothing step.

##### Initialization

Here, we relax the feasibility constraints, so we need to minimize the quadratic function subject to linear boundary constraints. To do it, we simply replace the linear constraints by 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.

More precisely, the new quadratic function is expressed in the form where is a matrix, our variable vector ( and ) and is a vector. The system of equations is constructed line-by-line:

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

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

 λX2i =λcos4θi λX2i+1 =λsin4θi,

where we set in our experiments.

From and , we just need to solve the linear system to obtain a minimizer of . Then from we can obtain an initialization of by projecting corresponding vectors on the set of feasible coefficients:

 ai←(X2i,X2i+1)⊤/∥(X2i,X2i+1)∥.
##### Smoothing iterations

Each smoothing iteration is similar to the initialization problem, except that we add to the objective function a new quadratic penalty term that corresponds to a linear approximation of the feasibility constraint as done in [Kowalski et al. (2012), p. 6]. As before it is expressed by a new set of linear equations when constructing and : for each vertex , we add one equation , where denotes the solution obtained at the previous iteration.

To solve linear systems we use OpenNL library [Lévy (Lévy)]: it automatically constructs from the set of linear equations and then solves it by the conjugate gradient method.

### 2.7 Toy problem revisited

This section explains our optimization approach on the toy problem already presented in § 2.5. As we mentioned before, at the initialization step we relax the constraints of feasibility of . Unfortunately we can not plot the corresponding energy since without the constraints it becomes four-dimensional. Figure 6: First row: initialization stage solution. Second row: corresponding functional interpretation. Third and fourth rows: frames obtained by the projection of initialization ai and after smoothing iterations.

Top row of Figure 6 shows the solution of the initialization stage. Intuitively, we allow the points not to be on the unit circle. Hence and lie on the chord between locked points and . Second row of Figure 6 shows the corresponding functions and the third row gives the frames obtained by projecting points on the circle of constraints.

Note that the initialization stage produces the correct sense of rotation (Figure 5). However it does not directly produce the optimum point . The reason being that the initialization stage produces points and (before normalization) in the way that all three chord segments are equal: . But after projecting the points onto the feasible circle (3rd row of Figure 6) the corresponding arc lengths are not equal. Therefore, we need a few smoothing iterations to reach the optimum (Figure 5—bottom row).

### 2.8 Results Figure 7: Evaluation of 2D frame field optimization algorithms: singular triangles are highlighted in red. Compared algorithms are: (left) steepest descent of the field curvature, initialized with an axis aligned frame field, (middle) smoothing iterations with our objective function E initialized with axis aligned frame field after 102, 103 and 106 iterations (from left to right), (right) initialization step alone (left) and the initialization step followed by 103 smoothing iterations (right).

Figure. 7–left shows the optimization of the field curvature by a gradient descent algorithm, initialized by an axis aligned frame field. We obtain a frame field that bends to match the boundary constraints, but its singularities remain on the boundary. The system is only able to reach the local minima corresponding to the initial field topology. The field curvature 444The value of the field curvature is relevant only for comparing different fields on the same mesh. is .

Figure 7–middle shows the optimization using only the smoothing iterations, initialized by an axis aligned frame field. We observe that smoothing iterations were able to slightly move the singularities away from the border and even to merge two singularities. It results in a better field curvature (equal to , and after , and iterations).

Figure 7–right shows that the initialization step alone finds a solution with a simple topology and a lower field curvature (). Smoothing iterations further decrease the field curvature (to ).

These observations suggest that our initialization step is mandatory, and smoothing iterations further improve the result. However, it is important to notice that each iteration takes approximately the same time as the initialization step.

##### Boundary constraints

When working only with feasible solutions a single equation (equation (3)) is sufficient to enforce each boundary constraint. Using it for the initialization step is wrong: for example, a normal constraint of angle forces but let free. As illustrated in Figure 8, it can produce very bad frames fields. Therefore we must use two separate equations (2).

There exists a similar issue in : Huang et al [Huang et al. (2011)] use a boundary condition that is not sufficient for the initialization step. It leads to a poor initialization for their smoothing algorithm, making it very slow, and getting possibly locked with a bad initial topology. This issue is discussed in details in §3.6. Figure 8: Boundary constraints for global optimization: if we only use one constraint (eq.3), the initialization step finds a constant frame field on a parallelogram (left). It is therefore mandatory to lock the frames to get the proper boundary constraints (right).

## 3 Optimization of 3d frame fields

Our objective is to extend to the optimization algorithm presented in previous section.

In , our framework allows to retrieve the representation vector that was the key to efficient optimization of direction fields. We now extend our framework to , find a generalization of this representation vector, express the boundary alignment condition with respect to this representation, and derive the optimization algorithm.

### 3.1 Problem settings

The problem is to define, inside a shape, a smooth frame field that is aligned with the boundary of the shape. We are working in discrete settings on a tet mesh. The problem to minimize the field curvature is defined as follows:

• The reference frame is the set of unit vectors forming normals of a cube aligned with coordinate axes (Figure 9).

• A frame is the reference frame rotated by a matrix : . Note that multiplying a matrix by a set is a slight abuse of notation, it means that we obtain a new set where each vector is rotated by the given matrix.

• A frame field is the definition of a frame for each vertex of the tet mesh.

• The boundary constraint: The frame of a boundary vertex must have one of its member vectors equal to the normal of the boundary.

• The rotation angle between two frames is the minimal angle of rotation that brings one frame to the other.

• The curvature of a frame field is the sum, over each edge, of the squared rotation angle between adjacent frames.

• A tet is called singular if any of its triangles is singular.555We can define what a singular tet is, but we are not able to characterize them by an equivalent of the index in . This fact is discussed in the supplemental material. The triangle is singular if and only if , where denotes the rotation matrix that brings the frame to the frame .

### 3.2 Frame representation

The reference frame is represented by the function , where is the real valued spherical harmonic of degree and order . 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]). Function is defined as , but we are interested by its restriction to the unit sphere .

Any other frame can be obtained as a rotation of by a matrix . It is represented by the function , where is a point of the unit sphere (Figure 9).

forms a functional basis over the unit sphere with an interesting property: applying a rotation to a spherical harmonic of degree produces another harmonic of degree . As a consequence, since we represent the reference frame by a sum of two spherical harmonics of degree 4, each frame function can be represented in the basis . Using it, we can rewrite the expression for the reference frame function as with . Any other frame can be represented by , where with being a rotation matrix acting on coefficients space. Appendix A.1 describes the construction of the rotation matrices.

A feasible coefficient vector is a vector that can be written as where is a rotation matrix that can be derived from a rotation. Geometrically, is constrained to be on a manifold of dimension embedded in the coefficient space.

At this point we can consider the coefficient vector as an extension of the representative vector used in the direction field literature. It is also the representation introduced in [Huang et al. (2011)]. Figure 9: A frame f is the reference frame ~f rotated by a 3×3 matrix R. The plots of the corresponding functions F and ~F are also rotated by R, and their coefficients vectors verify a=RB~a where RB is a 9×9 rotation matrix given in §A.1.

### 3.3 Objective function

As in , the objective function is the sum, over each edge , of the squared difference between frames located at the edges extremities. In our formalism, the difference between two frames is not the rotation angle, but the norm of the difference between the corresponding functions: . It gives the energy:

 E=∑ij∫2π0(Fj(α)−Fi(α))2dα

Here, the function basis is orthonormal, so the expression simplifies to:

 E=∑ij∥aj−ai∥2 (4)

### 3.4 Constraints

There are two types of constraints: each coefficient vector must be feasible, and boundary frames must have one vector aligned with the normal of the volume boundary. The first constraint is presented in the frame representation section, and will be enforced by our optimizer in a dedicated “projection” step (the counterpart of the normalization of in the case). Here we focus on the boundary constraint.

#### 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.

Case 1: the normal is equal to the axis.

Let us first consider the case where the fixed vector (the surface normal) is the axis . If we rotate around by angle , we obtain with being a rotation around . The simple structure of together with the null coefficients of gives the equation:

 a=(√512sin4θ,0,0,0,√712,0,0,0,√512cos4θ)⊤ (5)

As done in the construction of coefficient vectors in the case, we can get rid of the angle by replacing it by a vector .

 a = √712(0,0,0,0,1,0,0,0,0)⊤ (6) + c0(0,0,0,0,0,0,0,0,1)⊤ (7) + c1(1,0,0,0,0,0,0,0,0)⊤ (8)

With this equation, all frames having a vector equal to can be represented by the vector . As in the case it comes with a norm constraint: .

The variable defines the rotation of the frame around the surface normal i.e. a frame field. The optimization of this frame field using as variables is exactly what we did in by introducing the coefficient vector . Our solution restricted to the object boundary is therefore almost 666The boundary has curvature that was not assumed in our frame fields. equivalent to our solution (Figure 10). Figure 10: A 3D frame field produced on a (2D) disk (Left) produces the 2D frame field we could obtain from the 2D algorithm (Right).

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

To handle this more general case, we rotate the constraint. If we want the vector to be preserved, we first compute a rotation that brings axis to . From this rotation, we compute the corresponding rotation matrix , and derive the constraints:

 a = √712RB(0,0,0,0,1,0,0,0,0)⊤ (9) + c0RB(0,0,0,0,0,0,0,0,1)⊤ (10) + c1RB(1,0,0,0,0,0,0,0,0)⊤ (11)

This expression of the normal constraint gives us a set of linear equations per boundary vertex. It introduces two new variables and , and a quadratic constraint .

Note As in the case, the boundary constraint has a simpler expression [Huang et al. (2011)] : that is valid only if all are feasible. Consequently, it cannot be used safely during the initialization step (see Figure 19).

#### Non smooth vertices

Frames of vertices located on hard edges 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 to the first normal, and to the second normal. We compute the corresponding coefficient space rotation and fix the frame coefficient vector to .

### 3.5 Implementation

We have to minimize our objective function (eq. (4)) with linear equality constraints on boundary vertices eq. (9), quadratic equality constraints on boundary vertices, and the constraint that each is feasible.

As in the case (in §2.6), our minimization algorithm (Algo. 1) is formulated as a series of least squares problems (minimize ), where and are constructed without the feasibility constraint of at the first iteration (initialization), and with a linear approximation of it in the subsequent iterations (smoothing iterations).

• Initialization: Our variable vector must represent the representation vectors but also the variables introduced to express the boundary constraint. To do so, we first reorder vertices to have boundary vertices first 777It is possible to increase the performances by by doing a Hilbert sort in the boundary vertices block, and another one in the free vertices block. We can then organize the variable vector by blocks: , and where is the number of vertices.

As in the case, the matrix and the vector are constructed iteratively by adding new equations for the objective function (algorithm 3) and the boundary constraints (algorithm 2). In our approach, we do not explicitly enforce the feasibility of (), but it will be indirectly respected by the feasibility of .

The projection of on the set of feasible coefficient vectors is no longer a simple normalization. Instead we perform, for each vertex, a gradient descent (algorithm 5) initialized by . More precisely, starting with we rotate our current frame 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 norm induced by infinitesimal rotation matrices with Euler’s angles.

• Smoothing iteration: For the linearized feasibility constraint of , we must also add extra variables per vertex to our system. These variables account for the position in a local basis of the tangent space of the manifold of feasible . The introduction of these constraints in the matrix is detailed in algorithm 4.

This frame field design algorithm can be implemented without being expert in spherical harmonics. We give explicit construction of matrices in the appendix. The system to solve is simply a linear system with a positive definite matrix. We use the OpenNL library [Lévy (Lévy)] because the system can be directly constructed from the equations (lines of and elements of ).

In order to keep the algorithm easy to read, we did not detail how to lock frames for vertices with multiple normal constraints.

### 3.6 Results

It is impossible to compare frame field design algorithms only from the images and results presented in the state of the art papers. First of all, our implementation of [Huang et al. (2011)] has significantly better performances compared to what was presented in the original paper. Next, Li et al. [Li et al. (2012)] did not present any frame field results, but only hex meshes that was the main focus of their work. Therefore, we implemented both methods; there are few points worth noting:

Sampling: In previous works the frame fields were sampled either on each tet face or on each tet. Instead we sample it on vertices, otherwise we would not be able to compare corresponding energies.

Gimbal Lock: Both Huang and Li use Euler angles as variables in their L-BFGS optimization, which have numerical issues when the angles are close to the gimbal lock. Note that each frame can be represented by 48 triplets of equivalent Euler angles. In our implementation we choose the triplet maximizing the distance to the nearest gimbal lock.

Rendering: For rendering purposes, we rely on a combination of techniques (Figure 19) to show the spherical harmonics field, the frame field (locally and globally) and the field topology.

#### 3.6.1 Comparison with Huang’s method

Recall that Huang et al. proposed a method in two steps:

• find an initial frame field by solving a linear system and projecting the solution onto the manifold of feasible solutions

• represent each frame by a triplet of Euler angles and optimize the smoothness using an L-BFGS descent method.

Our implementation produces results very similar to those presented in [Huang et al. (2011)], but with significantly better timings. For example, the rockarm (Figure 11) with one million tetrahedra takes about 10 minutes on a single thread application on a Dell M6600 laptop compared to 155 minutes obtained by Huang et al. on a two-thread i7 processor.

Huang’s initialization is very similar to ours, however their boundary condition is not sufficient in this case (it requires the SH coefficients to be feasible). Moreover, they enforce the boundary condition with a penalty term that is very light ( weight). As a consequence, their initial fields are almost constant everywhere (it maximizes the smoothness), with a topology very far from being optimal. The smoothing iterations are performed with much higher weight () of the penalty term using L-BFGS.

After the initialization step we measured the deviation of the field from the given constraints on the rockarm model. Note that the penalty term being the sum of deviations over all vertices, we can conclude that deviation at a given vertex belongs to . Thus on the rockarm the initial frame field has the average deviation of , whereas the maximum frame deviation is .

If we use a much higher penalty weight to enforce the boundary constraint, we obtain an initial frame field with average deviation from constraints equal to and maximum deviation equal to . The field has better topology and L-BFGS converges faster for this initialization. The initialization provided by our method has average deviation from constraints equal to with maximum deviation of .

Figure 11 gives an illustration, it compares three methods: Huang’s algorithm (left image) Huang’s algorithm with much higher penalty weight (middle image) and our initialization followed by Huang’s smoothing iterations (right).

Huang’s paper is the pioneer work and the main contribution in [Huang et al. (2011)] was the introduction of the energy used in frame field optimization, however the initialization is not very good and smoothing stage was also outperformed by later works (see §3.6.3). Figure 11: Comparison of initialization of Huang’s algorithm with their weight of boundary penalty term (10−2) (left/blue), with a much higher weight (102) (middle/yellow), and our initialization (right/red). All smoothing iterations are performed by Huang’s algorithm. For the rocker model, the energies (we take the run with our initialization for 100%) are respectively 94,7%, 101,6% and (obviously) 100%. For the tet model, we obtain 91%, 89,7% and 100%.

#### 3.6.2 Comparison with Li’s method

Li’s initialization computes a frame field on the surface, then propagates it inside the volume by advancing front. As a consequence, the initial field perfectly matches boundary constraints, but is discontinuous across its medial axis.

The smoothing iterations are performed by L-BFGS. It acts on a new set of variables that characterizes, per vertex, the rotation that brings the reference frame to the current frame. For the frames located inside the object, variables are the Euler angles as in Huang’s method. For frames located on the object boundary, the rotation is characterized by a single rotation angle around the normal vector.

The frame field results presented in their paper were limited to an ellipsoid and a sphere (frame field design was not the primary objective). We guess that most of presented results were not fully automatically generated frame fields, as they wrote: “For instance, we use guiding boxes to modify the frames inside the narrow ears of Bunny (Figure 13-a) and the head of rock arm (Figure 13-b) to reduce singularities”.

Moreover, before implementing the method, we thought that their algorithm was strongly limited by the original field topology from the sentence: “However, our propagation-based frame field initialization likely generates singular edges around the medial axis of the volume, and most of them cannot be eliminated by frame optimization.” Surprisingly, our implementation of their method is able to generate smooth frame fields automatically in most situations, even when the topology of the initial field is complex close to the medial axis. Our implementation is slightly different:

• we initialize the field by our algorithm restricted to boundary vertices

• we sample the field on vertices

• and we prevent gimbal locks by a proper initialization of Euler angles.

The only real failure case we found using their method is due to the front propagation algorithm: when a boundary frame is copied to a large number on inner samples. In this case, the L-BFGS solver is sometimes locked with a bad topology (see Fig. 12).

On more complex examples, we have compared their algorithm against our initialization followed by their smoothing iterations. In most cases, we obtain an energy that is a bit better (Fig 13). We also observed that their topology often differs from ours (Fig 14), so we conjecture that our topology is somehow better. However, the quality of the field topology depends on the application, and is not well evaluated by the energy, even for topologies very far from being optimal (see e.g. Fig 11). Figure 12: Li’s algorithm (red) is compared to our algorithm (green) on a one-finger bowling ball. The hole has a huge impact on the initialization due to the advancing front approach (second column, inside the yellow box). As a result, their initialization provides a field with a poor topology and smoothing iterations are not sufficient to find the expected topology (like ours). Our final energy is 86% of theirs. Figure 13: Comparison of two fields generated by Li et al. smoothing iterations: using their initialization (red) or ours (blue). Our energy in percent of theirs is 99%, 99.87%, 100.5%, 99.97%, 99%, 99.6%. The difference is always very low (<1%) but always in our advantage except for the third model. Figure 14: Comparison of two fields generated by Li’s smoothing iterations: using their initialization (red) or ours (blue). Our energy in percent of theirs is 98.5%. Close-ups show the field where the singularity graphs diverge: one is inside the volume (leftmost) and others are on the object boundary. Our results are on the top row and theirs are on the bottom row with singularity encircled in white.

#### 3.6.3 Comparison of smoothing iterations

In the previous sections we have shown (Fig. 11 and 12) that our method provides the best initialization, however our smoothing iterations are not clearly better than others.

Figure 15 shows a comparison of three different smoothing strategies (our linearization, L-BFGS proposed by Huang et al. and L-BFGS proposed by Li et al.). In all three cases we use our method to initialize the field. L-BFGS smoothing proposed by Huang et al. is the slowest in all test cases. First of all, in our implementation the time to evaluate the energy and the gradient is four times slower for the method by Huang et al. than for the method by Li et al. Second, the cruicial difference between these two methods is the way to enforce the boundary alignment: Huang et al. use a penalty term, whereas Li et al. use the set of variables directly satisfying the boundary constraints. In our test we noticed that usage of penalty terms increases the number of iterations to converge and interferes with topological choices to be made, leading to inferior final fields.

We also noticed that the behaviour of linearization changes in function of how far the initialization is from the final solution.

• First row of figure 15 shows a simple case without topology changes, the smoothing iterations change the field geometry only. In this case the linearization of feasibility constraints works flawlessly, in two iterations the method converges, taking about the same time as the smoothing by Li et al.

• Middle row shows a second case, where few topology changes must be made. It slowes the linearization down, even if two iterations produce a reasonably good field.

• Finally, the bottom row shows the case where the initial topology is really bad. The linearization method fails on this model: two first iterations are still very far from the final solution and to reach the minimum it requires four times more time than the method by Li et al.

We conclude that the best solution is our initialization followed by the optimization of Li et al.. In practice, the implementation of our linearization smoothing iterations is almost free (incremental with respect to our initialization algorithm), whereas Li et al. smoothing algorithm is more difficult to implement. Moreover, in most cases few iterations of linearization steps suffice to obtain a fairly good field. As a consequence, we suggest starting with our smoothing iterations (almost free to implement), then possibly replace it with Li et al. smoothing algorithm if performances are not sufficient. Figure 15: Comparison of smoothing iteration algorithms combined with our initialization algorithm. We compare Huang’s method (blue), Li’s method (green), our method (orange), and our method limited to two iterations (red). Top row compares the initialization (left) with the field after two iterations of our algorithm (right). Middle and bottom rows compare the field with two iterations of our smoothing algorithm with other smoothing strategies. Singularity graphs reflect nicely the convergence of thee algorithms. We obtain energy (we take Li’s result for 100%) of resp. 99.98%, 100%, 99.98% and 100.5% on the sector, 101.5%, 100%, 100.4%, and 106.5% on the fandisk, and 116%, 100%, 109%, 185% on the one-finger bowling ball. Figure 16: Our results are shown using combinations of the following rendering techniques. We can plot for each vertex i its Fi (upper left), or its frame as a cube(upper right). We can show the singular tets (lower left), or a smoothed and refined version (lower middle) to better see it in 3D (thanks to the lighting). The field inside the volume can be rendered by curved french fries.

## Future works Figure 20: The thin object (left) frame field is suitable to produce a hex mesh: its singularity graph is basically two singularities of index 1/4 and −1/4 extruded in the z direction. The frame field of fat object (right) has a singularity graph that does not correspond to a hex mesh.

From an application point of view, our method makes it fast and easy to produce smooth frame fields. The smoothness is not necessary the optimal objective for applications such as hex remeshing (see Fig. 20), but can be used as a regulation term for more complex energies. It is also very easy to modify our method to add constraints inside a volume. Figure 19 shows a frame field constrained by faults in geological data. Note that the field is not constrained by the boundary of the model. Such a frame field is useful to steer element placement for fluid simulation used in oil exploration.

From a theoritical point of view, it would be interesting to better understand the shape of the feasible set of (the manifold embedded in ). It could help to find a better projection algorithm than our current gradient descent.

## Conclusion

This work unifies the frame field design problem in and . Both problems are formulated with a similar representation of frames, constraints and objective function. As a consequence, they can also be solved by similar algorithms.

From this analysis, we discovered that the best actual solution to produce smooth frame fields is to initialize it by our proposed extension of [Kowalski et al. (2012)], followed by smoothing iterations of [Li et al. (2012)]. The main drawback of this solution is requirement to implement two very different approaches (a sparse linear system solver and a L-BFGS descent). A fair alternative is to use our extension of [Kowalski et al. (2012)], it is simple to implement and requires a linear system solver only. In practice for our models we perform only two or three linearization iterations, however if the initialization is a bad guess (e.g. the sphere), it can be insufficient. With this approach we are able to generate (on a laptop) fields on the models up to few millions tetrahedra in less than 10 minutes, refer to Figure 19 for an illustration.

## Appendix A SH cookbook

### a.1 9d rotation

Let us denote by , and matrices of rotation around axes , and respectively. Any frame can be obtained as a composed rotation of the reference frame , where the reference frame is the set of 6 unit vectors aligned with coordinate axes:

 f=Rx(α)×Ry(β)×Rz(γ)×~f.

If are Euler angles of rotation between a frame and , the representation vector is calculated as , where , and are matrices of rotation defined as follows:

 RzB(γ)=⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣cos(4γ)0000000sin(4γ)0cos(3γ)00000sin(3γ)000cos(2γ)000sin(2γ)00000cos(γ)0sin(γ)000000010000000−sin(γ)0cos(γ)00000−sin(2γ)000cos(2γ)000−sin(3γ)00000cos(3γ)0−sin(4γ)0000000cos(4γ)⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦
 RxB(π/2)=⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣00000√14/40−√2/400−3/40√7/40000000000√2/40√14/400√7/403/40000000003/80√5/40√35/8−√14/40−√2/40000000000√5/401/20−√7/4√2/40−√14/40000000000√35/80−√7/401/8⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦
 RyB(β)=RxB(π/2)×RzB(β)×RxB(π/2)⊤
 RxB(α)=RyB(π/2)⊤×RzB(α)×RyB(π/2)

These matrices are called Wigner D-matrices and the literature on their construction is vast [Collado et al. (1989), Blanco et al. (1997), Choi et al. (1999), Ivanic and Ruedenberg (1996)]. However, as we are using degree 4 harmonics only, we performed the symbolic computation as follows.

The matrix is easy to compute, given a spherical harmonic we can rotate it around the -axis by the angle by evaluating . Then the element of the matrix is simply .

Matrices and are trickier to get, however we can use the fact that arbitrary rotation around the -axis can be decomposed into the rotation around -axis by followed by a rotation around -axis. To rotate a spherical harmonics around the -axis by we can perform the following substitution:

 (θ,ϕ)←(arccos(−sin(θ)sin(ϕ)),atan2(cos(θ),sin(θ)cos(ϕ))).

Then the matrix is calculated by evaluating double integrals of all possible products between basis functions and rotated functions.

### a.2 Linearization

In order to linearize the constraints in our linear solver we define matrices and as follows:

 ExB=⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣0000000−√20000000−√7/20−√200000−3/√20−√7/200000−√100−3/√200000√1000000003/√20000000√7/203/√200000√20√7/20000000√20000000⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦
 EyB=⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣0√20000000−√20√7/20000000−√7/203/√20000000−3/√200000000000−√100000000√100−3/√200000003/√20−√7/20000000√7/20−√20000000√20⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦
 EzB=⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣000000004000000030000000200000001000000000000000−10000000−20000000−30000000−400000000⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦

It is easy to verify that these matrices are chosen to verify the following equations for small rotations :

 RxB(α) =I9×9+αExB+o(|α|) RyB(β) =I9×9+βEyB+o(|β|) RzB(γ) =I9×9+γEzB+o(|γ|).

Finally, for small rotations the multiplication is commutative:

 RB(α,β,