Clothoid fitting and geometric Hermite subdivision

We consider geometric Hermite subdivision for planar curves, i.e., iteratively refining an input polygon with additional tangent or normal vector information sitting in the vertices. The building block for the (nonlinear) subdivision schemes we propose is based on clothoidal averaging, i.e., averaging w.r.t. locally interpolating clothoids, which are curves of linear curvature. To this end, we derive a new strategy to approximate Hermite interpolating clothoids. We employ the proposed approach to define the geometric Hermite analogues of the well-known Lane-Riesenfeld and four-point schemes. We present numerical results produced by the proposed schemes and discuss their features.


Introduction
Linear subdivision schemes are widely used in various areas such as geometric modeling, multiscale analysis and for solving PDEs, and are rather well studied; references are for instance [CDM91,DL02,Han18,PR08].
In the last two decades also the interest in nonlinear subdivision schemes has significantly increased.One class of such schemes deals with scalar real valued data, but employs nonlinear averaging/prediction/interpolation techniques, e.g., [DY00, KVD02, CDM03, GVS09], making the corresponding schemes nonlinear.Motivation may be to get more robust estimators for processing data or to be able to deal with discontinuities or nonuniform data.Another class of nonlinear schemes addresses data which live in a nonlinear space, such as a Riemannian manifold or a Lie group, see for instance [WD05,XY09,Gro10,Wei10,Moo16].A third class considers data living in Euclidean space, typically of dimension 2 or 3, but the averaging rules are nonlinear to take account of their geometric characteristics.Such rules may be formulated in terms of angles, perpendicular bisectors, interpolating circles, and the like, rather than acting on each component separately, as linear schemes do.References of such geometric schemes include [AEV03, SD04, MDL05, DFH09, CHR13].
In contrast to linear schemes, for which a rather well established analysis is available, nonlinear schemes are less well understood, and there are ongoing efforts to devise new and to improve existing tools.However, due to the nonlinear nature and the diversity of the proposed schemes, not as general results as in the linear case can be expected.The investigation of a particular class of nonlinear schemes is likely to require an additional particular analysis component not covered by a general theory.Papers providing an analysis framework for geometric curve subdivision are [DH12,ERS15].The first reference derives sufficient conditions for a convergent interpolatory planar subdivision scheme to produce tangent continuous limit curves.The second reference deals with subdivision schemes which are geometric in the sense that they commute with similarities and derives a framework to establish C 1,α -and C 2,α -regularity of the generated limit curves.
In this paper, we present a family of geometric Hermite subdivision schemes for the generation of planar curves where the data to be refined are point-vector pairs, the latter serving as information on tangents or normals.Schemes refining point-vector pairs of that type were already suggested in [CJ07,LD16].The basic idea of these two approaches is to locally fit circles to the data and then to sample new points from them, but the specific methods are different.Also the strategies for determining new vectors are not the same so that the two schemes have significantly different shape properties.In contrast, the approach proposed here relies on clothoids, which are curves of linear curvature, rather than circles.The top row of Figure 1 shows the refinement of initial data consisting of two point-normal pairs by the so-called clothoid average, as described in this paper.The resulting curve is S-shaped and interpolates the initial data.Also the scheme of Chalmoviansky and Jüttler [CJ07], as shown in the middle row, interpolates the initial data, but has three turning points, leading to a less natural shape and a less optimal distribution of curvature.The scheme of Dyn and Lipovetsky [LD16], shown in the bottom row, produces a straight line (as it would for any choice of parallel normals).This is likely not to address the designer's intent, and it also does not interpolate the normal directions specified at the endpoints.The good performance of our scheme is related to the fact that it is not based on the reproduction of circles, which can be seen as the geometric analogue of quadratic polynomials, but on the reproduction of clothoids, which can be seen as the geometric analogue of cubic polynomials.Thus, it is able to mimic a much larger variety of shapes.
The idea of geometric planar spline fitting, i.e., fitting splines based on clothoids as building blocks can be traced through the literature for more than 50 years [BdB65,Sto82,Meh74,MW92,HL92,Coo93].Further recent contributions are for instance [BF18] presenting an alternative to [Sto82] for the computation of a C 2 interpolating clothoid spline.Hermite interpolation problems w.r.t.clothoids are for instance the topic of [BF15].We also point out [MW09] where the authors employ so called e-curves as approximative substitutes for clothoids in the context of Hermite interpolation.
Maybe the two striking reason why clothoid splines and approximations of them are still a topic in the literature are as follows: (i) the corresponding clothoid splines are rather expensive two compute; (ii) there are plenty of applications.Let us discuss these points in more details.Concerning (ii), for a long time clothoids have been used by route designers as transitional curves between straight lines and circular arcs, and between circular arcs of different radii; see for instance [MW90,Baa84].Nowadays, they are further used in connection with path planning for autonomous vehicles, e.g., [Ber15, AGRA15, BLR + 12] in computer vision and image processing [KFP03,BL15], in curve editing for design purposes [HEWF13] or representing hand-drawn strokes sketched by a user [MS09,BLP10].Concerning (i), many of the above applications are rather time critical and it is often important to have algorithms which are as fast as possible at a given (sometimes moderate) approximation quality.Solving clothoidal spline fitting problems up to high precision can be done rather fast [Sto82] but even this is sometimes too expensive or not needed.Instead, frequently faster strategies typically using some approximation are employed.For instance, [BLP10] locally fits clothoids and arcs as primitives and then optimizes w.r.t. a certain graph to obtain a global fit.
The paper [HEWF13] uses a variational approach iteratively inserting control points and optimizing them such as to generate a polyline with linear discrete curvature (approximating the clothoid segment; for details see also [SK00].)We mention that this approach uses subdivision (explicitly the analogue of corner cutting) for clothoid blending.Other approaches replace clothoids by approximating curve segments which are easier to handle, e.g., [MW04,MW09].We point out that the clothoidal subdivision schemes suggested here may be used as a computationally rather cheap alternative to clothoid splines and that they may be employed in the various applications discussed above.
The paper is organized as follows: After presenting some basic facts about two-point Hermite interpolation with clothoids in the next section, we consider its approximate solution in Section 3. The formula we propose is explicit, fast to evaluate, and yields a very small error for a large range of input data.In Section 4, this result is used to define the so-called clothoid average of a pair of points and corresponding tangent directions, which in turn serves as a building block for new families of geometric Hermite subdivision schemes.In particular, we obtain geometric Hermite analogues of the Lane-Riesenfeld schemes and the four-point scheme.In Section 5, we present results produced by the proposed schemes and discuss their features to illustrate the potential of the new method.As a first theoretical result, Section 6 establishes convergence and G 1 -continuity of the Lane-Riesenfeld-type algorithm of degree 1. Concluding remarks and an outlook are given in Section 7.

Two-point interpolation
In the following, points in the plane, and in particular images of planar curves, are considered as complex numbers.In this sense, let p : [0, 1] → C be a twice differentiable Figure 1: Midpoint refinement using the proposed clothoid average (upper row), the circle average of [CJ07] (middle row), and the circle average of [LD16] (lower row).We always display the points p j together with the normals n j = i exp(iα j ) at level .function parametrizing a planar curve which is regular in the sense that the velocity v := |p | vanishes nowhere.It is called uniform if the velocity is constant.According to the lifting lemma, the tangent vector can be expressed in the form p = v exp(iα) with a differentiable function α : [0, 1] → R, called the tangent angle of p.We write α = arg p for brevity, and assume throughout that α is unrolled suitably so that jumps are avoided.The curvature of p is κ := α /v = Im p p /v 3 .In the uniform case, on which we focus below, curve, velocity, and tangent angle are related by the formula The integral appearing here is abbreviated by exp(iα(s)) ds, I(α) := I(α, 1), so that p = p(0) + vI(α, •).A curve q : [0, 1] → C starting at q 0 := q(0) = 0 and ending at q 1 := q(1) ∈ R + is said to be in normal position.To tell this from the general case, we use the letters w = |q |, β = arg q , and λ = β /w to denote the velocity, tangent angle, and curvature, respectively.We have the relation Curves in general and normal position are linked by similarity.Denoting the secant between the endpoints of p by d := p 1 − p 0 and its angle by ϕ := arg d, we have the relations Let J = {0, 1}.Given points p J = (p 0 , p 1 ) with d := p 1 −p 0 = 0 and angles α J = (α 0 , α 1 ), the corresponding two-point Hermite interpolation problem is to find a curve p such that A special case of the above problem is to find a curve q such that q(j) = j, β(j) = β j , j ∈ J, for given β J .We recall that we use letters q, β to indicate that the sought curve is in normal position.Figure 2 illustrates the setting.While it is simple to specify curves merely satisfying these constraints, the challenge is to find solutions which are fair in some sense.For instance, it is a classical task to solve (3) in the space of clothoids, which are curves characterized by a linear curvature profile.In principle, this nonlinear problem is well understood and various more or less complicated methods for its numerical treatment are described in the literature, see for instance [WM09,BF15] and the references therein.Before we present our own approximate approach in the next section, we introduce some notation and basic facts.
Using the symbol P n to denote the space of polynomials of degree at most n over the unit interval, we define as the set of all uniform curves p : [0, 1] → C with curvature in P n .The corresponding subset of such curves in normal position is denoted by K + n .In particular, K 0 contains straight lines and circular arcs, while curves in K 1 \ K 0 are segments of clothoids.The tangent angle α = κ of a clothoid is a quadratic polynomial.For clothoids, the integral appearing in (1a) can be transformed to the so-called Fresnel integral F (x) := x 0 exp(iu 2 ) du, which does not possess a finite representation with respect to elementary functions.For later use we state the following lemma.The proof follows immediately form the Tait-Kneser Theorem.
In the following, for simplicity of wording, the term clothoid addresses not only segments of true clothoids, but also segments of circles and straight lines.In this sense, we consider the solution of (3) with clothoids, i.e., with curves p ∈ K 1 .By similarity according to (2), this problem can be reduced to (4) with data β J = α J − arg d, which are the angles between the secant d and the boundary tangent angles of p.Using the quadratic Lagrange polynomials with respect to the break points 0, 1/2, 1, the tangent angle of q can be written as Hence, given β J , the sought solution is characterized by the value β 1/2 = β(1/2) and the velocity w by means of (1b).The task is to find these two values.Figure 3 demonstrates that the solution of (4) is not unique.More precisely, it is known that for each pair β J there exists a countable family of solutions, but in applications, one is typically interested in curves avoiding excess rotation.For instance, if the boundary data β J are small in modulus, also the overall maximum β ∞ of tangent angles should be small.The following theorem guarantees existence of such a solution.
Theorem 2.2 There exists a smooth function F : U → R 2 , defined on some neighborhood U = (−u, u) 2 of the origin, with and the following property: Let then the tangent angle β and the velocity w define a solution q = wI(β, •) ∈ K + 1 of (4).In particular, I(β) = 1/w is real and positive.
Actually, it is possible to choose the domain U = (−π, π) 2 for F , covering almost all possible pairs of boundary tangent angles, but this is not needed here.
Proof.The idea is to parametrize the set solutions of (4) for varying β J as a twodimensional surface in R 4 and then to apply the implicit function theorem.Given α J ∈ R 2 , define the tangent angle α := α 0 2 0 + α 1 2 1 and the clothoid p := I(α, •) ∈ K 1 .Then q := p/p(1) ∈ K + 1 connects q(0) = 0 and q(1) = 1.Its velocity is w = 1/|p(1)|, and its tangent angle β = α − arg p(1) has values β j := β(j) = α j − arg p(1), j ∈ J, and β 1/2 := β(1/2) = − arg p(1).With these data, we define the surface Φ(α 0 , α 1 ) := [β 1/2 , w, β 0 , β 1 ] T .It is well defined and smooth in a neighborhood of the origin since are value and derivative of Φ at the origin.The lower (2 × 2)-submatrix of DΦ(0, 0) has determinant 2/3.Hence, by the implicit function theorem, there exists a neighborhood U of the origin and a smooth function F : U → R 2 such that [F (β J ), β J ] T defines a point on the trace of Φ for all β J ∈ U , corresponding to the clothoid solving (4).Value and derivative of F at the origin are given by One might suspect that employing functions α of the general form α = α 0 2 0 + α 1/2 2 1/2 + α 1 2 1 would define even more solutions of (4).To see that this is not true, let α : 1 be an angle function as considered above.Then the corresponding clothoids p := I(α, •) and p := I(α, •) are related by p = exp(iα 1/2 )p so that their normal forms coincide, Let p J , α J be boundary data for the general problem (3) with d = p 1 − p 0 = 0.If q is the solution of (4) for data β J := α J − arg d according to the preceding theorem, then p := p 0 + dq solves (3).With q = wI(β, •), an equivalent expression is which is independent of the velocity w.Hence, it suffices to know the first coordinate function The function f has the symmetry properties Hence, the second order derivatives vanish at the origin and we obtain the expansion

Approximate solution
Computing accurate solutions of the nonlinear problem (4) is possible, but the determination of β 1/2 = f (β 0 , β 1 ) requests more or less elaborate and/or computationally expensive numerical methods.Instead, we propose a good approximation f for later use with subdivision algorithms.More precisely, we seek a function f with the following properties: i) f : R 2 → R is a cubic polynomial.This choice combines modest complexity with sufficient flexibility to achieve good global approximation.
Before presenting our solution, let us explain the meaning of the angle defect as a measure for the quality of the approximation f .Theorem 3.1 Given Hermite data p J , α J with d = p 1 − p 0 = 0, let β J := α J − arg d and β := β 0 2 0 + f (β 0 , β 1 ) 2 1/2 + β 1 2 1 , as above.Then, analogous to (5), defines a clothoid interpolating the point data, i.e., p(0) = p 0 , p(1) = p 1 .At the boundaries, its tangent angle α differs from the prescribed values by the angle defect, Proof.Interpolation of the point data is trivial.For the tangent angle, we obtain To construct a suitable function f , we adopt the idea used in the proof of Theorem 2.2.For a large collection of angles α i J , i ∈ I, we compute points Φ(α i 0 , α i 1 ) = [β i 1/2 , w i , β i 0 , β i 1 ] T representing clothoids in normal position.Points for which one of the angles β i 0 , β i 1/2 , β i 1 lies outside the interval B are discarded as they correspond to clothoids with too large tangent angles, which are of little relevance for most applications, e.g., for design purposes.Denoting the index set of the remaining points by Ĩ, the task is to determine f such that f (β i 0 , β i 1 ) ≈ β i 1/2 for all i ∈ Ĩ.The ansatz for the cubic polynomial f satisfying ii) with symmetry properties according to iii) is Now, we can determine the unknown parameters f 0 , f 1 by standard approximation methods.The following choice was found when striving for a good compromise between the maximal angle defect δ(β J ) ∞ and simplicity of the coefficients.
It remains to show that this qualitative behavior is inherited by δ.To this end, we define the function 1 implies that the integral is nonzero so that ∆ is a smooth function over a compact domain and hence Lipschitz with some constant L. We obtain for β J ∈ U .For β J ∈ U , continuity of δ shows that the inequality remains valid with a possibly enlarged constant c 2 .The given numerical value for c 1 and c 2 can be verified by evaluation over a fine grid on B, see Figure 4.
In many applications, an error of less than a tenth of a degree for the interpolation of tangent angles will be acceptable so that the given approximation can be used directly.Moreover, the theorem shows that the approximation is exact for the symmetric case f (β 0 , −β 0 ) = f (β 0 , −β 0 ) = 0, whose solution is a segment of a circle or a straight line.
If, in the general case, higher accuracy is required, one step of the Newton iteration reduces the maximal angle defect to less than 5×10 −8 , and a second step to less than 5×10 −16 .

Hermite subdivision by the clothoid average
If a whole sequence of Hermite data points is to be interpolated, one might join always two of them by a clothoid and connect the segments to form a single composite curve.If the clothoids are determined exactly, their contact is G 1 , meaning that points and tangent directions of neighboring clothoids coincide at the junctions.If the clothoid is only approximated, as described in the preceding section, the contact is still continuous, but not G 1 in a strict sense.In any case, curvature is discontinuous, what may be insufficient for many design applications.Below, we propose different subdivision strategies to generate a (visually) smooth curve from a sequence p 0 j of points in C and corresponding tangent angles α 0 j , j ∈ Z.As building block we define the (approximate) clothoid average as follows: Definition 4.1 Given a pair of Hermite couples h j := (p j , α j ), j ∈ J, with d := p 1 − p 0 = 0, let p ∈ K 1 be the clothoid with endpoints p(j) = p j and tangent angles α(j) ≈ α j constructed by the approximation F of F according to Theorem 3.2.Then the approximate clothoid average of h 0 and h 1 at t ∈ R is defined by evaluation of p and the corresponding angle function α at t, and written as Sequences of Hermite couples are denoted by H := (h j ) j∈Z .Now, we generate sequences H 0 , H 1 , H 2 , . . .from given initial data H 0 by means of a binary subdivision operator S : H → H +1 , the rules of which are based on the clothoid average.
The simplest subdivision operator of that type is inserting a new Hermite couple always between two given ones, Formally, S 1 corresponds to the Lane-Riesenfeld algorithm of degree 1. Defining the averaging operator we can proceed in that direction and define Lane-Riesenfeld-type algorithms of degree n by applying (n − 1) rounds of averaging to the output of S 1 , For instance, the Chaikin-type algorithm, obtained for n = 2, explicitly reads Needless to say that this symbolic representation of a nonlinear process cannot be simplified through commutativity, associativity, or distributivity.As a further example, we define a family of interpolatory four-point schemes S 4 ω with tension parameter ω < 0, Of course, many other subdivision schemes, including schemes of arbitrary arity, can be constructed in this spirit.Let H = S H 0 = (h j ) j∈Z be the sequence of Hermite couples h j = (p j , α j ) at level ∈ N 0 .A common feature of the algorithms presented above is the reproduction of circles.That is, if there exists a midpoint m ∈ C and a radius r > 0 such that p j = m − ir exp(iα j ) for all j, then p +1 j = m − ir exp(iα +1 j ) for all j.Further, clothoids are almost reproduced.That is, if p ∈ K 1 is a clothoid with angle function α, and if p j = p(t j ), α j = α(t j ), j ∈ Z, for certain parameters t j , then there exist parameters t +1 j such that The quality of the approximation is determined by the magnitude of the angle defect according to the preceding section.Concerning convergence, the examples in the next section suggest that all algorithms presented here are G 1 -convergent in the following sense: There exists a curve p in C with tangent angle α which, respectively, are the limits of points and angles generated by the algorithm.More precisely, if j( ) is a sequence of integers such that t = lim →∞ 2 − j( ), then lim →∞ p j( ) = p(t) and lim This is analogous to standard Hermite subdivision, where the slope of the limit must coincide with the limit of slopes, and that is why we suggest to call the procedures introduced here Geometric Hermite subdivision.A proof of G 1 -convergence of the Lane-Riesenfeld-type algorithm S 1 of degree 1 is given in Section 6; a more general theory is currently developed and beyond the scope of this paper.

Numerical Experiments
In this section, we present numerical examples illustrating the shape properties of the Hermite subdivision algorithms S n and S 4 ω as introduced in the preceding section.Throughout, we use the same set of initial data H 0 .To avoid a special treatment of boundaries, it is assumed to be periodic, h 0 j = h 0 j+8 , j ∈ Z.All figures are structured as follows: On the left hand side, we see the initial points p 0 j and normals n 0 j := i exp(iα 0 j ) together with the polygon p j as obtained after = 8 rounds of subdivision.The middle figure shows the points p j for = 5 together with the corresponding normals n j := i exp(iα j ).On the right hand side, estimated curvature values κ j , = 8, are plotted versus the normalized chord length where σ is chosen such that the length of one loop equals 1.Specifically, the curvature value κ j is computed as the reciprocal radius of the circle interpolating the points p j−1 , p j , p j+1 .The Fresnel-integrals I(α, •) appearing in the definition of the clothoid average are computed using Gauss-Legendre quadrature with three nodes.
Figure 5 shows the Lane-Riesenfeld-type algorithm S 1 .By construction, it is interpolatory.The plots suggests that the limit is G 1 , i.e., free of kinks, and that the tangent angle of the limit equals the limit of tangent angles.The same observation is true for all subsequent cases.Curvature looks piecewise linear, as it would be the case when connecting always two consecutive points by a clothoid.In fact, the pieces are not exact clothoids due to the angle defect, but the deviation is very small.The uneven distribution of spikes in the middle figure indicates that the standard parametrization t j = j2 − → p j does not converge to a differentiable limit.
Figure 6 shows the Lane-Riesenfeld-type algorithm S 2 .The scheme is no longer interpolatory, even though the initial data points p 0 j are very close to the generated curve.The same is true for all schemes S n , n ≥ 2. Curvature looks continuous, but it still has certain imperfections.Thanks to the averaging step, the standard parametrization is now smoothed out, and we conjecture that it is C 1 .
Figure 7 shows the Lane-Riesenfeld-type algorithm S 3 .Now, the curvature distribution is free of artifacts so that the generated curve can be rated Class A according to the conventions of the automotive industry.However, there are certain spots where curvature seems to be not differentiable with respect to arc length.That is, the limit is G 2 , but not G 3 .The same seems to be true also for Lane-Riesenfeld variants of even higher degree.
Figure 8 shows the four-point scheme with weight ω = −1/18.It is interpolatory and seems to generate a G 2 -limit.Extended experiments show that the value ω = −1/18  6 G 1 -convergence of the scheme S 1 While the focus of this paper is on the construction of geometric Hermite subdivision schemes, we also want to outline a proof for the G 1 -convergence of the Lane-Riesenfeldtype algorithm S 1 .It is based on the results in [DH12], but we want to remark that the notion of G 1 used there is special.It would be good to have a link between the approach of Dyn/Hormann and standard theory, which calls a curve G 1 if it has a C 1 -reparametrization.
Convergence and smoothness are local properties of S 1 in the sense that the Hermite couples h j depend only on h 0 k , h 0 k+1 for k ≤ j2 − ≤ k + 1, corresponding to the interval [k, k + 1] of the standard parametrization.That is, convergence and smoothness can be studied for initial data H 0 which are periodic in the sense that h 0 j = h 0 j+k , j ∈ Z, for some k ≥ 2. This assumption avoids technical problems with unboundedness of sequences.Throughout, as before, H = (h j ) j∈Z , h j = (p j , α j ).Further, we define the   vectors d j := p j+1 − p j and the pairs of angles β j,J = (β j,0 , β j,1 ) by β j,0 := α j − arg d j , β j,1 := α j+1 − arg d j .
The Euclidean disk in R 2 with radius 3π/4 is denoted by B * := {β J : β J 2 ≤ 3π/4}.Theorem 6.1 Let H 0 be periodic.If d 0 j = 0 and β 0 j,J ∈ B * for all j ∈ Z, then the iterates H := S 1 H 0 define a sequence of polygons (p j ) j∈Z converging to a G 1 -limit in the sense of [DH12].Moreover, the angles α j converge to the arguments of d j ,  First, the sequence of maximal secant lengths is summable, implying that the sequence of polygons (p j ) j∈Z , converges to a continuous limit, see Theorem 3 and Proposition 4 in [DH12].Second, to address G 1 -continuity, we consider the exterior angles  Hence, also the sequence of maximal exterior angles is summable, implying that the limit curve is G 1 , see Theorem 18 in [DH12].Third, the final statement of the theorem is a simple consequence of α j − arg(d j ) = β j,0 .

Conclusion and Outlook
In this paper, we have proposed geometric Hermite subdivision schemes and we have demonstrated that they are a reasonable means for designing curves with prescribed tangents or normals.More precisely, we have first proposed an explicit strategy to approximate Hermite interpolating clothoids and used it to define the clothoid averages.Then we have used clothoidal averaging to define geometric Hermite subdivision schemes.Particular instances considered were the geometric Hermite analogues of the Lane-Riesenfeld schemes and of the four-point scheme.Examples demonstrate that these schemes yield very convincing results.Finally, we have presented some first smoothness results.More precisely, we obtain smoothness in the sense of [DH12] for the first order geometric Lane-Riesenfeld Hermite scheme.However, the notion of [DH12] is not standard; standard theory calls a curve G 1 if it has a C 1 reparametrization.It would be desirable to obtain related G 1 -or even G 2 -smoothness results for wider classes of schemes.This is an interesting topic of ongoing research.

Figure 2 :
Figure 2: Interpolation problem in general position (left) and normal position (right).