Remarks on High Accuracy Geometric Hermite Interpolation

Robert Schaback

Abstract. In a recent paper deBoor, Höllig, and Sabin proposed a method for local sixth– order interpolation of convex planar curves by geometrically C^{2} piecewise cubic curves, using positional, directional, and curvature data at breakpoints. This note extends their results to curves with simple zeros of curvature and gives some hints to stabilize the numerical construction.

Keywords. Curves, Interpolation, Geometric Continuity, Convergence Order.

1 Hermite Interpolation

This paper follows the work [1] of deBoor, Höllig, and Sabin concerning Hermite interpolation of smooth curves f : IR → IR^{2}. Their piecewise cubic geometrically C^{2} interpolants use positional, directional, and curvature data. They are approximations of order O(h^{6}) with respect to the maximum distance h of two adjoining interpolation points, provided that the curvature κ_{f} of f does not vanish. The approximation order can drop to four in presence of double zeros of curvature, as is shown in [1] by an example.

However, numerical experience strongly suggests that the approximation order still is six when simple zeros of curvature are allowed, provided that the numerical solution of the problem is carried out with some care. More precisely,

Theorem 1.1 Let f be a smooth planar curve, parametrized by arclength around a point f(0) where the curvature κ_{f} of f has a simple zero. Then there is a neighborhood of zero in which Hermite interpolation with piecewise cubics in the sense of deBoor, Höllig, and Sabin [1] has a solution with approximation order six for sufficiently dense data.

Proof: Following [1] the angle θ(h) between tangents at f(0) and f(h) can be written as

θ(h) = ∑_{i=1}^{4} θ_{i}h^{i} + O(h^{5}).  (1)

Then the curvature of f at h is κ_{f}(h) = θ′(h), and κ_{f}(0) = θ_{1}, κ′_{f}(0) = 2θ_{2}. Now assume the data f(−h), f(h), f′(−h), f′(h), κ_{f}(−h), κ_{f}(h) to be given. A cubic polynomial in Bernstein–Bézier form is constructed, whose control points are f(−h), f(−h) + δ_{0}(h)f′(−h), f(h) − δ_{1}(h)f′(h), f(h), and where positive scalars δ_{0}(h) and δ_{1}(h) have to be found such that the given curvature values κ_{f}(−h) and κ_{f}(h) are attained. With the same notation as in [1], but explicitly stating the symmetry with respect to the substitution h → −h, the data

1File : usr/nam/rschaba/tex/fertig/dbhs.tex, TEXed May 18, 1994, Status : Not submitted

\begin{aligned}a(h) & := & f(h) − f(−h) = −a(−h) \\ d_{0}(h) & := & f′(−h) = d_{1}(−h) & d_{1}(h) & := & f′(h) = d_{0}(−h) \\ c(h) & := & d_{0}(h) × d_{1}(h) = −c(−h) \\ c_{0}(h) & := & d_{0}(h) × a(h) = c_{1}(−h) & c_{1}(h) & := & a(h) × d_{1}(h) = c_{0}(−h) \\ κ_{0}(h) & := & κ_{f}(−h) = κ_{1}(−h) & κ_{1}(h) & := & κ_{f}(h) = κ_{0}(−h), \end{aligned}  (2)

are calculated, where the cross–product is defined by (u_{x}, u_{y}) × (v_{x}, v_{y}) := u_{x}v_{y} − u_{y}v_{x}. Then the unknowns δ_{0} and δ_{1} have to satisfy the system

\begin{aligned}0 & = c(h)δ_{0}(h) − c_{1}(h) + \frac{3}{2}κ_{1}(h)δ_{1}(h)^{2} \\ 0 & = c(h)δ_{1}(h) − c_{0}(h) + \frac{3}{2}κ_{0}(h)δ_{0}(h)^{2}, \end{aligned}  (3)

which may have up to four solutions in general. The paper [1] of deBoor, Höllig, and Sabin handles all positive solutions, while here a specific solution branch as a function of h is selected.

Because of the symmetry of the problem all solutions should have rotational symmetry

δ_{0}(−h) = −δ_{1}(h)  (4)

around zero. Thus there is a single functional equation of the form

F (δ_{1}, h) := −c(h)δ_{1}(−h) − c_{1}(h) + \frac{3}{2}κ_{1}(h)δ_{1}(h)^{2} = 0  (5)

for the function δ_{1}(h) around h = 0. To investigate this equation, a good approximate solution of the special form

δ̂_{1}(h) = \frac{2}{3}h + d_{2}h^{2} + d_{3}h^{3}  (6)

is constructed to satisfy (5) up to terms of order six in h. The expansion (6) is motivated by a result of [1]: there are up to three positive solutions of (3) of the form δ_{i}(h) = 2h/3 + (−1)^{i}ϵ(h)h^{2} + O(h^{3}) with bounded ϵ(h) if κ_{f} does not vanish.

Using REDUCE (up to terms of order ten, for safety) equation (6) is put into (5) to get

F (δ̂_{1}, h) = h^{4}c_{4} + h^{5}c_{5} + O(h^{6}),

leading to

\begin{aligned}c_{4} := & −(4 θ_{1}^{3} − 27 θ_{1}d_{2}^{2} − 72 θ_{1}d_{3} − 72 θ_{2}d_{2} − 24 θ_{3})/18 = 0 \\ c_{5} := & (20 θ_{1}^{3}d_{2} + 12 θ_{1}^{2} θ_{2} + 45 θ_{1}d_{2}d_{3}+ \\ & 45 θ_{2}d_{2}^{2} + 60 θ_{2}d_{3} + 60 θ_{3}d_{2} + 16 θ_{4})/15 = 0. \end{aligned}  (7)

These equations are now solved for d_{2}, d_{3}. In case θ_{1} ≠ 0 it is convenient to solve c_{4} = 0 for d_{3} and c_{5} = 0 for d_{2}. However, in the situation of this paper it is reasonable to solve first for d_{3} from c_{5} = 0, because in the case θ_{1} = 0 ≠ θ_{2} the equation c_{4} = 0 completely determines d_{2}. This leads to

\begin{aligned}d_{3} := & −(20 θ_{1}^{3}d_{2} + 12 θ_{1}^{2} θ_{2} + 45 θ_{2}d_{2}^{2} + 60 θ_{3}d_{2} + 16 θ_{4})(15(3 θ_{1}d_{2} + 4 θ_{2}))^{−1} \\ c_{4} = & −(540 θ_{1}^{4}d_{2} + 368 θ_{1}^{3} θ_{2} − 405 θ_{1}^{2}d_{2}^{3} − 540 θ_{1} θ_{2}d_{2}^{2} + 1080 θ_{1} θ_{3}d_{2}+ \\ & 384 θ_{1} θ_{4} − 1440 θ_{2}^{2}d_{2} − 480 θ_{2} θ_{3})(90(3 θ_{1}d_{2} + 4 θ_{2}))^{−1} \end{aligned}  (8)

with denominators that will not vanish when θ_{1} ≈ 0 ≠ θ_{2}, and provided that there finally is a bounded solution d_{2} of c_{4} = 0. The numerator of c_{4} now is a cubic polynomial in d_{2}, and essentially the same polynomial results if c_{4} = 0 is solved for d_{3} first, inserting the result into c_{5} = 0. Of this cubic equation REDUCE gives three roots r_{1}, r_{2}, and r_{3}, where r_{3} does not explicitly contain i = \sqrt{−1}. The roots are complicated fractions containing square roots and involving θ_{1},…, θ_{4}. A simplified description is

z_{1} := −2025θ_{1}^{12} − 12150θ_{1}^{9}θ_{3} + 14076θ_{1}^{8}θ_{2}^{2} + 3456θ_{1}^{6}θ_{2}θ_{4} − 24300θ_{1}^{6}θ_{3}^{2}
+45360θ_{1}^{5}θ_{2}^{2}θ_{3} − 24880θ_{1}^{4}θ_{2}^{4} + 5184θ_{1}^{4}θ_{4}^{2} − 25920θ_{1}^{3}θ_{2}θ_{3}θ_{4} − 16200θ_{1}^{3}θ_{3}^{3}
+15360θ_{1}^{2}θ_{2}^{3}θ_{4} + 86400θ_{1}^{2}θ_{2}^{2}θ_{3}^{2} − 98400θ_{1}θ_{2}^{4}θ_{3} + 33600θ_{2}^{6}
z_{2} := 3 · \sqrt{z_{1}} + 72θ_{1}^{4}θ_{2} + 216θ_{1}^{2}θ_{4} − 540θ_{1}θ_{2}θ_{3} + 320θ_{2}^{3}
z_{3} := z_{2}^{1/3}
z_{4} := z_{3}^{2}
z_{5} := 5^{2/3}(θ_{1}^{4} + 18 · θ_{1}θ_{3} − 20 · θ_{2}^{2})
z_{6} := 9 · 5^{1/3}z_{3}θ_{1}
z_{7} := 4 · 5^{1/3}z_{3}θ_{2}
r_{1} := (−z_{4}(1 + i\sqrt{3}) + z_{5}(1 − i\sqrt{3}) − z_{7})/z_{6}
r_{2} := (−z_{4}(1 − i\sqrt{3}) + z_{5}(1 + i\sqrt{3}) − z_{7})/z_{6}
r_{3} := (2z_{4} − 2z_{5} − z_{7})/z_{6}.

Now the behaviour of r_{3} in case θ_{1} ≈ 0 ≠ θ_{2} is studied. The denominator of r_{3} has the form

18θ_{1}θ_{2}(5(15\sqrt{21} + 40))^{1/3} + O(θ_{1}^{2})

for θ_{1} → 0, and the numerator has the value

8 θ_{2}^{2}((15\sqrt{21} + 40)^{2/3} − (15\sqrt{21} + 40)^{1/3}5^{1/3} − 5^{5/3})

at θ_{1} = 0, as calculated by REDUCE. This result is disappointing, because it suggests a singularity of r_{3} at θ_{1} = 0. However, the complicated numerical factor in the numerator can be shown to be zero, since the solution x = (5^{1/3}(1 + \sqrt{21}))/2 of x^{2} − 5^{1/3}x = 5^{5/3} satisfies x^{3} = 15\sqrt{21} + 40.

Taking derivatives with REDUCE, the rule of de l’Hospital yields

r_{3} = −\frac{1}{3} \frac{θ_{3}}{θ_{2}} + O(θ_{1})

for θ_{1} → 0. If the root r_{3} is chosen to define d_{2} and if d_{3} is calculated from (8), there is a bounded solution of c_{4} = c_{5} = 0 with respect to θ_{1} → 0, if θ_{2} ≠ 0.

The construction up to here yields F (δ̂_{1}, h) = O(h^{6}) for small h. To compare this approximate solution with an exact solution, the system (3) is rescaled as

\begin{aligned}0 & = c(h)ϵ_{0}(h)h^{−1} − c_{1}(h)h^{−2} + \frac{3}{2}κ_{1}(h)ϵ_{1}(h)^{2} \\ 0 & = c(h)ϵ_{1}(h)h^{−1} − c_{0}(h)h^{−2} + \frac{3}{2}κ_{0}(h)ϵ_{0}(h)^{2}. \end{aligned}  (9)

in new unknowns ϵ_{i}(h) := δ_{i}(h)/h for i = 0, 1. If this system is written as

T_{h}(ϵ_{0}(h), ϵ_{1}(h)) = 0, T_{h} : IR^{2} → IR^{2},

the approximate solution satisfies

T_{h}(−δ̂_{1}(−h)/h, δ̂_{1}(h)/h) = O(h^{4}).

The derivative T′_{h} of T_{h} is uniformly Lipschitz continous for h → 0. Using REDUCE again, one gets the determinant as

det T′_{h} = 16h^{2}(θ_{2}^{2} + O(θ_{1})) + O(h^{4})

and ∥(T′_{h})^{−1}∥ ≤ O(h^{−2}) for θ_{1} ≈ 0 ≠ θ_{2} and h → 0.

Now a Newton–Kantorovitch–type theorem (see e.g. [3], p. 421) yields the existence of an exact solution ϵ_{0}(h), ϵ_{1}(h) of the system (9) satisfying the inequality

|ϵ_{i}(h) − δ̂_{i}(h)/h| ≤ O(h^{2}), i = 0, 1; δ̂_{0}(h) := −δ̂_{1}(−h),

when h is sufficiently small. This means there is a solution of the form

δ_{i}(h) = δ̂_{i}(h) + O(h^{3}), i = 0, 1

to the system (3). The rest follows as in [1].

Remark : The proof in [1] for the case θ_{1} ≠ 0 allowed all of the three possible solution branches of type (6) of the system (3) to be chosen. These solution branches of (3) are approximated by the solution branches of (5) for the approximation (6), and there is a useless fourth branch of (3) with the behaviour δ̂_{1}(h) = −2h+d_{2}h^{2} +d_{3}h^{3}. In the general case one has to make sure that the solution branch defined by r_{3} is taken whenever θ_{1} becomes small with respect to θ_{2}. Therefore the overall approximation order of the cubic Hermite interpolation method is six, provided that there are no double zeros of curvature and that the correct solution branch near inflection points is chosen. The next section will give an easy recipe to choose the proper solution branch.

2 Choosing the Solution Branch

Since the numerical treatment of the method of deBoor, Höllig, and Sabin [1] appears to be rather hazardous, some hints for the implementation seem to be necessary.

The coefficients of the system (3) are explicitly calculated from given data via (2), and the chordlength d = ∥f(h) − f(−h)∥_{2} = ∥a(h)∥_{2} can be used as an approximation of the arclength h. Then the system (3) should have solutions δ_{0}, δ_{1} near 2h/3.

Theorem 2.1 To get O(h^{6}) convergence for interpolation of any smooth curve without double zeros of curvature it suffices to choose any positive solution (δ_{0}, δ_{1}) of (3) which is not too far away from (2h/3, 2h/3).

Proof : Since (3) is equivalent to a fourth–order algebraic equation, there are at most four complex–valued solution branches. Constructing expansions with respect to h with REDUCE, one can find a single branch with the behaviour (−2h, −2h) + O(h^{2}). This branch yields (useless) negative values of δ_{0}, δ_{1} and is automatically discarded, because there are up to three branches of the type (2h/3, 2h/3)+O(h^{2}), which were already found in [1] under the assumption of nonvanishing curvature. Each of the branches is useful in case of data with curvature bounded away from zero, and it then does not matter which branch is chosen.

In case θ_{1} ≈ 0 ≠ θ_{2} the construction of the first section can be carried out. This gives three solution branches (two of which may be complex) of the system (7), and the Newton–Kantorovitch argument shows that for each of the three approximate solutions of type (6) there exists a neighbouring solution of (3) with the same local behaviour. Thus all interesting solutions of (3) are covered by this approach, but REDUCE yields two complex solution branches of (3) via the solutions

r_{1} = (0.5103 + i · 1.728)\frac{θ_{3}}{θ_{2}} + O(θ_{1})

r_{2} = (0.5103 − i · 1.728)\frac{θ_{3}}{θ_{2}} + O(θ_{1})

of (7) for θ_{1} → 0 with θ_{2} ≠ 0. Thus there is only one solution of the required type, if θ_{1} is small and θ_{2}, θ_{3} are bounded away from zero. This was already observed in [1] for cases without sign change of curvature, if curvature is sufficiently small. If in the situation of this paper also θ_{3} is small, the three roots tend to coalesce stably for θ_{3} → 0, and in the limit θ_{3} = 0 all roots satisfy d_{2} = 0, d_{3} = −\frac{4}{15}θ_{4}/θ_{2}. The suggested rule for the choice of solutions of (3) will thus work in all cases except when curvature and its derivative vanish simultaneously.

In the latter case the equations (7) can be satisfied only if θ_{3} = θ_{4} = 0, i.e.: the zeros of curvature of the given curve must be either simple or at least of order four to make the approach of section 1 possible. It is conjectured that (at least generically) the local convergence order of the method of deBoor, Höllig, and Sabin is six if and only if the zeros of curvature of the interpolated curve are either simple or of order at least four.

3 Solving the Local System

Since the system (3) normally has multiple solutions which may coalesce, standard approaches like Newton’s method or direct evaluation of the roots of the fourth–degree polynomial behind (3) are not stably applicable. As a compromise between efficiency, reliability, and simplicity it is suggested to write (3) as a single equation φ(x) = 0 which is solved by a straightforward bisection technique T that starts from a given point x_{0} and finds the nearest zero x = T (x_{0}, φ) of φ. This can be done with linear convergence. In a large number of cases the constants

R_{0} = \frac{3}{2} \frac{κ_{1}}{c_{1}} (\frac{c_{0}}{c} )^{2} , R_{1} = \frac{3}{2} \frac{κ_{0}}{c_{0}} (\frac{c_{1}}{c} )^{2} ,

can be calculated and lie in some interval [0, γ] ⊂ IR with a fixed γ ≥ 1. Then the method T is applied to the function φ(x) := x − 1 + R_{1}(1 − R_{0}x^{2})^{2} with respect to the variable x = cδ_{0}/c_{1}, where 0 and 1 are possibly swapped to ensure |R_{0}| ≥ |R_{1}|. This works in most cases where curvature is bounded away from zero, because R_{i} = \frac{3}{4} + O(h^{2}) for h → 0 and θ_{1} ≠ 0.

If curvature tends to zero, a term of type θ_{2}^{2}θ_{3}^{−2}h^{−1} in the R_{i} for θ_{1} = 0 spoils this approach, and R_{i}(h) goes to infinity for h → 0 in case θ_{1} = 0. Since data from smooth curves satisfy κ_{1}c_{1} = 2θ_{1}^{2}h^{2} + O(h^{3}) for θ_{1} ≠ 0 or κ_{1}c_{1} = \frac{8}{3}θ_{2}^{2}h^{4} + O(h^{5}) for θ_{1} = 0, the R_{i} will be asymptotically positive. Furthermore, R_{0}(h)/R_{1}(h) always tends to 1 for h → 0. Thus degeneration of φ will occur only in the form of R_{i} tending simultaneously to +∞, if there are no multiple zeros of curvature.

But then a different choice of φ automatically becomes feasible and stable:

φ(x) := x^{2} − S_{0} + T_{0}\sqrt{S_{1}x − T_{1}},

where x = δ_{0}/κ_{0} and

S_{0} = \frac{2}{3} \frac{c_{0}}{κ_{0}^{3}} , S_{1} = \frac{2}{3} \frac{c_{1}}{κ_{1}^{3}} , T_{0} = \frac{2}{3} \frac{cκ_{1}}{κ_{0}^{3}} , T_{1} = \frac{2}{3} \frac{cκ_{0}}{κ_{1}^{3}} .

Here S_{i} = (6θ_{2}^{2})^{−1} + O(h) and T_{i} = O(h) for h → 0 in case of θ_{1} = 0. This function does not degenerate for curvature tending to zero, if there is no double zero of curvature.

If both methods are not applicable or give negative results, the starting value δ_{0} = δ_{1} = 2h/3 (which still yields fourth order of convergence and fairly good pictures) is kept and a warning message is given. This can only happen for large h or for data from functions that are either nonsmooth or have multiple zeros of curvature.

The bisection method will always use δ_{0} = 2h/3 as a starting point and therefore it will automatically pick a solution δ_{0} of (3) near to 2h/3. The other component δ_{1} will also be of this size in case of dense data from a smooth curve without multiple zeros of curvature, because then all solution branches satisfy the symmetry property (4) for small h. Thus the solution should simply be rejected if it is in some sense too far away from (2h/3, 2h/3).

Explicit numerical examples are suppressed here, since they will not lead much beyond the examples of [1] and will always reproduce smooth curves with graphically invisible errors.

The paper [4] contains a method to estimate curvature data to high accuracy. This can be applied to make the method of deBoor, Höllig, and Sabin work for positional and tangent data only. High–accuracy Hermite interpolation of nonplanar curves by piecewise cubic rationals is treated in [2].

References

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

[2] Höllig, K., Algorithms for Rational Spline Curves, Technical Report # 710, Comp. Science Dept., Univ. of Wisconsin, 1987

[3] Ortega, J.M., and W.C. Rheinboldt, Iterative Solution of Nonlinear Equations in Several Variables, Academic Press 1970

[4] Schaback, R., Convergence of Planar Curve Interpolation Schemes, in: Approximation Theory VI, Volume 2, of C.K. Chui, L.L. Schumaker, and J.D. Ward (eds.), Academic Press, New York, 1989, 581–584