- Email: [email protected]
A scheme for interpolation with trigonometric spline curves
Journal of Computational and Applied Mathematics 263 (2014) 246–261
Contents lists available at ScienceDirect
Journal of Computational and Applied Mathematics journal homepage: www.elsevier.com/locate/cam
A scheme for interpolation with trigonometric spline curves Imre Juhász a , Ágoston Róth b,∗ a
Department of Descriptive Geometry, University of Miskolc, H-3515 Miskolc-Egyetemváros, Hungary
b
Department of Mathematics and Computer Science, Babeş–Bolyai University, RO-400084 Cluj-Napoca, Romania
article
abstract
info
We present a method for the interpolation of a given sequence of data points with C n continuous trigonometric spline curves of order n+1 (n ≥ 1) that are produced by blending elliptical arcs. Ready to use explicit formulae for the control points of the interpolating arcs are also provided. Each interpolating arc depends on a global parameter α ∈ (0, π) that can be used for global shape modification. Associating non-negative weights with data points, rational trigonometric interpolating spline curves can be obtained, where weights can be used for local shape modification. The proposed interpolation scheme is a generalization of the Overhauser spline, and it includes a C n Bézier spline interpolation method as the limiting case α → 0. © 2013 Elsevier B.V. All rights reserved.
Article history: Received 27 May 2013 Received in revised form 12 December 2013 Keywords: Interpolation Trigonometric spline Blending Bézier spline Overhauser spline Shape parameters
1. Introduction Nowadays, in CAD and CAGD the standard description form of curves is c(t ) =
n
Fi (t )di ,
t ∈ [a, b] ⊂ R
(1)
i =0
where di ∈ Rδ (δ ≥ 2) are called control points, and Fi : [a, b] → R are sufficiently smooth combining functions. The most well-known examples are Bézier, B-spline and NURBS curves. In general, curve (1) does not pass through the control points just approximates the shape of the control polygon determined by them. However, if functions {Fi (t ) : t ∈ [a, b]}ni=0 are linearly n independent, we can always solve the following interpolation problem. Given the sequence of data points pi ∈ Rδ i=0 along with associated strictly monotone parameter values {ti ∈ [a, b]}ni=0 , find those control points of (1) for which equalities c (ti ) = pi ,
i = 0, 1, . . . , n
are fulfilled. This problem can be reduced to solve the system of linear equations F0 (t0 ) F0 (t1 )
∗
.. . F0 (tn )
F1 (t0 ) F1 (t1 )
.. . F1 (tn )
··· ··· .. . ···
d0 Fn (t0 ) Fn (t1 ) d1
.. .. . . Fn (tn ) dn
p0
p1 =. . . pn
Corresponding author. Tel.: +40 264405300. E-mail addresses: [email protected] (I. Juhász), [email protected] (Á. Róth).
0377-0427/$ – see front matter © 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.cam.2013.12.034
I. Juhász, Á. Róth / Journal of Computational and Applied Mathematics 263 (2014) 246–261
247
Fig. 1. The arc ai (t ) of the Overhauser spline is a blend of parabolas ci (t ) and ci+1 (t ).
for the unknown control points {di }ni=0 . A drawback of this method is that the interpolating curve will be globally controlled by data points, i.e., the displacement of any pi results in the change of the shape of the whole curve, even if the combining functions are splines. This disadvantage can be overcome by applying a local interpolating scheme, that was initiated by Overhauser [1]. He considered a sequence of data points pi with associated parameter values ti and constructed a C 1 cubic interpolating spline curve, the arcs of which are linearly blended parabolas. The ith (i = 1, 2, . . . , n − 2) arc of his spline is of the form ai (t ) =
1−
t − ti ti+1 − ti
ci (t ) +
t − ti ti+1 − ti
ci+1 (t ),
t ∈ [ti , ti+1 ],
where ci (t ), t ∈ [ti−1 , ti+1 ] is the parabolic arc that passes through the points pi−1 , pi , pi+1 at ti−1 , ti , ti+1 , respectively. This construction is illustrated in Fig. 1. The blending concept of Overhauser has various generalizations. Pobegailo [2] blended line segments and circular arcs with polynomials of degree 2n + 1 (specified in power basis) to produce interpolating curve of variable smoothness at data points. Wiltsche [3] studied blending in a wider context, used polynomials like [2] but specified them in Bernstein basis to blend different types of parametric curves, in addition introduced shape parameters. Wenz [4] linearly blended line segments and circular arcs to produce G1 interpolating spline. Szilvási-Nagy and Vendel [5] applied trigonometric blending of line segments and circular arcs to produce G2 interpolating spline. In order to describe an interpolation subspline scheme of degree 2k − 1 and of C k−1 continuity, Röschel [6] uses B-spline basis functions of degree k − 1 (k ≥ 1) for blending Lagrange interpolants of degree k that interpolate k consecutive data points at given strictly increasing parameter values. His method also includes procedures like Overhauser or quintic interpolation schemes as special cases. Replacing the Lagrange by Hermite interpolants, Gfrerrer and Röschel [7] generalize the previous idea. In this case, the generated polynomial subspline interpolates the Hermite input data that consists of parameter values and corresponding zeroth or higher order derivatives, and the common order k of B-spline basis functions used for blending is greater than or equal to the maximal order of derivatives appearing in the interpolation input. These blended Hermite interpolants will also be C k−1 continuous. Schneider [8] uses shape parameters (i.e., tension-like or straightness factors) and special weight functions in order to pull the parabolic arcs of the classical Overhauser interpolation scheme towards the chords between consecutive data points and to ensure the smoothness of the obtained interpolating curves at their common joints, respectively. The advantage of these methods is that there is no need for the solution of linear systems; the interpolating curve can be computed locally, as a consequence of which the spline curve is locally modifiable. A drawback of these procedures is that we lose the control point representation when the arcs to be blended are specified by control points. Our proposed method combines the assets of the two options described above, i.e., it produces the interpolating spline locally and provides directly the control points of the interpolating arcs with arbitrary order of continuity at joints. We demonstrate the method for trigonometric curves, that include polynomial curves as a special case. The advantages of the proposed method are: treats both trigonometric and polynomial local interpolating spline curves in a unified way; provides a ready to use control point based description which is advantageous from data storage and transmission point of view and is suitable for nowadays CAD/CAM systems; it is intuitive and user friendly; possesses both global and local (its rational counterpart) shape parameters; offers a more versatile tool for shape design than polynomial methods (e.g. ensures higher order circular precision), since these trigonometric curves form a proper subset of rational curves. The outline of the rest of the paper is as follows. Section 2 recalls the normalized B-basis of the vector space α F2n = span {cos(kt ), sin(kt ) : t ∈ [0, α ]}nk=0
(2)
of truncated Fourier series of finite order n ≥ 1 (i.e., trigonometric polynomials of degree at most 2n), presents some of its properties that will be used latter and specifies the curve defined by this basis in the form (1). A subdivision formula for the special case n = 1 is also provided for these curves. Section 3 proposes a trigonometric blending function and shows its essential properties. Section 4 constructs and subdivides a single elliptical arc that interpolates three consecutive data points at specified parameter values. Section 5 uses the blending function introduced in Section 3 to produce the arcs of a C n interpolating spline curve of order n + 1, by blending two quadratic trigonometric arcs, the global shape parameter of which is α . Ready to use formulae for the control points of the interpolating arcs are provided. In Section 6 it is shown that when α → 0 trigonometric arcs degenerate to Bézier curves of degree 2n + 1. By assigning non-negative weights of rank 1
248
I. Juhász, Á. Róth / Journal of Computational and Applied Mathematics 263 (2014) 246–261
with data points, Section 7 formulates the rational counterpart of our local interpolating scheme with rational trigonometric spline curves. Conclusions section closes the paper. 2. Trigonometric basis functions and curve We will produce interpolating trigonometric curves, that are described by means of the combination of control points and basis functions, therefore we need a proper basis in the space of trigonometric polynomials. We will use a transformed version of the basis specified in Sánchez-Reyes [9]. Let α ∈ (0, π) be an arbitrarily fixed parameter and consider the function system
Aα2n,i (t ) : t ∈ [0, α ] i=0 =
2n
α 2n−i c2n ,i sin
α−t
2
sini
t
2
2n : t ∈ [0, α ] ,
(3)
i =0
where the normalizing positive coefficients i
α c2n ,i
2 α i−2r n i−r 1 , 2 cos = r 2 sin2n α2 r =0 i − r
i = 0, 1, . . . , 2n
fulfill the symmetry α α c2n ,i = c2n,2n−i , i = 0, 1, . . . , n. Function system (3) is the B-basis of the vector space (2). Therefore, 2n
(4)
Aα2n,i (t ) ≡ 1
(5)
i =0
and Aα2n,i (0) = Aα2n,i (α) =
1, 0,
for i = 0, otherwise,
1, 0,
(6)
for i = 2n, otherwise.
(7)
Basis functions (3) enjoy the symmetry Aα2n,i (α − t ) = Aα2n,2n−i (t ) ,
∀t ∈ [0, α ] , i = 0, 1, . . . , n
(8)
since α
α
A2n,i (α − t ) = c2n,i sin
2n−i
α − (α − t )
sin
2
α i = c2n ,2n−i sin
= Aα2n,2n−i (t ),
α−t
2
sin2n−i
i
α−t
2
t
2
∀t ∈ [0, α ] . The product of basis functions A2n,i (t ) and Aα2m,j (t ) can be expressed in the form α
α α c2n ,i c2m,j α Aα2n,i (t )Aα2m,j (t ) = α A2(n+m),i+j (t ), c2(n+m),i+j
∀t ∈ [0, α ]
(9)
since α
α
α
A2n,i (t )A2m,j (t ) = c2n,i sin
2n−i
α−t
2
α α 2(n+m)−i−j = c2n ,i c2m,j sin
=
α α c2n ,i c2m,j
c2α(n+m),i+j
sin
i
t
2
α−t 2
Aα2(n+m),i+j (t ),
α
2m−j
· c2m,j sin
α−t 2
j
sin
t
2
t
sini+j
2
∀t ∈ [0, α ] .
Using the linear reparametrization t (v) = αv, v ∈ [0, 1], one obtains from (3) the system Aα2n,i (αv) : v ∈ [0, 1] i=0 =
2n
α i c2n ,i sin
αv 2
sin2n−i
α (1 − v) 2
2n : v ∈ [0, 1] , i=0
by means of which it is easy to show that (3) is a generalization of the Bernstein basis of degree 2n.
I. Juhász, Á. Róth / Journal of Computational and Applied Mathematics 263 (2014) 246–261
249
Proposition 2.1 (Bernstein Polynomials of Degree 2n as Special Case). If α → 0 basis functions (3) converge to the Bernstein polynomials of degree 2n, i.e., lim Aα2n,i (αv) = B2n i (v),
∀v ∈ [0, 1] , i = 0, 1, . . . , 2n.
α→0
(10)
Proof. In Ref. [9] functions (3) are obtained by raising the identity Aα2,0 (αv) + Aα2,1 (αv) + Aα2,2 (αv) = 1,
∀v ∈ [0, 1]
to the power n, i.e., 2n
Aα2n,i (αv) = Aα2,0 (αv) + Aα2,1 (αv) + Aα2,2 (αv)
n
= 1,
∀v ∈ [0, 1]
i =0
which in the limiting case α → 0 takes the form
1
sin2 sin2 α2
lim
α→0
α (1 − v)
+
2
2 cos α2
sin sin2 α2
α (1 − v)
2
sin
αv 2
+
1
sin2 sin2 α2
αv n 2
n = (1 − v)2 + 2 (1 − v) v + v 2 = ((1 − v) + v)2n 2n = B2n i (v) i=0
for all values of v ∈ [0, 1]. Corollary 2.1 (Product Rule of Bernstein Polynomials as Special Case). The limiting case α → 0 of the linear reparametrization t (v) = αv, v ∈ [0, 1] of the product rule (9) results in the well-known identity
2n k
2m B2n k (v)Bℓ (v) =
2m
ℓ
2(n+m) k+ℓ
(n+m) B2k+ℓ (v),
∀v ∈ [0, 1]
of Bernstein polynomials for all k = 0, 1, . . . , 2n and ℓ = 0, 1, . . . , 2m. Definition 2.1 (Trigonometric Curves of Order n). Let n ≥ 1 be a finite order. By means of functions (3) and control points
2n
di ∈ Rδ i=0 the curve of order n (degree 2n)
gαn (t ) =
2n
Aα2n,i (t )di ,
t ∈ [0, α ] ,
(11)
i=0
can be defined, where α is a global shape parameter. Curves (11) form a proper subset of rational curves, i.e., any curve of type (11) has a rational parametrization. However, rational parametrizations have drawbacks, e.g., their higher order derivatives become very involved and the parametrization of the curve is often poor (cf. [10]), thus it is worth using trigonometric parametrization. Basis functions (3) form the normalized B-basis of function space (2), therefore there must be a recursive corner cutting algorithm for curves (11), cf. [11]. Unfortunately, the general formula is unknown as yet, there is a scheme only for a very specific arrangement of control points in [12]. We managed to solve the problem only in the special case n = 1 for any placement of control points. We will need this result in Section 4. Let us consider the quadratic curve gα1 (t ) =
2
Aα2,i (t )di ,
t ∈ [0, α ] ,
i=0
and the functions µ10 (·; α) , µ11 (·; α) , µ20 (·; α) : [0, α ] → [0, 1], where
µ (t ; α) = 1 − 1 0
µ (t ; α) = 1 1
, sin α2 cos 2t sin2
t
t , sin α2 sin α2 − sin α− cos 2t 2
µ (t ; α) = 1 − 2 0
t sin α− 2
2
t sin α− cos 2t 2 sin α2
.
(12)
250
I. Juhász, Á. Róth / Journal of Computational and Applied Mathematics 263 (2014) 246–261
For any arbitrarily chosen value β ∈ (0, α) we introduce the subdivision points d10 = 1 − µ10 (β; α) d0 + µ10 (β; α) d1 ,
d11 = 1 − µ11 (β; α) d1 + µ11 (β; α) d2
and d20 = 1 − µ20 (β; α) d10 + µ20 (β; α) d11 ,
using which we define curves β
β
β
β
α−β A2,0 (t )d20
α−β A2,1 (t )d10
g1,ℓ (t ) = A2,0 (t ) d0 + A2,1 (t )d10 + A2,2 (t )d20 , α−β g1,r (t )
=
+
+
t ∈ [0, β ] ,
α−β A2,2 (t )d2 ,
t ∈ [0, α − β ] .
(13) (14)
By means of settings and notations above, one can easily prove the next proposition. Proposition 2.2 (Subdivision of Quadratic Trigonometric Curves). Quadratic trigonometric curves (13) and (14) form two consecutive arcs of curve (12) and the continuity at their common point d20 is C ∞ , i.e., β g1,ℓ (t ) = gα1 (t ), α−β g1,r (t )
∀t ∈ [0, β ] ,
α
= g1 (β + t ) ,
(15)
∀t ∈ [0, α − β ]
(16)
and dz dt
β
g (t ) z 1,ℓ
=
t =β
dz
α−β
g (t ) dt z 1,r
,
∀z ∈ N.
t =0
3. Blending function and the blended curve We are going to describe an interpolation scheme that is based on curve blending, therefore first of all we specify an appropriate blending function, then we show an essential continuity property of the resulted blended curve. 3.1. Trigonometric blending function Definition 3.1 (Trigonometric Blending Function of Degree 2n). By means of basis functions (3) we define the function bα2n (t ) =
2n
1 α A2n,n (t ) + Aα2n,i (t ), 2 i=n+1
t ∈ [0, α ]
(17)
that we will use for blending curves. Blending function (17) is non-negative and strictly increasing, specifically bα2n (0) = 0
(18)
bα2n (α) = 1
(19)
and due to Eqs. (6) and (7). Proposition 3.1 (Symmetry Property of the Blending Function). Blending function (17) is symmetric in the sense bα2n (t ) = 1 − bα2n (α − t ) ,
∀t ∈ [0, α ]
i.e., its graph is symmetric with respect to the point of coordinates α2 , 12 . Proof. Subsequent application of Eqs. (8) and (5) yields 1 − bα2n (α − t ) = 1 −
2n
1 α A2n,n (α − t ) − Aα2n,i (α − t ) 2 i=n+1 1
= 1 − Aα2n,n (t ) − 2
Aα2n,i (t )
i=0 2n
=
n −1
1 α A2n,n (t ) + Aα2n,i (t ) = bα2n (t ) 2 i=n+1
for all values of t ∈ [0, α ].
I. Juhász, Á. Róth / Journal of Computational and Applied Mathematics 263 (2014) 246–261
251
Proposition 3.2 (Vanishing Higher Order Derivatives). Derivatives of blending function (17) vanish at the endpoints up to the order n − 1, i.e., dz dt
bα (t ) z 2n
= 0,
(20)
=0
(21)
t =0
dz
b2n (t ) dt z α
t =α
for all z = 1, 2, . . . , n − 1. Proof. Function system (3) is an extended Chebyshev system, thus result of Mazure [13] on the derivatives of B-basis in extended Chebyshev spaces can be applied, that results in the Eqs. dz
A2n,i (t ) z dt α
dz
=
dt
t =0
A (t ) z 2n,2n−i α
= 0,
0 ≤ z < i ≤ 2n.
(22)
t =α
This immediately yields dz dt
bα (t ) z 2n
= 0,
z = 1, 2, . . . , n − 1.
(23)
t =0
Basis functions (3) form a partition of unity, thus 2n dz α A2n,i (t ) = 0, z dt i=0
∀t ∈ [0, α ] , z ≥ 1.
Taking into account Eqs. (22), (23), we obtain
n −1 dz α A ( t ) 2n,i dt z i=0
+
t =0
1 dz 2 dt
A (t ) z 2n,n α
= 0,
z = 1, 2, . . . , n − 1.
(24)
t =0
Utilization of identities (8) and their derivatives dz dt
Aα (t ) z 2n,i
=
t =α
dz dt
Aα (α − t ) z 2n,2n−i
= (−1)z
t =α
dz dt
Aα (t ) z 2n,2n−i
t =0
α
for i = 0, 1, . . . , n yields the derivative of b2n (t ) at t = α in the form
A ( t ) 2n,i dt z i=n+1 2n dz
α
+
t =α
1 dz 2 dt
A ( t ) z 2n,n α
n −1 dz α 1 dz α = A2n,n (α − t ) A − t + (α ) 2n,i z z dt i=0 2 dt t =α t =α t =α
= (−1)
z
n dz α A2n,i (t ) z dt i=0
t =0
A2n,n (t ) + =0 z 2 dt t =0 1 dz
α
due to Eq. (24). 3.2. Curve blending Let us consider curves c1 (t ) and c2 (t ), t ∈ [0, α ] that fulfill conditions c1 (0) = c2 (0) = p1 ,
(25)
c1 (α) = c2 (α) = p2 . By means of the blending function (17) we generate the curve c(t ) = 1 − bα2n (t ) c1 (t ) + bα2n (t )c2 (t ),
t ∈ [0, α ] .
(26)
It is obvious, that c(t ) is a convex combination of curves c1 (t ) and c2 (t ), and the blended curve fulfills the endpoint interpolation properties
c(0) = p1 , c (α) = p2 .
252
I. Juhász, Á. Róth / Journal of Computational and Applied Mathematics 263 (2014) 246–261
Proposition 3.3 (Higher Order Hermite Conditions of the Blended Curve). If curves c1 (t ) and c2 (t ), t ∈ [0, α ] are at least n times continuously differentiable at t = 0 and t = α , then for the blended curve (26) equalities dz
dz c (t ) = z c1 (t ) , dt z dt t =0 t =0 dz dz c ( t ) c ( t ) = 2 z z dt
dt
t =α
t =α
hold for all z = 1, 2, . . . , n. Proof. The zth order derivative of c(t ) is dz dt
c(t ) = z
z k z d k=0
dt k
k
d z −k
1 − bα2n (t )
dt
c (t ) + z −k 1
dk α d z −k b ( t ) c2 (t ) . 2n dt k dt z −k
At first we assume that 0 ≤ z < n and evaluate this derivative at t = 0. In this case the only non-vanishing term of the above sum is
dz 1 − b2n (0) c (t ) z 1
α
α
dt
+ b2n (0)
t =0
dz dt
c (t ) z 2
t =0
(i.e., the one that belongs to k = 0) due to Proposition 3.2. Thus, we obtain dz dt
c (t ) z
dz dz α = 1 − b2n (0) + b2n (0) z c2 (t ) c 1 (t ) z dt dt t =0 t =0 dz = c (t ) z 1 α
t =0
dt
t =0
because of Eq. (18). Now, we consider the case z = n at t = 0 that yields dn
dn
c (t ) n dt
t =0
= + c1 (t ) n dt t =0 dn = − c1 (t ) dt n t =0 dn c (t ) = n 1 dt
n dn α c1 (0) + d bα (t ) c2 (0) 1 − b ( t ) 2n 2n n n dt dt t =0 t =0
n dn α p1 + d bα (t ) p1 b ( t ) 2n n 2n dt n dt t =0 t =0
t =0
using assumption (25). The equality for the derivative at the endpoint (t = α ) can be proved analogously. 4. An interpolating elliptical arc We are going to use elliptical arcs for local interpolation, therefore we have to describe the elliptical arc as a quadratic trigonometric curve of type (11) that interpolates three consecutive points at specified parameter values. The task is as follows. Given the sequence of data points p0 , p1 , p2 and associated parameter values u0 < u1 < u2 , with the constraint u2 − u0 < π , find those control points q0 , q1 , q2 for which the curve e(u) =
2
α +α1
A2,0j
(u)qj ,
u ∈ [0, α0 + α1 ]
(27)
j=0
passes through the given points at the given parameter values, i.e., equalities e(0) = p0 , e (α0 ) = p1 , e (α0 + α1 ) = p2 are fulfilled, where αi = ui+1 − ui , (i = 0, 1). The solutions are q0 = p0 , α +α1
q1 =: q1 0 q2 = p2 .
=
1 α +α1
A2,01
(α0 )
α +α1
p1 − A2,00
α +α1
(α0 ) p0 − A2,02
(α0 ) p2 ,
I. Juhász, Á. Róth / Journal of Computational and Applied Mathematics 263 (2014) 246–261
253
Fig. 2. The elliptical arc that interpolates data points p0 , p1 , p2 that is split into two arcs. Settings are α = π/2 and knots are chosen by centripetal parametrization.
Using the subdivision formula detailed in Proposition 2.2, we split the elliptical arc (27) at the parameter value u = α0 . Control points of the left and right arcs are as follows. The left arc eℓ (u) =
2
α
A2,0j (u) aj,ℓ ,
u ∈ [0, α0 ]
j =0
is defined by control points a0,ℓ = p0 , α
0 a1,ℓ =: a1,ℓ = 1 − µ10 (α0 ; α0 + α1 ) p0 + µ10 (α0 ; α0 + α1 ) q1 , a2,ℓ = p1 ,
while the right arc er (u) =
2
α
A2,1j (u)aj,r ,
u ∈ [0, α1 ]
j=0
is generated by a0,r = p1 , α
a1,r =: a1,1r = 1 − µ11 (α0 ; α0 + α1 ) q1 + µ11 (α0 ; α0 + α1 ) p2 , a2,r = p2 .
The interpolating elliptical arc and its subdivision are illustrated in Fig. 2. Such left and right elliptical arcs will be used in the forthcoming section to produce C n (n ≥ 1) interpolating trigonometric spline curves. Remark 4.1. If points p0 , p1 , p2 are collinear the process described above results in two straight line segments parametrized in the form (11). 5. Interpolating trigonometric spline curve Our objective is to generalize the Overhauser spline curve by solving the following interpolation problem. Given the sequence of data points {pi }m i=0 with associated parameter values t0 < t1 < · · · tm−1 < tm , moreover the shape parameter α ∈ (0, π). Find the control points of a spline curve of continuity C n (n ≥ 1) that interpolates the given data points, and that consists of arcs of type (11). We intend to use elliptical arcs (27) defined by consecutive triplets of data points pi−1 , pi , pi+1 (i = 1, 2, . . . , m − 1), where the maximum of domains ti+1 − ti−1 is equal to α , therefore we have to scale parameter values ti that results in the new knots uα0 = t0 , uαi = uαi−1 +
α dmax
(ti − ti−1 ) ,
i = 1, 2, . . . , m,
254
I. Juhász, Á. Róth / Journal of Computational and Applied Mathematics 263 (2014) 246–261
where dmax = max {ti+1 − ti−1 : i = 1, 2, . . . , m − 1}. Thus the domain of the interpolating elliptical arc will be αi−1 + αi , α +α with αi = uαi+1 − uαi , (i = 0, 1, . . . , m − 1). Its control points are pi−1 , qi i−1 i , pi+1 , where α
qi i−1
+αi
=
1 α 1 +αi A2,i− 1
α
1 pi − A2,i− 0
(αi−1 )
+αi
α
(αi−1 ) pi−1 − A2,i−2 1
+αi
(αi−1 ) pi+1 .
We split this arc into two (also quadratic) elliptical arcs, as it is described in Section 4, that results in the left and right arcs α
α
α
α
α
i−1 i−1 1 + A2,i−2 1 (u)pi , ei,ℓi−1 (u) = A2,i− 0 (u)pi−1 + A2,1 (u)ai,ℓ
u ∈ [0, αi−1 ]
and α
α
α
α
α
ei,ir (u) = A2,i 0 (u) pi + A2,i 1 (u)ai,ir + A2,i 2 (u)pi+1 , where α
u ∈ [0, αi ] , α
ai,ℓi−1 = 1 − µ10 (αi−1 ; αi−1 + αi ) pi−1 + µ10 (αi−1 ; αi−1 + αi ) qi i−1
α ai,ir
αi−1 +αi
= 1 − µ (αi−1 ; αi−1 + αi ) qi 1 1
Blending
α ei,ir (u)
αi gn+ 1,i (u)
and
α ei+i 1,ℓ (u) (i αi
= 1 − b2n (u)
+αi
,
+ µ (αi−1 ; αi−1 + αi ) pi+1 . 1 1
= 1, 2, . . . , m − 1) we obtain the sequence of arcs
α ei,ir (u)
α
α
+ b2ni (u)ei+i 1,ℓ (u),
u ∈ [0, αi ] , i = 1, 2, . . . , m − 1.
(28) α
α
i+1 i Proposition 5.1 (Continuity of the Interpolating Trigonometric Spline Curve of Order n+1). Arcs gn+ 1,i (u) and gn+1,i+1 (u)(i = 1, 2, . . . , m − 2) are C n continuous at their common joint pi+1 .
Proof. On account of Proposition 3.3 dz
αi gn+ 1,i (u) z du dz du
= u=αi
α ei+i 1,ℓ (u) z du
,
u=αi
α
dz
g i+1 (u) z n+1,i+1
=
u =0
dz du
α
e i+1 (u) z i+1,r
u =0
that are equal for all z ≥ 0 due to Proposition 2.2. α
m−1
i In what follows, we determine the control points of blended arcs gn+ 1,i (u) : u ∈ [0, αi ] i=1 . We have that
α
α
α
α
α
i i i i i gn+ 1,i (u) = 1 − b2n (u) ei,r (u) + b2n (u)ei+1,ℓ (u)
=
2n α 1 α α α α α 1 − A2ni ,n (u) − A2ni ,j (u) A2,i 0 (u)pi + A2,i 1 (u)ai,ir + A2,i 2 (u) pi+1 2 j =n +1
+ =
2n α 1 αi α α α α A2n,n (u) + A2ni ,j (u) A2,i 0 (u)pi + A2,i 1 (u)ai+i 1,ℓ + A2,i 2 (u)pi+1 2 j =n +1
n −1
α 1 α α α α A2n,j (u) + A2ni ,n (u) A2,i 0 (u)pi + A2,i 1 (u) ai,ir + A2,i 2 (u)pi+1 2 j =0
+
αi
2n α 1 αi α α α α A2n,n (u) + A2ni ,j (u) A2,i 0 (u)pi + A2,i 1 (u)ai+i 1,ℓ + A2,i 2 (u)pi+1 . 2 j =n +1
α
By rearranging the right hand side, using identity (9) and collecting the coefficients of basis functions A2(i n+1),j (u) : u ∈ [0, αi ]
2(n+1)
, we obtain control points
j=0
αi αi α α α α c2ni ,j−1 c2,i1 α c2ni ,j−2 c2,i2 c2n,j c2,0 i p + a + pi+1 , for j = 0, 1, . . . , n, i α α α i , r c2(in+1),j c2(in+1),j c2(in+1),j αi α α α α α c2n,n+1 c2,i0 c2ni ,n c2,i1 1 αi c2ni ,n−1 c2,i2 αi α p + a + a + pi+1 , for j = n + 1, dj i = i α α α i,r i+1,ℓ c2(in+1),n+1 c2(in+1),n+1 2 c2(in+1),n+1 αi αi αi αi αi αi c2n,j c2,0 pi + c2n,j−1 c2,1 aαi + c2n,j−2 c2,2 pi+1 , for j = n + 2, n + 3, . . . , 2 (n + 1) , α
c2(in+1),j
α
c2(in+1),j
i+1,ℓ
α
c2(in+1),j
I. Juhász, Á. Róth / Journal of Computational and Applied Mathematics 263 (2014) 246–261
α
α
α
255
α
Fig. 3. The arc g3,i i (u) of the C 2 interpolating spline is obtained by blending elliptical arcs ei,ir (u) and ei+i 1,ℓ (u) with b4 i (u); α = 3 and knots ti are chosen by chord length parametrization.
Fig. 4. The effect of the order of continuity on the shape of an interpolating plane trigonometric spline curve; α = π2 and the initial knots were chosen by centripetal parametrization. Subfigures (a), (b) and (c) show the curve and its signed curvature plot for continuities C 1 , C 2 and C 3 , respectively.
α
where the convention c2ni ,j = 0 if j < 0 or j > 2n is adopted. Fig. 3 shows the construction of an arc of a C 2 interpolating spline along with its control points. Fig. 4 shows interpolating plane trigonometric spline curves of different orders of continuity for a fixed sequence of data points and centripetal parametrization along with their curvature plot. The above described method does not provide the first and last arcs of the interpolating spline, i.e., arcs that connect data points p0 , p1 and pm−1 , pm , respectively. In order to specify these arcs, we can apply, e.g., pseudo data points p−1 , pm+1 , or α α0 we can use arcs e1,ℓ (u), emm−−11,r (u). Both methods guarantee the required order of continuity at p1 and pm−1 , respectively.
256
I. Juhász, Á. Róth / Journal of Computational and Applied Mathematics 263 (2014) 246–261
Fig. 5. The effect of the parametrization of data points (the choice of knots ti ) both on a plane C 3 interpolating spline curve and its signed curvature plot; α = π/2. Observe that the uniform parametrization furnishes unwanted undulations.
Data points have a local influence on the interpolating curve, since the repositioning of data point pi affects at most four
αj
arcs, namely gn+1,j (u) : u ∈ 0, αj
i+1
.
j=i−2
α
α
i i Remark 5.1 (Placement of Control Points). Due to the construction of the arc gn+ 1,i (u), control points dj
are in the plane spanned by data points pi−1 , pi , pi+1 and pi , pi+1 , pi+2 , respectively.
n j =0
α
and dj i
2(n+1) j =n +2
Fig. 5 illustrates the effect of the parametrization (the choice of knots ti ) of data points, while Fig. 6 shows the influence of the global shape parameter α on the shape of the interpolating plane spline curve. Remark 5.2 (Higher Order Circular Precision). Due to the construction and local controllability of the arcs (28), our method is also able to reproduce arcs with circular precision of order n if at least 4 consecutive data points are uniformly distributed on a circle and α is the double of the central angle determined by two consecutive data points on the circle. Fig. 7 illustrates a C 2 interpolating trigonometric spline curve that has four arcs with second order circular precision. This circular precision of order n is the circular counterpart of the linear precision of Bézier curves, i.e., equidistantly placed Bézier points on a line generate a linear polynomial. Circular precision is important in engineering applications as mentioned in articles [14–16] that produce G1 or G2 continuous polynomial interpolating arcs with zeroth order circular precision. 6. Asymptotic property of the interpolating trigonometric spline curve This section examines the limiting case α → 0 of the interpolating C n continuous trigonometric spline curve that consists α
m−1
i of arcs gn+ 1,i (u) : u ∈ [0, αi ] i=1 of order n + 1. Observe that
α
lim αi = lim uαi+1 − uαi = lim (ti − ti−1 ) = 0 α→0 α→0 α→0 dmax
I. Juhász, Á. Róth / Journal of Computational and Applied Mathematics 263 (2014) 246–261
257
Fig. 6. The impact of the global shape parameter α on the shape of a C 2 interpolating spline curve; knots ti are specified by chord length parametrization.
Fig. 7. Data points {pi }8i=2 are uniformly distributed on a circle with central angle π6 of arcs. Parameter settings are n = 2 and α = π3 . Observe that the
6
zeroth, first and second order derivatives of the plane C 2 interpolating trigonometric arcs gα3,i (u) : u ∈ [0, α ] i=3 coincide with the same order derivatives of the support circle. As expected, the signed curvature plot is constant on the interval [u3 , u7 ].
for all i = 0, 1, . . . , m − 1. Using parametrization u(t ) = (αi−1 + αi )
t − ti−1 ti+1 − ti−1
,
t ∈ [ti−1 , ti+1 ]
and Proposition 2.1, we obtain for α → 0, that the local interpolating elliptical arc α
ei i−1
+αi
(u) =
2
α
qj,ii−1
+αi αi−1 +αi A2,0 (u),
u ∈ [0, αi−1 + αi ]
j =0
of order one degenerates to the quadratic Bézier curve (i.e., to the parabolic arc) e0i (t ) :=
2 j =0
where q00,i = pi−1 ,
q0j,i B2j
t − ti−1 ti+1 − ti−1
,
t ∈ [ti−1 , ti+1 ] ,
(29)
258
I. Juhász, Á. Róth / Journal of Computational and Applied Mathematics 263 (2014) 246–261
1
q01,i =
B21
ti −ti−1 ti+1 −ti−1
ti − ti−1 ti − ti−1 pi − pi−1 B20 − pi+1 B22 , ti+1 − ti−1 ti+1 − ti−1
q02,i = pi+1 . Due to equalities lim µ
α→0
1 0
(αi−1 + αi )
t − ti−1 ti+1 − ti−1
; αi−1 + αi
= lim µ α→0
1 1
(αi−1 + αi )
= lim µ20 (αi−1 + αi ) α→0
=
t − ti−1 ti+1 − ti−1
t − ti−1 ti+1 − ti−1 t − ti−1 ti+1 − ti−1
∈ [0, 1] ,
; αi−1 + αi ; αi−1 + αi
∀t ∈ [ti−1 , ti+1 ] ,
the special subdivision scheme detailed in Proposition 2.2 coincides with the classical de Casteljau algorithm of quadratic Bézier curves when α → 0. Splitting the parabola (29) into two (also quadratic) Bézier arcs, one obtains the left and right arcs e0i,ℓ (t ) = B20 (vi−1 (t )) pi−1 + B21 (vi−1 (t )) a0i,ℓ + B22 (vi−1 (t )) pi ,
t ∈ [ti−1 , ti ]
and e0i,r (t ) = B20 (vi (t )) pi + B21 (vi (t )) a0i,r + B22 (vi (t )) pi+1 ,
t ∈ [ti , ti+1 ] ,
where a0i,ℓ = a0i,r =
1−
1−
ti − ti−1
pi−1 +
ti+1 − ti−1 ti − ti−1 ti+1 − ti−1
q01,i +
ti − ti−1 ti+1 − ti−1
q01,i ,
ti − ti−1 pi+1 , ti+1 − ti−1
and
vi ( t ) =
t − ti ti+1 − ti
,
t ∈ [ti , ti+1 ] , i = 1, 2, . . . , m − 1.
Due to Corollary 2.1, by blending parabolas e0i,r (t ) and e0i+1,ℓ (t ) (i = 1, 2, . . . , m − 1, t ∈ [ti , ti+1 ]) one obtains the sequence of Bézier arcs g0n+1,i (t ) = 1 − b02n (t ) e0i,r (t ) + b02n (t )e0i+1,ℓ (t ),
t ∈ [ti , ti+1 ] , i = 1, 2, . . . , m − 1
(30)
of degree 2 (n + 1), where α
b02n (t ) = lim b2ni (t ) α→0
= lim
α→0
=
1 2
2n
α 1 αi A2n,n (αi vi (t )) + A2ni ,k (αi vi (t )) 2 k=n+1
B2n n (vi (t )) +
2n
B2n k (vi (t )) ,
∀t ∈ [ti , ti+1 ]
k=n+1
denotes the blending function of degree 2n. Naturally, Bézier arcs g0n+1,i (t ) and g0n+1,i+1 (t ) (i = 1, 2, . . . , m − 2) are still C n continuous at their joint pi . Observe that the Bézier ordinates
y0 = 0, y1 = 0, . . . , yn−1 = 0, yn =
of the blending function b02n (t ) = of the Bézier function
2n
k=0
1 2
, yn+1 = 1, . . . , y2n = 1
yk B2n k (vi (t )) of apparently even degree 2n, can be obtained by elevating the degree
2n−1
ℓ=0
xℓ Bℓ2n−1 (vi (t )) ,
t ∈ [ti , ti+1 ]
of odd degree 2n − 1 having Bézier ordinates
(x0 = 0, x1 = 0, . . . , xn−1 = 0, xn = 1, xn+1 = 1, . . . , x2n−1 = 1) .
I. Juhász, Á. Róth / Journal of Computational and Applied Mathematics 263 (2014) 246–261
259
Indeed, applying the geometric interpretation of the degree elevation, one has that y0 = x0 = 0,
0, j j 1 yj = x j −1 + 1 − xj = , 2n 2n 2 1, y2n = x2n−1 = 1,
j = 1, 2, . . . , n − 1,
j = n, j = n + 1, n + 2, . . . , 2n − 1,
i.e., the next proposition is proved. Proposition 6.1 (Degree Reduction). The actual degree of the polynomial blending function b02n (t ), t ∈ [ti , ti+1 ] is 2n − 1. The special case n = 1 of the above polynomial is the linear blending function vi (t ) that appears in the classical cubic Overhauser interpolation scheme, therefore our proposed method is a generalization of the classical cubic Overhauser
m−1
interpolation scheme for C n continuity. Naturally, this also means that Bézier arcs g0n+1,i (t ) : t ∈ [ti , ti+1 ] i=1 are of degree 2n + 1 and the control points of the ith arc can be expressed as
2n−1 2 2n−1 2 2n−1 2 j−1 1 j − 2 j 0 2 0 pi + ai,r + pi+1 , for j = 0, 1, . . . , n, 2n+1 2n+1 2n+1 j j j d0j = 2n−1 2 2n−1 2 2n−1 2 j 0 j−1 1 j − 2 2 0 pi + ai+1,ℓ + pi+1 , for j = n + 1, n + 2, . . . , 2n + 1. 2n+1 2n+1 2n+1 j
j
j
ti m i =0
Moreover, in this special case, original knots { } the shape parameter α .
do not have to be rescaled, since the Bernstein polynomials do not possess
7. Interpolating rational trigonometric spline curve m We can assign non-negative weights {wi }m i=0 of rank 1 to data points {pi }i=0 that will serve as additional local shape parameters of the interpolating spline curve. This weighted interpolating spline curve can be described by means of the rational counterpart of basis functions (3), i.e., by ωi Aα2n,i (t ) Rα2n,i (t ) = , t ∈ [0, α ] , i = 0, 1, . . . , 2n, 2n
j =0
ωj Aα2n,j (t )
2n
where ωj j=0 are non-negative weights of rank 1 that are unknown at the moment. A rational curve in Rδ can be considered as the central projection of a curve in the δ + 1 dimensional space from the origin on to the δ dimensional hyperplane w = 1 (assuming that the last coordinate of space Rδ+1 is denoted by w ). The space Rδ+1 is called the pre-image space. It facilitates our calculations if we do them in the pre-image space, where coordinates of data points are
wi pi , wi
i = 0, 1, . . . , m.
In the pre-image space calculations go like in Section 5. Control points of the interpolating spline in the pre-image space are
α
dj i
αi
ωj
=
α
α
α
α
α
α
α
c2ni ,j−1 c2,i1 α c2ni ,j−2 c2,i2 i w p + a + wi+1 pi+1 i i α α α i , r i i c c2(n+1),j c2(in+1),j 2(n+1),j , αi αi αi αi αi αi c2n c c c c c ,j 2,0 2n,j−1 2,1 2n,j−2 2,2 αi wi + αi wi,r + αi wi+1 αi c2(n+1),j c2(n+1),j c2(n+1),j
c2ni ,j c2,i0
α
c2ni ,n+1 c2,i0
α
wi pi +
c αi 2(n+1),n+1 αi αi c2n ,n+1 c2,0
α
α
α
for j = 0, 1, . . . , n,
c2ni ,n−1 c2,i2 c2ni ,n c2,i1 1 αi α ai+1,ℓ + ai,ir + αi wi+1 pi+1 αi c2(n+1),n+1 2 c2(n+1),n+1 , α α α α c2ni ,n c2,i1 1 αi c2ni ,n−1 c2,i2 αi
for j = n + 1, wi + αi wi+1,ℓ + wi,r + αi wi+1 α c2(in+1),n+1 c2(n+1),n+1 2 c2(n+1),n+1 αi αi α α α α c2n,j c2,0 c2ni ,j−1 c2,i1 α c2ni ,j−2 c2,i2 i w p + a + w p i i i + 1 i + 1 α α i+1,ℓ c αi c2(in+1),j c2(in+1),j 2(n+1),j , for j = n + 2, n + 3, . . . , 2n + 2. α α α α α α c2ni ,j c2,i0 c2ni ,j−1 c2,i1 α c2ni ,j−2 c2,i2 i w i + αi wi+1,ℓ + αi wi+1 αi c2(n+1),j
c2(n+1),j
c2(n+1),j
260
I. Juhász, Á. Róth / Journal of Computational and Applied Mathematics 263 (2014) 246–261
Fig. 8. Local shape modification of a C 2 rational interpolating trigonometric curve, the weight of the data point marked with filled disc takes the values 0.3, 1, 2 while the rest of the weights is 1; α = π/2 and knots ti are specified by centripetal parametrization.
In the model space control points and associated weights are
1 αi
ωj
α
m
dj i j=0
α
and ωj i
m j=0
, respectively. In a design situation
the user has to specify just data points pi and associated non-negative weights wi , control points and weights of the rational interpolating arcs are calculated automatically. The user can adjust the curvature of the interpolating spline curve at a data point locally, by altering the corresponding weight. This kind of shape modification is shown in Fig. 8. 8. Conclusions A method is proposed that produces interpolating trigonometric C n continuous (n ≥ 1) spline curves to a given sequence of data points. Arcs of the spline curve are computed locally, i.e., there is no need for the solution of a system of linear equations, and explicit formulae are provided for the control points of the arcs. Arcs are determined by blending two quadratic trigonometric curves (two elliptical arcs) with an appropriate trigonometric function. The domain of each arc depends on α , where α ∈ (0, π) is a global shape parameter of the interpolating spline curve. It is also shown that associating non-negative weights with data points rational trigonometric interpolating splines can be obtained, that are locally modifiable. The method is derived for trigonometric curves, that also include polynomial curves as the special case α → 0. The proposed interpolation scheme is a genuine generalization of the Overhauser spline. The generalizations are: the order of continuity n ≥ 1, the more versatile shape and the additional global shape parameter. A further direction of the research is to extend this method to smooth surface interpolation. Acknowledgments Imre Juhász carried out his research in the framework of the Center of Excellence of Mechatronics and Logistics at the University of Miskolc. The research of Ágoston Róth was supported by the European Union and the State of Hungary, cofinanced by the European Social Fund in the framework of TÁMOP-4.2.4.A/2-11/1-2012-0001 ‘National Excellence Program’. All assets, infrastructure and personnel support of his research was contributed by the Romanian national grant CNCSUEFISCDI/PN-II-RU-TE-2011-3-0047. The authors also thank the reviewers for their useful comments. References [1] A.W. Overhauser, Analytic Definition of Curves and Surfaces by Parabolic Blending, Technical Report No. SL68-40, Ford Motor Company Scientific Laboratory, 1968. [2] A.P. Pobegailo, Local interpolation with weight functions for variable smoothness curve design, Comput.-Aided Des. 23 (8) (1991) 579–582. [3] A. Wiltsche, Blending curves, J. Geom. Graph. 9 (1) (2005) 67–75. [4] H.-J. Wenz, Interpolation of curve data by blended generalized circles, Comput. Aided Geom. Des. 13 (8) (1996) 673–680. [5] M. Szilvási-Nagy, T.P. Vendel, Generating curves and swept surfaces by blended circles, Comput. Aided Geom. Des. 17 (2) (2000) 197–206.
I. Juhász, Á. Róth / Journal of Computational and Applied Mathematics 263 (2014) 246–261
261
[6] O. Röschel, An interpolation subspline scheme related to B-spline techniques, in: Computer Graphics International, 1997. Proceedings, IEEE, 1997, pp. 131–136. [7] A. Gfrerrer, O. Röschel, Blended Hermite interpolants, Comput. Aided Geom. Des. 18 (9) (2001) 865–873. [8] W. Schneider, A simple technique for adding tension to parabolic blending interpolation, Comput. Math. Appl. 12 (11) (1986) 1155–1160. [9] J. Sánchez-Reyes, Harmonic rational Bézier curves, p-Bézier curves and trigonometric polynomials, Comput. Aided Geom. Des. 15 (9) (1998) 909–923. [10] L. Piegl, On NURBS: a survay, Comput. Graph. Appl. 11 (1) (1991) 55–71. [11] E. Mainar, J.M. Peña, Corner cutting algorithms associated with optimal shape preserving representations, Comput. Aided Geom. Des. 16 (9) (1999) 883–906. [12] J. Sánchez-Reyes, Single-valued curves in polar coordinates, Comput.-Aided Des. 22 (1) (1990) 19–26. [13] M.-L. Mazure, Chebyshev–Bernstein bases, Comput. Aided Geom. Des. 16 (7) (1999) 649–669. [14] G. Farin, Geometric Hermite interpolation with circular precision, Comput.-Aided Des. 40 (4) (2008) 476–479. [15] D.J. Walton, D.S. Meek, G2 Hermite interpolation with circular precision, Comput.-Aided Des. 42 (9) (2010) 749–758. [16] Y. Li, C. Deng, C-shaped G2 Hermite interpolation with circular precision based on cubic PH curve interpolation, Comput.-Aided Des. 44 (11) (2012) 1056–1061.
Contents lists available at ScienceDirect
Journal of Computational and Applied Mathematics journal homepage: www.elsevier.com/locate/cam
A scheme for interpolation with trigonometric spline curves Imre Juhász a , Ágoston Róth b,∗ a
Department of Descriptive Geometry, University of Miskolc, H-3515 Miskolc-Egyetemváros, Hungary
b
Department of Mathematics and Computer Science, Babeş–Bolyai University, RO-400084 Cluj-Napoca, Romania
article
abstract
info
We present a method for the interpolation of a given sequence of data points with C n continuous trigonometric spline curves of order n+1 (n ≥ 1) that are produced by blending elliptical arcs. Ready to use explicit formulae for the control points of the interpolating arcs are also provided. Each interpolating arc depends on a global parameter α ∈ (0, π) that can be used for global shape modification. Associating non-negative weights with data points, rational trigonometric interpolating spline curves can be obtained, where weights can be used for local shape modification. The proposed interpolation scheme is a generalization of the Overhauser spline, and it includes a C n Bézier spline interpolation method as the limiting case α → 0. © 2013 Elsevier B.V. All rights reserved.
Article history: Received 27 May 2013 Received in revised form 12 December 2013 Keywords: Interpolation Trigonometric spline Blending Bézier spline Overhauser spline Shape parameters
1. Introduction Nowadays, in CAD and CAGD the standard description form of curves is c(t ) =
n
Fi (t )di ,
t ∈ [a, b] ⊂ R
(1)
i =0
where di ∈ Rδ (δ ≥ 2) are called control points, and Fi : [a, b] → R are sufficiently smooth combining functions. The most well-known examples are Bézier, B-spline and NURBS curves. In general, curve (1) does not pass through the control points just approximates the shape of the control polygon determined by them. However, if functions {Fi (t ) : t ∈ [a, b]}ni=0 are linearly n independent, we can always solve the following interpolation problem. Given the sequence of data points pi ∈ Rδ i=0 along with associated strictly monotone parameter values {ti ∈ [a, b]}ni=0 , find those control points of (1) for which equalities c (ti ) = pi ,
i = 0, 1, . . . , n
are fulfilled. This problem can be reduced to solve the system of linear equations F0 (t0 ) F0 (t1 )
∗
.. . F0 (tn )
F1 (t0 ) F1 (t1 )
.. . F1 (tn )
··· ··· .. . ···
d0 Fn (t0 ) Fn (t1 ) d1
.. .. . . Fn (tn ) dn
p0
p1 =. . . pn
Corresponding author. Tel.: +40 264405300. E-mail addresses: [email protected] (I. Juhász), [email protected] (Á. Róth).
0377-0427/$ – see front matter © 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.cam.2013.12.034
I. Juhász, Á. Róth / Journal of Computational and Applied Mathematics 263 (2014) 246–261
247
Fig. 1. The arc ai (t ) of the Overhauser spline is a blend of parabolas ci (t ) and ci+1 (t ).
for the unknown control points {di }ni=0 . A drawback of this method is that the interpolating curve will be globally controlled by data points, i.e., the displacement of any pi results in the change of the shape of the whole curve, even if the combining functions are splines. This disadvantage can be overcome by applying a local interpolating scheme, that was initiated by Overhauser [1]. He considered a sequence of data points pi with associated parameter values ti and constructed a C 1 cubic interpolating spline curve, the arcs of which are linearly blended parabolas. The ith (i = 1, 2, . . . , n − 2) arc of his spline is of the form ai (t ) =
1−
t − ti ti+1 − ti
ci (t ) +
t − ti ti+1 − ti
ci+1 (t ),
t ∈ [ti , ti+1 ],
where ci (t ), t ∈ [ti−1 , ti+1 ] is the parabolic arc that passes through the points pi−1 , pi , pi+1 at ti−1 , ti , ti+1 , respectively. This construction is illustrated in Fig. 1. The blending concept of Overhauser has various generalizations. Pobegailo [2] blended line segments and circular arcs with polynomials of degree 2n + 1 (specified in power basis) to produce interpolating curve of variable smoothness at data points. Wiltsche [3] studied blending in a wider context, used polynomials like [2] but specified them in Bernstein basis to blend different types of parametric curves, in addition introduced shape parameters. Wenz [4] linearly blended line segments and circular arcs to produce G1 interpolating spline. Szilvási-Nagy and Vendel [5] applied trigonometric blending of line segments and circular arcs to produce G2 interpolating spline. In order to describe an interpolation subspline scheme of degree 2k − 1 and of C k−1 continuity, Röschel [6] uses B-spline basis functions of degree k − 1 (k ≥ 1) for blending Lagrange interpolants of degree k that interpolate k consecutive data points at given strictly increasing parameter values. His method also includes procedures like Overhauser or quintic interpolation schemes as special cases. Replacing the Lagrange by Hermite interpolants, Gfrerrer and Röschel [7] generalize the previous idea. In this case, the generated polynomial subspline interpolates the Hermite input data that consists of parameter values and corresponding zeroth or higher order derivatives, and the common order k of B-spline basis functions used for blending is greater than or equal to the maximal order of derivatives appearing in the interpolation input. These blended Hermite interpolants will also be C k−1 continuous. Schneider [8] uses shape parameters (i.e., tension-like or straightness factors) and special weight functions in order to pull the parabolic arcs of the classical Overhauser interpolation scheme towards the chords between consecutive data points and to ensure the smoothness of the obtained interpolating curves at their common joints, respectively. The advantage of these methods is that there is no need for the solution of linear systems; the interpolating curve can be computed locally, as a consequence of which the spline curve is locally modifiable. A drawback of these procedures is that we lose the control point representation when the arcs to be blended are specified by control points. Our proposed method combines the assets of the two options described above, i.e., it produces the interpolating spline locally and provides directly the control points of the interpolating arcs with arbitrary order of continuity at joints. We demonstrate the method for trigonometric curves, that include polynomial curves as a special case. The advantages of the proposed method are: treats both trigonometric and polynomial local interpolating spline curves in a unified way; provides a ready to use control point based description which is advantageous from data storage and transmission point of view and is suitable for nowadays CAD/CAM systems; it is intuitive and user friendly; possesses both global and local (its rational counterpart) shape parameters; offers a more versatile tool for shape design than polynomial methods (e.g. ensures higher order circular precision), since these trigonometric curves form a proper subset of rational curves. The outline of the rest of the paper is as follows. Section 2 recalls the normalized B-basis of the vector space α F2n = span {cos(kt ), sin(kt ) : t ∈ [0, α ]}nk=0
(2)
of truncated Fourier series of finite order n ≥ 1 (i.e., trigonometric polynomials of degree at most 2n), presents some of its properties that will be used latter and specifies the curve defined by this basis in the form (1). A subdivision formula for the special case n = 1 is also provided for these curves. Section 3 proposes a trigonometric blending function and shows its essential properties. Section 4 constructs and subdivides a single elliptical arc that interpolates three consecutive data points at specified parameter values. Section 5 uses the blending function introduced in Section 3 to produce the arcs of a C n interpolating spline curve of order n + 1, by blending two quadratic trigonometric arcs, the global shape parameter of which is α . Ready to use formulae for the control points of the interpolating arcs are provided. In Section 6 it is shown that when α → 0 trigonometric arcs degenerate to Bézier curves of degree 2n + 1. By assigning non-negative weights of rank 1
248
I. Juhász, Á. Róth / Journal of Computational and Applied Mathematics 263 (2014) 246–261
with data points, Section 7 formulates the rational counterpart of our local interpolating scheme with rational trigonometric spline curves. Conclusions section closes the paper. 2. Trigonometric basis functions and curve We will produce interpolating trigonometric curves, that are described by means of the combination of control points and basis functions, therefore we need a proper basis in the space of trigonometric polynomials. We will use a transformed version of the basis specified in Sánchez-Reyes [9]. Let α ∈ (0, π) be an arbitrarily fixed parameter and consider the function system
Aα2n,i (t ) : t ∈ [0, α ] i=0 =
2n
α 2n−i c2n ,i sin
α−t
2
sini
t
2
2n : t ∈ [0, α ] ,
(3)
i =0
where the normalizing positive coefficients i
α c2n ,i
2 α i−2r n i−r 1 , 2 cos = r 2 sin2n α2 r =0 i − r
i = 0, 1, . . . , 2n
fulfill the symmetry α α c2n ,i = c2n,2n−i , i = 0, 1, . . . , n. Function system (3) is the B-basis of the vector space (2). Therefore, 2n
(4)
Aα2n,i (t ) ≡ 1
(5)
i =0
and Aα2n,i (0) = Aα2n,i (α) =
1, 0,
for i = 0, otherwise,
1, 0,
(6)
for i = 2n, otherwise.
(7)
Basis functions (3) enjoy the symmetry Aα2n,i (α − t ) = Aα2n,2n−i (t ) ,
∀t ∈ [0, α ] , i = 0, 1, . . . , n
(8)
since α
α
A2n,i (α − t ) = c2n,i sin
2n−i
α − (α − t )
sin
2
α i = c2n ,2n−i sin
= Aα2n,2n−i (t ),
α−t
2
sin2n−i
i
α−t
2
t
2
∀t ∈ [0, α ] . The product of basis functions A2n,i (t ) and Aα2m,j (t ) can be expressed in the form α
α α c2n ,i c2m,j α Aα2n,i (t )Aα2m,j (t ) = α A2(n+m),i+j (t ), c2(n+m),i+j
∀t ∈ [0, α ]
(9)
since α
α
α
A2n,i (t )A2m,j (t ) = c2n,i sin
2n−i
α−t
2
α α 2(n+m)−i−j = c2n ,i c2m,j sin
=
α α c2n ,i c2m,j
c2α(n+m),i+j
sin
i
t
2
α−t 2
Aα2(n+m),i+j (t ),
α
2m−j
· c2m,j sin
α−t 2
j
sin
t
2
t
sini+j
2
∀t ∈ [0, α ] .
Using the linear reparametrization t (v) = αv, v ∈ [0, 1], one obtains from (3) the system Aα2n,i (αv) : v ∈ [0, 1] i=0 =
2n
α i c2n ,i sin
αv 2
sin2n−i
α (1 − v) 2
2n : v ∈ [0, 1] , i=0
by means of which it is easy to show that (3) is a generalization of the Bernstein basis of degree 2n.
I. Juhász, Á. Róth / Journal of Computational and Applied Mathematics 263 (2014) 246–261
249
Proposition 2.1 (Bernstein Polynomials of Degree 2n as Special Case). If α → 0 basis functions (3) converge to the Bernstein polynomials of degree 2n, i.e., lim Aα2n,i (αv) = B2n i (v),
∀v ∈ [0, 1] , i = 0, 1, . . . , 2n.
α→0
(10)
Proof. In Ref. [9] functions (3) are obtained by raising the identity Aα2,0 (αv) + Aα2,1 (αv) + Aα2,2 (αv) = 1,
∀v ∈ [0, 1]
to the power n, i.e., 2n
Aα2n,i (αv) = Aα2,0 (αv) + Aα2,1 (αv) + Aα2,2 (αv)
n
= 1,
∀v ∈ [0, 1]
i =0
which in the limiting case α → 0 takes the form
1
sin2 sin2 α2
lim
α→0
α (1 − v)
+
2
2 cos α2
sin sin2 α2
α (1 − v)
2
sin
αv 2
+
1
sin2 sin2 α2
αv n 2
n = (1 − v)2 + 2 (1 − v) v + v 2 = ((1 − v) + v)2n 2n = B2n i (v) i=0
for all values of v ∈ [0, 1]. Corollary 2.1 (Product Rule of Bernstein Polynomials as Special Case). The limiting case α → 0 of the linear reparametrization t (v) = αv, v ∈ [0, 1] of the product rule (9) results in the well-known identity
2n k
2m B2n k (v)Bℓ (v) =
2m
ℓ
2(n+m) k+ℓ
(n+m) B2k+ℓ (v),
∀v ∈ [0, 1]
of Bernstein polynomials for all k = 0, 1, . . . , 2n and ℓ = 0, 1, . . . , 2m. Definition 2.1 (Trigonometric Curves of Order n). Let n ≥ 1 be a finite order. By means of functions (3) and control points
2n
di ∈ Rδ i=0 the curve of order n (degree 2n)
gαn (t ) =
2n
Aα2n,i (t )di ,
t ∈ [0, α ] ,
(11)
i=0
can be defined, where α is a global shape parameter. Curves (11) form a proper subset of rational curves, i.e., any curve of type (11) has a rational parametrization. However, rational parametrizations have drawbacks, e.g., their higher order derivatives become very involved and the parametrization of the curve is often poor (cf. [10]), thus it is worth using trigonometric parametrization. Basis functions (3) form the normalized B-basis of function space (2), therefore there must be a recursive corner cutting algorithm for curves (11), cf. [11]. Unfortunately, the general formula is unknown as yet, there is a scheme only for a very specific arrangement of control points in [12]. We managed to solve the problem only in the special case n = 1 for any placement of control points. We will need this result in Section 4. Let us consider the quadratic curve gα1 (t ) =
2
Aα2,i (t )di ,
t ∈ [0, α ] ,
i=0
and the functions µ10 (·; α) , µ11 (·; α) , µ20 (·; α) : [0, α ] → [0, 1], where
µ (t ; α) = 1 − 1 0
µ (t ; α) = 1 1
, sin α2 cos 2t sin2
t
t , sin α2 sin α2 − sin α− cos 2t 2
µ (t ; α) = 1 − 2 0
t sin α− 2
2
t sin α− cos 2t 2 sin α2
.
(12)
250
I. Juhász, Á. Róth / Journal of Computational and Applied Mathematics 263 (2014) 246–261
For any arbitrarily chosen value β ∈ (0, α) we introduce the subdivision points d10 = 1 − µ10 (β; α) d0 + µ10 (β; α) d1 ,
d11 = 1 − µ11 (β; α) d1 + µ11 (β; α) d2
and d20 = 1 − µ20 (β; α) d10 + µ20 (β; α) d11 ,
using which we define curves β
β
β
β
α−β A2,0 (t )d20
α−β A2,1 (t )d10
g1,ℓ (t ) = A2,0 (t ) d0 + A2,1 (t )d10 + A2,2 (t )d20 , α−β g1,r (t )
=
+
+
t ∈ [0, β ] ,
α−β A2,2 (t )d2 ,
t ∈ [0, α − β ] .
(13) (14)
By means of settings and notations above, one can easily prove the next proposition. Proposition 2.2 (Subdivision of Quadratic Trigonometric Curves). Quadratic trigonometric curves (13) and (14) form two consecutive arcs of curve (12) and the continuity at their common point d20 is C ∞ , i.e., β g1,ℓ (t ) = gα1 (t ), α−β g1,r (t )
∀t ∈ [0, β ] ,
α
= g1 (β + t ) ,
(15)
∀t ∈ [0, α − β ]
(16)
and dz dt
β
g (t ) z 1,ℓ
=
t =β
dz
α−β
g (t ) dt z 1,r
,
∀z ∈ N.
t =0
3. Blending function and the blended curve We are going to describe an interpolation scheme that is based on curve blending, therefore first of all we specify an appropriate blending function, then we show an essential continuity property of the resulted blended curve. 3.1. Trigonometric blending function Definition 3.1 (Trigonometric Blending Function of Degree 2n). By means of basis functions (3) we define the function bα2n (t ) =
2n
1 α A2n,n (t ) + Aα2n,i (t ), 2 i=n+1
t ∈ [0, α ]
(17)
that we will use for blending curves. Blending function (17) is non-negative and strictly increasing, specifically bα2n (0) = 0
(18)
bα2n (α) = 1
(19)
and due to Eqs. (6) and (7). Proposition 3.1 (Symmetry Property of the Blending Function). Blending function (17) is symmetric in the sense bα2n (t ) = 1 − bα2n (α − t ) ,
∀t ∈ [0, α ]
i.e., its graph is symmetric with respect to the point of coordinates α2 , 12 . Proof. Subsequent application of Eqs. (8) and (5) yields 1 − bα2n (α − t ) = 1 −
2n
1 α A2n,n (α − t ) − Aα2n,i (α − t ) 2 i=n+1 1
= 1 − Aα2n,n (t ) − 2
Aα2n,i (t )
i=0 2n
=
n −1
1 α A2n,n (t ) + Aα2n,i (t ) = bα2n (t ) 2 i=n+1
for all values of t ∈ [0, α ].
I. Juhász, Á. Róth / Journal of Computational and Applied Mathematics 263 (2014) 246–261
251
Proposition 3.2 (Vanishing Higher Order Derivatives). Derivatives of blending function (17) vanish at the endpoints up to the order n − 1, i.e., dz dt
bα (t ) z 2n
= 0,
(20)
=0
(21)
t =0
dz
b2n (t ) dt z α
t =α
for all z = 1, 2, . . . , n − 1. Proof. Function system (3) is an extended Chebyshev system, thus result of Mazure [13] on the derivatives of B-basis in extended Chebyshev spaces can be applied, that results in the Eqs. dz
A2n,i (t ) z dt α
dz
=
dt
t =0
A (t ) z 2n,2n−i α
= 0,
0 ≤ z < i ≤ 2n.
(22)
t =α
This immediately yields dz dt
bα (t ) z 2n
= 0,
z = 1, 2, . . . , n − 1.
(23)
t =0
Basis functions (3) form a partition of unity, thus 2n dz α A2n,i (t ) = 0, z dt i=0
∀t ∈ [0, α ] , z ≥ 1.
Taking into account Eqs. (22), (23), we obtain
n −1 dz α A ( t ) 2n,i dt z i=0
+
t =0
1 dz 2 dt
A (t ) z 2n,n α
= 0,
z = 1, 2, . . . , n − 1.
(24)
t =0
Utilization of identities (8) and their derivatives dz dt
Aα (t ) z 2n,i
=
t =α
dz dt
Aα (α − t ) z 2n,2n−i
= (−1)z
t =α
dz dt
Aα (t ) z 2n,2n−i
t =0
α
for i = 0, 1, . . . , n yields the derivative of b2n (t ) at t = α in the form
A ( t ) 2n,i dt z i=n+1 2n dz
α
+
t =α
1 dz 2 dt
A ( t ) z 2n,n α
n −1 dz α 1 dz α = A2n,n (α − t ) A − t + (α ) 2n,i z z dt i=0 2 dt t =α t =α t =α
= (−1)
z
n dz α A2n,i (t ) z dt i=0
t =0
A2n,n (t ) + =0 z 2 dt t =0 1 dz
α
due to Eq. (24). 3.2. Curve blending Let us consider curves c1 (t ) and c2 (t ), t ∈ [0, α ] that fulfill conditions c1 (0) = c2 (0) = p1 ,
(25)
c1 (α) = c2 (α) = p2 . By means of the blending function (17) we generate the curve c(t ) = 1 − bα2n (t ) c1 (t ) + bα2n (t )c2 (t ),
t ∈ [0, α ] .
(26)
It is obvious, that c(t ) is a convex combination of curves c1 (t ) and c2 (t ), and the blended curve fulfills the endpoint interpolation properties
c(0) = p1 , c (α) = p2 .
252
I. Juhász, Á. Róth / Journal of Computational and Applied Mathematics 263 (2014) 246–261
Proposition 3.3 (Higher Order Hermite Conditions of the Blended Curve). If curves c1 (t ) and c2 (t ), t ∈ [0, α ] are at least n times continuously differentiable at t = 0 and t = α , then for the blended curve (26) equalities dz
dz c (t ) = z c1 (t ) , dt z dt t =0 t =0 dz dz c ( t ) c ( t ) = 2 z z dt
dt
t =α
t =α
hold for all z = 1, 2, . . . , n. Proof. The zth order derivative of c(t ) is dz dt
c(t ) = z
z k z d k=0
dt k
k
d z −k
1 − bα2n (t )
dt
c (t ) + z −k 1
dk α d z −k b ( t ) c2 (t ) . 2n dt k dt z −k
At first we assume that 0 ≤ z < n and evaluate this derivative at t = 0. In this case the only non-vanishing term of the above sum is
dz 1 − b2n (0) c (t ) z 1
α
α
dt
+ b2n (0)
t =0
dz dt
c (t ) z 2
t =0
(i.e., the one that belongs to k = 0) due to Proposition 3.2. Thus, we obtain dz dt
c (t ) z
dz dz α = 1 − b2n (0) + b2n (0) z c2 (t ) c 1 (t ) z dt dt t =0 t =0 dz = c (t ) z 1 α
t =0
dt
t =0
because of Eq. (18). Now, we consider the case z = n at t = 0 that yields dn
dn
c (t ) n dt
t =0
= + c1 (t ) n dt t =0 dn = − c1 (t ) dt n t =0 dn c (t ) = n 1 dt
n dn α c1 (0) + d bα (t ) c2 (0) 1 − b ( t ) 2n 2n n n dt dt t =0 t =0
n dn α p1 + d bα (t ) p1 b ( t ) 2n n 2n dt n dt t =0 t =0
t =0
using assumption (25). The equality for the derivative at the endpoint (t = α ) can be proved analogously. 4. An interpolating elliptical arc We are going to use elliptical arcs for local interpolation, therefore we have to describe the elliptical arc as a quadratic trigonometric curve of type (11) that interpolates three consecutive points at specified parameter values. The task is as follows. Given the sequence of data points p0 , p1 , p2 and associated parameter values u0 < u1 < u2 , with the constraint u2 − u0 < π , find those control points q0 , q1 , q2 for which the curve e(u) =
2
α +α1
A2,0j
(u)qj ,
u ∈ [0, α0 + α1 ]
(27)
j=0
passes through the given points at the given parameter values, i.e., equalities e(0) = p0 , e (α0 ) = p1 , e (α0 + α1 ) = p2 are fulfilled, where αi = ui+1 − ui , (i = 0, 1). The solutions are q0 = p0 , α +α1
q1 =: q1 0 q2 = p2 .
=
1 α +α1
A2,01
(α0 )
α +α1
p1 − A2,00
α +α1
(α0 ) p0 − A2,02
(α0 ) p2 ,
I. Juhász, Á. Róth / Journal of Computational and Applied Mathematics 263 (2014) 246–261
253
Fig. 2. The elliptical arc that interpolates data points p0 , p1 , p2 that is split into two arcs. Settings are α = π/2 and knots are chosen by centripetal parametrization.
Using the subdivision formula detailed in Proposition 2.2, we split the elliptical arc (27) at the parameter value u = α0 . Control points of the left and right arcs are as follows. The left arc eℓ (u) =
2
α
A2,0j (u) aj,ℓ ,
u ∈ [0, α0 ]
j =0
is defined by control points a0,ℓ = p0 , α
0 a1,ℓ =: a1,ℓ = 1 − µ10 (α0 ; α0 + α1 ) p0 + µ10 (α0 ; α0 + α1 ) q1 , a2,ℓ = p1 ,
while the right arc er (u) =
2
α
A2,1j (u)aj,r ,
u ∈ [0, α1 ]
j=0
is generated by a0,r = p1 , α
a1,r =: a1,1r = 1 − µ11 (α0 ; α0 + α1 ) q1 + µ11 (α0 ; α0 + α1 ) p2 , a2,r = p2 .
The interpolating elliptical arc and its subdivision are illustrated in Fig. 2. Such left and right elliptical arcs will be used in the forthcoming section to produce C n (n ≥ 1) interpolating trigonometric spline curves. Remark 4.1. If points p0 , p1 , p2 are collinear the process described above results in two straight line segments parametrized in the form (11). 5. Interpolating trigonometric spline curve Our objective is to generalize the Overhauser spline curve by solving the following interpolation problem. Given the sequence of data points {pi }m i=0 with associated parameter values t0 < t1 < · · · tm−1 < tm , moreover the shape parameter α ∈ (0, π). Find the control points of a spline curve of continuity C n (n ≥ 1) that interpolates the given data points, and that consists of arcs of type (11). We intend to use elliptical arcs (27) defined by consecutive triplets of data points pi−1 , pi , pi+1 (i = 1, 2, . . . , m − 1), where the maximum of domains ti+1 − ti−1 is equal to α , therefore we have to scale parameter values ti that results in the new knots uα0 = t0 , uαi = uαi−1 +
α dmax
(ti − ti−1 ) ,
i = 1, 2, . . . , m,
254
I. Juhász, Á. Róth / Journal of Computational and Applied Mathematics 263 (2014) 246–261
where dmax = max {ti+1 − ti−1 : i = 1, 2, . . . , m − 1}. Thus the domain of the interpolating elliptical arc will be αi−1 + αi , α +α with αi = uαi+1 − uαi , (i = 0, 1, . . . , m − 1). Its control points are pi−1 , qi i−1 i , pi+1 , where α
qi i−1
+αi
=
1 α 1 +αi A2,i− 1
α
1 pi − A2,i− 0
(αi−1 )
+αi
α
(αi−1 ) pi−1 − A2,i−2 1
+αi
(αi−1 ) pi+1 .
We split this arc into two (also quadratic) elliptical arcs, as it is described in Section 4, that results in the left and right arcs α
α
α
α
α
i−1 i−1 1 + A2,i−2 1 (u)pi , ei,ℓi−1 (u) = A2,i− 0 (u)pi−1 + A2,1 (u)ai,ℓ
u ∈ [0, αi−1 ]
and α
α
α
α
α
ei,ir (u) = A2,i 0 (u) pi + A2,i 1 (u)ai,ir + A2,i 2 (u)pi+1 , where α
u ∈ [0, αi ] , α
ai,ℓi−1 = 1 − µ10 (αi−1 ; αi−1 + αi ) pi−1 + µ10 (αi−1 ; αi−1 + αi ) qi i−1
α ai,ir
αi−1 +αi
= 1 − µ (αi−1 ; αi−1 + αi ) qi 1 1
Blending
α ei,ir (u)
αi gn+ 1,i (u)
and
α ei+i 1,ℓ (u) (i αi
= 1 − b2n (u)
+αi
,
+ µ (αi−1 ; αi−1 + αi ) pi+1 . 1 1
= 1, 2, . . . , m − 1) we obtain the sequence of arcs
α ei,ir (u)
α
α
+ b2ni (u)ei+i 1,ℓ (u),
u ∈ [0, αi ] , i = 1, 2, . . . , m − 1.
(28) α
α
i+1 i Proposition 5.1 (Continuity of the Interpolating Trigonometric Spline Curve of Order n+1). Arcs gn+ 1,i (u) and gn+1,i+1 (u)(i = 1, 2, . . . , m − 2) are C n continuous at their common joint pi+1 .
Proof. On account of Proposition 3.3 dz
αi gn+ 1,i (u) z du dz du
= u=αi
α ei+i 1,ℓ (u) z du
,
u=αi
α
dz
g i+1 (u) z n+1,i+1
=
u =0
dz du
α
e i+1 (u) z i+1,r
u =0
that are equal for all z ≥ 0 due to Proposition 2.2. α
m−1
i In what follows, we determine the control points of blended arcs gn+ 1,i (u) : u ∈ [0, αi ] i=1 . We have that
α
α
α
α
α
i i i i i gn+ 1,i (u) = 1 − b2n (u) ei,r (u) + b2n (u)ei+1,ℓ (u)
=
2n α 1 α α α α α 1 − A2ni ,n (u) − A2ni ,j (u) A2,i 0 (u)pi + A2,i 1 (u)ai,ir + A2,i 2 (u) pi+1 2 j =n +1
+ =
2n α 1 αi α α α α A2n,n (u) + A2ni ,j (u) A2,i 0 (u)pi + A2,i 1 (u)ai+i 1,ℓ + A2,i 2 (u)pi+1 2 j =n +1
n −1
α 1 α α α α A2n,j (u) + A2ni ,n (u) A2,i 0 (u)pi + A2,i 1 (u) ai,ir + A2,i 2 (u)pi+1 2 j =0
+
αi
2n α 1 αi α α α α A2n,n (u) + A2ni ,j (u) A2,i 0 (u)pi + A2,i 1 (u)ai+i 1,ℓ + A2,i 2 (u)pi+1 . 2 j =n +1
α
By rearranging the right hand side, using identity (9) and collecting the coefficients of basis functions A2(i n+1),j (u) : u ∈ [0, αi ]
2(n+1)
, we obtain control points
j=0
αi αi α α α α c2ni ,j−1 c2,i1 α c2ni ,j−2 c2,i2 c2n,j c2,0 i p + a + pi+1 , for j = 0, 1, . . . , n, i α α α i , r c2(in+1),j c2(in+1),j c2(in+1),j αi α α α α α c2n,n+1 c2,i0 c2ni ,n c2,i1 1 αi c2ni ,n−1 c2,i2 αi α p + a + a + pi+1 , for j = n + 1, dj i = i α α α i,r i+1,ℓ c2(in+1),n+1 c2(in+1),n+1 2 c2(in+1),n+1 αi αi αi αi αi αi c2n,j c2,0 pi + c2n,j−1 c2,1 aαi + c2n,j−2 c2,2 pi+1 , for j = n + 2, n + 3, . . . , 2 (n + 1) , α
c2(in+1),j
α
c2(in+1),j
i+1,ℓ
α
c2(in+1),j
I. Juhász, Á. Róth / Journal of Computational and Applied Mathematics 263 (2014) 246–261
α
α
α
255
α
Fig. 3. The arc g3,i i (u) of the C 2 interpolating spline is obtained by blending elliptical arcs ei,ir (u) and ei+i 1,ℓ (u) with b4 i (u); α = 3 and knots ti are chosen by chord length parametrization.
Fig. 4. The effect of the order of continuity on the shape of an interpolating plane trigonometric spline curve; α = π2 and the initial knots were chosen by centripetal parametrization. Subfigures (a), (b) and (c) show the curve and its signed curvature plot for continuities C 1 , C 2 and C 3 , respectively.
α
where the convention c2ni ,j = 0 if j < 0 or j > 2n is adopted. Fig. 3 shows the construction of an arc of a C 2 interpolating spline along with its control points. Fig. 4 shows interpolating plane trigonometric spline curves of different orders of continuity for a fixed sequence of data points and centripetal parametrization along with their curvature plot. The above described method does not provide the first and last arcs of the interpolating spline, i.e., arcs that connect data points p0 , p1 and pm−1 , pm , respectively. In order to specify these arcs, we can apply, e.g., pseudo data points p−1 , pm+1 , or α α0 we can use arcs e1,ℓ (u), emm−−11,r (u). Both methods guarantee the required order of continuity at p1 and pm−1 , respectively.
256
I. Juhász, Á. Róth / Journal of Computational and Applied Mathematics 263 (2014) 246–261
Fig. 5. The effect of the parametrization of data points (the choice of knots ti ) both on a plane C 3 interpolating spline curve and its signed curvature plot; α = π/2. Observe that the uniform parametrization furnishes unwanted undulations.
Data points have a local influence on the interpolating curve, since the repositioning of data point pi affects at most four
αj
arcs, namely gn+1,j (u) : u ∈ 0, αj
i+1
.
j=i−2
α
α
i i Remark 5.1 (Placement of Control Points). Due to the construction of the arc gn+ 1,i (u), control points dj
are in the plane spanned by data points pi−1 , pi , pi+1 and pi , pi+1 , pi+2 , respectively.
n j =0
α
and dj i
2(n+1) j =n +2
Fig. 5 illustrates the effect of the parametrization (the choice of knots ti ) of data points, while Fig. 6 shows the influence of the global shape parameter α on the shape of the interpolating plane spline curve. Remark 5.2 (Higher Order Circular Precision). Due to the construction and local controllability of the arcs (28), our method is also able to reproduce arcs with circular precision of order n if at least 4 consecutive data points are uniformly distributed on a circle and α is the double of the central angle determined by two consecutive data points on the circle. Fig. 7 illustrates a C 2 interpolating trigonometric spline curve that has four arcs with second order circular precision. This circular precision of order n is the circular counterpart of the linear precision of Bézier curves, i.e., equidistantly placed Bézier points on a line generate a linear polynomial. Circular precision is important in engineering applications as mentioned in articles [14–16] that produce G1 or G2 continuous polynomial interpolating arcs with zeroth order circular precision. 6. Asymptotic property of the interpolating trigonometric spline curve This section examines the limiting case α → 0 of the interpolating C n continuous trigonometric spline curve that consists α
m−1
i of arcs gn+ 1,i (u) : u ∈ [0, αi ] i=1 of order n + 1. Observe that
α
lim αi = lim uαi+1 − uαi = lim (ti − ti−1 ) = 0 α→0 α→0 α→0 dmax
I. Juhász, Á. Róth / Journal of Computational and Applied Mathematics 263 (2014) 246–261
257
Fig. 6. The impact of the global shape parameter α on the shape of a C 2 interpolating spline curve; knots ti are specified by chord length parametrization.
Fig. 7. Data points {pi }8i=2 are uniformly distributed on a circle with central angle π6 of arcs. Parameter settings are n = 2 and α = π3 . Observe that the
6
zeroth, first and second order derivatives of the plane C 2 interpolating trigonometric arcs gα3,i (u) : u ∈ [0, α ] i=3 coincide with the same order derivatives of the support circle. As expected, the signed curvature plot is constant on the interval [u3 , u7 ].
for all i = 0, 1, . . . , m − 1. Using parametrization u(t ) = (αi−1 + αi )
t − ti−1 ti+1 − ti−1
,
t ∈ [ti−1 , ti+1 ]
and Proposition 2.1, we obtain for α → 0, that the local interpolating elliptical arc α
ei i−1
+αi
(u) =
2
α
qj,ii−1
+αi αi−1 +αi A2,0 (u),
u ∈ [0, αi−1 + αi ]
j =0
of order one degenerates to the quadratic Bézier curve (i.e., to the parabolic arc) e0i (t ) :=
2 j =0
where q00,i = pi−1 ,
q0j,i B2j
t − ti−1 ti+1 − ti−1
,
t ∈ [ti−1 , ti+1 ] ,
(29)
258
I. Juhász, Á. Róth / Journal of Computational and Applied Mathematics 263 (2014) 246–261
1
q01,i =
B21
ti −ti−1 ti+1 −ti−1
ti − ti−1 ti − ti−1 pi − pi−1 B20 − pi+1 B22 , ti+1 − ti−1 ti+1 − ti−1
q02,i = pi+1 . Due to equalities lim µ
α→0
1 0
(αi−1 + αi )
t − ti−1 ti+1 − ti−1
; αi−1 + αi
= lim µ α→0
1 1
(αi−1 + αi )
= lim µ20 (αi−1 + αi ) α→0
=
t − ti−1 ti+1 − ti−1
t − ti−1 ti+1 − ti−1 t − ti−1 ti+1 − ti−1
∈ [0, 1] ,
; αi−1 + αi ; αi−1 + αi
∀t ∈ [ti−1 , ti+1 ] ,
the special subdivision scheme detailed in Proposition 2.2 coincides with the classical de Casteljau algorithm of quadratic Bézier curves when α → 0. Splitting the parabola (29) into two (also quadratic) Bézier arcs, one obtains the left and right arcs e0i,ℓ (t ) = B20 (vi−1 (t )) pi−1 + B21 (vi−1 (t )) a0i,ℓ + B22 (vi−1 (t )) pi ,
t ∈ [ti−1 , ti ]
and e0i,r (t ) = B20 (vi (t )) pi + B21 (vi (t )) a0i,r + B22 (vi (t )) pi+1 ,
t ∈ [ti , ti+1 ] ,
where a0i,ℓ = a0i,r =
1−
1−
ti − ti−1
pi−1 +
ti+1 − ti−1 ti − ti−1 ti+1 − ti−1
q01,i +
ti − ti−1 ti+1 − ti−1
q01,i ,
ti − ti−1 pi+1 , ti+1 − ti−1
and
vi ( t ) =
t − ti ti+1 − ti
,
t ∈ [ti , ti+1 ] , i = 1, 2, . . . , m − 1.
Due to Corollary 2.1, by blending parabolas e0i,r (t ) and e0i+1,ℓ (t ) (i = 1, 2, . . . , m − 1, t ∈ [ti , ti+1 ]) one obtains the sequence of Bézier arcs g0n+1,i (t ) = 1 − b02n (t ) e0i,r (t ) + b02n (t )e0i+1,ℓ (t ),
t ∈ [ti , ti+1 ] , i = 1, 2, . . . , m − 1
(30)
of degree 2 (n + 1), where α
b02n (t ) = lim b2ni (t ) α→0
= lim
α→0
=
1 2
2n
α 1 αi A2n,n (αi vi (t )) + A2ni ,k (αi vi (t )) 2 k=n+1
B2n n (vi (t )) +
2n
B2n k (vi (t )) ,
∀t ∈ [ti , ti+1 ]
k=n+1
denotes the blending function of degree 2n. Naturally, Bézier arcs g0n+1,i (t ) and g0n+1,i+1 (t ) (i = 1, 2, . . . , m − 2) are still C n continuous at their joint pi . Observe that the Bézier ordinates
y0 = 0, y1 = 0, . . . , yn−1 = 0, yn =
of the blending function b02n (t ) = of the Bézier function
2n
k=0
1 2
, yn+1 = 1, . . . , y2n = 1
yk B2n k (vi (t )) of apparently even degree 2n, can be obtained by elevating the degree
2n−1
ℓ=0
xℓ Bℓ2n−1 (vi (t )) ,
t ∈ [ti , ti+1 ]
of odd degree 2n − 1 having Bézier ordinates
(x0 = 0, x1 = 0, . . . , xn−1 = 0, xn = 1, xn+1 = 1, . . . , x2n−1 = 1) .
I. Juhász, Á. Róth / Journal of Computational and Applied Mathematics 263 (2014) 246–261
259
Indeed, applying the geometric interpretation of the degree elevation, one has that y0 = x0 = 0,
0, j j 1 yj = x j −1 + 1 − xj = , 2n 2n 2 1, y2n = x2n−1 = 1,
j = 1, 2, . . . , n − 1,
j = n, j = n + 1, n + 2, . . . , 2n − 1,
i.e., the next proposition is proved. Proposition 6.1 (Degree Reduction). The actual degree of the polynomial blending function b02n (t ), t ∈ [ti , ti+1 ] is 2n − 1. The special case n = 1 of the above polynomial is the linear blending function vi (t ) that appears in the classical cubic Overhauser interpolation scheme, therefore our proposed method is a generalization of the classical cubic Overhauser
m−1
interpolation scheme for C n continuity. Naturally, this also means that Bézier arcs g0n+1,i (t ) : t ∈ [ti , ti+1 ] i=1 are of degree 2n + 1 and the control points of the ith arc can be expressed as
2n−1 2 2n−1 2 2n−1 2 j−1 1 j − 2 j 0 2 0 pi + ai,r + pi+1 , for j = 0, 1, . . . , n, 2n+1 2n+1 2n+1 j j j d0j = 2n−1 2 2n−1 2 2n−1 2 j 0 j−1 1 j − 2 2 0 pi + ai+1,ℓ + pi+1 , for j = n + 1, n + 2, . . . , 2n + 1. 2n+1 2n+1 2n+1 j
j
j
ti m i =0
Moreover, in this special case, original knots { } the shape parameter α .
do not have to be rescaled, since the Bernstein polynomials do not possess
7. Interpolating rational trigonometric spline curve m We can assign non-negative weights {wi }m i=0 of rank 1 to data points {pi }i=0 that will serve as additional local shape parameters of the interpolating spline curve. This weighted interpolating spline curve can be described by means of the rational counterpart of basis functions (3), i.e., by ωi Aα2n,i (t ) Rα2n,i (t ) = , t ∈ [0, α ] , i = 0, 1, . . . , 2n, 2n
j =0
ωj Aα2n,j (t )
2n
where ωj j=0 are non-negative weights of rank 1 that are unknown at the moment. A rational curve in Rδ can be considered as the central projection of a curve in the δ + 1 dimensional space from the origin on to the δ dimensional hyperplane w = 1 (assuming that the last coordinate of space Rδ+1 is denoted by w ). The space Rδ+1 is called the pre-image space. It facilitates our calculations if we do them in the pre-image space, where coordinates of data points are
wi pi , wi
i = 0, 1, . . . , m.
In the pre-image space calculations go like in Section 5. Control points of the interpolating spline in the pre-image space are
α
dj i
αi
ωj
=
α
α
α
α
α
α
α
c2ni ,j−1 c2,i1 α c2ni ,j−2 c2,i2 i w p + a + wi+1 pi+1 i i α α α i , r i i c c2(n+1),j c2(in+1),j 2(n+1),j , αi αi αi αi αi αi c2n c c c c c ,j 2,0 2n,j−1 2,1 2n,j−2 2,2 αi wi + αi wi,r + αi wi+1 αi c2(n+1),j c2(n+1),j c2(n+1),j
c2ni ,j c2,i0
α
c2ni ,n+1 c2,i0
α
wi pi +
c αi 2(n+1),n+1 αi αi c2n ,n+1 c2,0
α
α
α
for j = 0, 1, . . . , n,
c2ni ,n−1 c2,i2 c2ni ,n c2,i1 1 αi α ai+1,ℓ + ai,ir + αi wi+1 pi+1 αi c2(n+1),n+1 2 c2(n+1),n+1 , α α α α c2ni ,n c2,i1 1 αi c2ni ,n−1 c2,i2 αi
for j = n + 1, wi + αi wi+1,ℓ + wi,r + αi wi+1 α c2(in+1),n+1 c2(n+1),n+1 2 c2(n+1),n+1 αi αi α α α α c2n,j c2,0 c2ni ,j−1 c2,i1 α c2ni ,j−2 c2,i2 i w p + a + w p i i i + 1 i + 1 α α i+1,ℓ c αi c2(in+1),j c2(in+1),j 2(n+1),j , for j = n + 2, n + 3, . . . , 2n + 2. α α α α α α c2ni ,j c2,i0 c2ni ,j−1 c2,i1 α c2ni ,j−2 c2,i2 i w i + αi wi+1,ℓ + αi wi+1 αi c2(n+1),j
c2(n+1),j
c2(n+1),j
260
I. Juhász, Á. Róth / Journal of Computational and Applied Mathematics 263 (2014) 246–261
Fig. 8. Local shape modification of a C 2 rational interpolating trigonometric curve, the weight of the data point marked with filled disc takes the values 0.3, 1, 2 while the rest of the weights is 1; α = π/2 and knots ti are specified by centripetal parametrization.
In the model space control points and associated weights are
1 αi
ωj
α
m
dj i j=0
α
and ωj i
m j=0
, respectively. In a design situation
the user has to specify just data points pi and associated non-negative weights wi , control points and weights of the rational interpolating arcs are calculated automatically. The user can adjust the curvature of the interpolating spline curve at a data point locally, by altering the corresponding weight. This kind of shape modification is shown in Fig. 8. 8. Conclusions A method is proposed that produces interpolating trigonometric C n continuous (n ≥ 1) spline curves to a given sequence of data points. Arcs of the spline curve are computed locally, i.e., there is no need for the solution of a system of linear equations, and explicit formulae are provided for the control points of the arcs. Arcs are determined by blending two quadratic trigonometric curves (two elliptical arcs) with an appropriate trigonometric function. The domain of each arc depends on α , where α ∈ (0, π) is a global shape parameter of the interpolating spline curve. It is also shown that associating non-negative weights with data points rational trigonometric interpolating splines can be obtained, that are locally modifiable. The method is derived for trigonometric curves, that also include polynomial curves as the special case α → 0. The proposed interpolation scheme is a genuine generalization of the Overhauser spline. The generalizations are: the order of continuity n ≥ 1, the more versatile shape and the additional global shape parameter. A further direction of the research is to extend this method to smooth surface interpolation. Acknowledgments Imre Juhász carried out his research in the framework of the Center of Excellence of Mechatronics and Logistics at the University of Miskolc. The research of Ágoston Róth was supported by the European Union and the State of Hungary, cofinanced by the European Social Fund in the framework of TÁMOP-4.2.4.A/2-11/1-2012-0001 ‘National Excellence Program’. All assets, infrastructure and personnel support of his research was contributed by the Romanian national grant CNCSUEFISCDI/PN-II-RU-TE-2011-3-0047. The authors also thank the reviewers for their useful comments. References [1] A.W. Overhauser, Analytic Definition of Curves and Surfaces by Parabolic Blending, Technical Report No. SL68-40, Ford Motor Company Scientific Laboratory, 1968. [2] A.P. Pobegailo, Local interpolation with weight functions for variable smoothness curve design, Comput.-Aided Des. 23 (8) (1991) 579–582. [3] A. Wiltsche, Blending curves, J. Geom. Graph. 9 (1) (2005) 67–75. [4] H.-J. Wenz, Interpolation of curve data by blended generalized circles, Comput. Aided Geom. Des. 13 (8) (1996) 673–680. [5] M. Szilvási-Nagy, T.P. Vendel, Generating curves and swept surfaces by blended circles, Comput. Aided Geom. Des. 17 (2) (2000) 197–206.
I. Juhász, Á. Róth / Journal of Computational and Applied Mathematics 263 (2014) 246–261
261
[6] O. Röschel, An interpolation subspline scheme related to B-spline techniques, in: Computer Graphics International, 1997. Proceedings, IEEE, 1997, pp. 131–136. [7] A. Gfrerrer, O. Röschel, Blended Hermite interpolants, Comput. Aided Geom. Des. 18 (9) (2001) 865–873. [8] W. Schneider, A simple technique for adding tension to parabolic blending interpolation, Comput. Math. Appl. 12 (11) (1986) 1155–1160. [9] J. Sánchez-Reyes, Harmonic rational Bézier curves, p-Bézier curves and trigonometric polynomials, Comput. Aided Geom. Des. 15 (9) (1998) 909–923. [10] L. Piegl, On NURBS: a survay, Comput. Graph. Appl. 11 (1) (1991) 55–71. [11] E. Mainar, J.M. Peña, Corner cutting algorithms associated with optimal shape preserving representations, Comput. Aided Geom. Des. 16 (9) (1999) 883–906. [12] J. Sánchez-Reyes, Single-valued curves in polar coordinates, Comput.-Aided Des. 22 (1) (1990) 19–26. [13] M.-L. Mazure, Chebyshev–Bernstein bases, Comput. Aided Geom. Des. 16 (7) (1999) 649–669. [14] G. Farin, Geometric Hermite interpolation with circular precision, Comput.-Aided Des. 40 (4) (2008) 476–479. [15] D.J. Walton, D.S. Meek, G2 Hermite interpolation with circular precision, Comput.-Aided Des. 42 (9) (2010) 749–758. [16] Y. Li, C. Deng, C-shaped G2 Hermite interpolation with circular precision based on cubic PH curve interpolation, Comput.-Aided Des. 44 (11) (2012) 1056–1061.