Constructing self-supporting surfaces with planar quadrilateral elements

We present a simple yet effective method for constructing 3D self-supporting surfaces with planar quadrilateral (PQ) elements. Starting with a triangular discretization of a self-supporting surface, we first compute the principal curvatures and directions of each triangular face using a new discrete differential geometry approach, yielding more accurate results than existing methods. Then, we smooth the principal direction field to reduce the number of singularities. Next, we partition all faces into two groups in terms of principal curvature difference. For each face with small curvature difference, we compute a stretch matrix that turns the principal directions into a pair of conjugate directions. For the remaining triangular faces, we simply keep their smoothed principal directions. Finally, applying a mixed-integer programming solver to the mixed principal and conjugate direction field, we obtain a planar quadrilateral mesh. Experimental results show that our method is computationally efficient and can yield high-quality PQ meshes that well approximate the geometry of the input surfaces and maintain their self-supporting properties.


Introduction
Self-supporting surfaces, as one of the most ancient and elegant techniques for building curved shapes, can be found in many architectural heritages, because of their advantageous structural properties and efficient use of material. A 3D surface is self-supporting if the internal stresses are along the tangent directions so that they balance the surface's weight, and there is no shear stress or bending moment within the surface [1][2][3]. Designing self-supporting surfaces is an overconstrained problem since there are 3 constraints in the equilibrium equation whereas the surface has only 2 degrees of freedom, which are the isotropic stress function and the z-coordinate.
Constructing self-supporting surfaces has received considerable attention in architectural geometry [4] this century, leading to several computational tools [5][6][7][8][9]. Most existing self-supporting surface generation methods target triangle meshes, the dominant representation in digital geometry processing thanks to its flexibility and ease of construction. We note that quadrilateral meshes are highly desired in architectural geometry since most smooth surfaces have two dominant local directions, associated with principal curvature directions, and aligning quad elements with given directions is a natural way to capture shape features. Using triangle meshes, however, necessitates an arbitrary choice of a third edge direction. Moreover, triangle meshes require more edges than quad meshes, so are heavier in architectural design [4].
In this paper, we are interested in modeling self-supporting surfaces with planar quadrilateral (PQ) meshes which are highly desired for glass structures [10][11][12]. Since PQ meshes are the discrete counterpart of conjugate curve networks [13,14], computing conjugate directions plays a critical role in PQ mesh construction [4,11]. Existing methods [11,12,15,16] apply global optimization to compute conjugate directions on polygonal meshes. Though they can yield smooth PQ meshes, these methods are computationally expensive, and it is not easy to find physically-valid solutions due to high non-linearity in the optimization.
To overcome the aforementioned challenges in existing methods, we propose a simple yet effective method for constructing self-supporting surfaces with planar quadrilateral elements. Our idea is to parameterize an existing triangle-mesh-based selfsupporting surface and then extract a PQ mesh from the parameterization. Our method takes a triangular discretization of a self-supporting surface as input, and the principal curvatures and directions for each triangular face are calculated using a new discrete differential geometry approach, yielding more accurate results than existing methods. Then, the principal direction field is smoothed to reduce the number of singularities. Next, all faces are partitioned into two groups in terms of principal curvature differences. We compute a stretch matrix for faces with similar principal curvatures; it turns the principal directions into a pair of conjugate directions. For faces with large curvature difference, we keep the smoothed principal directions. Finally, a mixed-integer programming solver is applied to the mixed principal directions and conjugate directions to compute a planar quadrilateral mesh.
Since our method does not require any non-linear optimization for computing conjugate direction fields (CDFs), it is easy to implement and computationally efficient. Experimental results on various selfsupporting surfaces (including both height and non-height fields) demonstrate the effectiveness and computational efficiency of our method when generating high-quality self-supporting surface with planar quadrilateral elements. Figure 1 illustrates an example generated by our method.

Self-supporting surfaces
Self-supporting surfaces are often designed by TNA, a popular computational tool developed by Block Fig. 1 Given a self-supporting surface as a triangular mesh, our method converts it into a planar quadrilateral mesh that approximates the input surface faithfully both geometrically and in physical performance. and his colleagues [8,17]. TNA approximates the stress field through uniaxial singular stresses and factors the over-constrained 3D equilibrium problem into horizontal and vertical sub-problems, making it easier to solve. Vouga et al. [18] applied TNA to approximate freeform shapes by self-supporting ones. Since their method solves an over-constrained linear system, it cannot guarantee convergence when vertices have high valence and/or large deformation is needed to make the reference mesh self-supporting.
Miki et al. [19] generalized the Airy stress function in the planar elasticity problem to the 3D thin shell problem and developed a method for computing self-supporting spline surfaces with high-order smoothness. However, the mean curvature of the thin shell causes the Airy stress function to be inconsistent with the true mechanical equilibrium. Therefore, this method can only compute approximate solutions.
Both Liu et al. [6] and de Goes et al. [5] projected 3D surfaces onto the 2D domain to simplify the problem. De Goes et al.'s method allows the user to specify the stress distribution, giving greater flexibility. Liu et al.'s method uses regular triangulation to ensure equilibrium in the xand y-directions.
Ma et al. [9] proved that the hyper-generatrix of a 4D minimal hyper-surface of revolution is a 3D self-supporting surface. As a result, constructing a 3D self-supporting surface is equivalent to 4D volume minimization, which can be further converted into a surface reconstruction problem with mean curvature constraint [20,21]. Ma et al. showed that their method computes solutions of the over-constrained equilibrium problem more accurately than the TNAbased approaches [5,6,18]. The input of our method is a 3D self-supporting surface represented by a triangle mesh, as generated by Ma et al.'s method [9].

Conjugate direction fields
Conjugate direction fields are at the core of constructing PQ meshes, and in discrete geometry as counterparts of conjugate curve networks. The first algorithm for designing PQ meshes, proposed by Liu et al. [11], optimized the fairness of mesh edges, with face planarity constraints. Zadravec et al. [12] computed smooth CDFs by optimizing the smoothness of the representation vector field. Their method can produce high-quality quad-dominant meshes with planar faces, but it is unable to handle ±k/4 (k ∈ Z 2 ) singularities, which are common for surfaces with convex corners such as a rounded cube. Liu et al.'s method [15] allows the presence of such singularities, but at the cost of solving a non-linear optimization problem.

Overview
Our method takes a self-supporting surface represented by triangle mesh as input. We first compute the principal curvatures and directions for each triangular face using a new method for computing the curvature tensor on triangle meshes. Then we group the triangular faces whose principal curvature difference is small. We call these faces free faces and the remaining faces fixed faces. Since the principal directions on fixed charts are already conjugate, our goal is to construct a pair of conjugate directions for each free chart. Towards this goal, we compute a stretching matrix that turns the principal directions of free charts into conjugate directions. The principal directions of the fixed charts and the conjugate directions of the free charts induce a conjugate direction field on the input mesh. Applying the mixed integer quadrangulation method [22] to the CDF, we parameterize the triangle mesh and extract the planar quadrilateral elements. We present the algorithmic details in Sections 3.2-3.6. We show the algorithmic pipeline of our method in Fig. 2 and pseudocode in Algorithm 1. To ease reading, we list major notation in Table 1 Fig. 2 Algorithmic pipeline. The input is a 3D self-supporting surface represented by a triangle mesh. First, we compute principal curvatures and principal directions for each triangular face. Then, we group the faces with similar principal curvatures, which are called free faces (green). We compute a smooth cross field on the free faces and generate a pair of conjugate directions by stretching the principal directions. Next, we parameterize the triangle mesh using the mixed principal direction field and conjugate direction field. Finally, we extract the planar quadrilateral mesh from the parameterization. Table 1 Notations   Item Description

Computing principal curvatures and directions
Curvature is the amount by which a curve deviates from a straight line, or a surface deviates from planar. For a smooth 3D surface r(u, v), the curvature tensor C satisfies Cr u = −n u and Cr v = −n v , where r u , r v are the tangent vectors of S and n u , n v the partial derivatives of unit normal n. Since we deal with 3D self-supporting surfaces represented by triangle meshes, we need to discretize C. The key idea is to consider r u , r v as the movement of observation points on the surface. Consider a triangular face ABC and its neighbors CBD, ACF , and BAE. Let n, n i (i = 1, 2, 3) be the normals of ABC, CBD, ACF , and AEB, and δn i be the derivatives of these normals.
Observe that each face normal is the cross product of two edge vectors, and each vertex normal is computed as the weighted average of face normals. Therefore, face normals are usually more accurate than vertex normals on triangle meshes. Our movement strategy is to use the centers O and O i (i = 1, 2, 3) of the circumscribed circles of ABC and its neighbors BCD, ACF , and BAE. We compute the derivatives as for i = 1, 2, 3. Since the curvature tensor C satisfies Cδr i = −δn i (i =1, 2, 3), we form an over-constrained linear system: (3) where e 1 and e 2 are the axes of the local coordinate system on ABC. Then we compute x i (i = 1, 2, 3) by solving a linear least squares problem. After that, we can form the curvature tensor C:

Smoothing principal directions
The computed principal directions are usually not smooth due to discretization. As a result, the principal direction fields contain many undesired singularities, which compromises the parameterization quality. To reduce singularities, we apply an iterative method to locally adjust singularities. Specifically, the method moves singularities based on their indices: singularities of the same index repel, and of opposite index attract. When a pair of opposite-index singularities meet, they cancel each other out, reducing the number of singularities by two. After all singularities stop moving, we obtain a smoothed principal direction field with fewer singularities than the original. Figure 4 shows such a principal curvature direction field before and after optimization. We observe that the number of singularities is significantly reduced using this local smoothing method.

Distinguishing fixed and free faces
On a self-supporting surface, some parts typically exist which are spherical or nearly so (see Fig. 5(a)). It is difficult to calculate principal directions on these parts, and the direction field produced in Section 3.2 is inaccurate there. However, using the ideas of Liu et al. [11,15], principal direction fields are not the only way to make a PQ mesh: conjugate direction fields can also be used for this purpose. Thus we need to manipulate these cross vectors into conjugate directions.
Regions with large principal curvature differences can be ignored, since we can directly use the computed principal directions, as they are already conjugate. We need to pay attention to regions with small principal curvature differences (i.e., almost spherical regions), on which principal directions are ill-defined. A face is classified as a free face when its principal curvatures satisfy where ε is a userspecified threshold. In our implementation, we set ε = √ 3. See Fig. 5 for two examples showing classification of fixed and free faces.

Computing the stretch matrix S
To generate conjugate directions for free faces, we compute a matrix S, which stretches the principal directions.
Let u and v be the tangent vectors at point p on the surface. If p, p + su, p + tv, and p + su + tv form a planar quadrilateral for some scalars s and t, they satisfy the following constraint: where C is the curvature tensor at p. This implies the tangent directions u and v satisfy: so are conjugate directions. If the edges of each quadrilateral face follow conjugate directions, its vertices are thus co-planar. Therefore, conjugate directions are often used in constructing planar quadrilateral meshes [11,12,15]. The surface curvature tensor C is where L, M , and N are the coefficients of the second fundamental form, and ⊗ is tensor product. It is well known that the principal directions d 1 Since d 1 and d 2 are perpendicular, multiplying by S 1 and S 2 does not influence each direction.
Multiplying S 1 and S 2 yields the stretch matrix S: where n is the unit normal. We can thus re-write S as Matrix S is symmetric. Now we show that the matrix S turns two perpendicular tangent directions into a pair of conjugate directions. Since the principal curvatures κ i and principal directions d i satisfy the following equations: Given two arbitrary perpendicular tangent vectors u and v, we can easily verify that the stretched vectors Su and Sv are conjugate because Fig. 6 for an illustration. Therefore, the stretch matrix S turns perpendicular tangent directions into conjugate directions. In our implementation, we stretch the cross vectors on free faces to generate conjugate directions. As a local operation, computing stretch matrix S is computationally efficient. Figure 7 shows examples of conjugate directions generated by stretching principal directions.

Parameterization & generating PQ meshes
Let (u, v, −u, −v) be the directional field, consisting of principal directions on fixed faces and conjugate directions on free faces. It is used as a guide to parameterize the triangle mesh and extract PQ elements.
We need to cut the input mesh into a topological disk so that the singularities of the directional field lie on the disk boundary. To do that, we take the singularities as seeds and apply Dijkstra's algorithm to compute a Voronoi diagram on the input mesh. The dual graph of the Voronoi diagram is a triangulation with singularities as vertices. We then segment the triangle mesh into patches along the edges of the Delaunay triangulation. Then, starting with an arbitrary patch, we iteratively glue the patches together by traversing them in breadth-firstsearch order. Since each patch is used only once, the resulting shape is a topological disk with singularities on its boundary.
To compute the parameterization, we adapt the global parametrization used in Liu et al.'s work [15]. Denote by (s, t) the parameters assigned to each vertex of the input triangle mesh. We minimize the following energy function: where A i is the area of face f i and ⊥ denotes 90 • rotation about the normal. To make the parameterization seamless, for each edge on the cut graph, we require the projections of parameter gradients on the two adjacent faces to match, leading to integer translations on the parameters. Equation (7) is a nonlinear energy function with integer variables introduced at the topology cut and singularities [22]. Like in Liu et al.'s method [15], we solve the mixed integer programming using a greedy heuristic. We perform three operations in each iteration. First, we relax the integer constraints and minimize the function using the L-BFGS method [24]. Second, we pick an integer variable that is closest to integer, round it off, and fix it. We repeat the first two steps until all integer variables are fixed.

Experimental results
We implemented our algorithm in C++ and evaluated it on a laptop with a 2.90 GHz Intel Core i7 CPU and 8 GB memory.

Principal curvatures and directions
To calculate the curvature tensor C, our method uses the observation point movement strategy, based on the centers of circumscribed circles. Here we compare our method to Rusinkiewicz's method [23] in terms of accuracy and robustness on a torus model: ((R + r cos v) cos u, (R + r cos v) sin u, r sin v) T setting R = √ 10 and r = 1. The principal curvatures are κ 1 = cos v/(R + r cos v) and κ 2 = 1/r, while the principal directions are d 1 = (− sin u, cos u, 0) T and d 2 = (− sin v cos u, − sin v sin u, cos v) T . We generated two sets of triangle meshes with different degrees of anisotropy τ and resolutions.
The anisotropy degree τ for a triangle f is defined as τ (f ) = P H/(2 √ 3A), where P , H, A are the half-perimeter, longest edge length, and area of f , respectively. τ 1 for all triangles and equality holds for equilateral triangles. The numerical results show that both methods converge with increasing number of resolutions. However, our method is more accurate and robust than theirs, especially for anisotropic meshes. See Fig. 8. Our method has better accuracy than Rusinkiewicz's method [23] since we adopt per-face normals in computing the curvature tensor C, while Rusinkiewicz's method [23] uses per-vertex normals, which are more sensitive to triangulation quality than per-face normals and often result in large accumulation errors on anisotropic meshes.

Evaluation
We evaluated our method and other approaches in terms of both physical and geometrical properties. Like Ref. [9], we adopt three physical measures: normal displacement θ, stress tensor error σ S , and differential stress σ D , to evaluate the stability of selfsupporting surfaces. We also measure the smoothness of parameter lines and the planarity of the resulting PQ meshes.
We define the smoothness measure by averaging the discrete geodesic curvatures of the parameter lines (see Fig. 9): where n is the number of vertices of the quad mesh. Given a PQ mesh with m patches, we define the planarity measure γ by averaging the volume-to-base ratio: Fig. 8 Comparison with Rusinkiewicz's method [23] in terms of accuracy and convergence plots. We generate sequences of isotropic (first two rows) and anisotropic meshes (last two rows) with increasing number of faces. The heat colour map shows the angle differences δθ between the computed principal directions and the ground-truth; warmer colors are larger errors. For each model, we show our results above and theirs below. τ is the anisotropy measure and |F | is the face count.
where V i is the volume of the tetrahedron formed by the four vertices of the i-th quad patch and A i is the area of the patch. Obviously γ = 0 for a completely planar quad mesh. Therefore, the lower the value γ, the better the quality of a PQ mesh in terms of planarity.

Comparison
We compared our method to the mixed integer quadrangulation (MIQ) method [22] using the above physical and planarity measures. Table 2 shows that the MIQ method produces quad meshes with higher self-supporting performance than that of our method. However, our method outperforms MIQ in terms of planarity. This comparison indicates that a conjugate direction field is important for constructing self-supporting surfaces with planar quadrilateral elements. Without a conjugate direction field, the resulting meshes are far from PQ meshes, and so difficult to construct using glass structures.
We also compared our method with Ma et al.'s  It is worth noting that our method and Liu et al.'s method [15] adopt different strategies to deal with faces with small principal curvature difference. Their method solves a global optimization problem to find a pair of conjugate directions v p and w p for each point p that satisfy κ p,1 (v p · e p,1 )(w p · e p,1 ) + κ p,2 (v p · e p,2 )(w p · e p,2 ) = 0 where κ p,1 and κ p,2 are the principal curvatures at p, and e p,1 , e p,2 are the corresponding principal directions. Since the principal directions are not accurate in sphere-like regions, the errors may accumulate as the iterative algorithm proceeds. To improve planarity, Liu et al. [15] introduced a planarity term in the objective function. However, the side effect of this term is larger shape approximation error, which in turn may compromise the equilibrium condition of the generated shape. In contrast, our method generates conjugate directions by manipulating the principal directions using a stretch matrix. We calculate smooth conjugate directions without planarity optimization. Therefore, our method is able to well approximate the input geometry and yields better results in terms of both geometric and physical properties. Figures 10 and 11 show results for height and nonheight fields, respectively. We visualize the physical measures using color maps; warmer colors indicate larger errors or distortion. We observe that our results are comparable to Ma et al.'s method [9] that produces triangle meshes and are better than Liu et al.'s method [15] that produces PQ meshes.
We may also qualitatively compare our method to Vouga et al.'s method [7]. Their method allows the user to design self-supporting surfaces using a reference surface. To generate PQ meshes, they solve a complicated non-linear optimization problem that takes equilibrium condition, shape smoothness, and planarity into consideration. Our method is based on local computation of a conjugate direction field and then constructs PQ meshes via mixed integer based parameterization. Our method aims to preserve the geometry of the input shape which is already selfsupporting, so does not need a reference surface. Our method (including both CDF computation and global parameterization) is more efficient than theirs.

Conclusions and future work
We have developed a simple yet effective method for generating 3D self-supporting surfaces with planar quadrilateral elements. Our method takes a triangular discretization of a 3D self-supporting surface as input and generates a PQ mesh by constructing a conjugate direction field. In contrast to existing methods that compute conjugate direction fields using global optimization, our method is a lightweight approach since it computes a local stretch matrix for each triangular face with small principal curvature difference and turns principal directions into conjugate directions. Computational results We visualize the quality of the results using 3 quantitative measures: normal displacement θ, stress tensor error σ S , and differential stress σ D . The lower the measures, the higher the quality of the resulting PQ meshes. show that the generated PQ meshes approximate the input shape very well, thereby satisfying the balance constraints. Comparisons with existing methods demonstrate the effectiveness of our method.
There are a few interesting directions for future work. First, our method requires the input shape to be self-supporting. From the application point of view, it is highly desirable to either construct PQ meshes directly from user-specified boundary curves or to take a reference surface as input and then optimize its geometry to fulfil the physical requirements of being self-supported. To make the method more flexible, we will also take the userspecified soft and/or hard directional constraints into consideration.