Dev2PQ: Planar Quadrilateral Strip Remeshing of Developable Surfaces

Developable surfaces are commonly used in architecture and product design due to the simplicity of their fabrication. Such surfaces are locally isometric to a planar domain, which means they can be manufactured by mere bending of sheet material, such as metal. Freeform developable surfaces form a rich and interesting shape space, but they are notoriously difficult to design due to their highly constrained nature. Therefore, in most cases, only simple forms are used, such as cylinders and cones.

The majority of methods for developable surface modeling use rulings-based representations (see, e.g., Tang et al. 2016 and Solomon et al. 2012) or isometry optimization (see, e.g., Burgoon et al. 2006 and Fröhlich and Botsch 2011). The recently proposed discrete orthogonal geodesic nets Rabinovich et al. 2018 and the checkerboard pattern isometries Jiang et al. 2020 represent developable surfaces without explicitly accounting for principal curvature directions or patch decomposition. A more recent approach to creating developable surfaces is through the approximation of existing shapes with developables. Proposed methods include cutting surfaces into developable pieces Sharp and Crane 2018, fitting developable patches to a surface Ion et al. 2020, and flowing existing surfaces to reduce a discrete developability measure Binninger et al. 2021 Kohlbrenner et al. 2021 Liu and Jacobson 2021 Sellán et al. 2020 Stein et al. 2018.

While these discrete representations of developable surfaces are excellently suited for creative exploration and design of freeform developable shapes, they fall short of providing a suitable final representation for manufacturing. For that purpose, it is especially important to have planar mesh faces aligned to principal curvature directions Alliez et al. 2003 Liu et al. 2006 Tang et al. 2016, that only the ruling-based representations provide. A curvature-line representation of a developable surface possesses the desired properties for fabrication once the shape is fixed, since the minimal curvature lines on a developable surface are rulings on which the normal is constant; as such, they can be easily tessellated into planar polygons that approximate the surface shape well. In fact, meshes composed of planar quadrilateral strips (with no interior vertices) constitute a well-known model for discrete developable surfaces, whose refinement and convergence properties have been studied Liu et al. 2006. Additionally, the planar quadrilateral strip representation allows for the shape to be fabricated by bending developable material, such as sheet metal, along the ruling edges of the mesh (as shown in Figure 22). General planar quadrilateral meshes that do not consist of strips, i.e., not aligned to minimal curvature, are not suitable for this fabrication process.

In this article, we develop a method to convert a triangle mesh representation of a developable surface into a discrete curvature-line representation in order to reap the benefits of both worlds: the support for developable shape creation provided by a representation of choice and the desirable properties for fabrication offered by the curvature-line representation. Our method produces strips of planar quadrilaterals (PQs) aligned to principal curvature directions that closely approximate the curved parts of a given input mesh, along with planar polygons representing the flat parts of the input (see Figure ). In particular, our method produces precisely straight rulings, modeled as individual edges in the output mesh.

The past decade has seen a highly active stream of fruitful research on field-aligned quad meshing, where principal curvature fields have naturally received special attention Bommes et al. 2012 Vaxman et al. 2016 de Goes et al. 2015. However, to the best of our knowledge, no existing general remeshing method is guaranteed to produce completely straight-edge sequences that are consistent with curvature lines, which are often difficult to obtain fully and faithfully for discrete meshes. In this work we exploit this specific constrained setting and the geometric structure of developable shapes to reproduce straight minimum-curvature lines, as well as automatically segment the input into curved and planar parts in a robust manner that is consistent with the structure dictated by developability.

Our method is based on fitting a scalar function on the input mesh, such that its isolines are straight and align as best as possible to the locally estimated rulings on non-planar regions (see Figure ). The condition that guarantees straight isolines on developable surfaces is simple to formulate: the normalized gradient of the scalar field needs to be divergence free. This non-linear and high-order condition is numerically difficult to enforce in a straightforward manner. We therefore devise a dedicated optimization scheme that factors the problem into a divergence-free and integrable directional field optimization that is subsequently integrated into a scalar function. This makes our problem tractable and practical to optimize. We extract the isolines of the obtained scalar field at the desired resolution and remesh the input into strips of planar quads whose chordal edges are the isolines, i.e., the rulings. We supplement the mesh by planar polygonal faces that represent the planar patches of the input surface. The flexibility of the field-to-function design allows for the automatic placement of singularities, flat regions, and curved folds without explicitly segmenting different curvature regions on the mesh.

We demonstrate the effectiveness of our approach on a variety of input developable shapes represented by general triangle meshes (and triangulated quad meshes) and show how our remeshing can be used side by side with freeform modeling of developables.

Remeshing general meshes into (planar) quad meshes is an active area of research. A comprehensive review is beyond the scope of this article, but we highlight the main features of existing approaches most closely related to our work.

As stated in the introduction, the quadrilaterals in freeform models of developable surfaces are usually non-planar, and typically neither triangle nor quad meshes are curvature aligned. Our goal is to obtain a curvature-aligned remeshing with planar faces. Planarization of general polygonal meshes has been explored in several works Alexa and Wardetzky 2011 Bouaziz et al. 2012 Tang et al. 2014 Diamanti et al. 2014 Poranne et al. 2013. These methods take arbitrary shapes as input and are not specifically targeted at developable surfaces. Typically, applying a general planarization method to developable surfaces leads to poor results in terms of curvature alignment and shape approximation (see Figure 2).

A different approach to obtaining PQ meshes from general developable input meshes is to utilize the fact that PQ meshes are a discrete model for conjugate nets and seek a remeshing that is aligned to ruling directions. Many curvature-aligned or just conjugate quad remeshing techniques for general shapes exist; see e.g. Bommes et al. 2012, Jakob et al. 2015, Diamanti et al. 2014, Liu et al. 2011, and Zadravec et al. 2010. Similar to our method, these techniques rely on numerical estimation of the principal curvature directions, but they do not guarantee exact alignment or straight edge sequences and may introduce unnecessary singularities on developable shapes. Their optimization process might fail to create precise, straight rulings on developable surfaces, unlike the algorithm we propose in this work (see Figure 4). From a fabrication viewpoint, a general PQ remeshing method is lacking since it merely produces a PQ mesh, rather than a PQ strip mesh, which allows for simpler construction.

A more promising approach to PQ meshing of developable shapes is a dedicated technique that utilizes their specific properties. Peternell [2004] converts a scan of a single torsal developable patch into a PQ mesh by thinning its tangent-space representation into a one-dimensional object (a simple curve). This approach is not immediately applicable to composite and possibly piecewise developable surfaces that consist of multiple torsal patches and planar regions. Kilian and colleagues [2008] compute a torsal patch decomposition for 3D scans of physical developable surfaces by estimating flat regions and ruling directions. This approach may struggle with developable meshes that are coarse in comparison to scans due to insufficient data density for reliable ruling fitting. Their method relies on the ruling estimates to initialize a planar mesh development, which is used in the subsequent optimization. The connectivity of this initial mesh cannot be changed during the optimization, and thus determines the approximation quality that can be obtained. Locally estimated rulings on developable meshes can be quite noisy and inaccurate, as we discuss in Section 4. We avoid a direct domain decomposition based on rulings employed in Kilian et al. 2008 and instead devise a global constrained optimization approach. As a result, our method is successful on coarse and noisy inputs. Additionally, in contrast to Kilian et al. 2008, our method can handle piecewise developable surfaces that are not curved folds constructed from a single sheet of material. Their method optimizes and deforms a surface to become developable as a whole, whereas our work takes a piecewise approach to remeshing the existing surface without deforming it. We compare our results to their work in Section 6.

Wang and colleagues [2019] use discrete parallel geodesic nets as a discrete model for developable surfaces. They require the geodesic strips to be of constant width for a surface to be developable but do not impose any requirements on the directions of these strips and as such do not have a ruling-aligned representation for developables. They also use parallel geodesic nets to approximate surfaces by piecewise developable strips. The individual geodesic strips of a parallel geodesic net are approximated by piecewise planar faces in a post-processing step, but this does not guarantee that compatibility between neighboring strips is preserved. In contrast, our method produces a complete, connected remeshing of the input developable surface with strips aligned to the rulings, rather than in the orthogonal direction.

While targeting geodesic fields, rather than planar quad remeshing, the works by Vekhter et al. [2019] and Pottmann et al. [2010] show parallels to our proposed method. They compute a unit curl-free field, while we compute a divergence-free field—these two kinds of fields are in fact duals. Nevertheless, our field has further constraints in terms of ruling alignment, which we take into account. In addition, our optimization strategy is different, interlacing integrability optimization with divergence reduction. We discuss this in further detail in Section 3.

We next summarize relevant facts about both continuous and discrete developable surfaces that inform our algorithm and offer a directional-field-based definition of developable surfaces. We provide a discrete setup for these fields in Section 4 and an optimization scheme in Section 5.

The input to our algorithm is a triangle mesh \(\mathcal {M} = \lbrace \mathcal {V},\mathcal {E},\mathcal {F}\rbrace\) representing the (piecewise) developable surface, where \(\mathcal {V}\) denotes the set of vertices, \(\mathcal {E}\) the set of edges, and \(\mathcal {F}\) the faces. To regularize the scale of surface curvature between different surfaces and our optimization parameters, we scale the input \(\mathcal {M}\) to have unit-length bounding box diagonal. We define \(u\) as a piecewise-linear vertex-based function \(u(v), \ v \in \mathcal {V}\), and consequently represent \(r\), \(r^\perp\), and \(\nabla u\) as face-based piecewise-constant tangent fields; we denote this space as \(\mathcal {X}\). We use the conforming discrete gradient \(G: (\mathcal {V}\rightarrow {\rm {I}\!R})\rightarrow \mathcal {X}\) and divergence \(D:\mathcal {X}\rightarrow (\mathcal {V}\rightarrow {\rm {I}\!R})\) operators, and the non-conforming discrete curl operator \(C:\mathcal {X}\rightarrow (\mathcal {E}\rightarrow {\rm {I}\!R})\). Their explicit expressions can be found in, e.g., Brandt et al. 2017.

Estimating rulings. We compute a ruling direction \(r(f),\ \forall f \in \mathcal {F}\), as the eigenvector corresponding to the minimal eigenvalue of the face-based shape operator \({S}(f)\), as defined in de Goes et al. 2020. Since we know the ruling only up to sign, we represent it unambiguously using a power representation Azencot et al. 2017 Knöppel et al. 2013: we first represent \(r(f)\) as a complex number in a local coordinate system and then square this complex number to have a representation that is invariant to the sign of the direction; i.e., we store \(R(f) =r^2(f)\). We also define \(R^{\perp }(f) = (r^{\perp }(f))^2\), the power representation of the ruling locally rotated by 90 degrees.

Confidence weights. A clean domain decomposition into planar and torsal regions would significantly simplify the fitting of individual developable patches. Unfortunately, we cannot obtain such a clean segmentation directly, because the curvature measure (like the ruling estimates) is noisy and does not delineate planar and torsal patches nicely (Figure 5). Therefore, we model on the assumption that the rulings are least reliable in planar or near-planar regions, and mostly consistent in strongly curved areas (see Figures and 6). We thus attach a relative confidence weight \(w(f)\) to each face \(f\in \mathcal {F}\), as a function of the discrete absolute max and min curvatures \(\kappa _1(f)\) and \(\kappa _2(f)\):For \(\kappa _1(f)\) and \(\kappa _2(f)\) we use the absolute largest and smallest eigenvalues of the shape operator \(\mathcal {S}(f)\) and set \(\theta _1=0.8\), \(\theta _2=-0.014\).

The confidence function is a logistic curve (see inset), facilitating a stronger distinction between confidence in planar and near-planar regions (albeit still small compared to stronger curved regions). The value of \(w(f)\) is capped at 0.8 by design, ensuring that we never fully rely on a ruling. We define \(\mathcal {V}_b\) to be all vertices on the boundary of \(\mathcal {M}\) and \(\mathcal {F}_b\) all faces that contain a vertex in \(\mathcal {V}_b\), and we set \(w(f)=0\) for these faces.

Creases. Our method requires as input the explicit identification of the set of crease edges \(\mathcal {E}_c\) that define curved-fold creases and boundaries to developable pieces. We define \(\mathcal {V}_c\) as all vertices that are incident on an edge in \(\mathcal {E}_c\), and from this we define the set of faces adjacent to them: \(\mathcal {F}_c\) is the set of all faces that have one or more vertices in \(\mathcal {V}_c\). We update the confidence weights by setting \(w(f)=0\) for all faces in \(\mathcal {F}_c\). Using \(\mathcal {V}_c\) and \(\mathcal {V}_b\), we initialize \(\mathcal {V}^{\ast }\) as the set of mesh vertices but with the boundary and crease vertices excluded, i.e. : \(\mathcal {V}^{\ast } = \mathcal {V} \setminus (\mathcal {V}_b\cup \mathcal {V}_c)\). The user can provide crease edges as input or we can make an initial guess for the crease edges based on the dihedral angle of adjacent faces and manually add edges belonging to softer creases. We require that the final set of crease edges divides the surface into smooth developable surfaces. The required seams for this segmentation are typically easily visually distinguishable (e.g., using reflection lines). Figures 10, 17, 19, 20, and 21 show examples of developable surfaces with creases, and Figure 19 explicitly shows them on a complex model. When creases end in the interior of a developable piece rather than on another crease or boundary, so-called open creases, we duplicate their interior vertices and define them as mesh boundaries. Examples of such open creases can be seen in Figures 10, 18, and 19 (the chair model).

We describe our approach to remeshing a (near-)developable input triangle mesh to a curvature-aligned, planar polygonal mesh consisting primarily of PQ strips.

We implemented our algorithm using libigl Jacobson et al. 2018 and Directional Vaxman et al. 2017a on a machine with i7-8569U CPU and 16 GB RAM. Our typical input mesh resolution is 1,800 faces, and for this approximate input size the vector field design part of our method takes 4–5 seconds, of which the majority of the time is spent in the ProjectCurlFree step, i.e., solving the convex optimization in Equations (22)–(24). We currently use CVX Grant and Boyd 2014, but this part can be optimized for better speed. Although we do not have a formal convergence guarantee for our alternating algorithm, we observe that it typically converges to our specified tolerance level within 10–20 iterations for vector fields without singularities and 40–50 iterations for shapes with planar parts that introduce singularities in the vector field. We also test our method on inputs of up to 160k faces, which does not cause problems for convergence. The parameterization part of our method takes approximately 10–15 seconds.

A variety of our results can be seen in Figure 23. Note that our method preserves the input boundary vertices, and therefore our output faces are quadrilateral-like higher-degree polygons, rather than actual quadrilaterals. Examples of our results with various boundary shapes and non-disk topologies are included in Figure 23. Our method is applicable to developables with curved folds, as seen in Figures, 1, and 21 (input models from Rabinovich et al. 2019) as well as in Figure 10. It can handle piecewise developable shapes, such as D-forms (shapes obtained by gluing together two planar domains with the same perimeter) from Jiang et al. 2020 and sphericons from Tang et al. 2016 (see Figure 17), as well as other shapes with creases from Tang et al. 2016 (see Figure 20). Works by Stein et al. [2018], Sellán et al. [2020], Ion et al. [2020], and others produce piecewise developable approximations, but they are not PQ meshes. Our work can be used to remesh those surfaces; see Figure 19 for an example. We successfully apply our method to glued constructions from Jiang et al. 2020, including point singularities; see Figure 18. We have physically fabricated some of our results, shown in Figure 22. Table 3 lists the most important statistics about our results.

Developable surface editing with dynamic connectivity. To demonstrate the utility of our approach, we use the point handle-based editing system of Rabinovich et al. 2018 to interactively deform an input discrete orthogonal geodesic net (DOG) and create a sequence of a few developable surfaces, on which we run our algorithm (offline) after trivial triangulation. See Figure 7 and the accompanying video for some examples of such editing sessions. Note the natural change in the combinatorics that our algorithm induces to model exact developability, which can change considerably even for small deformations in the input.

Planarity evaluation. Since our output meshes have no interior vertices inside the developable patches, the ultimate accuracy measure for the developability of our results is the planarity of the mesh faces. We measure planarity of each quadrilateral face by the ratio of the distance between the diagonals to their average length, in percent Liu et al. 2006. For higher-degree polygons, we compute the root-mean-square (RMS) error of all quads constructed from every four consecutive vertices in the polygon. An acceptable stringent tolerance for the planarity error is \(\le 1\%\) Vaxman et al. 2017b. It is generally not expected for parameterization-based methods to achieve planarity to more than first order, so that usually further planarization post-processing is needed. We show the raw maximum and mean planarity error values of our results without any post-processing in Table 3. Even though our output meshes are quite coarse, our planarity errors are typically very low without such planarity optimization, and very close to the tolerance. Our worst maximum planarity error is obtained on a mesh with very thin features (Figure 23, fifth row, middle column), and the output quad with this maximal planarity error is toward the end of the spiral. As can be seen in Figure 23 and Figure 9, the input triangulation is very coarse there. We planarize this example, which has the worst maximum planarity error (\(p=11.84\%\)), to zero planarity (\(p=0.0034\%\)) using ShapeUp Bouaziz et al. 2012 and reach a visually highly similar result; see Figure 8. This demonstrates the capability of our algorithm to utilize the information in the original mesh effectively. Interestingly, the triangulation of the input mesh is similar to a Schwartz lantern, and remeshing this input to a more uniform triangulation already drastically brings down the maximal planarity error of our algorithm to \(p=2.58\%\); see Figure 9.

Constraints and objective values. As described in Section 5.2, we have an iterative multi-step optimization algorithm. Because of its alternating nature, we can only guarantee that the converged solution will meet the curl-free constraint that is optimized for in the final step. We see that the divergence of the final (normalized) vector field is not fully zero everywhere, but it is nevertheless very low (see Table 1). A plot of the divergence values over the mesh shows that it mostly concentrates on planar regions (see inset).

As the divergence is close to zero in torsal regions, the resulting function isolines are very straight and similar to their simplified version there. The simplified isolines in planar regions possibly deviate more, but since any straight line in a planar region is a ruling, this does not pose a problem. Our final converged vector field may be far from unit-norm, but this does not raise any issues as long as the divergence of its normalized version is zero.

Comparison with state of the art. In Figure 10 we compare our results with the work by Kilian et al. [2008] through a visual comparison between (triangulated versions of) their output meshes and our outputs obtained from an isotropic remeshing of their output as our input. Because we require the input surface to already be developable, our method does not work directly on the surface scans that they use as input. We generated our input by explicitly preserving the crease edges and isotropically remeshing the patches in between. For the triangulated output meshes of Kilian et al. we determined a lower bound on the maximal and mean planarity errors by taking the global maximum and mean of the minimal planarity values attainable by pairing each face with each of its neighbors. We obtain visually similar results and also obtain comparable planarity error measures, while requiring less user input and manual tweaking. Specifically, we do not require the user to ensure that the initial mesh connectivity is a valid developable decomposition.

Effect of output resolution. We extract \(u\) isolines of varying isovalues to create output meshes of different resolutions; we then measure their planarity and approximation quality w.r.t. the input mesh in terms of Hausdorff distance (Figure 11). We note that the approximation quality and the planarity improve with higher resolution, although even for the coarsest resolution these metrics are already below tolerance.

Comparison with analytical principal-curvature directions. We test our method on an input mesh sampled from an analytical clothoid surface with varying resolution and compare the obtained vector field \(Y\) with the analytical max curvature directions. See Figure 12 and Table 2. Note that the input to our method is numerically estimated ruling directions \(r\), not their analytical values. As the data shows, upon refinement of the input mesh, our output field converges toward the analytical solution.

Robustness. Our method is robust with respect to the parameters \(\omega _a\) and \(\omega _s\), for which there is a range of values that leads to visually very similar results. As the relative weight of \(\omega _s\) with respect to \(\omega _a\) increases, the vector field turns into a more constant field, reducing alignment quality of the final output mesh. Because the smoothing step is the last one before the vector field is projected onto the divergence-free solution space, \(\omega _s\) is the dominating parameter in determining the rough layout of the final vector field. The parameter \(\omega _a\) (within reasonable range) has a smaller contribution but can affect the exact placement of singularities and the iterations required to reach convergence (see Figure 13). For noisy inputs, as in Figure 14, our method does not converge with our standard parameter settings, or it converges but generates a vector field with a large amount of singularities. For these cases, simply increasing \(\omega _s\) ensures that the optimization converges, although some small and noisy details may be lost (in Figure 14 (right) we use \(\omega _s = 0.15\)). This shows that our method with the help of the smoothness term manages to recover a principally aligned vector field even if the information from the input ruling directions is very weak. For optimal alignment the value of \(\omega _s\) should be chosen as small as possible but as high as necessary so to not introduce excessive singularities; e.g., for the cone in the second-to-last row of Figure 23 we use \(\omega _s = 0.00005\) to emphasize better alignment near the boundary.

Limitations. Our method is not entirely triangulation independent, as shown in Figure 15. If the input meshing does not allow for an accurate estimation of the principal curvature directions, this might lead to poor ruling estimation and diminished performance of our algorithm in terms of planarity error. This is most noticeable near the corners of the given input, where there is relatively little data for our algorithm to align to. In order to minimize bias introduced by the triangulation, it is advisable to triangulate polygonal input meshes with higher valences by inserting a new vertex at the face center and connecting it to the vertices of the original polygon in a triangle fan.

Furthermore, as we treat curved folds in an identical manner to creases (namely we assume they delineate a developable surface piece), there is no guarantee that the curved folded surface as a whole remains developable in one piece. An interesting direction for future work would be to incorporate known geometrical constraints at curved folds into our method.

If a cone apex is present in the mesh (see Section 3.1), the natural behavior encouraged by our algorithm is to place a singularity on the crease to compensate for the curvature of the seam (see the cylinder example in Figure 20 and Figure 16, bottom). Nevertheless, our algorithm may fail to put the singularity exactly on the seam, depending on discretization. As a result, our optimization might get stuck, oscillating between nearby solutions with different singularity configurations (Figure 16). Even though our method nominally fails to converge for the car example in Figure 16, a reasonable result close to the expected one is still obtained.

Our method may struggle with thin features, e.g., as part of a piecewise developable, as these can often provide no alignment information at all. In the future it would be interesting to see how the vector field on surrounding developable pieces can be used to add constraints to these thin features, since in the final meshing we wish to guarantee continuity throughout the pieces. As shown in Figure 12, our output quality with respect to Hausdorff distance, as well as mean and maximal planarity error, increases as the input resolution increases. Our method is therefore dependent on the input resolution but still performs well on low-resolution inputs.

Finally, we have no theoretical guarantees that the ruling-aligned edges in our output mesh do not intersect, although we never see this happen in our experiments. In torsal regions, the guiding field discourages overlapping behavior, as rulings on a developable are ordered. For planar regions, the guiding field contains more noise, but here the smoothness requirement for the vector field (and thus the corresponding parameterization) strongly discourages crossovers.

We presented an algorithm that converts developable surfaces represented by triangle meshes to a discrete curvature line parameterization, i.e., polygonal meshes with planar faces, where all interior edges correspond to rulings. This conversion from an unstructured triangulation of a developable surface to a curvature-aligned PQ mesh is an important step in the developable modeling and fabrication pipeline, for which thus far a robust practical solution was missing. Our method can only be expected to work well for nearly developable input shapes. We show a non-developable example in the inset, where the edge simplification step is omitted as it would degenerate the mesh. It would be interesting to see how we can adapt our remeshing algorithm to be usable for the approximation of non-developable surfaces by developable patches.

Another venue for further study would be the automatic tuning of the values \(\omega _s\) and \(\omega _a\) based on the noise levels of the estimated input rulings. Finally, it is conceivable to adapt the output mesh resolution based on the local curvature, allowing a denser representation in more curved areas.

The authors are grateful to Helmut Pottmann, Michael Rabinovich and Alexander Sorkine-Hornung for illuminating discussions and guidance, and to Kaan Baki for help with video production. We also thank the authors of Kilian et al. 2008, Tang et al. 2016, Stein et al. 2018 and Jiang et al. 2020 for providing us with models from their works.