A BRIEF SURVEY OF SPHERICAL INTERPOLATION AND APPROXIMATION METHODS FOR TEXTURE ANALYSIS

In texture analysis there are several instances when mathematical methods of spherical interpolation or approximation are required. Ad hoc adaptions of univariate or bivariate methods to the topology of spherical manifolds usually fail in one way or another. Therefore, this contribution will provide a brief survey of genuinely spherical methods.

In texture analysis there are several instances when mathematical methods of spherical interpolation or approximation are required. Ad hoc adaptions of univariate or bivariate methods to the topology of spherical manifolds usually fail in one way or another. Therefore, this contribution will provide a brief survey of genuinely spherical methods.

MOTIVATION AND INTRODUCTION
The most prominent functions of texture analysis, the orientation density function describing the probability distribution of proper rotations, and the pole density function describing the probability distribution of crystallographic lattice planes are genuinely spherical functions. The orientation density function is defined on the threedimensional projective space equivalent to the upper (lower) unit sphere H =-$4+ IRa, the pole density function is defined on the upper (lower) unit sphere H $3+ c IR3. Thus there are several instances when mathematical methods of spherical interpolation or approximation are required, e.g.
preprocessing of experimental pole figure data, postprocessing of calculated values of a pole density function, display of an orientation density function, density estimation of individual orientation data; Ad hoc adaptions of general multivariate methods (cf. Schumaker, 1976aSchumaker, , 1976bFranke and Schumaker, 1987) to the topology of spherical manifolds usually fail in one way or another. They may be successful provided the domain of the data and of the interpolant/approximant to be constructed is a sufficiently small subset of Sa, or if it suffices that smoothness constraints or boundary conditions are only approximately satisfied for a particular subset of Sd, cf. Dierckx (1984);Gmelig Meyling and Pfluger (1987). However, if the domain of the interpolating or smoothing function to be constructed contains the north (south) pole of the sphere and if smoothness constraints or boundary conditions have to be exactly satisfied genuine spherical methods (cf. Schumaker and Traas, 1991;Traas et al., 1993) are required, because there is no differentiable mapping of the entire sphere to a bounded hyperplane. Therefore, this contribution will provide a brief survey of genuinely spherical methods. References of the topic include Berens et al., 1968;Pawelke, 1972;Butzer et al., 1979;Freeden, 1981;Wahba, 1981Wahba, , 1984Wehrens, 1981;Lawson, 1984;Renka, 1984;Barnhill et al., 1987;Nielson and Ramaraj, 1987;Foley, 1990;Pottmann and Eck, 1990;Schumaker and Traas, 1991;Hoschek and Seemann, 1992;Taijeron et al., 1994. 2 BILINEAR SPHERICAL INTERPOLATION From the times of mechanical control of the texture goniometer two "equi-angular" grids to sample a pole density function have survived and are widely used. They are such that (i) data values Zp > 0 are sampled at sites (0i, Oj)i=l P1, j=l P2 $3+ c IR with Thus, the two grids are related to one another by a shift of half the angular step size in each spherical coordinate.
For a conversion of one equi-angular grid to the other bilinear spherical interpolation seems to provide the most simple and actually effective means. For (0,0)  Obviously, interpolated values are always nonnegative, and the resulting surface is C. Repeated conversion of one grid to the other by bilinear interpolation results in an ever smoothed interpolant with gradually diminishing modes. Spherical bilinear interpolation does not seem appropriate for scattered data, i.e. for data sites arbitrarily distributed on the sphere.  (cf. Butzer and Nessel, 1971) is to approximate a (given) function f by a function with "better" properties which is to be constructed by some smoothing operation on f itself. The smoothing operation may be provided by convolution, which reads for spherical functions The convolution integral (3) is called a singlular integral (cf. Aronszajn, 1950) if the sequence (Zp)taA is a kernel, or since it is to be applied to probability density functions more specifically a nonnegative approximate identity. For a fixed v and a variable u of S a c IR a the functions Zp(VU) usually display a bell-shaped graph rotationally symmetric with respect to v with normalized area under the curve whereby for/9--/90 the mode at v becomes larger and narrower in such a way that the area under the curve near v comes out to be 1. A kernel with some parameter set A, or more precisely a nonnegative approximate identity (Zta)pA is defined by its properties 0 < Zta(t) LI(-1, 1), The last two properties may be used alternatively. The nonnegativity of the kernel implies its uniform boundedness with respect to the L norm.
The name approximate identity originates in the property that the sequence {f, Zta tends uniformly to f as /9 /90 which can be seen from One of the most prominent features of the convolution integral f, Zta is that it preserves the "best" properties of each of its factors, i.e. if Zta is very smooth, then f, Zta will be so, too; and usually, Z; are chosen to be very smooth. Now it is presumed that the data values zp > O, p 1 P, have been sampled from a spherical function f g q(S3) at arbitrarily scattered data sites rp, p P, i.e. (rp, Zp) (rp, flrp)), p 1 P, provides the available data. Then an approximation It is emphasized that f(r; pe) is nonnegative and a density itself of the same class of functions as the approximate identity Zp. The approximant f does not interpolate the data. In fact, it is given as superposition of rotationally invariant decreasing functions each of which centered at a data site rp and multiplied by Zp f(rp). This approximation method may be referred to as "quasi-interpolation" because the interpolation^c onditions are thought of as being relaxed by the requirement that the dependence of f(r;pe) with respect to each particular f(rp) diminishes rapidly with increasing distance of r from rp. Of major interest in practical applications with finite data sets (rp, f(rp)), p P, is not the sequence f , Z, and its convergence but an individual member f f * Z,p with a fixed parameter Pe depending on the total number P of data. This fixed parameter Pe may be interpreted as the "window width" considering the data and controls the degree of smoothing. Therefore the choice of this tuning, parameter is crucial. Several mathematical methods exist to optimize the choice of this parameter. Whatever the choice is, the degree of smoothing is usually large because one generally chooses "very smooth" functions as kernel, e.g. spherical analogues of the Gaussian probability density function (Nikolayev and Ullemeyer, these proceedings). Therefore, repeated application of these methods ("iterated singular integrals") tends to make characteristics of the data disappear.
While the parameter Pe can be optimized for a given approximate identity, the choice of the approximate identity itself is not obvious. Early constructions of approximate identities which are optimal with respect to the rate of convergence are given by Bartlett (1963), Epachenikov (1969).
As an example for S c Ig without theoretical justification, the exponential function The reader may find CosiPoWi (Adam, 1989) a rather peculiar reference to the de la Vall6e Poussin kernel.
Obviously, approximation methods by singular convolution integral generalize easily to hyperspheres of any dimension.

INTERPOLATION BY RADIAL BASIS FUNCTIONS
Employing radial basis functions for the purposes of interpolation/approximation or estimation seems to originate in the earth sciences, cf. Krige, i951; Matheron, 1967;Crain and Bhattacharyya, 1967;Shepard, 1968. Radial interpolation methods are closely related to approximation by singular integrals (cf. Aronszajn, 1950;Butzer and Nessel, 1971) and splines as delivered by variational calculus (cf. Hogervorst, 1994). A recent survey is given by Kansa (1992).
The Gaussian reveals most obviously the close relation to approximation by the singular integral of Gauss-Weierstrass. All of them may be read as appropriately smoothed multivariate versions of univariate piecewise linear interpolation provided by Zp(x) Z(I x-xpl) which produces discontinuities of the first derivative at the data sites xp.
In case of the multiquadrics smoothing of the distance function is done by using arcs of hyperbolas that possess the r values as asymptotes, where c is the tuning parameter which must be supplied by the user. Considering the augmented data sites (x, 0) and (Xp, c) as elements of IRa/l, the multiquadric Zp(x) may be interpreted as the distance between these augmented data sites.
Next, generalizations of the multiquadric and inverse multiquadric approach for scattered spherical data are summarized as they are known for their excellent practical achievements (Foley, 1990;Pottmann and Eck, 1990).
Let gp(r) arccos (rrp) denote the geodesic distance of r from r. Then define f(r) XpZ(g(r)) (26) p=l where Z(t) is either the multiquadric or the inverse multiquadric. As usually, coefficients , p 1 P, are determined by solving the linear system of equations P /(re) _ c2B(g2(rp)) zp, p P p=l Each basis function Z(gp(r)) is a C function on the sphere except at the antipodal point-rp, where it is not differentiable. Thus the interpolant is C(Sa). Modifications of basis functions such that they are C on S a can be accomplished by piecewise quintic poynomial blending near the antipodal point. Yet another spherical generalization (Hardy and Goepfert, 1975;Pottmann and Eck, 1990) is defined as P f(r) _ Az,(r) (32) p=l where Zp(r) /1 + c -2c(rrp) (33) usually with 0 < c < 1, but c > 1 is feasible, too. c 3/8 yields Zp(rp) 5/8; for large P a smaller value of c is to be preferred; an optimal choice of c is unknown.
The linear system of equations to determine the coefficients has always a unique solution (Micchelli, 1986). However, the system may be ill conditioned, in particular for large P, and require preconditioning (cf. Hogervorst, 1994). The interpolant is not necessarily nonnegative. If the nonnegativity constraint is actually violated for one tuning parameter Cl another value c may apply.

APPROXIMATION BY SPHERICAL SPLINES
Polynomial splines may roughly be thought of as piecewise polynomials of a given polynomial order which are joined together to form a unique function at their defining knots. With polynomial interpolation/approximation in mind one may think of the decoupling of the total number of interpolation conditions (data) and the polynomial order of the interpolant as one major motivation of the development of polynomial splines. In accordance with the best known natural cubic splies the term spline is nowadays used for a function that is defined piecewise or satisfies a specific optimality constraint or both. Key references for splines in probability and statistics are provided by Boneva et al., 1971;Wegman and Wright, 1983;Silverman, 1989. This communication is confined to an approximation scheme provided by a tensor product of a polynomial and a periodic trigonometric spline of order rn n 3 (Schumaker and Traas, 1991;Traas et al., 1993).
It is again presumed that the data values zp > 0, p 1 P, have been sampled from a spherical function f q(S3) at arbitrarily scattered data sites r, p P, i.e. (rp, zp) (rp, f(rp)), p = 1 P, provides the available data. Let the function f(r) f(O, qg) be given in terms of spherical coordinates lattitude [0, n:] and longitude 99 [0, 2] of r S3.
Assume the spline approximation f(r) of f(r) is given as M N+2 f(r) f (, p) ,, cijN()Tn(cp) A solution can be obtained be the Householder transformation method. The approximate is C on S3.
This scheme of spline approximation generalizes to a better order of smoothness as well as to hyperspheres in a straightforward manner. Its major advantage is that repeated application does not change the initial spline approximation, because the spline approximation of a given spline is this spline itself. While B-splines are nonnegative functions any linear combination need not to be nonnegative. If the required nonnegativity of f is violated it may be corrected by appropriate subdivision (cf. Dyn et al., 1990) and final correction of corresponding control points or some other numerical heuristics. Promising results have been reported by Traas et al. (1993).

CONCLUSIONS
This communication is merely a brief review of methods of genuinely spherical interpolation and approximation. It is in no way complete, but rather a subjective selection of suggested readings. Any major omission is owed to the author's limited knowledge.