b/Surf: Interactive Bézier Splines on Surfaces
Claudio Mancinelli, Giacomo Nazzaro, Fabio Pellacini, Enrico Puppo
bb/Surf: Interactive Bézier Splines on Surfaces
CLAUDIO MANCINELLI ∗ , University of Genoa
GIACOMO NAZZARO ∗ , Sapienza University of Rome
FABIO PELLACINI,
Sapienza University of Rome
ENRICO PUPPO,
University of Genoa
Fig. 1. We propose algorithms to interactively edit Bézier splines on large meshes, including curve editing, curve transformations and import and editingcomplex SVG drawings. All computations occur in the intrinsic geodesic metric of the surface. All splines in this figure have been drawn interactively. Controlpoints and tangents of curves under editing are shown in the zoomed insets. Asian Dragon ∼ ∼ Bézier curves provide the basic building blocks of graphic design in 2D. Inthis paper, we port Bézier curves to manifolds. We support the interactivedrawing and editing Bézier splines on manifold meshes with millions oftriangles, by relying on just repeated manifold averages. We show that directextensions of the De Casteljau and Bernstein evaluation algorithms to themanifold setting are fragile, and prone to discontinuities when control poly-gons become large. Conversely, approaches based on subdivision are robustand can be implemented efficiently. We define Bézier curves on manifolds, byextending both the recursive De Casteljau bisection and a new open-uniformLane-Riesenfeld subdivision scheme, which provide curves with differentdegrees of smoothness. For both schemes, we present algorithms for curvetracing, point evaluation, and point insertion. We test our algorithms forrobustness and performance on all watertight, manifold, models from theThingi10k repository, without any pre-processing and with random controlpoints. For interactive editing, we port all the basic user interface inter-actions found in 2D tools directly to the mesh. We also support mappingcomplex SVG drawings to the mesh and their interactive editing.CCS Concepts: •
Computing methodologies → Mesh models ; Paramet-ric curve and surface models ; Graphics systems and interfaces .Additional Key Words and Phrases: geometric meshes, spline curves, userinterfaces, geometry processing ∗ Joint first authors
Vector graphics in 2D is consolidated since decades, as is supportedin many design applications, such as Adobe Illustrator [Adobe Inc.2019], and languages, like Scalable Vector Graphics (SVG) [W3C2010]. Bézier curves are the building blocks of most vector graphicspackages, since most other primitives can be converted into
Béziersplines (chains of Bézier curves) and edited as such [Farin 2001].In many design applications, it would be beneficial to edit vectorgraphics directly on surfaces, instead of relying on parametrizationor projections that have inherent distortions [Poerner et al. 2018].Yet, bringing vector graphics to surfaces is all but trivial, sincebasic rules of Euclidean geometry do not hold under the geodesicmetric on manifolds, and distances, shortest and straightest pathscannot be computed in closed form. In particular, in spite of severalattempts to defining curves under the geodesic metric, a completecomputational framework that supports their practical usage in aninteractive design setting is still missing.In this work, we port Bézier curves on surfaces. We begin byanalyzing the extensions of Bézier curves to the manifold setting,which are obtained by replacing straight lines with geodesic lines,and affine averages with the Riemannian center of mass. We showthat direct approaches to curve evaluation are fragile, and maygenerate discontinuous curves. On the contrary, approaches basedon subdivision and repeated averages are robust. In particular, thescheme based on recursive De Casteljau bisection lends itself to a r X i v : . [ c s . G R ] F e b anicinelli et al. Fig. 2. The GUI of our system. Curves on the eye and on the arm consist ofa single cubic segment each, while the two splines on the ears consist of twocubic segments each, with a sharp corner to the left and smooth junction tothe right. A control tangent (in blue) is depicted on the curve under editing. an efficient implementation. Inspired by recent theoretical results[Duchamp et al. 2018; Dyn et al. 2019], we propose a new schemebased on an open-uniform Lane-Riesenfeld subdivision, which war-rants higher smoothness, and can also be implemented efficiently.Focusing on such two schemes, we design algorithms for curvetracing, point evaluation, and point insertion. Our implementationrests on light data structures and just a few basic tools for geo-desic computations. The proposed algorithms remain interactiveon meshes made of millions of triangles such as the ones shown inFigure 1. To assess the robustness and performance of these algo-rithms, we trace curves on the more than five-thousands watertight,manifold, meshes of the Thingi10k repository [Zhou and Jacobson2016] with many randomly generated control polygons. We showthat our algorithms handle well all cases.We integrate curve tracing and point insertion algorithms ina user interface prototype that supports the interactive design ofBèzier splines on manifold meshes. We support all basic operationsthat 2D editors have, including: click-and-drag of control points andtangents; point insertion and deletion; and translation, rotation andscaling of curves. We also support mapping of 2D drawings ontothe surface. Figure 2 shows the interface of our system with somesimple curves traced on a model.All together, this work advances the state of the art with fivemain contributions: • We present a critical analysis of several definitions of Béziercurves in the manifold setting. We demonstrate the limitationsof direct methods, and show that subdivision schemes arealways well behaved. • We provide algorithms for curve tracing, point evaluation,and point insertion, for a recursive De Casteljau scheme, anda novel open-uniform Lane-Riesenfeld scheme that warrantshigher smoothness than previous definitions.
Fig. 3. Curves that wind about the object or require large control polygonsmay be challenging to draw with an approach based on parametrization. Thecollar and the curl consist each of a single cubic segment, while the spiralis a spline of four segments joined with smooth ( 𝐶 ) continuity. Controlpolygons are depicted in blue. • We provide a very fast implementation of such algorithmsthat runs at interactive rates for meshes of millions of trian-gles on single-core commodity CPUs. Upon publication, wewill release all source code with a permissive license. • We show that our algorithms are robust and interactive witha large test on Thingi10k models [Zhou and Jacobson 2016]. • We develop a prototype system that supports the interac-tive design of Bézier splines, supporting all operations com-monly found in 2D editors. Our editor remains interactivewith meshes of millions of triangles. We also show how toport existing 2D drawings onto surfaces.
The design of spline curves on manifolds has been addressed byseveral authors, both from a mathematical and from a computationalperspective. We review only methods addressing general surfaces.A traditional approach to circumvent the problems of the Rie-mannian metric consists of linearizing the manifold domain viaparametrization, designing curves in the parametric plane, and map-ping the result to the surface. Parametrization introduces seams,and drawing lines across them becomes problematic. Moreover, dis-tortions induced by parametrizations are hard to predict and control.The exponential map can provide a local parametrization on the flyfor the region of interest [Biermann et al. 2002; Herholz and Alexa2019; Schmidt 2013; Schmidt et al. 2006; Sun et al. 2013]. However,its radius of injectivity is small in regions of high curvature, whilecontrol polygons and curves may extend over large regions. Evencurves as simple as the ones depicted in Figure 3 may be hard tocontrol using either local or global parametrizations.Another approach consists of relaxing the manifold constraint,resolving the problem in a space that admits computations in closedform, and projecting the result back to the surface. [Wallner andPottmann 2006] computes weighted averages in 3D space and projectsthem to the nearest point on surface. [Panozzo et al. 2013] uses an /Surf: Interactive Bézier Splines on Surfaces embedding in a higher-dimensional Euclidean space, followed byPhong projection. These methods may support user interaction, butthey provide only approximate results, are prone to artifacts, andare hard to scale to large meshes.The design of curves can also be addressed as an optimizationproblem in a variational setting. [Noakes et al. 1989] and [Camarinhaet al. 1995] provide the basic variational theory of splines on man-ifolds. This approach is adopted in several other papers [Arnouldet al. 2015; Gousenbourger et al. 2018, 2014; Hofer and Pottmann2004; Jin et al. 2019; Pottmann and Hofer 2005; Samir et al. 2011].While most such works do not address implementation and perfor-mance, [Hofer and Pottmann 2004] and [Jin et al. 2019] eventuallyresort to projection methods. Overall, the variational approach istoo computationally expensive to support user interaction on largemeshes. Moreover, these curves are harder to control interactivelythan traditional Bézier splines.Concerning the specific case of Bézier curves, [Park and Ravani1995] first extended the De Casteljau algorithm to Riemannian man-ifolds, expressing geodesic lines of control polygons through theexponential map, without further developing the computationaldetails. Later on, the De Casteljau algorithm on surfaces has beenexplored by several other authors [Gousenbourger et al. 2018; Linand Walker 2001; Morera et al. 2008; Nava-Yazdani and Polthier2013; Popiel and Noakes 2007]. Among these, [Morera et al. 2008]extends the recursive De Casteljau bisection, and [Sharp and Crane2020] achieves performance on the same algorithm, by using a fastmethod for evaluating locally shortest geodesic paths. We adopt thesame outer structure of the algorithm of [Morera et al. 2008] forcurve tracing with the recursive De Casteljau bisection. [Absil et al.2016] defines Bézier curves both with the De Casteljau algorithmand with the Riemannian center of mass, and show that they mayproduce different results.Several authors have investigated the theoretical aspects of thesubdivision approach to splines in the manifold setting. [Noakes1998] proves that the recursive De Casteljau bisection converges andproduces a 𝐶 curve, while [Wallner 2006; Wallner and Dyn 2005]prove the convergence and smoothness of some other subdivisionschemes. Most recent results [Duchamp et al. 2018; Dyn et al. 2019;Dyn and Sharon 2017] focus on Lane-Riesenfeld schemes and showthat a scheme of order 𝑘 is convergent and 𝐶 𝑘 in the manifold andfunctional settings. These latter works motivate our approach tothe open-uniform Lane-Riesenfeld subdivision. In this section, we consider different constructions for Bézier curvesin the Euclidean setting, all of which produce the same curves.For each construction, we analyze the possible extensions to themanifold setting, showing that they may lead to different curves,some of which may be ill-defined.Here, we only provide the basics of each construction, and referthe readers to [Farin 2001; Salomon 2006] for further details inthe Euclidean setting. We assume that readers are familiar withgeodesic lines and shortest paths on manifolds, and provide a briefintroduction on the subject in Appendix A.
In the Euclidean setting, a Bézier curve is a polynomial parametricfunction of degree 𝑘 b 𝑘 : [ , ] −→ R 𝑑 , which is defined by means of a control polygon Π 𝑘 = ( 𝑃 , . . . , 𝑃 𝑘 ) ,where all 𝑃 𝑖 ∈ R 𝑑 . Curve b 𝑘 interpolates points 𝑃 and 𝑃 𝑘 , and it istangent to Π 𝑘 at them. When there is no ambiguity, we will omitthe subscript 𝑘 and will denote the control polygon simply by Π .All constructions of Bézier curves in the Euclidean setting relyon the computation of affine averages of points of the form ¯ 𝑃 = ℎ ∑︁ 𝑖 = 𝑤 𝑖 𝑃 𝑖 (1)where the 𝑤 𝑖 are non-negative weights satisfying the partition ofunity. For ℎ = , the affine average reduces to the linear interpolation ¯ 𝑃 = ( − 𝑤 ) 𝑃 + 𝑤𝑄. (2)By analogy with the Euclidean setting, a control polygon Π 𝑘 in themanifold setting consists of a polyline of shortest geodesic paths,connecting the control points that lie on a complete manifold 𝑀 .Affine averages are not available on manifolds, but they can besubstituted with the Riemannian center of mass as introduced by[Grove and Karcher 1973; Karcher 1977]. Given points 𝑃 , . . . , 𝑃 ℎ ∈ 𝑀 and weights 𝑤 , . . . , 𝑤 ℎ , as before, their Riemanninan Center ofMass (RCM) on 𝑀 is given by 𝑅𝐶𝑀 ( 𝑃 , . . . , 𝑃 ℎ ; 𝑤 , . . . , 𝑤 ℎ ) = argmin 𝑃 ∈ 𝑀 ℎ ∑︁ 𝑖 = 𝑤 𝑖 𝑑 ( 𝑃, 𝑃 𝑖 ) (3)where 𝑑 (· , ·) is the geodesic distance on 𝑀 . If 𝑀 is a Euclidean space,then the solution to Eq. 3 is the usual affine average of Eq. 1.The RCM requires that Eq. 3 has a unique minimizer. [Karcher1977] provides a condition of existence and uniqueness of the solu-tion, which requires all points 𝑃 𝑖 to be contained inside a stronglyconvex ball, whose maximum radius depends on the curvature of 𝑀 . In the following, we will refer to this condition as the Karchercondition . If such condition is satisfied, then the RCM is 𝐶 ∞ inboth the 𝑃 𝑖 ’s and the 𝑤 𝑖 ’s [Afsari 2009]. Unfortunately, the Karchercondition restricts the applicability of the RCM to relatively smallneighborhoods in the general case.For any two points 𝑃, 𝑄 ∈ 𝑀 , which are connected with a unique shortest path 𝛾 𝑃,𝑄 , with 𝛾 ( ) = 𝑃 and 𝛾 ( ) = 𝑄 , it is easy to showthat their RCM with weights ( − 𝑤 ) and 𝑤 , respectively, is alwaysdefined and lies at 𝛾 𝑃,𝑄 ( 𝑤 ) . This means that weighted averagesbetween pairs of points are defined and smooth, as long as suchpoints stay away from each other’s cut locus ∗ . We extend this binaryaverage to the cut locus, too, by picking one arbitrary shortest pathconnecting 𝑃 to 𝑄 , hence giving up continuity at the cut loci. Wethus define the manifold average between two points A : 𝑀 × 𝑀 × [ , ] −→ 𝑀 ; ( 𝑃, 𝑄 ; 𝑤 ) ↦→ 𝛾 𝑃,𝑄 ( 𝑤 ) (4)where 𝛾 𝑃,𝑄 is a shortest geodesic path joining 𝑃 to 𝑄 . We have that A( 𝑃, 𝑄 ; 𝑤 ) = 𝑅𝐶𝑀 ( 𝑃, 𝑄 ; ( − 𝑤 ) , 𝑤 ) as long as 𝑃 and 𝑄 do not lie ∗ For a definition of cut locus, normal ball and other terms related to the geodesic metricrefer to Appendix A.3 anicinelli et al.
Fig. 4. Example of a failure case of direct De Casteljau evaluation. The blackbullets at the discontinuities correspond to consecutive parameter valuesnear a critical value, and the blue/purple/pink lines provide the De Casteljauconstruction. Note how the pink line jumps from one side of the pole to theother as 𝑡 passes critical values, causing discontinuities. on each other’s cut locus. The averaging operator of Eq. 4 providesthe analogous of Eq. 2 in the manifold setting. The De Casteljau construction provides a recursive definition, whichevaluates a Bézier curve at each 𝑡 ∈ [ , ] as b 𝑘 ( 𝑡 ) = b 𝑘 ( 𝑡 ) , where b 𝑖 ( 𝑡 ) = 𝑃 𝑖 b 𝑟𝑖 ( 𝑡 ) = ( − 𝑡 ) b 𝑟 − 𝑖 ( 𝑡 ) + 𝑡 b 𝑟 − 𝑖 + ( 𝑡 ) (5)for 𝑟 = , . . . , 𝑘 and 𝑖 = , . . . , 𝑘 − 𝑟 . The construction for 𝑘 = and 𝑡 = . is exemplified in the left image of Fig. 6.This construction can be extended to the manifold setting in astraightforward way, by substituting the affine averages betweenpairs of points with the manifold average A defined above. Thisextension was proposed first by [Park and Ravani 1995]. As shownby [Popiel and Noakes 2007], if all consecutive pairs of controlpoints of the control polygon Π lie in a totally normal ball, thenthe resulting curve is 𝐶 ∞ . However, if the constraint is violated, theresulting curve can be discontinuous. In fact, even if all shortestgeodesics paths in Π are unique, some pairs of intermediate pointsinvolved in the construction may lie on each other’s cut locus, forsome value of the parameter 𝑡 . As 𝑡 passes such critical value, themanifold average A returns a discontinuous result, thus causing adiscontinuity in the curve. Figure 4 illustrates the construction nearfailure points, while Figure 5(a) provides another example of failure. A Bézier curve can be evaluated in closed form as an affine sum ofall its control points: b 𝑘 ( 𝑡 ) = 𝑘 ∑︁ 𝑖 = 𝐵 𝑘𝑖 ( 𝑡 ) 𝑃 𝑖 (6)where the 𝐵 𝑘𝑖 ( 𝑡 ) are the Bernstein basis polynomials of degree 𝑘𝐵 𝑘𝑖 ( 𝑡 ) = (cid:18) 𝑘𝑖 (cid:19) 𝑡 𝑖 ( − 𝑡 ) 𝑛 − 𝑖 . This expression can be rewritten for the manifold case as 𝔟 𝑘 ( 𝑡 ) = 𝑅𝐶𝑀 ( 𝑃 , . . . , 𝑃 𝑘 ; 𝐵 𝑘 ( 𝑡 ) , . . . , 𝐵 𝑘𝑘 ( 𝑡 )) . (7)This construction was mentioned, but not developed further, in[Absil et al. 2016; Panozzo et al. 2013].If the control points are close enough to fulfill the Karcher con-dition, then the resulting curve is 𝐶 ∞ , since both the RCM and theBernstein polynomials are 𝐶 ∞ . However, the Karcher condition is (b)(a)(c) (d) Fig. 5. An example of failure in tracing a curve with the direct De Casteljau(a) and the RCM evaluation (b). The same control polygon gives two smoothand nearly identical curves with the Recursive De Casteljau (c) and theOpen-uniform Lane-Riesenfeld schemes (d) described in Sections 3.5 and3.6, respectively. even more restrictive than the constraints required for the De Castel-jau construction, and it is not likely to be verified when the controlpoints lie far apart on a general surface. If the Karcher condition isnot fulfilled, Eq. 3 is no longer convex, and it might even have infin-itely many minima. In this case, the curve may be undetermined atsome intervals. Figure 5(b) provides an example of failure.
A Bézier curve of degree 𝑘 can be rewritten as a curve of degree 𝑘 + on a control polygon Π 𝑘 + with control points 𝑃 𝑘 + 𝑖 = 𝑖𝑘 + 𝑃 𝑘𝑖 − + (cid:18) − 𝑖𝑘 + (cid:19) 𝑃 𝑘𝑖 , 𝑖 = , . . . , 𝑘 + . (8)The repeated application of degree elevation generates a sequenceof control polygons Π 𝑛𝑑𝑒 that converges to the Bézier curve. Thesame stencil can be rewritten in the manifold setting as 𝑃 𝑘 + 𝑖 = A( 𝑃 𝑖 − , 𝑃 𝑖 ; 1 − 𝑖 /( 𝑘 + )) . (9)To the best of our knowledge, this scheme was never applied beforein the manifold setting.Note that, by construction, the control polygon Π 𝑛 + 𝑑𝑒 is formedby either portions of the segments of Π 𝑛𝑑𝑒 , or shortcuts betweentwo consecutive segments. By triangular inequality, this means thatthe polygon becomes progressively shorter and consecutive points 𝑃 𝑛𝑖 , 𝑃 𝑛𝑖 + become progressively closer as 𝑛 increases. In a finite num-ber of steps, each pair of consecutive points will be sufficiently closeto lie in a totally normal ball, hence, from [Popiel and Noakes 2007],the corresponding curve (of degree 𝑛 ) is 𝐶 ∞ and can be evaluatedwith the direct De Casteljau algorithm of Section 3.2.Unfortunately, the repeated application of Eq. 9 has a cost 𝑂 ( 𝑛 ) to generate a polygon with 𝑛 points, because all points are replaced /Surf: Interactive Bézier Splines on Surfaces at each iteration, while their number just increases by one. Also,point evaluation for a relatively large value of 𝑛 may become expen-sive. Hence, this approach is prohibitive for interactive usage. One step of the De Casteljau evaluation subdivides polygon Π intotwo control polygons Π 𝐿 and Π 𝑅 . See Fig. 6 (RDC). The junctionpoint of Π 𝐿 and Π 𝑅 lies on the curve. The recursive applicationof this procedure for 𝑡 = / defines a sequence of subdivisionpolygons Π 𝑛𝐷𝐶 , which converges to the Bézier curve.The extension to the manifold setting is straightforward, bymeans of the point evaluation procedure described in Section 3.2.In the following, we denote this scheme RDC for short. This exten-sion was studied first by [Noakes 1998], and implementations wereproposed in [Morera et al. 2008; Sharp and Crane 2020].Note that here the averaging operator A is not used to tracethe curve, but rather to estimate a discrete set of points, namelythe control points of the two sub-polygons Π 𝐿 and Π 𝑅 . As in theprevious case, the control polygons become progressively shorterand, in a finite number 𝑛 of subdivision steps, all control points ofeach single polygon will lie in a convex neighborhood.Let Π 𝑛𝐷𝐶 be the corresponding chain of control polygons. Fol-lowing [Noakes 1998], we know that each control polygon in Π 𝑛𝐷𝐶 will converge to a curve that is at least 𝐶 . In order to ensure 𝐶 continuity at all junction points between consecutive polygons, wenote that, by construction, the last segment of Π 𝐿 and the first seg-ment of Π 𝑅 form a geodesic line, and also that the first segmentof Π 𝐿 and the last segment of Π 𝑅 are portions of the first and thelast segment of their parent polygon, respectively. By induction onthe recursion tree, this ensures that the two segments incident ateach junction point in Π 𝑛𝐷𝐶 are branches of the same geodesic line,thus satisfying the condition stated in [Popiel and Noakes 2007] towarrant 𝐶 continuity at the junction point.We conclude that the limit curve of this subdivision exists, it is atleast 𝐶 , and, at the endpoints, it interpolates the control polygonand it is tangent to it. It remains an open question whether thissubdivision curve has a higher degree of smoothness. To achieve higher continuity, we propose a novel definition basedon subdivision. A Bézier curve of degree 𝑘 can be represented withan open-uniform B-spline † of degree 𝑘 (order 𝑘 + ), having thesame control polygon Π 𝑘 , and knot vector ( . . . . . . ) , wherethe and are repeated 𝑘 + times. Repeated knot insertion atthe midpoint of all non-zero intervals in the knot vector produces asequence of open uniform B-splines, all describing the same curve,whose control polygons Π 𝑛𝐿𝑅 converge to the curve itself.This subdivision process follows an open-uniform Lane-Riesenfeldscheme (OLR). The control points are nearly doubled at each levelof subdivision, by applying the standard even-odd stencils of theuniform Lane-Riesenfeld scheme “in the middle” [Lane and Riesen-feld 1980], while stencils for end conditions are applied near the † A B-spline is said to be open-uniform if it is uniform, except at its endpoints, whererepeated knots are inserted to make the curve interpolate the endpoints of its controlpolygon.
RDC OLR P
OLR (right) : One step ofsubdivision from Π (blue) to Π (red polygon). The even points 𝑃 𝑗 , aswell as the intermediate points 𝑄 𝑖 and 𝑅 𝑖 lie on segments of Π and areevaluated first. The evaluation of each odd point 𝑃 𝑗 + requires computingone shortest geodesic path (purple). This construction corresponds to onemidpoint subdivision followed by two steps of smoothing by averagingconsecutive points. endpoints. A number of 𝑘 − special stencils are needed at eachend of the polygon. The full constructions for the quadratic andcubic cases, as well as a sketch of the construction for a genericdegree 𝑘 > , are provided in Appendix B. One step of subdivisionfor 𝑘 = and 𝑛 = is exemplified in Figure 6 ( 𝑛 = is the first levelin which the stencils of the uniform LR subdivision apply).It is important to notice that, all affine averages necessary tocompute the stencils in this scheme can be factorized into weightedaverages between pairs of points, as shown in Appendix B. Thisfeature is intrinsic to the uniform L-R scheme, and it is easily gener-alized to the end conditions in the open-uniform scheme. Therefore,this scheme can be extended to the manifold setting, by substitut-ing each affine average with the corresponding application of themanifold average A . We omit the details for the sake of brevity. Tothe best of our knowledge, this scheme was never applied before inthe manifold setting.We show that this construction has strong continuity guarantees.Let Π 𝑛𝐿𝑅 be the subdivision polygon in the manifold setting, whichhas been obtained with 𝑛 levels of subdivision, starting at Π 𝐿𝑅 = Π 𝑘 .Let us consider two consecutive control points on Π 𝑛𝐿𝑅 , namely 𝑃 𝑛𝑗 , 𝑃 𝑛𝑗 + , with < 𝑗 < 𝑛 + , i.e. no such point is an endpoint. For 𝑛 sufficiently large, points 𝑃 𝑛𝑗 , 𝑃 𝑛𝑗 + are contained in a totally normalneighborhood and they are both far enough from the endpoints toundergo the standard LR scheme at the next level of subdivision.Applying the results of [Duchamp et al. 2018], this is sufficientto conclude that the sequence Π 𝑛𝐿𝑅 converges to a 𝐶 𝑘 − curve inthe open interval ( , ) . Concerning the endpoints of Π 𝑛𝐿𝑅 , the endconditions in Equations 10 and 11 guarantee that the limit curve anicinelli et al. Fig. 7. The geodesic line corresponding to the central segment of the controlpolygon can take two different routes at the cut locus of its endpoints, thusproducing two different curves. In practice, the curve will jump between thetwo configurations when dragging a control point across a cut locus. Thislimitation is intrinsic to the discontinuity of the geodesic metric and can beovercome by splitting a long curve into a spline with shorter segments. interpolates the endpoints of the initial control polygon Π 𝑘 , and itwill be tangent to the first and last segments of Π 𝑘 at its endpoints.Hence, the curve obtained with this method is 𝐶 𝑘 − , and providesenough control to weld it with 𝐶 continuity to other curves. The first main contribution of our work is the above analysis, fromwhich we can derive safe schemes for lifting Bézier curves to themanifold domain.The Bernstein point evaluation cannot be used safely on man-ifolds, since the Riemannian center of mass is fragile and can beused only “in the small”. Using direct De Casteljau evaluation isalso problematic. In fact, while the manifold average A for a pair ofpoints is defined over the whole domain, and smooth everywhereexcept at the cut loci, the discontinuity of operator A at the cut locimakes the direct De Casteljau evaluation fragile.Conversely, subdivision schemes can be safely defined with re-peated averages based on operator A and can work on any controlpolygon, since the arbitrary choices made with operator A at thecut loci do not affect convergence and smoothness.Of the subdivision schemes, the one based on degree elevationguarantees a 𝐶 ∞ curve, but it is computationally expensive andscales purely with curve complexity. At best, it can be used off-line,for instance, in order to obtain the fine evaluation of a curve thathas been designed interactively with one of the other schemes.This leaves two practical schemes for interactive Bézier curveson manifolds. The RDC scheme, based on the recursive De Casteljaubisection, guarantees 𝐶 continuity and it can be implemented easilyand efficiently via repeated geodesic averages. The OLR scheme,based on open-uniform Lane-Riesenfeld subdivision, guarantees 𝐶 𝑘 − smoothness for a curve of degree 𝑘 , and can be made to workefficiently, too. The arbitrary choices made with operator A atthe cut loci can in principle lead to different curves for the samecontrol polygon. In practice, since all our algorithms are determinis-tic, the same paths will always be chosen at the cut loci, for a given control polygon, thus returning the same curve. However, the curvemay jump to a different configuration for small displacements ofcontrol points, which make some of the paths in the constructioncross a cut locus. In Figure 7, a tiny displacement of one controlpoint takes one of the shortest paths in the control polygon to adrastically different route, resulting in a different curve. Note thatsuch jumps occur quite rarely, as the cut locus of each point coversa set of zero measure on the manifold. This fact is intrinsic to thediscontinuity of the manifold metrics and constitutes essentially theonly limitation of our approach. We now focus deriving practical algorithms for the RDC and OLRschemes, based on recursive De Casteljau bisection and on non-uniform LR subdivision, respectively. We provide algorithms forapproximating the curve with a geodesic polyline (curve tracing),evaluating a point on the curve for a given parameter value (pointevaluation), and splitting a curve at a given point into a spline withtwo segments (point insertion). We support surfaces represented astriangle meshes, and target interactivity for long curves and meshesof millions of triangles.In order to develop our algorithms, we assume to have proceduresfor (1) computing the point-to-point shortest path between pairsof points of 𝑀 ; (2) evaluating a point on a geodesic path at a givenparameter value; and (3) casting a geodesic path from a point in agiven direction. In all these cases, we consider generic points on thesurface, and not just the vertices of the mesh. The computational de-tails of such procedures, as well as additional algorithms to supportinteractive control, are provided in Section 5. We trace Bézier curves on surfaces by approx-imating them with a geodesic polygon. The tracing algorithm is arecursive subdivision that, at each step, takes a geodesic polygon Π and produces two sub-polygons Π 𝐿 and Π 𝑅 . Recursion is initial-ized by computing the 𝑘 shortest paths that constitute the polygonconnecting the initial control points 𝑃 , . . . , 𝑃 𝑘 .To split a control polygon, we compute a sequence of geodesicpolygons Π 𝑖 for 𝑖 = , . . . , 𝑘 , each containing 𝑘 − 𝑖 segments, where Π = Π , and Π 𝑘 degenerates to the midpoint of the curve. Since Π is known from recursion, this requires computing a total of 𝑘 ( 𝑘 − )/ further geodesic paths, each joining the midpoints ofthe segments in Π 𝑖 , in order to produce the points of Π 𝑖 + . Thepolygon Π 𝐿 is built by collecting all sub-paths that connect thestarting points of the polygons in the sequence, namely Π 𝑖 [ ] to Π 𝑖 + [ ] for 𝑖 = , . . . , 𝑘 − . The polygon Π 𝑅 is built likewise, byusing sub-paths connecting the end points of subsequent polygons.We support both uniform and adaptive subdivision. For uniformsubdivision, a maximum level of recursion is chosen either by theuser, or automatically computed on the basis of the total length 𝐿 ( Π ) of the initial polygon Π , and a threshold 𝛿 > . Since the pathsforming the subdivided polygon are shrinking through recursion,then after ⌈ log ( 𝐿 ( Π )/ 𝛿 )⌉ recursion levels, the length of a geodesicpath in the output will be no longer than 𝛿 . For adaptive subdivision, /Surf: Interactive Bézier Splines on Surfaces P
We encode triangle meshes 𝑀 with a simple indexed data structureconsisting of three arrays encoding the vertices, the triangles, andtriangles adjacencies, which provide the dual graph having thetriangles as nodes.We need to deal with generic points lying on the mesh, not justits vertices. A mesh point 𝑃 is encoded as a triple ( 𝑡, 𝛼, 𝛽 ) where 𝑡 isthe triangle index and 𝛼, 𝛽 are the barycentric coordinates of 𝑃 in 𝑡 . A vertex 𝑣 of 𝑀 can also be encoded as a generic mesh point, bymeans of any of its incident triangles.A geodesic path connecting two mesh points 𝑃 and 𝑄 is encodedwith a triangle strip ( 𝑡 , . . . , 𝑡 ℎ ) , where 𝑡 and 𝑡 ℎ contain 𝑃 and 𝑄 , respectively, and an array of real values 𝑙 , . . . , 𝑙 ℎ − , where 𝑙 𝑖 encodes the intercept of the path with the edge 𝑒 𝑖 common to 𝑡 𝑖 , 𝑡 𝑖 + ,parametrized along 𝑒 𝑖 . Point-to-point shortest path.
The literature offers several tech-niques for computing shortest paths [Crane et al. 2020]. We proposean algorithm that is derived by combining insights from the worksof [Lee and Preparata 1984; Xin and Wang 2007]. The algorithmconsists of three phases: (i) extraction of an initial strip; (ii) shortestpath in a strip; and (iii) strip straightening.Phase (i), which has been overlooked in several previous works, iscritical as it can become the bottleneck on large meshes. Given twomesh points 𝑃 and 𝑄 , we compute a strip of triangles that connectsthem, performing a search on the dual graph. We experienced a rel-evant speedup over the classical Dijkstra search by using a shortestpath algorithm based on the SLF and LLL heuristics [Bertsekas 1998],which do not require a priority queue, but just a double ended queue.The SLF and LLL heuristics govern the insertion and extraction ofweighted nodes in the queue. We weight each node as in a classicalA* search, with the sum of its current distance from the source plusits Euclidean 3D distance to the target. This heuristic prioritizesthe exploration of triangles closer to the destination in terms ofEuclidean distance, improving performance in most models.In phase (ii), the strip is unfolded in the 2D plane and the shortestpath within it is computed in linear time with the funnel algorithm[Lee and Preparata 1984]. See Fig. 9(a-b-c) for an example.In phase (iii), in order to obtain the locally shortest path on themesh, we remove reflex vertices from the strip where possible. Tothis aim, [Xin and Wang 2007] finds the reflex vertices that can beremoved by computing angles about a vertex inside and outside thestrip, respectively. However, in our experiments, the computationof angles slows down the algorithm, because the star of each reflexvertex is retrieved from a data structure that is not in cache memory.Instead, we select the reflex vertex 𝑣 that creates the largest turnin the polyline and, similarly to [Xin and Wang 2007], we updatethe strip by substituting the current semi-star of 𝑣 inside the stripwith its other semi-star. We perform the unfolding and the funnelalgorithm on the new strip: if 𝑣 still remains on the path, then it isfrozen; we repeat this procedure until all reflex vertices either areremoved or become frozen. See Fig. 9(d-e-f) for an example. Straightest geodesics.
A straightest geodesic is traced starting ata mesh point 𝑃 and following one given direction u in the tangent /Surf: Interactive Bézier Splines on Surfaces QP (a) (b) (c)(f)(e)(d) P PP PPQv Q QQQ
Fig. 9. Shortest path computation. Given a source 𝑃 and a target 𝑄 an initialstrip of triangles connecting them is found with a search on the dual graphof the mesh. (a) A shortest path within the strip is found by propagating afunnel, which is initialized with its apex at 𝑃 and its front at the first edgecrossing the strip. (b) The edges of the strip are processed one by one, totighten the front of the funnel. (c) When the funnel collapses, a new vertex,called a pseudo-source, is added to the path and the apex of the funnel ismoved to the pseudo-source. (d) When 𝑄 is reached, some reflex verticesmay still lie on the path. (e) Reflex vertices are analyzed for possible removal,starting at the vertex 𝑣 causing the sharpest turn. (f) The final path is foundwhen no more reflex vertices can be removed. plane 𝑇 𝑃 𝑀 for a given length. This is done by unfolding the trianglesthat are crossed by the line as it is being traced and intersectingtheir edges with the line in the 2D domain. The ending results is astraight 2D line that crosses a strip of unfolded 2D triangles, whichcan be mapped back into a geodesic surface path thanks to ourrepresentation. In case the path intersects a vertex 𝑣 , then we follow[Polthier and Schmies 1998], reflecting the incoming direction about 𝑣 in its tangent plane. Parallel transport.
We use an approach similar to [Knöppel et al.2013]. Each triangle has its own tangent 2D frame of reference. Adirection at a mesh point is represented by 2D coordinates withrespect to the frame of its containing triangle. When a directionmust be parallely transported from a mesh point 𝑃 to another 𝑄 , arotation must be applied to take into account the fact that 𝑃 and 𝑄 may lie on different triangles, with different frames of reference.To do so, the geodesic path from 𝑃 to 𝑄 is computed, then the stripcontaining the path is unfolded so that the coordinates of the axis ofthe frame of reference of 𝑄 ’s triangle can be expressed with respectto the frame of reference of 𝑃 ’s triangle, hence the rotation betweenthe two frames can be obtained. Leveraging the proposed algorithms, we developed a graphical ap-plication to allow users to interactively edit splines on meshes,imitating the same interaction of established 2D vector graphicstools. We focus on cubic curves since they are the most used in 2D.Our application supports the editing of curves by moving, adding,and deleting control points, and by translating, scaling and rotat-ing whole splines on the surface domain. All these operations aresupported by using the geodesic primitives just described. Here wedescribe the main editing feature.
Curve editing.
Borrowing the editing semantic from 2D tools,control points are distinguished in anchor points and handle points.Anchors are those points where two Bézier curves are joined, hencea spline passes through all of its anchor point. The preceding and fol-lowing control points of an anchor are its associated handle points.The handle points of an anchor determine two segments, both start-ing at the anchor itself. A spline is tangent to those segments at theanchor points.In the 2D setting, when an anchor is dragged, the two tangentsegments move with it and so do the associated handle points. Toobtain the same behavior on the surface, when moving an anchorpoint from 𝑃 to 𝑃 ′ , we find the two tangent directions of the tangentsegments at 𝑃 . Then, for each such segment, we trace a straightestgeodesics starting at 𝑃 ′ end for the same length of the segment, inthe direction of its tangent, rotated by the parallel transport from 𝑃 to 𝑃 ′ . The endpoint of each segment is the new position of thecorresponding handle point.In the 2D setting the user can impose an anchor to be "smooth",i.e. the two associated tangent segments are always colinear, whichautomatically ensure 𝐶 continuity at the anchor point. To providethe same functionality on the surface, whenever the handle point 𝑄 is moved, the opposite handle point 𝑄 is recomputed by tracinga straightest geodesic from the anchor 𝑃 along the tangent directiondefined from segment 𝑄 𝑃 to find the new position of handle 𝑄 . Rotation, Scaling and Translation.
Our application also supportstranslation, rotation and scaling of a whole spline. In the 2D set-tings these operations are obtained by just applying the same affinetransform to all control points of a spline. In the surface setting, wedefine the center of the transformation 𝐶 to be just the mesh pointunder the mouse pointer.To apply the transformation, the normal coordinates of the controlpoints are computed with respect to the center 𝐶 , in a sort of discreteexponential map. Then, the linear transformation is applied on these2D coordinates, which are finally converted back into mesh pointsby tracing straightest geodesic paths outward from 𝐶 , just like wedo for mapping 2D drawings.Translation needs special handling, as the center of the transfor-mation 𝐶 is dragged to a new position 𝐶 ′ . To compensate for thechange of reference frame, the normal coordinates are rotated bythe opposite angle of the parallel transport given by the tangentvector from 𝐶 to 𝐶 ′ .Note that, while the exponential map is not reliable to provide adense map, we apply normal coordinates just to a relatively smallset of sparse points. In this case, we can tolerate the distortionscaused by the curvature of the surface. Sometimes, it may be convenient to map a whole 2D vector drawing,made of several primitives in the Euclidean plane, to the surface.Note that, unlike standard methods based on parametrization, weare not mapping the result of the drawing, but rather its controlpoints: the final drawing is traced directly on the manifold, basedon its vector specification, and can be further edited after mapping.See Figures 1 and 10 for examples. anicinelli et al. Fig. 10. Example of importing a large SVG, made of 2269 splines, onto theArmadillo, consisting of 346K triangles. Our algorithm takes 60 millisecondsto import the control points and milliseconds to trace all curves.
This mapping is just meant to provide an initial placement of thecontrol points on the target surface, allowing the user to adjust andfine tune the drawing afterwards. Therefore, we can allow for somedistortion in the initial placement.Our method is analogous to [Biermann et al. 2002], and it is basedon the conversion between polar coordinated in 2D and normalcoordinates on the manifold. Each point of the SVG drawing isconverted into a mesh point by taking its polar coordinates, andtracing a geodesic from a center point in the given tangent direction,for the given distance.
We validate our work by tracing curves over a large number ofmeshes, by comparing it with a recent state-of-the-art solution,and by performing interactive editing sessions, as shown in theaccompanying video. In summary, our algorithms produce a validoutput in all trials, in a time compatible with interactive usage inover 99% of the trials (Table 1); and we beat the state-of-the-art byone order of magnitude in speed (Fig. 12).Concerning interactive usage, our system supports editing inall conditions for meshes of the order of one million triangles ona laptop computer. Interaction has been tested with meshes withseveral millions of triangles (see, e.g., Fig. 1) and interaction is stillsupported, provided that single curves do not span too large a frac-tion of the model. Such cases are rare in actual editing sessions, asreal designs are usually made of many splines, each consisting ofseveral small segments.
We tested our algorithms for robustness by running a large ex-periment on the Thingi10k repository [Zhou and Jacobson 2016].Our algorithm requires that the mesh is manifold and watertight,so we extracted the subset of meshes that have those properties,for a total of 5567 models. The models are used as is, without anypre-processing.For each model, we consider 100 random cubic curves. For eachcurve, we take the model in its standard pose, and pick points on algorithm percent of trials times at percentile< 0.001s < 0.1s 90% 99%RDC Uniform 71.3% 99.4% <0.0060 <0.0668OLR Uniform 73.1% 99.5% <0.0052 <0.0570RDC Adaptive 71.7% 99.4% <0.0058 <0.0661OLR Adaptive 64.8% 99.2% <0.0080 <0.0872
Table 1. Time performances of our algorithms in 556,700 trials. We reportthe percentage of trials in which tracing a curve takes less than 0.001 and 0.1seconds, and the running times at the 90th and 99th percentiles, respectively. it by casting random rays orthogonal to the view plane, until wefind four points that lie on the surface. These become the controlpoints of the spline. We place no restriction in the arrangement ofthe control points. This gives us a total of more than half millioncontrol polygons.For each test, we run both the RDC and the OLR tracing algo-rithms, in their uniform and adaptive configurations. The uniformRDC algorithm is expanded to 4 levels of recursion, which gener-ates a geodesic polyline consisting of 48 geodesic segments. Theuniform OLR algorithm is expanded to 6 levels of recursion, whichgenerates a geodesic polyline consisting of 64 geodesic segments. Infact, because of the different subdivision rules, we cannot generatethe same number of segments for both schemes. For the adaptivevariants, we set a threshold 𝜃 = ◦ for the maximum angle betweenconsecutive geodesic segments along the polyline. In this case, thenumber of geodesic segments in output is variable, depending onthe curve and on the method. Since all algorithms generate verysimilar curves, the final tessellated paths that approximate the curveon the mesh, consisting of one line segment per triangle crossed,have about the same number of segments in all four cases.Trials were executed on a Linux PC with an AMD Ryzen 5 2600xand 32GB memory, running on a single core in all experiments. Allour algorithms passed all the tests and generated curves that appearto be smooth in all cases.In Table 1 and Figure 11, we compare the timing performanceof the four algorithms. All algorithms perform quite similarly, andremain interactive in all cases, with roughly 70% of trials runningat less than 1 millisecond per curve, and 99% of the trials runningfaster than 0.1 second/curve. The few trials in which they take moretime are concerned, with very few exceptions, either with very longcurves on large meshes (>1M triangles), or with meshes containinga large fraction of extremely spiky triangles (which are numerousin the Thingi10K repository).There are small differences in the performances of the differentalgorithms. For uniform subdivision, the OLR algorithm resultsslightly faster than the RDC algorithm, beside generating a morerefined geodesic polyline. On the contrary, for adaptive subdivision,the RDC algorithm runs slightly faster than the OLR algorithm.These differences are probably due to the simpler structure of theOLR uniform algorithm in one case, and to the more involved struc-ture of the OLR adaptive algorithm in the other. In fact, both variantsof the RDC algorithm follow the same recursive pattern. On thecontrary, the uniform OLR algorithm expands the curve level by /Surf: Interactive Bézier Splines on Surfaces <0.1s<0.001s<0.1s<0.001s<0.1s<0.001s<0.1s<0.001s Fig. 11. The distributions of running times of our four algorithms for curve tracing in 556,700 trials on 5,567 models from the Thingi10k repository, tracing100 random curves on each model. All algorithms provide a valid output in all trials. The different algorithms have similar behavior and are compliant withinteraction (<0.1 seconds/curve) in more than 99% of the trials. For uniform subdivision, the OLR algorithm is slightly faster than the RDC algorithm; while foradaptive subdivision, the RDC algorithm performs slightly better than the OLR algorithm. Adaptive algorithms have slightly narrower distributions thanuniform algorithms. level, following a simpler pattern; while the OLR adaptive algorithmrequires a recursive pattern, with a slightly more involuted struc-ture than the RDC algorithms. Note that the adaptive algorithmsshow slightly narrower distributions, corresponding to a more stableperformance since the stopping criterion is related to curve shape.For the sake of brevity, we do not present here results on thealgorithms for curve tracing e point insertion, which run muchfaster than the tracing algorithms.In the previous experiments, the cost of computing a curve de-pends on both the length of the curve and the size of the mesh, withtrends that are not linear. Roughly speaking, the cost of finding theinitial path of geodesics depends on both the length of the curveand the size of the mesh, while the subsequent cost of finding theshortest path depends just on the length of the curve. As the relativelength of the curve grows, the cost of finding the initial path prevails,since it may requires exploring most of the mesh.
The flipOut algorithm was proposed recently [Sharp and Crane2020] as a fast solution to the computation of locally shortest geo-desic paths. On the basis of the flipOut algorithm, the same authorshave implemented the algorithm of [Morera et al. 2008], which usesthe same outer scheme of our RDC algorithm for curve tracing.We have used the implementation provided by the authors [Sharpet al. 2019] to run the same experiments of Sec. 6.1, with the sameparameter used for our RDC algorithm with uniform expansion.Unfortunately, the flipOut algorithm has some limitations thatprevent it from running on any control polygon. In particular, itrequires that the control polygon does not contain self-intersections,a case which is pretty common with cubic curves, and happens onaverage with 1/3 of randomly generated polygons. Therefore, weexcluded all the trials for which flipOut could not provide an output,keeping a total of 80,376 out of 556,700 trials.From a visual inspection of random samples of the results, itseems that both algorithms generate the same curves. In Figure12 we present a comparison between the performances of the twoalgorithms. Our RDC uniform algorithm exhibits a speedup of about
10x 2x1x14x
Fig. 12. The graph shows the distribution of the ratio of the running timesbetween our RDC uniform algorithm and the [Sharp and Crane 2020] basedalgorithm. Here we report only the 80,376 trials, out of 556,700, for which[Sharp and Crane 2020] could provide an output. The red vertical line showsthe average of the distribution at 0.073, corresponding to about 14x speedup;the remaining lines provide thresholds for 10x, 2x, and 1x speedups fromleft to right respectively.
14x on average; it runs faster than flipOut on 99.9% of the trials,being at least 2x on 99% and at least 10x on 77.7% of the trials. The flipOut method was faster than ours only on 64 over 80,376 trials, allof which, except one, on meshes with a size lower than 30K triangles.Our speedups seem for the most part due to phase (i) of the shortestpath algorithm described in Sec. 5.2.Besides, all our algorithms have no limitations, since they couldprovide a valid output in all 556,700 trials.
We have used extensively our system on a variety of models. Allediting sessions where performed on a MacBook laptop with a anicinelli et al. Fig. 13. A gallery of models from the Thingi10k repository, each with curves drawn with our method. The selected models span a wide range of shapes andthe sizes of meshes vary between about 120K and 5.7M triangles.
We propose methods for interactively drawing and editing of Béziercurves on meshes of millions of triangles, without any limitation onthe curve shape and extension of control polygons. Our algorithmsare robust, having been been tested on over five thousands shapeswith randomly generated control polygons, and they beat a state-of-the-art solution by one order of magnitude in speed. Our new Open-uniform Lane-Riesenfeld scheme provides the smoothest practicalsolution so far for Bézier curves in the manifold setting, while our De Casteljau variation is simple to implement, at the price of lesssmoothness.The main limitation of these methods lie in the discontinuitiesof the space of curves with respect to their control points: curvesare always smooth, but they may make jumps during editing. Sucha discontinuity is inherent of the geodesic metric, and it can beovercome by using a spline with shorter control polygons, insteadof a single large polygon, to define the curve. Our algorithms forpoint insertion greatly help in this task.In the future, we want to consider other types of splines. Anextension of our approach to B-splines seems straightforward. Anextension to interpolating splines seems easy, but it requires mani-fold extrapolation, which may become unstable. The most complexextension would be to handle NURBS, which at this point remainsunclear how to do. More generally, the smoothness analysis in thenon-uniform case needs a thorough investigation.
ACKNOWLEDGMENTS
We wish to thank Chiara Eva Catalano, Tom Duchamp, Kai Hor-mann, Daniele Panozzo, and Giuseppe Patané for helpful discussions;Marzia Riso for her help with experiments and figures; and MicheleSerpe for modeling a mesh for us. Models other than Thingi10k:Nefertiti is courtesy of Scan the World; Armadillo and Bunny arecourtesy of the Stanford 3D Scanning Repository. /Surf: Interactive Bézier Splines on Surfaces v i
In the following, we provide just a summary of basic concepts. Fora complete account on this subject, we refer to [do Carmo 1976](smooth setting) and [Crane et al. 2013] (discrete setting).
Smooth setting.
Let 𝑀 be a smooth surface embedded in R . Theembedding induces a Riemannian metric on 𝑀 , defining the length 𝐿 ( 𝛾 ) of any parametric curve 𝛾 on 𝑀 . Without lack of generality,we consider 𝛾 to be parametrized in [ , ] with constant speed.The geodesic distance 𝑑 ( 𝑝, 𝑞 ) between any two points 𝑝, 𝑞 ∈ 𝑀 isthe infimum of length 𝐿 over all curves 𝛾 such that 𝛾 ( ) = 𝑝 and 𝛾 ( ) = 𝑞 ; and one such curve ¯ 𝛾 satisfying 𝐿 ( ¯ 𝛾 ) = 𝑑 ( 𝑝, 𝑞 ) is called a shortest geodesic path between 𝑝 and 𝑞 .A shortest geodesic path may not be unique. For any 𝑝 ∈ 𝑀 , theset of points 𝑞 such that there exist more than one geodesic shortestpath connecting 𝑝 to 𝑞 belong to the cut locus of 𝑝 . If 𝑞 stays awayfrom the cut locus of 𝑝 , then both the geodesic distance 𝑑 ( 𝑝, 𝑞 ) andthe (unique) shortest path 𝛾 𝑝,𝑞 joining them vary smoothly with 𝑞 . On the contrary, if 𝑞 moves along a trajectory that crosses thecut locus of 𝑝 , then the distance 𝑑 ( 𝑝, 𝑞 ) changes with continuity,but not smoothly, at the cut locus; while 𝛾 𝑝,𝑞 may jump to a totallydifferent curve. These facts are relevant to assess the stability of themethods discussed in Section 3.A path 𝛾 𝑝,𝑞 is said to be locally shortest if there exist 𝛿 > suchthat for any interval [ 𝑡 − 𝛿, 𝑡 + 𝛿 ] ⊆ [ , ] the restriction of 𝛾 𝑝,𝑞 to [ 𝑡 − 𝛿, 𝑡 + 𝛿 ] is a shortest path. A curve connecting 𝑝 to 𝑞 canbe homotopically deformed to a locally shortest path; this is at thebasis of all algorithms for computing locally shortest paths.Geodesic curves can be also characterized by their straightness .In order to assess the curvature of lines in the intrinsic geometryof 𝑀 , one needs to introduce the covariant derivative, which weskip here for brevity. Intuitively, from an extrinsic point of view, ageodesic curve 𝛾 does not make any further turn except the strictlynecessary to follow the curvature of 𝑀 : it turns with 𝑀 , but it doesnot turn on 𝑀 . Thus, geodesics play the role of straight lines on 𝑀 .A geodesic curve is completely defined by its starting point andtangent vector. Following this observation, the exponential map exp 𝑝 : 𝑇 𝑝 𝑀 −→ 𝑀 maps vectors of the tangent plane 𝑇 𝑝 𝑀 at 𝑝 to points on the surface. In general, the exponential map is notinjective. A neighborhood 𝑈 of 𝑝 over which exp 𝑝 is invertible iscalled a normal neighborhood of 𝑝 ; the inverse of the exponential map on 𝑈 defines the logarithmic map log 𝑝 : 𝑈 → 𝑇 𝑝 𝑀 , whichprovides local coordinates around 𝑝 called normal coordinates .A set 𝑈 ⊂ 𝑀 is said to be a totally normal neighborhood if itis a normal neighborhood for all its points; and it is said to be strongly convex if it contains all shortest paths between pairs ofits points. The maximum radius that a region 𝑈 can be extendedto, while remaining a [totally] normal neighborhood, or a convexset, depends on the Gaussian curvature of 𝑀 and is not easy toassess. For this reason, methods based on exp and log maps cannotguarantee robustness in the general case. Discrete setting. If 𝑀 is a polyhedral manifold – which, w.l.o.g.,we can consider to be a triangle mesh – then it is no longer smooth.The concepts of shortest geodesic path, geodesic distance and lo-cally shortest path are nonetheless well defined. A geodesic path 𝛾 between two points 𝑝, 𝑞 on 𝑀 is a polyline, intersecting a stripof triangles of 𝑀 . As long as 𝛾 does not cross any vertex of 𝑀 , thesegments of this polyline can be obtained by flattening the strip oftriangles to the Euclidean plane.Geodesics that cross vertices are more complex to handle. A vertex 𝑣 𝑖 ∈ 𝑀 can be classified depending on the sign of its angle defect ,defined as the difference between 𝜋 and the total angle about vertex 𝑣 𝑖 (see Figure 14).It turns out that no (locally) shortest path can pass through avertex with positive angle defect; while infinitely many shortestpaths can reach a vertex with negative angle defect from a givendirection and take divergent directions to the other side of it (Figure14, left). This fact makes the definition of straightest geodesic morecomplicated than in the smooth setting. Following the definitionof [Polthier and Schmies 1998], a straightest geodesic intersectinga vertex 𝑣 is required to bisect the total angle about 𝑣 . With thisdefinition, a straightest geodesic 𝛾 is locally shortest if and only if itdoes not cross any vertex with positive angle defect. B OPEN-UNIFORM LANE-RIESENFELD SUBDIVISION
A Bézier curve with control polygon Π 𝑘 can be rewritten as an open-uniform B-spline of degree 𝑘 (order 𝑘 + ) with the same controlpolygon Π 𝑘 and knot vector ( . . . . . . ) , where the and are repeated 𝑘 + times [Salomon 2006].We apply the generalized Oslo algorithm [Goldman and Lyche1993] to repeatedly insert knots at the midpoints of all non-zerointervals in the knot vector, producing a sequence of open uniform B-splines, all describing the same curve, whose control polygons Π 𝑛𝐿𝑅 converge to curve itself. In the following, we sketch the constructionfor a generic value of 𝑘 , and we provide the complete solution forvalues 𝑘 = , , which are most relevant in the applications.Let Π 𝐿𝑅 = Π 𝑘 be our initial control polygon, and let Π 𝑛 + 𝐿𝑅 be thepolygon obtained from Π 𝑛𝐿𝑅 with one round of knot insertions. Theknot vectors associated to the first levels of the subdivision are: ( . . . . . . )( . . . . . . )( . . . . . . )( . . . . . . )· · · · · · anicinelli et al. where the first and the last node are always repeated 𝑘 + times, andnodes are renumbered at each level. Note that, at level 𝑛 , polygon Π 𝑛𝐿𝑅 contains 𝑛 + 𝑘 control points, and its associated knot vectorcontains 𝑛 + + 𝑘 nodes, with 𝑛 non-null intervals. Each pointof Π 𝑛 + 𝐿𝑅 is computed by applying the Oslo algorithm, as an affineaverage of at most 𝑘 + consecutive points of Π 𝑛𝐿𝑅 , with weightsdepending on 𝑘 + consecutive knots of the knot vector of level 𝑛 .The first few levels of the subdivision need special stencils, obtaineddirectly from knot insertion; as soon as 𝑛 ≥ ⌈ log ( 𝑘 + )⌉ , the knotvector contains at least 𝑘 + non-null intervals, and the subdivisionrules stabilize: we have the standard stencils of the uniform Lane-Riesenfeld subdivision [Lane and Riesenfeld 1980] for control pointsthat depend just on uniform nodes, plus a set of 𝑘 − end conditionsat each side of the curve.In the following, we give the stencils for the cases 𝑘 = , . Quadratic curves ( 𝑘 = ). We have Π = ( 𝑃 , 𝑃 , 𝑃 ) with initialnode vector ( ) . The polygons Π 𝐿𝑅 = ( 𝑃 , 𝑃 , 𝑃 , 𝑃 ) and Π 𝐿𝑅 = ( 𝑃 , ..., 𝑃 ) at the first two levels of subdivision are obtainedwith the special rules 𝑃 = 𝑃 , 𝑃 = 𝑃 + 𝑃 , 𝑃 = 𝑃 + 𝑃 , 𝑃 = 𝑃 and 𝑃 = 𝑃 , 𝑃 = 𝑃 + 𝑃 , 𝑃 = 𝑃 + 𝑃 ,𝑃 = 𝑃 + 𝑃 , 𝑃 = 𝑃 + 𝑃 , 𝑃 = 𝑃 From the second level on, we apply the following general rules: 𝑃 𝑛 + 𝑗 = 𝑃 𝑛𝑗 + 𝑃 𝑛𝑗 + 𝑃 𝑛 + 𝑗 + = 𝑃 𝑛𝑗 + 𝑃 𝑛𝑗 + 𝑗 = ... 𝑛 − 𝑃 𝑛 + = 𝑃 𝑛 𝑃 𝑛 + = 𝑃 𝑛 + 𝑃 𝑛 𝑃 𝑛 + 𝑛 + = 𝑃 𝑛 𝑛 + 𝑃 𝑛 𝑛 + 𝑃 𝑛 + 𝑛 + + = 𝑃 𝑛 𝑛 + (10)where the first two rows are the standard stencils of the Lane-Riesenfeld subdivision of order 1 (which coincides with the Chaikinsubdivision in the uniform case). Note that all stencils are affineaverages of at most two points. Cubic curves ( 𝑘 = ). We have Π = ( 𝑃 , 𝑃 , 𝑃 , 𝑃 ) and initialnode vector ( ) . The polygons at the first and second levelsare obtained with the special rules 𝑃 = 𝑃 , 𝑃 = 𝑃 + 𝑃 , 𝑃 = 𝑃 + 𝑃 , 𝑃 = 𝑃 + 𝑃 , 𝑃 = 𝑃 and 𝑃 = 𝑃 , 𝑃 = 𝑃 + 𝑃 , 𝑃 = 𝑃 + 𝑃 ,𝑃 = 𝑃 + 𝑃 + 𝑃 ,𝑃 = 𝑃 + 𝑃 , 𝑃 = 𝑃 + 𝑃 , 𝑃 = 𝑃 Note that the stencil to compute 𝑃 involves the weighted averageof three points, which can be conveniently factorized as: 𝑄 = 𝑃 + 𝑃 , 𝑄 = 𝑃 + 𝑃 , 𝑃 = 𝑄 + 𝑄 . From the second level on, we apply the following general rules: 𝑃 𝑛 + 𝑗 = 𝑃 𝑛𝑗 + 𝑃 𝑛𝑗 + 𝑗 = ... 𝑛 − 𝑃 𝑛 + 𝑗 + = 𝑃 𝑛 + 𝑗 + 𝑃 𝑛 + 𝑗 + + 𝑃 𝑛 + 𝑗 + 𝑗 = ... 𝑛 − 𝑃 𝑛 + = 𝑃 𝑃 𝑛 + = 𝑃 𝑛 + 𝑃 𝑛 𝑃 𝑛 + = 𝑃 𝑛 + 𝑃 𝑛 𝑃 𝑛 + = 𝑃 𝑛 + 𝑃 𝑛 + 𝑃 𝑛 (11)where, for the sake of brevity, we have omitted the end conditionsto the right end side, since they are symmetric to the end conditionsto the left end side. As in the previous case, the first two stencils inEq. 11 give the uniform LR scheme of order 2, which is applied toall central points. Again, the stencils involving the affine average ofthree points can be conveniently factorized as follows: 𝑄 𝑛 + 𝑗 = 𝑃 𝑛 + 𝑗 + 𝑃 𝑛 + 𝑗 + , 𝑄 𝑛 + 𝑗 + = 𝑃 𝑛 + 𝑗 + + 𝑃 𝑛 + 𝑗 + 𝑃 𝑛 + 𝑗 + = 𝑄 𝑛 + 𝑗 + 𝑄 𝑛 + 𝑗 + and 𝑅 𝑛 + = 𝑃 𝑛 + 𝑃 𝑛 , 𝑅 𝑛 + = 𝑃 𝑛 + 𝑃 𝑛 𝑃 𝑛 + = 𝑅 𝑛 + + 𝑅 𝑛 + . In conclusion, all computations for 𝑘 = , can be factorizedto involve just weighted averages between pairs of points. Thesame approach applies in the case 𝑘 > too, where the 𝑘 − special stencils at the boundaries can be computed through thegeneralized Oslo algorithm [Goldman and Lyche 1993], and theremaining (internal) points are obtained by applying the classical LRalgorithm. Factorization of operations as repeated averages of pairsof points applies for all values of 𝑘 by exploiting the nature of the LRscheme, which consists of one step of midpoint subdivision, followedby 𝑘 − steps of smoothing by averaging (see, e.g., [Goldman 2002],Chapter 7). C FROM B-SPLINE TO BÉZIER
Uniform case.
Given the control polygon ( 𝑃 , ..., 𝑃 𝑘 ) of a uniformB-spline segment b ( 𝑡 ) of degree 𝑘 , the control polygon ( 𝑄 , ..., 𝑄 𝑘 ) of a Bézier curve coincident with b ( 𝑡 ) is given by expression ( 𝑄 , 𝑄 , . . . , 𝑄 𝑘 ) 𝑇 = 𝑀 − 𝑏 𝑀 𝑠 ( 𝑃 , 𝑃 , . . . , 𝑃 𝑘 ) 𝑇 , where 𝑀 𝑠 and 𝑀 𝑏 are the matrices defining the 𝑘 -degree B-splineand the 𝑘 -degree Bézier curve, respectively. Both matrices are welldefined for an arbitrary degree 𝑘 , and their construction can befound in [Yamaguchi 1988].As an example, for 𝑘 = , the above equation leads to 𝑄 = ( 𝑃 + 𝑃 + 𝑃 ) 𝑄 = ( 𝑃 + 𝑃 ) 𝑄 = ( 𝑃 + 𝑃 ) 𝑄 = ( 𝑃 + 𝑃 + 𝑃 ) , where the first average (and similarly the last) can be factorized as ˜ 𝑄 = ( 𝑃 + 𝑃 ) ˜ 𝑄 = ( 𝑃 + 𝑃 ) 𝑄 = ( ˜ 𝑄 + ˜ 𝑄 ) . It is important to point out that the inverse computation of ( 𝑃 , ..., 𝑃 𝑘 ) from ( 𝑄 , ..., 𝑄 𝑘 ) leads to non-convex averages, i.e., withnegative weights. This explains why, in our case, we consider anopen-uniform LR subdivision, rather than computing the control /Surf: Interactive Bézier Splines on Surfaces polygon of a uniform B-spline, since this would have implied theextrapolation of long geodesic lines, which may become unstable. Non-uniform case.
In this case, the entries of 𝑀 𝑠 depend on thespacing between the entries of the knot vector of the B-spline. Theconstruction of 𝑀 𝑠 in the general case can be found in [Qin 2000];we report here, as an example, the case where the knot vector hasthe form ( . . . ) and b ( 𝑡 ) is the cubic B-spline segmentdefined in [ , ] . If ( 𝑃 , 𝑃 , 𝑃 , 𝑃 ) are the control points defining b ( 𝑡 ) , then we have 𝑄 = 𝑃 𝑄 = 𝑃 𝑄 = ( 𝑃 + 𝑃 ) 𝑄 = ( 𝑃 + 𝑃 + 𝑃 ) , where again we can factorize ˜ 𝑄 = ( 𝑃 + 𝑃 ) ˜ 𝑄 = ( 𝑃 + 𝑃 ) 𝑄 = ( ˜ 𝑄 + ˜ 𝑄 ) . REFERENCES
P A Absil, Pierre-Yves Gousenbourger, Paul Striewski, and Benedikt Wirth. 2016. Dif-ferentiable Piecewise-Bézier Surfaces on Riemannian Manifolds.
SIAM Journal onImaging Sciences
9, 4 (Jan. 2016), 1788–1828.Adobe Inc. 2019.
Adobe Illustrator . https://adobe.com/products/illustratorB. Afsari. 2009.
Means and Averaging on Riemannian Manifolds . Ph.D. Dissertation.https://books.google.it/books?id=HUAhnQAACAAJAuthor Anonymous. 2020. This reference is anoymized for the purpose of blind review.Antoine Arnould, Pierre-Yves Gousenbourger, Chafik Samir, Pierre-Antoine Absil, andMichel Canis. 2015. Fitting Smooth Paths on Riemannian Manifolds: EndometrialSurface Reconstruction and Preoperative MRI-Based Navigation. In
GeometricScience of Information . Springer, Cham, Cham, 491–498.Dimitri P. Bertsekas. 1998.
Network optimization: continuous and discrete models . AthenaScientific.H Biermann, I Martin, F Bernardini, and D Zorin. 2002. Cut-and-paste editing ofmultiresolution surfaces.
ACM Transactions on Graphics
21, 3 (July 2002), 312–321.M. Camarinha, F. Silva Leite, and P. Crouch. 1995. Splines of class 𝐶 𝑘 on non-euclideanspaces . IMA Journal of Mathematical Control and Information
12, 4 (1995), 399–410.Keenan Crane, Fernando de Goes, Mathieu Desbrun, and Peter Schröder. 2013. DigitalGeometry Processing with Discrete Exterior Calculus. In
ACM SIGGRAPH 2013courses (Anaheim, California) (SIGGRAPH ’13) . ACM, New York, NY, USA, 126.Keenan Crane, Marco Livesu, Enrico Puppo, and Yipeng Qin. 2020. A Survey of Algo-rithms for Geodesic Paths and Distances. arXiv.org (2020). arXiv:2007.10430 [cs.GR]M.P. do Carmo. 1976.
Differential Geometry of Curves and Surfaces . Prentice-Hall.Tom Duchamp, Gang Xie, and Thomas Yu. 2018. Smoothing nonlinear subdivisionschemes by averaging.
Numerical Algorithms
77, 2 (Feb. 2018), 361–379.Nira Dyn, Ron Goldman, and David Levin. 2019. High order smoothness of non-linearLane-Riesenfeld algorithms in the functional setting.
Computer Aided GeometricDesign
71 (May 2019), 119–129.Nira Dyn and Nir Sharon. 2017. Manifold-valued subdivision schemes based on geodesicinductive averaging.
J. Comput. Appl. Math.
311 (Feb. 2017), 54–67.Gerald Farin. 2001.
Curves and Surfaces for CAGD: A Practical Guide (5th ed.). MorganKaufmann Publishers Inc., San Francisco, CA, USA.Ron Goldman. 2002.
Pyramid Algorithms: A Dynamic Programming Approach to Curvesand Surfaces for Geometric Modeling (1st ed.). Morgan Kaufmann Publishers Inc.,San Francisco, CA, USA.R.N. Goldman and T. Lyche. 1993.
Knot Insertion and Deletion Algorithms for B-SplineCurves and Surfaces . Society for Industrial and Applied Mathematics.Pierre-Yves Gousenbourger, Estelle Massart, and P A Absil. 2018. Data Fitting onManifolds with Composite Bézier-Like Curves and Blended Cubic Splines.
Journalof Mathematical Imaging and Vision (Dec. 2018), 1–27.Pierre-Yves Gousenbourger, Chafik Samir, and P A Absil. 2014. Piecewise-Bezier C1Interpolation on Riemannian Manifolds with Application to 2D Shape Morphing. In . IEEE, 4086–4091.Karsten Grove and Hermann Karcher. 1973. How to Conjugate C1-Close Group Actions.
Mathematische Zeitschrift
132 (1973), 11–20. http://eudml.org/doc/171906Philipp Herholz and Marc Alexa. 2019. Efficient Computation of Smoothed ExponentialMaps.
Computer Graphics Forum (Jan. 2019).M Hofer and H Pottmann. 2004. Energy-minimizing splines in manifolds.
ACMTransactions on Graphics (2004), 284.Y Jin, D Song, T Wang, J Huang, Y Song, and L He. 2019. A shell space constrainedapproach for curve design on surface meshes.
Computer-Aided Design
113 (2019),24–34.H Karcher. 1977. Riemannian center of mass and mollifier smoothing.
Communicationson Pure and Applied Mathematics
30, 5 (Sept. 1977), 509–541. Felix Knöppel, Keenan Crane, Ulrich Pinkall, and Peter Schröder. 2013. Globally optimaldirection fields.
ACM Transactions on Graphics
32, 4 (July 2013), 1–14.J M Lane and RF Riesenfeld. 1980. A theoretical development for the computer gener-ation and display of piecewise polynomial surfaces.
IEEE Transactions on PatternAnalysis and Machine Intelligence
2, 1 (1980), 35–46.D T Lee and F P Preparata. 1984. Euclidean shortest paths in the presence of rectilinearbarriers.
Networks
14, 3 (Sept. 1984), 393–410.A Lin and M Walker. 2001. CAGD techniques for differentiable manifolds. In
Interna-tional Symposium on Algorithms for Approximation . 36–43.D M Morera, P C Carvalho, and L. Velho. 2008. Modeling on triangulations withgeodesic curves.
The Visual Computer
24 (2008), 1025–1037.E Nava-Yazdani and K. Polthier. 2013. De Casteljau’s algorithm on manifolds.
ComputerAided Geometric Design
30, 7 (2013), 722–732.Lyle Noakes. 1998. Nonlinear corner-cutting.
Advances in computational Mathematics
8, 3 (1998), 165–177.L Noakes, G Heinzinger, and B Paden. 1989. Cubic splines on curved spaces.
IMAJournal of Mathematical Control and Information ´ zier curves on Riemannian manifolds and Lie groupswith kinematics applications. Journal of Mechanical Design
ACM SIGGRAPH 2018 Talks . ACM Press, New York, New York, USA, 1–2.Konrad Polthier and Markus Schmies. 1998. Straightest geodesics on polyhedral surfaces.In
Mathematical Visualization . Springer-Verlag, New York, 135–150.Tomasz Popiel and Lyle Noakes. 2007. Bézier curves and C2 interpolation in Riemannianmanifolds.
Journal of Approximation Theory
Computer aided geometric design
22, 7 (Oct. 2005), 693–709.Kaihuai Qin. 2000. General matrix representations for B-splines.
The Visual Computer
16, 3 (2000), 177–186. https://doi.org/10.1007/s003710050206D Salomon. 2006. Curves and surfaces for computer graphics. Springer.Chafik Samir, P A Absil, Anuj Srivastava, and Eric Klassen. 2011. A Gradient-DescentMethod for Curve Fitting on Riemannian Manifolds.
Foundations of ComputationalMathematics
12, 1 (April 2011), 49–73.R Schmidt. 2013. Stroke Parameterization.
Computer Graphics Forum
32, 2pt2 (May2013), 255–263.R Schmidt, C Grimm, and B Wyvill. 2006. Interactive decal compositing with discreteexponential maps.
ACM Transactions on Graphics (TOG)
25, 3 (2006), 605–613.Nicholas Sharp and Keenan Crane. 2020. You Can Find Geodesic Paths in TriangleMeshes by Just Flipping Edges.
ACM Trans. Graph. (I3D ’13)
Constructive Approximation
24, 3 (Nov. 2006), 289–318.Johannes Wallner and Nira Dyn. 2005. Convergence and analysis of subdivision schemeson manifolds by proximity.
Computer Aided Geometric Design
22, 7 (2005), 593–622.Johannes Wallner and Helmut Pottmann. 2006. Intrinsic subdivision with smoothlimits for graphics and animation.
ACM Transactions on Graphics
25, 2 (April 2006),356–374.Shi-Qing Xin and Guo-Jin Wang. 2007. Efficiently Determining a Locally Exact ShortestPath on Polyhedral Surfaces.
Comput. Aided Des.
39, 12 (Dec. 2007), 1081–1090.https://doi.org/10.1016/j.cad.2007.08.001Fujio Yamaguchi. 1988.
Curves and Surfaces in Computer Aided Geometric Design (1 ed.).Springer-Verlag Berlin Heidelberg.Qingnan Zhou and Alec Jacobson. 2016. Thingi10K: A Dataset of 10,000 3D-PrintingModels. arXiv preprint arXiv:1605.04797arXiv preprint arXiv:1605.04797