Octahedral Frames for Feature-Aligned Cross-Fields

Cross fields ($n=4$) have been especially well-studied since their $\frac{\pi}{2}$-symmetry allows them to behave like local coordinate systems, resulting in intuitive seamless surface parameterizations.

Methods to compute intrinsically-smooth cross fields with alignment and singularity constraints were studied by Ray et al. [2008], Crane et al. [2010] and Knöppel et al. [2013]. More similar to our work, Jakob et al. [2015] instead formulate an extrinsic smoothness functional on cross fields, in an attempt to automatically align to surface features. Their method penalizes an extrinsic distance between neighboring crosses that does not use a shared tangent space or connection. The resulting energy is non-convex but is minimized to local optimality, often resulting in more singularities than necessary. Huang and Ju [2016] analyze this extrinsic energy, and find that it can be decomposed into an energy expressed in terms of intrinsic twisting and alignment to extrinsic curvature directions. We perform a similar analysis of the energy we introduce in the supplemental documents.

When using cross fields for quad mesh parameterization or processing, methods [Bommes et al., 2009; Campen et al., 2016; Knöppel et al., 2013; Brandt et al., 2018] often promote feature alignment by including a loss term penalizing disagreement with curvature directions. However, as we argue in the introduction and illustrate in Figure 2, alignment to curvature directions is often less important than alignment to sharp creases. Other parameterization methods such as [Bommes et al., 2009; Bommes et al., 2013; Campen et al., 2015] allow feature alignment but just assume that such feature curves are provided as input.

The three-dimensional generalization of a cross field is an octahedral field. Octahedral fields are often used in volumetric problems like hexahedral meshing [Nieser et al., 2011]. A single octahedral frame consists of three mutually-orthogonal vectors and their negations. Huang et al. [2011] introduced a particularly convenient representation of an octahedral frame as a rotation of the spherical function $g(x,y,z):(x,y,z)\in S^{2}\mapsto x^{4}+y^{4}+z^{4}$, encoded by coefficients in the spherical harmonic (SH) basis. Ray et al. [2016] used this representation to generate volumetric normal-aligned octahedral fields, and Solomon et al. [2017] combined the SH representation with the boundary element method to remove the need for a volumetric mesh. Both of these methods use normal-alignment constraints at the surface to enforce alignment of the octahedral frame with the volume’s boundary.

While there is no canonical three-dimensional generalization of arbitrary $n$-RoSy fields, the spherical harmonic representation allows for frames that mimic the symmetries of all platonic solids [Shen et al., 2016], including octahedral fields  [Solomon et al., 2017; Liu et al., 2018; Corman and Crane, 2019]. Algebraic characterization of the orbit of $g(x,y,z)$ under the space of rotations as a subset of all possible SH coefficients was presented by Palmer et al. [2019] and Chemin et al. [2018].

As introduced by Huang et al. [2011], the canonical axis-aligned octahedral frame can be represented by spherical harmonics as a function $g_{0}:\mathbb{S}^{2}\to\mathbb{R}$ written as $g_{0}=\sqrt{\frac{5}{12}}Y_{44}+\sqrt{\frac{7}{12}}Y_{40},$ where $Y_{lm}$ denotes the basis for real spherical harmonics. The function $g_{0}$ can be understood as the scaled projection of $x^{4}+y^{4}+z^{4}$ onto the fourth band ($l=4$) of spherical harmonics. Written differently, we can encode $g_{0}$ as a vector of coefficients in the full basis of fourth-band spherical harmonics $Y_{4(-4)},\ldots,Y_{44}$:

The space of octahedral frames can be described as all rotations of the canonical octahedral frame, that is, the orbit of $f_{0}$ under the group of 3D rotations $\mathrm{SO}(3)$ [Palmer et al., 2019, Definition 3.5]. We write this via exponentiation of the Lie algebra elements: the set of octahedral frames is

where $v\cdot L=v_{x}L_{x}+v_{y}L_{y}+v_{z}L_{z}$ and $L_{x},L_{y},L_{z}$ are the angular momentum operators expressed in the basis of band-four spherical harmonics. In this basis, $L_{x},L_{y},L_{z}$ are each $9\times 9$ matrices. The angular momentum operators are explicitly written in supplementary materials, §4. In this description, $v$ can be interpreted as an axis-angle representation of rotation, with corresponding rotation matrix $e^{[v]}$. $[v]$ denotes the skew-symmetric matrix that acts as $[v]u=v\times u$. Accordingly, $e^{v\cdot L}g_{0}$ encodes the octahedral frame whose directions are $\hat{x},\hat{y},\hat{z}$ rotated by $e^{[v]}$, where $\hat{}$ denotes normalization (see inset above).

Using such SH rotations, we can present an alternative interpretation of the octahedral frame $f_{0}$ as the sum of three orthogonal SH lobe-shaped functions. The $z$-aligned lobe is $l=[0,0,0,0,\sqrt{\frac{7}{12}},0,0,0,0]$ and is depicted in the inset. Lobes can be rotated in the same way that frames can i.e. by applying $e^{v\cdot L}$. The canonical octahedral frame can therefore be equivalently expressed as $f_{0}=l+e^{\frac{\pi}{2}L_{x}}l+e^{\frac{\pi}{2}L_{y}}l$.

The space of octahedral frames that are aligned to a unit vector $\hat{n}$ can be described by the set

where $v_{n}$ is any axis-angle rotation taking $\hat{z}$ to $\hat{n}$, (e.g., the vector parallel to $\hat{z}\times\hat{n}$ and has magnitude equal to the angle from $\hat{z}$ to $\hat{n}$) and $\theta$ encodes an additional twist of the frame about $\hat{n}$. The first rotation about $\hat{z}$ can be written in explicit form [Huang et al., 2011] as

The above allows us to formulate the set of all octahedral frames $g$ aligned to a given direction $\hat{n}$ in terms of two constraints:

where $W_{n}$ is the second through eighth rows of $e^{-v_{n}\cdot L}$. The linear constraint rotates the frame from normal alignment to $\hat{z}$ alignment, and the norm constraint ensures that the first and last components are of the appropriate form.

Lastly, we will make use of the projection operator $\pi_{\mathcal{V}}:\mathbb{R}^{9}\rightarrow\mathcal{V}$ onto the space of octahedral frames $\mathcal{V}$, as defined in [Palmer et al., 2019, §5.5].

We will later make use of a total variation energy (amongst others) to analyze the behavior of our cross fields on creased surfaces. Here, we introduce total variation and vectorial total variation definitions in $\mathbb{R}^{n}$ and provide intuition about their use. The extension to functions on a Riemannian manifold is straightforward, using the standard intrinsic gradient and divergence operators.

The total variation of a differentiable scalar function $h:\Omega\rightarrow\mathbb{R}$ is $TV[h]=\int_{\Omega}\|\nabla h\|_{2}~dA$ where $\Omega\subset\mathbb{R}^{n}$ [Ambrosio et al., 2000]. For non-differentiable $h$, the relevant definition is:

where $C^{1}_{c}$ denotes differentiable, compactly-supported vector fields. For smooth $h$, equivalence to $\int_{\Omega}\|\nabla h\|$ follows from integration by parts and Stokes’s theorem. In this case, the maximizing $\phi$ is $\nicefrac{{-\nabla h}}{{\|\nabla h\|}}$. If $h$ is the indicator function of a suitably regular (e.g., non-fractal) subset $A\subset\Omega$, then $TV[h]$ is the perimeter of $A$.

When $h:\Omega\rightarrow\mathbb{R}^{m}$ is vector-valued rather than scalar-valued, there are many different definitions for the vectorial total variation $VTV[h]$ [Sapiro, 1996]. We use one proposed by Di Zenzo [1986], which for differentiable $h$ is given by $VTV[h]=\int_{\Omega}\|\nabla h\|_{F}~dA$, where $\|\cdot\|_{F}$ is the Frobenius norm. More generally, we can take

where $h=(h_{1},h_{2},\ldots,h_{m})$, and $\phi=(\phi_{1},\phi_{2},\ldots\phi_{m})$ is a differentiable, compactly-supported $m$-tuple of vector fields. This definition is not equivalent to a sum of $m$ independent scalar total variations: The constraint on $\phi$ introduces nontrivial coupling between the dimensions. This definition of total vectorial variation is considered in the case where $\Omega$ is a surface in $\mathbb{R}^{3}$ by Bresson and Chan [2008], but without specific analysis for discontinuous $h$.

To calculate $\|\nabla f\|^{2}$, we first express it in an appropriate local coordinate system that simplifies the formulas in coordinates and better reveals the structure. Following the notation in §3.1, an octahedral field $f(r):\Omega\rightarrow\mathcal{V}\subset\mathbb{R}^{9}$ can be parameterized relative to a point $r^{*}$ by $v(r):\Omega\rightarrow\mathbb{R}^{3}$, where $v(r)$ is the axis-angle rotation from $f(r^{*})$ to $f(r)$. This implies that $v(r^{*})=[0,0,0]$. Without loss of generality, we rotate the surface so that the normal of $\Omega$ at $r^{*}$ is $\hat{z}$. We can then compute the gradient $\nabla f$ at the point $r^{*}$ from the formula $f(r)=e^{v(r)\cdot L}f(r^{*})$:

As the field $f(r)$ encodes an extrinsically embedded frame at each point, we take the gradient $\nabla$ to be the component-wise derivative of the field’s nine scalar functions rather than a covariant or Lie derivative along the surface in order to capture the extrinsic geometry of the surface. We use [Rossmann, 2002, §1.2.5] and the fact that $v(r^{*})=[0,0,0]$ to derive Equation (3).

By combining facts about the SH representation and standard results in differential geometry we show that the squared norm $\|\nabla f(r)\|^{2}$ at $r^{*}$ can then be expressed in the following more intuitive way:

We leave the full proof of this formula to supplementary materials. Proposition 4.1 gives a more intuitive form for Equation (3) and relates the spherical harmonic representation of an octahedral frame to properties of the frame it represents. Most notably, the Dirichlet energy of the SH representation can be effectively decoupled into extrinsic dependence of $\|\nabla f\|^{2}_{F}$ on the surface $\Omega$ and the intrinsic tangential twisting of the normal-aligned octahedral field $f(r)$. The values of $H$ and $K$ simply contribute a fixed quantity depending on $\Omega$ rather than the field. Therefore, the influence of $f$ on $\|\nabla f\|^{2}_{F}$ is just in $w$, the intrinsic twist of the cross field it represents. We stress that this behavior is quite different from the behavior of the component-wise derivative evaluated on vectors, as studied in [Huang and Ju, 2016], where their smoothness energy promotes alignment to extrinsic curvature directions.

Suppose we wish to measure smoothness of a normal-aligned octahedral field in the SH representation. We define the following class of convex smoothness energies using the $L^{p}$-norm of $\|\nabla f\|_{F}$ over the surface $\Omega$ for $p\geq 1$:

We now analyze the behavior of the $E_{p}$ energy for cross fields in several select cases.

It is sometimes beneficial to relax the normal alignment constraint, e.g., in cases where the mesh contains sliver triangles with unstable normal directions. In these cases, a smoother cross field can be obtained by deviating slightly from exact normal alignment. This relaxation changes the optimization problem from Equation (7) into the following:

This problem imposes a point-wise normal alignment constraint with tolerance $\epsilon$. When $\epsilon=0$, we recover the hard normal alignment formulation (7). On the opposite side of the spectrum, as $\epsilon\rightarrow\|u_{0}\|_{2}=\sqrt{\frac{7}{12}}\approx\;0.76$, the solution to (8) approaches a constant octahedral field. This is the case where normal alignment has relaxed so far that the octahedral frames are effectively unconstrained.

For values in between we perform the following experiment to obtain a rough correspondence between soft normal alignment parameter $\epsilon$ and maximum angle deviation from normal alignment: For each value of $\epsilon$ between $0$ and $.7$ (at intervals of $.05$), we sample $100000$ $\epsilon$-perturbations of a $\hat{z}-$aligned frame, extract the frame they represent, and compute its maximum angle deviation from the $\hat{z}$-axis. Results are shown in the inset.

We highlight that this parameter encodes a point-wise constraint uniformly applied over the mesh. As such its interpretation does not change with different meshes. Please see the supplemental materials for results on over 200 different meshes using a variety of values of $\epsilon$.

The benefit of soft normal alignment is demonstrated in Figure 5. Due to the influence of a sliver triangle in the buste mesh with unstable normal direction, the hard-normal-aligned cross field is forced to create a localized artifact. By using soft normal alignment, the sliver triangle’s unstable normal direction has less influence over the resulting cross field, therefore increasing the quality of the result. A similar benefit is demonstrated on the duck and armchair meshes shown in the supplementary materials.

Additionally we test soft normal alignment on a cube-mesh with artificial noise added in Figure 7. With hard normal alignment the cross fields exhibit undesirable alignment to noise which increases with $p$. With soft-normal alignment, the cross fields show significantly decreased sensitivity to noise.

Now we describe how to construct smooth cross-fields by numerically optimizing a discretization of $E_{p}$. We assume the surface $\Omega$ has been triangulated into a manifold mesh $\mathcal{M}=\left(V,E,F\right)$. Let $n_{t}$ be the normal direction of triangle $t\in F$. We represent a cross on $M$ as a normally-aligned octahedral frame $f_{t}\in\mathcal{V}\subset\mathbb{R}^{9}$ per triangle. We use the shorthand $f$ to denote the concatenation of all $f_{t}$ into a single $9|F|\times 1$ vector. $V$ is a $n_{v}\times 3$ matrix of vertex positions, where $n_{v}$ is the number of vertices. $E$ denotes the $n_{e}\times 2$ matrix of edges, where $n_{e}$ is the number of edges. The energy $E_{p}$ can be discretized as

where $t_{1}$ and $t_{2}$ are triangles adjacent to edge $e$, and $w_{e}$ are weights corresponding to the dual Laplacian. We use $w_{e}=\frac{\|e\|}{\|e^{*}\|}$, where $\|e\|$ is the length of edge $e$ and $\|e^{*}\|$ is the distance between barycenters of $t_{1}$ and $t_{2}$.

For $\epsilon=0$, the normal alignment constraint is discretized by the linear constraint $Wf=u$, where $W$ is a sparse block-diagonal matrix with a block $W_{n_{t}}$ for each triangle. It has dimensions $7|F|\times 9|F|$. The vector $u$ is a repetition of $u_{0}$ for each triangle, resulting in a $7|F|\times 1$ vector. For $\epsilon>0$, the normal alignment constraint is discretized by a second-order cone constraint: $\|W_{n_{t}}f_{t}-u_{0}\|_{2}\leq\epsilon$ per triangle.

For the case $p=1$ and a completely flat surface $\Omega$, our discretization agrees with the standard discretization of total variation in image processing [Rudin et al., 1992; Chambolle et al., 2010].

When the mesh contains boundaries, we do not enforce any boundary conditions; in PDE parlance, this choice corresponds to “natural boundary conditions.” While it is tempting to expect the natural boundary conditions for $p=2$ to imply zero Neumann boundary conditions [Stein et al., 2018], the SH representation vector is complicated by being constrained to a spatially varying linear subspace. We simply allow the cross on the boundary to be that which minimizes total energy. If desired, one can enforce a constraint that the cross field on the boundary be aligned to the boundary through the method described in §5.4.

To support manual guidance of the octahedral frame field, we can prescribe alignment of the frame field to streamlines. Streamline constraints combined with normal alignment result in a fully-determined frame. Therefore, prescribing streamlines is equivalent to prescribing the value of $f_{t}$ on a subset of triangles $T_{p}$. Denote the prescribed octahedral frame on triangle $t$ as $F_{t}$. We then add a new linear constraint that

This technique is demonstrated in Figure 6.

As a result of dropping the unit-norm constraint from Equation 7, we have no explicit guarantee that the tangential components of octahedral frames do not degenerate to zero. On a surface with a crease, however, the normal alignment constraint on one side of the crease imposes that the magnitude of the tangential component on the other side of the crease is close to one. As a result, we observe empirically that the vast majority of our octahedral frames do not degenerate.

In the case that octahedral frames do degenerate significantly, their norms can be too small to project robustly. We locate these by using the octahedral projection from [Palmer et al., 2019] to measure the distance from $f_{t}$ to the octahedral variety $\mathcal{V}$ : $d(f_{t})=\|\pi_{\mathcal{V}}(f_{t})-f_{t}\|_{2}$, and thresholding by $d(f_{t})>.665$. If such frames are found, we run the optimization again while holding non-degenerate frames to their projected values. In our experiments, just one round of re-solving results in 99.8% non-degenerate frames.

In its most general form, our problem formulation consists of minimizing a mixed-norm objective, with both linear and second-order cone constraints. This results in a convex problem that we solve with Mosek 9  [ApS, 2017]. The normal alignment constraint becomes $\left[\epsilon,(W_{n_{t}}f_{t}-u_{0})^{T}\right]\in\mathcal{L}^{8}$, where $\mathcal{L}^{8}$ is the 8-dimensional Lorentz cone. Likewise, the energy is formulated using a single $p$-norm cone. Our code is written in Matlab with a mex interface to Mosek; it builds cross-platform. Since our problem is convex, any dependence on initialization would entirely be due to non-unique solutions, which we do not observe in practice. Furthermore, we use the interior point method, which does not accept manual initialization. In the specific case of $\epsilon=0,p=2$, solving this optimization is equivalent to solving a linear system.