Algebraic Representations for Volumetric Frame Fields

TL;DR

This paper introduces a mathematically rigorous framework for representing and optimizing 3D volumetric frame fields, enabling more accurate and stable meshing in three-dimensional modeling.

Contribution

It develops a differential and algebraic geometric understanding of octahedral frame spaces and proposes geometry-aware optimization tools, including semidefinite relaxations, for volumetric frame fields.

Findings

Efficient algorithms produce high-quality, stable 3D frame fields.

Mathematically sound description of octahedral and odeco frames.

Generalization to anisotropic frames with independent axis scaling.

Abstract

Field-guided parametrization methods have proven effective for quad meshing of surfaces; these methods compute smooth cross fields to guide the meshing process and then integrate the fields to construct a discrete mesh. A key challenge in extending these methods to three dimensions, however, is representation of field values. Whereas cross fields can be represented by tangent vector fields that form a linear space, the 3D analog---an octahedral frame field---takes values in a nonlinear manifold. In this work, we describe the space of octahedral frames in the language of differential and algebraic geometry. With this understanding, we develop geometry-aware tools for optimization of octahedral fields, namely geodesic stepping and exact projection via semidefinite relaxation. Our algebraic approach not only provides an elegant and mathematically-sound description of the space of octahedral frames but also suggests a generalization to frames whose three axes scale independently, better capturing the singular behavior we expect to see in volumetric frame fields. These new odeco frames, so-called as they are represented by orthogonally decomposable tensors, also admit a semidefinite program--based projection operator. Our description of the spaces of octahedral and odeco frames suggests computing frame fields via manifold-based optimization algorithms; we show that these algorithms efficiently produce high-quality fields while maintaining stability and smoothness.

Equations

$$H=\frac{1}{|\mathrm{O}|}\sum _{o\in \mathrm{O}}\rho (o).$$

$$\rho (o^{\prime })H=\frac{1}{|\mathrm{O}|}\sum _{o\in \mathrm{O}}\rho (o^{\prime })\rho (o)=\frac{1}{|\mathrm{O}|}\sum _{o\in \mathrm{O}}\rho (o^{\prime }o)=\frac{1}{|\mathrm{O}|}\sum _{o\in \mathrm{O}}\rho (o)=H.\text{qed}$$

$$q_0={\left(0,0,0,0,\sqrt{\frac{7}{12}},0,0,0,\sqrt{\frac{5}{12}}\right)}^{\top }\in {\mathbb{R}}^9.$$

$$F=\rho (\mathrm{SO}(3))q_0=\{\rho (r)q_0\mid r\in \mathrm{SO}(3)\}.$$

$${\pi }_{\alpha }\circ \rho :\mathrm{SO}(3)\to F_{\alpha }$$

$$(D\rho {)}_I:\mathfrak{so}(3)=T_I\mathrm{SO}(3)\to \mathfrak{so}(9)\subset {\mathbb{R}}^{9\times 9}.$$

$$\{(D{R}_g{)}_I{l}_i{\}}_{i=1}^3.$$

$$D({\pi }_{\alpha }\circ \rho {)}_g(D{R}_g{)}_I{l}_i$$

$$=(D{\pi }_{\alpha })(D{R}_{\rho (g)}{)}_I(D\rho {)}_I{l}_i$$

$$=(D{\pi }_{\alpha })L_i\rho (g)$$

$$=L_i\rho (g)q_{\alpha },$$

$$⟨L_{i}q,L_{j}q⟩={\delta }_{ij}\forall q\in F_{\alpha }.$$

$$⟨L_i{q}_{\alpha },L_j{q}_{\alpha }⟩={\delta }_{ij}=⟨l_i,l_j⟩.$$

$$⟨L_i\rho (g)q_{\alpha },L_j\rho (g)q_{\alpha }⟩$$

$$=⟨D{\rho }_I(g^{\top }{l}_{i}g)q_{\alpha },D{\rho }_I(g^{\top }{l}_{j}g)q_{\alpha }⟩.$$

$$⟨L_i\rho (g)q_{\alpha },L_j\rho (g)q_{\alpha }⟩$$

$$=\sum _{r,s}{a}_{ir}{a}_{js}⟨L_r{q}_{\alpha },L_s{q}_{\alpha }⟩$$

$$=\sum _k{a}_{ik}{a}_{jk}={\delta }_{ij}$$

$$T=\sum _{i=1}^n{\lambda }_i(v_i{)}^{\otimes d}$$

$$T\in {\mathrm{Sym}}^d{\mathbb{R}}^n\mapsto p_T(x)=\frac{1}{d!}T(x,\dots ,x)\in \mathbb{R}[x].$$

$$p_T(x)=\sum _i{\lambda }_i(v_i^{\top }x{)}^d$$

$$p(x)=\sum _{d_1+\cdots +d_n=d}{u}_{d_1,\dots ,d_n}{x}_1^{d_1}\dots x_n^{d_n}$$

$$u^{\top }{A}_{i}u=0$$

$$q=(q_0,q_2,q_4)\in V_0\times V_2\times V_4$$

$$\Vert q_2{\Vert }^2=C_2\left(\sum _i{\lambda }_i^2-\sum _{i<j}{\lambda }_i{\lambda }_j\right).$$

$$\{q\in \tilde{F}\mid q_2=0\}$$

$${\left(\begin{array}{c}1\\ q_4\end{array}\right)}^{\top }{P}_i\left(\begin{array}{c}1\\ q_4\end{array}\right)=0,i=1,\dots ,15$$

$${\exp }_q(\mathfrak{v})=\rho (r^{\top }\exp (\Vert v\Vert l_3)r)q=\rho (r{)}^{\top }\exp (\Vert v\Vert L_3)\rho (r)q,$$

$$v=\Vert v\Vert (\cos \theta \sin \phi ,\sin \theta \sin \phi ,\cos \phi ).$$

$$r=\exp (-\varphi l_2)\exp (-\theta l_3)=r_{23}^{\top }\exp (-\varphi l_3)r_{23}\exp (-\theta l_3),$$

$$\rho (r)=R_{23}^{\top }\exp (-\varphi L_3)R_{23}\exp (-\theta L_3),$$

$${\Pi }_F(y)=\underset{q\in F}{\mathrm{argmin}}\Vert q-y{\Vert }^2.$$

$$\underset{Q\in {\mathbb{R}}^{10\times 10}}{\mathrm{argmin}}$$

Full text

Hexahedral meshing, octahedral frame fields, convex relaxations, convex algebraic geometry

††submissionid: TOG-19-0069.R1††copyright: rightsretained††journal: TOG††journalyear: 2020††journalvolume: 39††journalnumber: 2††article: 16††publicationmonth: 3††doi: 10.1145/3366786††ccs: Computing methodologies Computer graphics††ccs: Computing methodologies Mesh geometry models††ccs: Computing methodologies Volumetric models

1. Introduction

Inspired by the success of field–based approaches to quadrilateral meshing on surfaces (cf. [Vaxman et al., 2016]), recent research in applied geometry has focused on developing an analogous approach to hexahedral meshing. Motivated by applications in finite-element modeling, hexahedral meshing is the problem of dividing a given volume into hexahedral elements (deformed cubes) with minimal distortion and such that mesh boundary faces are aligned to the boundary of the volume.

Hexahedral meshing couples a geometry problem—minimizing distortion of mesh elements—to a combinatorial problem—placing mesh elements to achieve a desired connectivity structure. As in 2D, field-based approaches first ignore combinatorial constraints and solve for a frame field, which represents the local alignment and singular structure of a mesh (see Figure 1). Then, they integrate that field to guide the placement of hex elements [Nieser et al., 2011]. So as not to impose unnatural constraints, the space of frame fields must be expressive enough to represent the range of possible singularities that may appear in hexahedral meshes. These singularities are described by gluing relations restricted to the symmetries of a cube, i.e., the octahedral group.

One might hope that 2D field-based methods would extend easily to 3D. There are at least two obstructions to transferring ideas from the 2D case. First, the singular structure of a 3D frame field can be much more complicated than that of a cross field, comprising an embedded graph rather than a set of isolated points. Second, we must understand the geometry of the space of frames to measure and optimize smoothness of a frame field.

Complementing the recent results of Liu et al. [2018] on frame field singularities and necessary conditions for hex-meshability, in the present work we study the second challenge, namely characterizing the geometry of the space of frames. A 3D frame field is, intuitively, an assignment of three mutually orthogonal directions to each point in a volume. An orthonormal basis of three vectors—comprising an orthogonal matrix—is sufficient to specify a frame, but thanks to octahedral symmetry, multiple orthogonal matrices can specify the same frame. Formally, the space of octahedral frames can be viewed as the quotient of the group of 3D rotations $\mathrm{SO}(3)$ by the right action of the octahedral group $\mathrm{O}$ comprising the rotational symmetries of a cube (see §3). The non-commutativity of 3D rotations makes the geometry of this quotient more complicated than that of its 2D counterpart. Recent attempts to lift ideas and techniques from 2D have either ignored the geometry of the frame space entirely or treated it as a black box for nonlinear optimization.

Consider perhaps the simplest form of optimization on a space, projection—finding the closest point in the space to a given point in an ambient space. Previous work on frame fields has treated this projection problem as a nonlinear, nonconvex optimization problem over frames parametrized by Euler angles, with no guarantees on convergence or global optimality. Our description of the octahedral frame space as an algebraic variety suggests a different approach to projection based on semidefinite programming, which yields a certificate of global optimality in polynomial time. Our semidefinite relaxation of projection is exact in a neighborhood of the octahedral variety, and we conjecture—with strong empirical evidence—that it is so universally.

Even when conducting local optimization on the space of octahedral frames, parametrization by Euler angles may not be the best approach. We show that the map from $\mathrm{SO}(3)$ to the octahedral variety is a local isometry, enabling us to compute geodesics on the octahedral variety in closed form. Manifold-based optimization that moves along geodesics can then be used to accelerate local optimization dramatically.

Beyond precisely characterizing the space of octahedral frames, our algebraic approach admits a generalization to frames whose axes scale independently. This larger space better captures frame field geometry, e.g. allowing for a nonzero direction aligned to singular arcs even if the directions orthogonal to the arcs must vanish. We call these new objects odeco frames, thanks to their construction using orthogonally decomposable tensors, and we derive relevant projection operators.

Our experiments show how the theoretical objects we study enable volumetric frame field design in practice. In particular, we apply standard manifold-based optimization algorithms to field design, built on our differential and algebraic descriptions of octahedral and odeco frames (see Figure 2). The end result is an efficient suite of techniques for producing smooth fields that obey typical constraints for our target application.

Outline

In §3.1 we reintroduce the spherical harmonic representation of octahedral frames and demonstrate how this amounts to an equivariant isometric embedding of the quotient $\mathrm{SO}(3)/\mathrm{O}$ into $\mathbb{R}^{9}$. In §3.2 we introduce the odeco frames, whose axes can scale independently, and we exhibit the spaces of octahedral and odeco frames as nested varieties cut out by quadratic equations. In §4 we describe essential primitives for optimization over octahedral and odeco frames, namely geodesics and projection via semidefinite programming. In §5 we formulate the frame field optimization problem. In §6 we describe two algorithms for optimizing fields. §7 describes experiments using these algorithms. A discussion and conclusion follow in §8. Further experimental results are included in the supplemental material. In summary, our contributions are

• •

a proof of isometric embedding of $\mathrm{SO}(3)/O$ in $\mathbb{R}^{9}$;

• •

descriptions of the spaces of octahedral and more general odeco frames as nested algebraic varieties; and

• •

new state-of-the-art optimization techniques for volumetric frame fields valued in both varieties, featuring geodesics and SDP-based projection as primitives.

2. Related Work

2.1. 2D Frame Fields and Quadrilateral Meshing

Cross fields, and their application to quad meshing, have been studied extensively in geometry processing (cf. [Vaxman et al., 2016]). A useful insight from cross field research is that it is advantageous to replace field values defined up to some symmetry with a representation vector—that is, some function of the field value invariant under the relevant symmetry. In two dimensions, four vectors forming a right-angled cross can be represented in a unified manner by their common complex fourth power. This amounts to an embedding of a quotient manifold into a Euclidean space; we show that the octahedral variety generalizes this idea to 3D in a natural and isometric manner.

Recently, an effort has been made to unify and formalize the various algorithmic approaches to cross fields, borrowing from the Ginzburg–Landau theory from physics [Beaufort et al., 2017; Viertel and Osting, 2019; Osting and Wang, 2017]. This amounts to replacing an ill-posed unit norm constraint with a penalty term and taking the limit as the penalty parameter goes to infinity. Viertel and Osting [2019] show that local optima of such a procedure have isolated singularities with indices $\pm 1/4$, as appropriate for quad meshing. They propose a diffusion-generated algorithm to compute such local optima. Inspired by the work of Merriman, Bence, and Osher (MBO) [1992] on mean curvature flow, this algorithm alternates between finite-time diffusion and pointwise normalization. Osting and Wang [2017] study an analogous method applied to orthogonal matrix–valued fields. Our projection operators enable us to develop similar MBO diffusion-generated methods for optimization of octahedral and odeco fields. In §7 we demonstrate that a modified strategy, where the diffusion parameter is adjusted on-the-fly, frequently helps to avoid local minima.

2.2. 3D Frame Fields and Hexahedral Meshing

Huang et al. [2011] introduce a representation of frames as functions over the sphere that exhibit octahedral symmetry, parametrized by coefficients in the spherical harmonic basis. As an initialization step, they solve a Laplace equation, resulting in coefficients in the interior that do not correspond to octahedral frames. These must be projected via nonconvex optimization over frames parametrized by Euler angles. They further optimize the results by minimizing a discrete Dirichlet energy over Euler angles. Ray et al. [2016] refine the three stages of this approach, with improved boundary constraints in the initialization, an efficient projection technique, and an L-BFGS optimization algorithm.

Solomon et al. [2017] reformulate the initialization step of Ray et al. [2016] using the boundary element method (BEM). This provides a way to harmonically interpolate the boundary conditions exactly, ignoring the constraints in the interior. Finally, sampled interior values are approximately projected onto the constraints as in the previous methods. As the constraints are ignored at the Dirichlet energy–minimization stage, there is no sense in which the final frame fields are optimal.

A related but distinct problem is the computation of fields of symmetric second order tensors (i.e., symmetric matrices) [Palacios et al., 2017]. Every symmetric matrix has an orthonormal basis of eigenvectors, which corresponds to an octahedral frame. One might thus think symmetric matrix fields can be used to parametrize octahedral frame fields. While a symmetric matrix field corresponds to at least one frame field, the correspondence is not one-to-one—for example, a field of identity matrices corresponds to all frame fields. Moreover, symmetric matrix fields can only represent singularities whose indices are multiples of $\pm 1/2$, while indices $\pm 1/4$ are crucial for hex meshing [Liu et al., 2018]. Our fourth order octahedral and odeco fields are rich enough to represent all indices.

The work of Fang et al. [2016] on generalized polycubes imposes an even-more-extreme restriction that singularities may only appear on the boundary. After cutting handles, the regular frame field in the interior of the volume may be represented by a field of rotation (orthogonal) matrices. This is further relaxed to a field of matrices, and orthogonality is approximately enforced via a penalty term. The resulting field is used to construct a polycube map of the cut volume, through which a hex mesh is pulled back. By requiring regularity in the interior, this method may exclude fields that achieve lower distortion overall.

The main driving force for research on 3D frame fields has been hexahedral mesh generation. For a broader overview we refer the reader to the surveys of [Armstrong et al., 2015] and [Yu et al., 2015], and we limit the subsequent discussion to techniques involving frame fields. The idea of such methods is to construct a volumetric integer-grid map [Liu et al., 2018], through which portions of the Cartesian integer grid are pulled back to a structure-aligned hexahedral mesh in the input domain. Nieser et al. [2011] introduced a parametrization technique that turns a given frame field into an integer-grid map by solving a mixed-integer Poisson problem. Extraction of the hexahedral mesh from the map is hampered by degeneracies in the map, motivating the sanitization technique of Lyon et al. [2016] to improve robustness. Several refinements of the above ideas have been proposed, including different guidance for the frame field [Li et al., 2012] and heuristics to improve the singularity structure by decimation or splitting [Li et al., 2012; Jiang et al., 2014] to avoid degeneracies of the map.

While robust hexahedral meshing based on frame fields remains an open problem, recently several hex-dominant meshing algorithms have been proposed [Sokolov et al., 2016; Gao et al., 2017] that also use frame fields but circumvent the problem of non-meshable singularities. Gao et al. [2017] propose a hierarchical optimizer for frame fields that is based on local relaxation.

However, many practical applications demand pure hexahedral meshes—e.g. for the construction of volumetric spline spaces in isogeometric analysis—and consequently a full understanding of meshable field topology is required. To this end, Liu et al. [2018] enumerate the singular vertex types that may occur in a hex mesh with bounded edge valence; they also develop a topological index formula analogous to the Poincaré-Hopf formula for vector fields on surfaces. These local and global constraints being established, they propose an algorithm to generate a (meshable) frame field from a prescribed (meshable) singular structure. The theoretical portion of our work complements the topological work of Liu et al. with a closer study of the geometry of spaces of frames.

In concurrent work, Corman and Crane [2019] propose an alternative approach to computing a frame field with prescribed singular structure via a discretized connection on a frame bundle. This approach does not immediately extend to a method for de novo frame field design. However, the connection associated to a field is a natural object for the study of such properties as integrability. We leave to future work the study of connections associated to fields valued in the octahedral and odeco varieties and the extraction of such connections directly from the field coefficients.

2.3. Alternative Frame Representations

Complementing the spherical harmonic–based representation in the papers above on octahedral frame fields, Chemin et al. [2018] propose an equivalent representation of octahedral frames as certain symmetric tensors of order four. They introduce algebraic equations characterizing the octahedral frames among symmetric fourth-order tensors. These equations are equivalent under a change of basis to our defining equations for the octahedral variety. However, by using a basis suggested by the structure of $\mathrm{SO}(3)$, we are able to not only present the defining equations of the octahedral variety, but also illuminate how it is an isometric embedding of the quotient space $\mathrm{SO}(3)/\mathrm{O}$, enabling us to compute geodesics. Additionally, we use our algebraic description of the octahedral variety to introduce a semidefinite relaxation of projection and to place it in the context of the more general odeco variety.

In concurrent work, Golovaty et al. [2019] also represent frames by fourth-order symmetric tensors subject to some algebraic conditions, and they propose gradient flow on a Ginzburg-Landau–type energy to optimize for smooth fields.

An additional alternative representation is proposed by Gao et al. [2017], who use quaternions to encode octahedral frames. Their representation uses relatively few values, but a matching procedure has to be embedded in their optimization objective to account for the fact that $48$ quaternions correspond to the same frame. The concurrent work [Beaufort et al., 2019] uses quaternions to derive another parametrization of octahedral frames by points of a variety in three complex dimensions.

All previous work on 3D frame fields has only considered octahedral frames, which do not capture the unidirectional behavior of frame fields near singular curves (cf. Figure 3). In the algebraic geometry community, Robeva [2016] and Boralevi et al. [2017] have completely characterized a family of algebraic varieties of orthogonally decomposable (odeco) tensors. We show how the octahedral variety is embedded in one of the odeco varieties, and we introduce a technique for optimization over the odeco frames. For our purposes, odeco frames generalize octahedral frames by allowing independent scaling of the “axes” of a frame, including degeneration to unidirectionality at singular curves and to zero at singular nodes.

Finally, Shen et al. [2016] consider fields having different local discrete symmetry groups, such as the tetrahedral and icosahedral groups. They extend the methods of [Huang et al., 2011; Ray et al., 2016] to compute such fields.

2.4. Semidefinite Relaxations

Relaxation of algebraic optimization problems to semidefinite programs has been studied extensively in the field of real algebraic geometry and optimization. Blekherman et al. [2012] provide an introduction to this discipline. The efficacy of semidefinite relaxation in computer science was demonstrated dramatically in the seminal paper of Goemans and Williamson [1995] on the maximum cut problem. Since then, semidefinite relaxation has continued to be employed to solve both combinatorial and continuous optimization problems, such as angular synchronization [Singer, 2011]. In geometry processing and graphics, this machinery has been applied to such problems as correspondence [Kezurer et al., 2015], consistent mapping [Huang and Guibas, 2013], registration [Maron et al., 2016], camera motion estimation/calibration [Agrawal and Davis, 2003; Ozyesil et al., 2015], and deformation [Kovalsky et al., 2014].

Most relevant to the present work are relaxations of the Euclidean projection problem onto a variety defined by quadratic equations, an example of a quadratically-constrained quadratic program (QCQP). Cifuentes et al. [2017] have recently shown a stability result implying that the semidefinite relaxation of Euclidean projection onto a smooth, quadratically-defined variety is exact in a neighborhood of the variety. Cifuentes et al. [2018] have additionally shown that the region in which the relaxation is exact is a semialgebraic set, and they have provided a formula for the degree of its algebraic boundary. Unfortunately, computing this boundary is generally intractible for interesting varieties. A deeper theoretical understanding of when semidefinite relaxations of Euclidean projection are globally exact is still lacking.

3. Spaces of Frames

As previously studied in [Huang et al., 2011; Ray et al., 2016; Solomon et al., 2017; Chemin et al., 2018], the basic unknown in the volumetric frame field problem is a tuple of three mutually orthogonal directions at each point in a volume. These directions may be represented by an orthonormal basis of vectors, but the signs and order of the vectors are irrelevant thanks to octahedral symmetry. This redundancy makes detecting smoothness of a field of tuples difficult. Hence, it pays to use a unified representation invariant to the symmetries of the frame.

Below, we apply machinery from differential and algebraic geometry to derive a succinct description of this basic octahedral frame and show how it is related to rotations of the function $x_{1}^{4}+x_{2}^{4}+x_{3}^{4}$ on the unit sphere expressed in the spherical harmonic basis (as used to represent frames in [Huang et al., 2011; Ray et al., 2016; Solomon et al., 2017]) and to tensorial representations (as used in [Chemin et al., 2018]). Algebraic language not only provides a succinct description of previous representations but also suggests a means of generalizing to frames whose three axes scale independently (e.g. rotations of the function $\sum_{i}\lambda_{i}x_{i}^{4}$ for varying $\lambda\in\mathbb{R}^{3}$), whereas previous work requires them to have the same length ($\lambda_{1}=\lambda_{2}=\lambda_{3}$). This broader set better aligns with the realities of the frame field problem, since singular edges have nontrivial directionality that cannot be captured by existing representations. Relevant background material may be found in Appendix A.

3.1. Octahedral Variety

Intuitively, an octahedral frame field is a smooth assignment of (unoriented) orthonormal coordinate axes to each point $x$ in a region $\Omega\subset\mathbb{R}^{3}$. Such frames are called octahedral because they exhibit symmetries described by the octahedral group $\mathrm{O}$.111Technically, their stabilizers are conjugate to $\mathrm{O}$. The space of such octahedral frames can be identified with the quotient space $\mathcal{F}~{}=~{}\mathrm{SO}(3)/\mathrm{O}$ — that is, the quotient of the group of oriented rotations by the right action of the octahedral group. Since $\mathrm{O}$ is not a normal subgroup of $\mathrm{SO}(3)$ (indeed $\mathrm{SO}(3)$ is simple) $\mathcal{F}$ is not a group. However, as $\mathrm{O}$ is a finite group acting freely on $\mathrm{SO}(3)$, $\mathcal{F}$ is a manifold, and the surjective map $\mathrm{SO}(3)\to\mathcal{F}$ is a covering map. The universal cover of $\mathcal{F}$ is that of $\mathrm{SO}(3)$, namely $\mathrm{SU}(2)$, and so its fundamental group $\pi_{1}(\mathcal{F})$ is the lift of $\mathrm{O}$ to $\mathrm{SU}(2)$, the binary octahedral group $\mathrm{BO}$. In particular, $\mathrm{BO}$ classifies singular curves of octahedral fields [Mermin, 1979].

$\mathrm{SO}(3)$ admits a bi-invariant Riemannian metric (23), which means that the action of $\mathrm{O}$ by right-multiplication is isometric. Thus the Riemannian metric on $\mathrm{SO}(3)$ descends to $\mathcal{F}$, making it a Riemannian manifold. In particular, $\mathcal{F}$ has geodesics that lift to geodesics of $\mathrm{SO}(3)$. We will employ these geodesics to compute optimization sub-steps on $\mathcal{F}$ (cf. §4.1).

We have introduced $\mathcal{F}$ as an abstract smooth manifold, but an embedding of $\mathcal{F}$ in a Euclidean space is necessary to make it amenable to computation. Previous works have introduced this equivariant embedding via the action of $\mathrm{SO}(3)$ on polynomials and spherical harmonic coefficients. For completeness, we reintroduce it here, but in a slightly different way that highlights the role of representation theory and the orbit-stabilizer theorem. Later, we will show that the embedding of $\mathcal{F}$ in $\mathbb{R}^{9}$ is an equivariant isometric embedding (cf. [Moore, 1976]).

Consider the irreducible representation $\rho:\mathrm{SO}(3)\to\mathrm{SO}(9)$ corresponding to the fourth band of spherical harmonics. The $9\times 9$ orthogonal matrices $\rho(r)$ for $r\in\mathrm{SO}(3)$ are sometimes referred to as Wigner $D$-matrices. Form a linear operator $H\in\mathbb{R}^{9\times 9}$ as

[TABLE]

This $H$ is a projection operator onto the subspace of $\mathbb{R}^{9}$ invariant under all octahedral rotations:

Lemma 3.1.

For any $o^{\prime}\in\mathrm{O}$, $\rho(o^{\prime})H=H$.

Proof.

[TABLE]

Corollary 3.2.

For $q\in\mathbb{R}^{9}$, $\rho(\mathrm{O})\cdot q=q$ if and only if $q\in\mathrm{Im}(H)$.

Proof.

The forward implication follows by writing $q=Hq$. The reverse implication follows directly from Lemma 3.1. ∎

Because $\mathrm{O}$ is a maximal subgroup of $\mathrm{SO}(3)$, we have the following corollary. Here $\operatorname{stab}{q}$ denotes the stabilizer of $q$—that is, the subgroup of $\mathrm{SO}(3)$ leaving $q$ invariant (see Appendix A).

Corollary 3.3.

Let $q\in\mathrm{Im}(H)$ and nonzero. Then $\operatorname{stab}(q)=\mathrm{O}$.

It happens that $H\in\mathbb{R}^{9\times 9}$ has rank one, i.e., its image is one-dimensional, motivating the following definitions.

Definition 3.4.

The canonical octahedral frame is the normalized vector

[TABLE]

such that $H=q_{0}q_{0}^{\mathstrut\scriptscriptstyle{\top}}$.

Our canonical frame is the same as the normalized projection of the polynomial $\sum_{i}x_{i}^{4}$ into the fourth band of spherical harmonics, as in, e.g., [Solomon et al., 2017].

Definition 3.5.

The octahedral variety is the orbit of $q_{0}$ under the action of $\mathrm{SO}(3)$,

[TABLE]

To summarize, $F$ is the orbit of $q_{0}$, whose stabilizer is $\mathrm{O}$. A smooth version of the orbit-stabilizer theorem [Lee, 2012, Theorem 21.18] shows that there is an equivariant diffeomorphism $\phi:\mathcal{F}\to F$:

{\mathrm{SO}(3)}$${\mathrm{SO}(9)}$${\mathcal{F}}$${F}$$\scriptstyle{\rho}$$\scriptstyle{\cdot q_{0}}$$\scriptstyle{\phi}

Next, we will characterize the Riemannian geometry of $F$ by showing that $\phi$ is an isometry up to a uniform scale factor. To our knowledge, this observation has not appeared in previous work.

Proposition 3.6.

Let $\alpha=\sqrt{3/20}$ and $F_{\alpha}\coloneqq\alpha F=\{\alpha q:q\in F\}$. Let $\pi_{\alpha}:~{}\mathbb{R}^{9\times 9}~{}\to~{}\mathbb{R}^{9}$ denote matrix multiplication by the scaled canonical octahedral frame $q_{\alpha}\coloneqq\alpha q_{0}$. That is, $\pi_{\alpha}(A)=Aq_{\alpha}$. Then the map

[TABLE]

is a local isometry, making the induced diffeomorphism $\phi_{\alpha}:\mathcal{F}\to F_{\alpha}$ an isometry.

Proof.

Taking the differential of $\rho$ at the identity yields the associated Lie algebra representation

[TABLE]

$(D\rho)_{I}$ is characterized by $L_{i}\coloneqq(D\rho)_{I}(l_{i})$, the images of the Lie algebra generators $l_{i}$ (see Appendix A for definitions and supplemental material for explicit expressions). $\pi_{\alpha}$ is a linear map, so its differential is also multiplication by $q_{\alpha}$.

To see that $\pi_{\alpha}\circ\rho$ is a local isometry, first recall that the metric on $\mathrm{SO}(3)$ is bi-invariant. In particular, for each $g\in\mathrm{SO}(3)$, an orthonormal basis for $T_{g}\mathrm{SO}(3)$ is given by the right-translated Lie algebra generators

[TABLE]

So it suffices to prove that their images under $\pi_{\alpha}\circ\rho$ form an orthonormal basis in $T_{\rho(g)q_{\alpha}}F_{\alpha}$. But

[TABLE]

where the second equality follows by differentiating the representation property $\rho\circ\operatorname{R}_{g}=\operatorname{R}_{\rho(g)}\circ\rho$ at the identity. So the equation we need to prove is

[TABLE]

This isometry condition can be checked explicitly at $q_{\alpha}$:

[TABLE]

For a general $q=\rho(g)q_{\alpha}\in F_{\alpha}$, we compute

[TABLE]

The first equality follows because $\rho(g)\in\mathrm{SO}(9)$, and the second follows by Lemma A.1. Let $a_{ik}=a_{ik}(g)$ be the coefficients of the adjoint representation of $\mathrm{SO}(3)$ (see Appendix A), i.e., $g^{\mathstrut\scriptscriptstyle{\top}}l_{i}g\eqqcolon\sum_{k}a_{ik}l_{k}$; these coefficients form an orthogonal matrix. Now using linearity of $D\rho$ along with (3) and (4), we obtain

[TABLE]

as required. ∎

In summary, we have described the octahedral variety as an embedded submanifold in $\mathbb{R}^{9}$ isometric to $\mathcal{F}=\mathrm{SO}(3)/\mathrm{O}$. This isometry means that we can do manifold optimization over frames by working in the embedding (cf. §6.1). To show that $F$ is really an algebraic variety, we will need to exhibit equations that cut it out. We will delay this until §3.2.1, when we can give a unified description of the octahedral and odeco varieties.

3.2. Odeco Variety

The previous section provides more insight into the octahedral frames used in all previous work. Octahedral frames are suitable for representing singularity-free frame fields. However, frame fields commonly encountered in applications have singularities comprising an embedded graph. Consider the triangular prism shown in Figure 3. Near the singular curve, a smooth octahedral field would rotate infinitely quickly, and it would not have a well-defined value along the curve. This issue is analogous to the case of unit cross fields—note that for cross fields, the hairy ball theorem requires singularities on simply-connected surfaces.

One solution to this is to replace the hard constraint that the field values lie in the octahedral variety with a penalty term, as motivates the MBO method for octahedral fields detailed below (§6.2). Another solution is homogenization—i.e., allowing field values to scale, replacing singularities by zeros, as is considered for cross fields in [Knöppel et al., 2013]. But consider the triangular prism again: a scaled octahedral field would vanish completely at the singular curve since all three orthogonal axes must scale uniformly. This makes the octahedral representation unable to capture the alignment of the field to the singular curve. To cope with this problem and to show the value of our algebraic approach, we now describe a superset of the octahedral frames. This set allows the axes to scale independently; for instance, as shown in Figure 3, the frame axes orthogonal to the singular curve scale toward zero while the axis tangent to the singular curve remains nonzero.

The symmetric orthogonally decomposable (odeco) tensor varieties, introduced in [Robeva, 2016], parametrize symmetric tensors $T\in\operatorname{Sym}^{d}\mathbb{R}^{n}$ that can be written

[TABLE]

for some set of $n$ orthonormal vectors $v_{i}\in\mathbb{R}^{n}$, where $v^{\otimes d}$ denotes the $d$-wise tensor power of the vector $v$. That is, an odeco tensor encodes a set of orthogonal vectors up to permutation. If $d$ is even, $T$ is also invariant under sign changes to the $v_{i}$. Moreover, an odeco tensor assigns weights $\lambda_{i}$ to the vectors $v_{i}$. This property allows us to encode frames whose axes scale independently.

There is a one-to-one correspondence between symmetric tensors of order $d$ over $\mathbb{R}^{n}$ and homogeneous polynomials of degree $d$ in $n$ variables (cf [Robeva, 2016, §1.2]). This correspondence is given in one direction by taking a polynomial $p\in\mathbb{R}[x_{1},\dots,x_{n}]$ to its tensor of $d$th derivatives, and in the other direction by symbolic evaluation on the vector of formal variables $x=(x_{1},\dots,x_{n})$:

[TABLE]

This is a generalization of the correspondence between symmetric bilinear forms (i.e. symmetric $2$-tensors) and quadratic forms (i.e., homogeneous quadratic polynomials). Note that $T$ is odeco if and only if $p_{T}$ can be written as a sum of $d$th powers of linear forms

[TABLE]

where the $v_{i}$ are orthonormal as above. In this case, we also refer to $p_{T}$ as an odeco polynomial.

The defining equations of the odeco varieties are quadrics—homogeneous quadratic equations—in the tensor coefficients [Boralevi et al., 2017, Theorem 4] or equivalently in the coefficients of the associated polynomial $p_{T}$. That is, a homogeneous polynomial

[TABLE]

is odeco if and only if the coefficients $u$ satisfy

[TABLE]

for a finite set of symmetric matrices $A_{i}$. In the case relevant to us, where $n=3$ and $d=4$, there are exactly $27$ such defining equations. While dimension counting might suggest otherwise, these equations are not redundant—as can be seen by computing a Gröbner basis of the ideal they generate. The matrices $A_{i}$ are listed explicitly in the supplemental document. We will henceforth refer to this particular odeco variety simply as the odeco variety $\tilde{F}$. Figure 4 plots several fourth-order odeco polynomials over the unit sphere.

3.2.1. Octahedral Variety as a Subvariety of Odeco Variety

Recall that our octahedral frames were represented by coefficients in an irreducible representation of $\mathrm{SO}(3)$, while the odeco variety was defined using monomial or tensor coefficients. To see the relationship between the two varieties, it is beneficial to recast the odeco variety in the irreducible representation basis. This corresponds to looking at the coefficients of odeco polynomials in the basis of spherical harmonics.

The quartic polynomials comprising the odeco variety decompose as linear combinations of the spherical harmonics of bands [math], $2$, and $4$. Consider an odeco polynomial $\sum_{i=1}^{3}\lambda_{i}(v_{i}^{\mathstrut\scriptscriptstyle{\top}}x_{i})^{4}$, and let

[TABLE]

be its coefficients in this basis of even spherical harmonics. These $15$ coefficients give us a different representation of odeco frames in $\mathbb{R}^{15}$, where each band has a clear meaning. In particular, $q_{0}=C_{0}\sum_{i}\lambda_{i}$, where $C_{0}$ is a constant, and similarly

[TABLE]

The parenthesized expression is the squared distance between $(\lambda_{1},\lambda_{2},\lambda_{3})$ and the line $\lambda_{1}=\lambda_{2}=\lambda_{3}$. Thus $q_{2}=0$ if and only if $q_{4}$ is a scalar multiple of an octahedral frame. The set

[TABLE]

consists of scalar multiples of the octahedral variety indexed by $q_{0}$. Fix $q_{2}=0$ and let $q_{0}$ be a constant such that $\|q_{4}\|=1$. The octahedral variety is the intersection of this affine subspace with the odeco variety. Reducing the equations (6) of the odeco variety with respect to this subspace yields $15$ inhomogeneous quadratic equations

[TABLE]

cutting out the octahedral variety $F$. As was the case for the odeco variety, these equations are not redundant. The symmetric matrices $P_{i}$ are listed explicitly in the supplemental document.

4. Geodesics and Projection

We have introduced two spaces, the octahedral and odeco varieties, that can serve as target sets for frame fields. To compute such fields, we will have to solve optimization problems over products of many copies of these varieties. Naïvely, one might plug the quadratic constraints in (6) and (7) directly into a generic quadratic optimization solver. However, the $A_{i}$ or $P_{i}$ are not positive semidefinite, nor can their span be rewritten as the span of positive semidefinite matrices. So the constraints are nonconvex and challenging to enforce.

As an alternative, we use optimization algorithms that are tailored to the manifold-valued variables in our problem. These algorithms employ sub-steps such as geodesic traversal and projection. We derive these operations for a single frame below.

4.1. Octahedral Geodesics

By Proposition 3.6, the scaled octahedral variety $F_{\alpha}$ is locally isometric to $\mathrm{SO}(3)$ via the map $r\mapsto\rho(r)\cdot q_{\alpha}$. It follows that geodesics of $\mathrm{SO}(3)$ push forward to geodesics of $F_{\alpha}$. The relation (17) in Appendix A then allows us to compute geodesics on $F_{\alpha}$ in closed form without evaluating the representation map $\rho$ explicitly. Let $\mathfrak{v}\in T_{q}F_{\alpha}$. $\mathfrak{v}$ can be written in a basis induced by the $\mathrm{SO}(3)$ action, $\mathfrak{v}=\sum_{i=1}^{3}v_{i}L_{i}q$, Here the coefficient vector $v=(v_{1},v_{2},v_{3})\in\mathbb{R}^{3}$ is the “axis-angle” representation of a rotation, and the $\mathrm{SO}(3)$ exponential maps it to the corresponding rotation matrix. This can be computed by conjugation with a rotation $r$ taking the unit vector $v/\|v\|$ to $(0,0,1)$. Composing with the representation map, we have

[TABLE]

where (unsubscripted) $\exp$ denotes the ordinary matrix exponential.

To compute $\rho(r)$, define spherical coordinates $\theta,\varphi$ such that

[TABLE]

Then one possible choice for $r$ is

[TABLE]

where $r_{23}=\exp((\pi/2)l_{1})$. So

[TABLE]

where $R_{23}=\exp((\pi/2)L_{1})$. Combining (9) with (8), we can compute geodesics by products of two simple ingredients: $R_{23}$ and $\exp(tL_{3})$ for $t\in\mathbb{R}$. Closed-form expressions for both appear in the supplemental document, §2. Note that a formula similar to (9) is common in the graphics literature, e.g., for rotating BRDFs expressed in spherical harmonic coefficients (cf. [Kautz et al., 2002], Appendix).

4.2. Projection via Semidefinite Relaxation

$F$ and $\tilde{F}$ are both varieties defined by quadratic equations. $F\subset\mathbb{R}^{9}$ is cut out by $15$ inhomogeneous quadratic equations (7), while $\tilde{F}\subset\mathbb{R}^{15}$ is cut out by $27$ homogeneous quadratic equations (6). Consider the problem of finding the closest point in $F$ to a given point $y\in\mathbb{R}^{9}$:

[TABLE]

This problem is called Euclidean projection onto a quadratic variety. It is an example of a QCQP, and the general recipe for semidefinite relaxation detailed in §A.2 automatically applies. The SDP will have the form

[TABLE]

where $P_{i}$ are the symmetric matrices from (7), and

[TABLE]

Let $Q^{*}$ be an optimal solution to (SDPF). Since (SDPF) is a relaxation of (PF), $\langle Y,Q^{*}\rangle$ is a lower bound on the objective value of (PF). If it happens that $\operatorname{rank}(Q^{*})=1$, then $Q^{*}$ can be written in the form

[TABLE]

from which it follows that

[TABLE]

i.e., $q^{*}\in F$. In this case, $q^{*}$ is the globally optimal solution to (PF). This situation is called exact recovery. The foregoing discussion also applies, mutatis mutandis, to $\tilde{F}$.

In our MBO algorithm presented below (cf. §6.2), we alternate projection with stepping off the variety. We refer to the following theorem, which suggests that when taking small enough steps from smooth points of a variety, projection will be generically exact.

Theorem 4.1 (Cifuentes et al. [2017], Theorem 1.2).

Consider the Euclidean projection problem

[TABLE]

where $\mathcal{V}\subset\mathbb{R}^{n}$ is a real variety defined by quadratic equations $f_{1}(q)=\dots=f_{m}(q)=0$. Let $\bar{y}\in\mathcal{V}$ be a point at which the rank of the Jacobian $\nabla f(\bar{y})$ is equal to the codimension of $\mathcal{V}$. Then the semidefinite relaxation of projection is exact for $y\in\mathbb{R}^{n}$ sufficiently close to $\bar{y}$.

In addition to this theoretical motivation, we have ample empirical evidence that our relaxations are exact under much more general conditions. The blue histogram in Figure 5 shows the results of attempting to project $10^{6}$ random points onto $F$ via semidefinite relaxation. We use the ratio of the second largest eigenvalue $\lambda_{2}(Q^{*})$ to the largest eigenvalue $\lambda_{1}(Q^{*})$ as a proxy for exactness, as it indicates whether $Q^{*}$ was rank-one up to machine precision. For all octahedral projections, $\lambda_{2}(Q^{*})/\lambda_{1}(Q^{*})$ was $\approx 10^{-8}$ or less. Based on these results, we make the following conjecture.

Conjecture 4.2.

For a generic point $q_{0}\in\mathbb{R}^{9}$, the solution to (SDPF) has rank $1$, and therefore yields the exact projection $\Pi_{F}(q_{0})$.

We also tried projecting $10^{6}$ random quartic polynomials onto the odeco variety, as shown in the red histogram in Figure 5. The vast majority of odeco projections were also exact up to numerical precision: out of $10^{6}$ solutions, only $60$ had an eigenvalue ratio greater than $10^{-8}$. Theorem 4.1 gives us some intuition as to why this might happen. The stability result only holds near smooth points of the variety, but whereas the octahedral variety is a smooth manifold, the odeco variety has a singularity at the origin, separating polynomials of different signs.

Conjecture 4.3.

For a generic point $q_{0}\in\mathbb{R}^{15}$ representing a positive polynomial, the SDP relaxation yields the exact projection $\Pi_{\tilde{F}}(q_{0})$.

Indeed, the green histogram in Figure 5 shows the results of odeco projection on $10^{6}$ random positive quartic polynomials, generated as sums of squares of random quadratic polynomials. For all such positive initial points, the SDP solution had $\lambda_{2}(Q^{*})/\lambda_{1}(Q^{*})<10^{-8}$, supporting Conjecture 4.3.

For octahedral frames, we also compare our SDP-based projection to the previous state of the art method proposed by Ray et al. [2016]. Because Ray et al.’s method is based on nonconvex optimization, we would expect it to get stuck in local minima. Indeed, out of $100\,000$ trials, we found $600$ cases for which the result of Ray et al.’s method was at least $10^{-3}$ further from the initial point than our projection—a nearly $1\%$ error rate. Moreover, the difference between the computed projections can be substantial, as illustrated in Figure 6.

5. From Frames to Frame Fields

The overall aspiration of this work is to compute smooth frame fields in a region $\Omega\subseteq\mathbb{R}^{3}$ aligned to the normal $\boldsymbol{n}$ on its boundary $\partial\Omega$. We assume $\Omega$ to be compact with $\partial\Omega$ a union of smooth, embedded surfaces. We think of a frame field as a map $\phi:\Omega\to\mathcal{V}$ into some space of frames $\mathcal{V}$ satisfying alignment boundary conditions.

We’ve examined the geometry of two candidates for the fiber $\mathcal{V}$—the octahedral variety $F$ and the odeco variety $\tilde{F}$.222 The term “fiber” is intended to be suggestive. It may be fruitful to consider a bundle $\pi:P\to\Omega$ with fiber $\mathcal{V}$, whose restriction $\partial P\to\partial\Omega$ is nontrivial and encodes the boundary alignment condition. The map $\phi$ would be replaced with a section of $P$. It remains to describe the boundary conditions and to define an energy we want $\phi$ to minimize.

5.1. Boundary Conditions

In either $F$ and $\tilde{F}$, the frames aligned to a given direction $\boldsymbol{n}\in\mathbb{R}^{3}$ form the intersection of an affine subspace with the respective variety. In each case, this intersection is a lower-dimensional variety defined by a different set of quadratic equations. We can impose alignment boundary conditions by working over this lower-dimensional variety. For boundaries with sharp creases, we can optionally exclude creased vertices, where the normal direction is ill-defined, from the alignment constraint. Other constraints can also be imposed at creases—for example, alignment to the crease tangent.

5.1.1. Octahedral boundary conditions

Let $\rho:\mathrm{SO}(3)\to\mathrm{SO}(9)$ be the irreducible representation as in §3. First consider the octahedral frames aligned to the vertical $\boldsymbol{z}=(0,0,1)^{\mathstrut\scriptscriptstyle{\top}}$. The space of vertical-aligned frames $F_{\boldsymbol{z}}\coloneq\rho(A_{\boldsymbol{z}})q_{0}$, where $A_{\boldsymbol{z}}$ is the coaxial subgroup consisting of all rotations about $\boldsymbol{z}$, i.e.,

[TABLE]

and $q_{0}$ is the canonical octahedral frame. As derived in [Huang et al., 2011; Ray et al., 2016], such frames have the form

[TABLE]

where $q_{\boldsymbol{z}}=(0,0,0,0,\sqrt{7/12},0,0,0,0)^{\mathstrut\scriptscriptstyle{\top}}$, $s_{\boldsymbol{z}}=(\cos(4t),\sin(4t))^{\mathstrut\scriptscriptstyle{\top}}$, and $B_{\boldsymbol{z}}$ is the obvious $9\times 2$ matrix. Given $q\in\mathbb{R}^{9}$, the closest vertical-aligned frame is

[TABLE]

For a general unit normal $\boldsymbol{n}$, the aligned frames can be written

[TABLE]

where $r_{\boldsymbol{n}}\in\mathrm{SO}(3)$ is any rotation taking $\boldsymbol{z}$ to $\boldsymbol{n}$. The projection of $q\in\mathbb{R}^{9}$ onto the $\boldsymbol{n}$-aligned frames is then

[TABLE]

During optimization, this projection is used for boundary-aligned frames.

5.1.2. Odeco boundary conditions

Consider the standard odeco frame $\sum_{i}\lambda_{i}x_{i}^{4}$ rotated by some $r\in A_{\boldsymbol{z}}$, and fix $\lambda_{3}=1$. As described in §3.2.1, it has coefficients in the even-numbered irreducible representations (bands of spherical harmonics)

[TABLE]

Just as in the octahedral case, $q$ can be decomposed into a fixed part and a rotating part parametrized by a lower-dimensional vector

[TABLE]

where now $s_{\boldsymbol{z}}\in\mathbb{R}^{5}$ and $B_{\boldsymbol{z}}\in\mathbb{R}^{15\times 5}$. The equations (6) reduce to three quadratic equations in the coefficients of $s_{\boldsymbol{z}}$, which can be used to construct a semidefinite program to project onto the vertical-aligned odeco frames, as in §4.2. The details are given in the supplementary material. For frames aligned to a direction $\boldsymbol{n}$,

[TABLE]

where $\rho$ is now the direct sum representation on $V_{0}\times V_{2}\times V_{4}$. Projection onto the $\boldsymbol{n}$-aligned odeco frames can be constructed from projection onto the $\boldsymbol{z}$-aligned odeco frames as in (12).

5.2. Objective Function

As we are searching for smoothest fields, a natural choice for the energy is Dirichlet energy $\frac{1}{2}\int_{\Omega}\|\nabla\phi\|^{2}\;dx$, where the norm is taken with respect to a metric on $\mathcal{V}$. There are multiple problems this formulation will need to confront. As we have seen, $F$ is a smooth manifold. But being a quotient of the $3$-sphere, it has positive curvature, making mere existence of harmonic maps $\Omega\to F$ a hard problem. Even more fundamentally, it is not clear how to represent singularities of $\phi$—the map cannot be consistently defined along singular curves. $\tilde{F}$ attempts to solve this problem by representing singular frames with only one direction, while allowing the other two to vanish. Along a singular curve, we would expect an odeco field to take rank-one values of the form $\lambda v^{\otimes 4}$, where $v$ is tangent to the singular curve.

5.3. Discretization

Let $\mathcal{T}=(V,T)$ be a tetrahedral mesh of $\Omega$ with vertices $V$ and tetrahedra $T$. The set of boundary vertices will be denoted by $\partial V$. A discrete octahedral frame field on $\mathcal{T}$ is specified by a map $q:V\to\mathcal{V}=F$ or $\tilde{F}$, giving a frame $q_{i}$ for each vertex $i$. As the octahedral variety $F\subset\mathbb{R}^{9}$ and the odeco variety $\tilde{F}\subset\mathbb{R}^{15}$, we can think of $q$ as a $d\times n$ matrix, where $n=|V|$ and $d=9$ or $15$, respectively.

We then define a discretized Dirichlet energy using the Euclidean metric in the spherical harmonic basis to compare elements of $\mathcal{V}$. This is equivalent to measuring $L^{2}$ distance between the corresponding polynomials over the sphere, as employed in [Huang et al., 2011; Ray et al., 2016; Solomon et al., 2017]. Note that this would not be the case if we compared odeco frames in the monomial basis. The discrete energy is

[TABLE]

where $L$ is the linear finite element Laplacian on $\mathcal{T}$, with the usual cotangent weights.

6. Algorithms

We now have two variety-constrained optimization problems of the form

[TABLE]

where the variety $\mathcal{V}$ is either $F$ or $\tilde{F}$, $q_{i}$ are the columns of $q$, $\partial V$ denotes the set of boundary vertices, and $\mathcal{V}_{\boldsymbol{n}_{i}}$ is the variety of frames aligned to the normal at boundary vertex $i$.

6.1. Manifold Optimization Methods

In the case $\mathcal{V}=F$, (14) becomes a manifold-constrained optimization problem. The octahedral variety is a Riemannian manifold, and we can compute geodesics on it as described in §4.1. So we can apply the standard Riemannian trust-region (RTR) algorithm [Absil et al., 2007] as implemented in the Manopt toolbox for MATLAB [Boumal et al., 2014]. We consider the field to be a point in the product manifold

[TABLE]

In addition to the geodesics computed in §4.1, we also need a way to project vectors (e.g. gradients) in the ambient space $\mathbb{R}^{9}$ into the tangent space at a point $q$, $\operatorname{proj}_{q}:\mathbb{R}^{9}\to T_{q}F$. The vectors $\{L_{i}q\}_{i=1}^{3}$ form an orthogonal basis for this tangent space (cf. Proposition 3.6), so projection is just given by taking the inner product with each of these vectors.

The odeco variety $\tilde{F}$ is not a manifold, but it is smooth almost everywhere. We know of no way to compute exact geodesics on $\tilde{F}$, but we can implement a version of RTR by replacing the exact exponential map with a retraction, i.e., a first-order approximation (cf. [Absil et al., 2007, Definition 2.1]). In doing so, we give up the superlinear local convergence guarantees of RTR (cf. [Absil et al., 2007, Theorem 4.12]) but retain a global convergence guarantee. In practice, we find that odeco RTR converges at a reasonable rate.

At a smooth point of $\tilde{F}$, it is easy to project vectors onto the tangent space using the fact that $\tilde{F}$ is defined by quadratic equations. Let $q\in\tilde{F}$, and assume the coefficients of $q$ are expressed in spherical harmonic coefficients. Then differentiating the equations $q^{\mathstrut\scriptscriptstyle{\top}}\tilde{A}_{i}q=0$ shows that the normal space at $q$ is

[TABLE]

where $\tilde{A}_{i}$ are expressed in the spherical harmonic basis. Then $T_{q}\tilde{F}=(N_{q}\tilde{F})^{\perp}$, where the orthogonal complement is taken with respect to the standard metric under which the spherical harmonic functions are orthonormal.

We compute retractions as follows. The tangent space to the odeco variety decomposes into a rotational part and a scaling part:

[TABLE]

where the orthogonal complement is taken within $T_{q}\tilde{F}$. We then set

[TABLE]

where $v_{s}$ is the component of $v$ in $T^{s}_{q}\tilde{F}$, $v_{r}$ is the component of $v$ in $T^{r}_{q}\tilde{F}$, and

[TABLE]

where $v_{r}=\sum_{j}(v_{r})_{j}\tilde{L}_{j}q$. It is simple to verify that $\operatorname{retr}_{q}(v)\in\tilde{F}$ and that

[TABLE]

We compare the convergence behavior of octahedral RTR to that of the method of Ray et al. [2016] under two sets of conditions—Ray et al.’s initialization and combinatorial Laplacian (Figure 7); and our random initialization and linear finite element Laplacian (Figure 8). The quadratic local convergence of RTR stands in stark contrast to the slower linear convergence behavior of Ray et al.’s method. RTR converges faster but finds solutions of similar Dirichlet energy to previous work (see Table 1 in the supplemental document).

6.2. Generalized MBO Methods

As we have observed, it is possible to do optimization on the octahedral and odeco varieties by moving along curves that stay on the varieties exactly. However, this hard constraint sometimes causes pure manifold optimization to get stuck in local minima. To avoid such local minima, an approach that allows “tunneling” is required.

The method of Merriman, Bence, and Osher [1992] is a diffusion-generated method for computing (possibly singular) maps into a target space embedded in a Euclidean space. Following [Viertel and Osting, 2019; Osting and Wang, 2017], we first replace the hard constraint that the field values lie in the variety with a penalty term, yielding an energy of Ginzburg-Landau type,

[TABLE]

where $\mathcal{V}$ is our variety. Consider this energy in the limit $\epsilon\to 0$, as in [Viertel and Osting, 2019, §4]. The MBO method consists of alternating descent on the two terms of $E_{\epsilon}$. Gradient flow on the first term—Dirichlet energy—yields componentwise heat diffusion,

[TABLE]

with the boundary constraints given in §5.1. In practice, we do one step of implicit Euler integration per overall algorithm step. Gradient flow on the second term of (15) in the limit $\epsilon\to 0$ amounts to projection onto the variety. The overall algorithm is shown in Algorithm 1.

We follow the suggestion of [Viertel and Osting, 2019; van Gennip et al., 2014] to set $\tau_{0}$ relative to the inverse of the smallest nonzero eigenvalue of the Laplacian. In ordinary MBO, $\beta(k)=1$, i.e., $\tau$ does not change over the course of the algorithm. In Figure 9, we test the effects of different values of $\tau$. As we decrease $\tau$, the field is able to relax more and reduce the Dirichlet energy. On the other hand, larger values of $\tau$ allow more tunneling so that the field can escape local minima. This suggests a modified MBO scheme (mMBO) in which we start with a large $\tau$ and progressively reduce it over the course of optimization—analogously to what happens to the temperature in simulated annealing. Our mMBO uses an optimization schedule $\beta(k)=50k^{-3}$. This optimization schedule starts with a large diffusion time for robustness to random initialization, then sweeps $\tau$ through various orders of magnitude quickly. We have found that this heuristic produces a good balance between quick convergence and field quality.

Figure 10 depicts the convergence behavior of RTR, MBO with various $\tau$ values, and mMBO starting from a projection of an approximately harmonic field. While the energy value achieved by ordinary MBO is limited for each $\tau$, mMBO achieves a lower value by sweeping through various values of $\tau$. Note that the iteration counts shown do not reflect wall-clock time; in particular, RTR runs at least an order of magnitude faster than the other algorithms.

7. Experiments

We initialize our algorithms with vertexwise random octahedral fields, generated by—separately for each vertex—starting with the canonical frame, rotating it by a random angle between [math] and $2\pi$ about the positive $z$-axis, and then rotating the positive $z$-axis to a random direction. We do this even for optimization of odeco fields to avoid encountering odeco frames that have negative weights $\lambda_{i}$ (cf. Conjecture 4.3). When starting from octahedral initialization, the odeco frames we compute always have nonnegative weights.

In Figure 11, we demonstrate the robustness of MBO to random initialization. On the sphere, global rotations of a field do not change the objective value. Initializing odeco MBO with different fields of vertexwise random frames, we find that it converges to randomly rotated copies of a qualitatively identical, symmetric field.

Figure 12 illustrates the behavior of octahedral and odeco fields as the density of the underlying tetrahedral mesh increases. Due to the unit-norm constraint, the gradient of an octahedral field becomes unbounded close to its singularities. This leads to logarithmic divergence of the total energy as the tet mesh becomes finer. In contrast, the additional scaling degrees of freedom available to odeco fields allow renormalization of singularities, as shown in Figure 13. Thus, odeco field energy plateaus as mesh density increases. Note also the much smaller variance in energy between runs for odeco fields, quantitatively illustrating robustness to random initialization.

Table 1 (supplemental document) compares timings and energy values for our methods and the method of Ray et al. [2016], comprising $18$ types of experiments on $15$ different models, for $270$ total runs. Direct comparison of energy values with Ray et al. [2016] is not possible, as their method uses the graph Laplacian and initializes by solving a linear system, whereas our method uses the geometric (finite element) Laplacian and random initialization. To attempt a fair comparison, we report results of their method alone, their method followed by RTR with the geometric Laplacian, and their method substituted with the geometric Laplacian and random initialization. All experiments in the table were conducted on an Ubuntu workstation with a four-core, 3.6 GHz Intel Core i7-7700 and 32 GB of RAM. A MATLAB implementation of our algorithms accompanies this paper and also includes our implementation of [Ray et al., 2016].

The results in Table 1 show that RTR converges very quickly but, like previous work, easily gets stuck in local minima. In contrast, mMBO takes longer to converge but produces higher-quality fields that reflect the symmetries of the volume. Running RTR to polish the results of mMBO produces the best energies.

Despite sometimes getting stuck in local minima (see Figure 6), we observe that Ray et al.’s approximation of projection can be useful in practice. Table 1 (supplemental document) includes experiments with MBO and mMBO (see §6.2) substituting Ray et al.’s projection for the true projection. In most cases, the resulting energy is very similar, suggesting that the diffusion iterations in MBO smooth out any errors resulting from incorrect projection. This hybrid MBO can be a useful tradeoff between correctness and speed.

In Figures 14–17, we compare fields computed by octahedral mMBO + RTR to those computed by the method of Ray et al. [2016]. Our fields not only have lower Dirichlet energy, but also better singular structures. To visualize singular structures, we use the visualization technique of Liu et al. [2018]. To illustrate the importance of singular structure, we have generated hexahedral meshes from both sets of fields, following the methods laid out in [Nieser et al., 2011] and [Lyon et al., 2016]. The meshes are computed from the raw field data, with no preprocessing other than tetrahedral mesh refinement to resolve and localize singular curves. In particular, we do not “correct” singularities—thus, both sets of meshes show some degeneracies resulting from collapsed or flipped tetrahedra during the parametrization step. However, our fields yield meshes with fewer, smaller degeneracies, less distortion, and more regular structure.

Figures 18 and 19 compare fields produced by mMBO + RTR to fields produced by the method of Gao et al. [2017]. Our fields better respect the symmetries of the underlying models. We note that mMBO + RTR yields a higher-quality result without requiring a multiscale method like the one Gao et al. [2017] employ.

8. Discussion and Future Work

A stronger understanding of the unknowns in the volumetric frame field problem enables both theoretical and practical developments. On the theoretical side, our reexamination of the typical representation of octahedral frames not only yields useful geodesic formulas and projection operators but also suggests generalization to odeco frames. For both frame representations, optimization algorithms designed explicitly to traverse the manifold of unknowns yield gains in efficiency and in the quality of computed results.

Our work is intended not only to improve on existing frame field computation algorithms but also to inspire additional research into the structure of spaces of frames and to highlight the significance of this structure to practical aspects of field computation and meshing. Most immediately, a few open questions are left from our discussion. On the theoretical side, conjectures 4.2 and 4.3 remain to be proven. We anticipate that both can be embedded in a more general framework explaining when relaxations of projection problems are exact. On the algorithmic side, RTR seems to get stuck in local minima much more often than MBO, especially on denser meshes. We hypothesize that this is because RTR strictly moves along the manifold, while MBO is free to tunnel through the ambient space. On the other hand, RTR converges much faster and yields high-quality fields on coarse meshes. Given these observations, we anticipate that it may be possible to incorporate RTR into a multiscale method that leverages its efficiency while avoiding local minima that appear at the finest scales.

As with most existing frame field computation algorithms, even when mMBO+RTR converges to a smooth field, the topology of the field is not always hex-meshable. Figure 20 shows a case where the method of Ray et al. [2016] yields a field of higher Dirichlet energy, but which leads to fewer degeneracies than our field. In particular, the presence of a singular curve whose type changes from valence $3$ to $5$ leads to a degeneracy in the integer-grid map. While our method succeeds at finding fields of lower Dirichlet energy—the objective function of our method and that of Ray et al. [2016]—minimizing this energy is an incomplete proxy for our ultimate goal, namely to obtain smooth, meshable fields. We would like to investigate additional metrics, e.g., integrability of frame fields, that might make it possible to express meshability constraints rigorously.

To define such additional metrics, it might be fruitful to consider further alternative frame field representations. For example, given our use of Lie algebra representations, a logical next step might be to introduce an $\mathrm{SO}(3)$–principal bundle and work with connections on that bundle as variables. In this theory, quantities such as curvature, torsion, or the Chern-Simons functional might encode important features of frame fields. The discretization of connections in [Corman and Crane, 2019] represents a promising first step.

Finally, while our algorithms do not explicitly take account of symmetry, we find that fields computed by MBO consistently reproduce the symmetries of the volume. It would be interesting to develop a better theoretical understanding of this behavior and to develop machinery for explicitly promoting conformation to near-symmetry in deformed domains.

These and other challenges appear when extending well-known machinery from geometry processing on surfaces to volumetric domains, as demanded by applications in engineering, simulation, and other disciplines. Over the longer term, careful consideration of problems like the ones we consider in this paper will lead to a versatile and generalizable approach to geometry processing.

Appendix A Mathematical Background

Representation theory and algebraic geometry provide a concise language for the discussion of the octahedral and odeco varieties above. We begin with definitions of a few terms used throughout the paper. While a detailed introduction to either of these subjects is outside the scope of our discussion, we refer the reader to [Stillwell, 2008] for an introduction to matrix Lie groups and to [Blekherman et al., 2012] for a primer on real algebraic geometry and optimization.

A.1. Lie Groups and Representations

A Lie group is a group $G$ that is also a smooth manifold, and such that the multiplication $\cdot:G\times G\to G$ and inversion $\,{}^{-1}:G\to G$ are smooth maps. A matrix Lie group is a subgroup of the invertible matrices $\mathrm{GL}(n)\subset\mathbb{R}^{n\times n}$ that is simultaneously a submanifold. Examples include the orthogonal groups $\mathrm{O}(n)$ and the special orthogonal groups $\mathrm{SO}(n)$. The primary Lie group of interest to us in this paper will be $\mathrm{SO}(3)$.

For a matrix Lie group $G\subset\mathbb{R}^{n\times n}$, the Lie algebra is a linear subspace of $\mathbb{R}^{n\times n}$ identified with the tangent space to $G$ at the identity $I$, denoted $\mathfrak{g}=T_{I}G$. The inner automorphism $\operatorname{Ad}_{g}:G\to G$ taking $h\mapsto ghg^{-1}$ preserves the identity, and its derivative at the identity induces the Lie bracket $[\cdot,\cdot]:\mathfrak{g}\times\mathfrak{g}\to\mathfrak{g}$. In the case of a matrix Lie group, this is just the matrix commutator.

As left translation $\operatorname{L}_{g}:h\mapsto gh$ is a diffeomorphism, every element of $\mathfrak{g}$ is also associated to a left-invariant vector field on $G$. The exponential map $\exp:\mathfrak{g}\to G$ is given by starting at $I$ and integrating this vector field for one unit of time. For a matrix Lie group, $\exp$ is the ordinary matrix exponential.

A representation of a matrix Lie group $G$ is a smooth group homomorphism $\rho:G\to\mathrm{GL}(n)$ for some $n$. A representation of $G$ induces a representation of each subgroup $H\subset G$. A Lie group representation comes with a Lie algebra representation $(D\rho)_{I}:\mathfrak{g}\to\mathbb{R}^{n\times n}$, which preserves the Lie bracket. Moreover, representations commute with exponentials:

[TABLE]

An irreducible representation is one that cannot be decomposed as a direct sum of sub-representations. The irreducible representations of compact Lie groups such as $\mathrm{SO}(3)$ are finite-dimensional. The adjoint representation is the irreducible representation of $G$ on its own Lie algebra given by conjugation:

[TABLE]

Lemma A.1.

$\operatorname{Ad}$* commutes with representations of $G$:*

[TABLE]

Proof.

Using (17) and conjugating by $g$, we obtain

[TABLE]

Differentiating in $t$ at $t=0$ yields the result. ∎

The stabilizer $\operatorname{stab}(v)$ of a vector $v$ is the subgroup consisting of all elements $g\in G$ that preserve $v$—that is, such that $\rho(g)v=v$.

The Lie algebra $\mathfrak{so}(3)$ associated to $\mathrm{SO}(3)$ consists of the skew-symmetric $3\times 3$ matrices. $\mathfrak{so}(3)$ has a basis consisting of infinitesimal rotations about the three coordinate axes:

[TABLE]

Their Lie brackets are

[TABLE]

These generators might be familiar as the angular momentum operators from quantum mechanics. Indeed, $\mathrm{SO}(3)$ acts on functions on the sphere by rotating them, and its Lie algebra elements act as rotational derivatives. This action decomposes into irreducible representations. A basis for each irreducible representation of $\mathrm{SO}(3)$ is provided by spherical harmonics; each irreducible representation is referred to as a band of spherical harmonics. The real representations are odd-dimensional vector spaces, and the representation matrices $\rho(\cdot)$ are also called Wigner $D$-Matrices.

$\mathrm{SO}(3)$ admits a bi-invariant Riemannian metric whose value at the identity is given by

[TABLE]

Note that the generators $l_{i}$ form an orthonormal basis under this metric. Bi-invariance means that

[TABLE]

where $v,w\in\mathfrak{so}(3)$ and $\operatorname{R}_{g}$ is right translation by $g\in SO(3)$, defined analogously to left translation. (23) completely describes the Riemannian metric on $\mathrm{SO}(3)$. Under this metric, the matrix exponential is also the Riemannian exponential map at the identity, giving rise to geodesics.

The octahedral group $\mathrm{O}$ is a discrete subgroup of $\mathrm{SO}(3)$ generated by right-angle rotations about the three axes,

[TABLE]

$\mathrm{O}$ has $24$ elements comprising all the symmetries of the cube, or equivalently, of its dual, the octahedron.

A.2. Semidefinite Relaxations

A (real) algebraic variety is the set of solutions of a system of polynomial equations over $n$ real variables:

[TABLE]

where $p_{1},\dots,p_{k}$ are polynomials.

If $p_{1},\dots,p_{k}$ are quadratic in $x$, they may be written

[TABLE]

where $A_{i}$ are symmetric matrices, $\langle,\rangle$ denotes the usual inner product on matrices inducing the Frobenius norm, and

[TABLE]

A quadratically-constrained quadratic program (QCQP) is of the form

[TABLE]

where $f,p_{1},\dots,p_{k}$ are all quadratic. A QCQP may be rewritten in the more illuminating form

[TABLE]

Absent the rank constraint, (24) would be a semidefinite program (SDP). Semidefinite programs may be solved in polynomial time by interior point methods, implemented in common optimization software packages like Mosek [2020]. Thus it is natural to ignore the rank constraint and solve the associated SDP; this technique is called semidefinite relaxation. The optimal objective value of the relaxed problem is a lower bound for the globally optimal value of the QCQP (24), and if the recovered solution to the SDP has rank one, the solution is said to be exact, as it is the globally optimal solution of (24).

Acknowledgements.

The authors would like to thank Elina Robeva for helping us understand her work on the odeco varieties, Pablo Parrilo and Diego Cifuentes for sharing their insights into semidefinite relaxation, and Nicolas Boumal for helping us get started with manifold optimization. We thank Heng Liu for providing tetrahedral refinement, parametrization, and hexahedral mesh extraction code. We are grateful to Xifeng Gao and Daniele Panozzo for guiding us in using their code from [Gao et al., 2017]. Finally, we thank Mirela Ben-Chen, Amit Bermano, Edward Chien, Etienne Corman, Keenan Crane, Fernando de Goes, Steven Gortler, Leif Kobbelt, Braxton Osting, Nir Sharon, Nir Sochen, Amir Vaxman, Josh Vekhter, and Paul Zhang for valuable discussions. David Palmer appreciates the generous support of the Fannie and John Hertz Foundation through the Hertz Graduate Fellowship. David Bommes appreciates the generous support of the Sponsor European Research Council (ERC) https://erc.europa.eu under the European Union’s Horizon 2020 research and innovation program (AlgoHex, grant agreement no. Grant #853343). Justin Solomon acknowledges the generous support of Sponsor Army Research Office https://www.arl.army.mil/who-we-are/aro/ grant Grant #W911NF-12-R-0011, Sponsor National Science Foundation https://www.nsf.gov grant Grant #IIS-1838071, Sponsor Air Force Office of Scientific Research https://www.wpafb.af.mil/afrl/afosr/ award Grant #FA9550-19-1-0319, and a gift from Adobe Systems. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the authors and do not necessarily reflect the views of these organizations.