SPHERICAL DISPERSION WITH AN APPLICATION This paper appears in Anzeiger der Österreichischen Akademie der Wissenschaften, Mathematisch-naturwissenschaftliche Klasse, Abteilung II 132 (1995), 3–10.

TO POLYGONAL APPROXIMATION OF CURVES

Günter Rote and Robert F. Tichy

Abstract. Sphärische Dispersion mit einer Anwendung auf die polygonale Approximation von Kurven. Wir beschreiben ein Verfahren zur polygonalen Approximation von Raumkurven unter Verwendung ”gut verteilter“ Punkte auf der Kugeloberfläche und analysieren es. Dabei spielt die Dispersion der Punktmenge bezüglich ”Kugelspalten“ (Durchschnitten zweier Halbsphären) eine Rolle.

In this note we consider finite point sets A on the n-dimensional unit sphere S^{n}. The dispersion of A with respect to a given family R of subsets of S^{n}, called ranges, is defined by

d_{R}(A) := sup { μ(R) : R ∈ R, A ∩ R = ∅ } ,

where μ is the normalized surface measure on S^{n}. This notion extends the concept of dispersion on the unit interval which was introduced by Hlawka [Hl 76] and later investigated in more general form in Niederreiter [Ni 81]. In that setting, the dispersion is necessary for analyzing optimization algorithms based on uniformly distributed sequences, see Niederreiter [Ni 83], cf. also [NP 86] and chapter 6 in the monograph [Ni 92]. In [Ti 90], optimization problems on the sphere are considered. As a general reference on uniformly distributed sequences, see for example Kuipers and Niederreiter [KN 74] and Hlawka [Hl 79]. For applications of point sequences to numerical integration, we refer to Sobol [So 67].

The range space that is usually considered on the sphere is the range space C of spherical caps, i. e., intersections of the sphere with half-spaces. Recently, Blümlinger [Bl 91] proved a strong lower bound for the discrepancy with respect to the family S of spherical slices. A slice is the intersection of two half-spheres.

In the following we will consider the dispersion with respect to S, and we will explain an application of spherical slice dispersions to the piecewise linear approximation of curves in space, improving an algorithm of Schachinger [Sch 90] for such approximations.

Obviously, for any point set A on the sphere S^{n}, the following elementary relation holds between cap and slice dispersion:

(1)  d_{S}(A) ≤ c_{n} · d_{C}(A)^{1/n},

for some positive constant c_{n}. This is obvious from the fact that any slice contains a cap of corresponding volume.

Typeset by AMS-TEX

The second author was supported by the Austrian Science Foundation (Project P8274PHY)

We also have the trivial lower bound d_{S}(A) ≥ 1/N for every N-point set A. The following proposition states that this bound can be achieved, up to a constant factor.

Proposition. For every N there is a point set A with N points on the sphere S^{n} which has slice dispersion

d_{S}(A) = O(1/N).

Proof. W. l. o. g. we assume that S^{n} is the unit sphere in (n + 1)-dimensional space. We evenly distribute the N points over the \binom{n+1}{2} two-dimensional coordinate planes (on the sphere). On each of the “coordinate circles” C_{i}, i = 1, …, \binom{n+1}{2}, which are the intersections of the coordinate planes with the sphere, we place the corresponding points equidistantly.

Now let Z be an empty slice, i. e., a slice containing no point of A. Let U be the two-dimensional plane through the origin which is orthogonal to the two hyperplanes bounding Z. Let b_{i} be the projection factors for projections between U and the i-th coordinate plane, in the sense that a set of area a in the plane U is mapped to a set of area b_{i}a by an orthogonal projection onto the respective coordinate plane. By the generalized Pythagorean theorem (cf. [Ga 58], section IX.5, p. 223),

\begin{matrix}( & {n+1 \atop 2} & ) \\ ∑ \\ i=1 \end{matrix} b_{i}^{2} = 1.

Thus there is a coordinate plane whose projection factor b_{i} is at least \sqrt{1/\binom{n+1}{2}}.

The projection of the circle C_{i} onto the plane U is an ellipse with area πb_{i}, whose major axis is at most 1 and whose minor axis is therefore at least b_{i}. C_{i} contains at least ⌊N/\binom{n+1}{2}⌋ evenly spaced points. In order to see how these points restrict the opening angle of Z, let us imagine that a point moves on C_{i} at constant speed ω. The projection of this point onto U moves on the ellipse with speed at most ω, and since its distance from the origin is at least b_{i} its angular speed around the origin is at most ω/b_{i}. Since the angle between two adjacent points on C_{i} is at most 2π / ⌊N/\binom{n+1}{2}⌋, it follows that the opening angle of Z can be at most

2π / ⌊N/\binom{n+1}{2}⌋ / b_{i}, and its measure is μ(Z) ≤ 1 / ⌊N/\binom{n+1}{2}⌋ / b_{i} ≤ \sqrt{\binom{n+1}{2}} /

⌊N/\binom{n+1}{2}⌋ = O(n^{3}/N).

Remark. The above example shows that the converse of (1) does not hold since its cap dispersion is Ω(1).

In the remainder of this paper we will show how the slice dispersion on the sphere in three-dimensional space arises in a problem of piecewise linear approximation of curves in space.

For instance, in robotics it is an important problem to approximate a “general” curve by simple curves like straight lines, circles etc., because the arm of the robot can only run along such simple curves. The most important case is the approximation by a polygonal line.

Let us first consider a twice continuously differentiable spatial curve x⃗(s) parameterized by its arc length s. We want to construct a sequence of points x⃗(s_{0}), x⃗(s_{1}), …, x⃗(s_{M}) such that the polygonal line with vertices at these points is an ε-approximation of the curve segment C: { x⃗(s) : s_{0} ≤ s ≤ s_{M} }, where ε is an arbitrarily given (small) positive number. Let us recall here the definition of an εapproximation: A (closed) point set A is an ε-approximation of the (closed) set B if their Hausdorff-distance

δ(A, B) = max (max_{x∈A} d(x, B), max_{y∈B} d(y, A))

is not greater than ε (d denoting the Euclidean distance). We consider a line segment L_{k} with endpoints x⃗(s_{k}), x⃗(s_{k+1}), and denote the tangent vector by ˙ x⃗(s_{k}). Then δ(C, L_{k}) < ε is guaranteed provided that

∥x⃗(s_{k+1}) − x⃗(s_{k}) − ˙ x⃗(s_{k})(s_{k+1} − s_{k})∥ ≤ ε.

Applying Taylor’s formula yields

∥¨x⃗(s_{k})∥(s_{k+1} − s_{k})^{2} ≤ \frac{ε}{2}

as an approximate sufficient condition that the polygonal line is an ε-approximation. Thus we obtain the following iteration procedure for computing the vertices of the polygonal approximation:

s_{k+1} − s_{k} = \sqrt{\frac{ε}{2κ(s_{k})}},

where κ(s_{k}) denotes the curvature in the point x⃗(s_{k}), see McClure and Vitale [MV 75] and Müller [Mü 92]. This method, of course, has one disadvantage: one has to know arc length and curvature in advance. In the following we describe a different method which makes use of low-dispersed spherical point sequences.

Let us first consider a plane curve C: { x⃗(s) : 0 ≤ s ≤ σ }, which we want to approximate piecewise linearly. We assume that C is parameterized by the arc length s. We successively construct the vertices x⃗(s_{k}) for 0 = s_{0} < s_{1} < ⋯ < s_{M} = σ as follows:

Suppose that s_{k} has already been constructed. Let u_{k} be the largest value (s_{k} < u_{k} ≤ σ) for which there is a direction (i. e., a unit vector) w⃗_{k} such that for any s ∈ [s_{k}, u_{k}] there is a scalar λ with ∥x⃗(s) − x⃗(s_{k}) − λw⃗_{k}∥ ≤ ε. Then we set

(2)  s_{k+1} = max { s ≤ u_{k} : x⃗(s) = x⃗(s_{k}) + λw⃗_{k} for some λ }.

In other words, the curve C between s_{k} and s_{k+1} is contained in the infinitely long strip of width 2ε centered at the line through x⃗(s_{k}) and x⃗(s_{k+1}), and s_{k+1} is the largest possible value with this property. The maximum values u_{k} and the interpolation points x⃗(s_{k}) can be computed easily: We discretize the curve sufficiently fine and maintain the convex hull as we advance on the curve, see Imai and Iri [II 87]. Using the convex hull, it is straightforward to determine whether the curve still fits inside a strip of width 2ε centered at x⃗(s_{k}). Once this is no longer the case, we have found u_{k}. The maximum in (2) can be computed from the discretisation.

However, since the curve is only guaranteed to be contained in the infinitely long ε-strip, there are “pathological” examples where this procedure does not necessarily lead to an ε-approximation. For “reasonable” curves such as may be expected to arise in practice the procedure works. For example, requiring that the curve is smooth and has an upper bound less than 1/ε on the curvature suffices to ensure that the algorithm yields an ε-approximation. This would however exclude polygonal curves, which are especially important in practical applications, where one wants to approximate one such curve by another polygonal curve with fewer vertices.

What we need is a condition on the local “growth rate” of the arc length,

(3) d(x⃗(s), x⃗(s + Δs)) ≥ ρ · Δs, for all Δs ≤ λ and 0 ≤ s < s + Δs ≤ σ.

Here ρ ≤ 1 is the parameter determining the growth rate and the parameter λ makes the condition local. We will say that a curve satisfying (3) has λ-local minimum growth rate at least ρ.

Note that this local condition does not prevent the curve from crossing itself after making a “big” loop. If the curve should have cusps it must be subdivided at these points before applying the algorithm.

The following lemma states that this condition is sufficient for the correctness of our algorithm. Since we will need the spatial case we already formulate the lemma in arbitrary dimensions.

Lemma. Let C be a (continuous) curve from A to B whose (2ε/ρ)-local minimum growth rate is bigger than ρ, for some ρ > 0. If C is contained in the cylinder with axis AB and radius ε then the segment AB is an (ε\sqrt{1 + 1/ρ^{2}})-approximation of C.

Proof. Let us assume w. l. o. g. that A = x⃗(0) is the origin and B = x⃗(σ) lies on the positive x-axis. It is sufficient to show that the x-coordinate never goes below −ε/ρ and never exceeds the x-coordinate of B by more than ε/ρ. By symmetry, we just have to prove the first statement. Let x⃗(s_{0}) with s_{0} ∈ [0, σ] be a point on the curve with negative x-coordinate. Consider the two nearest points x⃗(s_{1}) and x⃗(s_{2}) on the curve with x-coordinates equal to 0: 0 ≤ s_{1} < s_{0} < s_{2} ≤ σ, x(s_{1}) = x(s_{2}) = 0, and x(s) < 0 for s_{1} < s < s_{2}. We show that s_{2} − s_{1} ≤ 2ε/ρ, from which x(s_{0}) ≥ −ε/ρ follows. Otherwise, assume that s_{2} − s_{1} > 2ε/ρ =: Δs. By continuity, we can find two intermediate points x⃗(s′) and x⃗(s′ + Δs) with s_{1} ≤ s′ < s′ + Δs ≤ s_{2} which have the same x-coordinate: If we let s′ vary from s_{1} to s_{2} − Δs we initially have 0 = x(s_{1}) > x(s_{1} +Δs), and at the end we have x(s_{2} −Δs) < x((s_{2} −Δ)+Δs) = 0; thus, there must be a crossover point s′. By (3), d(x⃗(s′), x⃗(s′ + Δs)) ≥ ρ′ · Δs = ρ′(2ε/ρ) > 2ε, where ρ′ > ρ is the local minimum growth rate. But then x⃗(s′) and x⃗(s′ +Δ) cannot both be contained in the cylinder with radius ε around the x-axis, a contradiction.

Remark. One way to ensure (3) is to require that the curve has locally increasing chords:

(4)  d(x⃗(s_{1}), x⃗(s_{4})) ≥ d(x⃗(s_{2}), x⃗(s_{3})),

for all 0 ≤ s_{1} ≤ s_{2} ≤ s_{3} ≤ s_{4} ≤ σ and s_{4} − s_{1} ≤ λ. The minimum growth rate of a plane curve with increasing chords is at least 3/(2π), as is shown in Rote [Ro 94], extending ideas of Larman and McMullen [LM 73]. This relation directly translates to the λ-local concepts. In three dimensions, the minimum growth rate of a curve with increasing chords is at least 0.034, and there is a positive lower bound for any dimension.

A stronger condition than (4) is to require that the curve locally has no angles sharper than π/2, i. e., for any three consecutive points x⃗(s_{1}), x⃗(s_{2}), and x⃗(s_{3}) with 0 ≤ s_{1} ≤ s_{2} ≤ s_{3} ≤ σ and s_{3} − s_{1} ≤ λ, the angle at x⃗(s_{2}) in the triangle x⃗(s_{1})x⃗(s_{2})x⃗(s_{3}) is at least π/2. Clearly, this implies the local increasing chords property.

A direct generalization of the above procedure to the three-dimensional case is not possible for computational reasons. Therefore Schachinger [Sch 90] used projections to reduce the three-dimensional case to the two-dimensional situation.

Let x^{i}(s), (i = 1, …, M) be orthogonal projections P^{i}(x⃗) of the spatial curve C: x⃗(s), 0 ≤ s ≤ σ, onto M suitably chosen planes E_{1}, …, E_{M}. We compute the interpolation points x⃗(s_{k}), k = 1, …, recursively and suppose that s_{0} = 0, s_{1}, …, s_{k} are known. For each projection i = 1, …, M, let S^{i} be the set of parameter values s′ > s_{k} such that there exists a scalar λ between 0 and 1 satisfying ∥x^{i}(s)−x^{i}(s_{k})− λ(x^{i}(s′) − x^{i}(s_{k}))∥ ≤ ε for any s ∈ (s_{k}, s′]. S^{i} is the set of possible parameter values for the next interpolation point s_{k+1} when seen in the i-th projection direction. Each set S^{i} is a finite union of intervals, and we define s_{k+1} as the largest value in the intersection of these M sets.

Now we connect the points x⃗(s_{k}) by line segments. In each projection i, the curve between x^{i}(s_{k}) and x^{i}(s_{k+1}) lies in the infinite 2ε-strip centered at the line through x^{i}(s_{k}) and x^{i}(s_{k+1}). Note that it is not sufficient to determine just the maximum possible parameter value s_{k+1} by (2) in each projection and take the minimum of these values, because this value might not be contained in each S^{i}.

Although the 2ε-strips cover the curve in each projection, the spatial curve between x⃗(s_{k}) and x⃗(s_{k+1}) does not necessarily lie in the infinite cylinder with radius ε centered at the line through x⃗(s_{k}) and x⃗(s_{k+1}).

We will now discuss how to choose the M projection planes to ensure that this holds for a cylinder with radius Aε, for a constant A > 1 which we want to be as small as possible.

Denoting the direction of the projection P^{i} by p_{i} we set for any straight line g

Z(p_{1}, …, p_{M}; g) = ⋂_{i=1}^{M} { x : δ(P^{i}(x), P^{i}(g)) ≤ 1 }.

This set is an intersection of M parallel slabs. The intersection of Z(p_{1}, …, p_{M}; g) with an orthogonal plane of g is a convex symmetric polygon Z, whose edges are parallel to the projections of p_{i} onto Z. The distance of the edges to the center of Z is 1. Setting A equal to the maximal distance of a vertex of Z to the center, the constant A fulfills the desired property if we take the line through the points x⃗(s_{k}) and x⃗(s_{k+1}) as g. Clearly,

A = \frac{1}{cos α/2},

where α is the maximal angle between two adjacent edges. If we consider the projection directions p_{1}, …, p_{M} and the line g as points on the sphere S^{2}, α is the opening angle of the largest empty slice with corners at the two points corresponding to g. Since we want α to be small for all directions g, we have to choose these M points exactly in such a way that the slice-dispersion is minimal. This can be achieved by taking the point set in the proposition for the two-dimensional sphere S^{2}.

Assuming a local smoothness property like (3), we can conclude by the lemma that the polygonal curve is an O(ε)-approximation of C. Acknowledgement. We are indebted to M. Blümlinger and W. Schachinger for valuable discussions and comments.

References

[Bl 91] M. Blümlinger, Slice discrepancy and irregularities of distribution on spheres, Mathe-

matika 38 (1991), 105–116.

[Ga 58] F. R. Gantmacher, Matrizenrechnung I, VEB Deutscher Verlag d. Wiss., Berlin, 1958. [Hl 76] E. Hlawka, Abschätzung von trigonometrischen Summen mittels diophantischer Approx-

imation, Österreich. Akad. Wiss. Math.-Natur. Kl. Sitzungsber. II 185 (1976), 43–50.

[Hl 79] E. Hlawka, Theorie der Gleichverteilung, Bibliographisches Institut, Mannheim, Wien,

Zürich, 1979.

[II 88] H. Imai and M. Iri, Polygonal approximations of a curve — formulations and algorithms, Computational Morphology — A Computational Geometric Approach to the Analysis of Form (G. T. Toussaint, ed.), North-Holland, Amsterdam-New-York, 1988, pp. 71-86.

[KN 74] L. Kuipers and H. Niederreiter, Uniform Distribution of Sequences, Wiley, New York,

1974.

[LM 73] D. G. Larman and P. McMullen, Arcs with increasing chords, Proc. Cambridge Philos.

Soc. 72 (1973), 205–207.

[MV 75] D. E. McClure and R. A. Vitale, Polygonal approximation of plane convex bodies, J.

Math. Anal. Appl. 51 (1975), 326–358.

[Mü 92] J. Müller, Step by step approximation of plane convex bodies, Arch. Math. (Basel) 58

(1992), 606–610.

[Ni 81] H. Niederreiter, A quasi Monte-Carlo method for the approximate evaluation of the extreme values of a function, Studies in Pure Mathematics (P. Erdös, ed.), Birkhäuser, Basel, 1981, pp. 523–530.

[Ni 83] H. Niederreiter, Quasi-Monte Carlo methods for global optimization, Mathematical Statistics and Applications, Vol. B, Proc. 4^{th} Pannonian Symp., Bad Tatzmannsdorf, 1983 (W. Grossmann, G. C. Pflug, I. Vincze, and W. Wertz, ed.), Reidel, Dordrecht, 1985, pp. 251–267.

[Ni 92] H. Niederreiter, Random Number Generation and Quasi-Monte Carlo Methods, SIAM,

Philadelphia, 1992.

[NP 86] H. Niederreiter and P. Peart, Localization of search in quasi-Monte Carlo methods for

global optimization, SIAM J. Sci. Stat. Comput. 7 (1986), 660–664.

[Ro 94] G. Rote, Curves with increasing chords, Math. Proc. Cambridge Philos. Soc. 115 (1994),

1–12.

[Sch 90] W. Schachinger, Dokumentation von Roboterbahnen, Projektpraktikum, Technische Uni-

versität Wien, 1990.

[So 67] I. M. Sobol, On the distribution of points in a cube and the approximate evaluation of integrals, Zh. Vychisl. Mat. i Mat. Fiz. 7 (1967), 784–802 (Russian); English translation in U.S.S.R. Comput. Math. and Math. Phys. 7 (1967), 86–112.

[Ti 90] R. F. Tichy, Random points on the sphere with applications to numerical analysis, Z.

angew. Math. Mech. 70 (1990), 642–646. Institut für Mathematik, Technische Universität Graz

Steyrergasse 30, A-8010 Graz, Austria