Computing cone singularities involves determining their positions, angles, and number to reduce parameterization distortion (Fig. 1). These singularities are typically constrained to mesh vertices on the input triangular mesh, making their positions discrete variables. To construct global seamless parameterizations, which can be obtained through rotationally seamless parameterizations [7], [17]- [19], the cone angles and holonomy angles for the non-contractible homology loops must strictly satisfy the constraint of being an integer multiple of π/2. Furthermore, it is generally desirable to minimize the number of singularities. Consequently, the discrete nature of all variables introduces significant challenges in achieving a favorable tradeoff among low parameterization distortion, small number of cones, and high efficiency.

Several methods have been proposed for this challenging problem [7], [12], [13]. [7] develop a progressive curvature concentration strategy to determine the positions and number of cone singularities, while constraining their angles to integer multiples of π/2 through iterative rounding. However, the two progressive mechanisms often produce suboptimal results with excessive cone points and parameterization distortion. [12] formulate the problem as a mixed-integer optimization solved via the Douglas-Rachford splitting algorithm. Since the number of integer-constrained variables equals the input vertex number, it suffers from high computational costs and a vast search space that frequently causes the optimization to converge to local optima, ultimately failing to achieve a favorable balance between competing objectives (Fig. 2). [13] rely on the enumeration of angles and positions to co-optimize singularity count and parameterization distortion, demonstrating improved tradeoff capability. However, this exhaustive enumer-(25, 0.237, 251.0) [12] (15, 0.230, 3.7) Ours Fig. 2: Comparison with [12]. We achieve fewer cones, lower distortion, and a much shorter time (in seconds, 60× faster).

(60, 0.244, 1804.3) [13] (58, 0.237, 9.9) Ours Fig. 3: Comparison with [13]. Fewer cones and lower distortion are obtained by our method in much less time (180× faster).

ation process dramatically increases computational complexity (Fig. 3). This paper proposes an efficient method for generating sparse cone singularities with integer constraints, enabling the construction of rotationally seamless conformal parameterizations with low distortion. Our method is inspired by [20], which adopts an alternating optimization approach for position, angle, and number of cones between three steps:

• Moving cones along the distortion energy gradient.

• Optimizing cone angles to reduce distortion.

• Inserting cones using a greedy branch-driven strategy. However, [20] only works effectively when cone angles are unconstrained. Adding integer constraints into their framework causes two major challenges:

• Solving a large-scale MIQP (with many continuous and integer variables) during cone angle optimization becomes impractical.

• Optimizing cone positions with integer-constrained angles fails to decrease the distortion using their strategy. We make two key contributions to improve efficiency:

• Efficient cone angle optimization. We formulate an MIQP for cone angle optimization and exploit structural properties of the constraint matrix for variable reduction, handling both genus-zero and high-genus cases, while further limiting variables by selecting a subset of cones. • Alternating optimization with fixed angles. We alternate between moving cones with fixed angles and re-optimizing angles. Moreover, we derive a new gradient expression to effectively decrease distortion. Our algorithm efficiently balances sparse integer-constrained cone placement with low distortion. We evaluate it on a dataset comprising 3885 models to demonstrate its capability and practicability. Compared with the state-of-the-art methods [12], [13], we achieve an order of magnitude speedup (30 times faster on average) for comparable cone density and distortion.

Cone singularities [8] first introduce cone singularities for conformal parameterizations. Based on their Gaussian curvature distribution, cone singularities can be classified into two types: cones with arbitrary angles and quantized cones. Cones with arbitrary angles are primarily employed to reduce area distortion in conformal parameterizations. Typically, cones are iteratively inserted at vertices exhibiting high distortion, and the log conformal factors are computed either by solving a convex optimization problem [6] or through a diffusion process [9]. [7] propose an iterative approach that progressively flattens regions with low Gaussian curvature, ultimately identifying cone singularities at isolated vertices with high curvature. By formulating an objective that balances area distortion and cone sparsity, global optimization methods have been developed to determine the optimal placement of cones, including convex ℓ 1 optimization [11] and non-convex ℓ 0 optimization [10]. [20] enhance efficiency by iteratively relocating cones using shape derivatives and introducing new cones to reduce parameterization distortion. Quantized cones are necessary for computing globally or rotationally seamless parameterizations, where cone angles are constrained to be integer multiples of π 2 . [7] progressively round the unconstrained cone angles to the nearest multiples of π 2 , which unavoidably leads to increased area distortion. To better control distortion, the quantized cone generation is instead formulated as a mixedinteger optimization problem, which can be solved using the Douglas-Rachford scheme [12] or a series of secondorder cone programming steps [13]. However, solving mixedinteger problems at scales comparable to the input mesh vertex number is expensive. Instead, we propose several explicit constructions that dramatically reduce the problem scale for high computational efficiency.

Field singularities Recent surveys [21], [22] on vector field design review various approaches for determining the number, indices, and placement of singularities that can also be used for conformal parameterizations. On discrete triangular meshes, vector fields can be represented in several ways, including face-based [23], [24], edge-based [25], [26], and vertexbased [27]- [29] formulations. Various discretizations have been proposed, such as angle-based [30], [31], complex-number [28], complex polynomials [32], and discrete 1-forms [25], [26], [33] representations. These methods define a smoothness energy optimized to determine vector fields and singularities, grounded in discrete analogies [34], [35] to parallel transport and trivial connections on smooth manifolds. Since the parameterization is often generated by tracing vector fields, no explicit correlation exists between parameterization distortion and singularity number. [36] employ Cartan’s method of moving frames to relate parametrization to field singularities, formulating a nonlinear least-squares problem that simultaneously places singularities and minimizes distortion. Although it allows for penalizing non-conformal distortion, the large number of frame variables and multiple constraints make the algorithm inefficient.

Integer-constrained optimization in geometric processing Integer-constrained optimization is common in geometric processing applications, e.g., vector field design, cone generation, and quad meshing [22], [37]. However, the integer-constrained optimization problem is generally an NP-hard problem [38], [39]. Therefore, rounding methods are mainly adopted in these methods. The direct method performs rounding on all variables simultaneously [40]. The iterative approach incrementally rounds one continuous variable to an integer at a time, updating the remaining variables after each step [7], [30]. Two-stage methods first reduce the problem size through an approximate equivalence transformation, and then apply branch-and-bound rounding to obtain an approximate solution [41]. Unlike these rounding-based strategies, the approach in [31] begins with a feasible solution and iteratively improves it using descent steps that strictly maintain integer constraints. [12] convert the integer constraints to binary constraints and then replace them with the intersection of a box set and a sphere to solve it in a continuous domain. [13] use the second-order cone programming to approximate and solve integer-constrained optimization. We develop a reduction framework to solve the integer-constrained optimization problem with high efficiency.

Discrete conformal parameterizations For a given Riemannian manifold, the change of Gaussian curvature under a conformal mapping can be described with the Yamabe equation [42]. We discuss the case of surfaces without boundaries. For surfaces with boundaries, we defer it to supplements. For genus-zero triangular mesh M with N vertices without boundaries, given the positive-semidefinite discrete Laplacian matrix L using cotangent weights, the FEM approximation of the Yamabe equation is [6], [8], [9], [42]

where u, k tar and k ori are N -dimensional vectors and their i-th entries u i , k tar i and k ori i represent the discrete log conformal factor, the target Gaussian curvature, and the original Gaussian curvature at each vertex v i , respectively.

Area distortion Since the Hencky energy [43] remains bounded for any solution to the Yamabe equation [11], we adopt it as the metric to quantify the conformal parameterization distortion, which is defined as E = M u 2 dA 1 2 . This l 2norm form can be discretized as E(u) = √ u T Au, where A is an area weighted matrix and its i-th diagonal element A ii is computed as A ii = Area(vi) N i=1 Area(vi) , where Area(v i ) = t∈Ωi Area(t)/3 and Ω i is the set of triangles incident to v i . The problem of optimizing cone configuration As a special kind of conformal mapping, cone parameterizations require k tar to be zero except at finite cone singularities [44]. Similar to [7], [12], [13], we aim to find N c vertices C = {c i } Nc i=1 as cone singularities for conformal parameterization under the following constraints and objectives. Methodology Determining the cone configuration involves computing their discrete positions, angles, and numbers, which is a challenging discrete optimization problem. We optimize each type of variable alternately until the distortion is less than the threshold ϵ tar or the iteration number reaches the maximum (1000 in the experiments) (Alg. 1). Specifically, given a set of cones, we optimize their angles without changing positions (Sec. III-B1) and reposition cones with fixed angles (Sec. III-B2). Neither step increases distortion. When freezing the angles and positions, we adaptively add (Sec. III-C1) and remove (Sec. III-C2) cones to change the number. The resulting cones induce a global seamless conformal parameterization. Alg. 1 shows the pseudocode, and Fig. 4 shows an example.

Formulation Given a set of cones, we perform the following constrained optimization to update angles and positions without changing the number:

where z is an N -dimensional vector whose each entry is an integer. Specifically, the iteration process of relocating cones is completely shown. We also report the number of zero-angle cones (denoted as N 0 ) and the runtime t. Five progressive results are shown, where the leftmost is the initialization and the rightmost is the result.

Key idea for high efficiency If we directly solve the problem (2), the number of integer variables is the same as the number N of vertices, which is time-consuming. In practice, reducing the number of integer variables is a promising way to improve efficiency. We observe that k tar is a sparse vector with known positions of non-zero entries under fixed positions of cones. Specifically, the number of non-zero entries is the same as that of cones. Thus, we can leverage the Yamabe equation to express u in terms of those non-zeros of k tar , thereby allowing us to reduce (2) to an equivalent problem, whose number of integer variables is similar to the cone number.

Solution of the Yamabe equation Given fixed positions of cones, we sequentially collect the non-zero entries of z into a vector z int = (z c1 , …, z c Nc ) T . We can use a linear transformation to connect z int and z: z = Tz int , where T is an N × N c -dimensional matrix defined by z int . Accordingly, the Yamabe equation becomes:

The solution to the Yamabe equation in ( 3) is unique up to a constant global scaling factor. Formally, since the rank of L is (N -1) and its null space is spanned by an N -dimensional all-ones vector 1 N , the general solution of (3) is

where u par is a particular solution and a is a free continuous variable. Without loss of generality, we select a mesh vertex v p that is not a cone singularity and set its log conformal factor u p to zero to find a particular solution. This gives:

where the non-singular matrix L differs from L only in p-th row and p-th column, with the diagonal element set to one while other entries are zeros and M is an identity matrix except that its p-th diagonal entry is zero. Then, by combining ( 4) and ( 5), we have

Optimization After substituting (6) into the problem (2), optimizing cone angles under fixed positions can be reduced to

This problem has N c integer variables and one continuous variable. To perform further simplification, we can express z c Nc as a linear combination of z c1 , . . . , z c Nc -1 using the Gauss-Bonnet theorem, i.e., Nc i=1 z ci = 8. Then, this linear topology constraint can be eliminated, and we convert the problem (7) to one without linear constraints. This allows us to apply existing mixed-integer solvers to optimize the objective and determine the cone singularity angles.

Active pruning As the number of cones is too large, solving (7) using an existing mixed-integer solver becomes inefficient. Therefore, if the number of cones is greater than N g (30 in our experiment), we apply active pruning by collecting N g cones to solve angles, while the angles of all other cones are kept fixed. Since we only solve angles after each adding cone step, it is necessary to always select the newly added Initialization Using ( 9) Using (10) Fig. 5: Different methods for computing h. After placing eight cones with angles of π 2 on the belly of the Santa model (left), cones are moved along the negative gradient directions computed using h from ( 9) and (10).

cones. As newly added cones primarily influence nearby cones, we select the cones closest to them to collect a total of N g cones.

2. Optimizing positions with fixed angles Given the number and angles of cones, we adopt a similar cone relocation strategy to [20], where cones are moved across the mesh to reduce area distortion progressively. In our setting, the cone position is constrained to mesh vertices.

A revisit of [20] According to [20], [43], to reduce area distortion under the Yamabe equation, a Lagrangian function is constructed to derive the gradient. The normal directional derivative ∇ n E of the area distortion with respect to the moving direction n of a varying-angle cone c is:

where h serves as the Lagrangian multiplier and is the solution of

Here, ∂δ c is the boundary of a small neighborhood of c. h, ∇ n E can be obtained by discretizing and solving (9) [20]. Therefore, the gradient grad c E for each cone c can be derived using ∇ n E (see Equation (11) in [20]). The cone is then moved to its adjacent vertex in the direction of -grad c E with the greatest magnitude, until E no longer decreases. Our case In [20], the cone angles are continuous variables that are continuously updated during movement. However, our angles remain unchanged during the movement, leading to a different Lagrangian function and h is the solution of ∆h = -2u, on M. (10) We add the detailed discussion in supplements. As shown in Fig. 5, if we still use (9) to derive the movement direction, some cones may get stuck at the initial positions, indicating that the solution h of ( 9) is no longer practical to our problem. Instead, using (10) to compute the directional derivative leads to effective movement. Moving cones We constrain each cone to only move toward its neighboring vertices along the computed direction, similar to [20]. To sufficiently reduce the distortion while ensuring high efficiency, we perform movement using the following steps in one iteration: (1) update all cones simultaneously; (2) if Step (1) does not decrease distortion, update the cones one by one. We iteratively conduct the movement until there is no distortion reduction.

3. Algorithm analysis Non-increasing distortion In these two steps, our algorithm does not increase distortion. After solving angles with fixed positions, the optimized solution will not lead to a higher distortion than the previous. If the distortion increases within repositioning, we reject the update and retain the previous solution. Hence, the distortion does not increase (Fig. 4).

Performance When solving angles with (7), we need to evaluate π 2 L -1 MT and L -1 Mk ori . To enable pre-factorization of L, we fix the vertex v p on the mesh. As a result, L remains fixed and can be pre-factorized. Additionally, since T has a small number of columns and the system is solved infrequently, solving angles is computationally efficient. During repositioning, the main cost comes from evaluating h. Unlike the discrete formulation of (9), our formulation of (10) has a fixed coefficient matrix, which allows for one-time factorization without the need for Cholesky updates in subsequent iterations [7], [20]. Thus, our method is overall highly efficient.

It is possible to increase the number of cone singularities to lower distortion. The effectiveness of introducing additional cones via a branch-driven strategy has been confirmed by [20], [45], and we adopt a similar approach.

A revisit of the branch-driven method Let f (v) be a scalar field defined on the mesh vertices, which serves as an indicator of parameterization distortion. We set a threshold f thres > 0, and identify all vertices v such that |f (v)| > f thres . Starting from each such vertex, we perform a depth-first search to find several vertex sets {S(v)} satisfying the following conditions:

Before removing After removing Fig. 7: Removing cones. After removing a pair of cones (in the orange circle), the distortion increases slightly from 0.1176 to 0.1181, while the number of cones decreases from 26 to 24.

• For each v ∈ S(v), we have |f (v)| > f thres ;

• Each set S(v) is simply connected on the mesh;

• All values of f (v) within a set share the same sign. This process yields a collection of simply connected vertex sets, called branches. The distortion of each branch is measured by the area-weighted norm of the function f , given by:

These branches are sorted in descending order of E branch . To introduce N a new cones, the vertex with the largest |f (v)| (that is not already a cone) is added from each of the top N a branches, as shown in Fig. 6.

Adding cones for initialization To determine the initial cone positions, we use Gaussian curvature as the scalar field f to apply the branch-driven strategy. The number N a of cones is set to eight.

Adding cones during iterations To determine new cones in each iteration, we follow [20] by setting the scalar field f as u. Unlike [20], our cone angles are restricted rather than arbitrary, so adding a single cone only leads to a slight reduction in distortion. However, adding too many cones at once may lead to redundancy. Thus, we insert more cones when the current distortion is large, and fewer when it is small. Specifically, we insert min{m, 10} cones if the distortion exceeds mϵ tar and the cone number is greater than N g . Otherwise, only one cone is added to avoid introducing unnecessary cones.

The newly inserted cones are initialized with zero angles before angle optimization. After optimization, some of them may still have zero angles, but we keep them as they may become effective in subsequent iterations.

As illustrated in Fig. 7, parts of cone pairs with opposite angles tend to be ineffective in reducing distortion when they are located too close to each other. Therefore, during the iteration, we detect such pairs and attempt to remove them, similar to [13]. Specifically, we collect cone pairs with opposite angles whose distance on the mesh (measured as the number of edges between them) is less than 5×10 -4 N , as shown in Fig. 7. If removing such a pair results in a relative distortion increase smaller than a threshold η, we remove the pair; otherwise, we retain it. Since removing cones may lead to an increase in energy, to guarantee the convergence of our algorithm, we dynamically update η during the optimization. Each time a pair of cones is successfully removed, we update η to 0.9η. In the limiting case, η approaches zero, at which point we stop removing cones. Hence, the convergence of our algorithm is still guaranteed. In the experiments, the initial η is 10%.

3. Gauss-Bonnet theorem When adding cones, the newly inserted cones are initialized with zero angles. When removing cones, we always remove pairs with opposite angles. Therefore, the total cone angle remains unchanged during both insertion and removal, ensuring the Gauss-Bonnet constraint is preserved throughout.

For a closed mesh with genus g > 0, a global seamless conformal parameterization must satisfy not only the previously mentioned constraints, but also an additional set of holonomy constraints: the holonomy angles for 2g non-contractible handles and tunnels must be integer multiples of π 2 [7]. In other words, the integral of geodesic curvature along each such loop must be a multiple of π 2 :

where γ i is the i-th non-contractible homology loop and k g is the geodesic curvature of the original mesh.

In the discrete case, it suffices to require that the integral of geodesic curvature on left side of the cut is a multiple of π 2 [12]. So if considering the i-th loop on the mesh, which contains vertices G i on the left side of the cut, the discrete form of the constraint (12)

where L l j denotes the Laplacian matrix computed using only the vertex v j and its 1-ring vertices on the left side of the loop and k g i,j is the geodesic curvature of the vertex in G i . The resulting system by combining ( 13) with (2) becomes over-constrained and has no solution, indicating that a global seamless conformal parameterization does not exist in general. Thus, we instead seek a rotationally seamless parameterization, in which the log conformal factors on the two sides of each loop are allowed to differ, same goal with [12], [13]. To achieve this, we cut the mesh along the loops corresponding to the handles and tunnels, and assign independent variables to the duplicated vertices.

Our pipeline contains two steps, similar to [13]. First, we ignore the holonomy angle constraints and solve (7) to determine the positions and angles of the cone singularities. Different from Sec. III, the topology constraint becomes Nc i=1 z ci = 8(1 -g) and the number of initial cones needs to be adjusted to |8(1 -g)| accordingly. In the second step, we cut the mesh along the non-contractible homology loops and solve for a new log conformal factor while keeping the Fig. 9: New variable assignment. For the red vertices on the cut boundaries, we introduce a nonzero δ u , while for the other green vertices, δ u = 0. cone configuration fixed, such that the additional constraints are satisfied and the difference of log conformal factors on either side of the cuts remains small (Fig. 8). Next, we will focus on discussing the second step.

The enumeration-driven rounding strategy [13] for computing the integer vector r = (r 1 , …, r 2g ) T is computationally expensive. Instead, we continue to use our variable reduction approach.

Null space of constraint matrix To enable this, we assign each vertex with a continuous variable u, assemble all the constraints (3) and ( 13) into a single linear system

and compute the null space Null(L h ) of the matrix L h . This step is indispensable, as the ability to express u in terms of the integer variables r = (r 1 , . . . , r 2g ) critically depends on the structure of Null(L h ).

To analyze the null space, we first examine its dimension. After cutting the mesh, each loop leads to two duplicated vertices on either side of the cut loops, except at the intersection points between handles and tunnels, where four vertices are produced. Denote that the number of vertices on the i-th handle and tunnel loops are G h i and G t i , respectively. Then, the number of newly introduced variables after cutting is

However, only 2g additional constraints are introduced, implying that the dimension of Null(L h ) is 1+

Numerical methods Numerical solutions for computing the null space of a matrix generally rely on matrix decomposition, which requires substantial computational resources and may suffer from numerical issues [46]- [48]. Therefore, we propose a practical method to explicitly construct the null space, which is more efficient and robust than numerical methods.

Explicit construction Following the approach used for the genus-zero case, the Yamabe equation imposes N -1 constraints, and we introduce N continuous variables. As a result, the constraint matrix is rank-deficient by one. This observation inspires us to treat the higher-genus case in a similar manner: due to the presence of 2g additional constraints, we likewise aim to introduce 2g additional variables so that the constraint matrix remains rank-deficient by one, ensuring its null space is trivial. To achieve this, we carefully design the variable assignment along the cuts, while preserving the assignments of all interior vertices. We set many split vertices along the cuts to share the same variable, and additionally introduce an offset between them, similar to [49].

As illustrated in Fig. 9, we assign the same variable to two vertices on either sides, except assigning three variables to four corners at the intersection by letting one pair of adjacent corners share the same variable. In other words, the conformal factors on both sides of each cut are represented by the same variable u i , except at the intersection points, where three variables (e.g., u 3 , u 6 , u 9 in Fig. 9) are used. Compared with the original uncut mesh, only two additional variables are introduced at each intersection point. Since there are g such intersections in total, the system therefore introduces 2g additional variables. However, since the conformal factors on the two sides of each cut are inherently different, we introduce a sparse variable vector δu as a modification to u. The entries of δu are set to zero at all interior vertices, on the left-side boundary vertices, and at the three corner vertices along the boundary. Hence, the constraints become:

1. For constraints on interior vertices v i :

, where L j i is the Laplacian matrix only considering the vertex v j i and its neighbours. 3) For constraints on boundary vertices v l i , v r i on opposite side:

, where L l i , L r i is the Laplacian matrix only considering the vertex v l i , v r i and their neighbours, respectively. 4) For constraints on holonomy angels on i-th loop:

.

By designing the variables in this way and using δu to express u, we only introduce 2g additional variables, which exactly match the number of added constraints. After assembling the above four parts into a complete linear system and moving all δu terms to the right hand side, we get:

where k g = ( j∈G1 k g 1,j , . . . , j∈G2g k g 2g,j ) T and L g , K are N ′ × N ′ matrices (N ′ = N + 2g). The rank of L g is N ′ -1 (we discuss it in supplements), and L g 1 N ′ = 0 holds, where 1 N ′ is the all-ones vector. Therefore, we follow the method mentioned in Sec. III-B1 and conclude that

Therefore, in this way, instead of computing Null(L h ) numerically, we can explicitly assemble K to construct the bases of the null space.

Formulation The energy term used to penalize the jump of u across cuts can be written as

where u l j and u r j are two opposite vertices on cuts. By substituting (17), u can be eliminated, and the expression is written as

where F j is a linear vector-valued function with respect to r, δu, a. Hence, our formulation is:

where λ d is a parameter to adjust the difference of the log conformal factor on either side of cuts. As illustrated in Fig. 10, our variable reduction method outperforms [13], achieving faster computation.

We implemented our algorithm in C++ and conducted all experiments on a computer with an i7-11700 2.5GHz CPU and 16GB of RAM. The code has been released at XXX.com. The integer optimizer used was Gurobi [50], while linear matrix factorization and solving were performed using Eigen [51]. For conformal parameterization, we utilized the Boundary First Flattening method (BFF) [52], and non-contractible homology loops were computed using the method described in [53]. For surfaces with boundaries, we adopt Dirichlet boundary conditions, that is, we set b = 0.

12.2 seconds [13] 2.9 seconds Ours Fig. 10: Comparison with [13] for solving the integer variables in (13). With the same cone configuration and λ d , the same results are obtained by both methods, but ours is faster.

Solving angles after each cone movement Solving angles after each cone movement may accelerate the reduction of distortion, which is adopted by [20]. We choose to move the cones until the distortion converges without updating angles. These two strategies are compared in Fig. 11. Our scheme is much faster and achieves identical results. We observe that the cone angles remain almost unchanged during cone movement; thus, our strategy is effective. Another reason is that the cost of solving angles is much higher than that of moving cones.

Active pruning strategy The number of feasible integer combinations to (7) is |N b | Ng-1 , indicating that the number of integer variables N g (the number of active cones) has a significant impact on the efficiency of solving (7). We evaluate how different N g s affect both the optimization efficiency and the solution quality in Fig. 12. If N g is too small, the pruned solution space becomes too restricted, and the target distortion cannot be reached after the maximum number of iterations. On the other hand, if no constraint is imposed on N g , lower distortion can be achieved, but at the cost of significantly Different parameters We test our algorithm under different target distortions ϵ tar in Fig. 13. As ϵ tar decreases, both the number of cones and the computation time increase.

Our approach remains effective over a broader range of integer values, as demonstrated in Fig. 14. For instance, on the Dolphin’s fin, cone singularities of π or even 3π/2 can be accommodated. Expanding the allowable range of integers reduces the number of singularities required to achieve a comparable level of distortion.

Different meshes Next, we present the results of our algorithm on various meshes to demonstrate its adaptability.

We evaluate our algorithm on the Camel model under four different types of tessellations (regular, irregular, nonuniform, and noisy) in Fig. 15, provided by [11]. For the regular, irregular, and nonuniform tessellations, the results are similar. In the noisy tessellation, the increased number of vertices with high curvature necessitates additional cones to adequately reduce the distortion. Nonetheless, there exists a lower bound on the achievable distortion, and our algorithm is unable to reduce it to 0.200. This reveals a fundamental difference between conformal parameterization with integer constraints and that without such constraints [10], [20].

We test the Human model (Fig. 16) with three different numbers of vertices. The resulting cone number N c and distortion E are similar, and the time concerning the number of vertices is approximately linear.

Our algorithm is capable of handling nonzero genus meshes, and the rotationally seamless conformal parameterizations generated exhibit small distortion differences across the cuts. Fig. 17 shows three models with nonzero genus. Dataset statistics Similar to [12], [13], we test our algorithm on the dataset containing 3885 mechanical and organic models provided by [54] to evaluate its effectiveness, robustness, and efficiency. The average number of the mesh vertices in the dataset is 6.8k.

Three different target distortions are used, including ϵ = 0.1, 0.2 and 0.4. We collect statistics on the number of cones N c , distortion E, and computation time t (in seconds), as shown in Fig. 18. Except for the group of ϵ = 0.1, our algorithm typically reaches the target distortion bound within 10 seconds. The success rates of reaching the target distortion in prescribed number of iterations (we set it to 1000) in the three groups are 62.16%, 99.10% and 99.74%, respectively.

Local injectivity We focus on efficiently computing sparse cones, rather than computing locally injective parameterizations. Theoretically, conformal parameterizations should not exhibit flips. We adopt the method of [52] for parameterizations where injectivity is not guaranteed. In our dataset (Fig. 18), when the target distortions are set to 0.1, 0.2, and 0.4, the percentages of parameterizations in the results that contain flipped triangles are 22.1%, 19.4%, and 18.8%, respectively. Within each case, the proportion of flipped triangles is 0.2%. However, if alternative parameterization methods are used, e.g., [55], local injectivity is guaranteed.

Competitors We compare with [12] and [13] on the dataset for generating integer-constrained cones. We use the distortion outputs of these two methods as the target distortion ϵ tar s for our algorithm. Since [12] and [13] adopt different integer bounds, we ensure fairness by adopting the same bounds when comparing with each of them. The number of cones produced by each method is denoted as N cmp c for the competitors and N our c for ours. Similarly, the resulting distortion for the competitor and ours are denoted as E cmp , E our , respectively.

In Fig. 19, (X cmp -X our )/X our are illustrated, where X denotes N c , E and t ind . Here, t cmp ind and t our ind are the runtimes for each individual model processed by the competitor’s and our algorithm.

To further investigate the statistics, we categorize the results based on the relationship between cone count and distortion. We also record the total computation time for each category of results, denoted as t cmp and t our , respectively. In Tab. I, II, we report the number and proportion of results in each category, as well as for all results combined. Additionally, we present the averages and standard deviations of N cmp c -N our c and E cmp -E our . For the four categories, we present some representative models in Figs. 20, 21, 22, 23, respectively. In terms of effectiveness, compared with [12], our method produces significantly fewer cones, with 87.1% of the cases achieving both fewer cones and lower distortion. Compared with [13], our method yields similar cone counts and distortion, with slightly better averages in both metrics. Regarding computational efficiency, our method achieves a substantial speed advantage over both methods, outperforming them by more than 30 times.

Our method can produce a seamless parameterization using BFF [52], followed by quad mesh extraction through (N cmp c -N our c )/N our c (E cmp -E our )/E our (t cmp ind t our ind )/t our ind [12] [13] Fig. 19: Box plots for (X cmp -X our )/X our , where X denotes N c , E and t ind . Here, t cmp ind and t our ind are the runtimes for each individual model processed by the competitor’s and our algorithm. libQEx [56]. For comparison, we also generate quad meshes using two alternative approaches. In the first approach, we replace our cone generation method with that of [13]; the second approach is to apply the pipeline of obtaining sparse cones (or integer-constrained cross field) using IOQ [31], followed by a seamless parameterization via libigl [57] and then generating a quad mesh using libQEx, as shown in Fig. 24.

We record the statistics of the quadrilateral meshes generated by different methods, as shown in Tab. III. Since our method is based on conformal parameterization, the interior angles are more uniform, although the edge lengths exhibit relatively larger variations. In terms of computational efficiency, our method demonstrates a significant advantage.

We propose an efficient method for generating sparse cone singularities, which can be used to derive low-distortion rotationally seamless conformal parameterizations. We iteratively optimize the cone number, angles, and positions in sequence to achieve this goal. When computing cone angles and holonomy angles for the non-contractible homology loops, we avoid the complexity of the original optimization by explicitly constructing an equivalent small-scale problem, which significantly improves efficiency. The efficiency and practical robustness of our algorithm are validated on more than 3885 models.

Global seamless conformal parameterizations Even on nonzero genus meshes, a globally seamless conformal parameterization exists if the Abel-Jacobi condition is satisfied [58], [59]. However, since our cone positions are restricted to mesh vertices, this condition typically cannot be fulfilled, leading to the general nonexistence of globally seamless conformal parameterizations. One potential solution is to dynamically refine the mesh by introducing new vertices, thereby enabling cone placement at more flexible locations and facilitating the construction of a globally seamless conformal parameterization.

Limitations The main limitation of our method arises when the target distortion is set very low (e.g., 0.1), which often leads to a large number of cones and significantly slows down the optimization. A potential remedy is to design more effective pruning strategies. Additionally, the existing branchdriven and greedy cone insertion strategies are less effective for CAD models. This is because CAD surfaces typically contain numerous vertices with Gaussian curvature close to ±π/2. Directly identifying and assigning cones at these locations could lead to higher efficiency and lower distortion. Finally, due to the discrete and combinatorial characteristics, we cannot theoretically guarantee that the final distortion is less than the specified threshold. Discussion on applications Conformal parameterization has its inherent limitations for wide range of applications. On the one hand, conformal methods are theoretically restricted and cannot achieve seamless parameterization if cones are restricted to vertices; on the other hand, imposing alignment constraints also presents theoretical challenges.

However, we still consider our research promising, because we can also serve as a tool to generate cones for other types of parameterizations. In geometry processing, the cone configuration has a significant impact on the parameterization result. However, most existing methods rely on simple heuristic strategies. [7] first proposes using low-area-distortion and strict-conformal parameterization as a proxy for low-distortion isometric parameterization, thereby enabling direct optimization of the cone configuration. Subsequent studies, including ours, continue to improve the practical effectiveness of this idea. Hence, our algorithm can be regarded as a novel approach for solving the cone configuration problem, which can be directly followed by an isometric parameterization stage. Moreover, if strict conformality is not required, various existing methods [7], [18], [19] to obtain seamless parameterizations can be employed; if alignment constraints are needed, method decribed in [60] can achieve this. Thus, we believe that our algorithm remains practically useful as a foundation for further parameterization tasks.

If the Neumann boundary conditions is given, the discrete Yamabe equation becomes:

where h denotes the directional derivative of the log conformal factor along the boundary normal. The matrix on the left-hand side of ( 12) is rank-deficient by one, so the approach described in Sec. III can be adopted to handle this case. Therefore, regardless of whether Dirichlet or Neumann boundary conditions are applied, u can still be expressed in terms of z int , allowing us to derive a reduced formulation. All techniques in Sec. III remain practical in these settings except that the topology constraint is removed, i.e., we do not enforce Nc i=1 z ci to be a constant value. Nonzero genus surfaces with boundaries For nonzero genus meshes with boundaries, the difference in treatment compared to boundary-free meshes is the same as in the genuszero case: we impose boundary conditions and remove the topological constraints. If the Dirichlet boundary conditions are applied, we similarly do not modify the Laplacian matrix.

Fig. 1 shows three models with boundaries handled by our method, where the Gauss-Bonnet constraint is no longer imposed, allowing positive and negative cones to be placed freely.

If, at each intersection, we use only a single variable instead of the three variables u 3 , u 6 , u 9 as in Fig. 9, then

with A being the Laplacian matrix and C = D = 0. In this case, rank(L g ) = N -1. However, when we introduce two variables at each intersection, two columns in C and D become nonzero, which increases the rank by 2. Since there are g intersections in total, the rank of L g becomes N + 2g -1.

c≥ N

c

the