Interpolation with piecewise quadratic visually C^{2} Bézier polynomials

Robert Schaback

Abstract. For data satisfying certain generalized convexity conditions, the existence and uniqueness of a piecewise quadratic curvature continuous parametric spline interpolant in IR^{2} is proved. Convexity is preserved by the interpolation, and the numerical construction of the interpolant can be carried out efficiently. The interpolation is 4–th order accurate, and a number of examples shows both the applicability and the limitations of this interpolation scheme.

1 Introduction

For given data points

f_{1},…, f_{n}, n ∈ IR^{2}, n ≥ 2, f_{i} ≠ f_{i+1} for 1 ≤ i ≤ n − 1  (1)

we want to construct a continuous interpolating curve

s : I → IR^{2}, I ⊂ IR a closed interval

whose pieces are of class C^{2} and whose tangent direction and curvature vary continuously. This property is called “visual” or “geometric” C^{2} continuity and is abbreviated by GC^{2} (see e.g. [Boehm et. al., ’84] and [Farin, ’82]). The interpolant should be convexity preserving in order to avoid unwanted oscillations. Since nonparametric strictly convex or concave C^{2} functions have no inflection points and nonvanishing curvature, we introduce a natural common generalization of convexity/concavity to the parametric situation:

Definition 1.1 A GC^{2} curve s in IR^{2} is (strictly) noninflecting, if the curvature of s nowhere vanishes.

If the cross–product for planar vectors is defined as

x × y = −x_{1}y_{2} + x_{2}y_{1} for x = ( \begin{matrix}x_{1} \\ x_{2} \end{matrix} ) , y = ( \begin{matrix}y_{1} \\ y_{2} \end{matrix} ) ,

the (signed) curvature of a planar GC^{2} curve f is f′(t) × f′′(t)/∥f′(t)∥^{3} at f(t). The sign of the curvature indicates whether the osculating circle lies to the “left” or “right” of the curve in the sense of a moving observer looking into the direction of the tangential vector. Reversing the parametrization will change the sign of this form of curvature.

The following theorem shows that parametric quadratic splines are appropriate tools for constructing noninflecting curves:

Theorem 1.1 A visual C^{2} parametric quadratic spline curve is strictly noninflecting or a straight line.

Proof. It is easy to verify that a quadratic polynomial in Bézier form either has identically zero curvature (in case it is a straight line) or has nonvanishing curvature.

The sign of the curvature of a quadratic polynomial s(t) in Bézier form with control points b_{0}, b_{1}, b_{2} in IR^{2} and

s(t) = b_{0}(1 − t)^{2} + 2b_{1}t(1 − t) + b_{2}t^{2}, t ∈ [0, 1]

is determined by the position of b_{1} relative to b_{0} and b_{2}. If b_{1} is collinear to b_{0} and b_{2}, the curve is linear. Otherwise, b_{1} lies in the “left” or “right” half–space defined by the oriented line L(b_{0}, b_{2}), thus determining the sign of the curvature.

2 Solvability conditions

In the following we shall construct the interpolating curve by patching together n−1 quadratic polynomials in Bézier form with control points f_{i} = b_{2i}, b_{2i+1}, b_{2i+2} = f_{i+1} for 1 ≤ i ≤ n−1. The interior control points b_{2i+1} will be determined by continuity requirements only. As a result of the above discussion, we can assume for the rest of this paper that the curve is always “turning to the right”. For “left–turning” data the indexing can be reversed to get “right–turning” data (positive curvature). Then the points b_{2i+1} have to be constructed on the “left” side of the oriented lines L(b_{2i}, b_{2i+2}) = L(f_{i}, f_{i+1}), except for the trivial case that all data points f_{i} are collinear, which we exclude from now on.

Since each quadratic piece has two degrees of freedom and visual C^{2} continuity imposes two conditions on each breakpoint f_{2}, f_{3},…, f_{n−1}, we have to fix two additional parameters via boundary conditions. This can be done by adding two points b_{0} ≠ b_{2} and b_{2n+2} ≠ b_{2n} defining tangent directions at b_{2} and b_{2n}, respectively. The complete interpolation problem then will be as follows:

Construct a visually C^{2} curve consisting piecewise of quadratic polynomials in Bézier form with control points f_{i} = b_{2i}, b_{2i+1}, b_{2i+2} = f_{i+1} 1 ≤ i ≤ n − 1, and satisfying boundary conditions

b_{3} − b_{2} = ρ_{1}(b_{2} − b_{0}),  ρ_{1} > 0 arbitrary,
b_{2n} − b_{2n−1} = ρ_{n}(b_{2n+2} − b_{2n}), ρ_{n} > 0 arbitrary,

by finding suitable control points b_{2i+1}, 1 ≤ i ≤ n − 1.

Since curvature continuity is a nonlinear condition and no additional degrees of freedom are available, the problem is equivalent to a nonlinear system of equations (see (2)). It is not always solvable (see Figure 1), even if the data form a piecewise linear curve that is a limit of noninflecting curves. We therefore must look for a weak sufficient condition for solvability. Furthermore, we shall see that unicity of the interpolant requires an even stronger condition than the one that will be used for the existence proof.

[figure: T_{2}; T_{1} = T_{3}; b_{2}; b_{0}; b_{8}; b_{6}; b_{4}; γ_{1}; γ_{3}; γ_{2}]

Figure 1: Unsolvable problem, n = 3

To discuss solvability conditions we introduce tangents T_{i} at interpolation points f_{i} = b_{2i}, 1 ≤ i ≤ n, where the boundary conditions force T_{1} and T_{n} to coincide with the lines L(b_{0}, b_{2}) and L(b_{2n}, b_{2n+2}). respectively. Then b_{2i+1} is the intercept of T_{i} with T_{i+1}. The problem posed in Figure 1 is unsolvable, because T_{2} has to intersect both T_{1} and T_{3} in finite points b_{3} and b_{5} which must lie on the “left” of the lines L(b_{2}, b_{4}) and L(b_{4}, b_{6}), respectively.

[figure: b_{10}; x; γ_{4}; b_{8}; 10; γ_{3}; b_{6}; 10; y; γ_{2}; b_{4}; 10; γ_{1}; b_{2}; b_{0}]

Figure 2: Partially unsolvable problem

Another example is given by Figure 2. This problem is unsolvable for \sqrt{x} + \sqrt{y} < 2\sqrt{10}. The piecewise linear interpolant of b_{0}, b_{2},…, b_{10} takes two successive 90^{∘} turns and fails to satisfy the second part of Definition 2.1 A data set b_{0}, b_{2},…, b_{2n+2} is called noninflecting, if all angles

γ_{i} := ^{}(b_{2i+2} − b_{2i}, b_{2i} − b_{2i−2}), (1 ≤ i ≤ n),

measured in counterclockwise direction, are in (0, 180^{∘}) or in (180, 360^{∘}) (i.e. the data are always “turning right” or always “turning left”). If the sum of two consecutive angles γ_{i} always is in (0, 180^{∘}) or in (180, 360^{∘}), the data set is called nonreflecting.

Note that any nonreflecting and noninflecting data set admits n−1 “accompanying triangles” formed by b_{2i}, b_{2i+2} and the intercept b̃_{2i+1} of the lines L(b_{2i−2}, b_{2i}) and L(b_{2i+2}, b_{2i+4}) for 1 ≤ i ≤ n − 1 (see Figure 4). Each accompanying triangle must contain a control point b_{2i+1} of the solution, if the latter exists. The main existence result of this paper will be

Theorem 2.1 For any noninflecting and nonreflecting data set there is at least one visually C^{2} interpolant by a parametric quadratic spline.

The proof is contained in the next two sections.

Note that Figure 2 can be made into a solvable problem by small perturbations of b_{2} and b_{8}. However, as Figure 3 shows, a perturbed problem with noninflecting and nonreflecting data may have multiple solutions. To get uniqueness, we therefore need additional conditions.

Definition 2.2 A nonreflecting and noninflecting data set is called nondegenerate, if all accompanying triangles have obtuse angles at the points b̃_{2i+1}.

An equivalent form of this condition is

Either all γ_{i} satisfy γ_{i} ∈ (0, 90^{∘}) and γ_{i} + γ_{i+1} < 90^{∘}
or all γ_{i} satisfy γ_{i} ∈ (270^{∘}, 360^{∘}) and 720^{∘} − γ_{i} − γ_{i+1} < 90^{∘},

and in short: A nondegenerate data set

• always turns to the left or always turns to the right,

• the direction changes only by acute angles,

• two successive direction changes sum up to an acute angle.

We remark that this form of nondegeneracy is a natural notion:

Theorem 2.2 An ordered and sufficiently dense sample of sufficiently many points from a noninflecting C^{2} curve will form a nondegenerate data set.

Our main result concerning uniqueness will be

Theorem 2.3 A nondegenerate data set has a unique visually C^{2} piecewise quadratic interpolant.

[figure: (2.5, −2.5); b_{0}; b_{2}; (0, 0); (1, 10); b_{4}; (10, 10); b_{6}; (10.1, 0); b_{8}; (7.6, −2.5); b_{10}]

Figure 3: Data with two solutions (a third solution is suppressed)

3 Principles of construction

Consider part of a noninflecting and nonreflecting data set as in Figure 5. We introduce tangents T_{i} at b_{2i} forming angles α_{i} with b_{2i+2} − b_{2i}. The intercept b_{2i+1} of T_{i} and T_{i+1} must lie in the accompanying triangle b_{2i}b_{2i+2}b̃_{2i+1} to avoid unwanted inflection points. Curvature at b_{2i} is (up to a constant factor) equal to

κ_{i}^{−} = \frac{dist (b_{2i−2}, T_{i})}{∥b_{2i−1} − b_{2i}∥_{2}^{2}} from the left and

κ_{i}^{+} = \frac{dist (b_{2i+2}, T_{i})}{∥b_{2i+1} − b_{2i}∥_{2}^{2}} from the right

(see e.g. [Boehm et. al., ’84, p. 11]).

Using the obvious identities

\frac{dist (b_{2i+2}, T_{i})}{h_{i}} = sin α_{i}

\frac{dist (b_{2i−2}, T_{i})}{h_{i−1}} = sin (γ_{i} − α_{i})

[figure: b_{0}; b_{2}; b_{4}; b̃_{3}; γ_{1}; b_{6}; b̃_{5}; γ_{2}; b_{8}; b̃_{7}; γ_{3}; b̃_{9}; b_{10}; γ_{4}; b_{12}; γ_{5}]

Figure 4: Accompanying triangles

\frac{dist (b_{2i−2}, T_{i})}{∥b_{2i−2} − b_{2i−1}∥_{2}} = sin δ_{i−1} = sin(γ_{i} − α_{i} + α_{i−1})
\frac{dist (b_{2i}, T_{i−1})}{∥b_{2i−1} − b_{2i}∥_{2}} = sin δ_{i−1}

we can express curvature continuity at b_{2i}, 2 ≤ i ≤ n − 1, by the equation

\frac{1}{h_{i}} sin α_{i} \frac{sin^{2}(γ_{i+1} − α_{i+1} + α_{i})}{sin^{2}(γ_{i+1} − α_{i+1})} = \frac{1}{h_{i−1}} sin(γ_{i} − α_{i}) \frac{sin^{2}(γ_{i} − α_{i} + α_{i−1})}{sin^{2} α_{i−1}}  (2)

involving only constants h_{j}, γ_{j} and the variables α_{j}. This gives n − 1 nonlinear equations for the unknowns α_{2},…, α_{n−1}; we use α_{1} = γ_{1} and α_{n} = 0 to satisfy the boundary conditions.

The approach via the system (2) is used later for efficient numerical treatment of the problem. It does not yield much insight for the proof of existence and uniqueness of a solution.

For the latter we theoretically investigate a numerically hazardous “shooting strategy”:

Start. Choose an arbitrary point b_{3} between b_{2} and b̃_{3} on T_{1}.

Iteration. For 2 ≤ i ≤ n − 1, if b_{2i−1} and T_{i} are given, construct b_{2i+1} on T_{i} by using curvature continuity at b_{2i} in the form

∥b_{2i+1} − b_{2i}∥_{2}^{2} = \frac{dist (b_{2i+2}, T_{i})}{dist (b_{2i−2}, T_{i})} ∥b_{2i−1} − b_{2i}∥_{2}^{2}.

If this quantity comes out to be too large, b_{2i+1} moves out of the accompanying triangle formed by b_{2i}, b_{2i+2}, b̃_{2i+1}, and the process can not be repeated. Otherwise the tangent T_{i+1} is fixed by b_{2i+1} and b_{2i+2}, and the next value of i can be used.

[figure: b_{2i−4}; b_{2i−2}; b_{2i}; h_{i−1}; b̃_{2i−1}; γ_{i−1}; T_{i}; T_{i−1}; α_{i−1}; b_{2i−1}; δ_{i−1}; b_{2i+2}; h_{i}; b̃_{2i+1}; γ_{i}; δ_{i}; b_{2i+1}; T_{i+1}; α_{i}; b_{2i+4}]

Figure 5: Local construction

This process tries to get through to the triangle b_{2n−2} b_{2n} b̃_{2n−1} and to hit the line L(b_{2n}, b̃_{2n−1}) to produce a curvature continuous solution. It is comparable to the “single shooting” method for solving a boundary value problem in ordinary differential equations. Its (theoretical) success is equivalent to the solvability of the problem, but experience shows that it is numerically extremely sensitive. There are simple nondegenerate data sets for which 32–bit precision does not suffice to find some b_{3} to shoot through to the last triangle. Replacement by a “multiple shooting” method is equivalent to solving the system (2), which is recommended for practical purposes (see section 6).

4 Existence proof

We now proceed to study the topological properties of a single step of the shooting method defined in the last section.

Let two lines L_{1} and L_{2} intersect in some point S with an angle α of less than 180 degrees. This angle forms two open cones C_{1} and C_{2}; let A_{1} ≠ S and A_{2} ≠ S be points on L_{1} and L_{2}, respectively, but on different sides of S (see Figure 6).

Definition 4.1 A continuous curve z : [0, 1] → IR^{2} is “transversal in C_{1}”, if (see Figure 6)

[figure: S; L_{2}; L_{1}; C_{1}; C_{2}; α; A_{1}; A_{2}; z(t)]

Figure 6: Transversal curve

a) z((0, 1)) ⊂ C_{1},

b) z(0) = A_{1},

c) z(1) ∈ L_{2},

d) z(1) ≠ S. Definition 4.2 On the configuration of Figure 6 we define a mapping F : C_{1} → C_{2} by

a) F (x) ∈ C_{2} lies on the line L(x) through x and S

b) \frac{dist (A_{1}, L(x))}{∥x − S∥_{2}^{2}} = \frac{dist (A_{2}, L(x))}{∥F (x) − S∥_{2}^{2}}

This definition makes sure that the piecewise quadratic polynomial defined in Bézier form by control points A_{1}, x, S and S, F (x), A_{2} is curvature continuous at S. Therefore F corresponds to a single step of the shooting method. If we introduce polar coordinates as in Figure 7, we get for the image F (r, φ) = (R(r, φ), φ) of (r, φ) the expression

R^{2}(r, φ) = \frac{h_{2}}{h_{1}} \frac{sin φ}{sin (α − φ)} r^{2}  (3)

which shows that F is a smooth function. This follows from

\frac{h_{1}}{r^{2}} = \frac{h_{2}}{R^{2}} , h_{1} = h_{1} sin(α − φ), h_{2} = h_{2} sin φ.

Lemma 4.1 Let z be a transversal curve in C_{1}. Then:

[figure: S; L_{2}; L_{1}; C_{1}; C_{2}; α; A_{1}; A_{2}; L(x); F (x); R; x; r; φ; h_{1}; h_{2}; h_{2}; h_{1}]

Figure 7: Construction of F (x) in polar coordinates

a) u(t) := F (z(t)) is a continuous curve on (0, 1) in C_{2}.

b) u(1) := S makes u continuous on (0, 1].

c) lim_{t→0} ∥u(t)∥_{2} = ∞.

d) lim_{t→0} dist(u(t), L_{1}) = 0.

Proof. Everything follows easily from elementary continuity arguments using the properties of z and F described in the polar coordinate system of (3) and Figure 7.

We now observe that in the configuration of Figure 8 two instances of Figure 7 are intertwined. A transversal curve z_{i−1} in the cone C_{i}^{−} will map by Lemma 4.1 into a curve z_{i} in the cone C_{i}^{+}. This curve starts in b_{2i} and has the line L(b_{2i−2}, b_{2i}) as an asymptote. Now the intersection of this curve with the accompanying triangle b_{2i} b_{2i+2} b̃_{2i+1} is transversal in the cone C_{i+1}^{−}. Applying Lemma 4.1 again we get a curve z_{i+1} in the cone C_{i+1}^{+} and so on. The shooting process starts with z_{1} being part of T_{1}.

To prove Theorem 2.1 we now use induction on the statement

“There is an open nonvoid interval on the line L(b_{2}, b̃_{3}) that maps uniquely onto a curve z_{i} in C_{i}^{+} that starts from b_{2i} and has the line L(b_{2i−2}, b_{2i}) as an asymptote.”

This statement follows from the previous lemma for i = 2, and the discussion illustrated by Figure 8 gives the induction step.

The resulting curve z_{n−1} in C_{n−1}^{+} necessarily intersects the line T_{n} given by the boundary condition at b_{2n} at least once. This proves existence of a solution.

Uniqueness is not guaranteed, because z_{n−1} might intersect T_{n} in several points. This can be actually observed by plotting the curves z_{i}, and this is how the solutions of Figure 3 were found.

[figure: b_{2i−2}; b_{2i}; b_{2i+2}; b_{2i+4}; C_{i}^{−}; C_{i+1}^{−}; C_{i}^{+}; C_{i+1}^{+}; z_{i−1}; b̃_{2i+1}; z_{i}; z_{i+1}; α_{i}; α_{i+1}]

Figure 8: Shooting method

5 Uniqueness

Definition 5.1 A transversal curve z in C_{1} is strictly monotonic, if in the polar coordinates of Figure 7, the radius r(φ) is strictly monotonic.

Lemma 5.1

a) For the image (R(r(φ), φ), φ) of a monotonic transversal curve z = (r(φ), φ) in C_{1} the radius R(r(φ), φ) is strictly monotonic.

b) If B_{2} ≠ S is a point on the line L_{1} and the boundary of C_{2} such that the angle γ = ^{}(A_{2}B_{2}S) is obtuse, the image curve F (z) meets the line L(B_{2}, A_{2}) in exactly one point (see Figure 9).

Proof. The first assertion is a direct consequence of (3). If F (z) meets \overline{B_{2}A_{2}} in two points P_{1}, P_{2} with φ_{1} < φ_{2} (see Figure 10), then R(r(φ), φ) is no monotonic function of φ.

For a nondegenerate data set we can apply Lemma 5.1 in each step of the shooting method, because all angles at b̃_{2i+1} are obtuse. Then the curves constructed in the previous section have unique intersections with the sides b̃_{2i+1} b_{2i+2} of the triangles b_{2i} b̃_{2i+1} b_{2i+2}. This also holds for i = n − 1, proving uniqueness.

[figure: S; L_{2}; L_{1}; C_{1}; C_{2}; α; A_{1}; A_{2}; F (z); z; γ; B_{2}]

Figure 9: Monotonicity

[figure: S; A_{1}; A_{2}; F^{−1}(P_{2}); F^{−1}(P_{1}); P_{1}; P_{2}; B_{2}; φ_{2}; φ_{1}; γ]

Figure 10: Uniqueness

6 Convergence

Let f : [a, b] ⟶ IR^{2} be a smooth planar curve with positive curvature, parametrized by arclength. This section analyzes the convergence of the interpolation process for data from f, i.e. for data density tending to zero. We can follow [deBoor et. al., ’87] to assume a local two–dimensional representation of f around a single point by

f(0) = ( \begin{matrix}0 \\ 0 \end{matrix} ) , f′(t) = ( \begin{matrix}cos θ(t) \\ sin θ(t) \end{matrix} )

θ(t) = θ_{1}t + θ_{2}t^{2} + O(t^{3})

without loss of generality. Curvature at f(t) is

κ(t) = θ′(t) = θ_{1} + \frac{1}{2}θ_{2}t + O(t^{2}).

The angle α(t) between the tangent at f(0) and the chord

f(t) − f(0) = ∫_{0}^{t} f′(σ)dσ = \frac{1}{6} ( \begin{matrix}6t − θ_{1}^{2}t^{3} \\ 3θ_{1}t^{2} + 2θ_{2}t^{3} \end{matrix} ) + O(t^{4})  (4)

can be expanded via

tan α(t) = \frac{3θ_{1}t^{2} + 2θ_{2}t^{3} + O(t^{4})}{6t − θ_{1}^{2}t^{3} + O(t^{4})}

as

α(t) = tan α(t) + O(t^{3}) = \frac{1}{2}θ_{1}t + \frac{1}{3}θ_{2}t^{2} + O(t^{3})  (5)

for t → 0. The chord length

h(t) = ∥f(t) − f(0)∥_{2} = t − \frac{1}{24}θ_{1}^{2}t^{3} + O(t^{4})

can replace arc length t up to terms of order O(t^{3}). This allows expansions like (5) to be written in terms of chord length h instead of arc length:

Lemma 6.1 The angle α = α(f, t_{0}, t_{1}) between the tangent to a smooth planar curve f at f(t_{0}) and a chord f(t_{1}) − f(t_{0}) of length h = ∥f(t_{1}) − f(t_{0})∥ is

α = \frac{1}{2}κ(t_{0})h + \frac{1}{3}κ′(t_{0})h^{2} + O(h^{3})  (6)

for t_{1} > t_{0}, t_{1} → t_{0}, where κ is the curvature of the curve.

The expansion (6) is independent of the parametrization except for the orientation hidden in the correspondence of κ′ and t_{1} > t_{0}. For t_{1} < t_{0}, the analogous expansion is

α = −\frac{1}{2}κ(t_{0})h + \frac{1}{3}κ′(t_{0})h^{2} + O(h^{3}),  (7)

obtained from Lemma 6.1 by reversing the parametrization and the sign of curvature. Note that in both cases the angle has to be taken as oriented from the tangent to the chord.

Similarly, one easily verifies expansions of angles between successive chords in terms of the curvature at the central or a non–central point:

Lemma 6.2 The angle γ = γ(f, t_{0}, t_{1}, t_{2}) between two chords f(t_{i+1}) − f(t_{i}), i = 0, 1, of a planar curve f with chord lengths h_{i} = ∥f(t_{i+1}) − f(t_{i})∥ is

\begin{aligned}γ & = & \frac{1}{2}κ(t_{1})(h_{0} + h_{1}) + \frac{1}{3}κ′(t_{1})(h_{1}^{2} − h_{0}^{2}) + O(max(h_{0}^{3}, h_{1}^{3})) \\ & = & \frac{1}{2}κ(t_{0})(h_{0} + h_{1}) + \frac{1}{3}κ′(t_{0})(2h_{0}^{2} + 3h_{0}h_{1} + h_{1}^{2}) + O(max(h_{0}^{3}, h_{1}^{3})) \end{aligned}  (8)

for t_{2} > t_{1} > t_{0}, t_{2} − t_{0} → 0, where κ is the curvature of the curve. The angle between the tangent at f(t_{1}) and the chord f(t_{1}) − f(t_{0}) is

α = \frac{1}{2}κ(t_{0})h_{0} + \frac{2}{3}κ′(t_{0})h_{0}^{2} + O(h_{0}^{3})  (9)

in terms of the curvature at t_{0}.

Proof: Expansion (9) simply follows from (8) for t_{2} = t_{1}, h_{1} = 0, and the first representation of γ in (6) can be written as the difference between (6) and (7). To obtain the non–central representation, the angle between the chord f(t_{2})−f(t_{1}) and the tangent at f(t_{0}) is evaluated as

\frac{1}{2}θ_{1}(h_{0} + h_{2}) + \frac{1}{3}θ_{2}(h_{0}^{2} + h_{0}h_{2} + h_{2}^{2}) + O(h_{2}^{3})

via an expansion of the chord like (4), using h_{2} = ∥f(t_{2}) − f(t_{0})∥. The rest follows easily, using h_{2} = h_{0} + h_{1} + O(max(h_{0}^{3}, h_{1}^{3})).

We now go over to the interpolation of f at points f_{i} = b_{2i} = f(t_{i}) for 1 ≤ i ≤ n with boundary conditions using the directions of f′(t_{1}) and f′(t_{n}). The density of knots is measured by

\begin{aligned}h & := & max_{1≤i≤n−1} h_{i}, \\ h_{i} & := & ∥f(t_{i+1}) − f(t_{i})∥, 1 ≤ i ≤ n − 1, h_{0} := h_{n} := 0, \end{aligned}

and we consider the asymptotic situation h → 0, n → ∞. The data–defined chord angles γ_{i} and the unknown exact tangent angles α_{i}^{∗} are representable as

\begin{aligned}γ_{i} & = & \frac{1}{2}κ_{i}(h_{i} + h_{i−1}) + \frac{1}{3}κ′_{i}(h_{i}^{2} − h_{i−1}^{2}) + O(max(h_{i−1}^{3}, h_{i}^{3})), \\ α_{i}^{∗} & = & \frac{1}{2}κ_{i}h_{i} + \frac{1}{3}κ′_{i}h_{i}^{2} + O(h_{i}^{3}), \end{aligned}

using Lemmas 6.1 and 6.2 and setting κ_{i} := κ(t_{i}), κ′_{i} := κ′(t_{i}). As a first application we eliminate the unknown curvature κ_{i} and get numerically available tangent angles

\begin{aligned}α̂_{i} & := & γ_{i}h_{i}/(h_{i−1} + h_{i}) \\ & = & \frac{1}{2}κ_{i}h_{i} + \frac{1}{3}κ′_{i}h_{i}(h_{i} − h_{i−1}) + O(h^{3}) \\ & = & α_{i}^{∗} + O(h^{2}) \end{aligned}  (10)

as good approximations of α_{i}^{∗}. We therefore recommend (10) for local methods and for starting Newton’s method to solve the system (2). The goal of this section is to prove

Theorem 6.1 Global GC^{2} piecewise quadratic interpolation of data from a smooth planar curve with nonvanishing curvature has an error of O(h^{4}) for data density h → 0.

The proof will be reduced to an application of

Theorem 6.2 If local Hermite interpolation of a smooth planar curve near a point with nonvanishing curvature is done by quadratic polynomials using exact positions and approximate tangent directions with errors O(h^{k}) for k = 2 or 3, the error will be O(h^{k+1}) in terms of the distance h between interpolation points.

The proof of this result can be omitted because it is just a slight variation of the argument given by deBoor, Höllig, and Sabin in [deBoor et. al., ’87] to prove their O(h^{6}) result for GC^{2} cubics.

Convergence of order O(h^{3}) for the convexity preserving local Hermite interpolation by GC^{1} piecewise quadratics using approximate tangent angles α̂_{i} defined in (10) follows immediately from Theorem 6.2. To prove O(h^{4}) convergence for the global GC^{2} method of this paper, we have to show

α̃_{i} = α_{i}^{∗} + O(h^{3}), h → 0  (11)

uniformly with respect to i for the solution α̃_{i} of the system (2). This will be done via Kantorovitch’s convergence theorem for Newton’s method (see e.g. [Ortega–Rheinboldt, ’70, p. 421]) applied to a scaled version of the system (2).

To evaluate (2) at the exact tangent angles α^{∗} we fix an index i and apply lemmas 6.1 and 6.2 to get the expansions

\begin{aligned}α_{i−1}^{∗} & = & \frac{1}{2}κ(t_{i})h_{i−1} − \frac{2}{3}κ′(t_{i})h_{i−1}^{2} + O(h^{3}) \\ α_{i+1}^{∗} & = & \frac{1}{2}κ(t_{i})h_{i} − \frac{2}{3}κ′(t_{i})h_{i} + O(h^{3}) \\ γ_{i+1} & = & \frac{1}{2}κ(t_{i})(h_{i} + h_{i+1}) + \frac{1}{3}κ′(t_{i})(2h_{i}^{2} + 3h_{i}h_{i+1} + h_{i+1}^{2}) + O(h^{3}) \\ γ_{i} − α_{i}^{∗} & = & \frac{1}{2}κ(t_{i})h_{i−1} − \frac{1}{3}κ′(t_{i})h_{i−1}^{2} + O(h^{3}) \\ γ_{i+1} − α_{i+1}^{∗} & = & \frac{1}{2}κ(t_{i})h_{i} + \frac{2}{3}κ′(t_{i})h_{i}^{2} + O(h^{3}) \\ γ_{i} − α_{i}^{∗} + α_{i−1}^{∗} & = & κ(t_{i})h_{i−1} − κ′(t_{i})h_{i−1}^{2} + O(h^{3}) \\ γ_{i+1} − α_{i+1}^{∗} + α_{i}^{∗} & = & κ(t_{i})h_{i} + κ′(t_{i})h_{i}^{2} + O(h^{3}). \end{aligned}

Putting these into (2) yields an O(h^{2}) expression after some calculations. Similar expansions can be derived for the angles α̂_{i} replacing α_{i}^{∗}, and (2) comes out to be of order O(h) at α̂. This is not surprising, because the approximate angles defined in (10) are exact up to O(h^{3}) for circular arcs.

To evaluate derivatives, we rescale (2) from the form F (α) = 0 into G(β) = 0 by setting β := T (α), G(β) := F (T^{−1}β) with T defined by

β_{j} := T_{j}(α) := \frac{α_{j}}{min(h_{j−1}, h_{j})}, 2 ≤ j ≤ n − 1.

Using β̂:= T (α̂) and β^{∗} := T (α^{∗}) we have

F (α̂) = G(β̂) = O(h), F (α^{∗}) = G(β^{∗}) = O(h^{2})

for the scaled version. Now we examine G′ at β^{∗} and get

\frac{∂G_{i}}{∂β_{i−1}} = 2\frac{sin^{2}(γ_{i} − α_{i}) sin(γ_{i} − α_{i} + α_{i−1})}{sin^{3} α_{i−1}} \frac{min(h_{i−1}, h_{i−2})}{h_{i−1}}
= (4 + O(h_{i−1}))\frac{min(h_{i−1}, h_{i−2})}{h_{i−1}} = 4 + O(h),
\frac{∂G_{i}}{∂β_{i+1}} = 2\frac{sin^{2}(α_{i}) sin(γ_{i+1} − α_{i+1} + α_{i})}{sin^{3}(γ_{i+1} − α_{i+1})} \frac{min(h_{i−1}, h_{i})}{h_{i}}
= (4 + O(h_{i}))\frac{min(h_{i−1}, h_{i})}{h_{i}} = 4 + O(h),
\frac{∂G_{i}}{∂β_{i}} = 2\frac{sin(γ_{i} − α_{i} + α_{i−1})}{sin^{2} α_{i−1}} \frac{min(h_{i−1}, h_{i−2})}{h_{i−1}}
·(sin(α_{i−1}) + 3 sin(γ_{i} − α_{i}) cos(γ_{i} − α_{i} + α_{i−1}))
+ 2\frac{sin(γ_{i+1} − α_{i+1} + α_{i})}{sin^{2}(γ_{i+1} − α_{i+1})} \frac{min(h_{i−1}, h_{i})}{h_{i}}
·(sin(γ_{i+1} − α_{i+1}) + 3 sin(α_{i}) cos(γ_{i+1} − α_{i+1} + α_{i}))
= (8 + O(h_{i−1}))\frac{min(h_{i−1}, h_{i−2})}{h_{i−1}} + (8 + O(h_{i}))\frac{min(h_{i−1}, h_{i})}{h_{i}} = 16 + O(h)
= 2\frac{∂G_{i}}{∂β_{i−1}} + 2\frac{∂G_{i}}{∂β_{i+1}} + O(h).

The row–sum norm of G′ at β^{∗} then is bounded by 32 + O(h), while the inverse has a norm bound of \frac{1}{4} + O(h). A similar calculation proves that all second derivatives of G_{i} at β^{∗} are bounded with respect to h → 0. This is where the special form of the scaling pays off.

Kantorovitch’s theorem now implies that there exist r > 0 and h_{0} > 0 such that there is a unique solution β̃ of G(β) = 0 for each interpolation problem for f with data density h less than h_{0}, satisfying ∥β̃− β^{∗}∥ < r in the sup–norm. Furthermore, the solution is reached with quadratic convergence from every starting point β with ∥β − β^{∗}∥ < r. Because of

|β̂_{i} − β_{i}^{∗}| = |α̂_{i} − α_{i}^{∗}|/ min(h_{i}, h_{i−1})
= (\frac{1}{3}|κ′(t_{i})|h_{i}h_{i−1} + O(max(h_{i}^{3}, h_{i−1}^{3}))/ min(h_{i}, h_{i−1})

the starting value β̂ will be successful for h → 0, if the ratio

max(h_{i}, h_{i−1})/ min(h_{i}, h_{i−1})  (12)

is uniformly bounded.

Since G(β^{∗}) = O(h^{2}), the radius of the existence domain can be (theoretically) chosen as an O(h^{2}) function, and therefore the exact solution β̃ satisfies ∥β̃ − β^{∗}∥ = O(h^{2}). But then

\begin{aligned}α̃_{i} & = & β̃_{i}(h_{i} + h_{i−1}) = β_{i}^{∗}(h_{i} + h_{i−1}) + O(h^{3}) \\ & = & α_{i}^{∗} + O(h^{3}), \end{aligned}

as required for O(h^{4}) convergence of the interpolation process.

7 Numerical treatment

The asymptotic considerations of the preceding section suggest a modified Newton–Raphson method for solving the nonlinear system (2). The Jacobians will be tridiagonal, and each iteration needs only O(n) operations. If (2) is written as F (α) = 0 with F : IR^{n−2} → IR^{n−2}, then a simple stepsize control is recommended to make sure that

a) ∥F (α)∥_{2} decreases at each step

b) α stays in the admissible region (e.g. 0 < α_{i} < γ_{i} for “right–turning” data or γ_{i} < α_{i} < 2π for “left–turning” data).

(see [Schaback, ’88] for a comparable situation). Good results are obtained for starting values (10), because these satisfy (2) up to O(h). A Gauss–Seidel iteration of (2) may be employed for an existence proof via Brouwer’s fixed–point theorem, but there are examples where the unique solution of (2) is a repelling fixed point of the Gauss–Seidel iteration.

Newton’s method also works well in cases with multiple solutions and in cases of reflecting data with existing solutions. In cases of nonexisting solutions the method creeps to the boundary of the admissible domain, using stepsizes tending to zero. There are nondegenerate data sets where the matrix F′ is not diagonally dominant at the solution, but no nondegenerate case is known where Newton’s method failed to converge. For sufficiently dense samples of data with density h → 0 and a bounded ratio (12), Newton’s method will converge quadratically to the solution when started at (10), as the convergence results of the preceding section show.

8 Examples and Concluding Remarks

Figures 11 and 12 show two examples.

[figure]

Figure 11: Logarithmic spiral, 8 points

For large numbers of data points, the solutions look very nice and closely approximate functions that supplied the data. If data points are far apart, the solution looks a little “too straight” at interior knots (see Figure 11). This is due to the fact that the interpolating parabolic pieces normally attain their curvature maximum in the interior of their domain.

Existence and uniqueness results extend to the case of piecewise rational quadratic curves with prescribed weights. This modifies the system (2) by some extra factors, but leaves the geometric proof for existence and uniqueness of solutions unchanged.

In a forthcoming paper we show how to place cubic polynomial pieces at positions where the data require an inflection point or a GC^{2} continuation into a straight line. Furthermore, the tangential boundary conditions can be replaced by “not–a–knot–” conditions.

References

[1] Boehm, W., Farin, G.F., Kahmann, J. (1984), A survey of curve and surface methods in CAGD, Computer Aided Geometric Design 1, 1 – 60

[2] deBoor, C., Höllig, K., and Sabin, M. (1987), High Accuracy Geometric Hermite Interpolation, Computer Aided Geometric Design 4, 269–278

[3] Farin, G.E. (1982), Visually C^{2} cubic splines, Computer Aided Design 14, 137 – 139

[4] Ortega, J.M., and W.C. Rheinboldt (1970), Iterative Solution of Nonlinear Equations in Several Variables, Academic Press [5] Schaback. R. (1988), Adaptive Rational Splines, NAM–Bericht 60, to appear in Con-

structive Approximation 1989.

Author’s address:

Prof. Dr. R. Schaback Institut für Numerische und Angewandte Mathematik Lotzestraße 16–18

D–3400 Göttingen

[figure]

Figure 12: Letter “O” formed by 15 points