<figure><img src="image_1.png" /></figure>

HAL
open science

# Rectangular Surface Parameterization

Etienne Corman, Keenan Crane

► To cite this version:

Etienne Corman, Keenan Crane. Rectangular Surface Parameterization. ACM Transactions on Graphics, 2025, 44 (4), pp.1-21. 10.1145/3731176. hal-05189504

HAL Id: hal-05189504

https://hal.science/hal-05189504v1

Submitted on 28 Jul 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Rectangular Surface Parameterization

ETIENNE CORMAN, Université de Lorraine, CNRS, LORIA, France

KEENAN CRANE, Carnegie Mellon University, USA

This paper describes a method for computing surface parameterizations that map infinitesimal axis-aligned squares in the plane to infinitesimal rectangles on the surface. Such rectangular parameterizations are needed for a broad range of tasks, from physical simulation to geometric modeling to computational fabrication. Our main contribution is a novel strategy for constructing frame fields that are perfectly orthogonal and exactly integrable, in the limit of mesh refinement. In contrast to past strategies for achieving integrability, we obtain maps that are less distorted and better preserve target field directions. The method supports user-defined distortion measures, sharp feature alignment, prescribed or automatic cone singularities, and direct control over boundary behavior (e.g., sizing or aspect ratio). By quantizing and contouring these maps we obtain high-quality anisotropic quad meshes, even without element-based optimization. Empirically, we outperform state-of-the-art research and commercial mesh generation algorithms in terms of element quality, accuracy, and asymptotic convergence rate in end-to-end simulation tasks, are competitive with the widely-used ZBrush package for automatic retopology, and provide Chebyshev nets of superior quality to methods specifically tailored to digital fabrication.

CCS Concepts: • Mathematics of computing → Mesh generation.

Additional Key Words and Phrases: surface parameterization, quad meshing

**ACM Reference Format:**

Etienne Corman and Keenan Crane. 2025. Rectangular Surface Parameterization. *ACM Trans. Graph.* 44, 4 (August 2025), 21 pages. https://doi.org/10.1145/3731176

## 1 Introduction

*Definition.* A surface parameterization $\phi(x, y) : U \to \mathbb{R}^n$ over a region $U \subset \mathbb{R}^2$ is rectangular if $\frac{\partial \phi}{\partial x} \cdot \frac{\partial \phi}{\partial y} = 0$, i.e., if the coordinate axes get mapped to orthogonal directions on the surface. In differential geometry, such a map is known as an *orthogonal parameterization* [Do Carmo 2016, Section 2.5]; we opt for “rectangular” since “orthogonal” is often overloaded to mean *orthonormal* (orthogonal and unit length). Equivalently, $\phi$ is rectangular if its Jacobian can be expressed at each point as $J_{\phi} = QD$ where $Q$ is orthogonal and $D$ is diagonal, describing an axis-aligned scaling followed by a rotation. At each point, the induced Riemannian metric $g$ is hence represented by a positive diagonal matrix $J_{\phi}^{T} J_{\phi} = D^{T} Q^{T} QD = D^{2}$, generalizing conformal maps, where it is a positive multiple of the identity. Importantly, not all parameterizations map (axis-aligned) squares to rectangles—not even infinitesimally (Figure 2).

Authors' Contact Information: Etienne Corman, etienne.corman@cnrs.fr, Université de Lorraine, CNRS, LORIA, Vandœuvre-lès-Nancy, France; Keenan Crane, kmcrane@cs.cmu.edu, Carnegie Mellon University, Pittsburgh, PA, USA.

© 2025 ACM.

This is the author's version of the work. It is posted here for your personal use. Not for redistribution. The definitive Version of Record was published in *ACM Transactions on Graphics*, https://doi.org/10.1145/3731176.

<figure><img src="image_3.png" /><figcaption>Fig. 1. Left: a rectangular surface parameterization (RSP) can adapt to long, skinny features, such as the folds of cloth in this bust. Right: such parameterizations lead to high-quality anisotropic quad meshes, even without any optimization or smoothing of mesh elements. Background: mesh quads sorted from most to least rectangular (top left to bottom right). Non-rectangular elements are largely due to irregular vertices, where more (or fewer) than four quads must fit around a vertex.</figcaption></figure>

*Motivation.* Numerous applications benefit from rectangular parameterization. In computational mechanics, it is a valuable starting point for quadrilateral mesh generation, since rectangular elements exhibit *asymptotically* faster convergence for quadrilateral FEM [Arnold et al. 2002] (see Figure 7). In computational architecture, special classes of orthogonal parameterizations correspond to *discrete nets* with essential properties for physical construction [Pottmann et al. 2007]. In digital fabrication, rectangular

<figure><img src="image_4.png" /><figcaption>Fig. 2. A map $\phi$ is *rectangular* if it exhibits no shearing of the coordinate grid (top), but can still apply unequal scaling along the coordinate axes.</figcaption></figure>

ACM Trans. Graph., Vol. 44, No. 4, Article. Publication date: August 2025.
<figure><img src="image_5.png" /><figcaption>Fig. 3. Past methods for making a field integrable can distort the input field, badly violating target properties such as orthogonality (bottom center), or incurring severe scale distortion (bottom left). Our method explores the full space of integrable rectangular fields, yielding a perfect rectangular parameterization with low metric distortion. Here for instance, our method (rectangular surface parameterization) applies just a gentle rotation of the input, yielding less-distorted and better-aligned results than two standard techniques used in nearly all prior work on field-aligned parameterization (curl correction and Hodge decomposition).</figcaption></figure>

networks are used to design piecewise developable surfaces [Rabinovich et al. 2018], optimize principal stress networks [Pellis and Pottmann 2018], and model textiles via Chebyshev nets [Sageman-Furnas et al. 2019]. Finally, in computer graphics, rectangular parameterizations help mimic structural and pigmentation patterns found in nature [Knöppel et al. 2015] (Figure 4); lines of principal curvature form a rectangular network commonly used for illustration and visualization [Hertzmann and Zorin 2000]; and automatic retopology guided by principal directions is commonly used for subdivision modeling [Rouca 2019].

**Integrability of Frame Fields.** Our method for rectangular surface parameterization (RSP) is based on field-guided mapping [Bommes et al. 2013b, Section 3.5], where one first solves for a frame field that encodes the principal directions of the mapping, then recovers the mapping itself. At the heart of RSP is a procedure for making a given frame field integrable, i.e., modifying it so that each frame axis locally corresponds to the gradient of a coordinate function. This seemingly small change has a big impact on the overall pipeline, since a nonintegrable field distorts the parameterization, which degrades mesh quality, and in turn diminishes performance of downstream applications (Section 6.5).

Though great care has been put into designing frame fields [Vaxman et al. 2016], there has been less intense study of how to preserve essential properties (such as orthogonality) when going from the field to a subsequent parameterization. Historically, two strategies are used almost universally to address nonintegrability:

* **Helmholtz-Hodge Decomposition** makes a field X integrable by extracting the curl-free or exact part, i.e., the component Y of X such that $Y = \nabla\phi$ for some scalar function $\phi$ [Kälberer et al. 2007, Section 4]. Note that Hodge decomposition is used implicitly in any method that seeks a “least squares” type solution to the equation $\nabla\phi = X$, since minimizing the residual $\|\nabla\phi - X\|^2$ yields a normal equation $\Delta\phi = \nabla \cdot X$ equivalent to the Poisson equation used in Hodge decomposition. For instance, all methods based on mixed

integer quadrangulation [Bommes et al. 2009, Section 5] effectively adopt the Hodge approach to integrability—even if not cast as a direct modification of the field itself.

* **Curl Correction** seeks a pointwise rescaling s of X such that $curl(sX) = 0$, i.e., such that the rescaled field Y is curl-free [Ray et al. 2006; Zhang et al. 2010]. Since scaling is isotropic, Y describes derivatives of a conformal mapping. Curl correction has been generalized to other settings, such as *polyvector fields* [Diamanti et al. 2015], and is representative of any method that scales field magnitude while keeping directions fixed [Zhang et al. 2010; Simons and Amenta 2024].

An attractive feature of these strategies is that they amount to sparse linear solves. A major drawback, however, is that in both methods the integrable field Y is “anchored” to any flaws in the input field X: Hodge decomposition seeks the closest integrable field Y (in the $L^2$ sense) to the input X; curl correction keeps the directions of X fixed, and is only able to modify magnitude. Hence, even if X exhibits erroneous directions or a poorly-designed global flow, the final field will retain some “memory” of this bad initial guess (Figure 3).

**Integrable Rectangular Fields.** We take a fundamentally different approach: we optimize any user-specified quality criterion (smoothness, alignment, etc.) over the full space of exactly integrable, rectangular frame fields, i.e., fields where axes are orthogonal and positively-oriented, but

need not have the same magnitude. The price we pay is that we must now solve a nonlinear problem rather than a linear one—yet the method remains efficient, and reliably yields optimal results (Figure 21). The benefit is that we often get fields that are better than those produced by linear methods in all respects: they are not only integrable, but are perfectly rectangular by construction, and generally lead to lower-distortion mappings than past methods (Figure 3, right).

<figure><img src="image_6.png" /></figure>
<figure><img src="image_7.png" /><figcaption>Fig. 4. A variety of natural phenomena exhibit rectangular patterns, where a continuously varying aspect ratio is needed to adapt to changes in curvature and pattern orientation. Our rectangular parameterizations can be used to synthesize patterns with a natural orthogonal structure, as seen here in the scales of *scleropages formosus* (Asian arowana).</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 5. The basic steps of our algorithm. The main cost is solving the optimization problem for the integrable rectangular frame. Optional steps needed for mesh extraction are shown in gray.</figcaption></figure>

**Method Overview.** How do we actually compute a rectangular parameterization? One idea is to simply assert that its derivatives must be orthogonal and positively oriented, but directly penalizing the map itself can result in severe violation of injectivity (Figure 10). Instead, we characterize a rectangular parameterization indirectly in terms of the frame field induced by the parameterization (Section 3). In particular, we ask, “if we already had a rectangular parameterization, what must be true about how the corresponding coordinate frame moves around on the surface?” Answering this question leads to our main integrability condition (Equation 10) on data encoding the frame at each point, namely, a rotation and positive scaling of an arbitrary reference frame. In turn, we obtain a positively oriented frame—and hence a locally injective map—by construction.

Algorithmically, we proceed in the opposite direction (Figure 5). First, we compute an orthonormal reference frame $X^0$ via any method. Second, we find rotation angles $\theta$ and positive scale factors $a, b$ that describe an integrable rectangular frame $X$ relative to $X^0$. Third, we solve a Poisson equation to recover the final map $f$ from $X$. Since $X$ is integrable, any deviation from orthogonality is due

<figure><img src="image_9.png" /><figcaption>Fig. 6. General-purpose parameterization methods like ARAP and SLIM minimize overall metric distortion, but exhibit severe shearing. Conversely, conformal methods like BFF give an orthogonal parameterization, but with extreme area distortion. A rectangular parameterization offers the best of both worlds: no shearing, yet far less area distortion than conformal maps.</figcaption></figure>

only to a small amount of discretization error, which vanishes under refinement (Figure 24). The final per-vertex map $f$ can then be used in any downstream application, such as quad meshing.

As with the metric scaling approach to conformal maps [Ben-Chen et al. 2008; Springborn et al. 2008], expressing our parameterization in terms of *log scale factors* (which necessarily induce strictly positive scaling) ensures it is flip-free, modulo discretization error. Overall, starting with a principled discretization of an integrable continuous formulation greatly improves results relative to past methods (Section 6).

*Outline.* After reviewing related work (Section 2), we derive our main integrability condition in the smooth setting (Section 3). We then discretize this condition on triangle meshes (Section 4), and describe how we optimize a frame for integrability, plus additional quality objectives (Section 5). Section 6 numerically validates our approach, via both synthetic metrics and end-to-end application tasks. Section 7 covers limitations and future work.

## 2 Related Work

Here we detail connections with methods from mesh parameterization (Section 2.1), field-guided mapping (Section 2.2), and quad meshing (Section 2.3). Overall, we find no past method able to optimize a given distortion metric directly over the space of rectangular maps. At one end of the spectrum, methods encourage orthogonality via soft penalties which still permit substantial shear. At the other end, conformal maps are exactly orthogonal, but exhibit major scale distortion. Moreover, no past anisotropic quad meshing scheme (to our knowledge) guarantees convergence to perfect rectangular elements under input mesh refinement.

### 2.1 Surface Parameterization

<table><thead><tr><th></th><th colspan="3">shear stretch scale</th></tr></thead><tbody><tr><th>conformal</th><td>X</td><td>X</td><td>✓</td></tr><tr><th>authalic</th><td>✓</td><td>✓</td><td>X</td></tr><tr><th>rectangular</th><td>X</td><td>✓</td><td>✓</td></tr></tbody></table>

Relative to traditional mesh parameterization, we break from the goal of finding a map that is “as isometric as possible,” instead seeking the

least-distorted rectangular map. In other words, unless requested by the user, we penalize neither uniform scale nor axis-aligned
<figure><img src="image_10.png" /><figcaption>Fig. 7. As shown by [Arnold et al. 2002], quadrilateral FEM converges asymptotically faster on rectangular meshes than on general quad meshes. Empirically, meshes produced by RSP achieve the theoretically-predicted optimal convergence rate on common problems, such as solving a Poisson equation (far right). In contrast, past curl correction schemes lead to highly non-rectangular elements, yielding suboptimal convergence.</figcaption></figure>

stretch—only shear is explicitly forbidden. In this way, rectangular parameterization complements area-preserving or *authalic* maps and angle-preserving or *conformal* maps (see inset). In fact, our approach provides a very natural generalization of the *log scale factors* used to describe conformal maps. In particular, we augment the conformal scale factor *u* with an additional authalic scale factor *v* (Section 3.3.1).

The full landscape of parameterization methods is too vast to detail here [Floater and Hormann 2005; Sheffer et al. 2007], but we can get a sense of tradeoffs by comparing to representative methods that penalize non-isometry, such as *ARAP* [Liu et al. 2008] or *SLIM* [Rabinovich et al. 2017], and to conformal methods like *BFF* [Sawhney and Crane 2017]. Here we use the implementations found in *libigl* example 503 [Jacobson et al. 2013]), *Blender* [2025], and by Sawhney and Crane [2019], *resp.* As seen in Figure 6, RSP achieves the best of both worlds: like a conformal map, it preserves right angles between principal axes, yet like non-conformal methods it does not exhibit extreme area distortion.

*Distortion Measures.* More broadly, many distortion metrics have been used for surface parameterization (Rabinovich et al. [2017] provide one list), but most are *isotropic* since they are expressed via scalar invariants of the Jacobian. In contrast, rectangular mapping demands an *anisotropic* distortion measure, where both scalar stretch and vectorial orientation are accounted for. In the field-guided setting, Bommes et al. [2009, Section 5.1] use a penalty that

promotes a prescribed anisotropy along a fixed direction field. Unlike RSP, this energy still permits large scale non-injectivity (Figure 10), and does not jointly optimize for the directions and anisotropy that best minimize distortion. Similarly, Levi [2023b] describe a *shear reduction* penalty that does not yield exact orthogonality (even under refinement), and requires that an injective map is already given as input. In contrast, RSP is orthogonal by construction, and jointly optimizes both stretch and orientation to minimize distortion.

*Local Injectivity.* A variety of recent methods offer hard guarantees on local injectivity [Rabinovich et al. 2017; Gillespie et al. 2021b; Fu et al. 2021]. In the smooth setting, the mapping defined by RSP is guaranteed to be locally injective, since positive scale factors $e^{u+\nu}$, $e^{u-\nu} > 0$ ensure the integrable frame (hence derivatives of the map) cannot invert. In the discrete setting, however, edge vectors mapped into the plane do fit together exactly (Section 4.2.1). Hence, our final Poisson step, which integrates edges in a least-squared sense, can in principle produce flipped triangles—though this rarely occurs in practice. Any such failure to close also vanishes under refinement (Figure 24), and is often mitigated by automatic cone placement, which can help account for a reference frame $X^0$ with poorly chosen global topology (Figure 20). Moreover, nothing prevents the *almost* perfectly integrable rectangular field generated by RSP from being used as input to a locally-injective field-guided parameterization method (e.g., [Myles et al. 2014]), rather than our simple Poisson integration scheme.
### 2.1.1 Conformal Mapping.
Mathematically, RSP is a natural generalization of the *metric scaling* approach to conformal flattening [Ben-Chen et al. 2008; Springborn et al. 2008]. Here, a Riemannian metric *g* is transformed into a flat metric *g*₀ by a positive scaling *g*₀ = e²*u* *g*, where *u*: *M* → ℝ is the *log conformal factor*. This transformation preserves angles between tangent vectors, but uniformly scales area. The factor *u* is obtained by solving the linear *Yamabe equation* [Aubin 2013; Bunin 2008; Sawhney and Crane 2017] or by minimizing a convex energy [Luo 2004; Springborn et al. 2008]. While this approach is efficient, it can induce extreme area distortion, and does not always admit a solution when seeking feature alignment [Myles and Zorin 2013]. RSP overcomes these limitations via an additional log factor *v*: *M* → ℝ that *stretches* the domain while preserving area. Our integrability condition (Equation 10) can be viewed as a generalization of [Myles and Zorin 2013, Proposition 3], which gives the relationship between the log conformal factor *u* and the rotation rate of the frame induced by a conformal mapping.

## 2.2 Field-Guided Parameterization

Past work on field-guided parameterization has focused heavily on *topological* questions, e.g., the conditions under which a field can be parameterized [Shen et al. 2022], or guarantees of local injectivity [Campen et al. 2019; Levi 2023a]. These properties are necessary to ensure a quad mesh can be generated—but do not by themselves guarantee that mesh elements will have good quality.

Geometric criteria are often secondary, and addressed indirectly through weak requirements like bounded distortion [Chien et al. 2016], or soft penalties on metric distortion [Capouellez and Zorin 2024]. Likewise, while orthogonality of the input field promotes orthogonality of the final map, this property is not achieved exactly unless the field is already integrable (Figure 3, center).

### 2.2.1 Integrable Fields.
Past methods for field integrability also treat geometric quality as an “afterthought”—largely prioritizing a topological property (integrability). For example, curl correction achieves integrability through significant isotropic [Ray et al. 2006] or anisotropic scaling [Zhang et al. 2010], but makes no attempt to optimize the direction of anisotropy in order to reduce geometric distortion. Diamanti et al. [2015] and Sageman-Furnas et al. [2019] develop integrability schemes for *polyvector fields*, Pluta et al. [2021] establishes an integrability condition for planar hexagonal meshing, and Coiffier and Corman [2023] use (as we do) moving frames to express integrability constraints for general fields—yet in all these schemes, only soft penalties are used to encourage geometric quality. Our approach to integrability guarantees that, under refinement, the parameterization always approaches a geometrically orthogonal one (Figure 24).

Less attention has been given to integrable *orthogonal fields*.SIMons and Amenta [2024] give conditions under which a frame with *fixed* orientation admits a rectangular parameterization after axis-aligned scaling. However, fixing frame directions can lead to severe scale distortion, as seen in Figure 7, left. Jezdimirović et al. [2023] study fields with fixed singularities and admissible rotations but fail to provide a practical algorithm both constraining orthogonality and minimizing distortion. Finally, Couplet et al. [2024] propose an integrability condition for orthogonal frame fields represented

<figure><img src="image_11.png" /><figcaption>Fig. 8. For field-guided mesh generation, starting with an integrable rectangular field both reduces distortion during the parameterization step, and helps elements adapt to geometric features. Each quad is colored according to the maximum deviation from 90° at any corner, and we give the median/mean/max such deviation across the entire mesh.</figcaption></figure>

with orthogonally decomposable (odeco) tensors [Palmer et al. 2020]. Here, integrability is imposed only weakly, and the method struggles to maintain orthogonality.
<figure><img src="image_12.png" /><figcaption>Fig. 9. Compared to recent learning-based methods, our method also produces meshes that more closely resemble models created by professional artists, with cleaner edge flow, fewer irregular regions, and better-shaped quads. Here for instance we compare against the *Meshtron* method of Hao et al [2024], which is directly trained on artist-generated models.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Fig. 10. Attempting to achieve orthogonality through penalty terms is unsatisfactory, since penalties may “fight” with other energy terms. *Left*: fixing field directions and promoting anisotropy during parameterization leads to large noninjective regions, and many angles far from 90°. *Right*: our method directly explores the space of rectangular maps, yielding injective, orthogonal results even for difficult constraints.</figcaption></figure>

## 2.3 Anisotropic Quad Meshing

Beyond field-guided parameterization, there are a variety of strategies for generating anisotropic rectangular meshes.

In 2D, *boundary layer meshes* concentrate anisotropic elements near the domain boundary to resolve characteristics of, e.g., flow equations, without globally inflating mesh resolution. Here, *paving* or *advancing front* methods incrementally add high-quality elements starting at the boundary, but must introduce highly irregular connectivity to resolve collisions between propagating wavefronts [Merhof et al. 2007]. Likewise, *indirect* methods first triangulate the domain, then merge triangles into quads [Remacle et al. 2013]—again resulting in many singularities and suboptimal, non-rectangular elements (Figure 28, *bottom, center-left*). Explicit control over boundary conditions enables RSP to produce such meshes (Figure 18, *bottom*), while

maintaining a regular, highly rectangular grid structure almost everywhere.

In 3D, early methods obtained anisotropic quads by dicing up surfaces along lines of principal curvature [Alliez et al. 2003; Marinov and Kobbelt 2004]. This local approach makes it difficult to obtain good global structure, resulting in many low-quality, non-quadrilateral elements. Field-guided methods addressed many of these issues, with limitations discussed in Section 2.2. Other methods modify the *metric* of the surface to encourage rectangular elements [Kovacs et al. 2010] or user-controlled anisotropy [Jiang et al. 2015]. Yet metric-based approaches still provide no guarantee that the field will be integrable with respect to the modified metric—which means that any subsequent parameterization can violate criteria like orthogonality.

*Learning-Based Methods.* Finally, some recent methods learn to generate mesh geometry and connectivity based on a given training set [Siddiqui et al. 2024; Son et al. 2024; Shen et al. 2024; Hao et al. 2024]. Here one could train on high-quality rectangular meshes made by artists or engineers. However, current methods leave much to be desired in terms of quality (Figure 9), ability to generalize to unseen object categories, and generation cost (e.g., on the order of minutes–hours for *Meshtron* [Hao et al. 2024], versus seconds–minutes for our method). In this context, RSP provides a way to generate high-quality training data for such networks; our integrability condition (Equation 10) may also prove useful as a regularizer for methods that generate fields [Dielen et al. 2021].

## 3 Smooth Formulation

We develop our method in the setting of continuous differential geometry. An introduction to relevant concepts (including moving frames) can be found in [O'Neill 2006]; further background on moving frames can be found in [Cartan et al. 2001]. Readers uninterested in derivation of the method may wish to skip to Section 5, where we describe the final algorithm.
In the smooth setting, our algorithm consists of a few basic steps:

I. Pick an orthonormal reference frame field $(X_1^0, X_2^0)$.

II. Solve for rotations $\theta$ and per-axis scale factors $a, b$ at each point that transform the reference field into an integrable rectangular field $(aX_1^\theta, bX_2^\theta)$.

III. Use the rectangular field to construct the derivatives $\mu$ of a rectangular parameterization.

IV. Integrate the derivative to obtain the parameterization $f$ itself, by solving a Poisson equation $\Delta f = \nabla \cdot \mu$.

This section is devoted to deriving the equation needed for Step II.

## 3.1 Basic Setup

Consider an oriented smooth surface $M$ with a Riemannian metric $g$ (Figure 11). Let $J_p : T_p M \to T_p M$ denote a counter-clockwise 90° rotation in each tangent space $T_p M$ compatible with the metric, i.e., $g_p(X|_p, JX|_p) = 0$ for all vector fields $X$, at all points $p \in M$. For brevity we drop subscripts $p$ throughout, implying that relationships hold at all points $p \in M$. We use $\Delta$ for the Laplace-Beltrami operator on $(M, g)$, $\nabla \cdot X$ for the divergence operator, and $\nabla_Y X$ for the covariant derivative of $X$ along the direction $Y$.

## 3.2 Moving Frames on Surfaces

A moving frame is any pair of vector fields $X_1, X_2$ on $M$ that are orthonormal with respect to $g$ and consistently oriented, i.e., $X_2 = JX_1$ and $g(X_i, X_j) = \delta_{ij}$, where $\delta_{ij}$ denotes the Kronecker delta. A classic example is the *Darboux frame*, where $X_1, X_2$ are the principal curvature directions—though in our case we do not make any special

<figure><img src="image_14.png" /><figcaption>Fig. 11. Given a surface $M$, we seek a parameterization $f : M \to \mathbb{R}^2$ that takes an orthonormal frame $(X_1, X_2)$ at each point to a rectangular orthogonal frame $(ae_1, be_2)$ in the plane.</figcaption></figure>

<table><caption>Fig. 12. Basic quantities used in our formulation.</caption><thead><tr><th>quantity</th><th>description</th></tr></thead><tbody><tr><td>$M \subset \mathbb{R}^3$</td><td>problem domain (surface)</td></tr><tr><td>$f : M \to \mathbb{R}^2$</td><td>surface parameterization</td></tr><tr><td>$df : TM \to \mathbb{R}^2$</td><td>differential / Jacobian of parameterization</td></tr><tr><td>$a, b : M \to \mathbb{R}_{>0}$</td><td>width, height of infinitesimal rectangles</td></tr><tr><td>$u, v : M \to \mathbb{R}$</td><td>logarithm of $a, b$</td></tr><tr><td>$e_1, e_2 \in \mathbb{R}^2$</td><td>standard basis in the 2D parameter plane</td></tr><tr><td>$X_1, X_2 : M \to TM$</td><td>axes of infinitesimal rectangles (moving frame)</td></tr><tr><td>$\omega_{12} : TM \to \mathbb{R}$</td><td>$\omega_{12}(Y)$ is rotation rate of $(X_1, X_2)$ in Y direction</td></tr><tr><td>$\sigma_1, \sigma_2 : TM \to \mathbb{R}$</td><td>coframe $\sigma_i(X_j) = \delta_{ij}$</td></tr></tbody></table>

choice of frame. We assume only that the frame varies smoothly away from a collection of isolated points, which will correspond to *cone singularities* in the final parameterization. Singularities are further discussed in Section 3.6; for now we consider a frame without singularities.

<figure><img src="image_15.png" /></figure>

**Coframe.** Every moving frame has a corresponding coframe, i.e., a pair of differential 1-forms $\sigma_1, \sigma_2 : TM \to \mathbb{R}$ such that

$$\sigma_i(X_j) = \delta_{ij}, \quad i, j \in \{1, 2\},$$

Intuitively, at each point $p \in M$ the coframe measures the extent of any given vector $Y$ along the frame directions, i.e., it gives the local coordinates of $Y$ in the basis $X_1, X_2$.

**Connection Form.** We also have differential 1-forms $\omega_{ij} : TM \to \mathbb{R}$ for $i, j \in \{1, 2\}$ defined by the relationship

$$\omega_{ij}(Y) := \sigma_i(\nabla_Y X_j),$$

for all tangent vector fields $Y$. Intuitively, $\omega_{ij}(Y)$ describes the change in the $i$th component of the $j$th frame direction as we move along $Y$, and hence expresses change in the frame relative to the moving basis of the frame itself. Since the frame

<figure><img src="image_16.png" /></figure>

is orthonormal, we always have $\omega_{11} = \omega_{22} = 0$, i.e., a unit vector cannot get longer or shorter. Moreover, $\omega_{12} = -\omega_{21}$: from point to point, the frame can only rotate, and $\omega_{12}(Y)$ captures the rate of rotation as we move in the $Y$ direction. We refer to $\omega_{12}$ as the *connection form* for the moving frame $X_1, X_2$, since it effectively tells us how local coordinate systems at nearby points “connect.”

**Structure Equations.** One can show that any $\sigma$ and $\omega$ that come from a frame $X$ must satisfy *Cartan's first structure equation*

$$\begin{aligned} d\sigma_1 &= -\omega_{12} \wedge \sigma_2, \\ d\sigma_2 &= -\omega_{21} \wedge \sigma_1, \end{aligned} \quad (1)$$

where $d$ and $\wedge$ are the exterior derivative and wedge product, resp. In fact, for a given $\sigma$, $\omega$ is the unique 1-form satisfying this equation [O'Neill 2014, Corollary 1.8.5]. Note that for our algorithm we will not need *Cartan's second structure equation*, which in 2D simply asserts that $\omega$ is closed ($d\omega = 0$).

## 3.3 Rectangular Parameterization

We now characterize rectangular parameterizations in terms of moving frames. Suppose in particular we want a parameterization $f : M \to \mathbb{R}^2$ that locally looks like a pure nonuniform scaling along some orthonormal frame $X$, i.e., no shearing. In particular, let $e_1 = (1, 0)$, $e_2 = (0, 1)$ be the standard basis for $\mathbb{R}^2$. Then we want our parameterization and moving frame to satisfy

$$\begin{aligned} df(X_1) &= ae_1, \\ df(X_2) &= be_2, \end{aligned} \quad (2)$$

where the differential of the parameterization $df : TM \to \mathbb{R}^2$ maps tangent vectors $X$ on the surface to the corresponding vectors in
<figure><img src="image_17.png" /><figcaption>Fig. 13. Any rectangular parameterization $f: M \to \mathbb{R}^2$ can be associated with an orthonormal frame $(X_1, X_2)$ and a pair of scale factors $a, b$ along the two axes. The value $\omega_{12}(Y)$ encodes the speed of rotation in any given direction $Y$, which can be thought of as the inner product of $Y$ with the representative vector for the connection form $\omega_{12}$. Here for instance the frame rotates in the $\theta$ direction, but remains parallel in the $r$ direction.</figcaption></figure>

the plane $\mathbb{R}^2$, and $a, b: M \to \mathbb{R}_{>0}$ are strictly positive functions that describe the scaling along each frame direction (Figure 11).

### 3.3.1 Log Scale Factors.
Alternatively, we can express our parameterization as

$$ 
\begin{aligned}
df(X_1) &= e^{u+\upsilon} e_1, \\
df(X_2) &= e^{u-\upsilon} e_2,
\end{aligned} 
\qquad (3) 
$$

for arbitrary (i.e., not necessarily positive) functions $u, v: M \to \mathbb{R}$, which we call the (log) conformal factor and (log) authalic factor, resp. In particular, when $v = 0$ we get $a = b = e^u$, i.e., the width and height of infinitesimal rectangles are scaled by the same amount, preserving angles. When $u = 0$ we get $a = e^v = 1/b$, i.e., the width and height are scaled by reciprocal amounts, preserving local area.

Log factors are natural for several reasons. First, exponentiation guarantees that scale factors $a, b$ are always nonnegative, avoiding overlap (and helping to prevent inverted mesh elements). Second, they yield an integrability condition linear in $u, v$, as discussed below. Finally, they make it straightforward to express objective functions and boundary conditions in simple geometric terms (e.g., angle or area preservation); see Sections 5.1.1 and 5.1.3. This setup directly generalizes the standard formulation of conformal maps in terms of the log conformal factor (Section 3.4.2).

## 3.4 Integrability Condition

We next derive an expression for the connection form $\omega_{12}$ associated with a rectangular parameterization $f$, i.e., for the rotation speed

of a rectangular frame, in terms of the log scale factors $u, v$. This relationship serves as an intermediate form of our final integrability condition, which will ultimately be expressed in terms of a rotation and scaling of a fixed reference frame (Section 3.4.1).

First note that the differential $df$ can be expressed via the coframe:

$$ 
\begin{aligned}
df = a\sigma_1 e_1 + b\sigma_2 e_2.
\end{aligned} 
\qquad (4) 
$$

Applying the exterior derivative to both sides of this equation (and noting that $e_1, e_2$ are constant) then yields

$$ 
\begin{aligned}
0 = ddf = (da \wedge \sigma_1 + ad\sigma_1)e_1 + (db \wedge \sigma_2 + bd\sigma_2)e_2.
\end{aligned} 
\qquad (5) 
$$

Since $e_1, e_2$ are linearly independent, we have

$$ da \wedge \sigma_1 + ad\sigma_1 = 0, $$

$$ db \wedge \sigma_2 + bd\sigma_2 = 0. $$

Applying the first structure equation (Equation 1) and making the substitutions $da = ad(u + v)$, $db = bd(u - v)$ yields

$$ 
\begin{aligned}
d(u + v) \wedge \sigma_1 &= \omega_{12} \wedge \sigma_2, \\
d(u - v) \wedge \sigma_2 &= \omega_{21} \wedge \sigma_1.
\end{aligned} 
\qquad (6) 
$$

Applying the 2-forms on both sides of these equations to the vectors $X_1, X_2$ (via the usual identity $\alpha \wedge \beta(X, Y) = \alpha(X)\beta(Y) - \alpha(Y)\beta(X)$), and recalling that $\sigma_i(X_j) = \delta_{ij}$, we get

$$ 
\begin{aligned}
d(u + v)(X_2) &= -\omega_{12}(X_1), \\
d(u - v)(X_1) &= -\omega_{21}(X_2).
\end{aligned} 
$$

Noting that $\omega_{12}$ can be expanded in the coframe as $\omega_{12} = \omega_{12}(X_1)\sigma_1 + \omega_{12}(X_2)\sigma_2$ (and recalling that $\omega_{12} = -\omega_{21}$), we then get an equation for the rotation speed of the rectangular frame $X_1, X_2$ in terms of the log scale factors $u$ and $v$:

$$ \omega_{12} = -D_{X_2}(u + v)\sigma_1 + D_{X_1}(u - v)\sigma_2. \qquad (7) $$

Here $D_X\phi$ denotes the directional derivative of a scalar function $\phi$ along the vector field $X$.

### Geometric Interpretation.
Geometrically, this equation tells us that in order for scale factors to vary over space, an infinitesimal rectangle must exhibit some amount of “bending” (see Figure 13 for an illustration). E.g., when $M \subset \mathbb{R}^2$, $\omega_{12} = 0$ only when the horizontal scaling $a$ is purely a function of $x$ and the vertical scaling $b$ is purely a function of $y$. In contrast, the inset example $f(x, y) := (\sqrt{x}, \frac{1}{2} - \frac{1}{2}\cos(\pi y))$ on the domain $M = [0, 1]^2$ exhibits independent scaling along each axis, and the frame does not rotate.

### 3.4.1 Final Integrability Condition.
We can express the frame $X$ in Equation 7 as a rotation of some fixed reference frame $X^0$ by a continuously-varying angle $\theta: M \to \mathbb{R}$:

$$ 
\begin{aligned}
X_1^\theta &:= \cos(\theta)X_1^0 + \sin(\theta)X_2^0, \\
X_2^\theta &:= -\sin(\theta)X_1^0 + \cos(\theta)X_2^0.
\end{aligned} 
\qquad (8) 
$$

Since $\theta$ is continuous, $X^0$ will have the same singularities as $X^0$ (see Section 3.6.2 for further discussion). Likewise, if $\sigma^0$ is the coframe of $X^0$, then

$$ 
\begin{aligned}
\sigma_1^\theta &:= \cos(\theta)\sigma_1^0 + \sin(\theta)\sigma_2^0, \\
\sigma_2^\theta &:= -\sin(\theta)\sigma_1^0 + \cos(\theta)\sigma_2^0.
\end{aligned} 
\qquad (9) 
$$
We then have the following relationship between the rotation speed of the original and rotated frames:

**LEMMA 3.1.** *The connection form $\omega_{12}^\theta$ of the rotated frame differs from the connection form $\omega_{12}^0$ of the reference frame by the differential $d\theta$ of the angle of rotation, i.e.,

$$ \omega_{12}^{\theta} = \omega_{12}^{0} - d\theta. $$

For example, for a constant reference frame on a flat region $M \subset \mathbb{R}^2$, we get just $\omega_{12}^\theta(Y) = d\theta(Y)$, i.e., the rate of rotation is just the change in angle as we move along any given direction $Y$. See Appendix A for a proof.

Substituting this expression for $\omega_{12}$ into Equation 7 yields our final integrability condition, which says how the rotation angles $\theta$ and log scale factors $u, v$ must be related in a rectangular map:

$$ \omega_{12}^{0} - d\theta = -D_{X_2}^{\theta}(u + v)\sigma_{1}^{\theta} + D_{X_1}^{\theta}(u - v)\sigma_{2}^{\theta}. \quad (10) $$

Note that since the frame $X^\theta$ and coframe $\sigma^\theta$ depend on $\theta$, this equation is linear in $u, v$ but nonlinear in $\theta$.

**3.4.2 Relationship to Conformal Maps.** A conformal map exhibits scaling but not shearing—in terms of scale factors, conformal maps correspond to the case $a = b$, or equivalently $v = 0$ (Equation 3). If we also assume a constant reference frame $X^0$ (so that $\omega_{12}^0 = 0$), Equation 10 reduces to just

$$ d\theta = (D_{X_2}u)\sigma_1 - (D_{X_1}u)\sigma_2, $$

or equivalently,

$$ d\theta = *du, \quad (11) $$

where * denotes the Hodge star (since for any 1-form $\alpha$ on a two-dimensional domain, $*\alpha(X) = -\alpha(JX)$ for all $X$). As discussed in Myles and Zorin [2013, Proposition 3], this equation gives the relationship between the scale factor $u$ and the connection 1-form $\omega_{12} = - * d\theta$ in a conformal parameterization.

## 3.5 Poisson Equation

Solving Equation 10 for $\theta, u, v$ provides data sufficient to determine the differential $df$ of a rectangular parameterization, which we can then integrate to obtain the parameterization $f$ itself. In particular, we know from Equation 4 that we can express the differential as

$$ \mu := a\sigma_1^\theta e_1 + b\sigma_2^\theta e_2, \quad (12) $$

i.e., as a 1-form that sends $X_1^\theta$ to $ae_1$ and $X_2^\theta$ to $be_2$. In the smooth setting the equation $df = \mu$ is exactly solvable, and we could in principle integrate $\mu$ directly to obtain $f$. In the discrete setting, however, there will be a small amount of nonintegrability arising from discretization error (i.e., the fact that we use a mesh with finitely many elements). Instead, as is standard practice in geometry processing, we seek the parameterization $f$ that minimizes the residual $\|df - \mu\|_{L^2}$, leading to a vector-valued Poisson equation

$$ \Delta f = \nabla \cdot \mu. \quad (13) $$

Solving this equation yields our final parameterization $f: M \to \mathbb{R}^2$.

**3.5.1 Local Injectivity.** Since (i) $a = e^{u+\nu}$ and $b = e^{u-\nu}$ are strictly positive, (ii) $(X_1^\theta, X_2^\theta)$ is a rotation of a positively-oriented frame, and (iii) $\mu$ is exactly integrable, the map $f$ defined by our procedure has a nondegenerate differential (i.e., Jacobian) by construction. Hence, in the smooth setting, it is always locally injective. In the discrete setting we no longer have a hard guarantee, due to discretization error. However, since this error goes to zero under refinement (Figure 24), we can reduce any overlap by simply refining the input mesh. In contrast, methods based on minimizing, e.g., an elastic [Liu et al. 2008] or other isometric energy [Rabinovich et al. 2017] have no such guarantee in the smooth setting—putting a far greater burden on discrete optimization to provide injectivity.

## 3.6 Cone Singularities

An important observation first made by Kharevych et al. [2006] in the context of conformal parameterization is that carefully-placed cone singularities can significantly reduce distortion. A complementary perspective is that such cones correspond to singularities in the frame field used for field-guided parameterization.

So far, we have considered a frame field $X$ that varies smoothly across the entire surface. Suppose instead that the field can discontinuously jump in angle by multiples of $\pi/2$ along a network of curves $\Gamma \subset M$. To find a parameterization $f$ aligned with $X$, we can now imagine cutting $M$ along $\Gamma$ before mapping it into the plane (Figure 14). If we pull back the metric of the plane under $f$, we get a metric on $M$ that is intrinsically flat, except at vertices $p \in \Gamma$ that map to a corner $f(p)$ of angle different from $2\pi$ (Figure 14, bottom right). Gluing such a corner back together along cut edges forms a cone—hence the name cone singularity.

In terms of our rectangular parameterization, the only thing we need to consider is the behavior of the log scale factors $u$ and $v$ as we cross the cut network $\Gamma$. From Equation 4 and the definition of $a, b$ and $\sigma_1, \sigma_2$, we have

$$ df(Y) = \langle e^{u+\nu}X_1, Y \rangle e_1 + \langle e^{u-\nu}X_2, Y \rangle e_2, $$

for any $Y$. Suppose that $X$ experiences a counter-clockwise 90-degree rotation as we cross $\Gamma$, i.e., $\tilde{X}_1 = X_2$ and $\tilde{X}_2 = -X_1$. Then on the other side of the cut we have

$$ d\tilde{f}(Y) = \langle e^{\tilde{u}+\tilde{v}}X_2, Y \rangle e_1 - \langle e^{\tilde{u}-\tilde{v}}X_1, Y \rangle e_2. $$

<figure><img src="image_18.png" /><figcaption>Fig. 14. Parametric distortion can be reduced by allowing cone singularities in the parameterization, corresponding to points $p$ around which the total rotation of the frame is nonzero (e.g., $\pm 90^\circ$). To do so, we cut the surface along a graph $\Gamma$ that passes through all singular points and turns the surface $M$ into a topological disk.</figcaption></figure>
Since $df(Y)$ and $d\tilde{f}(Y)$ are also related by a 90° rotation, we get

$$u = \tilde{u} \quad \text{and} \quad v = -\tilde{v}. \qquad (14)$$

I.e., conformal scaling (which is the same in all directions) is unchanged, whereas authalic scaling (which has reciprocal magnitude for $X_1$ and $X_2$) is applied to the opposite axis. A clockwise 90° rotation yields the same relationship. If the frame experiences a 180° rotation, we instead have $\tilde{u} = u$ and $\tilde{v} = v$, i.e., both scale factors are continuous at the cut. These relationships become linear conditions in our discrete formulation (Section 4).

Note that only the choice of singularities, and not the particular cut $\Gamma$, affects the quality of the final parameterization.

### 3.6.1 Singularity Index.
We can assign each point $p \in \Gamma$ an index, equal to the number of times the frame rotates as we walk around $p$. More precisely, recall that $\omega_{12}(Y)$ gives the angular velocity of the frame in the direction $Y$. Hence, if $\gamma_p$ is a small loop around $p$, then the index of $p$ is

$$\mathrm{ind}(p) := \frac{1}{2\pi} \int_{\gamma_p} \omega_{12}^0.$$

Since we allow only 90° rotations, all singularities in our theory will have index $k/4$, $k \in \mathbb{Z}$, though one typically wants only singularities of index $\pm 1/4$ to avoid significant distortion (corresponding to irregular vertices of degree beyond 3 or 5 in the resulting quad mesh).

### 3.6.2 Preservation of Singularities.
Recall that the frame $X^\theta$ used for rectangular parameterization differs from the reference frame $X^0$ by an angle $\theta$ at each point. Since $\theta$ varies continuously, the rotated field will (by construction) exhibit the same singularities as the reference field. In particular, applying Stokes' theorem, the new index around any point $p$ is given by

$$\frac{1}{2\pi} \int_{\gamma_p} (\omega_{12}^0 - d\theta) = \mathrm{ind}(p) - \frac{1}{2\pi} \int_{\partial\gamma_p} \theta = \mathrm{ind}(p).$$

<figure><img src="image_19.png" /><figcaption>Fig. 15. We adopt a cone placement strategy that does not penalize rectangular stretching, yet still yields lower area distortion as the number of cones increases (left to right). Even with many cones, the discrete parameterization remains almost perfectly rectangular. Inset histograms tabulate the angle between columns of the Jacobian of $f^{-1}$ across triangles.</figcaption></figure>

## 3.7 Discussion

Note that, due to Equation 5, our formulation implicitly relies on Poincaré's lemma, i.e., that on a simply connected domain, every closed differential form is also exact. However, our final algorithm still applies to nonsimply connected surfaces, since we effectively cut our surface into a disk (and apply boundary conditions that ensure seamless transitions across the cut); see Section 3.6.

It is also natural to wonder whether other, “simpler” equations established earlier in our derivation might provide a suitable starting point for rectangular parameterization. For instance, one could attempt to discretize Equation 2 and solve simultaneously for $f$, $X_1$, $X_2$, $a$, and $b$. However, it is unclear in this setting how to incorporate singularities, due to discontinuities in $X_1$, $X_2$ as well as $a$, $b$; moreover, one must still solve a system of nonlinear equations. Likewise, rather than work with a reference frame, one could try to solve Equation 7 directly, but would have to then incorporate the structure equations into the constraints, and would give no direct way to prescribe singularities. Our formulation, in terms of just an unknown rotation $\theta$ and log scale factors $u, v$, is much closer in spirit to classical conformal mapping, easily incorporates singularities, and leads to a robust discretization, as explored in the next section.

## 4 Discretization

We next describe how to discretize the two main equations needed for our algorithm: the integrability condition (Equation 10) and the Poisson equation (Equation 13). We consider a manifold oriented triangle mesh $M = (V, E, F)$, using indices $i \in V$ for vertices, pairs of indices $ij \in E$ for edges, and triples of indices $ijk \in F$ for triangles. We let $H$ be the set of oriented *halfedges*, using $\overline{ij} \neq \overline{ji}$ to denote two opposite orientations, and let $C$ be the set of all $3|F|$ triangle corners, using $\phi_i^{jk}$ to denote a quantity $\phi$ at corner $i$ of triangle $ijk$.

<figure><img src="image_20.png" /><figcaption>Fig. 16. The discrete quantities used in our formulation.</figcaption></figure>
Discrete quantities are illustrated in Figure 16. We assume triangles have nonzero corner angles $\alpha_i^{jk} \in (0, \pi)$ and areas $A_{ijk} \in \mathbb{R}^{>0}$, and use vertex areas $A_i := \frac{1}{3} \sum_{ijk} A_{ijk}$. We let $w_{ij}^k := \frac{1}{2} \cot \alpha_k^{ij}$ denote the cotan weight for halfedge $ij$ in triangle $ijk$ [MacNeal 1949, Equation 3.17]; the weight $w_{ij}$ for edge $ij$ is then the sum of cotan weights over its halvedges (two on the interior; one on the boundary). Nominally, we encode log scale factors as values $u_i, v_i$ at each vertex $i \in V$, though to handle singularities we more generally store values $v_i^{jk}$ at corners (Supplement, Appendix E).

The reference frame is encoded as a pair of orthonormal vectors $(X_1^0)_{ijk}$, $(X_2^0)_{ijk} \in \mathbb{R}^2$ relative to a local coordinate system in each triangle $ijk$, which by convention we align with the first edge $ij$. We use $\varphi_{ijk}$ to denote the angle of $(X_1^0)_{ijk}$ in this coordinate system; we use $\varphi_{ab}$ to denote the angle of this same vector relative to the direction of any of the three halvedges $\overrightarrow{ab} \in \{\overrightarrow{ij}, \overrightarrow{jk}, \overrightarrow{ki}\}$. We likewise use $\theta_{ijk}$ to denote the rotation of the current frame $X^\theta$ relative to the reference frame $X^0$. Hence, the angle of the rotated frame direction $X_1^\theta$ is given by $\eta_{ijk}^\theta := \varphi_{ijk} + \theta_{ijk}$. We can parallel transport any vector $Y_{ijk}$ to a neighboring triangle $jil$ to get the vector

$$Y_{ijk} \to Y_{ijk} := R_{\varphi_{ji}} - \varphi_{ij} + \pi Y_{ijk}, \quad (15)$$

where $R_\varphi$ denotes a counter-clockwise rotation by $\varphi$. Finally, the 1-form describing the rotation of the reference frame is encoded as a dual discrete 1-form, i.e., a value $\omega_{ij}^0$ at each (dual) edge that gives the counter-clockwise angle of rotation from $X_{ijk}^0$ to $X_{jil}^0$, where $ijk, jil$ are the two triangles containing edge $ij$.

## 4.1 Discretized Integrability Condition

<figure><img src="image_21.png" /></figure>

To discretize Equation 10, we integrate it along the dual edge $ab$ associated with each edge $ij \in E$. Throughout we assume each interior edge $ij$ is contained in two triangles $ijk, jil$, and let $a$ and $b$ be the circumcenters of $ijk$ and $jil$ (resp.), which we connect to the edge midpoint $m$ to form the dual edge. We use $T$ to denote the intrinsic unit vector tangent to dual edge $ab$. The left-hand side of Equation 10 is discretized via discrete exterior calculus [Desbrun et al. 2006], yielding

$$\omega_{ij}^0 - (\theta_{ijk} - \theta_{jil}). \quad (16)$$

To discretize the right-hand side, consider any facewise affine function interpolating values $\phi_i$ at vertices. Integrating over the dual edge $ab$, we have

$$\int_{ab} (D_{X_1^\theta} \phi) \sigma_2^\theta = \int_{ab} \langle X_1^\theta, \nabla \phi \rangle \langle X_2^\theta, T \rangle \, ds,$$

and similarly for $(D_{X_2^\theta} \phi) \sigma_1^\theta$. To obtain a more explicit formula, we perform integration piecewise over the two segments $am$ and $mb$. We then evaluate these expressions for $\phi = u - v$ and $\phi = u + v$, and add these terms to Equation 16. After some simplification, the contribution of triangle $ijk$ to the right-hand side is

$$\rho_{ij}^k(u, v, \theta) := \frac{1}{2} \cot \alpha_k^{ij} \left( (u_j - u_i) + \cos(2\eta_{ijk}(\theta))(v_j - v_i) + \sin(2\eta_{ijk}(\theta))(\cot \alpha_j^k(v_k - v_i) + \cot \alpha_i^k(v_k - v_j)) \right).$$

For any boundary edge $ij$, we let $\rho_{ji}^l = 0$. Also, since this formula assumes the reference frame is expressed relative to edge $ij$, we must re-express $\eta$ relative to edges $jk$ and $ki$ when evaluating terms for those edges. In particular, if $\eta_{ij} := \varphi_{ijk} + \theta_{ijk}$, then

$$\eta_{jk} := \eta_{ij} - (\pi - \alpha_j^k) \quad \text{and} \quad \eta_{ki} := \eta_{jk} - (\pi - \alpha_i^k).$$

To get the full residual corresponding to edge $ij$, we add the corresponding term $\rho_{ji}^l$ from the neighboring triangle $jil$, and subtract the discrete left-hand side, yielding our final discrete condition

$$F_{ij}(u, v, \theta) := \rho_{ij}^k(u, v, \theta) + \rho_{ji}^l(u, v, \theta) + (\theta_{ijk} - \theta_{jil}) - (\omega_{12}^0)_{ij} = 0, \quad (17)$$

which must hold for each edge $ij \in E$. Assembly of the overall residual and its Jacobian is detailed in Appendix E of the supplement.

When $v = 0$ (pure conformal scaling), Equation 17 reduces to a standard (cotan) discretization of Equation 11. Note also that Equation 17 is invariant to rotations $\eta \mapsto \eta + \pi$, correctly reflecting the fact that we only care about how much scaling there is along the horizontal/vertical axis, and not the particular direction.

## 4.2 Discretized Poisson Equation

We discretize the Laplace and divergence operators in Equation 13 using the usual cotan formula (see Appendix B of the supplement).

<figure><img src="image_22.png" /></figure>

### 4.2.1 Discrete Differential.
We discretize the term $\mu$ in Equation 13 by integrating the right-hand side of Equation 12 along each edge. Geometrically, $\mu_{ij} \in \mathbb{R}^2$ gives the vector along edge $ij$ in the parameter domain. Ideally, we could sum these vectors incrementally to get final coordinates $f_i$. Due to discretization error, however, they do not fit together perfectly—we instead solve a discrete Poisson equation (Supplement, Appendix B) that effectively finds a least-squares solution.

More explicitly, for each edge $ij$ we compute average scale factors

$$a_{ij} := e^{(u_i + u_j + v_i + v_j)/2}, \quad b_{ij} := e^{(u_i + u_j - v_i - v_j)/2},$$

which for $v = 0$ agree with the edge length ratios appearing in the theory of discrete conformal equivalence [Springborn et al. 2008, Equation 2]. Likewise, since $\sigma_i(Z) = \langle X_i^\theta, Z \rangle$, we take the average inner product of the vector $Z_{ij}$ along each edge with the frame directions $X_i^\theta$ in the two neighboring triangles—being careful to account for different coordinate systems. Overall, each triangle $ijk$ makes the following contributions to equations along its three edges:

$$\begin{aligned}
\mu_{ij}^k &= \ell_{ij} \begin{bmatrix} +a_{ij} \cos(\eta_{ijk}) \\ +b_{ij} \sin(\eta_{ijk}) \end{bmatrix}, \\
\mu_{jk}^l &= \ell_{jk} \begin{bmatrix} -a_{jk} \cos(\eta_{ijk} + \alpha_j^k) \\ -b_{jk} \sin(\eta_{ijk} + \alpha_j^k) \end{bmatrix}, \\
\mu_{ki}^j &= \ell_{ki} \begin{bmatrix} -a_{ki} \cos(\eta_{ijk} - \alpha_i^k) \\ -b_{ki} \sin(\eta_{ijk} - \alpha_i^k) \end{bmatrix}.
\end{aligned}$$

For each edge $ij$ contained in a pair of triangles $ijk, jil$ we then let

$$\mu_{ij} := \mu_{ij}^k + \mu_{ji}^l. \quad (18)$$

For edges $ij$ on the domain boundary, we let $\mu_{ji} = \mu_{ij}$, i.e., we effectively just use the frame from the single triangle containing $ij$.
<figure><img src="image_23.png" /><figcaption>Fig. 17. Since our fields are rectangular by construction, the shear-free condition does not “fight” with other objectives. We are therefore free to minimize any user-defined distortion metric $\Phi$ over the space of rectangular parameterizations. Here we minimize angle, metric, and area distortion (left to right).</figcaption></figure>

# 5 Algorithm

Our final algorithm is comprised of three basic steps:

(1) Compute a reference frame (or cross field) $X^0$, using any algorithm. Section 6.1 details algorithms used in this paper.

(2) Solve for scale factors $u, v$ and rotations $\theta$ that transform $X^0$ into an integrable rectangular field (Section 5.1).

(3) Recover the final parameterization from this field (Section 5.2).

A diagram of these steps is shown in Figure 5; detailed pseudocode can be found in the supplement (Appendix E).

## 5.1 Optimization

Given a reference frame $X^0$, we use constrained optimization to find log scale factors $u, v$ and angles $\theta$ that satisfy our integrability condition (Equation 10) while minimizing an objective $\Phi(u, v, \theta)$ that encodes quality criteria such as angle or area preservation (Section 5.1.1 gives several examples). Letting $F$ be our discretized integrability equation from Equation 17, we solve the problem

$$ \min_{u, v:V \to \mathbb{R} \atop \theta:F \to \mathbb{R}} \Phi(u, v, \theta) \quad \text{s.t.} \quad F(u, v, \theta) = 0. \quad (19) $$

### 5.1.1 Objective Functions.
The constraint $F = 0$ restricts our search to shear-free rectangular maps; we can then optimize any other objective $\Phi$ over this space. For instance:

* $\Phi_{\text{iso}}(u, v, \theta) := \sum_{i \in V} A_i(u_i^2 + v_i^2)$ — tries to make parameterization as isometric as possible

* $\Phi_{\text{area}}(u, v, \theta) := \sum_{i \in V} A_i u_i^2$ — encourages area preservation by minimizing scale distortion

* $\Phi_{\text{angle}}(u, v, \theta) := \sum_{i \in V} A_i v_i^2$ — encourages angle preservation by minimizing stretching

* $\Phi_{\text{align}}(u, v, \theta) := \sum_{ijk \in F} A_{ijk} \theta_{ijk}^2$ — encourages alignment with the given field $X^0$

Figure 17 demonstrates the effect of several different objectives. Since $u$ and $v$ are logarithmic scale factors, we obtain an infinite energy barrier to degenerate maps—further encouraging local in- jectivity beyond the reasons already discussed in Section 3.5.1. One can of course use other objectives, e.g., Section 6.6 describes how to find a parameterization close to a Chebyshev net.

### 5.1.2 Regularizer.
Equation 10 links the rotation speed of $X^0$ with derivatives of the log scale factors $u$ and $v$. In theory, the frame can rotate arbitrarily fast; in practice, however, rotation speed must be limited. For instance, our discretization does not distinguish between a frame that makes a full turn along a dual edge, and one that is not turning at all. This ambiguity introduces numerical instability, preventing convergence of the optimizer. Instead of constraining rotation along dual edges, we regularize variation in the authalic factor $v$ by adding a discrete Dirichlet energy to our objective:

$$ \Phi_{\text{reg}}(u, v, \theta) := \varepsilon v^\top L v. $$

As $\varepsilon \to \infty$ the authalic factor is constant, and Equation 10 reduces to the well-behaved linear integrability condition for conformal maps (Equation 11). By default we use $\varepsilon = 10^{-2}$, except in Section 6.6 where more variability is needed and we set $\varepsilon = 10^{-4}$.

### 5.1.3 Boundary Conditions.
Rectangular parameterizations accom- modate a variety of boundary conditions, naturally expressed as linear constraints on the log scale factors $u, v$ and frame rotations $\theta$.

**Feature Alignment.** To align the parameterization to boundary or feature curves (as in Figure 8), we simply constrain $\theta$ in triangles adjacent to feature edges, and omit integrability constraints $F_{ij} = 0$ for these edges. In our case, we start with a feature-aligned reference field $X^0$ and set $\theta = 0$.

**Element Shape.** We can also prescribe the shape of rectangles— and ultimately, mesh elements—near the boundary (Figure 18). Note
<figure><img src="image_24.png" /><figcaption>Fig. 18. Using scale factors as primary variables provides direct control over element sizing near the boundary. Here we prescribe values of u and v along the inner boundary, leaving the outer boundary free. Automatic cone singularities help adapt the rest of the field to these boundary conditions.</figcaption></figure>

that ab and a/b quantify local area and aspect ratio, resp., and that

$$u = \frac{1}{2} \log(ab) \quad \text{and} \quad v = \frac{1}{2} \log(a/b).$$

Assuming $X_1^0$ is tangent to the boundary, we then have:

* $u|_{\partial M} = p$ — prescribed area

* $v|_{\partial M} = p$ — prescribed aspect ratio

* $(u - v)|_{\partial M} = p$ — prescribed normal length

* $(u + v)|_{\partial M} = p$ — prescribed tangent length

for any prescribed function $p$ (defined at boundary vertices). As with conformal maps, not all boundary conditions are feasible [Sawhney and Crane 2017, Section 3.2.1], and we do not have a clear characterization of compatibility conditions—in practice the solver may simply fail if the boundary is overconstrained. Note also that extreme boundary conditions can induce distortion (Supplement, Figure 32), which can be mitigated via cone singularities (Figure 18).

### 5.1.4 Solver.
To solve Equation 19, we apply a constrained Newton method to the first order optimality conditions (see Supplement, Appendix C). In practice we always converge to a valid solution, despite the nonlinearity of this problem. The initial guess for $u, v$, and $\theta$ seems to have no impact on results (Section 6.2). We hence initialize with $u = v = \theta = 0$, i.e., no scaling or rotation of $X^0$.

## 5.2 Recovering Parameterization

After optimization, we have log scale factors $u, v$ and frame rotations $\theta$ that describe some rectangular parameterization $f$. However, we still need to recover $f$ itself, i.e., we need to solve a discrete Poisson equation (corresponding to Equation 13) for the unknown vertex coordinates $f: V \to \mathbb{R}^2$. This equation is solved independently for each scalar coordinate function, possibly after cutting the surface into a topological disk. As in many parameterization methods, the main implementation challenge is handling cuts through cones. See Supplement, Appendix Band D for further details.

## 5.3 Cone Singularities

To incorporate singularities, we start with a cross field [Vaxman et al. 2016, Section 2], and extract a frame field $X^0$ with $\pm 90^\circ$ or $180^\circ$ jumps. Since the initial field is free to rotate, only its singularities— and not the field quality—will have an impact on the solution (unless we explicitly encourage alignment with $X^0$).

### 5.3.1 Jump Curve and Sign Bit.
We first identify edges where $X^0$ is discontinuous, corresponding to the cut network $\Gamma$ from Section 3.6. In accordance with Equation 14, we also track the change in the sign of $v$ as we cross each edge. We mark edges in the cut network via a binary labeling $\Gamma: E \to \{0, 1\}$, and encode flips in the sign of $v$ via values $s: E \to \{-1, +1\}$.

To determine these values, we initialize $\Gamma_{ij} \leftarrow 0$ and $s_{ij} \leftarrow 1$ on all edges. We then perform breadth-first traversal: starting in any triangle $ijk$, we pick any of the four frames $X_{ijk}^0$ aligned with the cross field, and transport it to each neighboring triangle $jil$, via Equation 15. If $jil$ has not been visited, we set $X_{jil}^0$ to the cross direction closest to the transported frame $X_{ijk \to jil}^0$. Otherwise, if it has been visited, we let $\beta_{ij} \in [-\pi, \pi)$ be the angle from the transported frame to the existing value of $X_{jil}^0$. If $X_{jil}^0$ is not the closest of the four possible rotations (i.e., if $\beta_{ij}$ is greater than $\pm \pi/4$), there is a jump, and we set $\Gamma_{ij}$ to 1. Moreover, if the closest frame differs by a quarter rotation (i.e., if $\beta \in [\pi/4, 3\pi/4)$ or $\beta \in [-3\pi/4, -\pi/4)$), then $v$ requires a sign flip, and we set $s_{ij}$ to -1. Traversal continues until all faces have been visited.

Per-edge sign bits $s_{ij}$ in turn determine the sign of $v_i^{jk}$ at each triangle corner, relative to a single value $v_i$ stored at vertex $i \in V$. In particular, let $j_0, \dots, j_m$ be the neighbors of $i$ in counter-clockwise order. Then the value of $v$ at any corner around vertex $i$ is

<figure><img src="image_25.png" /></figure>

$v_i^{jn_{j+1}} := s_i^{jn_{j+1}} v_i$, where the sign $s$ for each corner is given by

$$v_i^{jn_{j+1}} := s_i^{jn_{j+1}} v_i,$$

$$s_i^{jn_{j+1}} := \prod_{p=0}^{n-1} s_{ij_{p+1}}. \qquad (20)$$

### 5.3.2 Automatic Singularities.
So far we have assumed cone singularities are fixed ahead of time. However, we can adapt the method of Myles and Zorin [2013] to find singularities that help minimize the objective $\Phi$. We start by cutting the mesh into a disk along edges (as
<figure><img src="image_26.png" /><figcaption>Fig. 19. Our greedy cone placement strategy for rectangular parameterization works better than even optimal cone placement strategies for conformal maps. Here we compare with the method of Li et al. [2022], yielding lower total root mean square stretch S [Sander et al. 2001, Section 3]. Colors plot local scaling $\sqrt{\zeta_1\zeta_2}$ and shearing $\zeta_1/\zeta_2$ on a log scale, where $\zeta_1, \zeta_2$ are singular values of the affine map from 3D to 2D triangles.</figcaption></figure>

<figure><img src="image_27.png" /><figcaption>Fig. 20. As with any field-guided method, solution quality depends on the global field topology. Here, for example, our default field generation strategy struggles on meshes with challenging limit cycles (left). Switching to automatic cone placement yields a flip-free parameterization (center) that is easily quantized and contoured to generate a quad mesh (right).</figcaption></figure>

shown in the inset), and flatten it using the procedure from Sections 5.1 and 5.2. We then glue back together the edge ij that incurs the smallest residual $F_{ij}$ (Equation 17). If the angle $\beta_{ij}$ between frames across the glued edge is closest to a nonzero multiple of $\pi/2$, we add this edge to the curve $\Gamma$. This process is repeated until all edges are glued back together, at which point $\Gamma$ determines the configuration of singularities (à la Section 5.3.1). To make this process more efficient, we glue batches of edges together after each solve, rather than just a single edge. In particular, we glue the largest batch of edges whose absolute error $|F_{ij}|$ sums to less than $\pi$.

<figure><img src="image_28.png" /></figure>

Like Myles et al., we do not make any claim of global optimality— but do obtain good results in practice, even for challenging field topology (Figure 20). Conversely, optimal cones for conformal parameterization can perform poorly for rectangular parameterization (Figure 19). As always, there is a tradeoff between number of singularities (which reduce field smoothness) and overall distortion. We control this tradeoff with a regularization term

$$ \Phi_{\text{smooth}} := \epsilon_{\text{smooth}} \sum_{ij \in E} w_{ij}((\omega_{12}^0)_{ij} + \theta_{ijk} - \theta_{jil})^2. $$

The effect of $\epsilon_{\text{smooth}}$ is illustrated in Figure 15.

# 6 Evaluation and Examples

## 6.1 Implementation

We implemented our algorithm in MATLAB, in double precision. Timings were measured a 2019 iMac with a 3.6GHz Intel Core i9 processor and 32GB of RAM. On our largest examples (about 50k input triangles), our method took no more than about a minute for the parameterization step (excluding the cost of reference field generation and mesh extraction).

### 6.1.1 Reference Field Generation.
To generate a reference field $X^0$ with prescribed singularities, we use the *trivial connections* method of Crane et al. [2010]. Alternatively, for organic surfaces, we start with a *globally optimal direction field à la* Knöppel et al. [2013], then improve this field via iterative *Ginzburg-Landau (GL)* smoothing [Viertel and Osting 2019]. Finally, for hard surface (e.g., CAD) models, we modify the definition of parallel transport near acute corners *à la* Desobry et al. [2022], identify sharp edges where the angle between face normals is greater than 40°, and again apply GL smoothing. For curvature-aligned examples, we use the method of Cohen-Steiner and Morvan [2003] (using a stencil of nine edges from four adjacent triangles), plus a small amount of GL smoothing.

### 6.1.2 Mesh Extraction.
For meshing examples, we need to (i) quantize the parameterization $f$ computed by our algorithm onto a regular grid, and (ii) contour the quantized parameterization to obtain mesh elements. For the first step we use the method of Coudert-Osmont et al. [2024]. For the second step we use the QEx method [Ebke et al. 2013], though since the method of Coudert-Osmont et al. [2024] already requires the parameterization to be locally injective, we do not rely on the robustness properties of this method (and simply use the QEx code for convenience).
<figure><img src="image_29.png" /><figcaption>Fig. 21. Empirically, we tend to get the same solution independent of how we initialize the solver—suggesting that our results may be globally optimal (even though the formulation is nonconvex). Here we plot the initial scale factors *u*, *v* and frame orientation *θ*, obtaining the same solution across constant, random, and structured initializations. We consistently observe this behavior for examples throughout the paper.</figcaption></figure>

## 6.2 Initialization, Optimality, and Robustness

Given the nonconvexity of our optimization problem, it is natural to wonder if we always get a high-quality solution. Importantly, to get a valid, integrable field, it is sufficient to find any feasible point satisfying Equation 17—we *do not need a globally optimal solution*. In practice, however, we do appear to obtain optimal solutions, since very different initialization strategies lead to identical results. E.g., in Figure 21 we initialize with (1) *u*, *v* = 0 (no scaling) and *θ* = 0 (no rotation of the reference frame), (2) *θ* = 0 and *u*, *v* sampled uniformly at random from [-1, 1], (3) *u*, *v* = 0 and uniformly random *θ*, (4) uniformly random *u*, *v*, and *θ*, (5) an intrinsic conformal parameterization where *v* = 0, *u* solves *Δu* = *dω*⁰, and *θ* solves *Δθ* = *d* *ω*⁰, (6) *u*, *v* and *θ* obtained by breadth-first integration of a nonintegrable reference field, (7) *u*, *v*, *θ* extracted from LSCM [Lévy et al. 2002], and (8) *u*, *v*, *θ* extracted from LSCM followed by ARAP [Liu et al. 2008]. If we intentionally try to break the solver (e.g., using large random initial values) it can get stuck near a saddle point, but we did not encounter this problem naturally on any of our models. In practice, we therefore default to just *u* = *v* = *θ* = 0.

As seen in Figure 22, RSP is also robust to difficult singularity configurations. In all cases we obtain a valid rectangular parameterization, and highly rectangular quads (away from singular points).

<figure><img src="image_30.png" /><figcaption>Fig. 22. Starting with an integrable rectangular field helps with robustness throughout the rest of the pipeline. Here we obtain valid, highly-rectangular meshes even for difficult cone configurations that are structured (top left), random (top right), or have increasingly large index (bottom row).</figcaption></figure>

## 6.3 Convergence

Under refinement of the input mesh *M*, we observe linear convergence to a smooth rectangular parameterization. Consider for instance Figure 24. Here, letting (*Y*₁, *Y*₂) be the Jacobian of the map *f*, *non-orthogonality* is the deviation of the angle between *Y*₁ and *Y*₂ from π/2; *integration error* is the Frobenius norm of the difference between (*X*₁⁰, *X*₂⁰) and (*Y*₁, *Y*₂); *angle deviation* is the angle between (*X*₁⁰, *X*₂⁰) and the rotation matrix *R* closest to (*Y*₁, *Y*₂); and *scale/stretch deviation* is the difference between *u*, *v* and the log magnitudes of (*Y*₁, *Y*₂), *resp.* In each case we plot the area-weighted *ℓ*₂ norm. Figure 23 again indicates that this error is quite small in absolute terms even on coarse meshes, and becomes smaller for fine meshes. Overall, Equation 17 appears to do an excellent job of characterizing integrability even at the discrete level.

<figure><img src="image_31.png" /><figcaption>Fig. 23. As predicted by the smooth theory, our discrete frame fields are almost perfectly integrable. Here we plot numerical integration error as a function of mesh size, for the 88 MAMBO models in our dataset; notice that error tends to get smaller as the mesh gets finer. As visualized on the left, integration error is tiny even for very coarse meshes.</figcaption></figure>
<figure><img src="image_32.png" /><figcaption>Fig. 24. Our smooth theory guarantees that fields are exactly orthogonal, and exactly integrable. In the discrete case, any failure of integrability is hence due only to discretization error—which vanishes under mesh refinement. Here we compare the gradient of the final piecewise linear parameterization (shown in the background) to the input field, achieving a linear reduction in error with respect to mean edge length of the input mesh.</figcaption></figure>

## 6.4 Quadrilateral Mesh Generation

A central motivation for orthogonal frame field generation is field-guided quadrilateral remeshing [Bommes et al. 2013b]. As noted in Section 1, rectangular elements are especially desirable in this setting due to, e.g., better asymptotic convergence rates for FEM [Arnold et al. 2002], superior mechanical properties [Pottmann et al. 2007], and optimal geometric approximation (Section 6.4.2). However, an orthogonal frame alone is not sufficient to obtain near-rectangular mesh elements: in general, the initial frame field will be distorted during the parameterization step (Figure 3), quantization of the parameterization (to “snap” it to a regular grid [Bommes et al. 2013a]), and extraction of quadrilateral elements [Ebke et al. 2013], where angles around irregular-degree vertices must deviate from 90°.

Though some distortion is inevitable, starting with a near-integrable rectangular field significantly improves mesh quality relative to state of the art methods. Here we compare against open source implementations of *Instant Meshes* [Jakob et al. 2015], *QuadriFlow* [Huang et al. 2018], and *QuadWild* [Pietroni et al. 2021], as well as *Exoside* [Rouca 2019], commercial software invoked by thousands of users via integration with *ZBrush*. For all methods we use the authors’ reference implementation with default parameters, targeting approximately the same element count.

### 6.4.1 Hard Surface Models.
We evaluated performance on hard surface meshing using 88 CAD models from the MAMBO dataset [Ledoux 2022] (Figure 8). Overall, RSP yields more rectangular elements than past methods, on all but a few examples where the mean difference in worst angle is around just 1° (Figure 25). (See supplemental material for mesh data and more detailed statistics.) We omit models where quantization via [Coudert-Osmont et al. 2024] failed—though our field optimization and parameterization

<figure><img src="image_33.png" /><figcaption>Fig. 25. Increase in error relative to RSP (on the MAMBO dataset), quantified as the mean deviation of corner angles from 90°. Overall, our method produces more rectangular elements than established alternatives, in turn providing better geometry (e.g., planarity and discrete conjugacy) and faster solver convergence (e.g., for quadrilateral FEM).</figcaption></figure>

steps succeed for all models. In a real implementation, one could achieve greater robustness with a less simplistic integration step and/or more robust quantization (see Section 7).

### 6.4.2 Organic Models.
Asymptotically, the best quadrilateral approximation of a positively-curved smooth surface aligns mesh edges with principal curvature directions [D’Azevedo 1999]. To get such a mesh, we first compute principal directions $E_1, E_2$ and curvatures $\kappa_1, \kappa_2$ in each triangle, via the method of Cohen-Steiner and Morvan [2003]. These directions are smoothed via a vector diffusion equation $\frac{\partial}{\partial t} E_i = a \Delta^\nabla E_i$ to obtain the reference field $X^0$, where $\Delta^\nabla$ is the connection Laplacian of Knöppel et al. [2013], and the variable diffusion rate $a := 1/(|K| + \varepsilon)$ helps adapt smoothing to local feature size [Sharp and Crane 2018, Section 5.6] (where K is Gaussian curvature and $\varepsilon = 10^{-3}$). We then minimize an objective

$$ \Phi_{\text{curv}}(u, v, \theta) := \Phi_{\text{align}}(u, v, \theta) + \|v - v^0\|_2^2, $$

where $\| \cdot \|_2$ is the area-weighted $\ell_2$ norm, and $v^0 := \frac{1}{2} \log(r)$ encourages rectangles to have a target aspect ratio of $r := |\kappa_1/\kappa_2|$, which is optimal for minimizing $L_\infty$ gradient error [Shewchuck 2002; Nebel and Chern 2023]. Examples are shown in Figure 27, and a comparison with other methods is provided in Figure 26.

## 6.5 Finite Element Simulation

Beyond standard synthetic metrics of element quality, we evaluated how overall mesh quality affects performance in end-to-end simulation of *microfluidic circuits*. These systems drive fluid through narrow channels to induce interactions between chemical compounds or biological species, enabling “lab on a chip” devices [Oh et al. 2012]. Since channels are shallow relative to overall circuit dimensions, they are well-modeled by two-dimensional low-Reynolds-number flows. To track mixing dynamics we use quadrilateral FEM
<figure><img src="image_34.png" /><figcaption>Fig. 26. For remeshing smooth, organic shapes our integrable rectangular approach leads to good normal approximation, and a natural flow of mesh topology that improves even on retopology tools in the commercial ZBrush modeling package (Exoside). Here we show both flat shaded (gray) and corresponding Catmull-Clark subdivision (yellow) surfaces.</figcaption></figure>

<figure><img src="image_35.png" /><figcaption>Fig. 27. Curvature-guided parameterizations (checkerboard pattern), showing cuts (yellow) and positive/negative cone points (red/blue), as well as the resulting quad mesh (visualized as Catmull-Clark surfaces). Histograms show the distribution of corner angles, which are tightly centered around 90°.</figcaption></figure>

for Stokes flow, and a low-order upwind finite volume scheme for advection-diffusion. In such problems, meshes cannot be tailored to a single solver, but must work well across different discretizations.

Figure 29 demonstrates a gradient generator, where compounds progressively flow through several serpentine mixers to achieve a uniform concentration gradient. Such mixers compensate for the lack of turbulent mixing in laminar flow. RSP provides an anisotropic quad mesh which, due to boundary alignment, naturally follows flow characteristics. We first solve a time-independent Stokes' equation

$$ \Delta \mathbf{u} + \nabla p = 0, \quad \nabla \cdot \mathbf{u} = 0, $$

for velocity **u** and pressure **p**. Positive Dirichlet conditions on **u** at inlets drives two species through the circuit. We use zero-Dirichlet conditions on **p** at the outlet, and zero-Neumann conditions everywhere else. We then solve a time-dependent advection equation

$$ \frac{d}{dt} c = -\nabla \cdot (c\mathbf{u}) $$

for the concentration $\phi$, until reaching equilibrium. At each step we update the concentration $c$ in each quad via

$$ c = c + \tau \frac{1}{A} \sum_{i=1}^{4} \ell_i \mathbf{n}_i \cdot \mathbf{u}_i c_i^+ $$

where $A$ is the quad area, $\ell_i$ is the length of the ith edge, $\mathbf{n}_i$ is its inward unit normal, $\mathbf{u}_i$ is the fluid velocity, and $c_i^+$ is the upwind
<figure><img src="image_36.png" /><figcaption>Fig. 28. Though 2D quad meshing is often viewed as a “solved problem,” it is deceptively hard to obtain correct results for end-to-end simulation. Here we test a variety of mesh generators on our microfluidic simulation. Most meshes led to INF/NaN values; the few that gave meaningful values still failed to reproduce the piecewise concentration observed in the real physics, due to numerical diffusion (even at higher element counts than RSP).</figcaption></figure>

concentration at edge i, equal to the value of c in either the quad or its neighbor depending on the sign of $n_i \cdot u_i$. The time step $\tau > 0$ is chosen according to a CFL condition, namely, $\tau$ must be smaller than the minimum over all quads of $A/\sum_{i=1}^4 |u_i \cdot l_i \cdot n_i|$.

As shown in Figure 29, simulation on the mesh generated by RSP yields not only the expected velocity and linear pressure profile, but also a concentration gradient stratified over the number of mixers in the final stage—just as in physical experiments. Figure 28 compares with other mesh generators, including those listed in Section 6.4, plus tools from MeshLab [Tarini et al. 2010] and Gmsh [Remacle et al. 2013; Baudouin et al. 2014]. Though some of these methods are “robust” in the sense that they reliably generate a topologically valid mesh, element quality can still be quite poor. As a result, most meshes lead to invalid solutions with INF/NaN values. Moreover,

even the successful examples (Exoside, Instant Meshes) exhibit non-physical diffusion—even with significantly more elements than the RSP mesh. Post-processing to improve element orthogonality will not solve this problem, since not all mesh connectivity admits a rectangular mesh aligned to flow characteristics. The fact that RSP provides explicit control over singularities and boundary behavior is hence a major advantage for real-world anisotropic mesh generation.

## 6.6 Chebyshev Nets

A Chebyshev net is a mapping $\phi(x, y) : U \to \mathbb{R}^n$ over a region $\Omega \subset \mathbb{R}^2$ such that

$$ \left| \frac{\partial \phi}{\partial x} \right| = \left| \frac{\partial \phi}{\partial y} \right| = d $$

for some constant $d > 0$, i.e., that preserve lengths in the two coordinate directions. Such nets provide a model of loosely woven textiles made from inextensible threads, which easily shear but do not scale or stretch (cheesecloth, gauze, fishnet, etc.). Finding a Chebyshev parameterization of a surface hence provides a cutting pattern for which the surface can be made out of a piece of cloth. For instance, a hemisphere can be perfectly covered by a single diamond-shaped piece of inextensible cloth, à la Figure 30, top center.

A slight tweak of RSP enables us to compute high-quality Chebyshev parameterizations. The basic observation is that infinitesimal squares in the plane are mapped to infinitesimal rectangles of constant diagonal, so long as the sum of squares of scale factors is constant (see inset). Since a, b are scale factors going from the surface to the plane, we can express this condition as

<figure><img src="image_37.png" /></figure>

$$ 1/a^2 + 1/b^2 = d^2. \qquad (21) $$

Contouring diagonals rather than isolines (equivalently: rotating $f$ by 45°) then yields a Chebyshev net (Figure 30).

For simplicity we let $d = 1$ (since any other $d$ is obtained via global rescaling), and encourage Equation 21 via a corresponding penalty, which is convex in log scale factors:

$$ \Phi_{\text{Chebyshev}}(u, v, \theta) := \sum_{ijk \in F} A_{ijk} \left( -2u_{ijk} + \log \left( e^{2v_{ijk}} + e^{-2v_{ijk}} \right) \right)^2, $$

where $u_{ijk} := \frac{1}{3}(u_i + u_j + u_k)$, and likewise for $v_{ijk}$.

Figure 31 shows that our maps are both smoother and closer to Chebyshev nets than those from the special-purpose algorithm of Oehri et al. [2024]. Like Oehri et al. [2024], and unlike Sageman-Furnas et al. [2019] or Liu et al. [2020], we can generate Chebyshev nets without any singularities—needed for fabricating garments from standard woven fabric (which has no irregular points).

## 7 Limitations and Future Work

Since Equation 19 is nonconvex, the solver could get stuck at a local minimum or infeasible point. In practice, however, we find our optimization strategy not only succeeds reliably, but produces
<figure><img src="image_38.png" /><figcaption>Fig. 29. Our method is suitable for generating high-quality 2D rectangular meshes of domains with complex geometry, as needed for numerical simulation. Here we simulate a microfluidic *gradient generator*, which uses a series of *serpentine mixers* to produce a more uniform concentration gradient of two input species. *Left:* schematic and photograph of a real gradient generator. *Center:* a rectangular mesh automatically generated via RSP captures flow characteristics, without any unwanted singularities. *Right:* simulation on this mesh accurately models the concentration gradient (qualitatively matching the photograph), even with a low-order numerical method.</figcaption></figure>

<figure><img src="image_39.png" /><figcaption>Fig. 30. Chebyshev nets are broadly used in architecture, textile design or, as shown here, to design literal *sheet bend* nets that wrap around curved surfaces. Our method can be used to generate high-quality Chebyshev nets by contouring diagonals of a rectangular parameterization rather than its isolines.</figcaption></figure>

<caption>Fig. 31. By simply reinterpreting our rectangular parameterizations, we reliably obtain lower-distortion Chebyshev nets than methods explicitly designed for this purpose. Here we measure error using the metric of Oehri et al. [2024, Equation 12].</caption>
the same high-quality minimizer independent of initialization (Sec-
tion 6.2). This situation is no different from existing field genera-
tion techniques, which are largely nonconvex [Viertel and Osting
2019; Sageman-Furnas et al. 2019; Jakob et al. 2015; Diamanti et al.
2015], apart from a few exceptions [Crane et al. 2010; Knöppel et al.
2013]. RSP is also of course more expensive than linear integrability
schemes—but well-worth it, given that the real-world engineering
bottleneck is typically human time spent on tweaking mesh genera-
tion and repairing defective meshes.

As discussed in Section 2.1, local injectivity can in principle fail
due to discretization error, though rarely happens in practice. If local
injectivity is crucial, one might replace our simple Poisson integra-
tion with a more robust approach [Myles et al. 2014]. An interesting
question is whether there is an exact, mimetic discretization of rect-
angular maps in the spirit of discrete conformal equivalence [Gillespie
et al. 2021b], which guarantees local injectivity.

Another question is whether existence of a rectangular parameter-
ization implies existence of a *seamless* rectangular parameterization
with the same cones/constraints. We have never found a counterex-
ample, and suspect the answer is “yes”, but lack a formal proof.

Our automatic cone placement strategy already works quite well (Section 5.3.2), but since singularity placement for *rectangular* parameterization is a new topic, there is of course room for improvement. Here one might build on recent insights from the conformal case [Soliman et al. 2018; Fang et al. 2021b; Li et al. 2022, 2023].

Since our method (like many in this space) discretizes a partial differential equation, its behavior will of course be influenced by the input mesh quality. However, since RSP is purely intrinsic—and we have been careful to express its discretization in intrinsic terms—it can be implemented without modification on *intrinsic triangulations* [Sharp et al. 2021], where *intrinsic Delaunay refinement* can be used to guarantee robustness on essentially arbitrarily bad input meshes [Sharp et al. 2019; Gillespie et al. 2021a].

Finally, our original motivation for this approach was establishing better tools for volumetric *hexahedral* meshing [Pietroni et al. 2022; Fang et al. 2021a]. Since our basic strategy (writing out derivatives of a moving rectangular frame) does not depend inherently on a 2D domain (or holomorphic functions), we are optimistic about its extension to 3D. However, significant open questions remain, *e.g.*, how to move from point-like singularities to 1D curve networks.

Acknowledgments

This work was supported by the French Agence Nationale de la
Recherche (ANR) under grant ANR-24-CE48-0895-01 (ORTHOMAP),
the National Science Foundation under award 1943123, a Packard
Foundation Fellowship, as well as gifts from Adobe and nTopology.

References

Pierre Alliez, David Cohen-Steiner, Olivier Devillers, Bruno Lévy, and Mathieu Desbrun. 2003. Anisotropic polygonal remeshing. In *ACM SIGGRAPH 2003*.

Douglas Arnold, Daniele Boffi, and Richard Falk. 2002. Approximation by quadrilateral finite elements. *Mathematics of computation* 71, 239 (2002).

Thierry Aubin. 2013. *Some nonlinear problems in Riemannian geometry*. Springer Science & Business Media.

Tristan Carrier Baudouin, Jean-François Remacle, Emilie Marchandise, François Hen-
rotte, and Christophe Geuzaine. 2014. A frontal approach to hex-dominant mesh
generation. *Adv. Mod. Sim. Eng. Sci.* 1 (2014), 1–30.

Mirela Ben-Chen, Craig Gotsman, and Guy Bunin. 2008. Conformal flattening by curvature prescription and metric scaling. In *Computer Graphics Forum*, Vol. 27.

David Bommes, Marcel Campen, Hans-Christian Ebke, Pierre Alliez, and Leif Kobbelt. 2013a. Integer-grid maps for reliable quad meshing. *ACM Trans. Graph*. 32, 4 (2013).

David Bommes, Bruno Lévy, Nico Pietroni, Enrico Puppo, Claudio Silva, Marco Tarini, and Denis Zorin. 2013b. Quad-mesh generation and processing: A survey. In *Computer Graphics Forum*, Vol. 32.

David Bommes, Henrik Zimmer, and Leif Kobbelt. 2009. Mixed-integer quadrangulation. *ACM Trans. Graph*. 28, 3 (2009).

Guy Bunin. 2008. A continuum theory for unstructured mesh generation in two dimensions. *Computer Aided Geometric Design* 25, 1 (2008).

Marcel Campen, Hanxiao Shen, Jiaran Zhou, and Denis Zorin. 2019. Seamless parametrization with arbitrary cones for arbitrary genus. *ACM Trans. Graph*. 39, 1 (2019).

Ryan Capouellez and Denis Zorin. 2024. Seamless Parametrization in Penner Coordinates. *ACM Trans. Graph*. 43, 4 (2024).

Élie Cartan et al. 2001. *Riemannian Geometry in an Orthogonal Frame: From Lectures Delivered by Lie Cartan at the Sorbonne in 1926-1927*. World Scientific.

Edward Chien, Zohar Levi, and Ofir Weber. 2016. Bounded distortion parametrization in the space of metrics. *ACM Trans. Graph*. 35, 6 (2016).

David Cohen-Steiner and Jean-Marie Morvan. 2003. Restricted delaunay triangulations and normal cycle. In *Symp. Comp. Geom*.

Guillaume Coiffier and Etienne Corman. 2023. The Method of Moving Frames for Surface Global Parametrization. *ACM Trans. Graph*. 42, 5 (2023).

Yoann Couderd-Osmont, David Desobry, Martin Heistermann, David Bommes, Nicolas Ray, and Dmitry Sokolov. 2024. Quad Mesh Quantization Without a T-Mesh. In *Computer Graphics Forum*, Vol. 43.

Matteo Couplet, Alexandre Chemin, and Jean-François Remacle. 2024. Integrable Frame Fields using Odeco Tensors. In *Proc. Int. Mesh. Roundtable*. SIAM.

Keenan Crane, Mathieu Desbrun, and Peter Schröder. 2010. Trivial Connections on Discrete Surfaces. *Computer Graphics Forum* (SGP) 29, 5 (2010).

Mathieu Desbrun, Eva Kanso, and Yiying Tong. 2006. Discrete differential forms for computational modeling. In *ACM SIGGRAPH 2006 Courses*.

David Desobry, François Protais, Nicolas Ray, Etienne Corman, and Dmitry Sokolov. 2022. Frame Fields for CAD models. *Lecture Notes in Computer Science* 13018 (Jan. 2022), 421–434. https://inria.hal.science/hal-03537852

Olga Diamanti, Amir Vaxman, Daniele Panozzo, and Olga Sorkine-Hornung. 2015. Integrable polyvector fields. *ACM Trans. Graph*. 34, 4 (2015).

Alexander Dielen, Isaak Lim, Max Lyon, and Leif Kobbelt. 2021. Learning direction fields for quad mesh generation. In *Computer Graphics Forum*, Vol. 40.

Manfredo P Do Carmo. 2016. *Differential geometry of curves and surfaces: revised and updated second edition*. Courier Dover Publications.

Eduardo D'Azevedo. 1999. On optimal bilinear quadrilateral meshes. *Engineering with Computers* 15 (1999), 219–227.

Hans-Christian Ebke, David Bommes, Marcel Campen, and Leif Kobbelt. 2013. QEx: Robust quad mesh extraction. *ACM Trans. Graph*. 32, 6 (2013).

Qing Fang, Wenqing Ouyang, Mo Li, Ligang Liu, and Xiao-Ming Fu. 2021b. Computing sparse cones with bounded distortion for conformal parameterizations. *ACM Trans. Graph*. 40, 6 (2021).

Xianzhong Fang, Jin Huang, Yiying Tong, and Hujun Bao. 2021a. Metric-driven 3D frame field generation. *IEEE Trans. Vis. and Comp. Graph*. 29, 4 (2021).

Michael S Floater and Kai Hormann. 2005. Surface parameterization: a tutorial and survey. *Advances in multiresolution for geometric modelling* (2005).

Blender Foundation. 2025. Blender: a 3D modelling and rendering package. https://www.blender.org. Version 3.6, accessed January 22, 2025.

Xiao-Ming Fu, Jian-Ping Su, Zheng-Yu Zhao, Qing Fang, Chunyang Ye, and Ligang Liu. 2021. Inversion-free geometric mapping construction: A survey. *Computational Visual Media* 7 (2021), 289–318.

Mark Gillespie, Nicholas Sharp, and Keenan Crane. 2021a. Integer Coordinates for Intrinsic Geometry Processing. *ACM Trans. Graph*. 40, 6 (2021).

Mark Gillespie, Boris Springborn, and Keenan Crane. 2021b. Discrete conformal equivalence of polyhedral surfaces. *ACM Trans. Graph*. 40, 4 (2021).

Zekun Hao, David W Romero, Tsung-Yi Lin, and Ming-Yu Liu. 2024. Meshtron: High-Fidelity, Artist-Like 3D Mesh Generation at Scale. *arXiv:2412.09548* (2024).

Aaron Hertzmann and Denis Zorin. 2000. Illustrating smooth surfaces. In *SIGGRAPH*.

Jingwei Huang, Yichao Zhou, Matthias Niessner, Jonathan Richard Shewchuk, and Leonidas J Guibas. 2018. Quadriflow: A scalable and robust method for quadrangulation. In *Computer Graphics Forum*, Vol. 37.

Alec Jacobson, Daniele Panozzo, C Schüller, Olga Diamanti, Qingnan Zhou, N Pietroni, et al. 2013. libigl: A simple C++ geometry processing library. (2013).

Wenzel Jakob, Marco Tarini, Daniele Panozzo, Olga Sorkine-Hornung, et al. 2015. Instant field-aligned meshes. *ACM Trans. Graph*. 34, 6 (2015).

Jovana Jezdimirović, Alexandre Chemin, and Jean-François Remacle. 2023. Integrable Cross-Field Generation Based on Imposed Singularity Configuration—The 2D Manifold Case. In *International Meshing Roundtable*. Springer.
Tengfei Jiang, Xianzhong Fang, Jin Huang, Hujun Bao, Yiying Tong, and Mathieu Desbrun. 2015. Frame field generation through metric customization. *ACM Trans. Graph.* 34, 4 (2015).

Felix Kälberer, Matthias Nieser, and Konrad Polthier. 2007. Quadcover-surface parametrization using branched coverings. In *Comp. Graph. For.*, Vol. 26.

Liliya Kharevych, Boris Springborn, and Peter Schröder. 2006. Discrete Conformal Mappings via Circle Patterns. *ACM Trans. Graph.* 25, 2 (apr 2006), 412–438.

Felix Knöppel, Keenan Crane, Ulrich Pinkall, and Peter Schröder. 2013. Globally optimal direction fields. *ACM Trans. Graph.* 32, 4 (2013).

Felix Knöppel, Keenan Crane, Ulrich Pinkall, and Peter Schröder. 2015. Stripe Patterns on Surfaces. *ACM Transactions on Graphics (TOG)* 34, 4 (2015), 1–11.

Denis Kovacs, Ashish Myles, and Denis Zorin. 2010. Anisotropic quadrangulation. In *Proc. Solid and Phys. Mod.*

Frank Ledoux. 2022. MAMBO. URL: https://gitlab.com/franck.ledoux/mambo (2022).

Zohar Levi. 2023a. Seamless Parametrization with Cone and Partial Loop Control. *ACM Trans. Graph.* 42, 5 (2023).

Zohar Levi. 2023b. Shear-reduced seamless parametrization. *Computer Aided Geometric Design* 101 (2023).

Bruno Lévy, Sylvain Petitjean, Nicolas Ray, and Jérôme Maillot. 2002. Least squares conformal maps for automatic texture atlas generation. *ACM Trans. Graph.* 21, 3 (2002).

Mo Li, Qing Fang, Wenqing Ouyang, Ligang Liu, and Xiao-Ming Fu. 2022. Computing sparse integer-constrained cones for conformal parameterizations. *ACM Trans. Graph.* 41, 4 (2022).

Mo Li, Qing Fang, Zheng Zhang, Ligang Liu, and Xiao-Ming Fu. 2023. Efficient Cone Singularity Construction for Conformal Parameterizations. *ACM Trans. Graph.* 42, 6 (2023).

Hao-Yu Liu, Zhong-Yuan Liu, Zheng-Yu Zhao, Ligang Liu, and Xiao-Ming Fu. 2020. Practical fabrication of discrete chebyshev nets. In *Computer Graphics Forum*, Vol. 39.

Ligang Liu, Lei Zhang, Yin Xu, Craig Gotsman, and Steven J Gortler. 2008. A local/global approach to mesh parameterization. In *Comp. Graph. For.*, Vol. 27.

Feng Luo. 2004. Combinatorial Yamabe flow on surfaces. *Communications in Contemporary Mathematics* 6, 05 (2004), 765–780.

Richard Henri MacNeal. 1949. The solution of partial differential equations by means of electrical networks. Ph.D. Dissertation. California Institute of Technology.

Martin Marinov and Leif Kobbelt. 2004. Direct anisotropic quad-dominant remeshing. In *Pac. Conf. Comp. Graph. App. IEEE*.

Dorit Merhof, Roberto Grosso, Udo Tremel, and Günther Greiner. 2007. Anisotropic quadrilateral mesh generation: an indirect approach. *Adv. Eng. Soft.* 38, 11-12 (2007).

Ashish Myles, Nico Pietroni, and Denis Zorin. 2014. Robust field-aligned global parametrization. *ACM Trans. Graph.* 33, 4 (2014).

Ashish Myles and Denis Zorin. 2013. Controlled-distortion constrained global parametrization. *ACM Trans. Graph.* 32, 4 (2013).

Nicolas Nebel and Albert Chern. 2023. Adaptive Surface Meshes from Harmonic Maps. *arXiv e-prints* (2023).

Annika Oehri, Aviv Segall, Jing Ren, and Olga Sorkine-Hornung. 2024. Chebyshev Parameterization for Woven Fabric Modeling. *ACM Trans. Graph.* 43, 6 (2024).

Kwang W Oh, Kangsun Lee, Byungwook Ahn, and Edward P Furlani. 2012. Design of pressure-driven microfluidic networks using electric circuit analogy. *Lab on a Chip* 12, 3 (2012).

Barrett O'Neill. 2006. *Elementary differential geometry*. Elsevier.

Barrett O'Neill. 2014. *The geometry of Kerr black holes*. Courier Corporation.

David Palmer, David Bommes, and Justin Solomon. 2020. Algebraic representations for volumetric frame fields. *ACM Trans. Graph.* 39, 2 (2020).

Davide Pellis and Helmut Pottmann. 2018. Aligning principal stress and curvature directions. In *Advances in Architectural Geometry*.

Nico Pietroni, Marcel Campen, Alla Sheffer, Gianmarco Cherchi, David Bommes, Xifeng Gao, Riccardo Scateni, Franck Ledoux, Jean Remacle, and Marco Livesu. 2022. Hex-mesh generation and processing: a survey. *ACM Trans. Graph.* 42, 2 (2022).

Nico Pietroni, Stefano Nuvoli, Thomas Alderighi, Paolo Cignoni, Marco Tarini, et al. 2021. Reliable feature-line driven quad-remeshing. *ACM Trans. Graph.* 40, 4 (2021).

Kacper Pluta, Michal Edelstein, Amir Vaxman, and Mirela Ben-Chen. 2021. PH-CPF: planar hexagonal meshing using coordinate power fields. *ACM Trans. Graph.* 40, 4 (2021).

Helmut Pottmann, Andreas Asperl, and Axel Kililan. 2007. *Architectural geometry*. SIAM.

Michael Rabinovich, Tim Hoffmann, and Olga Sorkine-Hornung. 2018. Discrete geodesic nets for modeling developable surfaces. *ACM Trans. Graph.* 37, 2 (2018).

Michael Rabinovich, Roi Poranne, Daniele Panozzo, and Olga Sorkine-Hornung. 2017. Scalable locally injective mappings. *ACM Trans. Graph.* 36, 4 (2017).

Nicolas Ray, Wan Chiu Li, Bruno Lévy, Alla Sheffer, and Pierre Alliez. 2006. Periodic global parameterization. *ACM Trans. Graph.* 25, 4 (2006).

J-F Remacle, François Henrotte, T Carrier-Baudouin, Eric Béchet, E Marchandise, Christophe Geuzaine, and Thibaud Mouton. 2013. A frontal Delaunay quad mesh generator using the L-∞ norm. *Int. J. Num. Meth. Eng.* 94, 5 (2013).

Maxime Rouca. 2019. Exoside Quad Remesher. https://exoside.com/ Automatic quad remeshing software, available as plugins for various 3D software, including Autodesk 3ds Max, Autodesk Maya, Autodesk Fusion 360, Modo, Blender, and Houdini.

Andrew O. Sageman-Furnas, Albert Chern, Mirela Ben-Chen, and Amir Vaxman. 2019. Chebyshev nets from commuting PolyVector fields. *ACM Trans. Graph.* 38, 6 (2019).

Pedro V Sander, John Snyder, Steven J Gortler, and Hugues Hoppe. 2001. Texture mapping progressive meshes. In *SIGGRAPH*.

Rohan Sawhney and Keenan Crane. 2017. Boundary First Flattening. *ACM Trans. Graph.* 37, 1, Article 5 (Dec. 2017), 14 pages.

Rohan Sawhney and Keenan Crane. 2019. boundary-first-flattening. https://github.com/GeometryCollective/boundary-first-flattening. Accessed January 22, 2025.

Nick Sharp and Keenan Crane. 2018. Variational Surface Cutting. *ACM Trans. Graph.* 37, 4 (2018).

Nicholas Sharp, Mark Gillespie, and Keenan Crane. 2021. Geometry Processing with Intrinsic Triangulations. (2021).

Nicholas Sharp, Yousuf Soliman, and Keenan Crane. 2019. Navigating Intrinsic Triangulations. *ACM Trans. Graph.* 38, 4 (2019).

Alla Sheffer, Emil Praun, Kenneth Rose, et al. 2007. Mesh parameterization methods and their applications. *Found. Trends Comp. Graph. Vis.* 2, 2 (2007).

Hanziao Shen, Leyi Zhu, Ryan Capouellez, Daniele Panozzo, Marcel Campen, and Denis Zorin. 2022. Which cross fields can be quadrangulated? Global parameterization from prescribed holonomy signatures. *ACM Trans. Graph.* 41, 4 (2022).

Tianchang Shen, Zhaoshuo Li, Marc Law, Matan Atzmon, Sanja Fidler, James Lucas, Jun Gao, and Nicholas Sharp. 2024. SpaceMesh: A Continuous Representation for Learning Manifold Surface Meshes. In *SIGGRAPH Asia 2024*.

J Shewchuck. 2002. What is a good linear finite element. *Interpolation, Conditioning, Anisotropy, and Quality Measures* (2002).

Yawar Siddiqui, Antonio Alliegro, Alexey Artemov, Tatiana Tommasi, Daniele Sirigatti, Vladislav Rosov, Angela Dai, and Matthias Nießner. 2024. MeshGPT: Generating triangle meshes with decoder-only transformers. In *Proc. Comp. Vis. Pat. Rec.*

Lance Simons and Nina Amenta. 2024. Anisotropy and Cross Fields. In *Computer Graphics Forum*, Vol. 43.

Yousuf Soliman, Dejan Slepčev, and Keenan Crane. 2018. Optimal Cone Singularities for Conformal Flattening. *ACM Trans. Graph.* 37, 4 (2018).

Sanghyun Son, Matheus Gadelha, Yang Zhou, Zexiang Xu, Ming C Lin, and Yi Zhou. 2024. DMesh: A Differentiable Representation for General Meshes. *arXiv e-prints* (2024).

Boris Springborn, Peter Schröder, and Ulrich Pinkall. 2008. Conformal equivalence of triangle meshes. In *ACM SIGGRAPH 2008 papers*.

Marco Tarini, Nico Pietroni, Paolo Cignoni, Daniele Panozzo, and Enrico Puppo. 2010. Practical quad mesh simplification. In *Comp. Graph. For.*, Vol. 29.

Amir Vaxman, Marcel Campen, Olga Diamanti, Daniele Panozzo, David Bommes, Klaus Hildebrandt, and Mirela Ben-Chen. 2016. Directional field synthesis, design, and processing. In *Comp. Graph. For.*, Vol. 35.

Ryan Viertel and Braxton Osting. 2019. An approach to quad meshing based on harmonic cross-valued maps and the Ginzburg-Landau theory. *SIAM Journal on Scientific Computing* 41, 1 (2019).

Muyang Zhang, Jin Huang, Xinguo Liu, and Hujun Bao. 2010. A wave-based anisotropic quadrangulation method. In *ACM SIGGRAPH 2010*.

## A Transformation of Connection Forms under Rotation

We here give a proof of Lemma 3.1, showing that the connection form for a frame rotated by an angle $\theta$ can be expressed as

$$ \omega_{12}^{\theta} = \omega_{12}^{0} - d\theta, $$

where $\omega_{12}^0$ is the connection form for $X^0$. First note that

$$ d\sigma_1^\theta = -\sin(\theta)d\theta \wedge \sigma_1^0 + \cos(\theta)d\sigma_1^0 + \cos(\theta)d\theta \wedge \sigma_2^0 + \sin(\theta)d\sigma_2^0. $$

Making substitutions $d\sigma_1^0 = -\omega_{12}^0 \wedge \sigma_1^0$ and $d\sigma_2^0 = -\omega_{21}^0 \wedge \sigma_2^0$ from the first structure equation (Equation 1) and collecting terms gives

$$ d\sigma_1^\theta = -(\omega_{12}^0 - d\theta) \wedge (-\sin(\theta)\sigma_1^0 + \cos(\theta)\sigma_2^0). $$

From Equation 9, we know that the second factor is equal to $\sigma_2^0$. Hence, we get $d\sigma_1^\theta = -(\omega_{12}^0 - d\theta) \wedge \sigma_2^0$. Again invoking the first structure equation $d\sigma_1^\theta = -\omega_{12}^0 \wedge \sigma_2^0$, we must therefore have $\omega_{12}^0 = \omega_{12}^0 - d\theta$, as desired.
# Rectangular Surface Parameterization—Supplemental Information

ETIENNE CORMAN, Université de Lorraine, CNRS, LORIA, France

KEENAN CRANE, Carnegie Mellon University, USA

**ACM Reference Format:**

Etienne Corman and Keenan Crane. 2025. Rectangular Surface Parameterization—Supplemental Information. *ACM Trans. Graph.* 44, 4 (August 2025), 6 pages. https://doi.org/10.1145/3731176

## B Discrete Poisson Equation

For completeness, we describe the standard cotan discretization of the Poisson equation [MacNeal 1949, Section 3.2], used to recover our final parameterization coordinates in Section 5.2.

### B.1 Discrete Laplacian

The cotan-Laplace matrix $L$ is then a sparse $|V| \times |V|$ matrix with nonzero entries

$$L_{i,j} = \begin{cases} -w_{ij}, & ij \in E, \\ \sum_{ik} w_{ik} & i = j. \end{cases} \qquad (1)$$

### B.2 Discrete Divergence

The discrete divergence operator is a sparse $|V| \times |E|$ matrix with nonzero entries

$$\mathrm{div}_{i,ij} = s_{ij} w_{ij}, \quad \forall i \in V, \quad ij \in E, \qquad (2)$$

where $s_{ij} = +1$ when edge $ij$ is oriented from $i$ to $j$, and equals $-1$ otherwise.

### B.3 Discrete Poisson Equation

Finally, we discretize the Poisson equation in Equation 13 as a scalar Poisson equation for each coordinate $f^1, f^2$ of the final rectangular parameterization $f: V \to \mathbb{R}^2$, namely

$$L f^p = \mathrm{div} \mu^p, \quad p = 1, 2. \qquad (3)$$

Here $L$, $div$, and $\mu$ are the discrete Laplacian, divergence, and differential given in Equations 1, 2, and 18, resp.

## C Underdetermined Newton Solver

To solve our main optimization problem (Equation 19), we need a method that can deal with underdetermined nonlinear constraints. We use the method of Polyak and Tremba [2020], detailed here. Similar to Newton's method, this method finds a zero of the first order optimality conditions by repeatedly solving linear systems, using a carefully-chosen line search.

<figure><img src="image_40.png" /><figcaption>Fig. 32. Here we show boundary conditions like those from Figure 18, but without any cone singularities. This example illustrates both the robustness of our method to extreme constraints—and also emphasizes the utility of singularities for obtaining low distortion/good element quality.</figcaption></figure>

More explicitly, let $x := (u, v, \theta) \in \mathbb{R}^{2|V|+|F|}$ be the primal variables and $\lambda \in \mathbb{R}^{|E|}$ be the Lagrange multiplier, the critical points of Equation 19 are solutions of the system of equations:

$$\begin{gathered} \nabla\Phi(x) + J_F(x)^\top\lambda = 0 \\ F(x) = 0, \end{gathered} \qquad (4)$$

where $\nabla\Phi$ denote gradient of the objective chosen in Section 5.1.1, and $J_F$ is the Jacobian of the constraint equation given in Section 4.1.

The nonlinear system of equations in Equation 4 is solved by taking Newton steps, using a line search that directly encourages satisfaction of this system. More explicitly, let $x_k$ and $\lambda_k$ denote the current guess at the $k$th step of optimization. Then we solve the linear system

$$\begin{bmatrix} \nabla^2 \Phi(x_k) + \nabla(J_F(x_k)^\top \lambda_k) & J_F(x_k)^\top \\ J_F(x_k) & 0_{|E| \times |E|} \end{bmatrix} \begin{bmatrix} y_k \\ \delta_k \end{bmatrix} = \begin{bmatrix} \nabla \Phi(x_k) + J_F(x_k)^\top \lambda_k \\ F(x_k) \end{bmatrix} \qquad (5)$$

<footnote>Authors' Contact Information: Etienne Corman, Université de Lorraine, CNRS, LORIA, Vandœuvre-lès-Nancy, France; Keenan Crane, Carnegie Mellon University, Pittsburgh, PA, USA.</footnote>
for the descent directions $y_k$ and $\delta_k$, where $\nabla^2\Phi$ denote the Hessian the objective function. The updated guesses are $x_{k+1} = x_k - \tau y_k$ and $\lambda_{k+1} = \lambda_k - \tau \delta_k$, where $\tau > 0$ is a step size. The step size is computed by performing backtracking line search until

$$ \|\nabla\Phi(x_k - \tau y_k) + J_F(x_k - \tau y_k)^\top(\lambda_k - \tau \delta_k)\| + \|F(x_k - \tau y_k)\| \le (1 - \frac{\tau}{2}) (\|\nabla\Phi(x_k) + J_F(x_k)^\top\lambda_k\| + \|F(x_k)\|). $$

This strategy guarantees that we converge to a local minimizer of the original optimization problem (Equation 19), at essentially the usual rate of Newton's method—see Polyak and Tremba [2020, Theorem 4.1]. In practice, we use a backtracking ratio of 0.9 for line search. Crucially, we find that Equation 5 always admits a solution and is only overconstrained when inadmissible boundary constraints are imposed (Section 5.1.3).

## D Singularities and Cuts

In the case where either (i) the field has singularities and/or (ii) the domain $M$ is not a topological disk, we must effectively cut along the curve $\Gamma$, plus additional topological cuts, before flattening the mesh. In practice, we compute topological cuts by considering the traversal tree used to construct $\Gamma$, and taking the complementary graph of primal edges not crossed by any edge in this dual tree. We then iteratively remove any degree-1 vertex from the primal graph, except for singular vertices, until no more vertices can be removed. The remaining primal edges define the cut. (This approach is similar in spirit to the tree-cotree strategy of Eppstein [2002].)

This problem of solving for a parameterization on a cut mesh comes up frequently in geometry processing, and is a perennial nuisance. Hence, we will first introduce a perspective that simplifies the situation both conceptually and in terms of implementation, before giving the details of our particular algorithm. We make just one simple change to our basic setup: rather than store values at vertices, we will store values at all triangle corners, using linear constraints to identify equivalent values as needed. In the absence of any cuts, our formulation ultimately yields the same solution as before, but with the additional flexibility to easily incorporate cuts.

In more detail, for any function $f_i$ at vertices, we now store values $f_i^{jk}$ at all triangle corners. A sparse matrix $U \in \mathbb{R}^{m \times |C|}$ encodes equivalence of values at different corners, where $m$ is the number of identifications, and $|C|$ is the number of corners. Two corners around a common vertex are equivalent if and only if they are not separated by a cut edge. For the $n$th equivalence $f_i^{jk} = f_i^{kl}$ of consecutive corners, $U$ will have nonzero entries $U_{n,ik} = 1$ and $U_{n,il} = -1$, i.e., the equation number determines the row index, and the corners are specified via the column indices.

Next, suppose we want to solve any linear equation $Af = b$, which in general may be under- or over-determined. We can express this system in terms of values on corners, plus the identification matrix $U$, as a least-squares problem subject to linear constraints:

$$ \begin{array}{ll} \min_{f \in \mathbb{R}^{|C|}} & \|Af - b\|_M^2 \\ \text{s.t.} & Uf = 0. \end{array} \qquad (6) $$

Here $|A|_M^2 := A^\top MA$ denotes the $\ell$-2 norm with a mass matrix $M \in \mathbb{R}^{|C|^2 \times |C|}$. We use in particular a diagonal matrix with nonzero

entries

$$ M_{jk, \ i, \ i}^{jk} = \frac{1}{2} \cot \alpha_i^{jk}, $$

i.e., we use the cotan weight associated with each corner. In the case where there is no cut, our problem amounts to solving the exact same cotan-Poisson problem in Equation 3.

Equation 6 can be solved using any available method—for instance, the method of Lagrange multipliers, we find the solution to this problem is given by the block linear system

$$ \begin{bmatrix} A^\top M A & U^\top \\ U & 0 \end{bmatrix} \begin{bmatrix} f \\ \lambda \end{bmatrix} = \begin{bmatrix} A^\top M b \\ 0 \end{bmatrix}, $$

where $\lambda \in \mathbb{R}^m$ are the Lagrange multipliers. By formulating the problem this way, rather than directly applying least-squares to the combined system $[AU]^\top f = [b0]$, we ensure that corner values are identified exactly, and only the original equations are approximated. Note that in general there may be more efficient ways to solve this kind of problem—we find that this approach is simple, and helps avoid implementation errors. For instance, one could directly eliminate variables by hand (though modern direct solvers are often quite good at effectively performing this same elimination); a general-purpose QP solver might also work fine depending on the use case.

In the case of our particular problem, we also have to modify Equation 3, splitting it into two different equations on the two sides of every edge. In particular, for each edge $ij$ in $\Gamma$ we modify the right-hand side so that both vectors $\mu_{ij}^k, \mu_{ji}^l$ are in the same coordinate system. In particular, consider the angles $\beta_{ij}$ from Section 5.3.1, which give the rotation between frames (after parallel transport) across any cut edge $ij$. Hence, to evaluate Equation 3 for an edge from corner $f_i^{jk}$ to corner $f_j^{ki}$, we use

$$ f_j^{ki} - f_i^{jk} = \frac{1}{2}(\mu_{ij}^k + R_{-\zeta_{ij}} \mu_{ji}^l), \qquad (7) $$

where $\zeta_{ij} := \frac{\pi}{2}[(\beta - \pi/4)/(\pi/2)]$ is the closest quarter-rotation taking the frame in triangle $ijk$ to the frame in triangle $jil$. In other words, we account for the jump in the frame across the cut. Likewise, on the other side of the edge, with endpoints at corners $i^{lj}$ and $j^{il}$, we have

$$ f_j^{il} - f_i^{lj} = \frac{1}{2}(R_{\zeta_{ij}} \mu_{ij}^k + \mu_{ji}^l). \qquad (8) $$

## References

* David Eppstein. 2002. Dynamic generators of topologically embedded graphs. *arXiv preprint cs/0207082* (2002).

* Richard Henri MacNeal. 1949. *The solution of partial differential equations by means of electrical networks*. Ph.D. Dissertation. California Institute of Technology.

* Boris Polyak and Andrey Tremba. 2020. New versions of Newton method: step-size choice, convergence domain and under-determined equations. *Optimization Methods and Software* 35, 6 (2020).

* B. Stellato, G. Banjac, P. Goulart, A. Bemporad, and S. Boyd. 2020. OSQP: an operator splitting solver for quadratic programs. *Mathematical Programming Computation* 12, 4 (2020), 637–672.
## E Pseudocode

Here we give detailed pseudocode for our rectangular parameterization algorithm. We omit pre-processing (namely, generation of the reference field) and post-processing (namely, quantization and contouring of the parameterization in the case of meshing), since these steps are handled by existing algorithms. All other methods are given explicitly in terms of basic operations on vectors and a standard halfedge mesh data structure, and should be easily translatable to real code. The only subroutines not spelled out here are:

* `UNIT(v)` — given a nonzero vector $v \in \mathbb{R}^3$, returns the unit vector in the same direction.

* `OBJECTIVE(M, A, u, v, θ)`, `GRADOBJECTIVE(M, A, u, v, θ)`, `HESSOBJECTIVE(M, A, u, v, θ)` — implements one of the objective functions from Section 5.1.1, and its gradient/Hessian, given the mesh connectivity $M$, the triangle areas $A$, log scale factors $u, v$ per vertex, frame rotation angles $θ$ per face.

* `SOLVE(A, b)` — returns $x$ solution of $Ax = b$

* `SOLVEQP(H, g, J, F)` — solves a quadratic program of the form $\min_y \frac{1}{2} y^T H y + y^T g$ subject to $Jy + F = 0$. (We use OSQP [Stelato et al. 2020].)

* `SPARSEFROMTRIPLETS((I, J, S), n, m)` — returns a sparse matrix $A$ of size $n \times m$ whose coefficient $A_{I_j J_i}$ is equal to $S_i$ for all $i$. If a coefficient appears several times in the lists $(I, J)$, their corresponding values in $S$ are summed.

* `APPEND(L, a)` — returns a new list appending the element $a$ to the list $L$.

* `ANGLE(x, y, z)` — returns the angle between the vector $y - x \in \mathbb{R}^3$ and $z - x \in \mathbb{R}^3$.

For simplicity we assume $M$ is without boundary.

### Algorithm 1 FLATTENMESH(M, x, W)

**Input:** A manifold, orientable, simply-connected triangle mesh $M = (V, E, F)$ with corners $C$, coordinates $x_i \in \mathbb{R}^3$ at each vertex $i \in V$, and a cross field given by any one of four representative unit vectors $W_{ijk} \in \mathbb{R}^3$ in each triangle $ijk \in F$

**Output:** Coordinates $f: C \to \mathbb{R}^2$ describing the flattening of each triangle into the 2D plane. Note that these values are stored at triangle corners rather than vertices to accommodate possible cuts.

1. **for each** $ij \in E$ **do** $\ell_{ij} \leftarrow |x_i - x_j|$ ▶ edge lengths

2. **for each** $j^k_i \in C$ **do** $\alpha_i^{j^k} \leftarrow \mathrm{ANGLE}(x_i, x_j, x_k)$ ▶ corner angles

3. **for each** $ijk \in F$ **do** ▶cross field angles

4. $N_1 \leftarrow (x_j - x_i) \times (x_k - x_i)$ ▶unnormalized normal

5. $A_{ijk} \leftarrow |N_1|/2$ ▶triangle area

6. $N_1 \leftarrow N_1/(2A_{ijk})$ ▶normalize

7. $T_1 \leftarrow \mathrm{UNIT}(x_j - x_i)$

8. $T_2 \leftarrow N_1 \times T_1$

9. $\xi_{ijk} \leftarrow \mathrm{atan2}(\langle T_2, W_{ijk} \rangle, \langle T_1, W_{ijk} \rangle)$ ▶angle relative to $ij$

10. $\Gamma, \zeta, s, \varphi, \omega^0 \leftarrow \mathrm{COMPUTECUTJUMPDATA}(M, \alpha, \xi)$

11. $u, v, \theta \leftarrow \mathrm{SOLVEOPTIMIZATIONPROBLEM}(M, A, \alpha, \varphi, \omega^0, s)$

12. $f \leftarrow \mathrm{RECOVERPARAMETERIZATION}(M, \Gamma, \zeta, \ell, \alpha, \varphi, \theta, u, v)$

13. return $f$

### Algorithm 2 COMPUTECUTJUMPDATA(M, α, ξ)

**Input:** A manifold, orientable triangle mesh $M = (V, E, F)$ with angles $\alpha_i^{jk}$ at each corner $i$ of each triangle $ijk \in F$, angles $\xi: F \to \mathbb{R}$ describing one of the four directions of the cross field in each triangle $ijk \in F$ as a rotation relative to the direction of the first edge $ij$.

**Output:** The curve $\Gamma: E \to \{0, 1\}$ cutting $M$ to a topological disk, the closest quarter-rotation $\zeta: E \to \mathbb{R}$ between two frame in adjacent triangle, sign bits $s: C \to \pm 1$ giving the sign of $v$ at each corner relative to the values stored at vertices, angles $\varphi: H \to \mathbb{R}$ which give the reference frame relative to each halfedge, and angles $\omega^0: E \to \mathbb{R}$ describing the rotation of the reference frame $X^0$ along each oriented dual edge.

1. **for each** $ijk \in F$ **do** ▶for convenience, express cross directions

2. $\xi \bar{i} \bar{j} \leftarrow \xi_{ijk}$

3. $\xi_j \bar{j} \leftarrow \xi \bar{i} - (\pi - \alpha_j^{ki})$

4. $\xi \bar{k} \bar{i} \leftarrow \xi_j \bar{j} - (\pi - \alpha_k^{ij})$

▶Perform breadth-first traversal (Section 5.3.1)

5. $ijk \leftarrow \mathrm{FIRST}(F)$ ▶get the first triangle in the mesh

6. **for each** $ijk \in F$ **do** $\varphi_{ijk} \leftarrow \infty$ ▶frame not yet set

7. $\varphi_{ijk} \leftarrow \xi_{ijk}$ ▶reference frame gets one of the four cross directions

8. $\mathrm{PUSH}(Q, \bar{i} \bar{j})$ ▶push any halfedge of $ijk$ onto a queue

9. **while** `NOTEMPTY`(Q) **do**

10. $\bar{i} \bar{j} \leftarrow \mathrm{POP}(Q)$ ▶go from $ijk$ to $jil$

11. $\varphi_{ij} \leftarrow \varphi_{ij} + \pi$ ▶parallel transport

12. $n^* \leftarrow \mathrm{MOD}(\mathrm{ROUND}(2\frac{\varphi_{ij} \to j - \xi_{ji}}{\pi}), 4)$ ▶closest cross index

13. $\xi_{ij} \leftarrow \frac{\pi}{2} n^*$ ▶jump angle across $ij$

14. $\xi^* \leftarrow \xi_{ij} + \xi_{ji}$ ▶closest cross

15. $\omega_{ij}^0 \leftarrow \varphi_{ij} - \xi^* + \pi$ ▶smallest rotation to neighboring cross

16. **if** $\varphi_{ji} == \infty$ **then** ▶neighboring frame not yet set

17. $\Gamma_{ij} \leftarrow 0$ ▶primal edge not in cut

18. $s_{ij} \leftarrow +1$ ▶no sign change

19. $\varphi_{ji} \leftarrow \xi^*$ ▶set reference frame in triangle

20. $\varphi \bar{j} \bar{k} \leftarrow \varphi \bar{i} \bar{j} - (\pi - \alpha_j^{ki})$

21. $\varphi \bar{k} \bar{i} \leftarrow \varphi \bar{j} \bar{k} - (\pi - \alpha_i^{ij})$

22. $\mathrm{PUSH}(Q, \bar{i} \bar{j})$ ▶visit other two neighbors

23. $\mathrm{PUSH}(Q, \bar{j} \bar{i})$ ▶visit other two neighbors

24. **else** ▶neighboring frame already set

25. $\Gamma_{ij} = 0$ ▶primal edge could be in cut

26. **if** `IsOdd`(n*) **then** $s_{ij} = -1$

27. **else** $s_{ij} = +1$

28. **end while** ▶Compute relative signs at corners

29. **for each** $i \in V$ **do** ▶iterate over vertices

30. $S = 1$ ▶cumulative product (Equation 20)

31. **for** $p \leftarrow 0, \dots, \mathrm{DEGREE}(M, i) - 1$ **do**

32. $s_i^{J_p J_{p+1}} \leftarrow S$

33. $S \leftarrow s_{ij}^{p+1} S$ ▶compute angle defect
34: $K_i \leftarrow 0$

35: **for each** $ijk \in F$ **do**

36: $K_i \leftarrow K_i + \alpha_i^{jk}$

37: $K_j \leftarrow K_j + \alpha_j^{ki}$

38: $K_k \leftarrow K_k + \alpha_k^{ij}$

39: $K \leftarrow 2\pi - K$

►compute cones indices at vertices

40: $c_i \leftarrow K$

41: **for each** $i \in V$ **do**

42: **for each** $j \in V$ such that $ij \in E$ **do**

43: $w \leftarrow \omega_{ij}^0$

44: **if** $i > j$ **then** $w \leftarrow -w$

45: $c_i \leftarrow c_i + w$

►edge incident to i

►accumulate all rotations

►compute the cut

46: $r \leftarrow TRUE$

47: **while** $r$ **do**

48: $d \leftarrow 0$

►initialize degree of unvisited edges

49: **for each** $\Gamma_{ij} == 1$ **do**

►iterate over unvisited edges

50: $d_i \leftarrow d_i + 1$

51: $d_j \leftarrow d_j + 1$

52: $r \leftarrow FALSE$

53: **for each** $i \in V$ **do**

►iterate over vertices

54: **if** $d_i == 1$ and $c_i == 0$ **then** ►if degree-1 and not a cone

55: **for each** $ij \in E$ **do** $\Gamma_{ij} = 0$

►remove edge

56: $r \leftarrow TRUE$

57: **end while**

58: **return** $(\Gamma, \zeta, s, \varphi, \omega^0)$

**Algorithm 3** SOLVEOPTIMIZATIONPROBLEM(M, A, α, φ, ω⁰, s)

**Input:** A manifold, orientable, simply-connected triangle mesh $M = (V, E, F)$, areas $A_{ijk}$ for each triangle $ijk \in F$, angles $\alpha_i^{jk}$ at each corner $j^k \in C$, angles $\varphi_{ij}$ describing the angle of the reference frame relative to each halfedge $\overline{ij} \in H$, angles $\omega_{ij}^0$ describing the rotation of the reference frame across each edge $ij \in E$, and values $s_i^{jk} \in \{\pm 1\}$ giving the sign of the value at each corner $j^k$ relative to a value stored at vertex $i$.

**Output:** Scale factors $u_i, v_i \in \mathbb{R}$ at each vertex $i \in V$, and angles $\theta_{ijk} \in \mathbb{R}$ for each triangle $ijk \in F$.

1: $u, v \leftarrow 0 \in \mathbb{R}^{|V|}$

►initialize scale factors

2: $\theta \leftarrow 0 \in \mathbb{R}^{|F|}$

►initialize frame angles

3: $\lambda \leftarrow 0 \in \mathbb{R}^{|E|}$

►initialize Lagrange multiplier

4: **while** TRUE **do**

5: $H \leftarrow HESSOBJECTIVE(M, A, u, v, \theta)$

6: $g \leftarrow GRADOBJECTIVE(M, A, u, v, \theta)$

7: $D \leftarrow BUILDHESSIAN(M, \alpha, u, v, \lambda, s, \varphi, \theta)$

►constraint

Hessian

8: $J \leftarrow BUILDJACOBIAN(M, \alpha, u, v, s, \varphi, \theta)$ ►constraint Jacobian

9: $F \leftarrow BUILDSYSTEM(M, \alpha, u, v, s, \varphi, \theta, \omega^0)$ ►constraint value

10: $x \leftarrow (u, v, \theta, \lambda)$

►current guess

11: $A \leftarrow \begin{bmatrix} H+D & J^T \\ J & 0_{|E| \times |E|} \end{bmatrix}$

►linear system

12: $b \leftarrow \begin{bmatrix} g + J^T \lambda \\ F \end{bmatrix}$

►right hand side

13: $y \leftarrow SOLVE(A, b)$

14: $\tau \leftarrow LINESEARCH(M, \alpha, s, \varphi, \omega^0, J, F, g, x, y)$

15: $(u, v, \theta, \lambda) \leftarrow x + \tau y$

►update guess

16: **if** $\|g + J^T \lambda\| + \|F\| < \varepsilon$ **then**

17: **break**

18: **end while**

19: **return** $(u, v, \theta)$

**Algorithm 4** LINESEARCH(M, α, s, φ, ω⁰, J, F, g, x, y)

**Input:** A manifold, orientable triangle mesh $M = (V, E, F)$ with angles $\alpha_i^{jk}$ at each corner $i$ of each triangle $ijk \in F$, sign bits $s: C \to \pm 1$ giving the sign of $v$ at each corner relative to the values stored at vertices, angles $\varphi: H \to \mathbb{R}$ giving the angle of the reference frame relative to each halfedge, angles $\omega^0: E \to \mathbb{R}$ describing the rotation of the reference frame $X^0$ across each edge, a sparse matrix $J \in \mathbb{R}^{|E| \times (2|V|+|F|)}$ giving the Jacobian of the constraint, values $F: E \to \mathbb{R}$ giving the discrete residual of Equation 17, a vector $g \in \mathbb{R}^{2|V|+|F|}$ giving the gradient of the objective function, two vectors $x, y \in \mathbb{R}^{2|V|+|F|}$ describing the current estimate and the descent direction

**Output:** Step size $\tau > 0$ decreasing the Lagrangian gradient norm

1: $(u, v, \theta, \lambda) \leftarrow x$

2: $\mathcal{E} \leftarrow \|g + J^T \lambda\| + \|F\|$

►norm of gradient Lagrangian

3: $\tau \leftarrow 1$

►step size

4: **while** TRUE **do**

5: $(u, v, \theta, \lambda) \leftarrow x + \tau y$

►update estimate

6: $g \leftarrow GRADOBJECTIVE(M, A, u, v, \theta)$

►gradient of objective

7: $J \leftarrow BUILDJACOBIAN(M, \alpha, u, v, s, \varphi, \theta)$

8: $F \leftarrow BUILDSYSTEM(M, \alpha, u, v, s, \varphi, \theta, \omega^0)$

9: $\mathcal{E}_{\text{new}} \leftarrow \|g + J^T \lambda\| + \|F\|$

►update gradient norm

10: **if** $\mathcal{E}_{\text{new}} \le (1 - \frac{\tau}{2}) \mathcal{E}$ **then**

►break

►stop when sufficient decrease

11: **else**

12: $\tau \leftarrow 0.9\tau$

►decrease step size

13: **end while**

14: **end while**

15: **return** $\tau$

**Algorithm 5** BUILDSYSTEM(M, α, u, v, s, φ, θ, ω⁰)

**Input:** A manifold, orientable triangle mesh $M = (V, E, F)$ with angles $\alpha_i^{jk}$ at each corner $i$ of each triangle $ijk \in F$, values $u, v: V \to \mathbb{R}$ at each vertex giving the log scale factors, sign bits $s: C \to \pm 1$ giving the sign of $v$ at each corner relative to the values stored at vertices, angles $\varphi: H \to \mathbb{R}$ giving the angle of the reference frame relative to each halfedge, angles $\theta: F \to \mathbb{R}$ describing the current frame $X_1^\theta$ as a rotation of the reference frame, angles $\omega^0: E \to \mathbb{R}$ describing the rotation of the reference frame $X^0$ across each edge.

**Output:** Values $\rho: E \to \mathbb{R}$ giving the residual of our discrete equation (Equation 17) on each edge.

1: **for each** $ijk \in E$ **do**

►left-hand side of Equation 10

2: $\rho_{ij} \leftarrow (\theta_{ijk} - \theta_{jil}) - \omega_{ij}^0$

►k and l are left/right of ij
3: **for each ijk $\in F$ do** ▶*right-hand side of Equation 10* ▶*add contribution to residual from dual edges inside this triangle*

4: $\rho_{ij} \leftarrow \rho_{ij} + \text{RESIDUAL}(i, j, k, \alpha, u, v, s, \varphi_{\vec{i}j} + \theta_{ijk})$

5: $\rho_{jk} \leftarrow \rho_{jk} + \text{RESIDUAL}(j, k, i, \alpha, u, v, s, \varphi_{\vec{j}k} + \theta_{ijk})$

6: $\rho_{ki} \leftarrow \rho_{ki} + \text{RESIDUAL}(k, i, j, \alpha, u, v, s, \varphi_{\vec{k}i} + \theta_{ijk})$

7: return $\rho$

**Algorithm 6 RESIDUAL(i, j, k, α, u, v, s, η)**

**Input:** Three vertices $i, j, k$ of a triangle, the angles $\alpha$ at triangle corners, the log scale factors $u, v$ at vertices, the signs $s$ at triangle corners, and the angle $\eta$ of the current frame $X_1^\theta$ relative to a coordinate system aligned with edge $ij$.

**Output:** The contribution to the residual of equation Equation 17 due to the piece of the edge dual to $ij$ within triangle $ijk$.

1: $(v_i^{jk}, v_j^{ki}, v_k^{ij}) \leftarrow (s_i^{jk} v_i, s_j^{ki} v_j, s_k^{ij} v_k)$ ▶*values at corners*

2: $\rho \leftarrow u_j - u_i$ ▶*conformal part*

3: $\rho \leftarrow \rho + \cos(2\eta)(v_j^{ki} - v_i^{jk})$ ▶*non-conformal part*

4: $\rho \leftarrow \rho + \sin(2\eta) \cot \alpha_j^{ki} (v_k^{ij} - v_i^{jk})$

5: $\rho \leftarrow \rho + \sin(2\eta) \cot \alpha_i^{jk} (v_k^{ij} - v_j^{ki})$

6: $\rho \leftarrow \frac{1}{2} \cot \alpha_k^{ij} \rho$

7: **if i > j then**

8: $\rho \leftarrow -\rho$

▸account for edge orientation

9: return $\rho$

**Algorithm 7 BUILDJACOBIAN(M, α, u, v, s, φ, θ)**

**Input:** A manifold, orientable triangle mesh $M = (V, E, F)$ with angles $\alpha_i^{jk}$ at each corner $i$ of each triangle $ijk \in F$, values $u, v: V \to \mathbb{R}$ at each vertex giving the log scale factors, signs $s: C \to \pm 1$ describing the values of $v$ at corners relative to the values stored at vertices, angles $\varphi: F \to \mathbb{R}$ describing the reference frame $X_1^0$ in each triangle $ijk \in F$ as a rotation relative to the direction of the first edge $ij$, and angles $\theta: F \to \mathbb{R}$ describing the current frame $X_1^\theta$ as a rotation of the reference frame.

**Output:** A sparse matrix $J \in \mathbb{R}^{|E|(2|V|+|F|)}$ giving the Jacobian of the constraint equations computed by *BUILDSYSTEM*.

1: $L_v \leftarrow ()$ ▶*list of nonzeros of Jacobian with respect to v*

2: $L_\theta \leftarrow ()$ ▶*list of nonzeros of Jacobian with respect to $\theta$*

3: **for each $ij \in D$ do** ▶*Jacobian of left-hand side of Equation 10*

4: $L_\theta \leftarrow \text{APPEND}(L_\theta, (ij, ijk, +1))$

5: $L_\theta \leftarrow \text{APPEND}(L_\theta, (ij, jil, -1))$

6: **for each $ijk \in F$ do** ▶*Jacobian of right-hand side of Equation 10*

7: $\eta_{ij} \leftarrow \varphi_{ijk} + \theta_{ijk}$ ▶*angle of $X_1^\theta$ relative to each edge*

8: $\eta_{jk} \leftarrow \eta_{ij} - (\pi - \alpha_i^{ki})$

9: $\eta_{ki} \leftarrow \eta_{jk} - (\pi - \alpha_j^{ki})$

10: ▶*append contribution to Jacobian from dual edges in this triangle*

11: $(L_u, L_v, L_\theta) \leftarrow \text{JACOBIAN}(i, j, k, \alpha, u, v, s, \eta_{ij}, L_u, L_v, L_\theta)$

12: $(L_u, L_v, L_\theta) \leftarrow \text{JACOBIAN}(j, k, i, \alpha, u, v, s, \eta_{jk}, L_u, L_v, L_\theta)$

13: $(L_u, L_v, L_\theta) \leftarrow \text{JACOBIAN}(k, i, j, \alpha, u, v, s, \eta_{ki}, L_u, L_v, L_\theta)$

14: ▶*Build Jacobian blocks from lists of nonzeros*

15: $J_u \leftarrow \text{SPARSEFROMTRIPLES}(L_u, |E|, |V|)$

16: $J_v \leftarrow \text{SPARSEFROMTRIPLES}(L_v, |E|, |V|)$

17: $J_\theta \leftarrow \text{SPARSEFROMTRIPLES}(L_\theta, |E|, |F|)$

18: return [$J_u$, $J_v$, $J_\theta$] ▶*assemble overall $|E| \times (2|V| + |F|)$ Jacobian*

**Algorithm 8 JACOBIAN(i, j, k, α, u, v, s, η, $L_u$, $L_v$, $L_\theta$)**

**Input:** The same inputs as RESIDUAL, plus the cumulative lists of entries $L_u, L_v, L_\theta$ of the three Jacobian blocks built so far.

**Output:** Modified lists $L_u, L_v, L_\theta$ that account for the derivatives (with respect to $u, v$ and $\theta$) of the residual due to the piece of the edge dual to $ij$ within triangle $ijk$.

1: $w \leftarrow \frac{1}{2} \cot \alpha_k^{ij}$ ▶*cotan weight*

2: **if i > j then**

3: $w \leftarrow -w$ ▶*account for edge orientation*

4: ▶*conformal part w.r.t. u*

5: $L_u \leftarrow \text{APPEND}(L_u, (ij, i, -w))$

6: $L_u \leftarrow \text{APPEND}(L_u, (ij, j, +w))$

7: ▶*non-conformal part w.r.t. v*

8: $L_v \leftarrow \text{APPEND}(L_v, (ij, i, -w s_i^{jk} (\cos(2\eta) + \sin(2\eta) \cot \alpha_j^{ki})))$

9: $L_v \leftarrow \text{APPEND}(L_v, (ij, j, +w s_j^{ki} (\cos(2\eta) - \sin(2\eta) \cot \alpha_i^{jk})))$

10: $L_v \leftarrow \text{APPEND}(L_v, (ij, k, +w s_k^{ij} \sin(2\eta)(\cot \alpha_i^{ki} + \cot \alpha_j^{jk})))$

11: ▶*non-conformal part w.r.t. $\theta$*

12: $(v_i^{jk}, v_j^{ki}, v_k^{ij}) \leftarrow (s_i^{jk} v_i, s_j^{ki} v_j, s_k^{ij} v_k)$ ▶*values at corners*

13: $L_\theta \leftarrow \text{APPEND}(L_\theta, (ij, ijk, -2w \sin(2\eta)(v_k^{ki} - v_i^{jk})))$

14: $L_\theta \leftarrow \text{APPEND}(L_\theta, (ij, ijk, +2w \cos(2\eta) \cot \alpha_j^{ki} (v_k^{ij} - v_i^{jk})))$

15: $L_\theta \leftarrow \text{APPEND}(L_\theta, (ij, ijk, +2w \cos(2\eta) \cot \alpha_i^{jk} (v_k^{ij} - v_j^{ki})))$

16: return $(L_u, L_v, L_\theta)$

**Algorithm 9 BUILDHESSIAN(M, α, u, v, λ, s, φ, θ)**

**Input:** A manifold, orientable triangle mesh $M = (V, E, F)$ with angles $\alpha_i^{jk}$ at each corner $i$ of each triangle $ijk \in F$, values $u, v: V \to \mathbb{R}$ at each vertex giving the log scale factors, values $\lambda: E \to \mathbb{R}$ at each dual-edge giving the Lagrange multiplier associated to the constraint, signs $s: C \to \pm 1$ describing the values of $v$ at corners relative to the values stored at vertices, angles $\varphi: F \to \mathbb{R}$ describing the reference frame $X_1^0$ in each triangle $ijk \in F$ as a rotation relative to the direction of the first edge $ij$, and angles $\theta: F \to \mathbb{R}$ describing the current frame $X_1^\theta$ as a rotation of the reference frame.

**Output:** A sparse matrix $H \in \mathbb{R}^{(2|V|+|F|) \times (2|V|+|F|)}$ giving the Jacobian of the matrix computed by JACOBIAN multiplied by the Lagrange multiplier $\lambda$.

1: $L_u \leftarrow ()$ ▶*list of nonzeros of Jacobian with respect to u*

2: $L_v \leftarrow ()$ ▶*list of nonzeros of Jacobian with respect to v*

3: $L_\theta \leftarrow ()$ ▶*list of nonzeros of Jacobian with respect to $\theta$*

4: **for each $ijk \in F$ do** ▶*Jacobian of right-hand side of Equation 10*

5: $\eta_{ij} \leftarrow \varphi_{ijk} + \theta_{ijk}$ ▶*angle of $X_1^\theta$ relative to each edge*

6: $\eta_{jk} \leftarrow \eta_{ij} - (\pi - \alpha_j^{ki})$

7: $\eta_{ki} \leftarrow \eta_{jk} - (\pi - \alpha_i^{ki})$

8: ▶*append contribution to Jacobian from dual edges in this triangle*

9: $(L_v, L_\theta) \leftarrow \text{HESSIAN}(i, j, k, \alpha, u, v, \lambda, s, \eta_{ij}, L_v, L_\theta)$

10: $(L_v, L_\theta) \leftarrow \text{HESSIAN}(j, k, i, \alpha, u, v, \lambda, s, \eta_{jk}, L_v, L_\theta)$
11: $(L_v, L_\theta) \leftarrow \text{HESSIAN}(k, i, j, \alpha, u, v, \lambda, s, \eta_{ki}, L_v, L_\theta)$

12: ▶Build Jacobian blocks from lists of nonzeros

13: $H_v \leftarrow \text{SPARSEFROMTRIPLETS}(L_v, |V|, |F|)$

14: $H_\theta \leftarrow \text{SPARSEFROMTRIPLETS}(L_\theta, |F|, |F|)$

15: return $\begin{bmatrix} 0|_V|x|_V & 0|_V|x|_V & 0|_V|x|_F \\ 0|_V|x|_V & 0|_V|x|_V & H_v^T \\ 0|_F|x|_V & H_v & H_\theta \end{bmatrix}$ ▶assemble overall
$(2|V| + |F|) \times (2|V| + |F|)$ Jacobian

**Algorithm 10 HESSIAN(i, j, k, α, u, v, λ, s, η, L_v, L_θ)**

**Input:** The same inputs as `BUILDHESSIAN`, plus the cumulative lists of entries $L_v, L_\theta$ of the two Hessian blocks built so far.

**Output:** Modified lists $L_v, L_\theta$ that account for the derivatives (with respect to $v$ and $\theta$) of the Jacobian due to the piece of the edge dual to $ij$ within triangle $ijk$.

1: $w \leftarrow \frac{1}{2} \cot \alpha_k^{ij}$ ▶cotan weight

2: **if** $i > j$ **then**

3: $w \leftarrow -w$ ▶account for edge orientation

4: ▶non-conformal part w.r.t $v$

5: $L_v \leftarrow \text{APPEND}(L_v, (i, ijk, 2s_i^{jk} \lambda_{ij} w(\sin(2\eta) - \cos(2\eta)) \cot \alpha_j^{ki}))$

6: $L_v \leftarrow \text{APPEND}(L_v, (j, ijk, 2s_j^{ki} \lambda_{ij} w(-\sin(2\eta) - \cos(2\eta)) \cot \alpha_i^{jk}))$

7: $L_v \leftarrow \text{APPEND}(L_v, (k, ijk, 2s_k^{ij} \lambda_{ij} w(\cos(2\eta)(\cot \alpha_i^{ki} + \cot \alpha_j^{jk})))$

8: ▶non-conformal part w.r.t $\theta$

9: $(v_i^{jk}, v_j^{ki}, v_k^{ij}) \leftarrow (s_i^{jk} v_i, s_j^{ki} v_j, s_k^{ij} v_k)$ ▶values at corners

10: $L_\theta \leftarrow \text{APPEND}(L_\theta, (ijk, ijk, -4\lambda_{ij} w \cos(2\eta)(v_j^{ki} - v_i^{jk})))$

11: $L_\theta \leftarrow \text{APPEND}(L_\theta, (ijk, ijk, -4\lambda_{ij} w \sin(2\eta) \cot \alpha_j^{ki}(v_j^{ij} - v_i^{jk})))$

12: $L_\theta \leftarrow \text{APPEND}(L_\theta, (ijk, ijk, -4\lambda_{ij} w \sin(2\eta) \cot \alpha_i^{jk}(v_k^{ij} - v_j^{ki})))$

13: return $(L_v, L_\theta)$

**Algorithm 11 RECOVERPARAMETERIZATION(M, Γ, ζ, ℓ, α, φ, θ, u, v)**

**Input:** The mesh connectivity $M = (V, E, F)$, the cut curve $\Gamma: E \to \{0, 1\}$ cutting the mesh to disk topology, the closest quarter-rotation $\zeta: E \to \mathbb{R}$ between two frame in adjacent triangle, the length $\ell_{ij} \in \mathbb{R}_{\ge 0}$ of each edge $ij \in E$, angles $\alpha_i^{jk}$ at each triangle corner, the angle $\varphi_{ijk} \in \mathbb{R}$ of the reference frame relative to edge $ij$ in each triangle $ijk$, the angle $\theta_{ijk}$ of the optimized frame relative to the reference frame, and the log scale/stretch factors $u_i, v_i$ at each vertex $i \in V$.

**Output:** $f: C \to \mathbb{R}^2$ giving the uv-coordinates per triangle corner.

1: **for each** $ij \in E$ **do** ▶edge scale from log-scale

2: $a_{ij} \leftarrow e^{(u_i+u_j+v_i+v_j)/2}$

3: $b_{ij} \leftarrow e^{(u_i+u_j-v_i-v_j)/2}$

4: **for each** $ijk \in F$ **do**

5: $\eta_{ij} \leftarrow \varphi_{ijk} + \theta_{ijk}$ ▶angle of $X_1^\theta$ relative to each edge

6: $\eta_{jk} \leftarrow \eta_{ij} - (\pi - \alpha_{jij}^{ki})$

7: $\eta_{ki} \leftarrow \eta_{jk} - (\pi - \alpha_k^{ki})$

8: ▶target edge vector

9: $\mu_{ij}^k \leftarrow \ell_{ij} + a_{ij} \cos(\eta_{ijk}) + b_{ij} \sin(\eta_{ijk})$

10: $\mu_{jk}^i \leftarrow \ell_{jk} (-a_{jk} \cos(\eta_{ijk} + \alpha_j^{ki}), -b_{jk} \sin(\eta_{ijk} + \alpha_j^{ki}))$

11: $\mu_{ki}^j \leftarrow \ell_{ki} (-a_{ki} \cos(\eta_{ijk} - \alpha_i^{jk}), -b_{ki} \sin(\eta_{ijk} - \alpha_i^{jk}))$

12: ▶right hand-side in Equation ??

13: **for each** $ij \in C$ **do** $b_k^{ij} \leftarrow \frac{1}{2}(R_{\xi_{ij}} \mu_{ij}^k + \mu_{ji}^l)$

14: 15: $L_A \leftarrow ()$

16: $L_U \leftarrow ()$

17: $p \leftarrow 0$

18: **for each** $ij \in C$ **do**

19: ▶build matrix A in Equation 6

20: $L_A \leftarrow \text{APPEND}(L_A, (\binom{ij}{i}, \binom{jk}{j}, +1))$

21: $L_A \leftarrow \text{APPEND}(L_A, (\binom{ij}{j}, \binom{ki}{k}, -1))$

22: ▶build matrix U in Equation 6

23: **if** $\Gamma_{ij} == 0$ **then** ▶if $ij$ not in the cut, add corner constraints

24: $L_U \leftarrow \text{APPEND}(L_U, (p, \binom{jk}{j}, +1))$

25: $L_U \leftarrow \text{APPEND}(L_U, (p, \binom{ij}{i}, -1))$

26: $L_U \leftarrow \text{APPEND}(L_U, (p + 1, \binom{ki}{k}, +1))$

27: $L_U \leftarrow \text{APPEND}(L_U, (p + 1, \binom{li}{l}, -1))$

28: $p \leftarrow p + 2$ ▶Build Jacobian blocks from lists of nonzeros

29: $A \leftarrow \text{SPARSEFROMTRIPLETS}(L_A, |C|, |C|)$

30: $U \leftarrow \text{SPARSEFROMTRIPLETS}(L_U, p, |C|)$

31: $f \leftarrow \text{SOLVEQP}(A^\top W A, -A^\top W b, U, 0)$ ▶corner coordinates