SSECOND ORDER SPIRAL SPLINES
LYLE NOAKES
Abstract.
Second order spiral splines are C unit-speed planar curves that can be used to interpolate a list Y of n + 1 pointsin R at times specified in some list T , where n ≥
2. Asymptotic methods are used to develop a fast algorithm, based on a pairof tridiagonal linear systems and standard software. The algorithm constructs a second order spiral spline interpolant for datathat is convex and sufficiently finely sampled. Introduction
Given a finite sequence T of real numbers T < T < . . . < T n , a pair V of vectors V , V n ∈ R m and a finite sequence Y of points Y , Y , . . . , Y n ∈ R m where m ≥
1, an interpolant of (
Y, V ) at T is defined to be a C curve y : [ T , T n ] → R m satisfying y ( T j ) = Y j for 0 ≤ j ≤ n , as well as the auxiliary conditions y (1) ( T ) = V and y (1) ( T n ) = V n on derivatives y (1) of y at T and T n . Then y is said to interpolate ( Y, V ) at T .A standard interpolant is the minimiser y of J ( y ) := (cid:90) T n T (cid:107) y (2) ( t ) (cid:107) dt where (cid:107) (cid:107) denotes the Euclidean norm, and y (2) is the second derivative. There is precisely one minimiser, namely the cubicpolynomial spline which is the C piecewise-cubic polynomial interpolant with knots at T . The cubic polynomial splineis easily calculated from T and ( Y, V ), by solving a tridiagonal system of linear equations [10]. It is much more difficult,both from a theoretical and practical point of view, to find C interpolants that are unit-speed , namely (cid:107) y (1) ( t ) (cid:107) = 1 for all t ∈ [ T , T n ].Unit-speed curves are significant for highway and railway design [19], [20], motion planning for robots [14], and path planningfor unmanned aerial vehicles and military aircraft [9], [18]. Whereas cubic polynomial splines always exist, an evidentnecessary condition for existence of a C unit-speed interpolant is(1) (cid:107) Y j − Y j − (cid:107) ≤ L j := T j − T j − for 1 ≤ j ≤ n. By an elastic spline we mean a critical point of the restriction of J to the space of unit-speed interpolants. An elastic splineis known to be precisely a C interpolant that is a track-sum y of elastica (elastic curves) y j : [ T j − , T j ] → R m . An elastica is defined as a critical point of J restricted to unit-speed curves y j where y j ( T j − ) , y (1) j ( T j − ) , y j ( T j ) , y (1) j ( T j ) are allprescribed in advance. Elastic splines are extensively studied, with well-developed algorithms for interpolation, as discussedin §
2. These algorithms for elastic splines require solutions of systems of nonlinear equations for variable parameters, andconsequently very significant effort, compared with the sparse linear systems needed for interpolation by (non-unit-speed)polynomial splines.Besides elastica there are simpler classes of unit-speed track-summands called polynomial spirals [14]. The well-studied classof clothoidal splines , whose curvatures are C and piecewise-affine, is insufficiently rich for interpolation when T is given inadvance. In the present paper we suppose that the inequalities (1) hold strictly, and consider the larger class of second orderspiral splines , namely interpolants that are C track-sums y θ of second order generalised Cornu spirals , whose curvatures are C and piecewise-quadratic in the parameter t .In the present paper the data ( Y, V ) is generated by unknown strictly convex curves x with suitably bounded derivatives,and sufficiently fine sampling specified by T . The setting is described in detail in §
3. In § § C cubic polynomial spline ˆ θ : [ t , T n ] → R . Our construction echoes the well-known algorithm forpolynomial cubic splines, where there is a single tridiagonal system [10]. As with the classical algorithm, our linear systemsare well-conditioned and very quick to solve. Our construction is more than twice as complicated as the classical algorithm,and requires more care to implement, but the description in § § § § § § (cid:15) > θ = θ + O ( (cid:15) ) where y θ is the desired interpolant. This immediately gives rise to some approximate interpolants by spiral splines, as illustratedin Examples 1, 2, 3. More significantly, as demonstrated later in §
10, availability of a suitable initial guess is absolutely keyto the nonlinear problem of determining θ and thereby y θ . The principle is illustrated in Examples 4, 5, 6. Date : January 10, 2019. The term elastic spline is sometimes used differently, to mean an interpolant where T is not specified in advance [5]. The length of y j is then T j − T j − . Sometimes the term elastica has a different meaning, where the length of y j is not prescribed in advance[5]. It seems that the plural of elastica is elastica. The nonconvex case is more complicated, with potentially 2 n interpolants. a r X i v : . [ m a t h . NA ] M a y LYLE NOAKES More Background on Unit-Speed Interpolants
A natural class of unit-speed interpolants is that of the unit-speed reparameterised polynomial splines , studied in [24], [11],[21] where there are algorithms for efficient computation. Whereas a polynomial spline can be found to interpolate (
Y, V )at T , the unit-speed reparameterisation interpolates ( Y, V ) at a different set of parameter values. This does not work wellwhen T is given in advance, as in our situation. Similar difficulties arise with clothoidal splines [23], [4], [3], [19], [8] andtheir C generalisations [18]. The algorithms in these papers are of significant interest, but do not apply when T is alreadyprescribed.Interpolation by unit-speed curves can also be performed using elastic splines [16], [12], namely critical points of therestriction of the functional J to the space of unit-speed C interpolants. An elastic spline turns out to be a C track-sum y of elastica y j : [ T j − , T j ] → R m , satisfying the auxiliary end conditions. Elastica (elastic curves) have a long and interestinghistory [15], [5], [1], and their study reduces to the case where m = 3. The elastica are completely known in terms ofelliptic functions [22], with simplifications for m = 2, which is the case of interest for the present paper. Although elasticsplines (sometimes called nonlinear splines or true splines ) are highly regarded as interpolants, they are less widely usedthan cubic polynomial splines. This is because the interpolation conditions for elastic splines require the solution of a systemof nonlinear equations that is even more complicated than for clothoidal splines. However, the relative sophistication ofelastica compared with cubic polynomials is not the main difficulty in computing elastic splines. As with most nonlinearproblems, indeed also with our present task of finding spiral spline interpolants, the main difficulty is the construction of asuitable initial guess.The algorithms in [16], [12] and [13] are for sliding elastic splines , where T is not given in advance. Another condition thatis sometimes imposed is that the total time T n − T should be fixed or, equivalently, the length of the spline is prescribedin advance. In the present paper the entire sequence T is assumed to be given in advance, by analogy with standardinterpolation by cubic polynomial splines.What is different from cubic polynomial splines is that our interpolants are required to be unit-speed. Because of thisadditional requirement, once T is given, fixing the rest of the sequence T is equivalent to prescribing the lengths L j of theinterpolant between all consecutive data points. In these circumstances, the main contribution to interpolation by elasticsplines is Theorem 4.1 of [17], where a nonlinear system of equations is given for determining elastic splines in great generality.As usual with nonlinear problems, an initial guess is required in order to get started. Unfortunately [17] says very little abouthow this might be constructed. The present paper contributes something in this respect, although our focus is primarilyon spiral splines, whereas [17] is concerned with elastic splines. Our guess hinges on the construction of ˆ θ in § §
5. In § θ to θ .We also mention the discrete elastic splines of [6], where a discrete analogue of J is optimised with respect to a variable finitesequence of points approximating y . In effect optimisation of J with respect to y is replaced by a large finite-dimensionaloptimisation, whose outcome depends critically on an unspecified initial guess in a high dimensional space.In summary, reparameterised polynomial splines and clothoidal splines are unsuitable for interpolating ( Y, V ) by T . Elasticsplines come with non-negligible computational issues, especially the generation of suitable initial guesses.3. Convex Generators and Admissible Samplings
Let n ≥ T < T n , 0 < C < B , and B i where 3 ≤ i ≤
5. A unit speed C ∞ curve x : [ T , T n ] → R is called a convex generator when C < (cid:107) x (2) (cid:107) ∞ ≤ B and (cid:107) x ( i ) (cid:107) ∞ ≤ B i for i = 3 , ,
5. For brevity we write X for the space X T ,T n ,C,B ,B ,B ,B of all convex generators.Let 0 < A m < A M be given, together with (cid:15) > C, B , B , B that define X . We define T (cid:15) to be the set of finite sequences T of reals T < T < T < . . . < T n ≤ b such that A m (cid:15) < L j < A M (cid:15) for all1 ≤ j ≤ n .A pair ( Y, V ) is called admissible with respect to T ∈ T (cid:15) when Y = ( x ( T ) , x ( T ) , . . . , x ( T n )) and V = ( x (1) ( T ) , x (1 )( T n ) forsome convex generator x ∈ X . Then, for 1 ≤ j ≤ n , set q j := ( Y j − Y j − ) /L j , r j := (cid:107) q j (cid:107) , and k j := (cid:113) − r j ) L j r j . Lemma 1.
For small (cid:15) > , r j is bounded away from , and k j ≥ C + O ( (cid:15) ) for all ≤ j ≤ n . Proof:
We have x ( T j ) − x ( T j − ) = x (1) ( T j − ) L j + x (2) ( T j − ) L j x (3) ( T j − ) L j O ( (cid:15) ) where x ∈ X . In particular L j r j =1 + O ( (cid:15) ), and r j is bounded away from 0 when (cid:15) > A clothoidal spline is a unit-speed C planar curve t (cid:55)→ y ( t ) whose curvature is C and piecewise-affine in t . Elastic splines and elastica mean different things depending on the context [17]. In the present paper, elastic splines are pinned and clamped ,and elastica are fixed-length . In [16], T is not given in advance and the elastic splines are sliding . On the other hand, the study of elastic splines does not reduce in this way. The objection raised at the end of [16] does not really apply when T is given. As noted on p.184 of [12], clothoidal splines are sometimes used to construct initial guesses for the computation of elastic splines. The convexity assumption built into admissibility is the reason that we are able to find a unique estimate ˆ θ (otherwise there would be many). ECOND ORDER SPIRAL SPLINES 3
Because (cid:107) x (1) (cid:107) = 1, we also have (cid:104) x (1) , x (2) (cid:105) = 0 and (cid:104) x (1) , x (3) (cid:105) = −(cid:107) x (2) (cid:107) , where (cid:104) , (cid:105) is the Euclidean inner product.Therefore, and by convexity, r j L j = L j − (cid:107) x (2) ( T j − ) (cid:107) L j
12 + O ( (cid:15) ) = ⇒ − r j = (cid:107) x (2) ( T j − ) (cid:107) L j
12 + O ( (cid:15) ) ≥ C L j
12 + O ( (cid:15) ) . So for small (cid:15) , 0 < r j < k j > − r j ) /L j ≥ C + O ( (cid:15) ). (cid:3) From admissibility the Y j are separated by no more than A M (cid:15) , and lie in the closed disc of radius nA M (cid:15) centred on Y .Assuming (cid:15) is sufficiently small for the conclusion of Lemma 1 to hold, write q j = r j (cos ω j , sin ω j ). We also require ω ∈ ( − π, π ), and ω j = ω j − + O ( (cid:15) ) for 2 ≤ j ≤ n , so that the ω j are uniquely determined. By admissibility, V = (cos ν , sin ν )and V n = (cos ν n , sin ν n ) where ν = ω + O ( (cid:15) ) and ν n = ω n + O ( (cid:15) ).Let θ : [ T , T n ] → R be a C cubic spline with knots T j such that θ ( T ) = ν and θ ( T n ) = ν n . For 1 ≤ j ≤ n and all s ∈ [0 , L j ]write θ j ( s ) := θ ( T j − + s ) = a j + b j s + c j s + d j s . Besides the two auxiliary end conditions,(2) a = ν , a n + b n L n + c n L n + d n L n = ν n , there are 2 n − θ to be C , namely a j +1 = a j + b j L j + c j L j + d j L j , (3) b j +1 = b j + 2 c j L j + 3 d j L j , (4)where 1 ≤ j ≤ n −
1. So we have 2 n affine equality constraints on 4 n coefficients a j , b j , c j , d j .The C spline θ determines a second order spiral spline y θ : [ T , T n ] → R given by(5) y θ ( t ) := Y + (cid:90) tT (cos θ ( s ) , sin θ ( s )) ds. We see that y θ is unit-speed with y (1) ( T ) = V and y (1) ( T n ) = V n . We are interested in finding θ such that y θ is aninterpolant of ( Y, V ) at T , namely for 1 ≤ j ≤ n ,(6) (cid:90) L j (cos θ j ( s ) , sin θ j ( s )) ds = Y j − Y j − . Condition (6) amounts to 2 n non-affine conditions on the 4 n coefficients, making 4 n conditions in total.4. Tridiagonal Systems
Given (
Y, V ) admissible with respect to T ∈ T (cid:15) where (cid:15) > n × n tridiagonal matrix ˆT := L − L . . . − L L + L ) − L . . . − L L + L ) − L . . .
00 0 − L L + L ) − L . . . − L n − L n − + L n − ) − L n − − L n − L n − + 8 L n and, if n ≥
3, the tridiagonal matrix ˜T := L L . . . L L + L ) L . . . L L + L ) L . . .
00 0 L L + L ) L . . . L n − L n − + L n − ) L n − L n − L n − + 3 L n / . Looking forward, admissibility is needed to guarantee the conclusions of Theorem 1, but may be troublesome to check in applications. Insuch cases (including the examples of the present paper) Theorem 1 may be taken to assert that when its conclusions are false, the sampling of x is too irregular or too sparse. The cure, which seems to be rarely needed, is to sample more regularly and more often from x . LYLE NOAKES
Solve(7) ˜T ˜ b = ˜ R := 6 ω − ν ω − ω ω − ω ω − ω ...... ω n − − ω n − (3 ω n − ν n ) / − ω n − for an n -dimensional column vector ˜ b .Taking σ to be the sign of the determinant of the 2 × V ... Y − Y ], set ρ j := σk j (cid:115) − k j L j − (˜ b j +1 − ˜ b j ) k j for 1 ≤ j ≤ n − , and(8) ρ n := σk n (cid:118)(cid:117)(cid:117)(cid:116) − k n L n − k n (cid:32) ω n − ν n L n + ˜ b n (cid:33) . (9)Then solve(10) ˆT ˆ b = ω − ν ) − L ρ ω − ω ) − L ρ + L ρ )24( ω − ω ) − L ρ + L ρ )24( ω − ω ) − L ρ + L ρ )......24( ω n − − ω n − ) − L n − ρ n − + L n − ρ n − )24(2 ω n + ν n − ω n − ) − L n − ρ n − + 4 L n ρ n ) for an n -dimensional column vector ˆ b . Setˆ c j := − b j − b j +1 L j + 5 ρ j L j for 1 ≤ j ≤ n − , (11) ˆ d j := 5(ˆ b j + ˆ b j +1 )6 L j − ρ j L j for 1 ≤ j ≤ n − , (12) ˆ c n := 6( ω n − ν n ) L n − b n L n + 5 ρ n L n , (13) ˆ d n := − ω n − ν n )3 L n + 10ˆ b n L n − ρ n L n , (14)and finally ˆ a := ν with(15) ˆ a j +1 := ˆ a j + ˆ b j L j + ˆ c j L j + ˆ d j L j for 1 ≤ j ≤ n − . Then, for all 1 ≤ j ≤ n , define ˆ θ j : [0 , L j ] → R by ˆ θ j ( s ) := ˆ a j +ˆ b j s + ˆ c j s + ˆ d j s . A C cubic polynomial spline ˆ θ : [ T , T n ] → R is then given by ˆ θ ( T n ) := ˆ θ n ( L n ), and ˆ θ ( t ) := ˆ θ j ( t − T j − ) for t ∈ [ T j − , T j ). Theorem 1. ˆ θ ( T ) = ν , ˆ θ ( T n ) = ν n and, for small (cid:15) > , (cid:107) θ − ˆ θ (cid:107) ∞ = O ( (cid:15) ) . Corollary 1. y ˆ θ satisfies the auxiliary conditions (2), and (cid:107) y θ − y ˆ θ (cid:107) ∞ = O ( (cid:15) ) . (cid:3) To prove Theorem 1 it suffices to show that, for all 1 ≤ j ≤ n , a j = ˆ a j + O ( (cid:15) )(16) b j = ˆ b j + O ( (cid:15) )(17) c j = ˆ c j + O ( (cid:15) )(18) d j = ˆ d j + O ( (cid:15) ) . (19)We now prove some asymptotic lemmas, for use in § § § ECOND ORDER SPIRAL SPLINES 5 Approximation and Interpolation
By Corollary 1, the C second order spiral spline y ˆ θ satisfies y ˆ θ ( T ) = Y , y (1)ˆ θ ( T ) = V , y (1)ˆ θ ( T n ) = V n and y ˆ θ )( T j ) = Y j + O ( (cid:15) ) for 1 ≤ j ≤ n . Usually y ˆ θ ( T j ) (cid:54) = Y j for 1 ≤ j ≤ n . So y ˆ θ is only approximately an interpolant of ( Y, V ).Because of the way y θ is constructed from θ , interpolation errors at T j tend to accumulate as j increases. This can be avoidedby using θ differently to define ¯ y θ : [ T , T j ] → R by ¯ y θ ( T ) = Y and for 1 ≤ j ≤ n ,¯ y θ ( t ) := Y j − + (cid:90) tT j − (cos θ ( s ) , sin θ ( s )) ds for T j − < t < T j . By construction ¯ y θ ( T ) = Y , and ¯ y θ ( T − j ) = Y j for all j , regardless of θ . However unless y θ interpolates ( Y, V ) exactly at T the left-interpolant ¯ y θ is not even a second order spiral spline, because ¯ y θ is only left-continuous but not continuous at T j for 1 ≤ j ≤ n − y ˆ θ ( T j ) = Y j + O ( (cid:15) ) for 1 ≤ j ≤ n , and ¯ y (1)ˆ θ ( T ) = V , ¯ y (1)ˆ θ ( T n ) = V n . So the left-interpolant ¯ y ˆ θ approximately interpolates ( Y, V ) at T . Because the errors in ¯ y ˆ θ ( T j ) ≈ Y j do not accumulate as j increases, ¯ y ˆ θ is better than y ˆ θ fordiagnostics. Example 1.
Set n = 10 with T = (0 , . , . , . , . , . , . , . , . , . , . , and let the convex generatorbe the second order spiral x : [0 , . → R given by x (0) = (0 , and x (1) ( t ) = (cos ψ ( t ) , sin ψ ( t )) with ψ ( t ) = t + t + t .This generates data ( Y, V ) where Y = ((0 , , (0 . , . , (0 . , . , (0 . , . , (0 . , . , (0 . , . , (0 . , . , (0 . , . , (0 . , . , (0 . , . , (0 . , . and V = ((1 , , ) , (0 . , . . Computation of the n coefficients for ˆ θ took 0.000918 seconds in Mathematica ona 2015 a 2.2GHZ MacBook Air with 8 GB RAM. The plot of y ˆ θ is shown in Figure 1, together with the data. Although ψ is just a cubic polynomial, the C cubic spline ˆ θ is not exactly ψ , as can be seen from the failures evident in Figure 1 of y ˆ θ to interpolate at T , T , T , T . However y ˆ θ is not a bad initial guess for an interpolant. Finer sampling would improve it.Coarser sampling would make it worse. (cid:3) In Example 1 the convex generator x is actually a second order spiral. Similar results are obtained when the data is generatedby a more complicated convex curve. Example 2.
Take T = (0 , . , . , . , . , . , . , . , . , . , . and let x : [0 , . → R be given by x (0) =(0 , and x (1) ( t ) = e t + t (1 + sin(5 t ) / . Using x to generate data ( Y, V ) , Mathematica takes 0.000998 seconds to compute ˆ θ . The plot of y ˆ θ is shown in Figure 2, together with the data. In Figure 2, y ˆ θ almost interpolates the data, except near T , T , T , T . (cid:3) Example 3.
Whereas in Example 2 y ˆ θ is nearly acceptable, the quality of the approximate interpolant falls away very quickly ifwe use define ( Y, V ) in the same way but over a larger interval [ T , T n ] . Here we take T = (0 , . , . , . , . , . , . , . , . , . , . ,and x is given by the same formula over the larger interval [0 , . . Computation of ˆ θ took . seconds, and the plotof y ˆ θ is shown in Figure 3, together with the data. We see that y ˆ θ is very far from interpolating, except at T , T , T . Eventhough the curve nearly passes through the terminal point Y , it arrives too soon and overshoots. This poor performance canbe corrected by more refined sampling or, as illustrated in Example 6, by the more sophisticated method of §
10 which startswith ˆ θ . (cid:3) Some Asymptotic Lemmas
Recall θ j ( s ) = a j + b j s + c j s + d j s where s ∈ [0 , L j ] for 1 ≤ j ≤ n , and Y, V, T are given with (
Y, V ) is admissible withrespect to T . Using Mathematica to calculate the Taylor polynomial ˆ F j ( s ) of degree 4 of F j ( s ) := (cid:82) s (cos u, sin u ) du ) about s = 0, we find Lemma 2.
For ≤ j ≤ n , ˆ F j ( L j ) = L j (cid:20) cos a j − sin a j sin a j cos a j (cid:21) (cid:20) α j β j (cid:21) where α j := 1 − b j L j − b j c j L j b j − c j − b j d j ) L j and β j := b j L j c j L j − ( b j − d j ) L j − b j c j L j . (cid:3) From Lemma 2, L j r j (cid:20) cos ω j sin ω j (cid:21) = L j q j = Y j − Y j − = F j ( L j ) = L j (cid:20) cos a j − sin a j sin a j cos a j (cid:21) (cid:20) α j β j (cid:21) + O ( (cid:15) ) = ⇒ We have chosen T in order to sample from x somewhat more frequently where curvature is large. Finer sampling always helps. LYLE NOAKES Y Y Y Y Y Y Y Y Y Y Y V V Figure 1. y ˆ θ for Example 1 r j = α j + β j + O ( (cid:15) )(20) (cos( a j − ω j ) , sin( a j − ω j )) = ( α j , β j ) (cid:107) ( α j , β j ) (cid:107) + O ( (cid:15) ) . (21) Lemma 3. ( b j + c j L j + 910 d j L j ) + c j L j
15 = k j (1 − k j L j
20 ) + O ( (cid:15) ) for ≤ j ≤ n, (22) a j + b j L j c j L j d j L j ω j + O ( (cid:15) ) for ≤ j ≤ n, (23) b j +1 L j +1 + b j L j c j +1 L j +1 + 2 c j L j d j +1 L j +1 + 3 d j L j ω j +1 − ω j + O ( (cid:15) ) for ≤ j ≤ n − . (24) Proof:
Using Mathematica to substitute for r j from (20) in k j and then expand in powers of L j , we obtain(25) ( b j + c j L j + 910 d j L j ) + ( b j
20 + c j
15 ) L j = k j . ECOND ORDER SPIRAL SPLINES 7 Y Y Y Y Y Y Y Y Y Y Y V V - - - Figure 2. y ˆ θ for Example 2where 1 ≤ j ≤ n . In particular, b j = k j + O ( (cid:15) ). Substituting b j = k j + O ( (cid:15) ) back in (25) gives (22). For (23) use Mathematicato Taylor expand the argument of the right hand side of (21) in powers of L j . For (24) substitute in (3) for a j , a j +1 from(23). (cid:3) Estimating the b j to O ( (cid:15) )From (24), for 1 ≤ j ≤ n − b j +1 L j +1 + b j L j c j +1 L j +1 + 2 c j L j ω j +1 − ω j + O ( (cid:15) ) . From (23) with j = 1, and from (2),(27) c L = 3( ω − ν ) L − b O ( (cid:15) ) . From (4) for 1 ≤ j ≤ n − c j L j = b j +1 − b j O ( (cid:15) ) . From (23) with j = n , and from (2),(28) ω n − b n L n − c n L n a n + O ( (cid:15) ) = ν n − b n L n − c n L n = ⇒ c n L n = − ω n − ν n )2 L n − b n O ( (cid:15) ) . LYLE NOAKES Y Y Y Y Y Y Y Y Y Y Y V V - - - Figure 3. y ˆ θ for Example 3Substituting for c L in (27), and for c L , c L , . . . , c n L n in (26) with 1 ≤ j ≤ n −
1, we obtain the system2 L b + L b = 6( ω − ν ) + O ( (cid:15) ) ,L j − b j − + 2( L j − + L j ) b j + L j +1 b j +1 = 6( ω j − ω j − ) + O ( (cid:15) ) for 2 ≤ j ≤ n − ,L n − b n − + (2 L n − + 3 L n b n = 3(3 ω n − ν n − ω n − ) + O ( (cid:15) ) , of n linear equations for b := ( b , b , . . . , b n ), namely(29) ˜ T b = ˜ R + O ( (cid:15) ) = ˜ T ˜ b + O ( (cid:15) )where ˜ T is the n × n tridiagonal matrix and ˜ R , ˜ b are the n -dimensional vectors all defined in § Lemma 4. b = ˜ b + O ( (cid:15) ) . Proof:
Let D ˜ T be the diagonal n × n matrix with the same diagonal entries as ˜ T . Because A m (cid:15) < L j < A M (cid:15) , we see that (cid:107) D − T (cid:107) ∞ ≤ / (2 A m (cid:15) ). By (29), ˜ T ( b − ˜ b ) = O ( (cid:15) ), and therefore (cid:107) D − T ˜ T ( b − ˜ b ) (cid:107) = O ( (cid:15) ) . Because ˜ T is strictly diagonally dominant by rows with dominance factor 1 / D − T ˜ T has condition number ≤ (cid:3) ECOND ORDER SPIRAL SPLINES 9 Proof of (17)
By (4) and Lemma 4, c j L j = b j +1 − b j O ( (cid:15) ) = ˜ b j +1 − ˜ b j O ( (cid:15) )for 1 ≤ j ≤ n −
1. Similarly by (28) and Lemma 4, c n L n = 3( ω n − ν n )2 L n − b n O ( (cid:15) ). Substituting for the c j L j in (22),( b j + c j L j + 910 d j L j ) = k j (1 − k j L j
20 ) − (˜ b j +1 − ˜ b j )
60 + O ( (cid:15) ) for 1 ≤ j ≤ n − , and( b n + c n L n + 910 d n L n ) = k n (1 − k n L n
20 ) − (cid:32) ω n − ν n L n + ˜ b n (cid:33) + O ( (cid:15) ) , where the right hand sides are known to O ( (cid:15) ). By Lemma 1, for (cid:15) small the k j are bounded away from 0 for all 1 ≤ j ≤ n .Therefore b j + c j L j + 910 d j L j = σ j | ρ j | + O ( (cid:15) ) = σ j k j + O ( (cid:15) )where σ j = ±
1, and the ρ j are defined as in (8), (9). In particular, by Lemma 1 σ j = sign( b j ) for 1 ≤ j ≤ n . By (4), b j = b j − + O ( (cid:15) ) for 2 ≤ j ≤ n . Because the k j are bounded away from 0 so are the | b j | . So the σ j are all equal. Recall thedefinition of σ in §
4, namely σ := sign(det[ V ... Y − Y ]). Lemma 5.
For (cid:15) > small, we have σ j = σ where ≤ j ≤ n . Proof:
By (4), b = b + O ( (cid:15) ). Then, from the first equation in system (7),3 b = 2 b + b + O ( (cid:15) ) = 2˜ b + ˜ b + O ( (cid:15) ) = 6( ω − ν ) L + O ( (cid:15) ) = ⇒ σ = sign( ω − ν ) , where we use σ = sign( b ). Now V = x (1) ( T ) and Y j = x ( T j ) for 0 ≤ j ≤ n , where x : [ T , T n ] → R is a convex generator.Write x (1) ( t ) = (cos ψ ( t ) , sin ψ ( t )) where ψ is C and ψ ( T ) = ν . By convexity, | ψ (1) ( T ) | = (cid:107) x (2) ( T ) (cid:107) > C >
0. Taylorexpanding x to second order in t about T , and recalling the definition of ω j , L r (cos ω , sin ω ) = Y − Y = L (cos ν , sin ν ) + L ψ (1) ( T )( − sin ν , cos ν ) + O ( (cid:15) ) = ⇒ L r sin( ω − ν ) = L ψ (1) ( T ) + O ( (cid:15) ) = det[ V ... Y − Y ] + O ( (cid:15) ) = ⇒ sign( ω − ν ) = σ. So σ n = σ n − = . . . = σ = σ . (cid:3) From Lemma 5 we now have, for 1 ≤ j ≤ n ,(30) b j + c j L j + 910 d j L j = ρ j + O ( (cid:15) ) . By (4), (30), for each 1 ≤ j ≤ n −
1, 2 c j L j + 3 d j L j = b j +1 − b j + O ( (cid:15) ) ,c j L j + 910 d j L j = ρ j − b j + O ( (cid:15) ) , where the right hand sides are presently known only with O ( (cid:15) ) errors. Nevertheless, solving these systems with 1 ≤ j ≤ n − c j L j = − b j − b j +1 ρ j O ( (cid:15) ) , (31) d j L j = 5( b j + b j +1 )6 − ρ j O ( (cid:15) ) . (32)By (23) with j = n and (2), and by (30) with j = n ,8 c n L n + 9 d n L n = 12( ν n − ω n ) L n − b n + O ( (cid:15) ) , c n L n + 9 d n L n = 10 ρ n − b n + O ( (cid:15) ) . Solving, c n L n = 6( ω n − ν n ) L n − b n + 5 ρ n + O ( (cid:15) ) , (33) d n L n = − ω n − ν n )3 L n + 10 b n − ρ n O ( (cid:15) ) . (34)In particular, substituting for c L , d L in (23) with j = 1,3 L b − L b = 24( ω − ν ) − L ρ + O ( (cid:15) ) . Substituting for c j L j , d j L j where 1 ≤ j ≤ n in (24), gives n − b , namely − L j b j + 3( L j + L j +1 ) b j +1 − L j +1 b j +2 = 24( ω j +1 − ω j ) − L j ρ j + L j +1 ρ j +1 ) + O ( (cid:15) ) for 1 ≤ j ≤ n − , and − L n − b n − + (9 L n − + 8 L n ) b n = 24(2 ω n + ν n − ω n − ) − L n − ρ n − + 4 L n ρ n ) + O ( (cid:15) ) . So ˆ T b = ˆ R + O ( (cid:15) ) and (17) follows, in the same way as for Lemma 4. (cid:3) Proof of Theorem 1
As noted at the end of §
4, now that (17) is proved it suffices to verify (16), (18), (19). Using (17) in (31), (35), then in(33), (35), c j L j = − b j − b j +1 ρ j O ( (cid:15) ) = − b j − b j +1 ρ j O ( (cid:15) ) = ˆ c j L j + O ( (cid:15) ) for 1 ≤ j ≤ n − ,d j L j = 5( b j + b j +1 )6 − ρ j O ( (cid:15) ) = 5(ˆ b j + ˆ b j +1 )6 − ρ j O ( (cid:15) ) = ˆ d j L j + O ( (cid:15) ) for 1 ≤ j ≤ n − ,c n L n = 6( ω n − ν n ) L n − b n + 5 ρ n + O ( (cid:15) ) = 6( ω n − ν n ) L n − b n + 5 ρ n + O ( (cid:15) ) = ˆ c n L n + O ( (cid:15) ) ,d n L n = − ω n − ν n )3 L n + 10 b n − ρ n O ( (cid:15) ) = − ω n − ν n )3 L n + 10ˆ b n − ρ n O ( (cid:15) ) = ˆ d n L n + O ( (cid:15) ) , according to definitions (11), (12), (13), (14) in §
4. This proves (18), (19). Finally by (2) a = ν = ˆ a , and by (3) a j +1 = a j + b j L j + c j L j + d j L j = ˆ a j + ˆ b j L j + ˆ c j L j + ˆ d j L j + O ( (cid:15) ) = ˆ a j +1 + O ( (cid:15) ) according to (15) for 1 ≤ j ≤ n −
1. So(16) is proved by induction. (cid:3)
Closing Gaps
As noted in §
5, the C unit-speed curve y ˆ θ does not interpolate ( Y, V ) exactly at T . Usually there are O ( (cid:15) ) errors inthe interpolation conditions. The left-interpolant ¯ y ˆ θ is not much better, except to illustrate the magnitudes and locations oferrors, because ¯ y ˆ θ is not even C .When sampling is sufficiently fine, namely for sufficiently small (cid:15) >
0, the curves y ˆ θ and ¯ y ˆ θ are almost indistinguishable, andeither of these unit-speed approximate interpolants might be used in practice. This is already true for some parts of y ˆ θ inExample 1 and for most of y ˆ θ in Example 2. However y ˆ θ is seriously problematic in Example 3.If finer sampling is inconvenient or impossible then, after finding ˆ θ , there is a second calculation we can do to correct thekinds of errors observed in Example 3, and also to a smaller extent in Examples 1, 2.We rewrite the conditions for a second order spiral spline interpolant as a nonlinear system of equations, then solve numericallyusing ˆ θ to generate an initial guess. This is easier than finding ˆ θ and simpler to code, but not so remarkably fast. It isunlikely to fail except when sampling is so coarse that ˆ θ does not give a good initial guess. Even in the case of Example 3the rough estimate ˆ θ is enough to get started.Given ( Y, V ) admissible with respect to T , first calculate ˆ θ as in §
4. Then, for u, v ∈ R n , define a j , b j , c j , d j ∈ R for 1 ≤ j ≤ n by c j := u j and d j := v j for all 1 ≤ j ≤ n , and a := ν . Using (3) and (4) j − ≤ j ≤ n , b j = b + 2 j − (cid:88) k =1 u k L k + 3 j − (cid:88) k =1 v k L k and(35) a j = ν + j − (cid:88) k =1 b k L k + j − (cid:88) k =1 u k L k + j − (cid:88) k =1 v k L k . (36)Using (35) to substitute for b k in (36),(37) a j = ν + ( T j − − T ) b + 2 j − (cid:88) k =1 k − (cid:88) i =1 u i L i L k + 3 j − (cid:88) k =1 k − (cid:88) i =1 v i L i L k + j − (cid:88) k =1 u k L k + j − (cid:88) k =1 v k L k where 2 ≤ j ≤ n . Similarly, using the second part of (2) in place of (3), ν n = ν + ( T n − T ) b + 2 n (cid:88) k =1 k − (cid:88) i =1 u i L i L k + 3 n (cid:88) k =1 k − (cid:88) i =1 v i L i L k + n (cid:88) k =1 u k L k + n (cid:88) k =1 v k L k = ⇒ b = ν n − ν − (cid:80) nk =1 (cid:80) k − i =1 u i L i L k − (cid:80) nk =1 (cid:80) k − i =1 v i L i L k − (cid:80) nk =1 u k L k − (cid:80) nk =1 v k L k T n − T . (38)Then a j and b j are also determined for 2 ≤ j ≤ n , by substitution for b in (37) and in (35). For 1 ≤ j ≤ n define θ j : [0 , L j ] → R by θ j ( t ) := a j + b j t + c j t + d j t . Define θ : [ T , T n ] → R by θ ( T ) := ν and θ ( t ) := θ j ( t − T j − ) for t ∈ ( T j − , T j ]. Because the a j , b j , c j , d j satisfy (3), (4), as well as the auxiliary conditions (2), the following lemma is easilyverified. Lemma 6. θ : [ T , T n ] → R is a cubic polynomial spline with θ ( T ) = ν and θ ( T n ) = ν n . (cid:3) ECOND ORDER SPIRAL SPLINES 11
So for any u, v ∈ R n , y θ : [ T , T n ] → R is a C second order spiral spline satisfying y (1) θ ( T ) = V and y (1) θ ( T n ) = V n . For1 ≤ j ≤ n define z j : R n × R n → R by(39) z j ( u, v ) := (cid:90) L j (cos θ j ( t ) , sin θ j ( t )) dt, and define z : R n × R n → ( R ) n by z ( u, v ) := ( z ( u, v ) , z ( u, v ) , . . . , z n ( u, v )).Then y θ interpolates ( Y, V ) at T precisely when z j ( u, v ) = Y j − Y j − for all 1 ≤ j ≤ n , namely when ( u, v ) ∈ R n × R n satisfiesthe nonlinear system of 2 n scalar equations in 2 n scalar unknowns(40) z ( u, v ) = z ∗ := ( Y − Y , Y − Y , . . . , Y n − Y n − ) ∈ ( R ) n whose right hand sides are found from Y . The key to solving such a system is a satisfactory initial guess (ˆ u, ˆ v ) to a solution( u, v ). Set ˆ u = ˆ c and ˆ v = ˆ d , where the ˆ a j , ˆ b j , ˆ c j , ˆ d j are found in §
4. By Corollary 1, z (ˆ u, ˆ v ) = z ∗ + O ( (cid:15) ).In practice the integrals on the right hand side of (39) are not easy to express as functions of u, v . So for numericalcomputations we replace the z j by approximations using the composite Simpson’s Rule. Then z becomes an explicit functionof ( u, v ) and (40) is solved using Mathematica’s FindRoot with (ˆ u, ˆ v ) as an initial guess. Example 4.
Using composite Simpson with intervals to estimate the z j , FindRoot took around . seconds to obtain anumerical solution θ of (40) for the data in Example 1. Figure 4 shows the interpolant y θ (yellow), together with the initialestimate y ˆ θ and data from Figure 1. (cid:3) Example 5.
Perhaps it is no surprise that y θ is successful for the data of Example 2, because y ˆ θ was already nearlyinterpolating. Figure 5 shows y θ (yellow), together with y ˆ θ and the data. FindRoot took some 0.3 seconds to find y θ . (cid:3) Example 6. ex6 More impressively, it takes some 0.3 seconds to improve y ˆ θ in Example 3 to the interpolant y θ shown (yellow)in Figure 6. Whereas y ˆ θ comprehensively failed to interpolate, the C cubic spline ˆ θ : [ T , T n ] → R retains considerableinformative power. This power, masked to some degree by the second order spiral spline y ˆ θ : [ T , T n ] → R , is key to themethod of the present section § (cid:3) Conclusion
Given interpolation data (
Y, V ) for convex unit-speed C curves in R , relative to some sequence T , we develop a methodfor interpolating ( Y, V ) at T by a C second order spiral spline y θ : [ T , T n ] → R . The method is developed using asymptoticarguments to define a pair of tridiagonal linear systems whose solutions define a real-valued C cubic polynomial spline ˆ θ with prescribed values at T , T n . For sufficiently finely sampled data, ˆ θ already defines a C second order spiral spline y ˆ θ that approximately interpolates ( Y, V ) at T . More importantly, ˆ θ is used to start a numerical method that finds another C cubic polynomial spline θ . Then y θ is the desired C second order interpolating spiral spline. References [1] Di Antonio L. The fabric of the universe is most perfect: Euler’s research on elastic curves, in Euler at, 300: an appreciation, 239–260,Mathematical Association of America (2007).[2] Audoly, B., Pomeau, Yves, Elasticity and Geometry: From Hair Curls to the Non-linear Response of Shells, Oxford University Press (2010).[3] Baran, Ilya, Lehtinen, Jaako and Popovic, Jovan, Sketching Clothoid Splines Using Shortest Paths, Eurographics 0 (1981).[4] Florence Bertails-Descoubes, Super-Clothoids, Eurographics 31 (2012)[5] G. Birkhoff, H. Burchard and D. Thomas, Nonlinear Interpolation by Splines, Pseudosplines and Elastica, General Motors Research LaboratoriesReport 468, (1965), 1–13.[6] Bruckstein, Alfred M., Holta Robert, J., and Netravalia, Arun, N., Discrete elastica, Applicable Analysis 78 (2001) 453–485.[7] Brunnett, Guido and Wendt, J¨org, A Univariate Method for Plane Elastic Curves, Computer Aided Geometric Design 14 (1997) 273–292.[8] Coope, Ian, D., Curve Fitting With Nonlinear Spiral Splines, Department of Mathematics, University of Canterbury NZ, No 63 (1991), 1–14.[9] Dai, Ran and Cochran, John E. Jr. (2010) Path Planning and State Estimation for Unmanned Aerial Vehicles in Hostile Environments, Journalof Guidance, Control and Dynamics 33 (2), 595–601.[10] de Boor, Carl, A Practical Guide to Splines, Applied mathematical Sciences 27, Springer (1978).[11] Eberly, David, Moving Along a Curve with Specified Speed, (2007) 1–15.[12] Edwards, J.A., Exact equations of the nonlinear spline, ACM Transactions on Mathematical Software 18 (1992), 174–192.[13] Golomb, Michael and Jerome, Joseph, Equilibria of the Curvature Functional and Manifolds of Interpolating Spline Curves, SIAM J. onMathematical Analysis 13 (3) (1982), 412–458.[14] Alonzo Kelly and Brian Nagy, Reactive Nonholonomic Trajectory Generation via Parametric Optimal Control, International Journal of RoboticsResearch 22 (2003), 583–601.[15] Levien, Raph, The Elastica: a Mathematical History, Technical Report No. UCB/EECS-2008-103
August 23, 2008.[16] E.H. Lee and G.E. Forsythe, Variational Study of Nonlinear Spline Curves, SIAM Review 15 (1973), 120–133.[17] Linners, Anders, Unified Representations of Nonlinear Splines, Journal of Approximation Theory 84 (1996), 315–350.[18] Looker, Jason R., Constant Speed Interpolating Paths, DSTO Defence Science and Technology Organisation, AR 014-939. DSTO-TN-0989,March 2011.[19] D.S. Meek and R.S.D. Thomas, A Guided Clothoid Spline, Computer Aided Geometric Design 8 (1991) 163–174.[20] D.S. Meek and D.J. Walton, An Arc Spline Approximation to a Clothoid, Journal of Computational and Applied Mathematics 170 (2004),59–77.[21] Peterson, John W., Arc-Length Parameterization of Spline Curves Interpolation errors accumulate along the trajectory. The related curve ¯ y ˆ θ of § Y Y Y Y Y Y Y Y Y Y Y V V Figure 4. y θ (yellow) and y ˆ θ for Example 4 [22] Singer, David A., Lectures on Elastic Curves and Rods, in AIP Conference Proceedings 1002, Curvature and Variational Modelling in Physicsand Biophysics, Santiago de Compostela, Spain, 17–18 September 2007 (eds ´Oscar J. Garay, Eduardo Garcia-R´ıo, Ram´on V´azquez-Lorenzo),3–32.[23] Stoer, Josef, Curve Fitting With Clothoidal Splines, J. Research of the National Bureau of Standards 87 No. 4 (1982) 318–346.[24] Arc-Length Parameterized Spline Curves for Real-Time Simulation, Wang, Hongling, Keaney, Joseph, Atkinson, Kendall, Curve and SurfaceDesign: Saint Malo (2002) Lyche, T., Mazure, M.-L., Schumacher, L. (eds) 387–396. [email protected] (Department of Mathematics & Statistics, The University of Western Australia, 35 Stirling High-way, Crawley, WA 6009, AUSTRALIA)
ECOND ORDER SPIRAL SPLINES 13 Y Y Y Y Y Y Y Y Y Y Y V V - - - Figure 5. y θ (yellow) and y ˆ θ in Example 5 Y Y Y Y Y Y Y Y Y Y Y V V - - - Figure 6. y θ (yellow) and y ˆ θθ